Phase separation of active Brownian particles on curved surfaces

The effect of curvature on an ensemble of repulsive active Brownian particles (ABPs) moving on a spherical surface is studied. Surface curvature strongly affects the dynamics of ABPs, as it introduces a new time scale $\tau=R/v_0$, with curvature radius $R$ and propulsion velocity $v_0$, in addition to the rotational diffusion time $\tau_r$. This implies that motility-induced phase separation (MIPS) disappears for small $R$. Furthermore, it causes a narrowing of the MIPS regime in the phase diagram of P{\'e}clet number $\text{Pe}$ and particle area fraction $\phi$. Also, the phase-separation boundary at low $\phi$ attains a turning point at small $R$, allowing for the possibility of a reentrant behavior. These results characterize the effect of curvature on ABP dynamics and MIPS, and will help to better understand the preferred occupation of certain niches by bacterial colonies in porous media.

In biophysical systems, active particles are often exposed to curved geometries and confinement.Examples include bacteria motion in porous media [27], cell migration on curved tissues of the gut [28], embryonic development [29], and actomyosin flows during cell division [30].Theoretical studies of single active particles indicate that for a tangential propulsion direction, their dynamics depends on the surface curvature [31,32], while for an unconstrained propulsion direction, particles predominantly accumulate in regions of higher curvature [33,34].Steric interactions among active elongated particles, which favor polar or nematic alignment, generate complex flow patters on spherical geometries, such as circulating band states [35][36][37].Further studies have shown topology-dependent collective dynamics of self-propelled rods [38], as well as segregation dynamics in binary mixtures of active and passive particles on spherical surfaces [39].
In this work, we study the dynamics and clustering behavior of repulsive ABPs, with a freely diffusing propulsion direction, but constrained to move on a curved surface in two/three spatial dimensions (2D/3D).We show that the confinement radius R introduces a new length scale that changes qualitatively the dynamics of ABPs on * Theoretical Physics of Living Matter, Institute of Biological Information Processing and Institute for Advanced Simulation, Forschungszentrum Jülich, 52425 Jülich, Germany Email: p.iyer@fz-juelich.de, r.winkler@fz-juelich.de, d.fedosov@fz-juelich.de, g.gompper@fz-juelich.de curved surfaces, such that their ballistic motion for large Péclet numbers is suppressed at distances smaller than R. The Péclet number is defined as Pe = v 0 τ r /σ, where v 0 is the propulsion velocity, τ r = D −1 r is the rotational diffusion time with a rotational diffusion D r , and σ is the ABP diameter.The diagram of phase separation on a sphere is constructed, and shows that curvature drastically changes the phase boundaries, and completely suppresses MIPS at small R. A simple model that considers the effective persistence length of particle motion on the sphere is then used to rationalize the observed effects of sphere curvature on MIPS.Furthermore, we study MIPS in a paradigmatic example of porous media represented by two connected spheres with unequal radii.
We consider first a model of an ABP as a disc of diameter σ in 2D, confined to a ring of radius R. While the ABP motion is restricted to one dimension, the propulsion vector e is free to rotate in 2D, see Fig. 1(a).We neglect the effects of translational noise and focus on rotational noise.Then, the equations of motion for the position r = (R cos θ, R sin θ) and propulsion direction e = (cos ψ, sin ψ) of the ABP are where Γ ψ is a Gaussian and Markovian random process with zero mean and Γ ψ (t)Γ ψ (t ) = δ(t−t ).In the limit of small misalignment angles δ ≡ θ − ψ with |δ| 1, Eq. ( 1) can be linearized and analytically solved [40], yielding the angular mean-squared displacement (MSD) Thus, irrespective of τ r , the crossover between the diffusive and ballistic regime is determined by the new time scale τ = R/v 0 .The shift of the onset of the diffusive behavior to earlier times for smaller radii R originates from a fast alignment of particle orientation with the surface normal r, as the ABP moves along the surface, upon which the translational motion of the particle nearly stops.Note that the MSD from Eq. ( 1) agrees well with simulation results for a single ABP on a ring [40].
The discussion above suggests that the misalignment angle δ plays an essential role.The corresponding Fokker-Planck equation [40] yields the stationary-state distribution P (δ) ∼ exp[(τ r /τ ) cos(δ)], from which we obtain δ 2 = τ /τ r = R/(σPe) for τ /τ r 1.For the tangential velocity v = v 0 | sin(δ)|, this implies that v v 0 τ /τ r .Thus, the ABP velocity slows down with increasing curvature (or decreasing τ ), which is also reflected in Eq. ( 2) where the MSD in the ballistic regime is given by R The behavior of ABPs and the corresponding MSD in 3D are more complex, as both time scales τ and τ r become important due to additional angular degrees of freedom.Figure 1(b) shows that the ratio τ r /τ = (σPe)/R determines particle dynamics.For τ r /τ 1, the ABP does not 'see' the effect of curvature and exhibits diffusive motion for times larger than τ r .However, when τ r /τ 1, the particle moves ballistically only up to time τ .For t > τ , diffusive motion due to sphere curvature sets in, and can be described as a stop-and-go motion due to fast alignment of e along r. Figure 1(b) also demonstrates the reduction of effective particle velocity with increasing Pe as the magnitude of the MSD in the ballistic regime drops.
We consider next an ensemble of N ABPs on a sphere with area packing fraction φ, to study how MIPS is affected by the sphere curvature 1/R.In the simulations [40], Pe is changed by varying τ r while v 0 is kept fixed.Figure 2 shows simulation snapshots for two different curvatures at Pe = 890, and demonstrates the absence of MIPS for the small radius R/σ = 16.1.The full phase diagram for different R values is presented in Fig. 3(a).Here, the binodal is constructed by measuring co-existing densities in the phase separated state, whereas the spinodal is obtained by computing the particle pressure [40,41].A sudden drop/change in pressure marks the transition from the homogeneous to the phaseseparated state [41][42][43][44].Two main effects of curvature can be seen in Fig. 3(a) for decreasing R: (i) the lower part of binodals and spinodals shifts to larger Pe, and (ii) the two-phase region becomes narrower and the slope of the left spinodal/binodal changes sign for large Pe. Figure 3(b) shows the variation of the critical Péclet number Pe c (R), where MIPS is first observed (binodal) with increasing Pe for a fixed initial density of φ = 0.5.Since the value of R sets the total number N of particles for a given φ, the system may inherently suffer from finite-size effects.To account for the possible effects of finite N , planar simulations with the same number of ABPs are performed, see Fig. 3(b).Pe c increases with decreasing N in both cases, see Fig. 3(b).However, this effect is much less pronounced for the planar systems, demonstrating that finite-size effects are sub-dominant.Furthermore, the width of the MIPS region becomes narrower with increasing Pe [see Fig. 3(a)], which explains the sudden disappearance of MIPS at R/σ ≈ 12.5 and φ = 0.5 in Fig. 3(b) (i.e., the spinodal has a turning point before φ = 0.5 is reached for R/σ 12.5).This clearly implies that the loss of MIPS at R/σ ≈ 12.5 is not due to finite-size effects.
MIPS occurs as a result of slowing down of ABPs due to crowding, which promotes a further reduction in velocity and clustering through a positive feedback mechanism [10,45].A requirement for MIPS is that the life time of small clusters is larger than the persistent travel time of ABPs [21,46].This means that the directed selfpropelled motion should dominate over diffusive motion on the length scale of particle diameter σ, i.e. σ/τ r v 0 or Pe = v 0 τ r /σ 1.For ABP motion on a curved surface, this argument has to be modified as follows.First, the propulsion velocity v 0 has to be replaced by a radiusdependent velocity v(R).In general, v(R) decreases with decreasing R, e.g., in 2D, v(R) = v 0 τ /τ r = v 0 R/Pe for R/Pe < 1.Second, on a curved surface, the time scale τ becomes relevant in addition to τ r .Thus, we have to distinguish the two cases ατ > τ r and ατ < τ r (α is a constant of order unity), which represent large and small radii R, respectively.In both cases, the shorter time scale τ min = min(ατ, τ r ) determines the dynamics [see Figure 1(b)].As a result, we can define an effective Péclet number Pe ef f = v(R)τ min /σ, which has to exceed the threshold Pe c for phase separation to occur.Hence, a larger bare Pe is required for a smaller R to compensate for the reduced effective surface velocity.
From this argument, all the trends observed in Fig. 3 can be understood.With decreasing R, Pe c first increases, because τ r is the relevant time scale and v(R) decreases, and this increase has to be compensated by a larger Pe.Note that Pe ef f increases with increasing τ r and v 0 only as long as τ r < ατ .When τ r = ατ , Pe ef f reaches a maximum as a function of Pe for a fixed R. A further increase of Pe only causes a decrease in the effective surface velocity v(R), without any increase in τ min = ατ .This leads to a decrease in Pe ef f and the turning of the low-φ branch of two-phase coexistence toward larger φ values in Fig. 3(a).
At low-to-intermediate particle densities (φ < 0.5) and smaller radii, with τ r > ατ , MIPS is absent for all Pe.This inversion of time scales and disappearance of MIPS occur when particle diffusion dominates over the mini-mum run length for cluster formation.In this case, an increase of τ r cannot lead to MIPS, because the slowing down due to translational ABP motion on a curved surface always precedes rotational diffusion.Furthermore, this argument indicates that for smaller R values, the occurrence of MIPS requires larger area fractions φ because Pe c is larger, as shown in Fig. 3(a) where MIPS for R/σ = 16.1 occurs only at φ 0.35.This curvature effect further lowers the range of R where MIPS is observed.
To verify that the effective Péclet number Pe ef f indeed controls phase separation on curved surfaces, average surface velocity v(R) for a single ABP is measured in simulations.Figure 3(c) shows a heat map of Pe ef f for various radii, where α = 6 is selected for a good fit of the simulation data for MIPS.Lower and upper boundaries of the MIPS region for a fixed φ = 0.5 agree well with the black dashed line for Pe ef f = 25.Furthermore, the heat map of Pe ef f nicely explains the loss of MIPS at small R. Figure 3(c) also shows that for a larger critical Pe c (or smaller R), the MIPS regime becomes narrower as a function of R, consistent with the onset of MIPS at larger φ.Noteworthy, Pe ef f reaches a maximum and then decreases as a function of Pe, which explains the turning of the phase boundary at large Pe.This supports the existence of a reentrant behavior (i.e., from homogeneous to MIPS and back to homogeneous density) with increasing Pe for a wide range of radii.Note that we have not observed any significant change in the right binodal/spinodal at large φ with curvature.This is due to high densities, at which inter-particle collisions are very frequent, so that the particle dynamics is significantly affected and the simple estimate based on single-particle Pe ef f is not valid.Following the characterization of the behavior of ABPs at a surface with a fixed curvature, we make a further step toward the understanding of their behavior in porous media.We consider a paradigmatic example of two spherical pores with unequal radii, which are connected by a small passage as shown in Fig. 4(a).For low area densities of ABPs, the solution of the Langevin equation in a convex, non-spherical confinement in 3D shows that the single-particle density is proportional to the squared local curvature of the boundary [33,34].For two connected pores in 3D, this implies that the ratio of particle area fractions φ 1 and φ 2 within the larger and smaller spheres with radii R 1 and R 2 is given by φ 1 /φ 2 = (R 2 /R 1 ) 2 .Therefore, for non-interacting ABPs at the steady state, their area fraction in the smaller pore is larger than that in the larger pore, while the average number N i of particles in each pore i = 1, 2 should be the same, i.e.N 1 = N 2 = N/2.As the (equilibrium) particle pressure with excluded-volume interactions increases faster with area fraction than for an ideal gas, the particle number N 2 in the smaller pore should be smaller than N/2, and correspondingly N 1 > N/2.Then, our results in Fig. 3(a) for MIPS in a single pore allow for predictions of the steady-state behavior of ABPs in the two pores.As N is increased for a given Pe, the particles in either sphere can phase separate only when their surface density exceeds φ c (R, Pe) for MIPS.The density in the smaller pore, which fills first, increases linearly with N , and MIPS is expected at N = 2φ c (R 2 , Pe)A 2 /A σ , where A 2 is the area of the smaller sphere minus the area of the passage between two pores, and A σ is the area occupied by an ABP.The number of particles in the smaller pore cannot exceed N 2,max = φ cp A 2 /A σ , where φ cp is the close-packing density.Therefore, the number of particles in the larger pore at the steady state is expected to be N 1 (N ) = N − min(N/2, N 2,max ).As a result, for a given Pe, phase separation in the larger sphere is first expected to occur near the value of N that satisfies the equality N 1 (N ) = φ c (R 1 , Pe)A 1 /A σ (A 1 is the area of the larger sphere minus the passage area).
To test our predictions, we consider two connected pores with R 1 /σ = 16.1 and R 2 /R 1 = 0.6, where MIPS in the larger pore is expected at N 2900 and N 3300 for Pe = 90 and Pe = 890, respectively.These predictions are nicely confirmed by the simulation results in Fig. 4(b), where the onset of phase separation as a function of N is characterized by a sudden rise in the fraction n cl of particles occupying clusters of size greater than N 1 /2.Note that the reduction of the effective Péclet number at large Pe implies that a larger N is required for MIPS.
Another interesting observation is the temporal evolution of this system after starting with a uniform equilibrium distribution of ABPs and approaching the transition regime.For sufficiently large N , the larger pore shows a 'dynamic' MIPS state.At short times, the initial area fraction is large enough to show MIPS, which eventually dissolves as ABPs are lost to the smaller pore with time.For low Péclet numbers, due to large fluctuations in the particle number in the larger sphere, the MIPS state can both be dynamically restored and lost in time [40].
In summary, if the propulsion direction of ABPs can vary diffusively (i.e., it is not aligned with the local tangent plane of a surface), the behavior of active particles on curved surfaces is very different in comparison to ABPs confined to a plane.Non-zero curvature results in a 'stop-and-go' motion, such that particles slide along the surface when their orientation is different from the local normal, and then stop after their orientation becomes perpendicular to the surface.This behavior governs motility-induced phase separation on curved surfaces, e.g., the MIPS region rapidly shrinks with increasing curvature and eventually disappears.Furthermore, curved surfaces lead to a possible reentrant behavior, where MIPS for a fixed surface density of ABPs first appears with increasing Pe, and then can disappear.The single-pore results also allow us to predict the dynamics of active particles in connected pores with distinct curvatures as in porous media.These results will help to better understand the preferred occupation of certain geometries and niches by bacterial colonies [47].

FIG. 1 :
FIG. 1: (a) Schematic diagram of an ABP confined to a ring.(b) MSD of a particle moving on a sphere for different Pe (or τ r ), with v 0 fixed.For τ r /τ < 1, the ballistic-to-diffusive transition occurs at time τ r /τ , whereas for τ r /τ > 1, the time scale τ determines the ballistic-to-diffusive transition irrespective of τ r .

FIG. 2 :
FIG. 2: Simulation snapshots for different sphere radii at Pe = 890 and φ = 0.5.(a) R/σ = 26.8,N = 5760, and (b) R/σ = 16.1,N = 2074.A planar system with the same number of particles as in (b) exhibits phase separation [40], indicating that the absence of MIPS for R/σ = 16.1 is not due to finite-size effects.See also movies S1 and S2.

FIG. 3 :
FIG. 3: (a) Pe-φ phase diagrams of motility-induced phase separation (MIPS) for three R values.Coexisting densities from the local density distributions (circles) and abrupt pressure drops (squares) are employed to identify the transition.The simulations for determining the coexisting densities are performed for an average area fraction of φ = 0.5.(b) Critical Péclet number Pe c (R) at which MIPS is first observed for increasing Pe at fixed φ = 0.5.The color-bar shows the variance of the local density distribution.Symbols mark the identification of no-MIPS (uni-modal P (φ loc )) (circles) and MIPS (bi-modal P (φ loc )) (triangles).The black line with bullets is obtained from the threshold σ 2 φ = 0.0305 of variance of the local density, and it follows well the boundary where the two peaks in the local density distribution merge.(c) Heat map of the effective Péclet number Pe ef f for a single ABP at the sphere surface as a function of Pe and R. The black and red dashed lines represent Pe ef f = 25 and Pe ef f = 35, and match well the lower and upper MIPS boundaries from the simulations.

FIG. 4 :
FIG. 4: (a) Snapshot of a system of two interconnected pores of radii R 1 /σ = 16.1 and R 2 = 0.6R 1 for N = 2500 at Pe = 90 (see also movie S3).(b) Fraction n cl of particles in the large sphere occupying clusters of size greater than N 1 /2 as a function of N .The sudden jump identifies MIPS.