Interacting particles in an activity landscape

We study interacting active Brownian particles (ABPs) with a space-dependent swim velocity via simulation and theory. We find that, although an equation of state exists, a mechanical equilibrium does not apply to ABPs in activity landscapes. The pressure difference originates in the flux of polar order and the gradient of swim velocity across the interface between regions of different activity. In contrast to motility-induced phase separation of ABPs with a homogeneous swim velocity, a critical point does not exist for an active-passive patch system, which continuously splits into a dense and a dilute phase with increasing activity. However, if the global density is so high that not all particles can be packed onto the inactive patch, then MIPS-like behavior is restored and the pressure is balanced again.


Introduction
Active particles, like motile microorganisms, live in a complex environment and are affected by various external fields like gravity, fluid flows or walls [1]. Mostly it is assumed that their intrinsic properties, like propulsion speed or tumbling rate, are constant in space [2]. This is generally not true because organisms often respond to a spatially varying stimulus such as light intensity or chemical concentration. For example, the bacteria Escherichia coli are chemotactic: they sense concentration gradients and adjust their tumbling rate to swim up or down a chemical gradient. Others do not respond to gradients but react in a non-directional way to the local stimulus intensity, a behaviour called kinesis. For instance, some cyanobacteria are photokinetic and their swim speed depends on the local light intensity [3]. In this context, it has been predicted theoretically that the local density of particles ρ(r) performing a run-and-tumble motion is inversely proportional to their local propulsion speed v(r) [4]. This fact was used to arrange millions of light-powered bacteria into a complex pattern, such as Leonardo da Vinci's Mona Lisa, via light fields [5,6]. Optical fields can also be used to create activity landscapes for synthetic microswimmers, such as thermophoretic Janus particles [7]. An activity landscape can have a technological application, such as traps [8], but it can also be used to study fundamental questions of active matter [9]. One such question is whether there is a coexistence criterion for active matter. The introduction of an active component of pressure [10,11], often called "swim pressure", allowed to formulate a mechanical equilibrium condition of equal pressures in phase separating purely repulsive active Brownian particles (ABPs) with a homogeneous propulsion speed [12]. The separation into a dense-inactive and a dilute-active phase is triggered by a slowdown during collisions and is named motility-induced phase separation (MIPS) [13]. The behaviour of ABPs in activity landscapes, like the aggregation of particles in slow regions, resembles in some sense MIPS, however, the differences and similarities are not well studied. Here we show that an activity pattern can still lead to pressure imbalance between regions of different activity even though an equation of state exists for ABPs [12]. We show that, in contrast to MIPS, a critical point does not exist for an active-passive patch system, which continuously splits into a dense and a dilute phase with increasing activity [14,15]. However, if the global density is so high that not all particles can be packed onto the inactive patch, then MIPS-like behaviour is restored and the pressure is balanced again. The particles in the low activity region exhibit all states of matter of a two dimensional system: a liquid, a hexatic phase and a crystal. The transition to a state with a high orientational order for a sufficiently dense and active system can be predicted theoretically.

Model
We consider N active Brownian particles (ABPs) in 2D with a space-dependent propulsion velocity v(r). The particles interact via a repulsive pair potential V (r) = k 2 (σ − r) 2 if r ≤ σ, i.e., the inter-particle distance r is smaller then the particle diameter σ = 2a, and V (r) = 0 otherwise [16]. The repulsion strength k is chosen such that the particle overlap is around 0.01σ. The positions r i = (x i , y i ) and the orientations e i = (cos θ i , sin θ i ) evolve according to the overdamped Langevin equations: where 1/µ t = γ t is the translational friction coefficient, is the force on i-th particle due to j-th particle, D t and D r are the translational and rotational diffusion constant and η i , ξ i are zero-mean unit-variance Gaussian white noises. For spherical Brownian particles it is D r = 3D t /σ 2 . We use a simulation box of size L x /σ = 160 and L y /σ = 40 with periodic boundary conditions in both directions and consider mainly a step-like activity landscape where Θ denotes the unit step function and α controls the extent of the active region (α = 0.5 means that a half of the box is active). Our length and time unit is the particle radius a and the rotational diffusion time τ r = 1/D r , respectively. Thus, the Péclet number, P e = v a /(aD r ), is a dimensionless measure of activity and φ 0 = πa 2 N/(L x L y ) is the global packing fraction.

Theory
First we ask whether a mechanical equilibrium exists or in other words: are pressures equal in neighboring regions of different mobility? Starting from (1) and (2) for the time evolution of the noise-averaged probability density ψ(r, where . . . denotes noise averages andρ = dθψ the fluctuating particle density. In the following we assume that the system is in the steady state, ∂ t ψ = 0, and is only inhomogeneous along the x-direction. Integrating (4) over θ gives with the particle flux where m = dθ cos(θ)ψ is the polarisation and I 1 = dr µ t f x (r − r) ρ(r )ρ(r) is a contribution due to the pair-potential. Similarly, multiplying (4) by cos(θ) and integrating over θ gives with the flux of polar order where Q = 1 2 dθ cos(2θ)ψ encodes nematic order along x-direction and I 2 = dr µ t f x (r − r) ρ(r )m(r) contains the effect of interactions. All setups considered here are flux-free steady states, i.e., J ρ (x) = 0. Thus inserting (7) into the vanishing particle flux J ρ in (6) yields Let us next consider two neighboring domains A and B with constant activity v A and v B , respectively. Then multiplying (9) with γ t = k B T /D t , integrating from the bulk of region A to the bulk of region B and using integration by parts gives We define the ideal passive pressure as p id = ρk B T , the interaction or "direct" pressure as ∂ x p D = −γ t I 1 [17,20] and the "swim" pressure as p S = γ t v 2 ρ/2D r +γ t vI 2 /D r [10,17,21]. For a homogeneous bulk phase in region A and B we can set Q = ∂ x m = 0 and obtain finally which is the main result of our study. Especially, for an activity step, where an ABP swims with v A for x < x if and abruptly propels with v B when crossing the interface at What does this mean? It means that the bulk pressures in both activity regions are not equal and that the pressure difference is governed by the activity profile v(x if ) and the flux of polar order J m (x if ) at the interface between this regions. Thus an inhomogeneous activity pattern can still lead to unequal pressures in both phases even though an equation of state exists for ABPs [12]. A similar expression as (11) have been obtained for underdamped quorum-sensing active particles undergoing motility-induced phase separation [22], which lack an equation of state for the pressure even for spatially uniform activity [12]. Next, we examine the case of non-interacting ABP's, which can be solved exactly in all details [9,23]. In particular, let us consider a setup where particles are inactive (v p = 0) for x < x if = 0 and are active (v a > 0) for x > 0. We non-dimensionalize (12) using πa 2 /k B T . The left hand side of (12) reads as using the bulk density ratio λ a /λ p = ρ a /ρ p [9,14,15,23,24], where λ α = (D r /D t + v 2 α /2D 2 t ) −1/2 denotes the polarisation decay length and φ = πa 2 ρ the packing fraction. In the following we use the polarisation profile obtained in [9,23], which in the passive region ( for our particular setup. In order to calculate the right hand side of Eq. (12) we use where J m (−∞) = 0 because of v(−∞) = ∂ x m(−∞) = 0, see (8) for the individual terms of J m . In total we get for the right hand side of (12) which is equivalent to (13) and thus proves the validity of (12) for the ideal system. Due to λ p > λ a , the pressure in the active region is larger than the pressure in the passive region.
With (12) we want to estimate the densities φ a and φ p of interacting particles in an active-passive patch system. An approximate equation of state for interacting ABP's in 2D [25] reads as where φ max = π/ √ 12 is the packing fraction of a close-packing of monodisperse disks. The "swim" pressure vanishes and the interaction pressure diverges at φ max , respectively. For passive hard disks we use [26] which, however, does not capture the existence of a "hexatic" phase at φ ∈ [0.7, 0.72], see [27]. We obtain one equation by inserting (16) and (17) into (12) using the interface contribution of non-interacting particles (15), however, with a density dependent activity v(φ) = v a (1 − φ/φ max ). And we get another equation by assuming a rectangular density profile such that φ 0 = αφ a + (1 − α)φ p . The numerical solution of these two equations give the densities φ a and φ p for each P e and φ 0 . Figure 1(a,b) shows snapshots of a system with a step-like activity pattern with α = 0.5, i.e., half of the box, namely [−L x /4, L x /4], is active and the rest is passive. In Figure 1(a,c,e) the global packing fraction is φ 0 = 0.4 < (1 − α) φ max ≈ 0.4534 and the maximum possible packing fraction in the passive patch is below close-packing. Figure 1(a) demonstrates the general behaviour of a system with a position dependent activity v(r), namely, that particles tend to accumulate in the less mobile region or in other words, ρ ∝ 1/v [4,28]. A strong positive polarization m (particles point towards the passive region) appears solely at the active-passive interface and is zero otherwise, see figure 1(c). Similarly to the behavior of ABP's near walls [29], only the particles that point towards the interface can also approach it and once they cross the interface they keep their orientation for τ r . In Figure 1(e) we show the normal component of the local "swim" p S (x) [30] and interaction pressure tensor p D (x) [31]. Surprisingly, the pressure is not equal in both regions, even though an equation of state for ABPs exists [12], and is larger in the dilute active region, where the "swim" pressure p S is the dominant contribution. There is a pressure gradient, but we do not observe particle flux. How to resolve this contradiction? One can take two points of view. If one argues that the "true" pressure consists only of the ideal p id and the "direct" pressure p D (neglecting the "swim" pressure p S ) then the pressure gradient is balanced by the swim force density γ t vm created by the polarization of the active particles [32] or in other words which follows from J ρ (x) = 0. For our case that means that the dense passive phase is held together by the rim of polarized particles in the active region. If one also considers p S then the pressure gradient is balanced by the term on the right hand side of (11) consisting of the flux of polar order J m and the gradient in the activity ∂ x v. For a dense system, φ 0 > (1 − α) φ max , the situation is different and is reminiscent of so-called motility-induced phase separation [16,17,25], as can be seen in figure 1(b,d,f). The dense phase extends into the active region and the polarized region is no longer located at the active-passive boundary but at the interface between the dense and the dilute phase. Moreover, this interface fluctuates heavily and the pressure is constant through the system. ABPs with a homogeneous propulsion speed separate only above a critical point P e cr ≈ 26.7, φ cr ≈ 0.597 into a dense-immobile and a dilute-mobile phase [33]. The densities of the coexisting phases are independent of the global density φ 0 at fixed activity P e. The situation is completely different for systems with a mobility pattern. In figure 2 the packing fractions φ a and φ p as a function of P e for two different φ 0 below overcrowding, φ 0 < (1 − α) φ max ≈ 0.4534, are shown. No signs of a critical point are visible, i.e., the system splits continuously, starting at P e = 0, into a dense and a dilute phase with increasing P e. In contrast to the ideal case, the densities φ a and φ p depend on the global density φ 0 at fixed P e and in the limit P e → ∞ all particles occupy the passive patch, φ p = φ 0 /(1 − α) and φ a = 0, as long as φ 0 < (1 − α) φ max . Moreover, figure 2 demonstrates that the theory, based on (12), agree very well with the simulation results for densities below overcrowding. Figure 3 indicates the interdependencies between the densities of the coexisting phases, the corresponding pressures and the polarisation at the interface between both phases. Independent of the global density, the density ratio behaves as for ideal particles, φ p /φ a ∝ P e, see figure 3(a), however, φ p /φ a decreases with φ 0 at fixed P e due to excluded volume effects as is evident from figure 3(b). The pressure ratio behaves similarly to the non-interacting case, p p /p a ∝ 1/P e, for φ 0 < (1 − α) φ max ≈ 0.4534, see figure 3(c). However, the pressure in both phases becomes equal, p a = p p , for the overcrowded case, φ 0 > (1 − α) φ max , see figure 3(d). The total interface polarisation m tot = p a dx m(x), which is proportional to the flux of polar order J m in the bulk of the active region [34], saturates with increasing P e, see figure 3(e). Besides the trivial dependency m tot ∝ φ 0 there is a pronounced change in growth of m tot as a function of φ 0 at (1 − α) φ max for large P e, see figure 3(f). Once the passive patch is fully occupied the dense phase have to extend into the active region, which means that the position of the interface between the dense and the dilute phase x if does not coincide anymore with the boundary of the active and passive patch αL x /2, see figure 4(a). Consequently, the interface start to fluctuate strongly and its width λ increases significantly above (1 − α) φ max , see figure 4(b), resembling the behaviour of an interface during motilityinduced phase separation [35]. We obtained x if and λ by fitting the density profile with a hyperbolic tangent.

Results
As clearly visible in figure 1(a), the less active patch has a much higher particle density than the other. However, monodisperse disks can only be packed up to φ max = π/ √ 12 arranged in a hexagonal lattice [36] and thus a crystalline order in the less mobile region is expected for a sufficiently dense and active system. Hexagonal packing has a long-range orientational (sixfold) order, which can be measured using the local hexatic order parameter Ψ(r j ) = N j k=1 exp (i6θ jk )/N j , where θ jk is the angle formed by the bond that connects the jth disk and its kth (out of N j ) nearest neighbor (found with a Voronoi tessellation algorithm) and the x axis [27]. For a perfect triangular lattice, all the angles 6θ jk are the same and |Ψ(r j )| = 1. Because of the polycrystalline character of the dense phase we use the average magnitude of the hexatic order parameter |Ψ| = N i=1 |Ψ(r i )|/N in order to get a global information on the order [37], which   For a perfect triangular lattice |Ψ| = 1. The isoline φ p (φ 0 , P e) = 0.716 (black dashed line) indicates the transition to a pure hexatic phase [27] and is a numerical result of the theory based on (12). however does not vanish in an isotropic fluid. In figure 5 the order parameter |Ψ| within the passive patch as function of φ 0 and P e is shown. |Ψ| displays a sharp increase beyond the isoline φ p (φ 0 , P e) = 0.716 (dashed line), which indicates the transition to a pure hexatic phase [27] and is a numerical result of the theory based on (12). The minimum global packing fraction necessary for a crystalline patch is φ 0 = (1 − α)0.716.

Summary
To summarize our main findings: we have theoretically justified the pressure imbalance in an activity landscape. The pressure difference originates in the flux of polar order and the gradient of swim velocity across the interface between regions of different activity. We found that although the density is lower in the more active area the corresponding pressure is higher. For not too dense systems, φ 0 < (1 − α)φ max , the densities of the coexisting phases can be predicted from the momentum balance equation (11). Moreover, we have studied the effect of interactions on ABPs in activity landscapes. Excluded volume effect become significant above a global packing fraction φ 0 = (1 − α)φ max , i.e., when, in the limit P e → ∞, the passive patch is completely filled with particles. In that case the system is reminiscent of a system with a homogeneous swim velocity undergoing MIPS and the pressure becomes equal in both phases. The dense phase in the less-active region exhibits a crystalline order for a sufficiently dense and active system, and the transition line can also be predicted from the momentum balance equation (11).