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.


Model and methods
Our two-dimensional model comprises two identical, impermeable, disklike colloidal inclusions of diameter σ c , 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 R 1 = −R 2 = −(d + σ c )x/2 for the 'left' ( R 1 ) and the 'right' ( R 2 ) inclusion in Fig. 1a, with x 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.
The rods self-propel due to active forces F i,SP (t) of fixed magnitude F SP = |F i,SP (t)| and fixed self-propulsion angle, θ , with i = 1, . . . , N labeling the individual rods. The angle θ is defined relative to the instantaneous orientation unit vector û i (t) ; i.e., û i (t) · F i,SP (t) = F SP cos θ . We define û i (t) = (cos φ i (t), sin φ i (t)) , where φ i (t) is the parametrizing polar angle measured from the x-axis; see Fig. 1b. The configuration space of the active rods www.nature.com/scientificreports/ is thus spanned by the center-of-mass position vectors {r i (t)} = { x i (t), y i (t) } and the orientation vectors {û i (t)} that obey overdamped Langevin dynamics 91,92 where Rû i =û i × ∇û i is the rotation operator 91,92 and ∇û i is the unconstrained partial differentiation operator with respect to the Cartesian components of û i . M T and M R are the relevant translational ( T ) and rotational ( R ) mobility tensors, where I is the unit tensor, and ζ , ζ ⊥ and ζ R are the (Stokes) friction coefficients associated with translation parallel ( ) and perpendicular ( ⊥ ) to the longitudinal axis of the rods and rotation ( R ) in the x-y plane, respectively. Also, in Eqs. (1) and (2), ξ i,T (t) = ξ i,� (t)û i (t) + ξ i,⊥ (t)û i (t) ×ẑ and ξ i,R (t) = ξ i,R (t)ẑ are, respectively, the translational and rotational, Gaussian, white noises of zero mean, �ξ i,� (t)� = �ξ i,⊥ (t)� = �ξ i,R (t)� = 0 , and two-point correlations �ξ i,� (t)ξ j,� (s)� = ζ � k B Tδ ij δ(t − s) , �ξ i,⊥ (t)ξ j,⊥ (s)� = ζ ⊥ k B Tδ ij δ(t − s) and �ξ i,R (t)ξ j,R (s)� = ζ R k B Tδ ij δ(t − s) , and ẑ is the out-of-plane (normal) unit vector. The friction coefficients depend on the length and aspect ratio p = l/σ of the rods and the medium viscosity according to 93 where the so-called correction coefficients ν , ν ⊥ and ν R depend only on the aspect ratio and are given (for 2 ≤ p ≤ 30 ) 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 n b = l/l b − 1 equidistant beads of diameter σ , with l b = σ/2 being the center-to-center distance between two adjacent beads in a given rod. In a given configuration of inclusions (labelled by α = 1, 2 ) and active rods (labelled by i, j = 1, . . . , N , with their constituent beads labelled further by a, b = 1, . . . , n b ), we have where r a i is the position of the ath bead of the ith rod, etc, and U WCA is the Weeks-Chandler-Andersen (WCA) pair potential 94 For the ensuing steric interaction between two individual beads (from two different rods), |r| is taken as the center-to-center distance of the beads and σ eff = σ . For the interaction between a given bead and either of the inclusions, |r| gives the center-to-center bead-inclusion distance and σ eff = (σ + σ c )/2 . 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 σ c /σ , the rescaled intersurface distance of the inclusions, d/σ , the area fraction of active rods , the self-propulsion angle θ , and the Péclet number (relative strength of the self-propelling force), Pe = σ F SP /(k B T).
Since our focus will be on the role of nonaxial self-propulsion, we fix the aspect ratio as p = 3 ( n b = 5 beads), the inclusion diameter as σ c = 12σ and the area fraction of active rods as = 0.2 (see below). The Péclet number is varied over the interval Pe ∈ [5,40] , that covers a realistic range of values; e.g., synthetic active rods of length 2 µm with self-propulsion speeds in the range 4-18 µm/s have been used in experiments 13 , yielding Péclet numbers that, using our definitions, fall in the range Pe ≃ 3-13. We restrict the self-propulsion angle to the first angular quadrant θ ∈ [0, π/2] , 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 (1) Scientific Reports | (2021) 11:23100 | https://doi.org/10.1038/s41598-021-02295-y www.nature.com/scientificreports/ 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 θ ∈ [0, π/2] , and a binary mixture of longitudinally ( θ = 0 ) and transversally ( θ = π/2 ) self-propelling rods with relative fractions 1 − n π/2 and 0 ≤ n π/2 ≤ 1 , respectively.
In the simulations, we use a rectangular box of periodic boundaries of side lengths L x = L y = 53.5σ and N = 200 active rods, giving the desirable area fraction � = We fix the WCA interaction strength as ǫ = 20k B T , 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 (3k B T/πησ 3 )δt = 10 −4 and run the simulations between 1.3 × 10 7 and 5 × 10 7 timesteps, with 5 × 10 6 and 10 7 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, ρ(r) , as a function of the radial distance, r = |r − R α | − σ c /2 , from the surface of an isolated inclusion (practically, either of the two inclusions α = 1, 2 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 d/σ = 14.5 ). Being a well-known trait in systems of self-propelled particles [41][42][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][51][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.
The θ-dependence is seen more clearly in the color-coded two-dimensional maps of the density field, ρ(x, y) , in the simulation box; see Fig. 2b ( θ = 0 ), c ( π/4 ) and d ( π/2 ). 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 π/2 at fixed Pe , and the transverse component of the self-propulsion velocity, Pe sin θ , 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 ( θ = 0 ), f ( π/4 ) and g ( π/2 ), 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 û = 1 , the local average û ≤ 1 ). In regions away from the inclusions, the active rods adopt random orientations on average ( û ≃ 0 ). When the self-propulsion is longitudinal ( θ = 0 , Fig. 2e), the rods point radially toward their respective (here, the right) inclusion on average. When the self-propulsion is nonaxial ( θ = π/4 and π/2 ; 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 (a) Rescaled number density of active rods as a function of the rescaled radial distance, r/σ , from the surface of an isolated inclusion for θ = 0, π/4 and π/2 . (b-d) show color-coded maps of the active-rod density field ρ(x, y) for θ = 0, π/4 and π/2(left to right). Panels (e-g) show the mean vector field, û , around the right inclusion for θ = 0 , π/4 and π/2 , respectively; here, the arrows show the average orientation of û and the colorcoding indicates û . In all cases, d/σ = 14.5 and Pe = 20.

Juxtaposed inclusions. When the two inclusions are brought into small intersurface separations, see
Fig. 3a-c with d/σ = 0.25 , 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) selfpropelling 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. While the spatial distributions of active rods for θ = 0 and π/2 (panels a and c) exhibit up-down symmetry, the situation turns out to be different for other values of θ (panel b for θ = π/4 ). 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, Q(y) = , embedding the right inclusion, where x 1 = 0 and x 2 = σ c + d/2 + l (recall the inclusion diameter of σ c /σ = 12 , the intersurface distance d/σ = 0.25 and length of rod l/σ = 3 in Fig. 3, giving x 2 /σ = 15.125 ). 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.
We emphasize that the asymmetric up-down population imbalance emerges only when the inclusions are juxtaposed and when the self-propulsion axis is tilted 0 < θ < π/2 . 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 π/2 < θ < π . For the case considered here ( 0 < θ < π/2 ), 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.   , c) show schematic views of the mechanism, causing asymmetric rod accumulation on the up versus down sides of the inclusions for 0 < θ < π/2 and π/2 < θ < π.

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][51][52] . To proceed, we focus on the inclusion shown on the right side ( α = 2 ) in Fig. 1  Special cases with θ = 0 and π/2. 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, Pe , and in the special case of longitudinally (axially) self-propelling rods with θ = 0 . The effective force in this case is found to be central; i.e., it only has a longitudinal component, F(d) = F(d)x . 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 F y = 0 , as will explicitly be verified later (see Fig. 6b). 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][51][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 Pe 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 ( θ = π/2 ), 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 θ = π/2 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 www.nature.com/scientificreports/ 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 ( 0 < θ < π/2 ), 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 ( F x ) and transverse ( F y ) components of the effective force on the inclusions are shown as functions of the intersurface distance, d, for Pe = 20 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 θ ∈ [0, π/2]. As seen in Fig. 6a, at small separations d/σ < 2 , the peak amplitudes of the longitudinal interaction force, F x , decrease as θ increases from 0 (pink squares) to π/2 (dark blue triangle-downs); in other words, F x 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 d/σ > 2 and F x becomes consistently larger as θ increases from 0 to π/2 ; this is because F x falls off more weakly with the separation at larger θ . The increase in the interaction range, as previously noted for θ = π/2 in Fig. 5b, thus occurs continually as θ increases. For the most part, the data for the uniformly distributed (random) case fall closer to θ = π/4 (green triangle-ups). Figure 6b confirms that F y indeed vanishes within the simulation error bars in the special cases of θ = 0 and π/2 over the whole range of intersurface separations. In other cases of θ , including the case of randomly distributed θ , F y 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., d = 0 (surface contact), and examining the detailed behavior of the force components over the θ axis. The aforementioned trends in F x (d) are clearly mirrored by the longitudinal contact force, F x (d = 0) , shown in Fig. 7a. As seen, F x (d = 0) exhibits consistently larger (smaller) repulsive strengths at larger Pe (larger θ ). While, at small Pe , we find a nearly monotonic decrease in F x (d = 0) down to zero, the expected  www.nature.com/scientificreports/ nonmonotonic alterations (due to stronger accumulation of active rods around the inclusions) emerge at large Pe , with a pronounced hump in the subinterval π/3 < θ < π/2 . The transverse contact force F y (d = 0) 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 π/3 < θ < π/2 at around θ ≃ 0.4π. 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, d ∞ , where the interaction force vanishes, into surface contact. We compute this quantity along the x-axis as W = d ∞ 0 F x (d) dx , by taking d ∞ /σ = 10 that gives vanishing interaction force within the simulation error bars. In an equilibrium setting with Pe = 0 , 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 θ ∈ [0, π/2] . As seen, for the parameters in the figure, W in the case with θ = 0 (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 θ ≃ π/6 and then turns up and matches the random sample (dashed line) at around θ ≃ 0.27π , growing to a global maximum at θ ≃ 5π/12 , before falling off as θ tends to π/2 . 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 ( θ = 0 ) and transversally ( θ = π/2 ) self-propelling rods with relative fractions 1 − n π/2 and n π/2 , respectively. W is a convex function of n π/2 , maximized in the monodisperse cases ( n π/2 = 0, 1 ) and a globally minimized at n π/2 ≃ 0.35 . 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 ( θ = 0 ) and transversally ( θ = π/2 ) self-propelling phantom rods, respectively. In these cases, the effective force is again central and F(d) = F(d)x.
Longitudinal self-propulsion (θ = 0). For θ = 0 , 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 interrod 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 ( θ = 0 ) 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, d < σ . The force-distance profiles in Fig. 8a increase up to a maximum www.nature.com/scientificreports/ 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 Pe 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 Pe = 20 and 40 in the figure). As Pe is increased ( Pe 20 ), 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 ( F/F SP ) acting on the inclusions at a fixed level.
Transverse self-propulsion (θ = π/2) : depletion-like attraction. For θ = π/2 , 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 θ = 0 as well (Fig. 8a). The shortranged effective force is found to be attractive in distances less than or comparable to the rod width, d < σ , 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, Pe , 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 d < σ , 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 ( Pe = 0 ). The attractive force behavior as a function of d in Fig. 8b can be estimated analytically, when d is sufficiently small ( d < σ ). We assume that Pe is sufficiently large ( Pe 20 ) 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 ∼ −F SP sin φ , where φ is limited by the gap exclusion as φ ∈ [−π/2 The longitudinal effective force is then estimated as up 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 Pe = 20 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 ( θ = 0 ) and transversally ( θ = π/2 ) 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 σ < d < l (between the width of the rod and its length). This effect will be present only at weak to intermediate values of Pe , as otherwise the strong tendency of the rods to tangentially adhere to the inclusions suppresses such rod-bridging configurations.
Tilted self-propulsion (0 < θ < π/2). Finally, for nonaxial self-propulsion angles 0 < θ < π/2 , 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, F x and F y , as functions of the intersurface distance d are shown in Fig. 9a,b for Pe = 20 and a few different values of θ . The longitudinal force profiles in Fig. 9a display similar features to those with θ = 0 , even though the overall force magnitude decreases as θ increases. On close approach to purely transverse self-propulsion with θ = π/2 , the force magnitudes significantly reduce and the force profiles change to those of Fig. 8b (for illustration purposes, we have excluded θ = π/2 , whose data become indiscernible in the scale of the graph). The transverse force component F y , Fig. 9b, again vanishes for the special cases θ = 0 and π/2 ; 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 F y 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  Fig. 7c, W is convex here and varies mostly monotonically as a function of θ , with a sharp drop as θ approaches π/2 . The work drops negative values in this limit, signifying the attractive depletion-like force profiles of Fig. 8b.

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, Pe . 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 ( θ = π/2 ) as compared with the case of longitudinally self-propelling rods ( θ = 0 ). 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 0 < θ < π/2 ( π/2 < θ < π ). 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) twodimensional 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][51][52][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][82][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).