Nonequilibrium dynamics of mixtures of active and passive colloidal particles

TL;DR

This paper develops a mesoscopic field theory to analyze the nonequilibrium dynamics of active and passive colloidal mixtures, predicting instabilities and cluster behaviors that align with simulations.

Contribution

It introduces a novel field theory for active-passive colloidal mixtures and identifies two types of spinodal instabilities, including a Hopf bifurcation leading to moving clusters.

Findings

Good agreement between theory and particle simulations.

Identification of two spinodal instability types.

Prediction of moving clusters in active-passive mixtures.

Abstract

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.

Equations

$$\displaystyle\begin{split}\dot{\vec{r}}^{\mu}_{i}(t)&=\beta D^{\mu}_{\mathrm{T}}\Big{(}\!\vec{F}^{\mu}_{\mathrm{int},i}(\{\vec{r}^{\mu}_{i}\})+F^{\mu}_{\mathrm{A}}\hat{u}(\phi^{\mu}_{i})\Big{)}+\vec{\xi}^{\mu}_{\mathrm{T},i}(t)\;,\end{split}$$

$$\begin{array}{rl}{\dot{\varphi }}_i^{\mu }(t)& ={\xi }_{\mathrm{R},i}^{\mu }(t)\end{array}$$

$$\begin{array}{r}{\vec{F}}_{\mathrm{int},i}^{\mu }(\{{\vec{r}}_i^{\mu }\})=-\underset{(\mu ,i)\ne (\nu ,j)}{\sum _{\nu =1}^n\sum _{j=1}^{N_{\nu }}}{\vec{\nabla }}_{{\vec{r}}_i^{\mu }}{U}_2^{(\mu \nu )}(\Vert {\vec{r}}_i^{\mu }-{\vec{r}}_j^{\nu }\Vert ).\end{array}$$

$$\begin{split}\dot{P}(\{\vec{r}^{\mu}_{i}\},\{\phi^{\mu}_{i}\},t)&=\sum^{n}_{\mu=1}\sum^{N_{\mu}}_{i=1}D^{\mu}_{\mathrm{T}}\vec{\nabla}_{\vec{r}^{\mu}_{i}}\!\cdot\!\Big{(}\vec{\nabla}_{\vec{r}^{\mu}_{i}}P(\{\vec{r}^{\mu}_{i}\},\{\phi^{\mu}_{i}\},t)-\beta\big{(}\vec{F}^{\mu}_{\mathrm{int},i}(\{\vec{r}^{\mu}_{i}\})+F^{\mu}_{\mathrm{A}}\hat{u}(\phi^{\mu}_{i})\big{)}P(\{\vec{r}^{\mu}_{i}\},\{\phi^{\mu}_{i}\},t)\!\Big{)}\\ &\quad\>\!+\sum^{n}_{\mu=1}\sum^{N_{\mu}}_{i=1}D^{\mu}_{\mathrm{R}}\partial^{2}_{\phi^{\mu}_{i}}P(\{\vec{r}^{\mu}_{i}\},\{\phi^{\mu}_{i}\},t)\end{split}$$

$$\begin{split}\rho_{\mu}(\vec{r}^{\mu}_{1},\phi^{\mu}_{1},t)&=N_{\mu}\Bigg{(}\!\underset{(\mu,1)\neq(\nu,j)}{\prod^{n}_{\nu=1}\prod^{N_{\nu}}_{j=1}}\!\!\int_{\mathbb{R}^{2}}\!\!\!\!\!\>\!\mathrm{d}^{2}r^{\nu}_{j}\!\int^{2\pi}_{0}\!\!\!\!\!\!\!\mathrm{d}\phi^{\nu}_{j}\Bigg{)}P(\{\vec{r}^{\mu}_{i}\},\{\phi^{\mu}_{i}\},t)\end{split}$$

$$\begin{split}\dot{\rho}_{\mu}(\vec{r},\phi,t)&=D^{\mu}_{\mathrm{T}}\boldsymbol{\triangle}_{\vec{r}}\rho_{\mu}(\vec{r},\phi,t)+D^{\mu}_{\mathrm{R}}\partial^{2}_{\phi}\rho_{\mu}(\vec{r},\phi,t)-v^{\mu}_{0}\vec{\nabla}_{\vec{r}}\!\cdot\!\big{(}\hat{u}(\phi)\rho_{\mu}(\vec{r},\phi,t)\big{)}\\ &\quad\>\!+\beta D^{\mu}_{\mathrm{T}}\sum^{n}_{\nu=1}\vec{\nabla}_{\vec{r}}\!\cdot\!\!\!\int_{\mathbb{R}^{2}}\!\!\!\!\!\>\!\mathrm{d}^{2}r^{\prime}\,U^{(\mu\nu)\prime}_{2}(\lVert\vec{r}-\vec{r}\hskip 1.0pt^{\prime}\rVert)\frac{\vec{r}-\vec{r}\hskip 1.0pt^{\prime}}{\lVert\vec{r}-\vec{r}\hskip 1.0pt^{\prime}\rVert}\int^{2\pi}_{0}\!\!\!\!\!\!\!\mathrm{d}\phi^{\prime}\,\rho^{(2)}_{\mu\nu}(\vec{r},\vec{r}\hskip 1.0pt^{\prime}\!,\phi,\phi^{\prime}\!,t)\end{split}$$

$$\begin{array}{rl}{\rho }_{\mu \nu }^{(2)}(\vec{r},\vec{r}{}^{\prime },\varphi ,{\varphi }^{\prime },t)& ={\rho }_{\mu }(\vec{r},\varphi ,t){\rho }_{\nu }(\vec{r}{}^{\prime },{\varphi }^{\prime },t)g_{\mu \nu }(\vec{r},\vec{r}{}^{\prime },\varphi ,{\varphi }^{\prime },t).\end{array}$$

$$\underset{\Vert \vec{r}-\vec{r}{}^{\prime }\Vert \to \infty }{\lim }{g}_{\mu \nu }(\vec{r},\vec{r}{}^{\prime },\varphi ,{\varphi }^{\prime },t)=1-\frac{{\delta }_{\mu \nu }}{N_{\mu }},$$

$$\begin{split}\dot{\rho}_{\mu}(\vec{r},\phi,t)&=D^{\mu}_{\mathrm{T}}\boldsymbol{\triangle}_{\vec{r}}\rho_{\mu}(\vec{r},\phi,t)+D^{\mu}_{\mathrm{R}}\partial^{2}_{\phi}\rho_{\mu}(\vec{r},\phi,t)-v^{\mu}_{0}\vec{\nabla}_{\vec{r}}\!\cdot\!\big{(}\hat{u}(\phi)\rho_{\mu}(\vec{r},\phi,t)\big{)}\\ &\quad\>\!+\beta D^{\mu}_{\mathrm{T}}\sum^{n}_{\nu=1}\vec{\nabla}_{\vec{r}}\!\cdot\!\bigg{(}\rho_{\mu}(\vec{r},\phi,t)\!\int_{\mathbb{R}^{2}}\!\!\!\!\!\>\!\mathrm{d}^{2}r^{\prime}\,U^{(\mu\nu)\prime}_{2}(\lVert\vec{r}-\vec{r}\hskip 1.0pt^{\prime}\rVert)\frac{\vec{r}-\vec{r}\hskip 1.0pt^{\prime}}{\lVert\vec{r}-\vec{r}\hskip 1.0pt^{\prime}\rVert}\!\>\!\int^{2\pi}_{0}\!\!\!\!\!\!\!\mathrm{d}\phi^{\prime}\,g_{\mu\nu}(\vec{r},\vec{r}\hskip 1.0pt^{\prime}\!,\phi,\phi^{\prime}\!,t)\rho_{\nu}(\vec{r}\hskip 1.0pt^{\prime}\!,\phi^{\prime}\!,t)\!\bigg{)}\>\!.\end{split}$$

$$\begin{split}\dot{\rho}_{\mu}(\vec{r},\phi,t)&=D^{\mu}_{\mathrm{T}}\boldsymbol{\triangle}_{\vec{r}}\rho_{\mu}(\vec{r},\phi,t)+D^{\mu}_{\mathrm{R}}\partial^{2}_{\phi}\rho_{\mu}(\vec{r},\phi,t)-\vec{\nabla}_{\vec{r}}\!\cdot\!\big{(}v_{\mu}(\vec{r},\phi,t)\hat{u}(\phi)\rho_{\mu}(\vec{r},\phi,t)\big{)}\\ &\quad\>\!+\beta D^{\mu}_{\mathrm{T}}\vec{\nabla}_{\vec{r}}\!\cdot\!\big{(}\rho_{\mu}(\vec{r},\phi,t)\vec{\nabla}_{\vec{r}}\Upsilon_{\mu}(\vec{r},\phi,t)\big{)}\end{split}$$

$$\begin{array}{r}{v}_{\mu }(\vec{r},\varphi ,t)=v_0^{\mu }-\beta D_{\mathrm{T}}^{\mu }{F}_{\mathrm{res}}^{\mu }(\vec{r},\varphi ,t)\end{array}$$

$$\begin{split}F^{\mu}_{\mathrm{res}}(\vec{r},\phi,t)&=\sum^{n}_{\nu=1}\sum^{\infty}_{m=1}\!\!\!\sum^{\lfloor\frac{m-1}{2}\rfloor}_{k=0}\!\!\sum^{k}_{l=0}\bigg{(}\frac{A^{(\mu\nu)}_{m,2k}}{(m-2k)!(2k)!}-\frac{A^{(\mu\nu)}_{m,2k+2}}{(m-2k-1)!(2k+1)!}\bigg{)}\frac{(-1)^{l}k!}{(k-l)!l!}\boldsymbol{\triangle}_{\vec{r}}^{k-l}\big{(}\hat{u}(\phi)\!\cdot\!\vec{\nabla}_{\vec{r}}\big{)}^{m-2(k-l)}\rho_{\nu}(\vec{r},t)\\ &\quad\>\!+\sum^{n}_{\nu=1}\sum^{\infty}_{m=0}\sum^{m}_{l=0}\frac{A^{(\mu\nu)}_{2m,2m}}{(2m)!}\frac{(-1)^{l}m!}{(m-l)!l!}\boldsymbol{\triangle}_{\vec{r}}^{m-l}\big{(}\hat{u}(\phi)\!\cdot\!\vec{\nabla}_{\vec{r}}\big{)}^{2l}\rho_{\nu}(\vec{r},t)\end{split}$$

$$\begin{split}\Upsilon_{\mu}(\vec{r},\phi,t)&=\sum^{n}_{\nu=1}\sum^{\infty}_{m=1}\!\!\!\sum^{\lfloor\frac{m-1}{2}\rfloor}_{k=0}\!\!\sum^{k}_{l=0}\frac{A^{(\mu\nu)}_{m,2k+2}}{(m-2k-1)!(2k+1)!}\frac{(-1)^{l}k!}{(k-l)!l!}\boldsymbol{\triangle}_{\vec{r}}^{k-l}\big{(}\hat{u}(\phi)\!\cdot\!\vec{\nabla}_{\vec{r}}\big{)}^{m-2(k-l)-1}\rho_{\nu}(\vec{r},t)\end{split}$$

$$\begin{array}{rl}{A}_{m,k}^{(\mu \nu )}& =-{\int }_0^{\infty }\mathrm{d}r{r}^{m+1}{U}_2^{(\mu \nu )\prime }(r){\int }_0^{2\pi }\mathrm{d}{\varphi }_{\mathrm{R}}{g}_{\mu \nu }(r,{\varphi }_{\mathrm{R}})\cos ({\varphi }_{\mathrm{R}}{)}^{m-k+1}\sin ({\varphi }_{\mathrm{R}}{)}^k\end{array}$$

$$\dot{\rho}_{\mu}=\sum^{n}_{\nu=1}\partial_{i}\Big{(}D^{(\mu\nu)}(\{\rho_{\mu}\})\partial_{i}\rho_{\nu}\Big{)}$$

$$\begin{split}D^{(\mu\nu)}(\{\rho_{\mu}\})&=\bigg{(}D^{\mu}_{\mathrm{T}}+\frac{v^{2}_{\mu}}{2D^{\mu}_{\mathrm{R}}}\bigg{)}\delta_{\mu\nu}+\frac{\rho_{\mu}}{2}\bigg{(}a^{(\mu\nu)}_{1}-a^{(\mu\nu)}_{0}\frac{v^{\mu}_{0}v_{\mu}}{D^{\mu}_{\mathrm{R}}}\bigg{)},\end{split}$$

$$v_{\mu}(\{\rho_{\mu}\})=v^{\mu}_{0}\bigg{(}1-\sum^{n}_{\nu=1}a^{(\mu\nu)}_{0}\rho_{\nu}\bigg{)},$$

$$\begin{array}{r}{a}_0^{(\mu \nu )}=\beta D_{\mathrm{T}}^{\mu }\frac{A_{0,0}^{(\mu \nu )}}{v_0^{\mu }},\end{array}$$

$$\begin{array}{r}{a}_1^{(\mu \nu )}=\beta D_{\mathrm{T}}^{\mu }(A_{1,0}^{(\mu \nu )}+A_{1,2}^{(\mu \nu )}).\end{array}$$

$$\begin{array}{rl}{a}_0^{(\mu \nu )}& =-\frac{1}{v_0^{\mu }}\frac{\partial v_{\mu }}{\partial {\rho }_{\nu }},\end{array}$$

$$\displaystyle\begin{split}a^{(\mu\nu)}_{1}&=\frac{2}{\rho_{\mu}}\bigg{(}D^{(\mu\nu)}-\bigg{(}D^{\mu}_{\mathrm{T}}+\frac{v^{2}_{\mu}}{2D^{\mu}_{\mathrm{R}}}\bigg{)}\delta_{\mu\nu}\bigg{)}-\frac{v_{\mu}}{D^{\mu}_{\mathrm{R}}}\frac{\partial v_{\mu}}{\partial\rho_{\nu}}\;.\end{split}$$

$${\dot{\delta \rho }}_{\mu }=\sum _{\nu =1}^n\mathrm{bar}{D}_{\mu \nu }△\delta {\rho }_{\nu }$$

$${\dot{\widetilde{\delta \rho }}}_{\mu }=-k^2\sum _{\nu =1}^n\mathrm{bar}{D}_{\mu \nu }{\widetilde{\delta \rho }}_{\nu },$$

$$\vec{\widetilde{\delta\rho}}(\vec{k},t)=\exp\!\big{(}\!-k^{2}\boldsymbol{\mathrm{bar}{\mathrm{D}}}t\big{)}\vec{\widetilde{\delta\rho}}(\vec{k},0)$$

$$\mathrm{Re}({\lambda }_i)⩾0\forall i\in \{1,\dots ,n\},$$

$$\mathrm{bar}\mathrm{D}=\left(\begin{array}{cc}\mathrm{bar}{D}_{\mathrm{AA}}& \mathrm{bar}{D}_{\mathrm{AB}}\\ \mathrm{bar}{D}_{\mathrm{BA}}& \mathrm{bar}{D}_{\mathrm{BB}}\end{array}\right)$$

$$\mathrm{bar}{D}_{\mathrm{AA}}=-\mathrm{bar}{D}_{\mathrm{BB}}\wedge \mathrm{bar}{D}_{\mathrm{AA}}^2<-\mathrm{bar}{D}_{\mathrm{AB}}\mathrm{bar}{D}_{\mathrm{BA}}$$

$$\mathrm{bar}{D}_{\mathrm{AA}}\mathrm{bar}{D}_{\mathrm{BB}}=\mathrm{bar}{D}_{\mathrm{AB}}\mathrm{bar}{D}_{\mathrm{BA}},$$

$$\begin{array}{rl}{\lambda }_1& =\frac{1}{2}(\mathrm{bar}{D}_{\mathrm{AA}}+\mathrm{bar}{D}_{\mathrm{BB}})-\frac{1}{2}\sqrt{(\mathrm{bar}{D}_{\mathrm{AA}}-\mathrm{bar}{D}_{\mathrm{BB}}{)}^2+4\mathrm{bar}{D}_{\mathrm{AB}}\mathrm{bar}{D}_{\mathrm{BA}}},\end{array}$$

$$\begin{array}{rl}{\lambda }_2& =\frac{1}{2}(\mathrm{bar}{D}_{\mathrm{AA}}+\mathrm{bar}{D}_{\mathrm{BB}})+\frac{1}{2}\sqrt{(\mathrm{bar}{D}_{\mathrm{AA}}-\mathrm{bar}{D}_{\mathrm{BB}}{)}^2+4\mathrm{bar}{D}_{\mathrm{AB}}\mathrm{bar}{D}_{\mathrm{BA}}}\end{array}$$

$${\vec{w}}_1=\left(\begin{array}{c}{\lambda }_1-\mathrm{bar}{D}_{\mathrm{BB}}\\ \mathrm{bar}{D}_{\mathrm{BA}}\end{array}\right),{\vec{w}}_2=\left(\begin{array}{c}{\lambda }_2-\mathrm{bar}{D}_{\mathrm{BB}}\\ \mathrm{bar}{D}_{\mathrm{BA}}\end{array}\right)$$

$$\overrightarrow{\delta \rho }(\vec{r},t)=\sum _{i=1}^2{A}_i{e}^{\mathrm{i}({\vec{k}}_i\cdot \vec{r}-{\omega }_{i}t)}{\vec{w}}_i$$

Full text

active colloidal particles, active-passive mixtures, motility-induced instability, mesoscopic field theory, particle-resolved simulations

I 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 Leptos et al. (2009); Lin et al. (2011); Jepson et al. (2013). 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 active-particle suspension Angelani et al. (2011); Ni et al. (2015); Tanaka et al. (2017), crystallization and melting of hard-sphere glasses by doping with active particles Ni et al. (2014); Kümmel et al. (2015), and mesoscale turbulence mediated by passive particles Hinz et al. (2014). Furthermore, the intriguing phenomenon of motility-induced phase separation (MIPS) Cates and Tailleur (2015), 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 Stenhammar et al. (2015); Wysocki et al. (2016). For purely active systems, where MIPS has already been extensively studied Tailleur and Cates (2008); Fily and Marchetti (2012); Cates and Tailleur (2013); Bialké et al. (2013); Buttinoni et al. (2013); Stenhammar et al. (2013); Redner et al. (2013); Stenhammar et al. (2014); Speck et al. (2014); Wittkowski et al. (2014); Wysocki et al. (2014); Takatori et al. (2014); Solon et al. (2015); Tiribocchi et al. (2015), 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 Stenhammar et al. (2015); Wysocki et al. (2016), 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-correlation 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 one-component 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 Stenhammar et al. (2015); Wysocki et al. (2016).

The article in organized as follows: in Sec. II we derive dynamic equations for multicomponent mixtures of particles with different motilities as well as useful approximations of these equations. Based on these, in Sec. III 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 Sec. IV. We finally conclude in Sec. V.

II Dynamic equations for mixtures

We consider an $n$-component 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” Bechinger et al. (2016)) undergo translational and rotational diffusion, respectively. In the following, $n$ denotes the number of different species of active or passive colloidal particles and $N_{\mu}$ is the total number of particles of species $\mu\in\{1,\dotsc,n\}$. The translational and rotational diffusion coefficients of the particles are $D^{\mu}_{\mathrm{T}}$ and $D^{\mu}_{\mathrm{R}}$, respectively, and the particles’ motilities are set by active driving forces $F^{\mu}_{\mathrm{A}}\geqslant 0$ ten Hagen et al. (2015), which are zero for passive particles. $\vec{r}^{\mu}_{i}(t)$ is the position and $\phi^{\mu}_{i}(t)$ is the orientation of the $i$th particle of species $\mu$ at time $t$. The interactions between a particle of species $\mu$ and a particle of species $\nu$ are described by the pair-interaction potentials $U^{(\mu\nu)}_{2}(\lVert\vec{r}^{\mu}_{i}\!-\vec{r}^{\nu}_{j}\rVert)$, where $U^{(\mu\nu)}_{2}(r)$ and $U^{(\nu\mu)}_{2}(r)$ can be different111This is the case especially for colloidal particles that catalyze chemical reactions on their surfaces and breaks 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 Soto and Golestanian (2014). However, in this article we do not further study such systems. Our simulations described in Sec. IV use the same interaction potential for all particles. and $\lVert\,\cdot\,\rVert$ is the Euclidean norm.

II.1 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 equations222In the limit of a one-component system, these Langevin equations reduce to the simpler ones that have been used in Refs. Redner et al. (2013); Bialké et al. (2013).

[TABLE]

with the interaction forces

[TABLE]

Here, an overdot denotes differentiation with respect to time and $\beta=1/(k_{\mathrm{B}}T)$ is the inverse thermal energy with the Boltzmann constant $k_{\mathrm{B}}$ and the absolute temperature $T$ of the implicit solvent. Furthermore, $\hat{u}(\phi)=(\cos(\phi),\sin(\phi))^{\mathrm{T}}$ is a normalized orientation vector and $\vec{\nabla}_{\vec{r}^{\mu}_{i}}=(\partial_{x^{\mu}_{i}},\partial_{y^{\mu}_{i}})$ is the del symbol that involves partial derivatives with respect to the elements of the vector $\vec{r}^{\mu}_{i}=(x^{\mu}_{i},y^{\mu}_{i})^{\mathrm{T}}$. The elements of the translational noise $\vec{\xi}^{\mu}_{\mathrm{T},i}(t)$ and the rotational noise $\xi^{\mu}_{\mathrm{R},i}(t)$ in the Langevin equations are statistically independent Gaussian white noises with mean values $\langle\vec{\xi}^{\mu}_{\mathrm{T},i}(t)\rangle=\vec{0}$ and $\langle\xi^{\mu}_{\mathrm{R},i}(t)\rangle=0$, and correlations $\langle\vec{\xi}^{\mu}_{\mathrm{T},i}(t_{1})\!\otimes\!\vec{\xi}^{\nu}_{\mathrm{T},j}(t_{2})\rangle=2D^{\mu}_{\mathrm{T}}\mathbf{1}\delta_{ij}\delta_{\mu\nu}\delta(t_{1}-t_{2})$ and $\langle\xi^{\mu}_{\mathrm{R},i}(t_{1})\xi^{\nu}_{\mathrm{R},j}(t_{2})\rangle=2D^{\mu}_{\mathrm{R}}\delta_{ij}\delta_{\mu\nu}\delta(t_{1}-t_{2})$, where $\otimes$ denotes the dyadic product and $\mathbf{1}$ is the identity matrix.

The collective dynamics of the $n$-component mixture can equivalently be described by the Smoluchowski equation (see Refs. Dhont (1996); Risken (1996) for details)

[TABLE]

with the many-particle probability density $P(\{\vec{r}^{\mu}_{i}\},\{\phi^{\mu}_{i}\},t)$. Here, the operator $\partial_{\phi^{\mu}_{i}}$ means partial differentiation with respect to $\phi^{\mu}_{i}$. By integrating away all degrees of freedom except for $\vec{r}^{\mu}_{1}$ and $\phi^{\mu}_{1}$, the one-particle densities Wittkowski and Löwen (2011)

[TABLE]

can be obtained from $P(\{\vec{r}^{\mu}_{i}\},\{\phi^{\mu}_{i}\},t)$. Many-particle densities like the two-particle densities $\rho^{(2)}_{\mu\nu}(\vec{r},\vec{r}\hskip 1.0pt^{\prime}\!,\phi,\phi^{\prime}\!,t)$ can be obtained analogously Wittkowski et al. (2012). Applying the integration (5) to the Smoluchowski equation (4) and renaming $\vec{r}^{\mu}_{1}\to\vec{r}=(x,y)^{\mathrm{T}}$ and $\phi^{\mu}_{1}\to\phi$ leads to the dynamic equations

[TABLE]

for the one-particle densities $\rho_{\mu}(\vec{r},\phi,t)$. Here, $v^{\mu}_{0}=\beta D^{\mu}_{\mathrm{T}}F^{\mu}_{\mathrm{A}}$ denotes the speed of a free particle of species $\mu$ and we used the notation $U^{(\mu\nu)\prime}_{2}(r)=\mathrm{d}U^{(\mu\nu)}_{2}(r)/\mathrm{d}r$. The dynamics of the one-particle densities $\rho_{\mu}(\vec{r},\phi,t)$ depends on the two-particle densities $\rho^{(2)}_{\mu\nu}(\vec{r},\vec{r}\hskip 1.0pt^{\prime}\!,\phi,\phi^{\prime}\!,t)$, which can be decomposed into the one-particle densities $\rho_{\mu}(\vec{r},\phi,t)$ and the pair-correlation functions $g_{\mu\nu}(\vec{r},\vec{r}\hskip 1.0pt^{\prime}\!,\phi,\phi^{\prime}\!,t)$:

[TABLE]

Notice that this standard definition of $g_{\mu\nu}(\vec{r},\vec{r}\hskip 1.0pt^{\prime}\!,\phi,\phi^{\prime}\!,t)$ implies that

[TABLE]

where the term $\delta_{\mu\nu}/N_{\mu}$ is usually negligibly small. As a further consequence of Eq. (7), the conventional radial distribution functions cannot be simply obtained by an angular integration and a time average of $g_{\mu\nu}(\vec{r},\vec{r}\hskip 1.0pt^{\prime}\!,\phi,\phi^{\prime}\!,t)$. Instead, they need to be calculated from $\rho^{(2)}_{\mu\nu}(\vec{r},\vec{r}\hskip 1.0pt^{\prime}\!,\phi,\phi^{\prime}\!,t)$.

Inserting Eq. (7) into Eq. (6) results in the exact dynamic equations for the one-particle densities

[TABLE]

The pair-correlation functions $g_{\mu\nu}(\vec{r},\vec{r}\hskip 1.0pt^{\prime}\!,\phi,\phi^{\prime}\!,t)$, however, are unknown and have to be approximated appropriately.

II.2 Approximate dynamic equations

In order to approximate the pair-correlation functions, we first assume global translational invariance of the system, leading to $g_{\mu\nu}(\vec{r},\vec{r}\hskip 1.0pt^{\prime}\!,\phi,\phi^{\prime}\!,t)=g_{\mu\nu}(\vec{r}\hskip 1.0pt^{\prime}\!-\vec{r},\phi,\phi^{\prime}\!,t)$. With the parametrization $\vec{r}\hskip 1.0pt^{\prime}-\vec{r}=r\hat{u}(\phi_{\mathrm{R}})$ (see Fig. 1) this can be rewritten as $g_{\mu\nu}(r,\phi_{\mathrm{R}},\phi,\phi^{\prime}\!,t)$.

Assuming also invariance of the system with respect to global rigid rotations, the number of variables of the pair-correlation functions can be reduced further, yielding $g_{\mu\nu}(r,\phi_{\mathrm{R}}-\phi,\phi^{\prime}\!-\phi,t)$. In addition, we assume that the pair-correlation functions are not explicitly time-dependent: $g_{\mu\nu}(r,\phi_{\mathrm{R}}-\phi,\phi^{\prime}\!-\phi)$. These common approximations become exact when the system is in a homogeneous steady state. In accordance with Ref. Bialké et al. (2013), we now further approximate $g_{\mu\nu}(r,\phi_{\mathrm{R}}-\phi,\phi^{\prime}\!-\phi)$ by $g_{\mu\nu}(r,\phi_{\mathrm{R}}-\phi)$ with the symmetry property $g_{\mu\nu}(r,-\theta)=g_{\mu\nu}(r,\theta)$. Notice that the pair-correlation functions $g_{\mu\nu}(r,\phi_{\mathrm{R}}-\phi)$ still depend on the driving forces $F^{\mu}_{\mathrm{A}}$ of the active particles and on the constant average particle number densities $\mathrm{bar}{\rho}_{\mu}=\langle\rho_{\mu}(\vec{r},t)\rangle_{\vec{r}}$, where $\rho_{\mu}(\vec{r},t)=\int^{2\pi}_{0}\!\!\mathrm{d}\phi\,\rho_{\mu}(\vec{r},\phi,t)$ are the orientation-independent translational densities (i.e., concentration fields) and $\langle\,\cdot\,\rangle_{\vec{r}}$ denotes a spatial average. Regarding the pair-interaction potentials $U^{(\mu\nu)}_{2}(r)$ we assume short-range interactions between the particles. This allows us to simplify the final term of Eq. (9) through a gradient expansion Yang et al. (1976); Evans (1979); Emmerich et al. (2012).

With these approximations, Eq. (9) simplifies to the dynamic equations for $n$-component mixtures of active and passive particles

[TABLE]

with the effective speeds of the active particles

[TABLE]

that describe the slow-down of active particles due to collisions Cates and Tailleur (2013); Bialké et al. (2013). Here, the resistive forces

[TABLE]

and the functions

[TABLE]

with the coefficients

[TABLE]

are gradient expansions of the translational densities $\rho_{\mu}(\vec{r},t)$. The four terms on the right-hand side of Eq. (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 $\Upsilon_{\mu}(\vec{r},\phi,t)$ reduces to the excess chemical potential corresponding to the particles of species $\mu$. 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-correlation 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 $g_{\mu\nu}(r,\phi_{\mathrm{R}}-\phi,\phi^{\prime}\!-\phi)$ on $\phi^{\prime}\!-\phi$. 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 $\phi^{\prime}\!-\phi$-dependence of the pair-correlation functions, which depend only on $r$ when the motilities of all particles vanish. With additional assumptions on the pair-correlation functions, one can reduce the number of independent coefficients $A^{(\mu\nu)}_{m,k}$ in Eqs. (12) and (13). These possible further approximations of the pair-correlation functions are described in Appendix A, but not used in this article.

We proceed by expanding the orientation-dependent one-particle densities $\rho_{\mu}(\vec{r},\phi,t)$ with respect to the orientation $\phi$, 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

[TABLE]

for the translational densities $\rho_{\mu}(\vec{r},t)$. Here and in the following, $\partial_{i}$ denotes the $i$th element of the del symbol $\vec{\nabla}_{\vec{r}}=(\partial_{1},\partial_{2})=(\partial_{x},\partial_{y})$ and Einstein’s summation convention is used for Latin but not for Greek indices. Furthermore, the density-dependent diffusion tensor $\boldsymbol{\mathrm{D}}(\{\rho_{\mu}\})$ has the elements

[TABLE]

the density-dependent swim speeds $v_{\mu}(\{\rho_{\mu}\})$ are given by

[TABLE]

and the coefficients $a^{(\mu\nu)}_{0}$ and $a^{(\mu\nu)}_{1}$ are defined as333For $v^{\mu}_{0}\to 0$ the pair-correlation functions $g_{\mu\nu}(r,\phi)$ become independent of $\phi$ and the coefficients $A^{(\mu\nu)}_{0,0}$ vanish. Therefore, the coefficients $a^{(\mu\nu)}_{0}$ do not diverge when $v^{\mu}_{0}\to 0$.

[TABLE]

The division of $A^{(\mu\nu)}_{0,0}$ by $v^{\mu}_{0}$ in the definition of the coefficients $a^{(\mu\nu)}_{0}$ is motivated by the fact that the effective swim speeds $v_{\mu}(\{\rho_{\mu}\})$ are proportional to $v^{\mu}_{0}$ in mixtures of active and passive particles Stenhammar et al. (2015). Equations (15) and (16) constitute the second main result of this article.

On the right-hand side of Eq. (16), the term $D^{\mu}_{\mathrm{T}}\delta_{\mu\nu}$ describes the diffusivity of a free passive particle, the term $a^{(\mu\nu)}_{1}\rho_{\mu}/2$ takes interactions with particles of the same or other species into account, and the two remaining contributions describe the effect of the particles’ motilities on the collective diffusion. Since the diffusion tensor $\boldsymbol{\mathrm{D}}(\{\rho_{\mu}\})$ and the swim speeds $v_{\mu}(\{\rho_{\mu}\})$ can directly be obtained from particle-resolved computer simulations, Eqs. (16) and (17) allow to determine the coefficients $a^{(\mu\nu)}_{0}$ and $a^{(\mu\nu)}_{1}$:

[TABLE]

Especially for one-component systems, this constitutes an alternative to determining the pair-correlation functions in simulations and calculating $a^{(\mu\nu)}_{0}$ and $a^{(\mu\nu)}_{1}$ using Eqs. (14), (18), and (19).

III Dynamical instability

For certain combinations of the free swim speeds $v^{\mu}_{0}$ and the average particle number densities $\mathrm{bar}{\rho}_{\mu}$ the mixture becomes unstable, particle clusters form, and the system phase separates Stenhammar et al. (2015). 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.

III.1 Multicomponent mixtures

To derive the spinodal, we consider a system with slightly perturbed translational densities $\rho_{\mu}(\vec{r},t)=\mathrm{bar}{\rho}_{\mu}+\delta\rho_{\mu}(\vec{r},t)$, where $\delta\rho_{\mu}(\vec{r},t)$ are the small perturbations. Inserting this perturbation ansatz into Eq. (15) and linearizing the resulting expression in $\delta\rho_{\mu}(\vec{r},t)$ leads to the dynamic equations for the perturbations

[TABLE]

with the constant diffusion coefficients $\mathrm{bar}{D}_{\mu\nu}=D^{(\mu\nu)}(\{\mathrm{bar}{\rho}_{\mu}\})$ corresponding to the average densities $\mathrm{bar}{\rho}_{1},\dotsc,\mathrm{bar}{\rho}_{n}$. Fourier transforming Eq. (22) results in

[TABLE]

where $\widetilde{\delta\rho}_{\mu}(\vec{k},t)$ denotes the spatial Fourier transform of $\delta\rho_{\mu}(\vec{r},t)$ and $k=\lVert\vec{k}\rVert$ is the wave number corresponding to the wave vector $\vec{k}$. If we define the perturbation vector $\vec{\widetilde{\delta\rho}}=(\widetilde{\delta\rho}_{1},\dotsc,\widetilde{\delta\rho}_{n})^{\mathrm{T}}$ and the constant diffusion tensor $\boldsymbol{\mathrm{bar}{\mathrm{D}}}=\boldsymbol{\mathrm{D}}(\{\mathrm{bar}{\rho}_{\mu}\})=(\mathrm{bar}{D}_{\mu\nu})_{\mu,\nu=1,\dotsc,n}$, the solution of Eq. (23) can be written as

[TABLE]

with the matrix exponential $\exp(\!-k^{2}\boldsymbol{\mathrm{bar}{\mathrm{D}}}t)$. An $n$-component mixture of colloidal particles with different motilities is therefore stable if the real parts of all eigenvalues $\lambda_{1},\dotsc,\lambda_{n}$ of $\boldsymbol{\mathrm{bar}{\mathrm{D}}}$ are nonnegative, i.e.,

[TABLE]

and unstable otherwise (see Appendix C for details). The spinodal of the dynamical instability is the outer surface that encloses the surfaces which are defined by the implicit functions $\operatorname{Re}(\lambda_{i})=0$ with $i\in\{1,\dotsc,n\}$. In the next subsection, we apply the stability criterion (25) to binary mixtures.

III.2 Binary mixtures

We now focus on binary mixtures with two species $\mu\in\{\mathrm{A},\mathrm{B}\}$ of particles. The constant diffusion tensor $\boldsymbol{\mathrm{bar}{\mathrm{D}}}$ is then a real $2\!\times\!2$-dimensional matrix

[TABLE]

and the spinodal is the outer surface that encloses the two nonoverlapping444Combining Eqs. (27) and (28) leads to the unfulfillable condition $\mathrm{bar}{D}^{2}_{\mathrm{AA}}<\mathrm{bar}{D}^{2}_{\mathrm{AA}}$. This shows that the surfaces described by Eqs. (27) and (28) have no points in common. surfaces given by

[TABLE]

and

[TABLE]

respectively. Equations (27) and (28) describe the stability of binary mixtures of colloidal particles with different motilities and constitute the third main result of this article.

As can be seen from the eigenvalues

[TABLE]

and corresponding eigenvectors

[TABLE]

of $\boldsymbol{\mathrm{bar}{\mathrm{D}}}$, respectively, Eqs. (27) and (28) correspond to instabilities of different types. Under the condition (27), the eigenvalues $\lambda_{1}$ and $\lambda_{2}$ are imaginary and complex conjugates of each other. The surface described by Eq. (27), where a complex conjugate pair of eigenvalues of $\boldsymbol{\mathrm{bar}{\mathrm{D}}}$ passes through the imaginary axis, is associated with a Hopf bifurcation. 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 (28), $\lambda_{1}$ and $\lambda_{2}$ are real and one of these eigenvalues is zero. The surface described by Eq. (28), where one of the two real eigenvalues of $\boldsymbol{\mathrm{bar}{\mathrm{D}}}$ passes through zero, is associated with a stationary bifurcation. This type of bifurcation is possible also for one-component systems. Hence, there is a fundamental difference between the already widely studied motility-induced phase separation 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.

To study the features of these bifurcations, we now insert the plane-wave ansatz

[TABLE]

with the perturbation vector $\vec{\delta\rho}=(\delta\rho_{\mathrm{A}},\delta\rho_{\mathrm{B}})^{\mathrm{T}}$, constants $A_{1},A_{2}\in\mathbb{C}$, imaginary unit $\mathrm{i}$, wave vectors $\vec{k}_{1},\vec{k}_{2}\in\mathbb{R}^{2}$, and angular frequencies $\omega_{1},\omega_{2}\in\mathbb{C}$ into Eq. (22).555Note that in the plane-wave ansatz (32) the perturbation vector $\vec{\delta\rho}$ is complex. The perturbations that could be observed in simulations or experiments correspond to the real part of $\vec{\delta\rho}$. This results in the dispersion relations

[TABLE]

with $k_{1}=\lVert\vec{k}_{1}\rVert$ and $k_{2}=\lVert\vec{k}_{2}\rVert$. When decomposing the eigenvalues $\lambda_{i}=\operatorname{Re}(\lambda_{i})+\mathrm{i}\operatorname{Im}(\lambda_{i})$ with $i\in\{1,2\}$ into their real parts $\operatorname{Re}(\lambda_{i})$ and imaginary parts $\operatorname{Im}(\lambda_{i})$, where $\operatorname{Im}(\lambda_{1})=-\operatorname{Im}(\lambda_{2})$, and inserting the dispersion relations (33) and (34) into the ansatz (32), we obtain

[TABLE]

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 of $\boldsymbol{\mathrm{bar}{\mathrm{D}}}$ has a negative real part. Furthermore, we see that the perturbations are traveling plane waves for $\operatorname{Im}(\lambda_{1})\neq 0$, i.e., for the Hopf bifurcation associated with Eq. (27), whereas the perturbations are static for $\operatorname{Im}(\lambda_{1})=0$, as holds for the stationary bifurcation associated with Eq. (28). 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 $\operatorname{Im}(\lambda_{2})k$. Note that this expression holds only for small $k$, since the continuum description by the field theory (15) breaks down for $k\gtrsim 2\pi/(R_{\mathrm{A}}+R_{\mathrm{B}})$, where $R_{\mathrm{A}}$ and $R_{\mathrm{B}}$ are the radii of particles of species A and B, respectively.

III.2.1 Active-passive mixtures

In the following, we assume that all particles have the same radius $R$ and thus identical translational diffusion coefficients $D_{\mathrm{T}}=D^{\mathrm{A}}_{\mathrm{T}}=D^{\mathrm{B}}_{\mathrm{T}}$ as well as identical rotational diffusion coefficients $D_{\mathrm{R}}=D^{\mathrm{A}}_{\mathrm{R}}=D^{\mathrm{B}}_{\mathrm{R}}$. The elements of the diffusion tensor (26) then simplify to

[TABLE]

with the effective swim speeds

[TABLE]

In the special case of a mixture of active (A) and passive (P) particles with $\mu\in\{\mathrm{A},\mathrm{P}\}$, Eqs. (36)-(41) further simplify due to $\mathrm{bar}{v}_{\mathrm{P}}=v^{\mathrm{P}}_{0}=0$. With the notation $\mathrm{bar}{v}=\mathrm{bar}{v}_{\mathrm{A}}$ and $v_{0}=v^{\mathrm{A}}_{0}$ this results in

[TABLE]

and

[TABLE]

Equations (42)-(45) show that for purely repulsive interactions, where $a^{(\mathrm{PA})}_{1}>0$ and $a^{(\mathrm{PP})}_{1}>0$, the passive particles are always stable and that the clustering is triggered by a dynamical instability of the active particles. According to Eqs. (27) and (28), the spinodal for the instability of the active-passive mixture is defined by

[TABLE]

and

[TABLE]

III.2.2 Large $v_{0}$

We now consider large values $v_{0}\to\infty$ of the self-propulsion speed of a free particle $v_{0}$, but not too large particle densities $\mathrm{bar}{\rho}_{\mathrm{A}}$ and $\mathrm{bar}{\rho}_{\mathrm{P}}$ so that $\mathrm{bar}{v}\to\infty$. Assuming that the coefficients $a^{(\mathrm{AA})}_{0}$, $a^{(\mathrm{AP})}_{0}$, $a^{(\mathrm{AA})}_{1}$, $a^{(\mathrm{AP})}_{1}$, $a^{(\mathrm{PA})}_{1}$, and $a^{(\mathrm{PP})}_{1}$ neither vanish nor diverge when $v_{0}\to\infty$, Eq. (47) then no longer has a solution, whereas Eq. (48) simplifies to

[TABLE]

This means that only the stationary bifurcation but not the Hopf bifurcation occurs in the limit $v_{0}\to\infty$. If we set $a^{(\mathrm{PA})}_{1}=a^{(\mathrm{PP})}_{1}=0$ in Eq. (49), it further reduces to

[TABLE]

which is a simple condition for the low-density branch of the spinodal, as found previously in Ref. Stenhammar et al. (2015).

III.3 One-component systems

Another limiting case of Eqs. (47) and (48) is that of a one-component system of active particles where $\mathrm{bar}{\rho}_{\mathrm{P}}=0$. Also in this case, Eq. (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 Eq. (48) the spinodal condition

[TABLE]

with the effective speed of the active particles

[TABLE]

When choosing $a^{(\mathrm{AA})}_{1}=0$, Eq. (51) becomes equivalent to the simpler spinodal condition that has been proposed in Ref. Bialké et al. (2013). This approximation is, however, in general not applicable. The term proportional to $a^{(\mathrm{AA})}_{1}$ in the diffusion coefficient $\mathrm{bar}{D}_{\mathrm{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}$. This value depends on $\mathrm{bar}{\rho}_{\mathrm{A}}$: For small values of $\mathrm{bar}{\rho}_{\mathrm{A}}$, the instability cannot occur at all. Considering $a^{(\mathrm{AA})}_{0}$ and $a^{(\mathrm{AA})}_{1}$ as nonzero finite constants, which is – as we will see in Sec. IV.2 – only an approximation, the minimal swim speed $v_{0,\mathrm{min}}$ of a free particle that allows for cluster formation is given by

[TABLE]

and corresponds to the active-particle density

[TABLE]

This result for the minimal value of $v_{0}$ at the spinodal generalizes the previous result $v_{0,\mathrm{min}}=4\sqrt{D_{\mathrm{T}}D_{\mathrm{R}}}$ of Ref. Bialké et al. (2013), which is obtained from Eq. (53) when choosing $a^{(\mathrm{AA})}_{1}=0$ and has been shown to lead to too small values for $v_{0,\mathrm{min}}$ compared to particle-resolved simulations Stenhammar et al. (2013, 2014). Note that, within the approximations mentioned above, $(\mathrm{bar}{\rho}_{\mathrm{A,min}},v_{0,\mathrm{min}})$ is the critical point of the one-component system of active particles.

Next, we consider the spinodal condition (51) in the limit $\mathrm{bar}{v}\to\infty$ and assume that $a^{(\mathrm{AA})}_{0}$ and $a^{(\mathrm{AA})}_{1}$ 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 $\mathrm{bar}{\rho}_{\mathrm{A}}=1/(2a^{(\mathrm{AA})}_{0})$ Cates and Tailleur (2013).

IV Simulations

To confirm our analytical results from Sec. III, we carried out Brownian dynamics computer simulations of binary mixtures of isometric active and passive colloidal particles in two spatial dimensions. These particle-resolved simulations were based on solving the Langevin equations (1) and (2) numerically and made use of the LAMMPS molecular dynamics package Plimpton (1995). All particles in our simulations had the same diameter $\sigma=2R$, diffusion coefficients $D_{\mathrm{T}}$ and $D_{\mathrm{R}}=3D_{\mathrm{T}}/\sigma^{2}$, and repulsive interaction potential $U_{2}(r)$. For the latter, we used the Weeks-Chandler-Andersen potential, i.e., the truncated and shifted Lennard-Jones potential

[TABLE]

Here, $\varepsilon$ determines the interaction strength as well as the Lennard-Jones time scale $\tau_{\mathrm{LJ}}=\sigma^{2}/(\varepsilon\beta D_{\mathrm{T}})$. For convenience, we chose $\sigma$, $\tau_{\mathrm{LJ}}$, and $\varepsilon$ as units for length, time, and energy, respectively. The active particles had the low-density self-propulsion speed $v_{0}=\beta D_{\mathrm{T}}F_{\mathrm{A}}=24\sigma/\tau_{\mathrm{LJ}}$ with the active driving force $F_{\mathrm{A}}=24\varepsilon/\sigma$. We studied active-passive mixtures with various total average area fractions $\Phi=\Phi_{\mathrm{A}}+\Phi_{\mathrm{P}}$, fractions of particles that are active $\chi_{\mathrm{A}}=\Phi_{\mathrm{A}}/\Phi$, and Péclet numbers $\mathrm{Pe}=v_{0}\sigma/D_{\mathrm{T}}$, where $\Phi_{\mathrm{A}}=\mathrm{bar}{\rho}_{\mathrm{A}}\pi\sigma^{2}/4$ and $\Phi_{\mathrm{P}}=\mathrm{bar}{\rho}_{\mathrm{P}}\pi\sigma^{2}/4$ 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 $\mathrm{Pe}$ by changing $D_{\mathrm{T}}$ for all particles, while we kept $v_{0}$ constant (see Refs. Stenhammar et al. (2014, 2016) for details). In the limiting case $\chi_{\mathrm{A}}=0$, which corresponds to a system of only passive particles, we chose $D_{\mathrm{T}}=\sigma^{2}/\tau_{\mathrm{LJ}}$.

Our simulations were carried out for particles with homogeneously distributed random initial positions and orientations in a quadratic simulation domain with size $150\sigma\times 150\sigma$ and periodic boundary conditions. The number of particles was $5730$ for $\Phi=0.2$ and $22918$ for $\Phi=0.8$. Furthermore, the simulations comprised $5\cdot 10^{7}$ time steps of length $\Delta t=5\cdot 10^{-5}\tau_{\mathrm{LJ}}$, of which we discarded the first $2\cdot 10^{6}$ time steps. The duration of a simulation was thus $2500\tau_{\mathrm{LJ}}$. For $\mathrm{Pe}=100$, where the Brownian time $\tau_{\mathrm{B}}=\sigma^{2}/D_{\mathrm{T}}$ is $\tau_{\mathrm{B}}\approx 4.2\tau_{\mathrm{LJ}}$, this equals $600\tau_{\mathrm{B}}$. When the particles have diameter $\sigma=1\,\mathrm{\mathrm{\SIUnitSymbolMicro}\mathrm{m}}$ and translational diffusion coefficient $D_{\mathrm{T}}=k_{\mathrm{B}}T/(3\pi\eta\sigma)$ Einstein (1905), and are dispersed in water with dynamic viscosity $\eta=10^{-3}\,\mathrm{\mathrm{Pa}\,\mathrm{s}}$ and temperature $T=293\,\mathrm{\mathrm{K}}$, the Brownian time scale is $\tau_{\mathrm{B}}\approx 2.3\,\mathrm{\mathrm{s}}$ and the duration of the simulations corresponds to about $1400\,\mathrm{\mathrm{s}}$.

IV.1 Pair-correlation functions

Since the pair-correlation functions $g_{\mathrm{AA}}(r,\phi_{\mathrm{R}}-\phi)$, $g_{\mathrm{AP}}(r,\phi_{\mathrm{R}}-\phi)$, $g_{\mathrm{PA}}(r,\phi_{\mathrm{R}}-\phi)$, and $g_{\mathrm{PP}}(r,\phi_{\mathrm{R}}-\phi)$ are needed to evaluate the coefficients $a^{(\mathrm{AA})}_{0}$, $a^{(\mathrm{AP})}_{0}$, $a^{(\mathrm{AA})}_{1}$, $a^{(\mathrm{AP})}_{1}$, $a^{(\mathrm{PA})}_{1}$, and $a^{(\mathrm{PP})}_{1}$ via Eqs. (14), (18), and (19), we used the simulations described above to determine the pair-correlation functions for $\chi_{\mathrm{A}}=0.5$ and various parameter combinations of $\mathrm{Pe}\in[0,500]$ and $\Phi\in[0,0.8]$. However, we did not include parameter combinations where the isotropic state of the simulated system was not stable. The resulting pair-correlation functions for $\mathrm{Pe}=50$ and $\Phi=0.6$ are shown in Fig. 2.

For fixed $\phi_{\mathrm{R}}-\phi$, the shapes of these functions of $r$ are qualitatively similar to the typical shape of a radial distribution function for repulsive disks or spheres. However, the maxima and minima of $g_{\mathrm{AA}}(r,\phi_{\mathrm{R}}-\phi)$ and $g_{\mathrm{AP}}(r,\phi_{\mathrm{R}}-\phi)$ as functions of $r$ depend strongly on $\phi_{\mathrm{R}}-\phi$. They are most strongly pronounced for $\phi_{\mathrm{R}}-\phi=0$, corresponding to the active (reference) particle colliding with other particles right in front of it, and become increasingly less pronounced when $\lvert\phi_{\mathrm{R}}-\phi\rvert\leqslant\pi$ increases (see Figs. 2a and 2b). Furthermore, the positions of the maxima and minima move to slightly larger particle separations $r$ for growing $\lvert\phi_{\mathrm{R}}-\phi\rvert\leqslant\pi$. This is consistent with the bow wave formation in front of a moving active particle reported for one-component systems of active particles in Ref. Bialké et al. (2013). In contrast to $g_{\mathrm{AA}}(r,\phi_{\mathrm{R}}-\phi)$ and $g_{\mathrm{AP}}(r,\phi_{\mathrm{R}}-\phi)$, the pair-correlation functions $g_{\mathrm{PA}}(r,\phi_{\mathrm{R}}-\phi)$ and $g_{\mathrm{PP}}(r,\phi_{\mathrm{R}}-\phi)$ do not depend on $\phi_{\mathrm{R}}-\phi$ so that it is here sufficient to focus on the case $\phi_{\mathrm{R}}-\phi=0$.

When comparing the four pair-correlation functions for $\phi=\phi_{\mathrm{R}}$, it becomes obvious that $g_{\mathrm{AA}}(r,0)$ and $g_{\mathrm{AP}}(r,0)$, on the one hand, as well as $g_{\mathrm{PA}}(r,0)$ and $g_{\mathrm{PP}}(r,0)$, on the other hand, are very similar. The main difference between these two pairs of correlation functions is that the maxima of $g_{\mathrm{AA}}(r,0)$ and $g_{\mathrm{AP}}(r,0)$ are significantly larger and shifted to smaller $r$ compared to $g_{\mathrm{PA}}(r,0)$ and $g_{\mathrm{PP}}(r,0)$ (see Fig. 2c). This is a consequence of the different type of the reference particle, which is active for $g_{\mathrm{AA}}(r,0)$ and $g_{\mathrm{AP}}(r,0)$ but passive for $g_{\mathrm{PA}}(r,0)$ and $g_{\mathrm{PP}}(r,0)$. When the reference particle is active, it pushes other particle ahead of itself so that in front of it (i.e., at $\phi=\phi_{\mathrm{R}}$) on average over time the concentration of particles is larger and the particle separations are smaller than around a passive reference particle.

IV.2 Coefficients

Using Eq. (14) as well as our results for the pair-correlation functions for $\chi_{\mathrm{A}}=0.5$ and various choices of $\mathrm{Pe}$ and $\Phi$, we calculated the corresponding values of the coefficients $a^{(\mathrm{AA})}_{0}$, $a^{(\mathrm{AP})}_{0}$, $a^{(\mathrm{AA})}_{1}$, $a^{(\mathrm{AP})}_{1}$, $a^{(\mathrm{PA})}_{1}$, and $a^{(\mathrm{PP})}_{1}$. While for small $\mathrm{Pe}$ we observed a strong dependence of these coefficients on $\mathrm{Pe}$, we found that their dependence on $\mathrm{Pe}$ is negligible for $\mathrm{Pe}>50$. Under the condition $\mathrm{Pe}>50$, the calculated coefficients are approximately given by the expressions

[TABLE]

and

[TABLE]

Interestingly, the coefficient $a^{(\mathrm{AA})}_{0}$ is independent of the total particle area fraction $\Phi$ and $a^{(\mathrm{AP})}_{0}$ has only a moderate dependence on $\Phi$, whereas the other coefficients depend strongly on $\Phi$.

With Eqs. (56) and (57) the effective swim speed of an active particle (46) becomes $\mathrm{bar}{v}=v_{0}(1-0.86\mathrm{bar}{\rho}_{\mathrm{A}}-(1.04-0.67(\Phi-0.5)^{2})\mathrm{bar}{\rho}_{\mathrm{P}})$ or, equivalently, $\mathrm{bar}{v}=v_{0}(1-1.09\Phi_{\mathrm{A}}-(1.32-0.85(\Phi-0.5)^{2})\Phi_{\mathrm{P}})$. This is close666For $\Phi_{\mathrm{A}},\Phi_{\mathrm{P}}\geqslant 0$ and $\Phi_{\mathrm{A}}+\Phi_{\mathrm{P}}\leqslant\Phi_{\mathrm{cp}}$ with the close-packing area fraction $\Phi_{\mathrm{cp}}=\pi/(2\sqrt{3})$, the maximal difference of both expressions for $\mathrm{bar}{v}$ is less than $0.061v_{0}$. to the result $\mathrm{bar}{v}=v_{0}(1-1.08\Phi_{\mathrm{A}}-1.21\Phi_{\mathrm{P}})$ from Ref. Stenhammar et al. (2015), which has been obtained from particle-resolved simulations for different $\Phi_{\mathrm{A}}$ and $\Phi_{\mathrm{P}}$ by calculating $\mathrm{bar}{v}(\Phi_{\mathrm{A}},\Phi_{\mathrm{P}})$ directly as the time-averaged speed of an active particle. For $\Phi_{\mathrm{P}}=0$, our expression reduces to $\mathrm{bar}{v}=v_{0}(1-1.09\Phi_{\mathrm{A}})$, which is similar to the corresponding result $\mathrm{bar}{v}=v_{0}(1-1.05\Phi_{\mathrm{A}})$ from Refs. Stenhammar et al. (2013, 2014).

The equality of $a^{(\mathrm{AP})}_{1}$ and $a^{(\mathrm{PA})}_{1}$ in Eq. (60) is a consequence of using the same interaction potential for all particles and considering only systems with $\chi_{\mathrm{A}}=0.5$ when calculating the pair-correlation functions. These special conditions lead to $\int^{2\pi}_{0}\!\!\!\mathrm{d}\phi\,g_{\mathrm{AP}}(r,\phi)=\int^{2\pi}_{0}\!\!\!\mathrm{d}\phi\,g_{\mathrm{PA}}(r,\phi)$ and thus $a^{(\mathrm{AP})}_{1}=a^{(\mathrm{PA})}_{1}$ here, but they do not hold in general.

IV.3 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 $\mathrm{Pe}$ and the average area fraction of the active particles $\Phi_{\mathrm{A}}$ are sufficiently large, the homogeneous state becomes unstable and particle clusters form Stenhammar et al. (2015).

Figure 3 shows the dynamical state diagram of a one-component system of passive particles ($\chi_{\mathrm{A}}=0$), a binary mixture with the same number of active and passive particles ($\chi_{\mathrm{A}}=0.5$), and a one-component system of active particles ($\chi_{\mathrm{A}}=1$) in the $\mathrm{Pe}$-$\Phi$ plane, where the stable and unstable regions were determined by visual inspection.

As expected, the system is homogeneous for all $\mathrm{Pe}$ and $\Phi$ when there are no active particles (see Fig. 3a), whereas the state diagram of systems with half of the particles or all particles being active contains both stable and unstable regions (see Figs. 3b and 3c). When $\chi_{\mathrm{A}}>0$, Eq. (46) predicts $\mathrm{bar}{v}\leqslant 0$ for the largest particle packing fractions. The corresponding region in the state diagram is in good agreement with the region where the total particle area fraction $\Phi$ is near the close-packing area fraction $\Phi_{\mathrm{cp}}=\pi/(2\sqrt{3})\approx 0.907$ or even larger so that effective motion of the active particles is not possible.

In the case of an active-passive mixture with $\chi_{\mathrm{A}}=0.5$, clustering of the particles requires $\mathrm{Pe}>50$ and $\Phi>0.4$ (see Fig. 3b). Evidently, the predicted spinodal for the onset of a dynamical instability in the mixture, which is in general given by Eqs. (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 Eq. (47) describes the spinodal curve, whereas Eq. (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 ($\chi_{\mathrm{A}}=1$) is qualitatively similar to that of active-passive mixtures with $\chi_{\mathrm{A}}=0.5$ (see Fig. 3c). Also in the case of solely active particles, the occurrence of a motility-induced instability requires $\mathrm{Pe}>50$, but in this case an instability can already be observed for $\Phi$ down to $\Phi\approx 0.4$ Stenhammar et al. (2014). For $\chi_{\mathrm{A}}=1$, Eq. (47) has no solution and the spinodal curve is entirely described by Eq. (48) or, equivalently, Eq. (51). The predicted spinodal is in rather good agreement with that directly obtained from the simulations. Especially the minimal values of $\mathrm{Pe}$ required for the instability to occur, that are associated with the critical points of the spinodals, are very close together. When $a^{(\mathrm{AA})}_{1}=0$ is chosen in Eq. (51), the predicted spinodal includes too small values of $\mathrm{Pe}$ down to $4\sqrt{3}\approx 6.9$ and 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 $a^{(\mathrm{AA})}_{1}=a^{(\mathrm{AP})}_{1}=a^{(\mathrm{PA})}_{1}=a^{(\mathrm{PP})}_{1}=0$ in the case of a mixture with $\chi_{\mathrm{A}}=0.5$ leads to an even stronger deterioration of the predictions. Then the predicted spinodal not only includes regions of too small $\mathrm{Pe}$, but also corresponds to the wrong type of instability.

Remarkably, through Eqs. (47) and (48), which describe the spinodals for the active-passive mixtures with $\chi_{\mathrm{A}}=0.5$ and for the purely active systems with $\chi_{\mathrm{A}}=1$, 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, than that of purely active systems, where clusters only move through particle diffusion Stenhammar et al. (2015); Wysocki et al. (2016). Given the different nature of the instabilities for $\chi_{\mathrm{A}}=0.5$ and $\chi_{\mathrm{A}}=1$, 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. This could, at least partly, explain, why the agreement of the predicted spinodals with the simulation results is better for $\chi_{\mathrm{A}}=1$ than for $\chi_{\mathrm{A}}=0.5$.

V 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 self-propulsion 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 one-component 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.

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.

Acknowledgements.

We thank Edgar Knobloch, Hartmut Löwen, and Uwe Thiele for helpful discussions. R.W. is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – WI 4170/3-1. J.S. acknowledges funding by a Project Grant from the Swedish Research Council (2015-05449). M.E.C. holds a Royal Society Research Professorship.

Appendix A Possible further approximations of the pair-correlation functions

The number of independent coefficients $A^{(\mu\nu)}_{m,k}$ in Eqs. (12) and (13) can be reduced by taking into account the orientational dependence of the pair-correlation functions $g_{\mu\nu}(r,\phi)$ only up to second order, i.e., by approximating

[TABLE]

and by considering rather hard particles with very short-range interactions where the interaction lengths are much smaller than the particle radii $R_{\mu}$. With these approximations, the expressions for the coefficients (14) simplify to

[TABLE]

with

[TABLE]

When the particles are hard disks or spheres in a plane, the coefficients $\alpha^{(\mu\nu)}_{i}$ can be expressed in a simpler form. Following a procedure that is a standard trick when deriving an expression for the pressure of a passive hard-sphere fluid Hansen and McDonald (2009), we write $g^{(\mu\nu)}_{i}(r)$ in Eq. (64) as $\exp(-\beta U^{(\mu\nu)}_{2}(r))\exp(\beta U^{(\mu\nu)}_{2}(r))g^{(\mu\nu)}_{i}(r)$ and afterwards use $U^{(\mu\nu)\prime}_{2}(r)\exp(-\beta U^{(\mu\nu)}_{2}(r))=-\beta^{-1}\mathrm{d}\exp(-\beta U^{(\mu\nu)}_{2}(r))/\mathrm{d}r$. For hard-sphere interactions, the function $\exp(-\beta U^{(\mu\nu)}_{2}(r))$ becomes the Heaviside step function $\mathrm{H}(r-R_{\mu}-R_{\nu})$. Taking into account that $\mathrm{d}\mathrm{H}(r)/\mathrm{d}r=\delta(r)$, the coefficients $\alpha^{(\mu\nu)}_{i}$ can thus be written as

[TABLE]

with $g^{(\mu\nu)}_{i}(R_{\mu}+R_{\nu})^{+}=\lim_{\epsilon\to 0^{+}}g^{(\mu\nu)}_{i}(R_{\mu}+R_{\nu}+\epsilon)$ for hard particles. What now remains to determine are the values of $g^{(\mu\nu)}_{i}(R_{\mu}+R_{\nu})^{+}$. They can still depend on the speeds $v^{\mu}_{0}$ and the average particle number densities $\mathrm{bar}{\rho}_{\mu}$, but are not functions of $r$ or $\phi$. Furthermore, due to their close relationship to the pair-correlation functions $g_{\mu\nu}(r,\phi)$ it is possible to estimate at least lower and upper bounds for the values of $g^{(\mu\nu)}_{i}(R_{\mu}+R_{\nu})^{+}$. It should also be possible to make good estimates for their dependence on $v^{\mu}_{0}$ and $\mathrm{bar}{\rho}_{\mu}$.

Appendix B Details on the derivation of Eq. (15)

Following the common notation in Ginzburg-Landau theories for liquid crystals de Gennes and Prost (1995); Wittkowski et al. (2010, 2011a, 2011b), it is useful to approximate the orientation-dependent one-particle densities

[TABLE]

in the dynamic equations (10) by the order-parameter fields

[TABLE]

These order-parameter fields are the translational densities $\rho_{\mu}(\vec{r},t)$, the local polarizations $\vec{\mathcal{P}}_{\mu}(\vec{r},t)$, and the symmetric and traceless nematic tensors $\boldsymbol{\mathcal{Q}}_{\mu}(\vec{r},t)$. The corresponding dynamic equations for the order-parameter fields are777In the following, the Einstein summation convention is used for Latin indices. There is, however, no implicit summation over Greek indices.

[TABLE]

with the currents $\vec{J}^{\rho_{\mu}}(\vec{r},t)$ and the quasi-currents888Notice that the translational densities $\rho_{\mu}(\vec{r},t)$ are conserved quantities, whereas the local polarizations $\vec{\mathcal{P}}_{\mu}(\vec{r},t)$ and the nematic tensors $\boldsymbol{\mathcal{Q}}_{\mu}(\vec{r},t)$ are not conserved. $\vec{\Phi}^{\mathcal{P}_{\mu}}(\vec{r},t)$ and $\boldsymbol{\Phi}^{\mathcal{Q}_{\mu}}(\vec{r},t)$. Up to second order in gradients, these currents and quasi-currents are given by

[TABLE]

with the density-dependent swim speeds

[TABLE]

the tensor $\chi_{ijkl}=\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}-\delta_{ij}\delta_{kl}$, and $\partial_{i}$ being the $i$th element of the del symbol $\vec{\nabla}_{\vec{r}}=(\partial_{1},\partial_{2})=(\partial_{x},\partial_{y})$. When considering the limiting case of a one-component system ($n=1$) and setting $A^{(\mu\nu)}_{1,0}=A^{(\mu\nu)}_{1,2}=0$, which means neglecting all effects of the particle interactions except for the slow-down of active particles due to collisions that is taken into account through the density-dependent effective swim speeds (76), our Eqs. (70)-(75) reduce to the corresponding dynamic equations for active Brownian particles in two spatial dimensions of Ref. Cates and Tailleur (2013).999Setting $n=1$ and $A^{(\mu\nu)}_{1,0}=A^{(\mu\nu)}_{1,2}=0$ in our Eqs. (70)-(75) results in dynamic equations that are equivalent to those that are obtained by setting $\varphi=\rho/(2\pi)$, $d=2$, $\alpha=0$, $D_{t}=\text{const.}$, $D_{r}=\text{const.}$, and $\chi_{abc}=0$ in Eqs. (4)-(6) of Ref. Cates and Tailleur (2013).

Since only the translational densities $\rho_{\mu}(\vec{r},t)$ are conserved quantities, their relaxation times are much larger than the relaxation times of $\vec{\mathcal{P}}_{\mu}(\vec{r},t)$ and $\boldsymbol{\mathcal{Q}}_{\mu}(\vec{r},t)$. When the nematic tensors $\boldsymbol{\mathcal{Q}}_{\mu}(\vec{r},t)$ are the fastest-relaxing order-parameter fields, we can describe the system on a larger time scale where $\dot{\mathcal{Q}}_{\mu,ij}=0$. Equation (75) can then be used to express $\boldsymbol{\mathcal{Q}}_{\mu}(\vec{r},t)$ in terms of gradients of the order-parameter fields. Inserting this expression recursively into Eqs. (73)-(75) and still neglecting all terms of third or higher order in gradients results in the quasi-stationary approximation

[TABLE]

which involves dynamic equations for $\rho_{\mu}(\vec{r},t)$ and $\vec{\mathcal{P}}_{\mu}(\vec{r},t)$ and a constitutive equation that expresses $\boldsymbol{\mathcal{Q}}_{\mu}(\vec{r},t)$ in terms of $\rho_{\mu}(\vec{r},t)$ and $\vec{\mathcal{P}}_{\mu}(\vec{r},t)$. On the even larger time scale of the relaxation times of the translational densities $\rho_{\mu}(\vec{r},t)$, we can use $\dot{\mathcal{P}}_{\mu,i}=0$ to further simplify our dynamic equations (77)-(79). Proceeding with $\vec{\mathcal{P}}_{\mu}(\vec{r},t)$ in a similar way as with $\boldsymbol{\mathcal{Q}}_{\mu}(\vec{r},t)$ leads to the full quasi-stationary approximation

[TABLE]

which involves only a dynamic equation for $\rho_{\mu}(\vec{r},t)$ and constitutive equations for the other order-parameter fields.

The first term on the right-hand side of Eq. (80) is proportional to the local polarizations (81). This shows that gradients of $v_{\mu}(\{\rho_{\mu}\})\rho_{\mu}(\vec{r},t)$ lead to local polarizations that drive currents of the translational densities. The other two terms on the right-hand side of Eq. (80) describe diffusive currents, where the first one represents the translational diffusion of noninteracting particles and the second one takes interactions with particles of the same and other species into account. As can be seen from the first term on the right-hand side of Eq. (82), gradients of $v_{\mu}(\{\rho_{\mu}\})\rho_{\mu}(\vec{r},t)$ also cause local nematic order. Like the local polarizations described by Eq. (81), the first term on the right-hand side of Eq. (82) vanishes for passive particles. In contrast, the second term on the right-hand side of Eq. (82) does not in general vanish for passive particles. When we set $n=1$ and $A^{(\mu\nu)}_{1,0}=A^{(\mu\nu)}_{1,2}=0$, our Eqs. (80)-(82) corresponds to Eqs. (8) and (9) and the unnumbered equation before Eq. (8) in Ref. Cates and Tailleur (2013). Furthermore, in this limit our Eq. (80) reduces to the simpler diffusion equation that is given by Eq. (20) in Ref. Bialké et al. (2013).

Appendix C Stability criterion for mixtures

From Eq. (24) it follows that an $n$-component mixture of colloidal particles with different motilities is stable, if

[TABLE]

where $\lVert\,\cdot\,\rVert_{\mathrm{F}}$ is the Frobenius norm, and unstable otherwise. The constant diffusion tensor $\boldsymbol{\mathrm{bar}{\mathrm{D}}}=\boldsymbol{\mathrm{bar}{\mathrm{D}}}_{0}+\boldsymbol{\mathrm{bar}{\mathrm{D}}}_{1}$ can be uniquely decomposed into a diagonalizable matrix $\boldsymbol{\mathrm{bar}{\mathrm{D}}}_{0}$ with the same eigenvalues as $\boldsymbol{\mathrm{bar}{\mathrm{D}}}$ and a nilpotent matrix $\boldsymbol{\mathrm{bar}{\mathrm{D}}}_{1}$ that commutes with $\boldsymbol{\mathrm{bar}{\mathrm{D}}}_{0}$ (Jordan-Chevalley decomposition). Equation (83) can therefore be written as

[TABLE]

Since $\boldsymbol{\mathrm{bar}{\mathrm{D}}}_{1}$ is nilpotent, the matrix exponential $\exp(\!-\boldsymbol{\mathrm{bar}{\mathrm{D}}}_{1}t)$ is a finite-order polynomial in $-\boldsymbol{\mathrm{bar}{\mathrm{D}}}_{1}t$. With special exceptions, which constitute a null set among the set of all possible realizations of $\boldsymbol{\mathrm{bar}{\mathrm{D}}}$, $\exp(\!-\boldsymbol{\mathrm{bar}{\mathrm{D}}}_{1}t)$ can therefore be neglected in Eq. (84).

Furthermore, we can take advantage of the fact that the Frobenius norm is invariant under unitary transformations and replace $\boldsymbol{\mathrm{bar}{\mathrm{D}}}_{0}$ by a corresponding diagonal matrix $\boldsymbol{\Lambda}=\operatorname{diag}(\lambda_{1},\dotsc,\lambda_{n})$, where $\lambda_{1},\dotsc,\lambda_{n}$ are the eigenvalues of $\boldsymbol{\mathrm{bar}{\mathrm{D}}}_{0}$ and thus of $\boldsymbol{\mathrm{bar}{\mathrm{D}}}$. The stability criterion (84) then becomes

[TABLE]

This means that the mixture is stable, if the real parts of all eigenvalues of $\boldsymbol{\mathrm{bar}{\mathrm{D}}}$ are nonnegative, i.e.,

[TABLE]

When one of the eigenvalues of $\boldsymbol{\mathrm{bar}{\mathrm{D}}}$ has a negative real part, the mixture is unstable.