Effective interactions mediated between two permeable disks in an active fluid

We study steady-state properties of a bath of active Brownian particles (ABPs) in two dimensions in the presence of two fixed, permeable (hollow) disklike inclusions, whose interior and exterior regions can exhibit mismatching motility (self-propulsion) strengths for the ABPs. We show that such a discontinuous motility field strongly affects spatial distribution of ABPs and thus also the effective interaction mediated between the inclusions through the active bath. Such net interactions arise from soft interfacial repulsions between ABPs that sterically interact with and/or pass through permeable membranes assumed to enclose the inclusions. Both regimes of repulsion and attractive (albeit with different mechanisms) are reported and summarized in overall phase diagrams.


Effective interactions mediated between two permeable disks in an active fluid Mahmoud Sebtosheikh * & Ali Naji
We study steady-state properties of a bath of active Brownian particles (ABPs) in two dimensions in the presence of two fixed, permeable (hollow) disklike inclusions, whose interior and exterior regions can exhibit mismatching motility (self-propulsion) strengths for the ABPs. We show that such a discontinuous motility field strongly affects spatial distribution of ABPs and thus also the effective interaction mediated between the inclusions through the active bath. Such net interactions arise from soft interfacial repulsions between ABPs that sterically interact with and/or pass through permeable membranes assumed to enclose the inclusions. Both regimes of repulsion and attractive (albeit with different mechanisms) are reported and summarized in overall phase diagrams.
Effective interactions between colloidal inclusions in fluid media is of crucial importance in predicting phase behaviors of colloidal suspensions and manipulating them for technological applications. Depletion forces are one of the primary types of effective interactions that emerge in size-asymmetric suspensions of (large) inclusions and (small) depletant particles 1,2 . They arise due to steric exclusion of depletants from the proximity of inclusions, are typically short-ranged and attractive.

Model and methods
We consider a two-dimensional minimal model of disklike active Brownian particles (ABPs) 4 with constant self-propulsion speeds and diameter σ in a base fluid, involving two fixed fluid enclosures, or hollow disklike inclusions of effective diameter σ c ; see Fig. 1. Each inclusion has a permeable enclosing membrane of thickness w = σ , modeled by a soft steric potential (see below) that acts within the annulus of inner/outer radii (σ c ∓ w)/2 . The inclusion centers are at R 1 and R 2 along the x-axis at outer surface-to-surface distance d (we label the inclusion on the left/right by k = 1, 2 , respectively).
The ABPs exhibit different self-propulsion speeds, v c and v m , inside and outside the inclusions, respectively, expressed through the discontinuous motility field 58,60 with r = (x, y) being the spatial coordinates and �(·) the Heaviside step function. The ABPs (labeled by i = 1 . . . , N ) are characterized by their position and self-propulsion orientation vectors {r i (t)} = {(x i (t), y i (t))} and {n i (t)} = {(cos θ i (t), sin θ i (t))} , respectively. Their stochastic dynamics are given by the Langevin equations with angular coordinates {θ i } measured from the x-axis, µ T being the translational mobility, U = U({r j }, {R k }) the sum of interaction potentials between the constituent particles (see below), and η i (t) and ζ i (t) being translational and rotational white noises, respectively. The noise distributions are assumed to be Gaussian with zero mean, where D T and D R are the translational and rotational diffusivities, respectively. We adopt the Einstein-Smoluchowski-Sutherland relation, D T = µ T k B T , and Stokes diffusivities for the no-slip spheres, yielding D R = 3D T /σ 265 .
The ABPs are assumed to interact via a Weeks-Chandler-Andersen (WCA) pair potential, giving the energy of interaction between the ith and the jth ABP as www.nature.com/scientificreports/ where r ij = |r i − r j | . As they go through the interfacial regions (membranes) of the inclusions, the ABPs experience a soft, repulsive, WCA-type potential, which we adopt from Ref. 66 , giving the energy of interaction between the ith ABP and the kth inclusion as , and −1 = 1 + (α/σ ′ ) 2 . We use representative values of ǫ = 10k B T and ǫ ′ = 0.0127k B T.
The Langevin equations (2) are numerically solved using standard Brownian Dynamics simulations, and the behavior of the system will be characterized in terms of the key dimensionless parameters given by the interior (c) and exterior (m) Péclet numbers, the rescaled center-to-center distance of the inclusions, d/σ , the rescaled inclusion size, σ c /σ , and the area fraction of ABPs, φ . In the simulations, we use N = 200 − 800 ABPs at fixed area fractions φ = {0.2, 0.4} , with φ = Nπσ 2 /(4L x L y ) defined in the absence of the inclusions. L x and L y are lateral dimensions of the simulation box with periodic boundary conditions. Strong correlations of ABPs with inclusion boundaries suppress motility-induced clustering/phase separation of ABPs 67 that can otherwise arise in the bulk. We fix the inclusion size as σ c = {5σ , 10σ } with L x = L y = {28σ , 56σ } , respectively (we use L x = 2L y = 56σ in Fig. 8b to capture the full range of the interaction force without undesired boundary effects). The area fraction of the inclusions φ c = πσ 2 c /(2L x L y ) = 0.05 (0.025 in Fig. 8b). The simulations are performed for (2 − 6) × 10 7 timesteps of rescaled size D T δt/σ 2 = 1.33 × 10 −5 . The first 10 7 steps are used for relaxation purposes, and the steady-state averages are calculated over 3-20 independent simulations.

Results
Spatial distribution of ABPs. When the inclusions are placed at relatively large surface separations, the spatial distribution of ABPs around each of them exhibits radial symmetry, with a typical radial number density profile, ρ(r) , from the given inclusion center as shown in Fig. 2a for fixed d/σ = 8 , Pe m = 40 and different Pe c . As seen, for Pe c = 0 (nonactive interior; red solid curve), there is a larger concentration of ABPs inside the inclusions relative to the cases with Pe c > 0 (blue and green solid curves), in accord with previous findings in inhomogeneous systems [56][57][58] . The oscillatory behavior of ρ(r) reflects ringlike ABP layers, forming due to steric repulsions between the ABPs inside the inclusions, as depicted also in Fig. 2b. As Pe c is increased, ABPs at the central regions inside the inclusions are more strongly depleted than those near the interfacial regions ( 2 ≤ r/σ ≤ 3 ), as ABPs are slowed down by the enclosing membrane potential and, thus, create larger steadystate densities near the interfacial regions.
The radial symmetry of the ABP distribution is broken when the two inclusions are placed nearby. This is shown in Fig. 2c, where multiple intersections between individual rings formed around each inclusion are illustrated; they are discerned more clearly in the closeup view in Fig. 2d. It is this asymmetric distribution of ABPs in and around the inclusions that causes an effective interaction force between them, which we shall explore in the following section.
Before proceeding further, we briefly examine the internal area fraction of ABPs within the inclusions as a function of Pe c and Pe m ; see Fig. 3. Being defined as φ in = N in σ 2 /(2σ 2 c ) , where N in is the number of particles trapped inside both inclusions, φ in increases by increasing Pe m or decreasing Pe c , which, as noted before, is due to the relatively longer times ABPs spend inside the inclusions. However, φ in cannot exceed a certain value due to steric repulsions between ABPs 58,61 , leading to a finite saturation level (dark red colors in the figure).

Effective interactions and force-distance profiles.
Due to the up-down symmetry of the system, only the x-component (along the center-to-center line) of the effective force acting on the inclusion is nonzero. It can be obtained, upon appropriate steady-state averaging, from the simulations using where the expression inside the brackets is the x-component of the instantaneous force imparted by the ith ABP on the k = 2 inclusion (on the right in Fig. 1) and U (ik) sWCA is defined by Eq. (4). Evidently, F(d) > 0 ( F(d) < 0 ) represents a repulsive (attractive) interaction force. We find identical results for the averaged force on the left inclusion within the simulation error bars. www.nature.com/scientificreports/ Figure 4a shows the effective force between the permeable inclusions as a function of their surface-to-surface distance, d, at fixed exterior Péclet number, Pe m = 20 , and for the interior Péclet number increased from Pe c = 0 up to Pe c = 160 . As a general trend, the force profiles indicate repulsive interactions with a characteristic oscillatory behavior, exhibiting a number of successive local maxima and minima, as the overall magnitude of the force decreases with d. These local minima and maxima can be understood directly based on the overlaps between the rings of ABPs that form around each of the inclusions (and also the intersections of the rings associated with one inclusion with the bounding surface of the other inclusion) as d is varied. These mechanisms have been elucidated in detail in Ref. 18 . As seen in Figure 4b, a similar behavior is found for the effective force, when the interior Péclet number, Pe c , is kept fixed (here, Pe c = 0 ) and Pe m is varied. The dependence of the force profiles in the two above-mentioned cases, however, turns out to be distinctly different, as we shall explore next. Fig. 4a show that, as the interior Péclet number, Pe c , is increased at fixed exterior Péclet number, Pe m , the magnitude of the force acting on the inclusions decreases, converging to a limiting curve already for Pe c = 80 . This behavior is illustrated in Fig. 5a, where we concentrate on the contact force F 0 = F(d = 0) as a function of Pe c at different fixed values of Pe m . The interaction force on the inclusions drops smoothly as Pe c increases. When Pe m is small (see, e.g., the data with Pe m = 25 in Fig. 5a, or Pe m = 20 in Fig. 4a), the force decreases monotonically and tends to a nonvanishing constant as Pe c is increased. When Pe m is sufficiently large (e.g., Pe m = 50 ), the force behavior with Pe c exhibits a weak and broad maximum before it drops to zero. This behavior can be understood as follows.

Role of interior/exterior motility strengths. Our data in
The effective force acting on the inclusions can be viewed as a resultant effect of two separate force components imparted by interior and exterior ABPs on the enclosing membrane of the inclusions, to be referred to as internal and external forces, respectively. To understand the dependence of the effective force on particle motility, one needs to understand how these force components vary with the Péclet numbers. The force components are directly influenced by the concentration of ABPs accumulated near the boundary regions of the inclusions and, specifically, also within their interior regions.
In the cases where Pe m is fixed (Fig. 5a), the external force depends only on the effective hardness of the inclusions and is directly related to the internal ABP concentration. It is nevertheless important to note that in addition to the hardness produced by the enclosing membrane, external ABPs also encounter resistance from In the cases where Pe c is kept fixed (Fig. 5b), there are several factors that control (e.g., enhance or suppress) the effective force imparted on the inclusions. By increasing Pe m , the ABPs present in the exterior regions more deeply overlap with the enclosing membranes. The internal concentration of ABPs increases and this is the factor that enhances the effective force. On the other hand, by increasing Pe m , the average time of interaction between external ABPs and the enclosing membranes decreases and also the ringlike structures of ABPs outside and inside the inclusions are weakened, as the exterior ABPs become more strongly motile and make the interior ABPs more strongly motile due to steric repulsion between ABPs. These latter factors are the ones that suppress the effective force produced on the inclusions. The exact contribution due to each of these factors cannot be determined as they are inter-related and they are not expected to contribute linearly to the effective force. At fixed Pe c = 0 , an initial increase in Pe m gives a dominant increase in the overlapping of external ABPs www.nature.com/scientificreports/ with the enclosing membranes and their internal concentration, and thus an increase in the effective force, but these factor are eventually overtaken by the factors that create effective force suppression. At fixed Pe c = 50 , by increasing Pe m in the range [0, 50], where the internal concentration is diluted (see Fig. 3), there appears to be a competition between the increase in the overlapping of external ABPs with the enclosing membranes and the decrease in the average time of the interaction between external ABPs and the enclosing membranes. As Pe m is increased, the former factor becomes dominant and, as a result, the effective force increases but, on further increase of Pe m , the latter factor dominates and, hence, the effective force decreases. In the range Pe m > 50 , the internal ABP concentration increases with Pe m ; this, being the dominant factor in this regime of parameters, causes an increase in both internal and external forces and, as a result, an increase in the effective force. The increase in the effective force continues until a secondary hump appears, as seen in Fig. 5b (blue curve). The said increase stops after the internal concentration is saturated, in which case, the decrease in the interaction time becomes dominant again and, thus, the effective force decreases. At Pe c = 100 , the same mechanisms hold except that the second hump does not appear in the range Pe m < 200 , because the internal concentration only weakly increases with Pe m .  www.nature.com/scientificreports/ Figure 6b,c show that a net attraction originates in the internal component, F in 0 , of the force acting on the inclusions, as F ext 0 never becomes attractive. Panel b also indicates that F in 0 is attractive typically when Pe c 60 , unless Pe m goes to zero, in which case the threshold Pe c (above which F in 0 becomes attractive) also tends to zero. Further insight can be obtained by comparing the regions of attraction and repulsion in Fig. 6a-c with the regions of the parameter space, where ABPs are found to be depleted from or saturated within the inclusions; see Fig. 3. Hence, it follows immediately that a net repulsion occurs when the ABPs are predominantly depleted from or saturated within the inclusions.

Regimes of attraction and repulsion. Even
The signs of the internal and external forces depend on the distribution of ABPs interacting with the enclosing membranes of the inclusions in their interior and exterior regions. Based on our simulation results for the ABP distributions and the resulting force components, we can divide the (Pe m , Pe c ) plane to three parametric regions I-III, as shown in Fig. 6d. The boundary lines here are determined by numerical interpolation of data obtained by fixing Pe c ( Pe m ) and scanning the Pe m ( Pe c ) axis at the resolution of P m = 5 ( P c = 5 ). In all of these three parametric regions, ABPs are more strongly concentrated at the narrow and wedge-shaped exterior region, intervening the two inclusions, where ABP trapping effects are predominant. This is the primary source of the all-repulsive external component of the force F ext 0 . However, as schematically shown in Fig. 7, in region I, the ABPs that are found within the inclusions show stronger steric overlaps with the enclosing membranes at the distal parts of the interior region (the ABPs found at the proximity of the intervening wedge-shaped gap between the inclusions show weaker steric overlaps with the enclosing membranes mainly because the inside ABPs also experience stronger inward-pointing repulsions from the outside ABPs that are accumulated in the outside inter-inclusion gap). As such, the inside ABPs tend to push the inclusions apart, making the internal force, F in 0 , repulsive as well. Region I thus corresponds to the situation, where F in 0 , F ext 0 and F 0 > 0 . In regions II and III, the ABPs that are found within the inclusions are more strongly concentrated at the proximity of the inter-inclusion gap and induce an attractive internal force. In region II, because of low concentration or low motility of ABPs inside the inclusions, such an attractive internal force cannot dominate the repulsive force due to external ABPs and, as a result, the total effective force turns out to be repulsive. This region corresponds to the situation, where F in 0 < 0 , F ext 0 > 0 (that is, |F in 0 | < |F ext 0 | ) and F 0 > 0 . In region III, the inside ABPs acquire sufficient concentration and motility to induce attractive internal force to dominate the repulsive external force.   Fig. 3). In region III, F in 0 < 0 , F ext 0 > 0 (with |F in 0 | > |F ext 0 | ) and F 0 < 0.

Role of inclusion size and ABP area fraction.
The general behavior of the force profiles with the surface-to-surface distance between the inclusions remains the same as the inclusion size is increased; see Fig. 8a. The larger the inclusion size the larger will be the magnitude and the range of the force experienced from the ABPs. This is expected because the increased perimeter of the enclosing membrane of the inclusions and also its reduced curvature lead to larger residence times for the ABPs near the inclusions and, hence, larger proximal ABP concentrations 29 ; this is further supported by the effective force divided by the perimeter of inclusions in the inset of Fig. 8a.
In systems with larger area fractions, a relatively larger fraction of ABPs accumulate inside and around the inclusions, which in effect leads to a larger and longer-ranged force profile as shown Fig. 8b. In this case, a larger number of ABP rings form around the inclusions as reflected by the larger number of peaks in the force profiles.

Concluding remarks
We investigate spatial distribution of ABPs and effective interactions mediated by them between permeable (hollow) dislike inclusions with mismatching ABP motility (self-propulsion) strengths in their interior/exterior regions. The inclusions are enclosed in permeable membranes modeled using soft, repulsive interfacial potentials. The ABPs are shown to distribute in ringlike high-density layers inside and outside the inclusions, which strongly influence the effective interactions mediated between the inclusions through the active bath. We determine the regimes of effective repulsive/attractive and elucidate their underlying mechanisms by decomposing the net force between the inclusions to external and internal force components; while the external force turns out to be repulsive, the internal force is found to be repulsive (attractive) below (above) a certain threshold interior Péclet number, indicating that a net attraction stems from a dominant internal force.
Our study builds on a basic model of ABPs interacting only through steric repulsions and demonstrates how a relatively complex nonequilibrium phenomenology emerges from a minimal description of inclusion permeability and spatial inhomogeneity in ABP motility. Our model can be generalized to include other types of spatial inhomogeneity in the system such as mismatching interior/exterior values for particle diffusivity (temperature), fluid viscosity (especially when hydrodynamic stresses are relevant), and stimulus/field intensities. The latter can be realized by subjecting the interior/exterior regions to external stimuli (e.g., activating light 59,68 ) or fields (e.g., magnetic/gravitational fields) that impart forces and/or torques of differing magnitudes on the ABPs inside and outside the inclusions. In all such cases, however, particle activity is expected to emerge as a crucial factor. Indeed, switching off the self-propulsion mechanism in the present model also leads to complete elimination of the phenomenology reported.
Our model assumes that the inclusions are placed at fixed positions within the active bath, an assumption that has frequently been used in the literature 8,9,[11][12][13][14][15][16]18,20,21,[23][24][25] , as it facilitates direct analogies with textbook examples of equilibrium depletion forces between colloidal inclusions fixed within a bath of nonactive depletants 1,2 , where the bath-mediated interactions represent the corresponding potentials of mean force between the inclusions. Several studies involving mobile (hard) inclusions in an active bath have also been published 10,11,19,22,24,26 . Effective interactions obtained in the case of fixed inclusions are expected to be different from those obtained in the case of mobile inclusions in an active bath in a similar setting 24 , and possible relations between the results in the two cases remain to be understood. It is thus interesting to study the role of inclusion mobility in the context of permeable inclusions as well, where many-body effects 18 can also be examined by assuming a suspension of many such inclusions within the active bath.
Possible real-world examples of permeable inclusions may be furnished by fluid enclosures such as lipid and polymer vesicles [30][31][32] , cells, soft tissues as well as immiscible/stabilized emulsion drops and active droplets [33][34][35][36][37]55 . Polymersomes (polymer vesicles) fabricated with controlled permeability 41 and polymer stomatocytes derived from controlled deformation of polymersomes 39 have recently been used to selectively entrap catalytically active platinum nanoparticles. In the biological context, cells can ingest and internalize foreign particles through phagocytosis; a process recently utilized in controlled internalization and cell-assisted assembly of nonactive polystyrene microparticles and crystallites within fibroblast cells 38 . Our study should motivate further studies along these lines using active particles as phagocytosed components. Active particles have also emerged as efficient agents for targeted drug and cargo delivery [40][41][42][43][44][45] . Soft biological tissues such as tumors can exhibit enhanced permeation to suitably designed active particles such as platinum-sputtered polymersome nanomotors 40 , and biological microswimmers such as magneto-aerotactic bacteria 45,46 .
Recent studies of semipermeable vesicles filled with active particles 37,[49][50][51][52][53][54] can also be generalized straightforwardly to incorporate permeation of active particles through the enclosing vesicle membrane 47,48 . Such extensions can provide insight into the interplay between permeability and shape fluctuations/deformability, brining the model vesicles closer to realistic examples of soft permeable inclusions.