Multi-species dynamical density functional theory for microswimmers: derivation, orientational ordering, trapping potentials, and shear cells

TL;DR

This paper develops a dynamical density functional theory for multi-species microswimmer systems, accounting for external potentials and interactions, and analyzes their collective behaviors, stability, and responses in various confined and shear flow environments.

Contribution

It extends existing theories to multi-species microswimmers, incorporating hydrodynamic interactions and external potentials, and applies the framework to analyze collective phenomena and stability.

Findings

Majority species influence minority behavior in traps.

Thresholds for puller fraction and propulsion strength induce collective motion.

Confined active swimmers exhibit circle-like trajectories under shear flows.

Abstract

Microswimmers typically operate in complex environments. In biological systems, often diverse species are simultaneously present and interact with each other. Here, we derive a (time-dependent) particle-scale statistical description, namely a dynamical density functional theory, for such multi-species systems, extending existing works on one-component microswimmer suspensions. In particular, our theory incorporates the effect of external potentials, but also steric and hydrodynamic interactions between swimmers. For the latter, a previously introduced force-dipole-based minimal (pusher or puller) microswimmer model is used. As a limiting case of our theory, mixtures of hydrodynamically interacting active and passive particles are captured as well. After deriving the theory, we apply it to different planar swimmer configurations. First, these are binary pusher--puller mixtures in external traps. In the considered situations, we find that the majority species imposes its behavior on the minority species. Second, for unconfined binary pusher--puller mixtures, the linear stability of an orientationally disordered state against the emergence of global polar orientational order (and thus emergent collective motion) is tested analytically. Our statistical approach predicts, qualitatively in line with previous particle-based computer simulations, a threshold for the fraction of pullers and for their propulsion strength that lets overall collective motion arise. Third, we let driven passive colloidal particles form the boundaries of a shear cell, with confined active microswimmers on their inside. Driving the passive particles then effectively imposes shear flows, which persistently acts on the inside microswimmers. Their resulting behavior reminds of the one of circle swimmers, though with varying swimming radii.

Figures

Equations

$${\mathbf{f}}_{i\pm }=\pm f_i{\hat{\mathbf{n}}}_i$$

$${\mathbf{r}}_{i+}=$$

$${\mathbf{r}}_{i-}=$$

$$\displaystyle\left[\begin{array}[]{c}{\mathbf{v}}_{i}\\[2.84544pt] {\bm{\omega}}_{i}\end{array}\right]=\sum_{j=1}^{N}\Bigg{(}\left[\begin{array}[]{cc}{\bm{\mu}}^{\mathrm{tt}}_{ij}&{\bm{\mu}}^{\mathrm{tr}}_{ij}\\[2.84544pt] {\bm{\mu}}^{\mathrm{rt}}_{ij}&{\bm{\mu}}^{\mathrm{rr}}_{ij}\end{array}\right]$$

$${\mathbf{\mu }}_{ij}^{tt}$$

$$\displaystyle{}+\frac{a^{2}_{i}+a^{2}_{j}}{4}\Big{(}\frac{1}{r_{ij}}\Big{)}^{3}\Big{(}{\bf{1}}-3\,{{\mathbf{\hat{r}}_{ij}\mathbf{\hat{r}}_{ij}}}\Big{)}\bigg{)},$$

$${\mathbf{\mu }}_{ij}^{rr}$$

$${\mathbf{\mu }}_{ij}^{tr}$$

$${\mathbf{\mu }}_{ii}^{\mathrm{tt}}={\mu }_i^{\mathrm{t}}1,{\mathbf{\mu }}_{ii}^{\mathrm{rr}}={\mu }_i^{\mathrm{r}}1,{\mathbf{\mu }}_{ii}^{\mathrm{tr}}={\mathbf{\mu }}_{ii}^{\mathrm{rt}}=0,$$

$${\mu }_i^{\mathrm{t}}=1/(6\pi \eta a_i),{\mu }_i^{\mathrm{r}}=1/(8\pi \eta a_i^3).$$

$${\mathbf{\Lambda }}_{ij}^{\mathrm{tt}}$$

$${\mathbf{\Lambda }}_{ij}^{\mathrm{rt}}$$

$${\mathbf{\mu }}_{ij}^{\mathrm{tt}\pm }$$

$${\mathbf{\mu }}_{ij}^{\mathrm{rt}\pm }$$

$${\mathbf{r}}_{ij}^{+}$$

$${\mathbf{r}}_{ij}^{-}$$

$${\mathbf{v}}_{0i}$$

$$U({\mathbf{r}}_1,\dots ,{\mathbf{r}}_N)=\sum _{k=1}^N{u}_{\mathrm{ext}}^k({\mathbf{r}}_k)+\sum _{k,l=1;k<l}^N{u}^{kl}(|{\mathbf{r}}_k-{\mathbf{r}}_l|).$$

$$u^{kl}(|{\mathbf{r}}_k-{\mathbf{r}}_l|)={ϵ}_0^{kl}\exp \left(-{\left(\frac{|{\mathbf{r}}_k-{\mathbf{r}}_l|}{{\sigma }_{kl}}\right)}^4\right),$$

$${\mathbf{F}}_j=-k_{\mathrm{B}}T{\nabla }_{{\mathbf{r}}_j}\ln P-{\nabla }_{{\mathbf{r}}_j}U({\mathbf{r}}_1,\dots ,{\mathbf{r}}_N),$$

$${\mathbf{T}}_j=-k_{\mathrm{B}}T{\hat{\mathbf{n}}}_j\times {\nabla }_{{\hat{\mathbf{n}}}_j}\ln P.$$

$$\frac{\partial P}{\partial t}={}-\sum_{i=1}^{N}\Big{(}\nabla_{\!{\mathbf{r}}_{i}}\cdot\left(\mathbf{v}_{i}P\right)+\left({\mathbf{\hat{n}}}_{i}\times\nabla_{\!{\mathbf{\hat{n}}}_{i}}\right)\cdot\left(\bm{\omega}_{i}P\right)\Big{)}$$

$${\rho }^{(m,n)}({\mathbf{X}}_{\mathrm{A}}^m,{\mathbf{X}}_{\mathrm{B}}^n,t)=$$

$$\frac{\partial {\rho }_{\mathrm{A}}(\mathbf{X},t)}{\partial t}=$$

$$\displaystyle+\!{\bm{\mathcal{K}}}_{\!\mathrm{AA}}^{\mathrm{tr}}\!+\!{\bm{\mathcal{K}}}_{\!\mathrm{AA}}^{\mathrm{ta}}\!+\!{\bm{\mathcal{K}}}_{\!\mathrm{AB}}^{\mathrm{tt}}\!+\!{\bm{\mathcal{K}}}_{\!\mathrm{AB}}^{\mathrm{tr}}\!+\!{\bm{\mathcal{K}}}_{\!\mathrm{AB}}^{\mathrm{ta}}\Big{)}$$

$$\displaystyle-(\mathbf{\hat{n}}\times\nabla_{\mathbf{\hat{n}}})\cdot\Big{(}{\bm{\mathcal{J}}}_{\!\!\mathrm{A}}^{\mathrm{rt}}\!+\!{\bm{\mathcal{J}}}_{\!\!\mathrm{A}}^{\mathrm{rr}}\!+\!{\bm{\mathcal{J}}}_{\!\!\mathrm{A}}^{\mathrm{ra}}\!+\!{\bm{\mathcal{K}}}_{\!\mathrm{AA}}^{\mathrm{rt}}$$

$$\displaystyle+\!{\bm{\mathcal{K}}}_{\!\mathrm{AA}}^{\mathrm{rr}}\!+\!{\bm{\mathcal{K}}}_{\!\mathrm{AA}}^{\mathrm{ra}}\!+\!{\bm{\mathcal{K}}}_{\!\mathrm{AB}}^{\mathrm{rt}}\!+\!{\bm{\mathcal{K}}}_{\!\mathrm{AB}}^{\mathrm{rr}}\!+\!{\bm{\mathcal{K}}}_{\!\mathrm{AB}}^{\mathrm{ra}}\Big{)}.$$

$${\mathcal{J}}_{\mathrm{A}}^{\mathrm{tt}}=$$

$$+\int \mathrm{d}{\mathbf{X}}^{\prime }{\rho }^{(2,0)}(\mathbf{X},{\mathbf{X}}^{\prime },t){\nabla }_{\mathbf{r}}{u}^{\mathrm{AA}}(|\mathbf{r}-{\mathbf{r}}^{\prime }|)$$

$$\displaystyle+\!\int\!\!\mathrm{d}{\mathbf{X}}^{\prime}\rho^{(1,1)}({\mathbf{X}},{\mathbf{X}}^{\prime},t)\nabla_{{\mathbf{r}}}u^{\mathrm{AB}}({|{\mathbf{r}}-{\mathbf{r}}^{\prime}|})\Big{)},$$

$${\mathcal{J}}_{\mathrm{A}}^{\mathrm{ta}}=$$

$${\mathcal{J}}_{\mathrm{A}}^{\mathrm{rr}}=$$

$${\mathcal{J}}_{\mathrm{A}}^{\mathrm{tr}}=$$

$${\mathcal{K}}_{\mathrm{AA}}^{\mathrm{tt}}=$$

$$+{\rho }^{(2,0)}(\mathbf{X},{\mathbf{X}}^{\prime },t){\nabla }_{{\mathbf{r}}^{\prime }}\left(u_{\mathrm{ext}}^{\mathrm{A}}({\mathbf{r}}^{\prime })+u^{\mathrm{AA}}(|\mathbf{r}-{\mathbf{r}}^{\prime }|)\right)$$

$$+\int \mathrm{d}{\mathbf{X}}^{\prime \prime }{\rho }^{(2,1)}(\mathbf{X},{\mathbf{X}}^{\prime },{\mathbf{X}}^{\prime \prime },t){\nabla }_{{\mathbf{r}}^{\prime }}{u}^{\mathrm{AB}}(|{\mathbf{r}}^{\prime }-{\mathbf{r}}^{\prime \prime }|)$$

$$\displaystyle+\!\int\!\!\mathrm{d}{\mathbf{X}}^{\prime\prime}\rho^{(3,0)}({\mathbf{X}},{\mathbf{X}}^{\prime},{\mathbf{X}}^{\prime\prime},t)\nabla_{{\mathbf{r}}^{\prime}}u^{\mathrm{AA}}({|{\mathbf{r}}^{\prime}-{\mathbf{r}}^{\prime\prime}|})\bigg{)},$$

$${\mathcal{K}}_{\mathrm{AA}}^{\mathrm{tr}}=$$

$$=$$

$${\mathcal{K}}_{\mathrm{AA}}^{\mathrm{ta}}=$$

$${\mathcal{K}}_{\mathrm{AA}}^{\mathrm{rt}}=$$

Full text

I Introduction

From a fundamental point of view, the study of active microswimmers Purcell (1977); Lauga and Powers (2009); Elgeti et al. (2015); Menzel (2015a); Zöttl and Stark (2016); Bechinger et al. (2016) — i.e., micron-sized self-propelling particles suspended in a fluid — is interesting already because of the inherent non-equilibrium nature of self-propelling particles.Ramaswamy (2010); Cates (2012); Marchetti et al. (2013); Menzel (2016) Unusual collective behavior arises from this feature, e.g., motility-induced phase separation (MIPS)Cates and Tailleur (2013); Buttinoni et al. (2013); Bialké et al. (2013, 2015); Cates and Tailleur (2015); Solon et al. (2018); Digregorio et al. (2018) and laning.Wensink and Löwen (2012); Menzel and Ohta (2012); Menzel (2013); Kogler and Klapp (2015); Romanczuk et al. (2016); Menzel (2016) Moreover, on the applied side, natural biological microswimmers Purcell (1977); Kessler (1985); Eisenbach and Giojalas (2006); Engstler et al. (2007); Berg (2008); Polin et al. (2009); Durham et al. (2009); Mussler et al. (2013); Goldstein (2015) occur in almost all locations on Earth, including the human body, and artificial microswimmers Paxton et al. (2004); Howse et al. (2007); Buttinoni et al. (2012); Walther and Müller (2013); Samin and van Roij (2015); Moran and Posner (2017) may in the near future be used in medical and technical applications on the microscale, e.g., for precise drug delivery,Wang and Gao (2012); Patra et al. (2013); Abdelmohsen et al. (2014); Ma et al. (2015); Demirörs et al. (2018) for non-invasive surgery,Nelson et al. (2010); Xi et al. (2013); Abdelmohsen et al. (2014) when guiding sperm cells, Magdanz et al. (2013) and to power microengines.Di Leonardo et al. (2010); Sokolov et al. (2010); Maggi et al. (2016)

Both biological and artificial microswimmers typically operate under complex conditions. Bechinger et al. (2016) For example, the complexity can arise from steric confinement of the swimmers Wensink and Löwen (2008); Elgeti and Gompper (2009); Kaiser et al. (2012); Licata et al. (2016); Janssen et al. (2017); Praetorius et al. (2018) or be induced by a complex dispersion medium.Gomez-Solano et al. (2016); Narinder et al. (2018); Liebchen et al. (2018); Ferreiro-Córdova et al. (2018); Daddi-Moussa-Ider and Menzel (2018); Puljiz and Menzel (2019) Here, we consider the complementing case of complexity caused by interactions between different swimmer species, as can occur in a diverse set of situations.

In medical contexts, active multi-species systems (including both active–active and active–passive mixtures) appear when active agents, e.g., pathogenic bacteria or cargo-delivering microrobots, interact with (similar-sized) human cells. Furthermore, real-world microorganisms can change between motile and non-motile (i.e., active and passive in our notation) behavior during their life, with the organization in many-particle biofilms Flemming et al. (2016); Vidakovic et al. (2018) and active carpets Berke et al. (2008); Mathijssen et al. (2018) as examples for extreme cases. Also, different mutant lines of the same bacterial species can show different motility properties, see, e.g., motile and non-motile strains of E. coli bacteria.Berg (2008) More in general, subgroups of swimmers may be identified, if a strong polydispersity, e.g., of swimming speeds, is present inside a system. Finally, at least two species of swimmers are necessary to construct “heteronuclear” (i.e., composed of different building blocks) microswimmer molecules.Babel et al. (2016); Küchler et al. (2016); Löwen (2018)

Despite these manifold possible applications, studies on mixtures of microswimmers (and active particles in general) are still relatively rare. The problems regarded thus far include predator–prey dynamics,Mecholsky et al. (2010); Sengupta et al. (2011) mixtures of active rotors with opposite senses of rotationNguyen et al. (2014); Yeo et al. (2015) (see also the corresponding macroscale equivalent in Ref. Scholz et al., 2018), transport of passive V-shaped cargo particles by active rods in the bulk Angelani and Di Leonardo (2010); Kaiser et al. (2014, 2015a, 2015b) and by circle swimmers in channels Liao et al. (2018), depletion interactions between passive particles induced by an active bath,Ray et al. (2014); Ni et al. (2015) segregation effects in mixtures of Taylor-line swimmers propelling by self-deformation,Agrawal and Babu (2018) mixtures in which the activity is introduced by an effective colored noise,Wittmann et al. (2018) mesoscale transport phenomena in multi-species microorganism systems,Ben-Jacob et al. (2016) and MIPS-like phase separation in active–passive mixtures.Stenhammar et al. (2015); Kümmel et al. (2015); Wysocki et al. (2016); Wittkowski et al. (2017); Alaimo and Voigt (2018) Furthermore, collective behavior in mixtures of straight-propelling particlesMenzel (2012); Guisandez et al. (2017) and those migrating on circular trajectoriesLevis et al. (2018); Levis and Liebchen (2019) has been studied assuming Vicsek-type Vicsek et al. (1995); Toner and Tu (1995); Ihle (2011) effective alignment interactions between the swimmers. In addition to that, particle-based computer simulations of binary mixtures of microswimmers with different types of propulsion mechanisms, subject to mutual hydrodynamic interactions, have been performed to quantify the effect on the overall collective alignment behavior.Pessot et al. (2018)

In the present work, we cover multi-component microswimmer suspensions subject to external potentials. Different species here are mutually interacting, both via steric interaction potentials and via (far-field) hydrodynamic interactions. The latter may, following classical statistical mechanics (for passive particles), affect the dynamic behavior but, in general, not the appearance of static equilibrium systems. Microswimmer suspensions, however, are inherently out of equilibrium so that even steady states may be significantly altered by hydrodynamic interactions, calling for their incorporation in the physical description. Additionally, interesting phenomena can appear when hydrodynamic effects interplay with, e.g., magnetic interactions.Babel et al. (2016); Guzmán-Lastra et al. (2016)

Generally, supplementing experiments and many-body particle-based simulations with statistical descriptions, e.g., density-field equations, allows for thorough theoretical analysis. Ideally, the observed phenomena are explained in this way and new types of behavior are predicted, leading to a better understanding of the underlying physical effects. A well-established way of finding such density-field equations in non-equilibrium colloidal systems is dynamical density functional theory (DDFT).Evans (1979, 1992); Marconi and Tarazona (1999, 2000); Archer and Evans (2004); Chan and Finken (2005); Wensink and Löwen (2008); Español and Löwen (2009); Evans (2010); Löwen (2010); Wittkowski and Löwen (2011); Menzel et al. (2016); Hoell et al. (2017, 2018) Accordingly, we successfully derived a DDFT for one-species microswimmer systems and applied it to several example situations in previous works.Menzel et al. (2016); Hoell et al. (2017, 2018) In other contexts, DDFTs for mixtures of passive colloidal particles have been developed before.Archer (2005); Goddard et al. (2013); Härtel et al. (2015); Babel et al. (2018); Somerville et al. (2018) Here, we combine these two approaches and explicitly allow for different species of active microswimmers (and / or passive particles). In addition to the applications listed above, such a DDFT might in the future help to find dynamic correlation functions in one-component systems via “test-particle” methods.Archer et al. (2007); Hopkins et al. (2010); Hoell et al. (2018) We remark that multi-species DDFT approaches can also be used to describe the dynamics of other kinds of active matter, e.g., the growth of tumors in cell tissues.Al-Saedi et al. (2018)

Below, the employed microswimmer model — introduced in previous works Menzel et al. (2016); Hoell et al. (2017); Pessot et al. (2018); Hoell et al. (2018) — and its implications for hydrodynamic interactions are overviewed in Sec. II. It is then used in Sec. III as an input to derive the statistical theory, namely the multi-species dynamical density functional theory for microswimmers. Subsequently, several applications of the theory are discussed in Sec. IV, where we confine ourselves to planar arrangements within three-dimensional fluids. First, extending the one-component case analyzed previously,Menzel et al. (2016); Hoell et al. (2017) we explore binary mixtures of microswimmers in an external trap and find additional steady states resulting from interspecies interactions. Second, the possibility of emergent overall orientational order due to hydrodynamic interactions in binary mixtures of microswimmers is discussed. Third, microswimmers confined inside an externally driven ring of passive colloidal particles are investigated. The passive particles induce a shear flow that the enclosed active swimmers are exposed to. Finally, a short summary and an outlook are given in Sec. V.

II Swimmer model and the resulting hydrodynamic interactions

Before a particle-scale statistical description can be developed in Sec. III, an appropriate discretized swimmer model must first be defined. In particular, the hydrodynamic interactions between individual swimmers are specified below. For this purpose, we briefly review the previously-introduced minimal swimmer model.Menzel et al. (2016); Hoell et al. (2017); Pessot et al. (2018); Hoell et al. (2018)

Since a microswimmer cannot exert a net force on the surrounding liquid,Purcell (1977); ten Hagen et al. (2015) the far-field fluid flow around a swimmer (to lowest order) can typically be described as if it were caused by a force dipole acting on the fluid. (Exceptions are “neutral-type” swimmers with a vanishing time-averaged force-dipole contribution,Najafi and Golestanian (2004); Zargar et al. (2009); Daddi-Moussa-Ider et al. (2018a, b) which only feature higher-order multipole terms in the far-field flow caused, e.g., by an effective force quadrupole.) Here, we explicitly resolve the force dipole by two oppositely-oriented forces of equal magnitude.

Depending on whether the forces push out or pull in the fluid along the axis of self-propulsion, one distinguishes between pusher (extensile) and puller (contractile) microswimmers.Underhill et al. (2008); Baskaran and Marchetti (2009) Consequently, pushers draw in the fluid from the transverse directions, while pullers expel it along them. Our model can cover both cases, as detailed below.

Low Reynolds numbers — as are typical for microswimmers Purcell (1977) — and incompressibility of the fluid are henceforth assumed. Particularly, this means that the response of the fluid to a force is linear, overdamped, and instantaneous. In the bulk, the analytically-known Oseen tensor then explicitly connects the effect of a point-like force center to the resulting fluid flow.Happel and Brenner (2012); Dhont (1996); Kim and Karrila (2013) For finite-sized spherical particles subject to net forces and torques, the way to find (approximate) expressions for the induced hydrodynamic interactions between them is well-established.Dhont (1996); Kim and Karrila (2013)

This said, we now detail our minimal microswimmer model, see Fig. 1, referring to one swimmer labeled by $i$. In this model, a no-slip boundary encloses the spherical swimmer body, the latter being centered at position ${\mathbf{r}}_{i}$ and being of radius $a_{i}$. Below, ${\mathbf{v}}_{i}$ and ${\bm{\omega}}_{i}$ denote the velocity and angular velocity of the sphere, respectively.

Additionally, two oppositely-oriented forces

[TABLE]

of equal magnitude are exerted by the swimmer onto the surrounding fluid at positions

[TABLE]

respectively, relative to its body center. They move and rotate along with the sphere and create the fluid flow that (self-)propels the swimmer. Here, ${\mathbf{\hat{n}}}_{i}$ is the unit vector describing the orientation of the axially symmetric swimmer, $L_{i}>2a_{i}$ is the distance between the two force centers, and $|f_{i}|$ sets the magnitude of the forces. Depending on the sign of $f_{i}$, either pusher ($f_{i}>0$) or puller ($f_{i}<0$) microswimmers are constructed. Furthermore, the real number $\alpha_{i}$ (with $a_{i}/L_{i}<\alpha_{i}<1/2$) is a geometric parameter, see Fig. 1, that quantifies the breaking of the front–rear symmetry, which implies self-propulsion. The swimmer self-propels in the direction of ${\mathbf{\hat{n}}}_{i}$ for pushers, see Fig. 1(a), and ${}-{\mathbf{\hat{n}}}_{i}$ for pullers, see Fig. 1(b).

Moreover, an isotropic steric interaction between the swimmers is assumed that avoids unphysical overlap between force centers and bodies of different swimmers. As indicated in Fig. 1 and further detailed later, the effective center-to-center range of the steric interactions is denoted by $\sigma_{i}$.

Next, we specify the hydrodynamic interactions in a system of $N$ potentially different such model swimmers, labeled by $i=1,\dots,N$. For shorter notation, we furthermore define the phase-space coordinate ${\mathbf{X}}_{i}=\left\{{\mathbf{r}}_{i},{\mathbf{\hat{n}}}_{i}\right\}$ of each swimmer $i$. In our overdamped system of microswimmers in suspension,  ${\mathbf{v}}_{i}$ and ${\bm{\omega}}_{i}$ follow instantaneously from the microstate ${\mathbf{X}}^{N}=\left\{{\mathbf{X}}_{1},\dots,{\mathbf{X}}_{N}\right\}$.

In principle, hydrodynamic interactions are many-body interactions.Happel and Brenner (2012); Dhont (1996); Kim and Karrila (2013) Yet, already the lowest-order contributions beyond pairwise interactions are of fourth order in the ratio of body size to swimmer distanceDhont (1996) and can be neglected when one is primarily interested in the effect of far-field hydrodynamic interactions, e.g., in semi-dilute suspensions.Beenakker and Mazur (1983, 1984); Nägele (1996); Nägele and Baur (1997); Banchio et al. (2006) This is further supported by our use of repulsive steric interactions between swimmers, as detailed below, that keep them at distances from each other that are significantly larger than their hydrodynamic radii $a_{i}$, see also Fig. 1. Thus, here we only account for pairwise interactions and restrict ourselves to an expansion up to (including) the third order, also known as the Rotne-Prager level.Rotne and Prager (1969); Reichert and Stark (2004)

Following this idea, ${\mathbf{v}}_{i}$ and ${\bm{\omega}}_{i}$ are connected to the (non-hydrodynamic) forces ${\mathbf{F}}_{j}$ and torques ${\mathbf{T}}_{j}$ acting on the swimmer bodies $j=1,\dots,N$ and the self-propulsion forces that the swimmers exert via Menzel et al. (2016)

[TABLE]

Here, the mobility tensors representing the passive hydrodynamic interactions between two swimmer bodies $i\neq j$ are given by Dhont (1996); Rotne and Prager (1969); Reichert and Stark (2004); Menzel et al. (2016)

[TABLE]

where $\eta$ is the dynamic viscosity of the fluid, “$\times$” denotes the outer vector product, $\bm{1}$ represents the identity matrix, ${\mathbf{r}}_{ij}={\mathbf{r}}_{j}-{\mathbf{r}}_{i}$ is the distance vector, $r_{ij}=|{\mathbf{r}}_{ij}|$ denotes its absolute value, and ${\mathbf{\hat{r}}}_{ij}={\mathbf{r}}_{ij}/r_{ij}$ is the corresponding unit vector. Additionally, the passive “self” (i.e., $i=j$) mobilities read (no summation over repeated indices in these expressions)

[TABLE]

with

[TABLE]

Next, the active contribution to Eq. (14) is given by the tensors Menzel et al. (2016)

[TABLE]

with

[TABLE]

and

[TABLE]

As can be seen, there is only little change to the one-species case ($a_{i}=a_{j}\equiv a$) Menzel et al. (2016) at this order of the expansion in $a_{k}/r_{ij}$, $k=i,j$, namely only in Eq. (15).

Setting $i=j$ in Eqs. (20) and (21), the velocity and angular velocity of a free swimmer $i$ are obtained as Pessot et al. (2018)

[TABLE]

and, respectively, ${\bm{\omega}}_{0i}={\mathbf{0}}$. Thus, in the absence of thermal noise and outer influences, this kind of swimmer self-propels along a straight trajectory. Corresponding circle swimmers of axial asymmetry and a non-vanishing ${\bm{\omega}}_{0i}$ were considered in a previous work.Hoell et al. (2017)

We remark that, for simplicity, we here do not account for the distortions caused by the finite spherical swimmer bodies on the flow field induced by the active force centers.Kim and Karrila (2013); Adhyapak and Jabbari-Farouji (2017) That is, when discussing the active mobility tensors ${\bm{\Lambda}}_{ij}^{\mathrm{tt}}$ and ${\bm{\Lambda}}_{ij}^{\mathrm{rt}}$ for $i\neq j$, in effect we only consider terms in $a_{j}/L_{j}$ to leading order.

Finally, the forces and torques in Eq. (14) remain to be defined. First, we set the overall potential in our system of $N$ swimmers as

[TABLE]

Here, the external potentials $u^{k}_{\mathrm{ext}}$ can differ for different particles $k$. Additionally, pairwise steric interactions are introduced via $u^{kl}({|{\mathbf{r}}_{k}-{\mathbf{r}}_{l}|})$, which we specify for the applications in Sec. IV as the GEM-4 potential Mladek et al. (2006); Archer et al. (2014)

[TABLE]

with the potential strength $\epsilon^{kl}_{0}$ and the effective diameter $\sigma_{kl}=(\sigma_{k}+\sigma_{l})/2$ being set for each pair $k$ and $l$.

The forces ${\mathbf{F}}_{j}$ in Eq. (14) then read

[TABLE]

where the effect of thermal forces enters via the first term based on the effective entropic potential,Doi and Edwards (1988) which involves the microstate probability density $P=P({\mathbf{X}}^{N},t)$, the Boltzmann factor $k_{\mathrm{B}}$, and the temperature $T$ of the system. This expression ensures that the correct (translational) diffusion terms eventually appear in the statistical description in Sec. III.

Similarly, the torques in Eq. (14) are given by

[TABLE]

Again, this expression correctly reproduces (rotational) diffusion in the statistical description, see Sec. III.

III Derivation of the dynamical density functional theory

In this section, we derive the partial differential equations describing the dynamical microscopic statistics of a multi-component microswimmer system via dynamical density functional theory (DDFT), building on the derivations of the one-component case.Menzel et al. (2016); Hoell et al. (2017) For this purpose, the hydrodynamic swimmer model overviewed in Sec. II is used as an input. The resulting theory covers, combines, and extends several previously considered theories for systems of, e.g., one-species microswimmer suspensions,Menzel et al. (2016) “dry” — i.e., not hydrodynamically-interacting — self-propelled particles,Wensink and Löwen (2008); Wittkowski and Löwen (2011) hydrodynamically interacting passive colloidal particles,Rex and Löwen (2009) and binary mixtures of dry passive colloidal particles.Archer (2005)

First, we specify our system, which contains two different species of microswimmers suspended in a surrounding bulk fluid. For these species, the number of corresponding swimmers in the system is given by $N_{\mathrm{A}}$ and $N_{\mathrm{B}}$, respectively, adding up to a total of $N=N_{\mathrm{A}}+N_{\mathrm{B}}$ swimmers. Here, we order the swimmers by species, such that swimmers $1,\dots,N_{\mathrm{A}}$ belong to the first species and swimmers $N_{\mathrm{A}}+1,\dots,N$ to the second species. Additionally, a constant temperature $T$ of the fluid and a constant volume of the system are assumed. We adhere to the swimmer model introduced in Sec. II, with all swimmers of species $\nu\in\{\mathrm{A},\mathrm{B}\}$ featuring the same parameters $a_{\nu}$, $f_{\nu}$, $\alpha_{\nu}$, $L_{\nu}$, and $\sigma_{\nu}$. Setting $f_{\nu}=0$, also passive particles can be described accordingly, i.e., active–passive mixtures are covered by our theory as well.

Our starting point to derive the statistical description is the many-body Smoluchowski equation Doi and Edwards (1988)

[TABLE]

for the overdamped dynamics of our microswimmers. Here, $P=P({\mathbf{X}}_{1},...,{\mathbf{X}}_{N},t)$ denotes the microstate probability density of the corresponding configuration at time $t$. The velocities ${\mathbf{v}}_{i}$ and the angular velocities ${\bm{\omega}}_{i}$ are again related to the forces, torques, and the self-propulsion mechanisms via Eq. (14).

Next, we introduce ${\mathbf{X}}^{m}_{\mathrm{A}}=\{{\mathbf{X}}_{1},\dots,{\mathbf{X}}_{m}\}$ and ${\mathbf{X}}^{n}_{\mathrm{B}}=\{{\mathbf{X}}_{N_{\mathrm{A}}+1},\dots,{\mathbf{X}}_{N_{\mathrm{A}}+n}\}$ as short notations for the sets containing the phase-space coordinates of the first $m$ swimmers of species A and, respectively, the first $n$ swimmers of species B in the system. Since all swimmers are identical, we now define, for $m\leq N_{\mathrm{A}}$ and $n\leq N_{\mathrm{B}}$, the reduced $(m,n)$-swimmer density $\rho^{(m,n)}({\mathbf{X}}^{m}_{\mathrm{A}},{\mathbf{X}}^{n}_{\mathrm{B}},t)$ of finding (any) $m$ swimmers of species A and (any) $n$ swimmers of species B at the coordinates indicated in the argument. It is obtained from the full probability distribution $P({\mathbf{X}}^{N_{\mathrm{A}}}_{\mathrm{A}},{\mathbf{X}}^{N_{\mathrm{B}}}_{\mathrm{B}},t)$ by integrating out the degrees of freedom of $N_{\mathrm{A}}-m$ swimmers of species A and $N_{\mathrm{B}}-n$ swimmers of species B, reading

[TABLE]

Here, the prefactors result from the considered indistinguishability between swimmers of the same species. Particularly, we define the one-swimmer densities $\rho_{\mathrm{A}}({\mathbf{X}},t):=\rho^{(1,0)}({\mathbf{X}},t)$ and $\rho_{\mathrm{B}}({\mathbf{X}},t):=\rho^{(0,1)}({\mathbf{X}},t)$. Instead of referring to one specific swimmer, the coordinates ${\mathbf{X}}$ now simply identify “a swimmer” of the corresponding species. Furthermore, reduced densities with $m+n=2$ ($m+n=3$) will be referred to as two-swimmer (three-swimmer) densities below.

Our aim is to derive a physically well-grounded, closed set of coupled dynamical equations for the two one-swimmer densities. Thus, eventually, there shall be no remaining explicit dependence on the (generally unknown) higher-order densities. The starting point for our derivation is the many-body Smoluchowski equation given in Eq. (31). We first integrate out all swimmer coordinates except for those of one swimmer of species A. Second, we integrate out in the initial Eq. (31) all swimmer coordinates except for those of one swimmer of species B. As a result, we obtain one dynamical equation for $\rho_{\mathrm{A}}({\mathbf{X}},t)$ and one for $\rho_{\mathrm{B}}({\mathbf{X}},t)$, respectively. These equations (given below) form a coupled set, but at this point still contain higher-order densities and thus require an additional closure, as will be addressed afterwards via methods of dynamical density functional theory.

The corresponding equation for species A reads

[TABLE]

In this expression, the current densities labeled as ${\bm{\mathcal{J}}}_{\!\!\cdot}^{\cdot\cdot}$ do not involve hydrodynamic interactions between swimmers. These current densities are given by

[TABLE]

In contrast to that, the current densities involving hydrodynamic interactions between pairs of swimmers of species A follow as

[TABLE]

Third, the current densities associated with hydrodynamic effects of swimmers of species B on swimmers of species A are

[TABLE]

Here, the tensors ${\bm{\mu}}^{\cdot\cdot}_{\cdot\cdot}$ and ${\bm{\Lambda}}^{\cdot\cdot}_{\cdot\cdot}$ follow from the definitions in Eqs. (15)–(25) by inserting the parameters corresponding to the (phase-space) coordinates given in the subscripts and the combination of species referred to in the superscripts. The current densities ${\bm{\mathcal{K}}}_{\!\mathrm{AA}}^{\mathrm{tr}}$, ${\bm{\mathcal{K}}}_{\!\mathrm{AA}}^{\mathrm{rr}}$, ${\bm{\mathcal{K}}}_{\!\mathrm{AB}}^{\mathrm{tr}}$, and ${\bm{\mathcal{K}}}_{\!\mathrm{AB}}^{\mathrm{rr}}$ vanish for spherical swimmer bodies because the corresponding mobility tensors are independent of ${\mathbf{\hat{n}}}^{\prime}$, see Eqs. (16) and (17). Integrating the remaining gradient expressions over the closed surface of the unit sphere yields zero in each case. For non-spherical swimmer bodies, however, these current densities (as well as ${\bm{\mathcal{J}}}_{\!\!\mathrm{A}}^{\mathrm{tt}}$, ${\bm{\mathcal{J}}}_{\!\!\mathrm{A}}^{\mathrm{rt}}$, and ${\bm{\mathcal{J}}}_{\!\!\mathrm{A}}^{\mathrm{ra}}$) could be non-zero. Moreover, we remark that all ${\bm{\mathcal{K}}}_{\!}$’s become zero if hydrodynamic interactions are neglected.

An analogous dynamical equation for $\rho_{\mathrm{B}}({\mathbf{X}},t)$ follows by replacing $\mathrm{A}\to\mathrm{B}$, $\mathrm{B}\to\mathrm{A}$, and $\rho^{(m,n)}\to\rho^{(n,m)}$. Moreover, because of our convention of ordering species coordinates by first A and then B, we need to replace

[TABLE]

Obviously, the above equations (34)–(49) depend on unknown higher-order densities. In principle, one can now find dynamical equations for these quantities by applying corresponding integral operations on Eq. (31), but the resulting equations again contain unknown densities of even higher order. This escalating loop is typical for BBGKY-like hierarchiesHansen and McDonald (1990) and must be truncated and closed by appropriate approximations of the higher-order densities, e.g., as functions of the one-swimmer densities. In the following, DDFT methods will be employed for this purpose.

The main step in DDFT Evans (1979, 1992); Marconi and Tarazona (1999, 2000); Archer and Evans (2004); Chan and Finken (2005); Español and Löwen (2009); Evans (2010); Löwen (2010); Wittkowski and Löwen (2011) is the adiabatic approximation. It transfers equilibrium closure relations established in (classical) density functional theory (DFT) to the non-equilibrium case. Particularly, DDFTs imply that the higher-order densities relax faster than the one-swimmer densities,Español and Löwen (2009) as is conceivable for typical overdamped systems of colloidal particles (i.e., at low Reynolds numbers) and thus also for microswimmers.Menzel et al. (2016)

In equilibrium, DFT states that each observed density profile results from exactly one, uniquely specified external potential working on the corresponding particles.Singh (1991); Evans (1992); Chan and Finken (2005); Löwen (2010); Marconi and Tarazona (1999, 2000); Archer and Evans (2004) We call these potentials $\Phi^{\nu}_{\mathrm{ext}}({\mathbf{X}})$, $\nu=\mathrm{A,B}$, for the two species in our case. DDFT assumes these relations to hold at any time $t$. Thus, the external DFT potentials become time-dependent, and we denote them by $\Phi^{\nu}_{\mathrm{ext}}({\mathbf{X}},t)$. We remark that the equilibrium relations strictly hold only for $f_{\nu}=0$, $\nu=\mathrm{A,B}$, i.e., for passive particles. This limits the applicability of the theory when activity-induced correlation effects in the higher-order densities dominate the behavior of the system. Nevertheless, the overdamped nature of the systems favors the DDFT approach. Previously, bulk swimmer–swimmer pair distribution functions have been determined Hoell et al. (2018) by combining DDFT with a Percus-like Percus (1962) test-particle protocol.

We now discuss the above-introduced virtual external potentials, which may (and generally will) differ for the two species. In contrast to the “real” external potential introduced in Eq. (27), a dependence on the orientations of the swimmers here is allowed, and indeed even needed when the distributions of the orientations become non-uniform.

It must be stressed that these virtual potentials do not need to be determined explicitly. Repeating usual steps in derivations of DDFTs, we will in the following show two different ways of expressing $\Phi^{\nu}_{\mathrm{ext}}({\mathbf{X}},t)$ so that they can be eliminated from the mathematical description. Accordingly, we obtain expressions that help us to close the above BBGKY-like set of equations.

We start from the equilibrium grand potential as a functional of the one-swimmer densities, which is minimal for the equilibrium density distributions. The general ansatz for this functional can be written asArcher (2005)

[TABLE]

Here, all terms on the right-hand side except for the last one are known analytically. Namely,

[TABLE]

$\nu=\mathrm{A,B}$, is the ideal gas part, with $\lambda_{\nu}$ the corresponding thermal de Broglie wavelength $\lambda_{\nu}$. The contributions due to the external DFT potentials read

[TABLE]

$\nu=\mathrm{A,B}$. For our purposes, we may assume the chemical potentials to be combined with the external potentials.

Finally, the third contribution $\mathcal{F}_{\!\!\mathrm{exc}}$ includes interactions between the particles and represents the excess free energy beyond the ideal gas part. In almost all situations, an exact expression for $\mathcal{F}_{\!\!\mathrm{exc}}$ is not known analytically, and it must be approximated by an appropriate functional depending on the case at hand. Typically, this assumption needs to be carefully tested against experimental and simulational results. Nevertheless, the general theoretical framework up to this point applies to any interaction potential, here independent of the orientations of the swimmers (in principle, this restriction could be lifted, e.g., when describing rod-like active particlesRex et al. (2007); Wensink and Löwen (2008)).

In equilibrium, the actual magnitude of the grand potential is found by minimizing the grand potential functional over all possible density distributions. Thus, the equilibrium density fields $\rho_{\nu}^{\mathrm{eq}}({\mathbf{X}})$ satisfy

[TABLE]

for $\nu=\mathrm{A},\mathrm{B}$. Inserting Eqs. (50)–(52) leads to

[TABLE]

for $\nu=\mathrm{A,B}$.

Second, we employ standard equilibrium statistical mechanics.Gray and Gubbins (1984) In equilibrium, the static system properties are set completely by the temperature and the overall potential $U=U({\mathbf{X}}_{1},\dots,{\mathbf{X}}_{N})$ as defined in Eq. (27), writing $\Phi^{\nu}_{\mathrm{ext}}({\mathbf{X}})$ instead of $u^{\nu}_{\mathrm{ext}}({\mathbf{r}})$. Thus, the microstate probability density is given by

[TABLE]

where $\beta=(k_{\mathrm{B}}T)^{-1}$.

Applying the gradient with respect to the position of the first swimmer, which is of species A, leads to

[TABLE]

Since swimmers of the same species are considered to be identical and indistinguishable, we may write

[TABLE]

after integrating over the coordinates of all but the first swimmer of species A and using Eq. (32). This constitutes a lowest-order member of the binary-mixture translational Yvon-Born-Green (YBG) relations.Hansen and McDonald (1990); Gray and Gubbins (1984) Combining Eqs. (54) and (57), $\Phi^{\mathrm{A}}_{\mathrm{ext}}({\mathbf{X}})$ is eliminated and

[TABLE]

is obtained. Here, we now applied the adiabatic approximation and also switched to a time-dependent description. This equation is inserted into Eq. (34) on our way of closing our dynamical equations.

Based on Eqs. (32), (54), and (55), i.e., again applying the adiabatic approximation, we find two further helpful relations, namely

[TABLE]

and

[TABLE]

Analogues for species B follow after applying to Eqs. (58)–(60) the replacements listed below Eq. (49).

Inserting the above relations into Eqs. (34), (38), (41), (44), and (47) yields

[TABLE]

respectively, where we defined the vector fields

[TABLE]

This way, the two-swimmer density in Eq. (34) and all three-swimmer densities have been eliminated. Again, analogous relations apply to the dynamical equation for $\rho_{\mathrm{B}}({\mathbf{x}},t)$ and are obtained by considering the replacements introduced below Eq. (49).

Still, the remaining two-swimmer densities in the ${\bm{\mathcal{K}}}_{\!\,\cdot}^{\cdot\cdot}$ current densities must be addressed. For this purpose, as in a previous work,Hoell et al. (2017) we employ the Onsager-type Onsager (1949) approximations

[TABLE]

Here, for $|{\mathbf{r}}-{\mathbf{r}}^{\prime}|$ smaller than the sum of the radii of the involved swimmer bodies, we furthermore set the pair densities to zero to avoid the otherwise-appearing unphysical hydrodynamic divergences. Strictly speaking, this leads to a discontinuity, but typically the jump is vanishingly small, e.g., $\exp(-5\exp(-1/16))\approx 0.009\ll 1$ for $\epsilon_{0}^{\cdot\cdot}=5k_{\mathrm{B}}T$ and $a_{\cdot}=\sigma_{\cdot}/4$, see Eq. (28). This order of magnitude is sufficiently low to treat the function as basically “smooth” in the numerical evaluation.

Equations (67)–(69) implicitly assume $g_{\mu\nu}({\mathbf{X}},{\mathbf{X}}^{\prime},t)\approx\exp\left(-\beta u^{\mu\nu}(|{\mathbf{r}}-{\mathbf{r}}^{\prime}|)\right)$ for the pair distribution functions, with $\mu,\nu\in\{\mathrm{A,B}\}$. Using these relations is exact for passive equilibrium systems in the low-density limit,Hansen and McDonald (1990) as the expressions are based on the assumption that the two involved particles interact only with each other (and with no third particles). Adapting these relations to describe semi-dilute active suspensions thus constitutes a reasonable first-order approximation beyond assuming a constant pair distribution function. More generally, one could at this point also insert another reasonable approximation for the pair distribution function.

Similarly, our (pairwise) treatment of hydrodynamic interactions between the swimmers, see Eqs. (14)–(25), requires sufficiently large distances between the swimmer bodies. First, this is ensured by the steric interaction between the swimmers when half of its effective range, i.e., $\sigma_{\mu\nu}/2$ in Eq. (28), is larger than $a_{\kappa}$, $\alpha_{\kappa}L_{\kappa}$, and $(1-\alpha_{\kappa})L_{\kappa}$, with $\mu,\nu\in\{\mathrm{A,B}\}$ and $\kappa\in\{\mu,\nu\}$. The larger the mean distances are between the swimmers, the higher the accuracy of our description of hydrodynamic interactions will be. Together with the assumptions involved in Eqs. (67)–(69), we thus expect our DDFT for multi-species systems of microswimmers to perform best for (semi-)dilute suspensions of swimmers, within which our steric interaction potentials maintain a significant distance between the swimmer bodies, even when they are heading for collisions.

Finally, the excess functional ${\mathcal{F}}_{\!\!\mathrm{exc}}$ involving the effective steric interactions between the swimmers needs to be specified. As appropriate for GEM potentials,Archer et al. (2014) we from now on use a mean-field approximation, here for our case of binary mixtures, reading

[TABLE]

with $\mu,\nu\in\{\mathrm{A,B}\}$ and summing over repeated indices. In this way, our set of coupled dynamical equations for $\rho_{\mathrm{A}}({\mathbf{X}},t)$ and $\rho_{\mathrm{B}}({\mathbf{X}},t)$ is closed. We remark that, along the same lines, a theory for more than two different species can be derived as well, leading to a correspondingly further increased number of terms. Here, we continue by applying the above theory to concrete example situations in Sec. IV.

IV Applications

In this section, the DDFT derived in Sec. III is applied to several illustrative cases. Specifically, for simplicity, these will be setups in which the positions and orientations of the swimmers are constricted to the $xy$-plane. Still, a surrounding bulk fluid is considered with the planar swimmer ensemble embedded therein, allowing for three-dimensional fluid flows. Possible methods to experimentally realize this situation could be the confinement of microswimmers to the interface between two immiscible fluids of identical viscosity $\eta$, or the use of optical trapping fields.

In such a setup, the orientation of a swimmer is described by a single angle $\phi$ (measured from the $x$-axis) via ${\mathbf{\hat{n}}}=(\cos\phi,\sin\phi)$. The orientational gradient operator then reduces to ${\mathbf{\hat{n}}}\times\nabla_{{\mathbf{\hat{n}}}}={\mathbf{\hat{z}}}\partial_{\phi}$, where ${\mathbf{\hat{z}}}$ is the oriented Cartesian unit vector pointing (upwards) out of the $xy$-plane. Furthermore, the phase-space coordinate ${\mathbf{X}}$ in this situation becomes ${\mathbf{X}}=\{x,y,\phi\}$.

The numerical solution of the coupled set of partial differential equations derived in Sec. III is then performed on an equidistant $N_{x}\times N_{y}\times N_{\phi}$ grid using the finite-volume-method solver FiPyGuyer et al. (2009). Formally, numerical periodic boundary conditions are imposed on all coordinates $x$, $y$, and $\phi$, but hydrodynamic and steric interactions are cut at a distance chosen such that no (unphysical) interactions across the boundaries occur. As nevertheless all physical interactions inside the system should, of course, be accounted for, we further always set the length of the simulation box in both spatial directions to at least twice the largest relevant interparticle distance.

Since the orientation-dependent densities $\rho_{\nu}({\mathbf{X}},t)$ at time $t$ are still a function of $x$, $y$, and $\phi$, they cannot be easily plotted even for our planar configurations. For displaying our results, we thus further define the (orientation-integrated) spatial swimmer densities

[TABLE]

and the orientational vector fields

[TABLE]

where $\nu\in\{\mathrm{A,B}\}$. Moreover, the overall (average) one-species densities are described by $\bar{\rho}_{\nu}=A^{-1}\!\!\int_{A}\!\mathrm{d}{\mathbf{r}}\,\rho_{\nu}({\mathbf{r}},t)$, where $A$ is the area of the regarded system.

IV.1 Trapped binary swimmer system

While restricting the binary microswimmer configuration to two spatial dimensions as detailed above, we now additionally introduce radially-symmetric quartic trapping potentials given by

[TABLE]

with potential strengths $V^{\nu}_{0}$, distance $r=|{\mathbf{r}}|$ to the center of the trap, and $\nu=\mathrm{A,B}$. As in previous works,Menzel et al. (2016); Hoell et al. (2017) we use a quartic potential — instead of, e.g., a harmonic one ($\propto r^{2}$) — to observe more pronounced differences between activity-induced off-center density distributions (see below) and center-heavy equilibrium distributions for passive particles. Previously reported results for harmonic trapsNash et al. (2010); Hennes et al. (2014) showed qualitative agreement with our results for a quartic potential.Menzel et al. (2016); Hoell et al. (2017) For simplicity, we furthermore from now on assume that all species-related parameters are the same for both species, except for $f_{\mathrm{A}}=-f_{\mathrm{B}}>0$. Thus, species A is formed by pushers and species B represents pullers (of the same strength).

In analogous one-component suspensions,Menzel et al. (2016); Hoell et al. (2017) without any active drive, the external potential leads to center-heavy distributions following standard equilibrium statistics. When the active drive is switched on in the one-component systems, but hydrodynamic interactions are still neglected, the self-propelled particles start forming a radially symmetric high-density ring, along which the outward self-propulsion is balanced by the restoring trapping force.Menzel (2015b); Yan and Brady (2015); Menzel et al. (2016) With hydrodynamic interactions incorporated, this ring of microswimmers can become unstable against collapsing to one spot on this ring, which is induced by the hydrodynamic coupling through the resulting fluid flows.Nash et al. (2010); Hennes et al. (2014); Menzel et al. (2016); Hoell et al. (2017) In parts of the parameter space, pushers and pullers were observed to behave quite differently, with pushers showing a significantly more pronounced destabilization of the high-density ring and formation of a high-density spot, while pullers showed a much weaker density variation along the ring.Hoell et al. (2017)

We are now interested in pusher–puller mixtures. There is a crucial competition between hydrodynamic effects resulting from the external potential acting on the swimmer bodies and from the actively introduced forces exerted by the microswimmers themselves. We here concentrate on a parameter range for which the hydrodynamic interactions induced by the self-propulsion mechanism dominate those induced by the external potential force. Concerning our current densities, we thus always check that $|{\bm{\mathcal{K}}}_{\!\,\cdot\cdot}^{\mathrm{ra}}|>|{\bm{\mathcal{K}}}_{\!\,\cdot\cdot}^{\mathrm{rt}}|$, see, e.g., Eqs. (47) and (49), for our chosen parameters.

Numerical results for (steady-state) distributions of pusher–puller mixtures are shown in Fig. 2, for varying overall densities of the two species. In strong contrast to the corresponding one-component systems, for which the (steady-state) distributions strongly differed between pure pusher and pure puller systems,Hoell et al. (2017) we here frequently observe the same qualitative behavior when both species are present simultaneously. For instance, in Fig. 2 (a), pushers transfer their “spot-forming” tendency onto the pullers, which in the absence of the pushers would show a ring-like arrangement instead of the spot. However, the plots in Fig. 2 indicate the rough relation $\rho_{\mathrm{A}}({\mathbf{r}},{\mathbf{\hat{n}}},t)/\bar{\rho}_{\mathrm{A}}\approx\rho_{\mathrm{B}}({\mathbf{r}},{}-{\mathbf{\hat{n}}},t)/\bar{\rho}_{\mathrm{B}}$. Choosing, e.g., $|f_{\mathrm{A}}|\neq|f_{\mathrm{B}}|$, this approximate relation breaks down as the two species aggregate at different distances from the origin, but for sufficiently small deviations, we still observe a qualitatively similar collective behavior for both species.

In Fig. 2, the overall density $\bar{\rho}_{\mathrm{B}}$ of pullers increases from left to right, while the overall density $\bar{\rho}_{\mathrm{A}}$ for pushers decreases from the top row to the bottom row. We observe clear spot formation in Fig. 2 (a) and (b), while (c) shows less-pronounced instabilities of the high-density ring. Thus, we may conclude, that the dominating species imposes its behavior onto the other species.

For Fig. 2 (c) and (e), where $\bar{\rho}_{\mathrm{A}}=\bar{\rho}_{\mathrm{B}}$ and therefore $\rho_{\mathrm{A}}({\mathbf{r}},{\mathbf{\hat{n}}},t)\approx\rho_{\mathrm{B}}({\mathbf{r}},-{\mathbf{\hat{n}}},t)$ holds, the probability currents associated with the rotation due to the active forces approximately cancel each other by symmetry, e.g., ${\bm{\mathcal{K}}}_{\!\mathrm{AA}}^{\mathrm{ra}}\approx{}-{\bm{\mathcal{K}}}_{\!\mathrm{AB}}^{\mathrm{ra}}$, so that only the currents ${\bm{\mathcal{K}}}_{\!\,\cdot\cdot}^{\mathrm{rt}}$ can lead to spot formation. The latter starts to outperform the rotational diffusion for the case depicted in Fig. 2 (c), but not for the lower overall densities in Fig. 2 (e). The instability of the ring here seems to be a question of high-enough overall density because, e.g., $|{\bm{\mathcal{K}}}_{\!\mathrm{AA}}^{\mathrm{rt}}|\propto\bar{\rho}_{\mathrm{A}}^{2}$ and $|{\bm{\mathcal{J}}}_{\!\!\mathrm{A}}^{\mathrm{rr}}|\propto\bar{\rho}_{\mathrm{A}}$.

The bottom row of Fig. 2 shows the corresponding density distributions for a smaller $\bar{\rho}_{\mathrm{A}}$. Thus, a decreased density of pushers leads to an increased stability of the high-density ring against aggregation in one spot. When (significantly) more pullers than pushers are in the system, as in Fig. 2 (f), they dominate the overall behavior and re-stabilize the high-density ring.

In summary, the majority species seems to dominate the overall behavior of the system. A similar conclusion has recently been drawn for the unconfined motion in pusher–puller mixtures,Pessot et al. (2018) which we will treat as the next example using our theoretical approach.

At this point, we include a short remark on the performance of our theory. We can remove the second species from our DDFT equations derived in Sec. III by setting $\rho_{\mathrm{B}}({\mathbf{X}},t)\equiv 0$. Then, the present set of equations reduces to the previous DDFT for monodisperse microswimmers.Menzel et al. (2016) In that case, likewise, the statistical theory was evaluated by exposing the system of swimmers to a radial external trapping potential, in analogy to the above consideration for a pusher–puller mixture. There, hydrodynamic interactions lead to the formation a high-density spot of aligned swimmers as well, resulting in overall flow fields.Menzel et al. (2016); Hoell et al. (2017) This “hydrodynamic fluid pump” had previously been reported in particle-based computer simulations,Nash et al. (2010); Hennes et al. (2014) using different swimmer models. Thus, a qualitative comparison shows that our DDFT reproduces corresponding general phenomena. Adding another microswimmer species to the same framework, we expect a similarly successful performance of the present theory. Direct quantitative comparison could be carried out in the future by implementing a suitable particle picture into many-swimmer computer simulations including hydrodynamic interactions and thermal fluctuations, e.g., via multiparticle collision dynamicsMalevanets and Kapral (1999); Kapral (2008); Gompper et al. (2009); Zöttl and Stark (2012, 2014); Rühle et al. (2018) / stochastic rotation dynamics.Ihle and Kroll (2001); Downton and Stark (2009) Then, also higher swimmer densities could be addressed numerically. Another way to explicitly take into account the induced hydrodynamic fluid flows in computer simulations could be Lattice-Boltzmann methods.McNamara and Zanetti (1988); Nash et al. (2010); Alarcón and Pagonabarraga (2013); Kuron et al. (2019); Bardfalvy et al. (2019)

IV.2 Emergence of polar orientational order and collective motion in pusher–puller mixtures

In the absence of the spherical trapping potential considered in Sec. IV.1, previous particle-based computer simulations of planar arrangements of microswimmers with periodic boundary conditions and using the same swimmer model have identified a tendency of puller microswimmers to develop (global) collective polar orientational order.Pessot et al. (2018) Related observations were made in simulations of analogous three-dimensional configurations of squirmer microswimmers.Alarcón and Pagonabarraga (2013) Such order in the swimmer orientations naturally leads to collective motion, maintaining a common average propulsion direction. Moreover, we have performed a corresponding linear stability analysis of our DDFT for planar pure (one-species) pusher or puller systems, with spontaneous ordering identified beyond a threshold active drive for pullers,Hoell et al. (2018) in contrast to pushers. We now address the corresponding two-species situation. In related computer simulations for mixtures of pushers and pullers using the same swimmer model,Pessot et al. (2018) it was found that collective orientational order only develops if the fraction of pushers is sufficiently small. As we demonstrate, our DDFT reproduces these results and leads to a more quantitative insight.

For this purpose, the external potential in our planar arrangement is now set to $u_{\mathrm{ext}}({\mathbf{r}})\equiv 0$. For simplicity, we assume that the one-swimmer densities are spatially homogeneous, i.e., $\rho_{\nu}({\mathbf{X}},t)=\rho_{\nu}(\phi,t)/A$, with $A$ denoting the area (considered to be large) of the periodic plane containing the swimmers and $\nu\in\{\mathrm{A,B}\}$. Then, integrating Eq. (33) over all positions ${\mathbf{r}}$ in the periodic box leads to

[TABLE]

with the probability current densities defined in Eqs. (34)–(49). Following Ref. Hoell et al., 2018, the current densities ${\bm{\mathcal{K}}}_{\!\cdot\cdot}^{\mathrm{rt}}$ are neglected for sufficiently dilute suspensions, as all the contained non-vanishing terms scale with three-swimmer densities. Thus, Eq. (74) reduces to

[TABLE]

Here, the two-swimmer densities are related to the pair distribution functions via

[TABLE]

Thus, Eq. (75) becomes

[TABLE]

where the hydrodynamic interactions are comprised in the coupling functions

[TABLE]

with $\mu,\nu\in\{\mathrm{A,B}\}$. An analogous dynamical equation for species B is obtained by replacing $\mathrm{A}\to\mathrm{B}$ and $\mathrm{B}\to\mathrm{A}$. In the following, species A again represents pushers, and species B pullers.

To allow for further analytical treatment, we include additional simplifying assumptions. Considering systems in which all active agents propel with the same amplitude of the active drive and further are identical in all other microscopic parameters, the coupling and pair distribution functions, see Eq. (79), were determined in Ref. Hoell et al., 2018 by a modified Percus test-particle method. For this purpose, hydrodynamic interactions were neglected and only the interplay of self-propulsion and steric interactions was evaluated. As a result, we had extracted and approximated the basic functional form as Hoell et al. (2018)

[TABLE]

where $\tilde{C}_{\mathrm{AA}}=\tilde{C}_{\mathrm{BB}}=\tilde{C}/A>0$ is positive for same-species coupling, and $\tilde{C}_{\mathrm{AB}}=\tilde{C}_{\mathrm{BA}}={}-\tilde{C}/A$. This distinction follows from the fact of our puller microswimmers propelling into the direction of ${}-{\mathbf{\hat{n}}}$ and / or ${}-{\mathbf{\hat{n}}}^{\prime}$, see Fig. 1. Since $\phi$ and $\phi^{\prime}$ parameterize the orientations of ${\mathbf{\hat{n}}}$ and ${\mathbf{\hat{n}}}^{\prime}$, respectively, the swimming direction of a puller is shifted by an additional angle $\pi$ relatively to $\phi$ and / or $\phi^{\prime}$. If only one of the angles $\phi$ and $\phi^{\prime}$ refers to a puller, the additional shift of $\phi-\phi^{\prime}$ by $\pi$ requires a minus sign in the prefactor of $\sin(\phi-\phi^{\prime})$ in Eq. (80).

The value of $\tilde{C}>0$ generally depends on the overall density and the microscopic parameters. (Some further positive constant parameters are here incorporated by the coefficients $\tilde{C}$ when compared to the amplitude $C$ in Ref. Hoell et al., 2018.) Since a similarly simple analytically treatable expression is still missing for hydrodynamic interactions included on the level of pair distribution functions, we use Eq. (80) as an input for our further calculations.

We assume that, if collective order arises, there is only one common direction of polar ordering, i.e., that in this case species A and B collectively propel along a common direction. This assumption is motivated by previous simulation results.Pessot et al. (2018) We now test the linear stability of the uniform distributions $\rho_{\nu}(\phi,t)\equiv N_{\nu}/(2\pi)$ against the emergence of collective orientational ordering. To this end, the ansatz $\rho_{\mathrm{A}}(\phi,t)=N_{\mathrm{A}}/(2\pi)+\epsilon_{\mathrm{A}}(t)\cos(\phi-\phi_{0})$ and $\rho_{\mathrm{B}}(\phi,t)=N_{\mathrm{B}}/(2\pi)+\epsilon_{\mathrm{B}}(t)\cos(\phi-\phi_{0}+\pi)$ , with an arbitrary angle $\phi_{0}$ and $|\epsilon_{\nu}(t)|\ll N_{\nu}$ for $\nu\in\{\mathrm{A,B}\}$, is inserted into Eq. (78) and the equivalent equation for species B. This leads to the coupled ordinary differential equations

[TABLE]

with the coefficient matrix

[TABLE]

where $m_{\mu\nu}:=N_{\mu}f_{\nu}\tilde{C}/(2A)$.

We recall that species A (pushers) and species B (pullers) are considered to have the same amplitude of their active drive, i.e., $f_{\mathrm{A}}=-f_{\mathrm{B}}>0$. Additionally, we keep $N_{\mathrm{A}}+N_{\mathrm{B}}=N$ constant, i.e., only the ratio of pushers to pullers is varied. Moreover, all other parameters are assumed to be identical for the two species. Then, the eigenvalues of ${\mathbf{M}}$ are determined as $\left(\!{}-k_{\mathrm{B}}T\mu^{\mathrm{r,A}},{}-k_{\mathrm{B}}T\mu^{\mathrm{r,A}}+f_{\mathrm{A}}\tilde{C}N(\chi_{\mathrm{B}}-1/2)/A\right)$. Here, the first eigenvalue is always negative, but the second one becomes positive if

[TABLE]

with $\chi_{\mathrm{B}}:=N_{\mathrm{B}}/N$ denoting the fraction of pullers. The corresponding eigenvector is ($N_{\mathrm{A}},N_{\mathrm{B}})$.

Our system can thus be linearly unstable against polar orientational ordering only if the right-hand side of Eq. (89) is positive. Since $f_{\mathrm{A}}>0$, this implies that the pullers must outnumber the pushers ($\chi_{\mathrm{B}}>1/2$). If this condition is satisfied, the active drive additionally needs to be strong enough, i.e., Eq. (89) sets a threshold strength for $f_{\mathrm{A}}=-f_{\mathrm{B}}$. Particularly, the effect of the active drive and the hydrodynamic interactions need to outperform rotational diffusion. Furthermore, as indicated by the corresponding eigenvector $(N_{\mathrm{A}},N_{\mathrm{B}})$, if orientational order arises, it does so simultaneously for both species.

Our results roughly agree with those in the previous simulation study.Pessot et al. (2018) We stress that our theory only tests linear instability with respect to polar orientational ordering and that the above approximations were involved. In particular, the influence of hydrodynamic interactions on the pair distribution function was neglected. To address this question, possibly the results of particle-based computer simulations could be used as an input to the theory in the future.Pessot et al. (2018); Schwarzendahl and Mazza (2019) Since our previous theoretical analysis for single-species systems indicated polar orientational ordering for puller suspensions but not for pushers,Hoell et al. (2018) we again find that the majority species imposes its behavior onto the minority species, as observed already for the confined (trapped) mixtures in Sec. IV.1.

IV.3 Shear cell

As a third example, we now address a planar circular configuration which effectively represents a shear cell. We compose this shear cell of passive colloidal particles forming an effective circular rim and active microswimmers trapped inside. The passive particles are continuously driven along the circular rim of the trap, inducing a shear-like circular fluid flow inside. In a very loose analogy, this geometry is similar to setups of Taylor-Couette flow,Taylor (1923) but, of course, here in the limit of low Reynolds numbers. In fact, driving passive colloidal particles along ring-like trajectories can be realized experimentally via optical trapping potentials.Williams et al. (2016)

Considering the driven particles (that hydrodynamically interact with the interior microswimmers) as one component of a binary mixture naturally induces fluid flows to which the enclosed microswimmers are exposed. This avoids explicitly imposing such flows as an external flow field.Rauscher et al. (2007); Scacchi et al. (2017); Scacchi and Brader (2018); Stopper and Roth (2018) However, we do not account in the present work for, e.g., possible effects of shear banding, which have been addressed in the context of DDFT as well.Scacchi et al. (2017); Scacchi and Brader (2018); Stopper and Roth (2018) Our one-body density, particularly for passive particles within the cell, remains basically unchanged by the translational effects of the shear flow, as expected in the limits of our current theory regarding shear.Brader and Krüger (2011); Aerov and Krüger (2014) Instead, for active microswimmers within the cell, the induced rotation of the swimmer orientations, coupling to the directions of self-propulsion, can lead to changes in the spatial density.

In the context of our theory, the active microswimmers represent the first species A, while the driven colloidal particles are treated as species B. Consequently, $f_{\mathrm{B}}=0$, but we also define an effective potential of confinement

[TABLE]

for the passive particles, based on the error function $\mathrm{erf}(s)=(2/\sqrt{\pi})\int_{0}^{s}\mathrm{d}u\exp(-u^{2})$. For $V^{\mathrm{B}}_{0}\gg k_{\mathrm{B}}T$ and $R_{0}\gg\sigma_{\mathrm{R}}$, this potential effectively anchors the particles on a (small-width) ring of radius $R_{0}$. Additionally, the non-conservative driving force

[TABLE]

is taken into account to describe the continuous circular driving of the passive particles. Technically, we include it by adding $-{\mathbf{F}}_{\mathrm{d}}({\mathbf{r}})$ to $\nabla_{{\mathbf{r}}}u^{\mathrm{B}}_{\mathrm{ext}}({\mathbf{r}})$ in the corresponding equations. Here, $\omega_{\mathrm{d}}$ is the (signed) magnitude of the imposed (spatial) angular velocity with which the passive particles are driven along the ring.

For species A, we again choose the external trapping potential defined in Eq. (73), but take care when adjusting the potential strength that (even with $f_{\mathrm{A}}\neq 0$) it at all times hinders the majority of the swimmers from reaching the passive particles on the outer ring. This way, species A and B mainly interact with each other hydrodynamically, as described by, e.g., the current densities in Eqs. (64) and (65).

The driven ring of passive colloidal particles of species B is shown in Fig. 3. For typical parameters, (a) the corresponding density profile and (b) the hydrodynamic influences on the microswimmers of species A are depicted. For the latter, we define for species A the contribution to the velocity resulting from the fluid flows induced by species B as

[TABLE]

and the corresponding contribution to the $z$-component of the angular velocity as

[TABLE]

Here, the current densities, as defined in Eqs. (64) and (65) in combination with Eqs. (66) and (68), are proportional to $\rho_{\mathrm{A}}({\mathbf{X}},t)$ so that the above expressions do not diverge when the denominator vanishes.

The resulting density distribution of species B depicted in Fig. 3 is basically circularly symmetric and after initial equilibration does not vary over time any longer. Still, it represents the moving passive particles driven by the (tangential) external force defined in Eq. (91). The latter is the main source of the fluid flows induced by particles of species B. Resulting flow fields can be approximated inside the cell by evaluating the corresponding terms in Eqs. (92) and (93) under the assumption of $\rho_{\mathrm{B}}({\mathbf{X}}^{\prime},t)\equiv N_{\mathrm{B}}(2\pi)^{-2}R_{0}^{-1}\delta(r^{\prime}-R_{0})$. Considering the contribution of ${\mathbf{F}}_{\mathrm{d}}({\mathbf{r}}^{\prime})$ as dominant, ignoring steric interactions between species A and B, and introducing $b=r/R_{0}<1$, we obtain from Eq. (93)

[TABLE]

for the angular velocity. As shown in Fig. 3 (c), there is good quantitative agreement between this approximation [the integral expression in Eq. (94) is plotted as the dashed line] and the full numerical solution (solid line). For positions close to the outer ring of the driven particles of species B, the curve drops, most likely because of the decreased probability of finding the swimmers and the driven particles within close distances from each other, formally introduced by the Onsager-like terms in Eqs. (67)–(69). To leading order in $a$, the flow field induced by the driven species B can be similarly obtained as

[TABLE]

We now concentrate on species A that is confined inside the shear cell. For $f_{\mathrm{A}}=0$, passive particles are recovered. As seen in the steady states shown in Fig. 4 (a) and (b), the distribution of the inner passive particles remains virtually unaffected by the external driving of the outer passive particles, except for possible small deviations that cannot be resolved within the precision of our numerical discretization scheme. But when the enclosed swimmers are active ($f_{\mathrm{A}}\neq 0$), see Fig. 4 (c)–(j), the effects of the induced shear flows become significant. Figure 4 (c) and (d) show the situation of the enclosed swimmers for pushers and pullers without the external drive, i.e., $\omega_{\mathrm{d}}=0$. Here, for the chosen parameters, the microswimmers form high-density rings with average orientations tilted relatively to the outward direction for pushers (c) Hoell et al. (2017) and radially oriented for pullers (d). The directional sense of the tilt for pushers is spontaneously chosen by the system as either counterclockwise or clockwise (depicted here), depending on the initialization of our numerical evaluation. In contrast to these cases of vanishing external driving of species B, Fig. 4 (e)–(j) demonstrate that for $\omega_{\mathrm{d}}\neq 0$, the externally induced shear flows can lead to a collapse of the steady-state density distributions towards the center of the confinement. Moreover, with increasing external driving $\omega_{\mathrm{d}}$, both pushers and pullers furthermore show an increasing tendency of their locally averaged swimming direction to be reoriented by the externally imposed fluid flow [see Fig. 3 (b)]. This explains the different sense indicated by the white arrows for increased $\omega_{\mathrm{d}}$ from Fig. 4 (c) to Fig. 4 (e).

As a source of this behavior, the shear flow induced by the external driving of the outer ring persistently rotates the orientations of the internal swimmers so that the latter are hindered from efficiently swimming radially outwards against the trapping force. In this way, the behavior of species A becomes comparable to that of circle swimmers, i.e., self-propelled particles that additionally feature an active self-rotation.van Teeffelen and Löwen (2008); Fily et al. (2012); Kümmel et al. (2013); Löwen (2016); Sevilla (2016) Actually, we have observed a similar phenomenology as in Fig. 4 for increasing inherent curvature of the trajectories of circle swimmers in Ref. Hoell et al., 2017. In the present case, however, the (externally induced) rotation varies with the distance $r$ from the origin, so that the local radius of induced circle-swimming $R_{\mathrm{cs}}(r):=|{\mathbf{v}}_{0\mathrm{A}}/\omega_{\mathrm{AB}}(r)|$, determined from the definitions in Eqs. (26) and (94), is non-constant. It reaches a maximum at the origin and decreases with increasing $r$. For Fig. 4 (e)–(h), the length scale of $R_{\mathrm{cs}}(r)$ is comparable to the radius of the effective trap so that a high-density ring is still visible. However, the average orientations are significantly tilted from the radial direction (especially for pusher microswimmers). Further increasing the external driving strength, see Fig. 4 (i) and (j), leads to more localized density profiles and circling around the center of confinement.

The increasing localization can be quantified by the (negative) radial component of the averaged external trapping force experienced by the microswimmer ensemble, as given in Fig. 5 for the same (steady-state) data as in Fig. 4 (c)–(j). For vanishing angular driving speed $\omega_{\mathrm{d}}$ of the outer passive particles, we find a higher value for pullers (blue circles) than for pushers (red squares), corresponding to the more off-center density distribution of pullers caused by their stronger tendency to show radial orientation. Both curves drop for increasing $\omega_{\mathrm{d}}$. The reason is again the induced shear flow increasingly hindering the swimmers from self-propelling efficiently against the radial external trapping potential. The drop is somewhat delayed for our pullers, in accordance with a similar effect previously seen for circle swimmers in an external trap, where the pullers also showed a stronger tendency of maintaining a ring of outward-oriented swimmers.Hoell et al. (2017)

In related works, rosette-like trajectories have been reported for (single) circle swimmers with explicitly time-dependent self-propulsion velocities.Babel et al. (2014); Jahanshahi et al. (2017) Beyond the scope of the present work, when genuine circle swimmers are confined (as species A) in our setup, again high-density rings might be observed with average swimmer orientations along the local radial direction. For this purpose, the induced rotation $\omega^{\mathrm{AB}}$ should balance the inherent self-rotation of the circle swimmers.

V Conclusions

In this work, we have presented a dynamical density functional theory (DDFT) for multi-species suspensions of microswimmers. We have included (pairwise) hydrodynamic and effective steric interactions between swimmers. In this way, we conceptually extended the previous one-component equivalent.Menzel et al. (2016); Hoell et al. (2017, 2018) The theory is based on a discrete force-dipole minimal microswimmer model, which has already been used successfully in several previous works.Menzel et al. (2016); Hoell et al. (2017); Pessot et al. (2018); Hoell et al. (2018); Daddi-Moussa-Ider and Menzel (2018) We then applied our theory to three illustrative example situations of planar swimmer configurations inside a three-dimensional bulk fluid.

First, binary pusher–puller mixtures in external spherically symmetric trapping potentials have been discussed. For the two species only differing in their pusher / puller signature, we found that the majority species imposes its behavior on the minority species. For example, pushers at the considered propulsion strength on their own tend to form concentrated spots on the rim of the trap. Therefore, if pushers represent the majority species, this spot formation is conveyed to simultaneously present pullers. Conversely, pullers by themselves rather tend to form a roughly spherically symmetric high-density ring on the rim of the trap. Thus, if they represent the majority in a pusher–puller mixture, also pushers tend to organize themselves in a corresponding ring structure.

Second, in the absence of any external trapping potential, pusher–puller mixtures in large periodic boxes have been considered. In an analytical treatment analogous to the previously studied one-component case,Hoell et al. (2018) pullers are found to be able to establish the onset of collective polar orientational order of the whole mixture. Accordingly, pullers can induce oriented collective motion. For this purpose, they need to represent the majority species and show a sufficiently large magnitude of their active drive. Our results are qualitatively in line with previous agent-based computer simulations.Pessot et al. (2018)

Third, a microswimmer species is confined inside a circular ring of externally driven passive particles. The induced shear flow persistently rotates the confined swimmers, and thus can hinder them from forming the high-density rings that are typically observed for sufficiently quick self-propelled particles in radial trapping potentials. Instead, the swimmer densities tend to collapse towards the center of the confinement. Similar mechanisms have previously been found for circle swimmers (featuring an inherent self-rotation) without externally induced shear flows. One future task could be to focus further on the role of shear flows in our statistical theory.Brader and Krüger (2011); Aerov and Krüger (2014)

In the numerical examples, we have restricted our evaluations for hydrodynamically interacting microswimmers to small confined systems that suitably fit into the corresponding simulation box. Nevertheless, in the future, our set of partial differential equations could be solved numerically as well for (basically infinitely extended) bulk situations. For this purpose, (true) periodic boundary conditions are applied to a finite simulation box. Then, because of the long-range nature of the hydrodynamic interactions, the influence of all periodic images on the density distribution in the simulation box must be accounted for. Mathematically, this can be achieved by applying Ewald summation techniques Ewald (1921) to the mobility tensors. Corresponding results have been derived for passive particles,Beenakker (1986); Brady et al. (1988); Rex and Löwen (2008) but could, in principle, also be calculated for our active microswimmers, as has recently been demonstrated for a similar force-dipole-based microswimmer model.Adhyapak and Jabbari-Farouji (2018) For quantitative tests of our theory in the future and for extensions to higher densities, particle-based computer simulations (that include hydrodynamic flows of the surrounding fluid and thermal fluctuations explicitly) can be performed.

One very interesting question is whether our DDFT could be extended to describe the aforementioned motility-induced phase separation. In this context, existing statistical theories involved a density-dependent effective swimming speed and/or an anisotropic pair distribution function as additional inputs.Cates and Tailleur (2013); Bialké et al. (2013); Speck et al. (2014, 2015); Wittkowski et al. (2017) It would thus be interesting to study in the future the effect of at least one similar activity-induced term in our theory. Another promising statistical approach beyond the adiabatic approximation of DDFT is the power functional theory for “dry” self-propelled particles, which has recently been formulated and evaluated semi-analytically.Krinninger et al. (2016, 2017); Krinninger and Schmidt (2019)

The present work derives the multi-species DDFT for the case of uniaxial straight-swimming microswimmers with spherical bodies. However, only a few changes transfer it to the case of (inherently biaxial) circle swimmers.Hoell et al. (2017) Even more generally, changes will allow to describe swimmers with less-symmetric body shapes, e.g., rod-like, bodies. Nevertheless, we remark that more work is needed in the future regarding situations of still higher complexity. Examples are cases in which, for instance, additional phoretic chemical- or temperature-based interactions between swimmers become significant.Liebchen and Löwen (2019) Moreover, effects of the fluctuations of the propulsion mechanism itself could be taken into account.Pietzonka and Seifert (2017)

Beyond the direct numerical evaluations performed in this work, DDFTs can serve as a foundation to derive corresponding phase-field-crystal modelsvan Teeffelen et al. (2009); Menzel and Löwen (2013); Menzel et al. (2014); Wittkowski et al. (2014); Alaimo and Voigt (2018) and more macroscopic continuum theoriesHeidenreich et al. (2016); Reinken et al. (2018, 2019) of microswimmer suspensions. The latter allow for connections towards still-larger length scales of theoretical descriptions. Altogether, we thus expect our DDFT to provide a powerful tool for the statistical characterization of dynamic multi-species systems of suspended microswimmers, of future relevance both in fundamental physics and concerning the corresponding biological, technical, and medical applications.

Acknowledgements.

The authors thank Urs Zimmermann, Soudeh Jahanshahi, and Giorgio Pessot for helpful comments, as well as the Deutsche Forschungsgemeinschaft for support of this work through the SPP 1726 on microswimmers, Grant Nos. LO 418/17-2 and ME 3571/2-2.