Effective interactions between inclusions in an active bath

TL;DR

This study investigates how active particles mediate effective interactions between colloidal inclusions, revealing non-monotonic forces, the role of particle layering, and the influence of chirality on attraction and repulsion in two-dimensional active baths.

Contribution

It provides detailed insights into the mechanisms of multi-body interactions and the effects of chirality in active matter systems with colloidal inclusions.

Findings

Non-monotonic two-body force profiles with maxima and humps.

Formation of active particle rings around inclusions influences interactions.

Chirality induces depletion attractions between inclusions.

Abstract

We study effective two- and three-body interactions between non-active colloidal inclusions in an active bath of chiral or non-chiral particles, using Brownian Dynamics simulations within a standard, two-dimensional model of disk-shaped inclusions and active particles. In a non-chiral active bath, we first corroborate previous findings on effective two-body repulsion mediated between the inclusions by elucidating the detailed non-monotonic features of the two-body force profiles, including a primary maximum, and a secondary hump at larger separations that was not previously reported. We then show that these features arise directly from the formation, and sequential overlaps, of circular layers (or 'rings') of active particles around the inclusions, as the latter are brought to small surface separations. These rings extend to radial distances of a few active-particle radii from the surface of inclusions, giving the hard-core inclusions relatively thick, soft, repulsive 'shoulders', whose multiple overlaps then enable significant (non-pairwise) three-body forces in both non-chiral and chiral active baths. The resulting three-body forces can even exceed the two-body forces in magnitude and display distinct repulsive and attractive regimes at intermediate to large self-propulsion strengths. In a chiral active bath, we show that, while active particles still tend to accumulate at the immediate vicinity of the inclusions, they exhibit strong depletion from the intervening region between the inclusions, and partial depletion from relatively thick, circular, zones further away from the inclusions. In this case, the effective, predominantly repulsive, interactions between the inclusions turn to active, chirality-induced, depletion-type attractions, acting over an extended range of separations.

Equations

$${\dot{\mathbf{r}}}_i$$

$${\dot{\theta }}_i$$

$$V_{\mathrm{WCA}}({\mathbf{r}})=\left\{\begin{array}[]{l l}\!4\epsilon\!\left[\left(\frac{\sigma_{\mathrm{eff}}}{|{\mathbf{r}}|}\right)^{12}\!-2\left(\frac{\sigma_{\mathrm{eff}}}{|{\mathbf{r}}|}\right)^{6}+1\right]&\,\,\,|{\mathbf{r}}|\leq\sigma_{\mathrm{eff}},\\ \!0&\,\,\,|{\mathbf{r}}|>\sigma_{\mathrm{eff}},\end{array}\right.$$

$$\tilde{x}=\frac{x}{a},\tilde{y}=\frac{y}{a},\tilde{t}=\frac{D_{T}t}{a^2}.$$

$${\tilde{x}}_i(\tilde{t}+\Delta \tilde{t})={\tilde{x}}_i(\tilde{t})+[P{e}_s\cos {\theta }_i(\tilde{t})+{\tilde{f}}_i^x(\tilde{t})]\Delta \tilde{t}+\sqrt{2\Delta \tilde{t}}{R}_i^x$$

$${\tilde{y}}_i(\tilde{t}+\Delta \tilde{t})={\tilde{y}}_i(\tilde{t})+[P{e}_s\sin {\theta }_i(\tilde{t})+{\tilde{f}}_i^y(\tilde{t})]\Delta \tilde{t}+\sqrt{2\Delta \tilde{t}}{R}_i^y$$

$${\theta }_i(\tilde{t}+\Delta \tilde{t})={\theta }_i(\tilde{t})+\chi S_{\Omega }\Gamma \Delta \tilde{t}+\sqrt{2\chi \Delta \tilde{t}}{R}_i^{\theta },$$

$$\Gamma =\frac{|\Omega |}{D_R},$$

$$P{e}_s=\frac{a{V}_s}{D_T}=\frac{3V_s}{4D_{R}a}.$$

$$F_3\equiv F-2F_2\cos \theta ,$$

Full text

I Introduction

Active systems comprising self-propelled particles have emerged as a novel class of non-equilibrium systems, exhibiting unusual phases with diverse structural, dynamical and mechanical properties out of equilibrium Purcell (1977); Berg (2003); Lauga and Powers (2009); Ramaswamy (2010); Golestanian et al. (2011); Lauga (2011); Romanczuk et al. (2012); Marchetti et al. (2013); Yeomans et al. (2014); Elgeti et al. (2015); Cates and Tailleur (2015); Goldstein (2015); Lauga (2016); Zöttl and Stark (2016); Bechinger et al. (2016); Vicsek et al. (1995); Helbing (2001); Aditi Simha and Ramaswamy (2002); Shenoy et al. (2007); Kudrolli et al. (2008); Drocco et al. (2012). Typical examples of active particles include asymmetrically coated, or Janus, micro-/nano-particles that exhibit catalytic surface reactions and, hence, self-propulsion, in specific solvent mixtures Walther and Müller (2013, 2008); Jiang et al. (2010a); Perro et al. (2005); Jiang et al. (2010b); Valadares et al. (2010); Howse et al. (2007); Golestanian et al. (2005, 2007); Douglass et al. (2012); Volpe et al. (2014a); Bianchi et al. (2016); Buttinoni et al. (2013); Mano and Heller (2005). These synthetic particles have enabled detailed study of different aspects of active processes involving self-propulsion Walther and Müller (2013, 2008); Jiang et al. (2010a); Perro et al. (2005); Becker et al. (2003); Zöttl and Stark (2014); Volpe et al. (2011); Jiang et al. (2010b); Valadares et al. (2010); Howse et al. (2007); Golestanian et al. (2005); Najafi and Golestanian (2005); Golestanian et al. (2007); Najafi and Golestanian (2004); Douglass et al. (2012); Volpe et al. (2014a); Bianchi et al. (2016); Buttinoni et al. (2013); Mano and Heller (2005); Romanczuk and Schimansky-Geier (2011); Groß mann et al. (2013, 2012); Romanczuk et al. (2015); Uspal et al. (2015); Simmchen et al. (2016); Uspal et al. (2016). Self-propelled particles are also abundant in biology with examples including many of the living organisms, such as multiflagellate bacterium E. coli Berg (2003), biflagellate alga C. reinhardtii Goldstein (2015) and uniflagellated sperm cells Shack et al. (1974); Woolley (2003), which contain engines of their own to self-propel in fluid media and in the absence of external forces Purcell (1977); Berg (2003); Lauga and Powers (2009); Ramaswamy (2010); Golestanian et al. (2011); Lauga (2011); Romanczuk et al. (2012); Marchetti et al. (2013); Yeomans et al. (2014); Elgeti et al. (2015); Cates and Tailleur (2015); Goldstein (2015); Lauga (2016); Zöttl and Stark (2016); Bechinger et al. (2016).

In the most common models of self-propelled particles, it is assumed that the particles possess uniaxial symmetry (with examples ranging from Janus micro-/nano-spheres and rods to E. coli) and that the direction of the self-propulsion is along the axis of symmetry of the particle. Such active particles thus tend to orient along the self-propelling force and move in a straight line, although their translational and rotational motion may also be subject to passive and/or active noise. While passive noise (due, e.g., to thermal ambient fluctuations) is independent of the direction of motion, active noise (due, e.g., to fluctuations in the internal self-propulsion mechanism) correlates with that direction Romanczuk et al. (2012); Romanczuk and Schimansky-Geier (2011); Groß mann et al. (2013, 2012); Romanczuk et al. (2015).

In general, asymmetries in self-propulsion mechanism can generate a more complex behavior, including, e.g., situations with an unbalanced torque acting on the particle; in this case, the direction of motion and that of the active force are no longer aligned and the active particles tend to execute circular motion, with a characteristic radius inversely proportional to the magnitude of the torque (or, equivalently, their intrinsic angular velocity). Such chiral self-propelled particles have attracted mounting interest over the last few years van Teeffelen and Löwen (2008); Ai et al. (2015); Reichhardt and Olson Reichhardt (2013); Mijalkov and Volpe (2013); Breier et al. (2016); Ao et al. (2015, 2014); Romanczuk et al. (2012); Volpe and Volpe (2013); Volpe et al. (2014b); Nourhani et al. (2013, 2015); Li et al. (2014); Kümmel et al. (2013); Löwen (2016); Weber et al. (2011); Takagi et al. (2013); Friedrich and Jülicher (2008, 2009); Ebbens et al. (2010); Mirkovic et al. (2010); Boymelgreen et al. (2014); Marine et al. (2013); Wheat et al. (2010); Su et al. (2012); Di Leonardo et al. (2011). Depending on the sign of the torque, one can consider two types of left- and right-handed chiral self-propelled particles. Separation strategies for the two different types of active particles have been discussed in recent works Volpe and Volpe (2013); Volpe et al. (2014b); Ai et al. (2015). In addition, chiral self-propelled particles have shown the potential to be used as carries of non-chiral particles in ratchet channels Ai (2016). Another interesting, and yet unexplored, aspect of chiral self-propelled particles, to be addressed among other problems in this paper, is the non-equilibrium effective (possibly, depletion-type) interactions that an active bath of such particles can mediate between external, non-active, colloidal objects.

Depletion interactions have been studied extensively in the context of non-active macromolecular mixtures in equilibrium due both to their fundamental importance in the theory of colloidal stability Lekkerkerker and Tuinier (2011); Likos (2001) and their numerous technological applications as, e.g., in protein crystallization Jia et al. (2005) and phase separations of colloid/polymer mixtures Tuinier et al. (2003). In their most basic form, equilibrium depletion interactions can be established as effective, short-ranged, attractive forces between two hard spheres (inclusions) immersed in a bath of thermalized hard spheres of smaller size (depletants) Lekkerkerker and Tuinier (2011); Likos (2001). In this case, steric repulsions lead to exclusion of small spheres from the immediate vicinity of the inclusions, forming thin depletion layers (of thickness equal to small sphere radius) around them. The overlap between these layers at small surface separations creates a narrow intervening region of complete depletion between the inclusions and, hence, a net attractive force between them.

Depletion-type interactions have also been studied in non-equilibrium systems involving colloidal inclusions immersed in a bath of active particles Dzubiella et al. (2003); Angelani et al. (2011); Schwarz-Linek et al. (2012); Harder et al. (2014); Ni et al. (2015); Leite et al. (2016); Smallenburg and Löwen (2015); Angelani et al. (2011); Schwarz-Linek et al. (2012); Mallory and Cacciuto (2016); Tung et al. (2016), and their relation to Casimir-type forces has been discussed Ray et al. (2014). Studies of the phase behavior of suspensions of non-motile and motile E. coli, used as non-active and active colloids, respectively, in the presence of non-adsorbing polymers as depletant, show that activeness suppresses the depletion-driven phase separation of E. coli Schwarz-Linek et al. (2012). Harder et al. Harder et al. (2014) used Brownian Dynamics simulations to study the steady-state, effective two-body interactions mediated between large non-active colloids in a bath of small, non-chiral, active Brownian particles in two dimensions. For disk-shaped inclusions, they showed that the effective two-body force mediated between the inclusions can be strongly repulsive, in stark contrast with the non-active case, where, as noted above, depletion forces are attractive. The repulsive forces were shown to result from surface accumulation of active particles near the inclusions (a common trait for biological/colloidal active particles at confining boundaries Elgeti and Gompper (2016); Berke et al. (2008); Lauga et al. (2006); Spagnolie and Lauga (2012); Li et al. (2011); Li and Tang (2009); Li and Ardekani (2014); Hernandez-Ortiz et al. (2005); Schaar et al. (2015); Hu et al. (2015); Elgeti and Gompper (2009, 2013); Kantsler et al. (2005); Rusconi et al. (2014); Elgeti and Gompper (2015); Costanzo et al. (2012); Nash et al. (2010); Ezhilan and Saintillan (2015); Sabass and Seifert (2010)) and, in particular, in the intervening region between the inclusions as they are brought to small surface separations Harder et al. (2014).

Here, we first revisit the problem of effective two-body interactions between disk-shaped inclusions in an active bath within a commonly used, two-dimensional model of self-propelled particles consistent with the one studied in Ref. Harder et al. (2014). Using Brownian Dynamics simulations, and by considering both chiral and non-chiral active baths, we show how the salient, and previously unexplored, aspects of the spatial distribution of active particles, including formation of circular layers (or ‘rings’) of active particles around the inclusions, determine the qualitative nature (e.g., repulsive versus attractive) as well as quantitative features (e.g., non-monotonic behavior with distance) of the two-body interaction force profiles. The sequential overlaps of active-particle rings, which in a way create soft, repulsive, ‘shoulders’ Glaser et al. (2007); Ziherl and Kamien (2011) around the inclusions, generate the distinct features of the force profiles as the inclusions are brought to small surface separations. Active-particle chirality leads to suppression of particle rings and partial depletion of active particles from the intervening region and the farther proximity of the inclusions. We thus show how particle chirality can engender a crossover in the parameter space between two characteristically different limits of pure active repulsion and pure depletion attraction; hence, linking them as two complementary limits for the system behavior. We also argue that, through their spatial structuring around the inclusions, active particles can also cause significant non-pairwise effects in the interaction forces between the inclusions, for which we provide direct evidence by investigating these effects on the leading three-body level.

The paper is organized as follows: In Section II, we introduce our model and details of our simulations. The distribution of non-chiral and chiral active particles and their influence on the two- and three-body interactions between the inclusions are discussed in Sections III and IV, respectively. The paper is concluded in Section V.

II Model and methods

Our model consists of two or three identical, non-active and nearly hard colloidal inclusions of radius $a_{c}$ placed at fixed positions in a bath of identical, self-propelled, Brownian particles of smaller radius $a<a_{c}$ and area fraction $\phi$ in two spatial dimensions. Such two-dimensional models have widely been used in the literature as they enable computationally efficient simulations, while capturing key features of active systems Lauga and Powers (2009); Marchetti et al. (2013); Yeomans et al. (2014); Elgeti et al. (2015); Lauga (2016); Golestanian et al. (2011); Ramaswamy (2010); Cates and Tailleur (2015); Romanczuk et al. (2012); Zöttl and Stark (2016); Bechinger et al. (2016).

Active particles are assumed to self-propel at constant speed $V_{s}$ along their preferred direction of motion ${\mathbf{n}}=(\cos\theta,\sin\theta)$. The instantaneous configuration of active particles (labelled by $i$) can thus be described by the set of position vectors ${\mathbf{r}}_{i}(t)=(x_{i}(t),y_{i}(t))$, and orientation angles $\theta_{i}(t)$ (see Fig. 1). These degrees of freedom are assumed to evolve in time according to the overdamped Langevin equations (see, e.g., Refs. Marchetti et al. (2013); Elgeti et al. (2015); Romanczuk et al. (2012); Volpe et al. (2014b))

[TABLE]

where $D_{T}$ and $D_{R}$ are the (bare) translational and rotational diffusion coefficients of active particles, and $\Omega$ their intrinsic angular velocity, whose sign $S_{\Omega}=\pm 1$ gives the chirality type, with positive (negative) sign corresponding to counterclockwise (clockwise) angular velocity. Here, ${\bm{\eta}}_{i}(t)$ and $\zeta_{i}(t)$ are independent, white, Gaussian translational and rotational noises, respectively. They have zero mean, $\langle{\eta}_{i}^{\alpha}(t)\rangle=\langle\zeta_{i}(t)\rangle=0$, and two-point time correlations $\langle{\eta}_{i}^{\alpha}(t){\eta}_{j}^{\beta}(t^{\prime})\rangle=\delta_{ij}\delta_{\alpha\beta}\delta(t-t^{\prime})$ and $\langle\zeta_{i}(t)\zeta_{j}(t^{\prime})\rangle=\delta(t-t^{\prime})$, with $i,j$ denoting the active particle labels and $\alpha,\beta$ the two Cartesian directions $x,y$.

We assume the Einstein-Smoluchowski-Sutherland relation between translational mobility and diffusion coefficients $D_{T}=\mu_{T}k_{\mathrm{B}}T$, where $k_{\mathrm{B}}$ is the Boltzmann constant and $T$ the ambient temperature. This assumption and the choice of thermal (passive) noises are made to ensure that the system approaches its appropriate Boltzmann-weighted equilibrium determined by the potential $U(\{{\mathbf{r}}_{j}\})$ in the steady state, when the active self-propulsion is switched off. This enables direct comparisons between non-equilibrium effective interactions generated in an active bath and their equilibrium counterpart when bath particles are non-active. Our model thus excludes athermal noises associated with fluctuations in the self-propulsion and angular velocities Romanczuk et al. (2012); Romanczuk and Schimansky-Geier (2011); Groß mann et al. (2013, 2012); Romanczuk et al. (2015).

The potential energy $U$ in Eq. (1) gives the sum of pair potentials between all particles (including active particles and inclusions) in the system. The pair potential is taken in a Weeks-Chandler-Andersen (WCA) form as

[TABLE]

where $|{\mathbf{r}}|$ is the center-to-center distance between the considered pair of particles with $\sigma_{\mathrm{eff}}=2a$, when the particles considered are both of active particles, $\sigma_{\mathrm{eff}}=2a_{c}$, when the pair are both of inclusions, and $\sigma_{\mathrm{eff}}=a+a_{c}$, when one of the particles is an active particle and the other one is an inclusion. In all these cases, the interaction energy strength is taken to be the same and equal to $\epsilon$.

II.1 Simulations: Methods and parameters

In order to proceed, we use a dimensionless representation by rescaling the units of length and time as

[TABLE]

Equations (1) and (2) can then be solved numerically using Brownian Dynamics methods by rewriting them in dimensionless and discrete form for time evolution over a sufficiently small time step $\Delta\tilde{t}$ as

[TABLE]

where $\tilde{f}^{x}_{i}=-\partial\tilde{U}/\partial\tilde{x}_{i}$ and $\tilde{f}^{y}_{i}=-\partial\tilde{U}/\partial\tilde{y}_{i}$ are the Cartesian components of the dimensionless force derived from the rescaled potential $\tilde{U}=U/(k_{\mathrm{B}}T)$, and $R^{x}_{i}$, $R^{y}_{i}$ and $R^{\theta}_{i}$ are independent Gaussian-distributed random numbers with zero mean and unit variance. In Eq. (7), we have defined $\chi=a^{2}D_{R}/D_{T}$, which, using the fact that in the low-Reynolds-number (Stokes) regime the translational and rotational diffusion coefficients for no-slip spherical particles satisfy the relation $D_{R}=3D_{T}/4a^{2}$ Happel and Brenner (1983), can be fixed conventionally as $\chi=3/4$. We have also defined the chirality strength parameter as the dimensionless magnitude of the (intrinsic) particle angular velocity,

[TABLE]

and the so-called swim Péclet number as

[TABLE]

This latter quantity gives the ratio of the characteristic timescale of the translational (or rotational) diffusion of the active particles, $a^{2}/D_{T}$ (or $1/D_{R}$), to the characteristic timescale of their active self-propulsion (swim), $a/V_{s}$. Note also that, in the present context, the sign of the angular velocity of active particles is irrelevant and only its magnitude, $|\Omega|$, will play a role.

In rescaled units, the system is described by the size ratio $a_{c}/a$, the rescaled area fraction $\phi a^{2}$, the swim Péclet number $Pe_{s}$, the chirality strength parameter $\Gamma$, and the rescaled center-to-center distance between adjacent inclusions $\tilde{d}=d/a$ or, equivalently, their rescaled surface-to-surface distance $\tilde{\Delta}=\Delta/a$ with $\Delta=d-2a_{c}$ (see Fig. 1). Our focus will be on the roles of self-propulsion and chirality strengths as our main control parameters; we will then fix other parameter values as $\epsilon/(k_{\textup{B}}T)=10$, $a_{c}/a=5$ and $\phi a^{2}=0.1$, consistent with those chosen in Ref. Harder et al. (2014) (the effects due to varying size ratio and area fraction for active particles and inclusions will be discussed in the more general context of active mixtures elsewhere Zaeifi Yamchi and Naji (2017)). In our simulations, $Pe_{s}$ is increased from 0 up to around 45 and $\Gamma$ from 0 up to around 19, spanning a wide range of experimentally accessible, actual parameter values; e.g., the cases with $Pe_{s}=22.8$ and $\Gamma=19.1$, which will be considered later, can be mapped to $a=1\,\mu{\mathrm{m}}$, $a_{c}=5\,\mu{\mathrm{m}}$, $V_{s}\simeq 5\,\mu{\mathrm{m}}\cdot{\mathrm{s}}^{-1}$, $D_{T}\simeq 0.22\,\mu{\mathrm{m}}^{2}\cdot{\mathrm{s}}^{-1}$, $D_{R}\simeq 0.16\,{\mathrm{s}}^{-1}$ and $|\Omega|\simeq\pi\,\,{\mathrm{s}}^{-1}$ in an aqueous medium with shear viscosity $\eta=0.001\,{\mathrm{Pa\cdot s}}$ Volpe et al. (2014b). It is also worth mentioning that active particles can exhibit a wide range of chirality strengths in experiments; examples include Janus doublets ($\Gamma\simeq 14-27$) Ebbens et al. (2010), self-propelled rods ($\Gamma\simeq 8$) Takagi et al. (2013) and active L-shaped particles ($\Gamma\simeq 200$) Kümmel et al. (2013). Therefore, while our focus will primarily be on the generic aspects of the considered model, the parameters are varied within a realistic range of values.

Our simulations typically run for $10^{6}-10^{8}$ time steps (with typical time-step size chosen as $\Delta\tilde{t}\simeq 10^{-4}$) with $10^{6}$ steps used initially for relaxation purposes and the rest used for computing averaged quantities; the results are further averaged over about 20 statistically independent samples after the system has reached a steady state. For the most part, the simulations are performed using one hundred active bath particles distributed in random initial positions in a square box with periodic boundary conditions, whose lateral size is adjusted according to the given area fraction (test simulations for systems with up to 300-400 particles show only quantitative differences of at most 10%, while qualitative aspects of our results remain unchanged). In the plots shown later, computed error bars are typically smaller than the size of symbols.

II.2 Forces acting on inclusions

The net force acting on a fixed inclusion (in the absence or presence of other fixed inclusions in the bath) follows from the averaged sum of instantaneous forces exerted on it by bath particle collisions. To introduce some of the terminology that we shall use later, we describe here the method of calculating the two-body forces to be analyzed in Section III. The method of evaluating three-body forces will be discussed in Section IV.

Since the active bath is homogeneous and isotropic, the net force acting on a single inclusion due to active particle collisions in the bath turns out to be zero on average (within our simulation margin of error). Hence, in the case of two fixed inclusions, the net force, which is found to be non-zero and to act on each of the inclusions with equal magnitude and in opposite directions along the $x$ (center-to-center) axis, can be interpreted as the effective, two-body, interaction force mediated between them by the bath. We conventionally denote the rescaled force acting on the inclusion placed on the right in Fig. 1 as $\tilde{\mathbf{F}}_{2}=\tilde{F}_{2}\hat{\mathbf{x}}$; thus, an attractive (repulsive) interaction force is represented by a negative (positive) force amplitude $\tilde{F}_{2}$. We calculate this interaction force using $\tilde{\mathbf{F}}_{2}=\sum_{i}\langle\tilde{\mathbf{f}}_{i}\rangle$, where the brackets $\langle\cdots\rangle$ denote the average taken over many simulated configurations and $\tilde{\mathbf{f}_{i}}=(\tilde{f}_{i}^{x},\tilde{f}_{i}^{y})$ denote instantaneous force components imparted on the mentioned inclusion from the $i$th active particle; these force components follow from the respective WCA interaction potentials discussed in Section II.1.

III Effective two-body interactions

III.1 Non-chiral active bath

We start our analysis of effective two-body interactions by considering a system of two fixed inclusions in a bath of non-chiral active particles. Figure 2 shows the results for the rescaled two-body force amplitude $\tilde{F}_{2}$, as defined in Section II.2, as a function of the rescaled surface-to-surface distance between the inclusions, $\tilde{\Delta}$. We show the results for a few different values of the swim Péclet number, $Pe_{s}$, including the non-active case with $Pe_{s}=0$. For the sake of demonstration, the results are divided by $Pe_{s}+1$ to display the quantity $\hat{F}_{2}\equiv\tilde{F}_{2}/(Pe_{s}+1)$.

As seen in the figure, the two-body interaction force turns from its typical non-active (equilibrium) form, representing a relatively short-ranged, attractive depletion force (black filled circles) to a much stronger repulsive force with a longer range of action (of a few active-particle radii in surface separations) as $Pe_{s}$, or the self-propulsion strength, is increased. This is in agreement with the findings in Ref. Harder et al. (2014). Our data, however, resolve the behavior of the force in more detail around the first peak and reveals the presence of a secondary hump (which appears more like a plateau region for the give parameter values in the figure) at larger separations, when the swim Péclet number is sufficiently large. The first peak occurs at rescaled surface-to-surface distances $\tilde{\Delta}\lesssim 2$ (or, equivalently, $\Delta\lesssim 2a$) and its location shifts to smaller values as $Pe_{s}$ is increased. The somewhat peculiar non-monotonic behavior of the two-body force profiles and the presence of a secondary hump, which is not reported in Ref. Harder et al. (2014), can be understood by looking more closely at the steady-state distribution of active particles around the inclusions, giving insight into the mechanism underlying the salient features of the force profiles.

We define $g(\tilde{r})$ as the probability of finding an active particle at a distance $\tilde{r}$ from the center of one of the inclusions (e.g., the inclusion on the right in Fig. 3, inset) along the center-to-center axis in the outward direction pointing toward the bulk (this particular direction is chosen as the influence of the other inclusion on $g(\tilde{r})$ along this direction is negligible). For the sake of brevity, we refer to $g(\tilde{r})$ simply as the active particle-inclusion pair distribution function, or the pdf. As seen in Fig. 3, this quantity exhibits two peaks, indicating a layered structure for the spatial distribution of active particles at sufficiently large $Pe_{s}$. The first peak in $g(\tilde{r})$ indicates a high-density layer, or the primary ring, of active particles at a close distance from the inclusion surface; this distance tends to roughly one active-particle radius as $Pe_{s}$ is increased, while the second peak in $g(\tilde{r})$, indicating a less populated secondary ring of active particles, is found at a larger distance (specifically, for $Pe_{s}=45.8$, the primary and secondary rings are found at separations $\tilde{r}_{1}-a_{c}/a\simeq 0.85$ and $\tilde{r}_{2}-a_{c}/a\simeq 2.8$ from the inclusion surface, respectively). As noted previously, active particles are known to exhibit strong accumulation near surface boundaries Elgeti and Gompper (2016); Berke et al. (2008); Lauga et al. (2006); Spagnolie and Lauga (2012); Li et al. (2011); Li and Tang (2009); Li and Ardekani (2014); Hernandez-Ortiz et al. (2005); Schaar et al. (2015); Hu et al. (2015); Elgeti and Gompper (2009, 2013); Kantsler et al. (2005); Rusconi et al. (2014); Elgeti and Gompper (2015); Costanzo et al. (2012); Nash et al. (2010); Ezhilan and Saintillan (2015); Sabass and Seifert (2010); this tendency is opposed, in the present context, by the steric repulsions between the particles, giving rise to their mentioned layer ordering in the high density regions around the inclusions.

The angular dependence of the simulated active particle distribution is shown in Fig. 4, which complements the information obtained from $g(\tilde{r})$ by presenting the two-dimensional view of the active-particle density profiles around and in between the two inclusions. In this figure, the circular areas excluded by the inclusions themselves are shown as green disks. The thin (dark blue) circular areas (of thickness equal to one active-particle radius) that are attached to the inclusion disks are regions from which the active particles are further excluded due to their finite radius. The white areas with high active-particle density correspond to the primary rings around the inclusions and also the regions of overlap between the secondary rings, appearing themselves in yellow. It is worth mentioning here that each inclusion with its two active-particle rings can be thought of as a single colloidal unit with a hard core and a soft repulsive shoulder. Qualitatively similar scenarios, incorporating hard-core/soft-shoulder pair potentials, have been considered in the equilibrium context of non-active colloids and shown to result in a diverse phase diagram Glaser et al. (2007); Ziherl and Kamien (2011). Those findings may, however, be inapplicable to the present case, in which the circular soft-ring zones due to active-particle layering form only in a non-equilibrium steady state.

The non-monotonic behavior of the two-body interaction force profiles, and their increased range of action at elevated $Pe_{s}$ (Fig. 2), can be understood based on sequential overlaps (or intersections) occurring between the active-particle rings as the two inclusions are brought together. This is shown schematically in Fig. 5, where we re-plot the interaction force profile in the exemplary case of $Pe_{s}=45.8$ together with schematic illustrations of two apposed complexes comprising a central inclusion and its two active-particle rings at four indicative surface-to-surface distances. At large separations, the two complexes are decoupled and they begin to interact with a non-vanishing effective force only when their secondary rings come into contact; this corresponds to configuration I as indicated on the graph (occurring at $\tilde{\Delta}\simeq 5.6$). As $\tilde{\Delta}$ is further decreased, the effective repulsion between the inclusions is expected to increase as the secondary ring associated with one inclusions intersects the primary ring of another inclusion (configuration II; $\tilde{\Delta}\simeq 3.65$). This results from an increased build-up of active particle density at regions of ring intersection between the inclusions (see also the left panel in Fig. 4) and, as a result, increased collision forces from active particles pushing the inclusions apart. After this point, one can expect only a weak change in the magnitude of the repulsion (leading to the secondary hump or the plateau-like region in the force profile) at about the separation, where the secondary ring of one inclusion comes into contact with the surface of another inclusion (configuration III; $\tilde{\Delta}\simeq 2.8$). The repulsive force is expected to sharply increase as $\tilde{\Delta}$ is squeezed down to smaller separations, where the primary ring of one inclusion also comes into contact with the surface of another inclusion (configuration IV; $\tilde{\Delta}\simeq 0.85$), creating even larger collision forces from active particles pushing the inclusions apart. This will in fact be the global maximum in the effective force profile since active particles are excluded from the narrow intervening region between the inclusions as $\tilde{\Delta}$ is decreased even further.

III.2 Chiral active bath

We now turn to the case of chiral active particles in a system with two fixed colloidal inclusions. Chirality of the active particles adds a new characteristic length scale to the problem $R_{\Omega}=V_{s}/|\Omega|$ (see Fig. 1), representing the characteristic radius of a typical circular arc traversed by the particles (note that, in the present context, active particles do not move on perfect circular trajectories due to collisions with other particles and the rotational Brownian noise). Chirality effects are expected to dominate when $R_{\Omega}$ becomes smaller than the swim (run) length $\ell_{\textrm{run}}=V_{s}/D_{R}$ Ao et al. (2015, 2014), or equivalently, when the timescale associate with active-particle chirality $|\Omega|^{-1}$ is smaller than that of its rotational diffusion $1/D_{R}$, as, otherwise, the chirality effects are masked by the rotational diffusion process. These conditions translate to having a large chirality strength parameter $\Gamma\gg 1$. To gain a more systematic understanding of the chirality effects in the present context, we run Brownian Dynamics simulations similar to those reported in Section III.1 but assuming here that $\Gamma$ is finite (ranging from 0 up to around 19, corresponding to angular velocity magnitudes $|\Omega|$ from 0 up to around $\pi\,{\mathrm{s}}^{-1}$, when $D_{R}\simeq 0.16\,{\mathrm{s}}^{-1}$ Volpe et al. (2014b)).

As seen in Fig. 6, the peaks in the active particle-inclusion pdf are largely suppressed upon increasing the chirality strength parameter. Interestingly, however, we find an extended range of separations (Fig. 6, inset) over which the pdf drops below its bulk value, representing a region of partial depletion right after the primary ring (the width of this depletion region is of the order of the rescaled radius of curvature, i.e., $R_{\Omega}/a=4Pe_{s}/(3\Gamma)\simeq 3.2$ for the parameter values in the inset). The density map of Fig. 7 gives a clearer view of the primary (yellow) rings at the vicinity of the inclusions and the relatively thick, circular, partial depletion zones (in darker hazy colors) at their farther proximity; the overlap between these latter zones gives a narrow intervening region of strong particle depletion (in dark blue) between the inclusions, where their primary rings disappear, too.

As a consequence of the changes in the spatial distribution of active particles, the effective two-body force between the inclusions varies drastically as $\Gamma$ is increased; see Fig. 8. The repulsive small-distance peak in the effective force is suppressed and so is the force magnitude over the plateau-like region at larger surface-to-surface distances $\tilde{\Delta}$. The plateau-like region is also moved to smaller separations and, at sufficiently large chirality strength (data for $\Gamma\geq 6.1$ in the graph), it is followed by a relatively sharp drop to a larger-distance region, where the interaction force between the inclusions becomes attractive (negative). A closer view of this behavior is shown in the inset of Fig. 8.

The effective attraction found at sufficiently large $\Gamma$ thus appears to display some of the qualitative features of non-active (equilibrium) depletion attractions Lekkerkerker and Tuinier (2011). In fact, one can show using systematic arguments (albeit in a different context Jamali and Naji (2017)) that, when the limit $\Gamma\rightarrow\infty$ is taken, the steady-state properties of a system of self-propelled Brownian particles reduce to their corresponding equilibrium values (as obtained by setting the self-propulsion strength equal to zero). It should, however, be noted that, within our model, chirality effects will be seen only when the particles are active. Also, the resulting effective attraction at finite, even though large chirality strength, still shows significant deviations from its non-active limiting form; these deviations can be traced back to the structure of particle depletion zones, which remain different from their limiting non-active form (thin, dark blue, circular areas attached to the inclusions Fig. 7). Hence, in addition to the fact that the active, chirality-induced, depletion attraction reported here represents a qualitatively distinct state of non-equilibrium, it also exhibits a range of action a few times larger than its characteristically short-ranged, equilibrium counterpart Lekkerkerker and Tuinier (2011) (compare data in Fig. 8, inset, with $Pe_{s}=0$ in Fig. 2).

Our data also shows a non-monotonic trend in the behavior of the attractive force as a function of $\Gamma$ (Fig. 8, inset): Increasing $\Gamma$, after the force mediated between the inclusions turns attractive, slightly decreases the force magnitude and, as noted above, gradually shifts the attractive region to smaller $\tilde{\Delta}$ as $\Gamma$ is further increased.

IV Three-body interactions

IV.1 Non-chiral active bath

In non-active systems, depletion interactions between colloidal inclusions in a bath of (smaller) particles are generally known to be non-pairwise, especially, at relatively high volume fractions, where such (typically short-ranged) interactions become important Lekkerkerker and Tuinier (2011). Another important factor is the thickness of the depletion layer around the inclusions: As illustrated in Ref. Lekkerkerker and Tuinier (2011) (see Fig. 3.8 therein), non-pairwise effects become important when depletion layers are thick enough (relative to the inclusion size) to exhibit multiple overlap regions as more than two inclusions are brought to small surface separations from one another; in this case, the effective (net) interaction forces acting on the inclusions will involve significant multi-body contributions. The lowest order of multiple overlaps occurs in the case of three inclusions placed at equal surface-to-surface distances creating an equilateral-triangle formation Lekkerkerker and Tuinier (2011); this is indeed the case that we shall consider in what follows by first assuming that the active particles are non-chiral ($\Gamma=0$).

The two-dimensional active-particle density profile around and in between the inclusions in this case is shown in Fig. 9 by taking the surface-to-surface distance of the inclusions as $\tilde{\Delta}=1.5$ and the Péclet number of the active particles as $Pe_{s}=45.8$. The interesting aspect of the problem in this case is that the relatively thick shoulders comprising two rings of active particles enable large multiple-overlap areas over a range of surface separations comparable to a few active-particle radii. Hence, in analogy with the multiple-overlap mechanism described above Lekkerkerker and Tuinier (2011), we may as well anticipate significant non-pairwise three-body contributions to occur due to such overlaps in the present context.

The effective net force acting on individual inclusions in a triangular configuration similar to the one shown in Fig. 9 can be evaluated as a function of their (equal) surface-to-surface distance $\tilde{\Delta}$, and the swim Péclet number $Pe_{s}$. The net forces experienced by the inclusions are equal in size and are aligned with the bisectors of the angles formed by the triangular configuration of the inclusion centers. The results, in rescaled units and divided also by $Pe_{s}+1$ (see Section III.1), are shown in Fig. 10, where negative (positive) values represent forces pointing toward (away from) the center of the triangle, respectively. Comparing this figure with Fig. 2, the net force acting on a single inclusion is found to show roughly similar qualitative features in both the two-inclusion and the three-inclusion configurations. There are, however, remarkable differences between the results in the two cases at a given $Pe_{s}$: The magnitude of the net force in the three-inclusion configuration shows an overall increase by a factor of more than two; the secondary hump (which appeared more like a plateau region in the force profiles of the two-inclusion system) is more pronounced in the three-inclusion setting; finally, the net force in the non-active case ($Pe_{s}=0$) displays a repulsive hump in the three-inclusion system, one that is absent in the two-inclusion case (compare Figs. 10 and 2). These differences all point to the presence of sizable three-body effects in the system as we shall discuss below.

If we denote the effective net force acting on a given inclusion by $F$ and the two-body force acting on this reference inclusion by either of its two neighboring inclusions by $F_{2}$, then the three-body force follows as

[TABLE]

where, in the equilateral-triangle arrangement considered here $\theta=\pi/6$. We have also used the fact that, by construction, the two-body force $F_{2}$ acts along the line connecting the centers of the inclusions in question and that its values can be taken from our results in Fig. 2.

The resulting three-body force contribution is shown in rescaled units (divided also by $Pe_{s}+1$) in Fig. 11. (Note that, in our simulations, the set of values for the variable $\tilde{\Delta}$ at which the net force $\hat{F}$, Fig. 10, and the two-body force $\hat{F}_{2}$, Fig. 2, are calculated, do not necessarily coincide and, thus, $\hat{F}_{3}$ is obtained by cubic-spline interpolation of the corresponding data sets in those two figures; hence, reported in Fig. 11 are continuous and accurately representing curves rather than individual symbols.)

The simulated three-body force exhibits rapid, sizable variations in sign and magnitude as a function of $\tilde{\Delta}$ and $Pe_{s}$, in stark contrast with the behavior found in the case of the two-body interaction force in Fig. 2. The three-body force becomes strongly repulsive at small separations $\tilde{\Delta}\lesssim 1$ with magnitudes exceeding those of the two-body force; e.g., in this regime, we find the force ratio $|\hat{F}_{3}|/|\hat{F}_{2}|\lesssim 2$ for $Pe_{s}=22.8$ and $|\hat{F}_{3}|/|\hat{F}_{2}|\lesssim 1.3$ for $Pe_{s}=45.8$. Our data clearly indicates a non-monotonic behavior for the dependence of $|\hat{F}_{3}|$ on $Pe_{s}$ at small separations between the inclusions.

The three-body force profiles as a function of the surface-to-surface distance between the inclusions (Fig. 11) also appear to show two major maxima, whose locations appear to be roughly consistent with the locations of the primary and secondary active-particle rings around the inclusions (for both $Pe_{s}=22.8$ and 45.8, these two maxima are found in the intervals $\tilde{\Delta}\lesssim 1$ and $2\lesssim\tilde{\Delta}\lesssim 3$). This correspondence is, however, not precise and the relation between the variations in sign and magnitude of the three-body force and the ring-structure around the inclusions remains to be understood.

Finally, our data for non-active bath particles ($Pe_{s}=0$; Fig. 11) indicate that the repulsive hump, mentioned in our discussion of the net force in Fig. 10, is entirely due to three-body effects. In the two-inclusion system, as the surface separation of inclusions is decreased down to the particle diameter ($\tilde{\Delta}\simeq 2$), particles are freely depleted from the narrow inter-surface gap between the inclusions, leading to an attractive depletion force without encountering any possible repulsive barrier on the way (Fig. 2). This is not the case for three inclusions being brought to small surface separations, where the wedge-shaped, central void formed in the intervening region between the inclusions can accommodate a finite number of bath particles; these entrapped particles produce a counteracting outward force barrier, before it is overcome and the entrapped particles are squeezed out upon further packing of the three inclusions together.

IV.2 Chiral active bath

We proceed by considering the three-body interactions in the case of chiral active particles. The density map of active particles in the three-inclusion system for relatively high values of the swim Péclet number and the chirality strength parameter ($Pe_{s}=45.8$ and $\Gamma=19.1$) again shows a primary (yellow) ring of active particles in close vicinity and thick, circular, zones (in darker hazy colors) of partial depletion at the farther proximity of the inclusions, in addition to a central region (in dark blue) of strong depletion. Hence, the effective net force acting on individual inclusions is expected to weaken in strength as $\Gamma$ is increased as confirmed by our data in Fig. 13.

Comparing the net force acting on a single inclusion in the two- and three-inclusion configurations (Figs. 8 and 13), we reach similar conclusions as we did in the non-chiral case in Section IV.1; that is, the magnitude of the net force in the three-inclusion configuration shows an overall increase by a factor of more than two and the secondary hump in the force profile becomes more pronounced, pointing again to significant non-pairwise three-body effects in the system.

The three-body force, which contributes to the effective net force (Fig. 13) acting on individual inclusions in a bath of chiral active particles, can be calculated using Eq. (10). The results, in rescaled units and divided also by $Pe_{s}+1$, are shown in Fig. 14 for $Pe_{s}=45.8$ and different values of $\Gamma$. Note that the black solid curve in this figure coincides with the blue dashed curve in Fig. 11, representing the non-chiral case ($\Gamma=0$) with its characteristic variations in sign and magnitude as discussed in the previous section. As the chirality strength parameter is increased, these features are suppressed and the three-body force profiles become increasingly more smooth. Thus, particle chirality consistently also weakens the non-pairwise three-body contribution to the effective net force in the three-inclusion configuration. The three-body force becomes attractive over an extended of range surface separations (up to several active-particle radii; see Fig. 14, inset) for $\Gamma\geq 6.1$, where it can, therefore, be considered as an active three-body depletion attraction induced by the particle chirality.

V Conclusion and Discussion

In this paper, we use Brownian Dynamics simulations to study effective interactions generated by an active bath of non-chiral or chiral active particles between large, non-active, colloidal inclusions in the bath. The main emphasis of our study is on the role of the spatial distributions of active particles around the inclusions, and the way they influence the qualitative nature (such as repulsive versus attractive) as well as the quantitative features (such as non-monotonic behavior with distance) of the two- and three-body interactions that the active bath particles mediate between the inclusions.

In the case of a non-chiral active bath, the effective two-body force profiles are found to exhibit a non-monotonic behavior with finer details than previously reported Harder et al. (2014). We provide new insight into the behavior of the force profiles by establishing the formation, and the sequential overlaps, of rings of active particles around the inclusions as the root cause for the salient features of the force profiles that include a primary maximum and a secondary hump (or plateau-like region). The active-particle rings, on the other hand, create relatively thick, soft shoulders around the hard-core inclusions, enabling them to generate multiple-overlap regions as the inclusions are brought together; hence, leading to non-pairwise interactions between the inclusions. This is in line with the general arguments used in the equilibrium context of non-active colloids in a bath of small non-active depletants, where multiple overlaps between depletion layers (occurring when the depletion layers are thick enough relative to the inclusion size) are described as the mechanism underlying non-pairwise many-body interactions between the inclusions Lekkerkerker and Tuinier (2011). The lowest order of such multiple overlaps occurs in the case of three inclusions at equal surface-to-surface distances in an equilateral-triangle formation (Fig. 3.8 in Ref. Lekkerkerker and Tuinier (2011)), which is also the case considered in our analysis of three-body forces in an active bath. Our results thus show that, while an active bath can mediate strong repulsive two-body forces between inclusions, it can also produce comparably strong three-body forces with a more complex profile, displaying distinct repulsive and attractive regimes. Our analyses also indicate that multiple overlaps may be used as a measure for the significance of non-pairwise many-body contributions to the lowest order in the present non-equilibrium context with active particles. This point, however, remains to be confirmed by more systematic arguments.

In the case of chiral active bath particles, the spatial distribution of active particles and the effective two- and three-body interactions mediated by them between the inclusions are found to be strongly dependent on the particle chirality: The layered structure of active particles is largely suppressed at elevated chirality strengths, giving rise to strong depletion of active particles from the narrow intervening region between the inclusions, and also their partial depletion from relatively thick, circular, zones further away from the inclusions. Consequently, the repulsive interaction forces mediated between the inclusions are weakened and, eventually, and for both the two- and the three-body forces, turn to relatively long-ranged, active and chirality-induced, depletion-type attractions at large enough particle chirality strength.

The results reported in our work signify the important role of non-equilibrium active-particle layering and the resulting many-body interactions in active mixtures and also provide insight into the role of chirality. While in the equilibrium (non-active) context, standard theories (such as the free-volume and/or scaled particle theories in the case of hard spheres Lekkerkerker and Tuinier (2011); Likos (2001)) can be used to determine the phase behavior of the system, systematic approaches to study many-body effects and the phase behavior in active systems have been developed only quite recently (see Ref. Cates and Tailleur (2015) for a recent review). The importance of contributions from multi-body interactions in active particle/colloid mixtures can, nevertheless, be determined by resorting to numerical simulations. An efficient strategy is to develop direct, quantitative comparisons between implicit-active-particle simulations (in which only colloidal inclusions are explicitly modeled and they are then assumed to interact through effective forces mediated by active particles, such as those reported here, up to a given order in multi-body interactions) and explicit-active-particle simulations (in which both active particles and colloidal inclusions are modeled explicitly) Zaeifi Yamchi and Naji (2017).

On the other hand, with recent progress in design and manipulation of synthetic Janus particles Walther and Müller (2013, 2008); Jiang et al. (2010a); Perro et al. (2005); Jiang et al. (2010b); Valadares et al. (2010); Howse et al. (2007); Golestanian et al. (2005, 2007); Douglass et al. (2012); Volpe et al. (2014a); Bianchi et al. (2016); Buttinoni et al. (2013); Mano and Heller (2005), it is possible to develop active particles with different specifics (such as size, shape, and chirality); hence, exploring the non-equilibrium phase behavior of active particle/colloid mixtures by including these and other factors (such as inter-particle hydrodynamic coupling; see, e.g., Refs. Baskaran and Marchetti (2009); Evans et al. (2011); Molina et al. (2013); Zöttl and Stark (2014); Berke et al. (2008); Lauga et al. (2006); Spagnolie and Lauga (2012); Li and Ardekani (2014); Nash et al. (2010); Schaar et al. (2015); Elgeti and Gompper (2009); Hu et al. (2015); Hernandez-Ortiz et al. (2005); Elgeti and Gompper (2016)) appear as interesting potential directions for future research, where the question of effective interactions mediated by active constituents between inclusions is expected to play a central role.

Other interesting problems that can be studied within the present context include the role of active noise (represented by fluctuating self-propulsion and angular velocities) in the dynamics of active particles Romanczuk et al. (2012); Romanczuk and Schimansky-Geier (2011); Groß mann et al. (2013, 2012); Romanczuk et al. (2015); these effects can result in features distinct from those obtained with passive, thermal, noise considered here.

VI Acknowledgements

We thank the High Performance Computing Center (School of Computer Science, IPM) for computation time and Iran Science Elites Federation (ISEF) for partial support. A.N. acknowledges partial (Regular Associate) support from The Abdus Salam International Centre for Theoretical Physics (Trieste, Italy). We acknowledge useful discussions with M. Sebtosheikh and T. Jamali.