Noncentral forces mediated between two inclusions in a bath of active Brownian rods.

Using Brownian Dynamics simulations, we study effective interactions mediated between two identical and impermeable disks (inclusions) immersed in a bath of identical, active (self-propelled), Brownian rods in two spatial dimensions, by assuming that the self-propulsion axis of the rods may generally deviate from their longitudinal axis. When the self-propulsion is transverse (perpendicular to the rod axis), the accumulation of active rods around the inclusions is significantly enhanced, causing a more expansive steric layering (ring formation) of the rods around the inclusions, as compared with the reference case of longitudinally self-propelling rods. As a result, the transversally self-propelling rods also mediate a significantly longer ranged effective interaction between the inclusions. The bath-mediated interaction arises due to the overlaps between the active-rod rings formed around the inclusions, as they are brought into small separations. When the self-propulsion axis is tilted relative to the rod axis, we find an asymmetric imbalance of active-rod accumulation around the inclusion dimer. This leads to a noncentral interaction, featuring an anti-parallel pair of transverse force components and, hence, a bath-mediated torque on the dimer.

Introduction

Active particles are characterized by their ability to take up ambient free energy and convert it to self-propelled motion[1,2]. They constitute a large class of particles, including molecular motors[3-6] that are powered by the hydrolysis of any available adenosine-tri-phosphate molecules, motile microorganisms such as bacteria and sperm cells[7-12] that self-propel via cellular organelles such as cilia and flagella, and artificial nano and microswimmers[13-23] such as Janus particles and functionalized colloids that undergo active motions due to catalyzing surface reactions or external stimuli such as light. The self-propelled motion drives any system of particles suspended in a base fluid out of thermodynamic equilibrium, engendering a host of intriguing effects, including collective dynamic self-organization and nonequilibrium pattern formation[1,2,24-31].

Active particles are often subject to thermal (Brownian) and athermal (active) sources of noise, arising typically from ambient fluctuations and internal mechanisms, respectively[32-36]. Among various models of active motion[24-26], active Brownian particles (ABPs)—featuring a constant self-propulsion speed subject to rotational Brownian noise—have emerged as a minimal model, capable of capturing key aspects of active particulate systems, with the greater part of theoretical works on this model being devoted to computer simulations in two spatial dimensions[24-27]. These include motility-induced phase separation (see, e.g., Refs.[37-39] and references therein), nonequilibrium wall accumulation, boundary layering and barrier trapping (see, e.g., Refs.[40-52] and references therein), and collective self-organized motion (see, e.g., Refs.[26,29,30,53,54] and references therein). Although these phenomena are enabled and sustained by particle activity, they can arise only in the presence of interparticle and particle-wall interactions, including Vicsek[26,29,55] and Vicsek-like[56-59] alignment interactions and hydrodynamic interactions[57,60]. Most of the aforementioned phenomena can emerge in active systems with only short-ranged steric interactions between the constituent particles, signifying their generic nature.

Models of ABPs with short-ranged steric interactions are particularly useful in the study of active colloidal suspensions in a base molecular liquid. When ABP dynamics is supplemented by the translational Brownian noise, the properties of the (bulk or confined) suspension can be tuned continuously from its reference equilibrium state, when particle self-propulsion is set to zero, to an arbitrary nonequilibrium steady state, merely by adjusting a single control parameter, which is known as Péclet number and is defined as the ratio of the self-propulsion and rotational diffusivity timescales. The ABP framework thus facilitates direct comparisons between analogously designed active and passive systems.

In equilibrium colloid science, bath-mediated effective interactions between colloidal inclusions have been of primary interest, given that such interactions directly determine the phase behavior of the suspended colloids. The inclusions can be of a whole range of shapes and dimensions (from nano to microscale), while the bath can, for instance, be a polar molecular solvent (as is often the case in the study of Casimir-van-der-Waals interactions) or a colloidal or polymeric suspension of smaller-size particles itself (as is often the case in the study of depletion forces)[61]. In the latter case, the effective bath-mediated interactions stem from the steric depletion of smaller bath particles (or depletants) from the vicinity of the inclusions, yielding a short-ranged and attractive, depletion, interaction force between the inclusions[62,63].

When colloidal inclusions are suspended in an active bath of ABPs[31,47-52,64-78], the nature of bath-mediated interactions is altered significantly, as both attractive and repulsive forces of significantly larger magnitude and range are possible, depending on the size and concentration of the bath particles and, more importantly, also their activity strength. Active bath-mediated interactions between suspended inclusions can substantially be affected also by the presence of any confining boundaries[51,77], shape, mobility and concentrations[48,51,52,68,75,78] and, if present, dissimilarity[49,69] and permeability[52,79] of the inclusions. On the other hand, the shape of active particles itself plays a crucial role in determining their bulk properties[37,80-83], near-wall accumulations at rigid[40-43,56,84] and deformable[85-87] boundaries and, hence, the active (or ‘swim’) pressure they exert on these boundaries[56,64,79,88,89], and the effective interactions they mediate between external objects[47,66,68,69]. For instance, the nonequilibrium wall accumulation of self-propelled rods turns out to be weaker than that of self-propelled disks[41,42]. Active rods are often assumed to self-propel along their longitudinal (or long) axis and the torques experienced by the rods due to collisions with the wall reduces their near-wall ‘detention’ time. The effects of active-particle shape on their clustering, segregation and collective motion have been explored for a range of different shapes such as active semiflexible filaments[90], self-propelled triangles[81], squares[82], crescents and other elongated different shapes with fore-aft asymmetry[83]. Self-propelled rods with longitudinal and transverse self-propulsion direction (being parallel and perpendicular to the longitudinal axis of the rods, respectively) have been realized experimentally using, e.g., bimetallic Janus nano/microrods[13-17]. In particular, phase behaviors of transversally self-propelling rods have numerically and experimentally been investigated[17]. Despite these recent works, systems of active rods with transverse or, more generally, nonaxial (tilted) self-propulsion direction remain largely unexplored.

Here, we focus on a minimal model of active Brownian rods with constant self-propulsion speed and, using Brownian Dynamics simulations, explore the effects of nonaxial self-propulsion of the rods on their spatial distributions and the effective interactions they mediate between two colloidal inclusions. The paper is organized as follows. We introduce our model and methods in “Model and methods” and discuss our simulation results for the distribution of the active Brownian rods in “Results for spatial distribution of active rods”, followed by an analysis of the bath-mediated interaction force on the inclusions in “Results for effective interaction force”. The paper is summarized in “Summary”.

Model and methods

Our two-dimensional model comprises two identical, impermeable, disklike colloidal inclusions of diameter , placed at fixed positions at intersurface separation d within a bath of N identical, stiff, active Brownian rods of width and length l, moving in a still background medium of viscosity ; see Fig. 1. The x-axis is taken along the center-to-center axis of the inclusions with the origin placed in the middle, yielding the center position vectors for the ‘left’ () and the ‘right’ () inclusion in Fig. 1a, with being the x unit vector. We consider two separate cases of interacting and noninteracting rods. In the former case, active rods repel one another via nearly hard and short-ranged, steric pair potentials, as will be specified later, causing ‘contact’ repulsive forces and torques between the rods. In the latter case, the active rods act as ‘phantom’ particles with no interactions among themselves. In either case, they always interact sterically with the disklike inclusions.

Figure 1

(a) Schematic view of two fixed inclusions of diameter placed at intersurface distance d in a bath of identical active Brownian rods. (b) Active rods are of width and length l and self-propel due to active forces of fixed magnitude and self-propulsion angle relative to the orientation unit vector along the longitudinal axis of the rods.

The rods self-propel due to active forces of fixed magnitude and fixed self-propulsion angle, , with labeling the individual rods. The angle is defined relative to the instantaneous orientation unit vector ; i.e., . We define , where is the parametrizing polar angle measured from the x-axis; see Fig. 1b. The configuration space of the active rods is thus spanned by the center-of-mass position vectors and the orientation vectors that obey overdamped Langevin dynamics[91,92]where is the rotation operator[91,92] and is the unconstrained partial differentiation operator with respect to the Cartesian components of . and are the relevant translational () and rotational () mobility tensors,where is the unit tensor, and , and are the (Stokes) friction coefficients associated with translation parallel () and perpendicular () to the longitudinal axis of the rods and rotation () in the x–y plane, respectively. Also, in Eqs. (1) and (2), and are, respectively, the translational and rotational, Gaussian, white noises of zero mean, , and two-point correlations , and , and is the out-of-plane (normal) unit vector. The friction coefficients depend on the length and aspect ratio of the rods and the medium viscosity according to[93]where the so-called correction coefficients , and depend only on the aspect ratio and are given (for ) by[93]

Finally, in Eqs. (1) and (2), U gives the sum of the interaction potentials between all ‘particles’ in the system. We discretize each of the rods to a linear array of equidistant beads of diameter , with being the center-to-center distance between two adjacent beads in a given rod. In a given configuration of inclusions (labelled by ) and active rods (labelled by , with their constituent beads labelled further by ), we havewhere is the position of the ath bead of the ith rod, etc, and is the Weeks-Chandler-Andersen (WCA) pair potential[94]

For the ensuing steric interaction between two individual beads (from two different rods), is taken as the center-to-center distance of the beads and . For the interaction between a given bead and either of the inclusions, gives the center-to-center bead-inclusion distance and . In either case, we use the same interaction strength .

Simulation methods and parameters

We use standard Brownian dynamics methods[95] to numerically solve the governing Eqs. (1) and (2). The relevant parameter space of the system can be spanned by the following set of dimensionless parameters: The rescaled inclusion diameter , the rescaled intersurface distance of the inclusions, , the area fraction of active rods , the self-propulsion angle , and the Péclet number (relative strength of the self-propelling force), .

Since our focus will be on the role of nonaxial self-propulsion, we fix the aspect ratio as ( beads), the inclusion diameter as and the area fraction of active rods as (see below). The Péclet number is varied over the interval , that covers a realistic range of values; e.g., synthetic active rods of length 2  with self-propulsion speeds in the range 4–18  have been used in experiments[13], yielding Péclet numbers that, using our definitions, fall in the range –13. We restrict the self-propulsion angle to the first angular quadrant , as the results for other -quadrants can be reproduced (within the simulation error bars) using appropriate symmetry arguments; e.g., the results for a given self-propulsion angle, , in the first quadrant can be mapped to those of the second quadrant angle – using a left-right reflection (relative to the y-axis). For the sake of comparison with a monodisperse system of rods, we shall later consider two additional cases: a random mixture of active rods with self-propulsion angles being extracted from a uniform probability distribution over the interval , and a binary mixture of longitudinally () and transversally () self-propelling rods with relative fractions and , respectively.

In the simulations, we use a rectangular box of periodic boundaries of side lengths and active rods, giving the desirable area fraction . We fix the WCA interaction strength as , ensuring sufficiently impermeable particles (note that strongly self-propelling rods can nevertheless partially penetrate into the inclusions, an effect that we shall return to in “Noninteracting active Brownian rods”). In the simulations, we use the timesteps of rescaled size and run the simulations between and timesteps, with and initial steps used for relaxation to the steady state in the cases of interacting and noninteracting rods, respectively. The averaged quantities are then calculated by further averaging them over 36 statistically independent simulations.

Results for spatial distribution of active rods

Isolated inclusions

We begin our discussion by examining the distributions of active rods around the disklike inclusions in the steady state. Figure 2a shows the number density of active rods, , as a function of the radial distance, , from the surface of an isolated inclusion (practically, either of the two inclusions can be treated as an isolated one, provided that the inclusions are placed at sufficiently large intersurface distances to make them decoupled; here, we take ). Being a well-known trait in systems of self-propelled particles[41-43], the active rods in the present context also exhibit pronounced accumulations in the close proximity of the inclusions. This is indicated by the near-surface peaks (occurring at spacings nearly equal to the rod width, ) in the density profiles of Fig. 2a. The peaks thus portray not only the excessive nonequilibrium accumulation of the rods near the inclusions, but also their steric layering, akin to that of disklike active particles, as extensively explored in Refs.[50-52]. In the absence of hydrodynamic effects and/or particle-specific features (such as swimmer tails), surface accumulation of active particles arises generically due to their prolonged detention near hard boundaries: upon surface contact, the normal-to-surface component of the self-propulsion velocity pushes and holds the active particles against the surface, while its parallel-to-surface component persistently drives the particles along the boundary contour, until the rotational diffusion sets in, enabling particle escape from the boundary regions. In the case of active rods, this mechanism brings into play the nonaxial self-propulsion angle, , as already indicated by the data in Fig. 2a.

Figure 2

(a) Rescaled number density of active rods as a function of the rescaled radial distance, , from the surface of an isolated inclusion for and . (b–d) show color-coded maps of the active-rod density field for and (left to right). Panels (e–g) show the mean vector field, , around the right inclusion for , and , respectively; here, the arrows show the average orientation of and the color-coding indicates . In all cases, and .

The -dependence is seen more clearly in the color-coded two-dimensional maps of the density field, , in the simulation box; see Fig. 2b (), c () and d (). The density peaks of Fig. 2a appear here as distinct high-density (dark yellow to red) layers, or ‘rings’, around the inclusions. As increases from 0 to at fixed , and the transverse component of the self-propulsion velocity, , increases, the active rods show a significantly stronger surface accumulation (denser rings) over a wider range of radial distances from the inclusions (more rings). This picture can be corroborated by the vector field, , in Fig. 2e (), f () and g (), where we show the results in the regions around the right inclusion (with the results being indistinguishable from those obtained for the left inclusion; not shown). The arrows show the mean orientation of (see Fig. 1b) and the color-coding indicates (note that, while , the local average ). In regions away from the inclusions, the active rods adopt random orientations on average (). When the self-propulsion is longitudinal (, Fig. 2e), the rods point radially toward their respective (here, the right) inclusion on average. When the self-propulsion is nonaxial ( and ; Fig. 2f,g, respectively), the mean rod orientations form a counterclockwise vortex within a finite radial distance from the inclusions. This effect is stronger for purely transverse self-propulsion, as the motion of rods within the vortex is strongly hindered. For tilted self-propulsion (), the nonvanshing longitudinal component of self-propulsion causes chiral motion of rods around the inclusions in counterclockwise direction. As we shall see later, this chiral motion acts as the source of noncentrality in the effective force between the inclusions, when they are placed at small separations.

Juxtaposed inclusions

When the two inclusions are brought into small intersurface separations, see Fig. 3a–c with , the active rods are exposed to both types of convex and concave boundary curvatures. These types of curvatures occur, respectively, on the far sides of the inclusions and in the intervening region (the upper and lower wedge-shaped gaps) formed between the inclusions, respectively. In general, surface accumulation of active particles is stronger on concave boundaries as compared with convex ones[43,64]. This effect enhances particle densities in the intervening gaps as opposed to those on the far sides of the inclusions. The relative accumulation in the gaps is strongest (weakest) in the case of longitudinally (transversally) self-propelling rods (panels a and c), making the spatial distribution of the rods around the inclusion dimer strongly (weakly) anisotropic. We shall discuss the ramifications of this strong versus weak anisotropy on the effective bath-mediated interactions in following sections.

Figure 3

(a–c) show color-coded maps of the spatial distribution of active rods for and , respectively with fixed and .

While the spatial distributions of active rods for and (panels a and c) exhibit up-down symmetry, the situation turns out to be different for other values of (panel b for ). For the parameters in Fig. 3, we find an excess (deficit) of active rods above (below) the right inclusion, with the situation being reversed in the case of the left inclusion. The asymmetric density imbalance can be seen more clearly in Fig. 4a, where we plot the conventionally defined cumulative density, , of active rods in a vertical stripe , embedding the right inclusion, where and (recall the inclusion diameter of , the intersurface distance and length of rod in Fig. 3, giving ). The alternating peaks are again due to the active-rod ring formation around the inclusion, and the excess (deficit) of the rods above (below) the inclusion are mirrored by the larger (smaller) peaks in the plot.

Figure 4

(a) Rescaled cumulative density of active rods defined along the y-axis and relative to the right inclusion (see the text for definitions) for , and two juxtaposed inclusions at intersurface distance . The dashed lines indicate the y-coordinates at the upper and lowermost points of the inclusion. (b, c) show schematic views of the mechanism, causing asymmetric rod accumulation on the up versus down sides of the inclusions for and .

We emphasize that the asymmetric up-down population imbalance emerges only when the inclusions are juxtaposed and when the self-propulsion axis is tilted . The underlying mechanism for this effect can be understood by noting that, when the active rods come into contact with the inclusions, they adopt circular trajectories, going in counterclockwise direction around the inclusions. This type of chiral surface motion is driven by the longitudinal component of the self-propulsion, while the transverse component keeps the rods in contact with the inclusions. The effect is shown schematically in Fig. 4b. For comparison, Fig. 4c depicts the clockwise motion that would arise for . For the case considered here (), the counterclockwise surface motion proceeds around the right (left) inclusion and is hindered, when the active rods reach the upper (lower) wedge-shaped gaps. Hence, the prolonged detention time and imbalanced crowding of the active rods above (below) the right (left) inclusion.

Results for effective interaction force

Since the distribution of active particles around isolated inclusions is rotationally symmetric (Fig. 2), the instantaneous forces imparted on the inclusions by the steric interactions (collisions) with the active particles average out to zero in the steady state. Any net force on the inclusions that may arise at small intersurface separations can thus be interpreted as the effective interaction force mediated by the active bath between the two inclusions. In the present context, such interactions originate from the steric overlaps between the active-particle rings formed around individual inclusions, as they are bought to small intersurface separations[50-52]. To proceed, we focus on the inclusion shown on the right side () in Fig. 1 and denote the effective force experienced by this inclusion by . This quantity is computed in our simulations using the relationwhere the brackets denote the steady-state average and the expression inside is the instantaneous forces imparted by the constituent beads of the rods on the said inclusion. depends on various system parameters, with the dependence on the intersurface distance, d, explicitly indicated. For brevity, we may refer to as the effective (interaction) force between or acting on the inclusions without further specifying it as the force on the right inclusion. One should bear in mind that the force on the left inclusion is of equal magnitude and opposite direction, as explicitly verified in the simulations (with the error bars). Thus, by definition, () identifies a ‘repulsive’ (‘attractive’) longitudinal force component, and () identifies an ‘upward’ (‘downward’) transverse force component. When , the effective interaction between the inclusions will be noncentral. In this case, there is an anti-parallel transverse force component of the same magnitude acting on the left inclusion and, hence, a counterclockwise (clockwise) torque on the inclusion dimer.

Special cases with and

Figure 5a shows the simulated effective interaction force, Eq. (10), as a function of the intersurface distance, d, for different values of the Péclet number, , and in the special case of longitudinally (axially) self-propelling rods with . The effective force in this case is found to be central; i.e., it only has a longitudinal component, . Its root cause is the broken left-right symmetry in the active-particle distribution around each individual inclusion. The preserved up-down symmetry in this case (Fig. 3a) ensures , as will explicitly be verified later (see Fig. 6b).

Figure 5

Rescaled effective force acting on the inclusions as a function of their rescaled intersurface distance, , for (a) longitudinally and (b) transversally self-propelling rods and different values of , as indicated by the legends in (a).

Figure 6

(a) Longitudinal () and (b) transverse () components of the rescaled effective force acting on the inclusions as a function of the rescaled intersurface distance, , for and different values of , as indicated by the legends in (a).

The effective force F(d) in Fig. 5a is repulsive because of the previously discussed preferential accumulation of active rods in the wedge-shaped gaps between the inclusions (“Juxtaposed inclusions”). It exhibits a characteristic alternating profile due to the successive intersections that occur between the active-particle rings of the two inclusions, as they are brought into surface contact (hence, the interpeak spacings in the force-distance profiles also exhibit a rather intricate pattern, even though they remains of the order of the rod width, σ). These features of the force profiles in Fig. 5 are similar to those reported for active Brownian disks in Refs.[50-52]. The self-propelled rods, however, experience inter-particle and particle-inclusion torques that exacerbate their trapping and accumulation in the wedge-shaped gaps between juxtaposed inclusions; an effect that is not present in the case of disklike active particles. Increasing enhances the magnitude of the effective force because of stronger rod-inclusion collisions, but suppresses its oscillatory behavior. The latter is due to increased run lengths of active rods that, upon exceeding the inclusion size, lead to broadening and dilution of the circular rings they form around the inclusions.

The force-distance profiles in the case of transversally self-propelling rods (), Fig. 5b, show similar features to those of longitudinally self-propelling rods. However, the effective force here is of relatively smaller magnitude and of larger range of action (compare Fig. 5a,b). The extended range of interaction can be understood by noting that the active-particle layering around the inclusions for appears over a respectively larger range of radial distances from each inclusion (Fig. 3, compare panels a and c). But since the distribution of transversally self-propelling rods around the inclusions is relatively more isotropic, as compared with the case of longitudinally self-propelling rods, the repulsive force due to the excess accumulation of active rods in the gaps between the inclusions is also relatively weaker.

Tilted (nonaxial) self-propulsion

When the self-propulsion axis is tilted (), the imbalance of active-rod population on the upper (lower) side of the right (left) inclusion, see “Juxtaposed inclusions”, yields a nonvanishing transverse force component on the inclusions. In Fig. 6a,b, we show the longitudinal () and transverse () components of the effective force on the inclusions are shown as functions of the intersurface distance, d, for and selected values of the self-propulsion angle . We also include the results for a random sample (cyan symbols) with in the active mixture being distributed uniformly within the interval .

As seen in Fig. 6a, at small separations , the peak amplitudes of the longitudinal interaction force, , decrease as increases from 0 (pink squares) to (dark blue triangle-downs); in other words, decreases as the self-propulsion direction becomes increasingly more tilted away from the rod axis. This is because, upon increasing , active-rod accumulation around the inclusions becomes increasingly more isotropic (i.e., more active rods accumulate on the far sides of the inclusions), reducing the repulsive longitudinal force. The situation is reversed at large separations and becomes consistently larger as increases from 0 to ; this is because falls off more weakly with the separation at larger . The increase in the interaction range, as previously noted for in Fig. 5b, thus occurs continually as increases. For the most part, the data for the uniformly distributed (random) case fall closer to (green triangle-ups). Figure 6b confirms that indeed vanishes within the simulation error bars in the special cases of and over the whole range of intersurface separations. In other cases of , including the case of randomly distributed , is negative (i.e., according to our conventions, it pushes the left/right inclusion up/down, causing a clockwise effective torque on the inclusion dimer) and shows a distinctly long range.

Further insights into the -dependent behavior of the interaction force can be obtained by fixing the intersurface distance as, e.g., (surface contact), and examining the detailed behavior of the force components over the axis. The aforementioned trends in are clearly mirrored by the longitudinal contact force, , shown in Fig. 7a. As seen, exhibits consistently larger (smaller) repulsive strengths at larger (larger ). While, at small , we find a nearly monotonic decrease in down to zero, the expected nonmonotonic alterations (due to stronger accumulation of active rods around the inclusions) emerge at large , with a pronounced hump in the subinterval . The transverse contact force in Fig. 7b expectedly vanishes at the end points of the depicted -interval but the figure here reveals a pronounced global minimum (corresponding to a maximal force magnitude) in the subinterval at around .

Figure 7

(a) Longitudinal () and (b) transverse () components of the rescaled effective force acting on the inclusions at surface contact as a function of for different values of , as indicated by the legends in (a). (c) The effective work, W, as defined in the text, plotted as a function of for . Symbols show simulation data for a monodisperse system of active rods with given , and the dashed line shows W for a random dispersion with uniformly distributed . Inset shows W for a binary mixture of longitudinally () and transversally () self-propelling rods as a function of the fraction of transversally self-propelling rods, .

Another useful quantity that can complement the foregoing ones is the effective work required to bring the inclusions, in a quasistationary fashion from a sufficiently large intersurface distance, , where the interaction force vanishes, into surface contact. We compute this quantity along the x-axis as , by taking that gives vanishing interaction force within the simulation error bars. In an equilibrium setting with , W corresponds to the so-called potential of mean force between the inclusions[62]. W is shown in Fig. 7c, main set, as a function of the nonaxial self-propulsion angle, . The symbols show the results for a monodisperse system of active rods with given , while the dashed horizontal gives the corresponding value for a random dispersion of active rods, with uniformly distributed . As seen, for the parameters in the figure, W in the case with (longitudinal self-propulsion) nearly matches that in the case of randomly distributed . W then drops from this value and decreases down to a global minimum at and then turns up and matches the random sample (dashed line) at around , growing to a global maximum at , before falling off as tends to . These extrema signify the situations, where the overall repulsion between the inclusions is minimized/maximized.

As defined, W integrates the complex pattern of interaction force profiles (see Fig. 6) into a single quantity. Even though the outcome (Fig. 7c) gives a seemingly more tangible representation of the -dependent behavior of the effective interactions, the detailed features of W as a function (such as its limiting and extremal values) remain to be understood. In the inset of Fig. 7c, we show W in the case of a binary mixture of longitudinally () and transversally () self-propelling rods with relative fractions and , respectively. W is a convex function of , maximized in the monodisperse cases () and a globally minimized at . This indicates that the overall repulsion between the inclusions can also be suppressed by producing a binary mixture of active rods.

Noninteracting active Brownian rods

As a useful baseline to further elucidate the role of interparticle interactions in the present context, we consider an analogous system of noninteracting (‘phantom’) active Brownian rods by switching off (only) the steric inter-rod repulsions. Figure 8a,b show the corresponding interaction forces acting on the inclusions for longitudinally () and transversally () self-propelling phantom rods, respectively. In these cases, the effective force is again central and .

Figure 8

Rescaled effective force acting on the inclusions due to noninteracting (phantom) active rods as a function of the rescaled intersurface distance, , for (a) longitudinally and (b) transversally self-propelling rods and different values of , as indicated by the legends in (a). (c) Analytically predicted d-dependence of the active depletion-like force, Eq. (11), is fitted (lines) to the simulation data for , and 40 (blue and red symbols, respectively) and .

Longitudinal self-propulsion

For , the effective force is repulsive, as in the interacting case, but it shows qualitatively different features: it is of significantly shorter range, of significantly larger magnitude, and with only a single sharp peak, when compared with the case of interacting rods (compare Fig. 8a and Fig. 5a). These differences reflect the differences that occur in the spatial distribution of active rods, when the steric inter-rod repulsions are turned off. For the sake of brevity (and given that the noninteracting model is an idealistic example and not at the focus of the present work), here we merely summarize some of the key aspects of the simulated phantom rod distributions and present plots of spatial distributions of active rods in the “Supplementary Information”.

The most notable difference between spatial distributions of phantom active rods and their interacting counterparts in “Results for spatial distribution of active rods” is that phantom rods do not form the steric rings around the inclusions. Thus, in the case of longitudinally () self-propelling phantom rods, as the inclusions are brought into sufficiently small separations, the preferential accumulation of the rods in the concave gaps formed between the inclusions changes into their excessive overcrowding in those intervening regions (see Supplementary Fig. S1a online). This leads to the much stronger effective repulsion between the inclusions with force strengths that can be larger by more than an order of magnitude (compare Fig. 8a and Fig. 5a), albeit in a significantly shorter range of separations, . The force-distance profiles in Fig. 8a increase up to a maximum then suddenly drop to zero. The peak location indicates the typical gap width d between the inclusions for which the active rods can start squeezing through the gap. This is enabled by partial penetration of the rods into the interfacial regions of the inclusions, as the rods self-propel from and to the top/bottom fluid regions connected by the intervening gap between the inclusions (note that the phantom rods can still interact sterically with the inclusions). These processes release the longitudinal stress and reduce the repulsive force, which appears as a sharp drop beyond the peak force. Hence, as the Péclet number increases, the peak location in Fig. 8a shifts to the left, as the stronger self-propulsion force helps the active rods to squeeze through even narrower intervening gaps (of widths down, e.g., to 45% or 30% of the rod width; see and 40 in the figure). As is increased (), the relative changes in the peak amplitude remains small, as total trapping of active rods in the intervening gap between the inclusions maintains the maximum longitudinal rescaled force () acting on the inclusions at a fixed level.

Transverse self-propulsion : depletion-like attraction

For , the absence of steric inter-rod repulsions and, primarily, the ring formation, of phantom rods around the inclusions, also leads to significant modifications in the resulting interaction force as compared with the interacting case (compare Fig. 8b and Fig. 5b). The features turn out to be quite different from those in the case of as well (Fig. 8a). The short-ranged effective force is found to be attractive in distances less than or comparable to the rod width, , followed by a weak repulsive hump, before leveling off to zero (Fig. 8b).

The effective force in Fig. 8b resembles the typical force profiles in the context of equilibrium depletion interactions between colloidal particles[62,63]. This finding might appear unexpected, but it reflects the tendency of transversally self-propelling phantom rods to adhere tangentially (from their side edges as pushed by the self-propulsion force) to the inclusions and to act similar to the so-called depletants in the said context. The key difference with the case of longitudinally self-propelling phantom rods is that—instead of overcrowding the intervening gap regions—transversally self-propelling phantom rods distribute rather isotropically around the circular boundaries of the inclusions (see Supplementary Fig. S1a,c online), especially as the Péclet number, , is sufficiently large (see below). Such a tendency was noted in the case of interacting rods as well (“Juxtaposed inclusions”), but the effect is exacerbated here by the absence of steric inter-rod repulsions between phantom rods. This leads to the formation of a thin interfacial layer of phantom rods around the inclusions with a thickness comparable to the rod width. When the inclusions are brought into small separations , the noted interfacial layer is sterically depleted from the intervening gap between the inclusions (see Supplementary Fig. S1c online), engendering a steric depletion mechanism similar to the traditional one in the equilibrium context[62,63] and, hence, a short-range attraction.

It is important to note that the resulting depletion-like attraction remains of predominantly active nature because of the finite Péclet number, and is thus of much greater magnitude than its equilibrium value (). The attractive force behavior as a function of d in Fig. 8b can be estimated analytically, when d is sufficiently small (). We assume that is sufficiently large () to warrant tangential and isotropic adherence of active rods around the inclusions (except in the intervening gap). The longitudinal force component due to a single rod of orientation angle , adhered to the right inclusion and pushing against it, can be estimated as , where is limited by the gap exclusion as , with . The longitudinal effective force is then estimated asup to a prefactor that cannot be determined within the above analysis. The analytically predicted d-dependence of the active depletion-like force in Eq. (11) fits the simulation data for and 40 very well (compare line and symbols in Fig. 8c).

The forgoing discussions elucidate the main qualitative and quantitative differences (in sign, magnitude and range) that we find between the effective forces mediated by longitudinally () and transversally () self-propelling phantom rods between the inclusions (Fig. 8a,b). Our numerical inspections however indicate that the repulsive humps in the force profiles in Fig. 8b arise as a result of yet another effect; namely, the transversally self-propelling phantom rods can get horizontally trapped in the wedge-shaped gaps between the inclusions, when the gap size is in the range (between the width of the rod and its length). This effect will be present only at weak to intermediate values of , as otherwise the strong tendency of the rods to tangentially adhere to the inclusions suppresses such rod-bridging configurations.

Tilted self-propulsion

Finally, for nonaxial self-propulsion angles , the effective interaction between the inclusions becomes noncentral as in the case of interacting rods (“Juxtaposed inclusions”). The x and y-components of the resulting force, and , as functions of the intersurface distance d are shown in Fig. 9a,b for and a few different values of . The longitudinal force profiles in Fig. 9a display similar features to those with , even though the overall force magnitude decreases as increases. On close approach to purely transverse self-propulsion with , the force magnitudes significantly reduce and the force profiles change to those of Fig. 8b (for illustration purposes, we have excluded , whose data become indiscernible in the scale of the graph). The transverse force component , Fig. 9b, again vanishes for the special cases and ; however, for other self-propulsions angles, the transverse force profiles due to phantom active rods also differ significantly from those found in the case interacting rods (Fig. 6b), even though is still negative and creates an effective clockwise torque on the inclusion dimer. The qualitative differences seen in the force profiles in the case of phantom rods as opposed to sterically interacting rods are reflected by the work, W, done on the inclusions as they are brought into contact along x-axis from far enough distances; see Fig. 9c (here, we take ). Unlike the case of interacting rods in Fig. 7c, W is convex here and varies mostly monotonically as a function of , with a sharp drop as approaches . The work drops negative values in this limit, signifying the attractive depletion-like force profiles of Fig. 8b.

Figure 9

(a) Longitudinal () and (b) transverse () components of the rescaled effective force acting on the inclusions due to phantom active rods as a function of the rescaled intersurface distance, , for different values of , as indicated by the legends in (a). (c) Same as Fig. 7c but here we plot the effective work, W, for phantom active rods. In all cases, . The dashed line here indicates the reference value .

Summary

Using extensive Brownian Dynamics simulations of a minimal model of active Brownian rods with constant self-propulsion speed in two dimensions, we have investigated the role of nonaxial self-propulsion, identified by a tilt angle , on (1) the steady-state distributions of the active rods around two fixed, impermeable, disklike inclusions immersed in the bath, and (2) the effective interaction forces mediated by the active bath between the inclusions. Because the active rods are assumed to be rigid and interact with steric repulsions, they form layered structures (rings) around the inclusions at sufficiently high Péclet number, . This type of nonequilibrium particle layering occurs over a significantly larger range of radial distances from the inclusion boundaries in the case of transversally self-propelling rods () as compared with the case of longitudinally self-propelling rods (). This in turn leads to effective interaction forces of respectively longer range of action between the inclusions in the former case. In both these cases, the interaction force is central, with only a longitudinal component along the center-to-center axis of the two inclusions. It is also repulsive because of a preferential accumulation of rods in the intervening wedge-shaped gaps formed between the inclusions at sufficiently small separations, even though the effect is weaker in the case of transversally self-propelling rods, as they tend to adhere tangentially to the inclusions, displaying a relatively more isotropic spatial distribution around them. When the self-propulsion axis is tilted relative to the rod axis, the effective force acquires a finite transverse component and becomes noncentral. This effect is due an asymmetric imbalance of active-rod accumulation on the up and down sides of the inclusion dimer, and translates into a finite clockwise (counterclockwise) torque on the dimer for (). For comparison purposes, we also presented simulations result for a random mixture of active rods with uniformly distributed self-propulsion angles, a binary mixture of longitudinally and transversally self-propelling rods, and the ideal case of phantom active rods, lacking steric inter-rod repulsions.

Noncentral interactions have been reported in simulations of two hard inclusions embedded in a (dry) two-dimensional bed of granular disks, activated by a shaking substrate with in-plane oscillations[72]. In this latter case, nonaxial self-propulsion is induced along a fixed external axis determined by the in-plane direction of substrate oscillations, contrasting the inherent self-propulsion tilt assigned to individual particles relative to their body axis in the present model.

Our study is based on a customarily used minimal model of active particles in two dimensions (see, e.g., Refs.[37,40,44,45,47,48,50-53,83,84] and reference therein), which neglects hydrodynamic interactions and other types of nonsteric interactions such as alignment (e.g., Vicsek) interactions[26,29,55]. It will thus be useful to study possible roles of such interparticle interactions in future works, especially since (quasi-)two-dimensional realizations of the current model would be feasible in actual settings involving swarming bacteria[96], shaken (mono-layered) granular matter[46,54,97,98], microswimmers confined to fluid-fluid interfaces[99], and microswimmers confined by acoustic tweezers[100]. A careful analysis of hydrodynamic effects due to a fluid background will be of particular interest in (quasi-)two-dimensional settings. From this perspective, minimal models such as the one adopted here (see Ref.[101] for a review) should be viewed as first-step models that revolve around the role of particle activity. The present model can be extended to incorporate other tractable features of real systems, such as the confinement effects[51,77], deformability[90] and complex shape of self-propelled particles[81-83], and also mobility[47,66,72,75,78], shape dissimilarity[49,69] and permeability[52,79] of the inclusions. As we have focused on two-body interactions between the inclusions, it will be of interest to study three-body[50] and many-body effects in clusters[78] and linear arrays of inclusions. In the latter case, the noncentrality in the effective interactions is expected to be suppressed by increasing the array size and to vanish in the limit of infinite arrays (where each inclusion has similar neighbors on both sides).

Supplementary Information.