Dynamics and density distribution of strongly confined noninteracting nonaligning self-propelled particles in a nonconvex boundary
DDynamics and density distribution of strongly confinednoninteracting nonaligning self-propelled particles in a nonconvex boundary
Yaouen Fily, Aparna Baskaran, Michael F. Hagan
Martin Fisher School of Physics, Brandeis University, Waltham, MA 02453, USA (Dated: October 17, 2018)We study the dynamics of non-aligning, non-interacting self-propelled particles confined to abox in two dimensions. In the strong confinement limit, when the persistence length of the activeparticles is much larger than the size of the box, particles stay on the boundary and align withthe local boundary normal. It is then possible to derive the steady-state density on the boundaryfor arbitrary box shapes. In non-convex boxes, the non-uniqueness of the boundary normal resultsin hysteretic dynamics and the density is non-local, i.e. it depends on the global geometry of thebox. These findings establish a general connection between the geometry of a confining box and thebehavior of an ideal active gas it confines, thus providing a powerful tool to understand and designsuch confinements.
I. INTRODUCTION
Active systems are non-equilibrium systems whoseconstituent units consume energy to generate motion ormechanical forces. Originally inspired by biology, e.g. bacterial colonies [1, 2], healing tissues [3, 4], or flock-ing animals [5], the field now encompasses a variety ofartificial systems that share this ability to inject energyat the microscopic level and emulate the unique prop-erties of their biological counterparts, from flocking tospontaneous aggregation [6–15]. Beyond biomimetism,the study of active matter has led to new applicationsnot found in nature, such as bacteria-powered micro-gears [16–18].One distinctive yet relatively unexplored property ofactive systems is their sensitivity to boundary effects.Striking macroscopic effects may be obtained by pattern-ing confining walls on the micro-scale, as exemplified bythe rectification phenomenon [19–24]. More generally,any real-world system must have boundaries, and un-derstanding their role is paramount to designing activematter based devices. Whether the boundaries are onlypresent by necessity or designed as an integral compo-nent of an active system, it is important to note thatboundary effects are not merely size effects: the exactshape of the boundary is crucial. However, most existingstudies are only concerned with one among a handful ofspecific geometries [21, 25–38], and little is known abouthow the shape of a boundary affects the behavior of theactive system it confines.In this paper, we focus on non-aligning self-propelledparticles, a model that has recently attracted attentionas a minimal model for self-propelled matter [28, 39–46].In particular, we neglect alignment interactions such asthose that would arise from hydrodynamic couplings ina fluid environment. Our results thus apply to systemsin which such coupling torques are weak. Furthermore,we restrict ourselves to the “ideal active gas” limit inwhich particles interact with the wall, but not with eachother [28]. We recently showed, for such a system, ananalytic relationship between the density and pressure of the active gas and the shape of the box for a gen-eral class of box shapes [47]. In the strong confinementregime, where the persistence length of the active parti-cles is much larger than the size of the box, and when thebox is convex, we showed that particles never leave theboundary and always align their self-propulsion directionwith the local boundary normal. Furthermore, the den-sity and the pressure on the boundary are proportionalto the local boundary curvature. It is then possible topredict the density and pressure on the boundary of anyconvex box, regardless of the details of its shape. How-ever, existing applications suggest that active devices aremost effective when their boundaries have both convexand concave regions.In this paper, we extend the theoretical framework in-troduced in Ref. [47] to the case of non-convex boxes.The presence of concave regions is a significant compli-cation, as it implies that the same normal is found atmultiple locations on the boundary, leading to multi-stability and hysteresis. Furthermore, particles withinconcave regions undergo complex, accelerated dynamicsthat sometimes launches them off the wall. Nonetheless,we demonstrate that in the strong confinement regime:(i) this complex dynamics can be understood in termsof non-local “jumps”, (ii) the density of particles withinconcave regions vanishes and (iii) it is possible to predictthe density everywhere on the boundary. We present ageneral algorithm to obtain this relationship and we testour predictions against the results of molecular dynamicssimulations in a family of boxes with both concave andconvex regions. Finally, we discuss the role of interac-tions and the limits of the ideal gas approximation.The paper is arranged as follows. Section II intro-duces the model. Section III explores the particle dy-namics on the boundary and shows how the accelerateddynamics over concave regions can be recast as instanta-neous jumps between disparate convex regions (see alsoappendices A and B). Section IV presents a theory forthe density on the boundary in the strong confinementregime, and shows how to obtain the steady-state den-sity. Section V presents the results of molecular dynamicssimulations and compares them against the predictions a r X i v : . [ c ond - m a t . s o f t ] J a n s ˆ n ˆ ν θ ψφ xy FIG. 1. (Color online) Notations for a particle at the wall.The particle is characterized by its arclength s along the walland its orientation ˆ ν = cos θ ˆ x + sin θ ˆ y . The local normal tothe wall is ˆ n = cos ψ ˆ x + sin ψ ˆ y . of sections III and IV. Section VI discusses the scope ofour model and the role of convexity in confined activegases. II. MODEL
We consider a collection of confined, overdamped,self-propelled particles in two dimensions. Each par-ticle is characterized by its position r and orientationˆ ν = cos θ ˆ x + sin θ ˆ y . The dynamics obeys the followingequations of motion:˙ r = v ˆ ν + µ F w , ˙ θ = ξ ( t ) (1)where v is the self-propulsion speed, µ is the mobility, ξ is a white Gaussian noise with zero mean and correlations (cid:104) ξ ( t ) ξ ( t (cid:48) ) (cid:105) = 2 D r δ ( t − t (cid:48) ), and over-dots indicate timederivatives. The medium on which the self-propulsionforce is exerted is treated as a momentum sink; in par-ticular, we neglect hydrodynamic interactions. The con-fining walls are hard and frictionless; whenever the ve-locity of a particle is such that it would drive the particleinto the wall, its component normal to the wall is can-celled by the wall force. More precisely, the wall force F w is zero if the particle is not at the wall, or if itis at the wall but pointing away from it; otherwise itis equal to − (cid:16) v µ ˆ ν i · ˆ n (cid:17) ˆ n , where ˆ n = cos ψ ˆ x + sin ψ ˆ y is the local normal to the wall pointing outwards [48].This is the simplest choice of wall potential consistentwith overdamped dynamics. It neglects alignment termsthat can arise when particles are anisotropic or experi-ence hydrodynamic effects, and thusis appropriate whenthe re-orientation rate induced by such torques is slowin comparison to particle translation (this point is dis-cussed further in section VI). Finally, since particles arenon-interacting, we may restrict the discussion to a singleparticle.When a particle is at the wall, its configuration is char-acterized by its arclength s ∈ [0 , L ) along the boundary,where L is the box perimeter, and its orientation relativeto the local boundary normal φ = θ − ψ (see Fig. 1).There are two important lengths scales in the system:the active persistence length v /D r , i.e. the typical dis-tance a free (unconfined) particle travels before its orien-tation decorrelates, and the global size of the confining box. However, for a boundary with nonuniform curva-ture, the variations in the local radius of curvature leadto additional length scales. Which one is most relevantdepends strongly on the geometry of the box, as discussedin Ref. [47] and in the rest of the paper.The regime we study in this paper is the strong confine-ment regime , obtained when the persistence length v /D r is much larger than the size of the box. It is obtained atlarge self-propulsion, small angular noise or small boxsize, and we will often refer to it as the small angularnoise regime , or simply the small noise regime (the an-gular noise is the only noise in our model). III. DYNAMICS AT THE WALL
We first look at the dynamics of a single particle mov-ing along the wall. In particular, we explore the fun-damental difference between convex and concave regionsand show how the latter cause fast jumps and bi-stabilityin the low noise regime. To this end, we project the equa-tions of motion (1) onto the tangent to the wall:˙ s = v sin φ , ˙ φ = ξ ( t ) − v R ( s ) sin φ (2)where s is the arclength along the wall, φ = θ − ψ isthe angle between the particle’s orientation ˆ ν and theboundary normal ˆ n (see Fig. 1), and R ( s ) = ds/dψ isthe local radius of curvature. Eqs. (2) remain valid aslong as | φ | ≤ π/
2; as soon as | φ | > π/
2, the particleleaves the boundary.
A. Dynamics at zero angular noise
In the absence of noise ( D r = 0), the orientation ˆ ν =cos θ ˆ x + sin θ ˆ y of the particle is constant and its glidingvelocity only depends on its location (see Fig. 2):˙ s = v sin( θ − ψ ( s )) (3)When | θ − ψ | > π/ ψ ( s ) = θ ( i.e. theparticle is aligned with the normal) act as fixed points.They are stable in convex regions ( R ( s ) >
0) and un-stable in concave regions ( R ( s ) < s of the fixed point: ddt ( s − s ) = − v R ( s ) ( s − s ) + O (cid:0) ( s − s ) (cid:1) (4)Graphically, the fixed points are located at the in-tersection(s) of the curve y = ψ ( s ) with the horizontaldashed line y = θ in the top right panel of Fig. 2. In con-vex boxes, ψ ( s ) is monotonic and the fixed point corre-sponding to the orientation θ is always stable and unique. ABCD s s θπψ AB CD ˆ ν − v v Gliding velocity
FIG. 2. (Color online) Dynamics in a confining box in theabsence of angular noise. Top left: Real space representation.The concave region (between B and C ) is shown in orange(grey). The arclength s is counted counter-clockwise. Thepoints A and D have the same normal (shown as a straightarrow) as C and B , respectively. Top right: Normal angle ψ as a function of arclength s . The dynamics of a parti-cle with constant orientation θ is controlled by the distance φ = θ − ψ between the curve and the horizontal dashed line.The fixed points as shown as filled (stable) and empty (un-stable) squares. Bottom left: Gliding velocity ˙ s = v sin φ fora particle with orientation ˆ ν (shown inside the box). The ori-entation and fixed points are the same as those shown in thetop right panel. Arrow heads indicate the gliding direction.The colorbar is shown on the right. Conversely, in the presence of concavity there are multi-ple locations with the same normal, and multiple fixedpoints for some values of θ . B. Quasi-static dynamics
We now let the orientation θ vary slowly. When therate of change of θ is small enough, the particle spendsmost of its time at a fixed point, and only a small fractionof its time travelling between fixed points. This quasi-static regime is obtained when the angular noise D r issmall, i.e. in the strong confinement regime. The particleis then confined to the boundary: the fixed point condi-tion θ = ψ implies that the particle always points towardthe boundary, whereas leaving the boundary would re-quire pointing away from it ( | θ − ψ | > π/ dθ in the orientation θ causesa small displacement ds = R ( s ) dθ of the correspondingfixed point. In a convex region, the particle relaxes ex-ponentially toward the new fixed point. The quasi-staticregime is then obtained when the corresponding relax-ation time R ( s ) /v is much shorter than the reorientationtime D − .In a concave region, on the other hand, an infinitesimalchange in the orientation θ can trigger a large displace-ment, which we refer to as a jump . Consider the boxshown in Fig. 2, and a particle at point B with orientation θ = ψ B + dθ where dθ > ψ B + dθ in thevicinity of B , the particle has to travel to the next convexlocation with normal angle ψ B , i.e. point D . Concretely,the perturbation dθ sends the particle into the concaveregion where its gliding speed continuously increases. Itonly starts to decelerate once it reaches the end C of theconcave region, and eventually comes to a stop at point D . The quasi-static regime is obtained when the jumpfrom B to D is much faster than the reorientation time D − . In appendix B we show that this is the case forsmall angular noise. We also discuss the possibility andimplications of particles leaving the boundary during ajump. Within the quasi-static approximation, however,the details of a jump are irrelevant: it is considered in-stantaneous, and only its landing point matters.In summary, the presence of concavity causes non-trivial dynamics in the quasi-static regime. A particlereaching the end of a convex region experiences an in-stantaneous jump to a new convex location with the samenormal angle. As a result, the vicinity of a concave re-gion exhibits bi-stability and hysteresis (in Fig. 2, jumpsfrom B to D and from C to A create an hysteresis looparound ABCD ).Finally, these results have important consequences forthe steady-state density (see section IV). First, the in-stantaneousness of the jumps over concave regions im-plies that those regions are empty. Second, the fact thatjumps do not stop at the end of the concave region butcontinue into the next convex region causes non-localdensity fluxes within the system. As we show next, therequirement that these fluxes cancel at steady-state en-ables predicting the density profile everywhere on theboundary.
IV. QUASI-STATIC STEADY-STATE DENSITY
We now use the results of section III to predict thesteady-state density of a particle in a box of arbitraryshape in the quasi-static regime. Our starting point isthe assumption that the particle is always at a stablefixed point, i.e. at a convex point where the particle’sorientation is aligned with the boundary normal ( θ = ψ ).This has three important consequences. First, the parti-cle always points toward the boundary and never leavesit. As a result, the density is zero in the bulk and theproblem is effectively limited to the (one-dimensional)boundary. Second, the particle never visits concave re-gions, where fixed points are always unstable. The den-sity on the boundary therefore vanishes in those regions.Third, the normal angle at the location of the particlefollows a simple random walk: ˙ ψ = ˙ θ = ξ ( t ) where ξ isthe same noise as in Eqs. (1) and (2). Thus, the densityon the boundary in ψ space (unit circle) obeys the usualdiffusion equation: ∂ t ρ ψ = D r ∂ ψ ρ ψ . (5)whose steady-state solution is given by ρ ψ ( ψ ) = 12 π . A. Role of convexity
In a convex box, R ( s ) is positive everywhere and ψ ( s )is monotonic. The density of particles per unit lengthof boundary is then obtained by making the change ofvariable ψ → s . At steady-state, this yields [47] ρ ( s ) = ρ ψ ( ψ ) dψds = 12 πR ( s ) (6)The density on the boundary is thus proportional to thelocal boundary curvature.In a non-convex box, on the other hand, there are mul-tiple locations on the boundary with the same normal an-gle ψ (see Fig. 2). Unlike the density ρ ( s ), the normal an-gle density ρ ψ ( ψ ) does not discriminate between those lo-cations. Thus, ρ ( s ) cannot be inferred from ρ ψ ( ψ ) alonefor a box with concave boundary regions. B. Formulating the problem
To retain all the information contained in ρ ( s ) whileworking in normal angle space, where the dynamics issimply diffusive, we number the convex regions 1 to n andintroduce the normal angle density ρ ψi in region i . Region i is delimited by the two inflexion points A i − and A i (see Fig. 3). Inflexion point A i has normal angle ψ i andarclength s i . Each ρ ψi is defined over the entire [0 , π )interval but only takes non zero values between ψ i − and ψ i , so that we can write ρ ψ = (cid:80) i ρ ψi . The concaveregions, where the density is assumed to be zero, are notexplicitly described, but manifest themselves through theboundary conditions ρ ψi ( ψ i − ) = ρ ψi ( ψ i ) = 0. Inflexionpoints act as one-way teleportation devices that send theparticle to a new convex location with the same normalangle ψ . These instantaneous jumps often, but not al-ways, connect consecutive regions (see appendix A).From the point of view of each ρ ψi , the start of a jumpis a particle sink and its end a particle source. Apartfrom jumps, the dynamics in a convex region is indistin-guishable from that in a convex box and ρ ψi ( ψ ) obeys thesame diffusion equation as ρ ψ ( ψ ). The result is a set ofcoupled diffusion equations: ∂ t ρ ψi = D r ∂ ψ ρ ψi + (cid:88) k (cid:15) ik J k δ ( ψ − ψ k ) (7)where δ is the Dirac delta function and the sum is overjumps. Jump k occurs at normal angle ψ k and carries acurrent J k > A k to its landingpoint B k . (cid:15) ik encodes the relationship between jump k and region i . If jump k starts in region i , (cid:15) ik = − k lands in region i , (cid:15) ik = 1 and the corresponding A B A B A B A B A B A B s s s s s s s s s s s s L s0 π πψ A B A B A B A B A B A B πψ A B A A B B A A B B A B FIG. 3. (Color online) Three representations of a non-convexbox showing the concave regions in orange (grey) and thejumps over the concave regions (curved arrows). Convex re-gions are indexed by a number 1 to 3. Region i is delimitedby the two inflexion points A i − and A i . Each A i is thestarting point of a jump over the neighboring concave regionthat lands at B i , which has the same normal angle ψ i as A i .Top: Shape of the box. The ‘x’ on the right is the arclengthorigin, s = 0. Bottom left: Normal angle vs. arclength. Bot-tom right: Normal angle representation. The vertical axislabels the convex region corresponding to each interval. Re-gion 1 appears split into two parts due to periodic boundaryconditions. term in Eq. (7) is a source. Finally, if jump k does notinvolve region i then (cid:15) ik = 0; in other words the sum inEq. (7) is restricted to jumps that involve region i [49].In order for ∂ t ρ ψi to remain finite, each jump must createa discontinuity in ∂ ψ ρ ψi proportional to J k : ∂ ψ ρ ψi ( ψ + k ) − ∂ ψ ρ ψi ( ψ − k ) = − J k /D r where the superscripts ± symbolizeone-sided limits. The currents J k may then be eliminatedfrom Eq. (7) in favor of the densities ρ ψi .Once the density in normal angle space is known inevery convex region, the linear density on the boundary ρ ( s ) is obtained using the change of variable ψ → s ,which is monotonic within each convex region: ρ ( s ) = ρ ψi ( ψ ( s )) R ( s ) if s is in convex region i s is in a concave region (8)where R ( s ) is the radius of curvature. C. Steady-state
From the form of Eq. (7), it is clear that the steady-state density in each convex interval is piecewise linear,with a change of slope at the location ψ k of every jumpthat starts or lands in the interval. We also require ρ ψi tobe continuous, and to vanish at the ends of the interval(beyond which the boundary is concave and thus empty): ρ ψi ( ψ i − ) = ρ ψi ( ψ i ) = 0. As a result, the entire set ofdensity functions { ρ ψi ( ψ ) } ≤ i ≤ n is fully determined byits 2 n values at the location of every jump landing B k (see Fig. 3). Let x k = ρ ψi ( k ) ( ψ k ) , k ∈ [1 , n ] be thoseunknowns, with i ( k ) the index of the convex region inwhich jump k lands.On the other hand, Eq. (5) implies that the steady-state total density ρ ψ = (cid:80) i ρ ψi is equal to 1 / (2 π ).Using the piecewise linearity of ρ ψi , ρ ψi ( ψ ) can alwaysbe expressed as a linear combination of x k ’s. Writing ρ ψ ( ψ ) = (2 π ) − at 2 n distinct values of ψ then leads toa 2 n × n linear system that can be solved to obtain the x k ’s. In order for the system not to be degenerate, the2 n values of ψ need to be spread across all subregions ofall convex regions so that no x k is left out. A convenientchoice is to use the locations ψ k of the jumps.Once the ρ ψi ’s are known, the density per unit lengthof the boundary is given by Eq. (8).The entire process can be automated, i.e. it is possi-ble to write a program that, starting with a parametriza-tion of the boundary, identifies the convex regions andthe jumps, writes the linear system, solves for the x k ’s,and generates the functions ρ ψi and ρ . We used such aprogram to create most of our figures and analyze oursimulation results (section V).Below we go through the method step-by-step in sev-eral situations of interest. In sections IV C 1 and IV C 2,we consider boxes in which the number of convex loca-tions with the same normal is limited to 2 (section IV C 1)and 3 (section IV C 2). In both cases, we show thatthe density can be expressed in a simple form. In sec-tion IV C 3 we give a detailed description of the algorithmthat leads to the steady-state density in the general case.
1. Multiplicity no greater than We start with the box pictured in Fig. 3. In order towrite ρ ψ ( ψ k ) = (2 π ) − in terms of the unknowns { x j } ,we look at the normal angle representation (bottom rightpanel of Fig. 3) and take a vertical slice at ψ k . Lookingat, e.g. , the slice through B , we see that A and B arethe only two locations on the boundary where the normalangle is equal to ψ . At B , the density is by definition x , while at A it is zero because it is the entrance ofa concave region. Therefore, the total density at ψ issimply x , and the equation for that slice is x = (2 π ) − .The situation is the same at every jump location, and ψ ψ ψ ψ ψ ψ π πρ ψ πρ ψ πρ ψ FIG. 4. (Color online) Density of particles in normal anglespace as a function of the normal angle ψ in each of the threeconvex lobes of the box shown in Fig. 3. ψ i is the normal angleof the inflexion point A i as well as that of its correspondinglanding point B i (see Fig. 3). Ls s s s s s s s s s s s s00.4 ρ ρ FIG. 5. (Color online) Density of particles per unit lengthof boundary as a function of arclength s for the box shownin Fig. 3. s i and s (cid:48) i are the arclengths of the inflexion points A i and corresponding landing points B i (see Fig. 3)). Rightpanel: Heat map for the density along the boundary, shownin real space. the corresponding linear system is trivial: ∀ k ∈ [1 , n ] , x k = 12 π (9)Recalling the piecewise linearity and boundary conditions ρ ψi ( ψ i − ) = ρ ψi ( ψ i ) = 0, the resulting partial densities ρ ψi in ψ space then take the form shown in Fig. 4.Finally, Fig. 5 shows the density in real space and s space, where the overlaps between convex regions disap-pear and the empty concave regions reappear.More generally, let m ( ψ ) be the multiplicity, i.e. thenumber of convex locations that have normal angle ψ .Graphically, in the normal angle representation (bottomright panel of Fig. 3), m ( ψ ) is the number of convexintervals that intersect the vertical line at ψ . Since ajump connects two regions, m ( ψ k ) ≥
2. If m ( ψ ) neverexceeds 2 over the entire boundary, then A k and B k arethe only convex locations with normal angle ψ k , and theequation ρ ψ ( ψ k ) = (2 π ) − always takes the trivial formgiven by Eq. (9), regardless of the number n of convexregions or their sizes. A B A B A B A B A B A B A B A B
123 40 Ls0 π πψ A A B A B A B A B A A B A B πψ B A B A B A A B B A B A B B FIG. 6. (Color online) Three representations of a non-convexbox showing the concave regions in orange (grey) and thejumps over the concave regions (curved arrows). Convex re-gions are indexed by a number 1 to 4. Top: Shape of thebox. The origin of arclengths is the rightmost point on theboundary, halfway between B and B . Bottom left: Normalangle vs. arclength. Bottom right: Normal angle representa-tion. The vertical axis labels the convex region correspondingto each interval. Region 1 appears split into two parts due toperiodic boundary conditions.
2. An example with multiplicity We now consider the box shown in Fig. 6, defined inpolar coordinates by r ( θ ) = 1 + 0 . θ ).There are n = 4 convex regions, delimited by the 8 in-flexion points A to A . The jump starting at A i ends at B i , which is located in a neighboring lobe. Although thejumps involve leaving the boundary and flying straightthrough the bulk for a short period of time (see sec-tion A), these flights do not alter the location of the land-ing points B i and thus are irrelevant to the steady-state.In real space, the most important difference betweenthis box and the box of Fig. 3 is the order of the landingpoints. For example, B comes before B when movingcounter-clockwise. This “inversion” is a signature of mul-tiplicities higher than 2, as can be seen by comparing thebottom panels of Figs. 3 and 6.Following the method outlined in previous sections, weset out to write ρ ψ ( ψ k ) = (cid:80) i ρ ψi ( ψ k ) = (2 π ) − at eachjump location ψ k in terms of the variables x k . We startwith jump 7 and draw a virtual vertical line through A and B in the normal angle representation (bottom right ψ ψ ψ ψ ψ ψ ψ ψ π πρ ψ πρ ψ πρ ψ πρ ψ FIG. 7. (Color online) Density of particles in normal anglespace as a function of the normal angle ψ in each of the fourconvex lobes of the box shown in Fig. 6. panel of Fig. 6). This line intersects region 4 at A , region3 at B , and region 2 at a point C located between B and A . We can therefore write(2 π ) − = ρ ψ ( ψ ) + ρ ψ ( ψ ) + ρ ψ ( ψ )= x + ψ − ψ ψ − ψ ρ ψ ( ψ ) + ψ − ψ ψ − ψ ρ ψ ( ψ )= x + ψ − ψ ψ − ψ x (10)The second line is obtained by using the boundary con-dition ρ ψ ( ψ ) = 0, the definition ρ ψ ( ψ ) = x , and thelinearity of ρ ψ between B and A . The third line is ob-tained by using the boundary condition ρ ψ ( ψ ) = 0 andthe definition ρ ψ ( ψ ) = x . Applying the method to eachjump yields the system2 π α α α α
00 0 0 0 1 0 0 α α α α · x x x x x x x x = (11)where α k = ψ k − − ψ k ψ k − − ψ k − (indices are defined modulo2 n ). The solution to this system is not in general com-pact; however, here the four-fold symmetry implies that α k ≡ α is independent of k ,and the solution simply reads: ∀ k, x k = 12 π (1 + α ) (12)The corresponding densities in normal angle space, inarclength space and in real space are shown in Figs. 7and 8. s s s s s s s s s s s s s s s s L s0 . . ρ ρ FIG. 8. (Color online) Density of particles per unit lengthof boundary as a function of arclength s for the box shownin Fig. 6. Right panel: Heat map for the density along theboundary, shown in real space.
3. General case
The method used above to compute the steady-statedensity in simple boxes can be readily extended to arbi-trary shapes. We now write down the steps leading tothe solution in the general case.As a preliminary step, in each region i we iden-tify the set K i = { k ∈ [1 , n ] | s (cid:48) k ∈ [ s i − , s i ] } of jumps that land in the region, with s (cid:48) k the ar-clength of landing point B k . We also define the set E i = { A i − , A i } ∪ { B k | k ∈ K i } of jump ends inregion i , i.e. the two ends of the region plus any jumplandings. Then, for each jump landing B k , we performthe following sequence of operations:1. Identify the region i that contains B k . By defini-tion ρ ψi ( ψ k ) = x k .2. Pick a region j other than i . In the set E j , identifythe two points C and D whose normal angles areclosest to ψ k on each side: ψ C = max M ∈E j { ψ M | ψ M ≤ ψ k } ψ D = min M ∈E j { ψ M | ψ M ≥ ψ k } By construction, ρ ψj is linear between C and D ;therefore ρ ψj ( ψ k ) = ψ D − ψ k ψ D − ψ C ρ ψj ( ψ C ) + ψ k − ψ C ψ D − ψ C ρ ψj ( ψ D ) . (13)Furthermore, ρ ψj ( ψ C/D ) is either 0 if
C/D is A j − or A j , or x l if C/D is B l ; therefore the right-handside of Eq. (13) is a linear combination of 0, 1 or 2of the x l ’s (if the region has no jump landing, or if ψ k is outside of the region, the right-hand side issimply zero).3. Repeat the previous step until every region hasbeen considered, then sum Eq. (13) over j . The left-hand side is (by definition) ρ ψ ( ψ k ) = (2 π ) − ,while the right-hand side is a linear combination of x l ’s.At the end of step 3, a linear system of the form (cid:80) j a ij x j = (2 π ) − is obtained, whose coefficients a ij de-pend on the ψ k ’s. After inverting the system to get the x j ’s and thus the normal angle space densities { ρ ψi ( ψ ) } ,“unfolding” normal angle space onto arclength space anddividing by the radius of curvature yields the density ρ ( s )per unit length of boundary (see Eq. (8)).As pointed out at the beginning of section IV, theprocess can be automated. This is particularly usefulwhen the number of non-zero elements in a ij is large,which happens when the multiplicity m is large. On theother hand, a small multiplicity leads to a sparse ma-trix. In particular, m = 2 yields a diagonal matrix (seesection IV C 1). V. SIMULATIONS
We test the results of sections III and IV by performingmolecular dynamics simulations of Eqs. (1) in the familyof boxes defined in polar coordinates by r ( θ ) = 1 + r sin ( kθ ) (14)The boxes shown in Figs. 3 and 6 belong to this family,with ( k, r ) = (2 , .
5) and (4 , .
6) respectively.
A. Angular distribution
We first test our core assumption, that the deviation ofthe particle orientation from the boundary normal van-ishes: φ ≈
0, and the reasoning that led to it (see sec-tion III), by characterizing the statistics of φ at variouslocations along the boundary. In convex regions, the the-ory predicts that the distribution of φ will exhibit a nar-row peak centered around φ = 0. On the other hand,particles jumping over a concave region should generatesecondary peaks, much smaller than the primary peak,and centered around a position-dependent non-zero value φ = ψ − ψ where ψ and ψ are the normal angles at theinflexion point where the particle entered the jump andat the current location, respectively. In concave regions,there should be no central peak, only “secondary” ones.To compare these predictions with the simulation re-sults, we measure the distribution P ( φ ) at regularlyspaced locations along the boundary. At each point, wedetermine the heights of the distribution’s peaks { P max } and their corresponding orientations { φ max } . In Fig. 9we show the peak orientations as a function of arclengthin a single convex region, as well as the prediction asso-ciated with each jump involving the region. Every de-tected peak falls on one of the predicted branches: φ = 0or φ = ψ i − ψ with i ∈ , , , A , A , A , or A . − L/ L/ s − π/ π/ φ max s = + L s = − L A B A A A B FIG. 9. (Color online) Positions of peaks in the distributionof the orientations P ( φ ) relative to the normal as a functionof the arclength s over one lobe of the box defined by Eq. (14)with k = 4 and r = 0 .
5. The lobe geometry is shown in theright panel. The angular diffusion constant is D r = 10 − ,deep in the strong confinement regime. The simulation datais shown as circles. The theoretical prediction is shown aslines. Colors denote the origin of the particles: blue for parti-cles jumping from A (upper gray curve), orange for particlesjumping from A (lower gray curve), and black for particlesfrom the lobe under consideration. Dotted lines represent thepaths of particles that fly through the bulk. Fig. 9 also illustrates particles leaving the boundary tofly straight through the bulk (see middle row of Fig. 12):since φ is only defined at the boundary, the peak corre-sponding to these particles is absent between their leav-ing the boundary and their hitting it again (dotted linesin Fig. 9). B. Steady-state density
We now assess the accuracy of the predictions made insection IV by plotting in Fig. 10 the observed boundarydensity ρ as a function of the predicted density ρ pr forvarious box shapes and angular diffusion constants D r .The plot can be interpreted in terms of two linearasymptotes. In boundary locations with moderate andhigh density, we observe good agreement between theoryand simulation, i.e. ρ ≈ ρ pr , up to D r ∼ − . Deep inthe strong confinement regime ( D r = 10 − ), the agree-ment is excellent and persists over two decades. In loca-tions with low density, on the other hand, the predictionunderestimates the density: the predicted density van-ishes altogether over regions of zero or negative curvaturewhile the observed density is finite everywhere, resultingin a horizontal asymptote whose position depends on thebox’s shape and the angular noise’s strength. Note thatbecause of the logarithmic scale, regions where the pre-dicted density is zero do not appear in Fig. 10; insteadthe data visible on the left of the plot comes from weaklyconvex areas near the inflexion points.As discussed in section III B and appendix B, the den-sity in flat and concave regions is controlled by the ra-tio of the time spent crossing them to the reorientation − − − Predicted Density10 − − − D e n s i t y D r − − − − FIG. 10. (Color online) The density observed in simulationsat various positions along the boundary is plotted againstthe predicted density for six boxes from the family definedby Eq. (14) and for four values of the angular diffusion con-stant D r . The dashed line corresponds to a density equal tothe predicted density. The symbol associated with each boxgeometry is shown at the center of that box on the right side. time D − . This ratio only vanishes in the D r → D r . Since the crossing time is largest inflat regions, this is where the deviations from the quasi-static theory are most prominent. Additionally, the timeit takes to cross the vicinity of an inflexion point growswith its “flatness”, as inferred from the second derivativeof the normal angle with respect to the arclength ( i.e. the derivative of the curvature). This explains the po-tentially counter-intuitive observation that the accuracyof our predicted density is better in “strongly concave”boxes [ e.g. the box denoted by left-pointing triangles ( (cid:47) )in Fig. 10] than in “weakly concave” boxes [ e.g. the boxdenoted by squares ( (cid:3) )]. Despite their weaker concavity,the latter exhibit larger “flat” regions. ψ ψ ψ ψ ψ πρ ψ D r − − − − FIG. 11. (Color online) Density in normal angle space inthe convex region shown in Fig. 9 for several values of theangular diffusion constant D r from deep in the strong confine-ment regime ( D r = 10 − ) to outside of it ( D r = 10 − ). Thedashed line is the theoretical prediction from section IV C 2(see Eq. (12)), corresponding to D r → Finally, we plot in Fig. 11 the density ρ ψ = Rρ in nor-mal angle ( ψ ) space. There, each convex region must betreated separately; we consider the region labelled 1 inFig. 6. The predicted density is piecewise linear with atrapezoidal shape. At very small noise ( D r = 10 − ), thedensity observed in simulations closely matches the pre-diction. As the noise is increased, the trapezoidal shapegets smoothed out and some of the density is transferredfrom the tip to the base (as well as the neighboring con-cave region, not shown in this representation), suggest-ing that finite noise effects may be treated as a pertur-bation to our zero-noise theory. On the other hand, at D r = 10 − the predicted form of the density is not rec-ognizable anymore and a different approach is required.Note that this description over-emphasizes the finitenessof the observed density at the ends of the convex intervalwhere the radius of curvature, and thus ρ ψ , diverges. VI. DISCUSSION
In summary, we have presented a systematic approachto predict the density of a non-aligning ideal active gasin a small box of arbitrary shape, thus establishing aconnection between the geometry of a confining box andthe properties of the active gas it confines. Our resultshold as long as the persistence length (the distance a freeparticle travels before it loses its orientation) is muchlarger than the size of the box.In the special case of convex boxes, there is a strikinglysimple relationship between the density and the bound-ary geometry: the density is zero in the bulk and on theboundary it is proportional to the local boundary cur-vature [47]. Here, we have shown that boundaries withconcave regions lead to a much richer particle dynamics,including multi-stability, hysteretic dynamics, and parti-cles flying through the bulk of the box between disparateboundary locations. However, we showed that the par-ticle density still vanishes in the bulk and we describedan algorithm to calculate the steady-state density profileon the boundary of a 2D box with any shape. The pre-dicted particle density vanishes in concave regions, whilein convex regions it can be written as the product of thelocal curvature and a “splitting factor” which obeys thefollowing property: given a unit vector ˆ n , the sum of thesplitting factor over all the locations on the boundarywhere the normal is ˆ n is equal to one. In other words,boundary points that share the same normal also sharethe same “pool” of particles with the corresponding ori-entation.Despite the complexity intrinsic to concave regions,understanding non-convex shapes is essential to ratio-nally design active micro-devices with specific function-alities. This is nicely illustrated by the micro-gearsused in Refs. [16–18], which trap particles in sharp cor-ners where they exert torques that make the gear ro-tate [50]. Effectively trapping active particles requiressharp corners [29, 30, 33, 47], and the total torque on thegear is maximized by having several such trapping sites.This cannot, however, be achieved with convex shapes in which the number of sharp corners is limited to two.It is then clear that understanding and using non-convexconfinements is a necessary step toward designing a broadclass of active devices. Scope of the model.
Finally, we consider limitationsand avenues for extension of the model. Firstly, since weneglect interparticle interactions our results are limitedto dilute systems. For example, above a threshold pack-ing fraction, steric effects will prevent all particles fromresiding on the boundary. The effect of steric interactionswill be discussed in a future publication, but preliminarysimulations confirm that our results apply at least quali-tatively at finite particle densities. Secondly, our resultsapply to the strong confinement limit, in which particlescircumnavigate the box faster than they reorient and thustend to align with the boundary normal. Within thislimit, we expect the results to remain valid regardless ofthe reorientation mechanism, including angular diffusionor aligning interactions with the wall such as may arisedue to hydrodynamics. However, if particle-wall interac-tions drive particles to align with the wall, or divert awayfrom it, on timescales comparable to the circumnaviga-tion time ( ∼ v /R with R the boxsize), then a differentapproach is required. Appendix A: Types of jumps
Depending on the geometry of the concave region theycross, jumps may lead to particles leaving the boundaryto travel in the interior of the box. This situation, whichis illustrated in Fig. 12, occurs when the angle betweenthe normals at the two ends of the concave region ( B and C ) is larger than π/
2. As a result, there exists a point E between B and C where φ = ψ B − ψ E = π/
2. At E ,the particle’s orientation is aligned with the tangent andit leaves the boundary to travel in a near-straight linethrough the box, until it hits the boundary again. Thethree possible jump scenarios are explained in Fig. 12. Appendix B: Jump duration
To evaluate the duration of a jump, we consider thebox shown in the top panel of Fig. 12 and a particle atthe entrance B of the concave region with orientation θ = ψ B + dθ where dθ > D r = 0), the gliding speed ˙ s = v sin φ iscontrolled by the relative angle φ = θ − ψ ≈ ψ B − ψ ,shown in the top right panel of Fig. 12.In particular, we want to show that in the small angu-lar noise limit D r →
0, the duration of the jump from B to D is negligible compared with the reorientation time τ ≡ D − .To this end, we divide the jump into three regions num-bered (i), (ii) and (iii) (see top right panel of Fig. 12),0 BCD EFD BCD EF BCD s s π/ ψ B − ψ B C D
E FD s π/ ψ B − ψ B C D
E Fs π/ ψ B − ψ B C D(i) (ii) (iii)
FIG. 12. (Color online) Types of jumps over a concaveregion. The right column shows the orientation φ = ψ B − ψ relative to the boundary as a function of the arclength s (counted counter-clockwise). Concave regions are in orange(grey). Trajectories through the interior of the box are shownas arrows. B and C are the inflexion points delimiting theconcave region. D is the first location after C with the samenormal as B . Top row: the particle follows the boundaryand stops at D . See section B for the definition of regions(i), (ii), (iii). The dashed (in green) and the dot-dashed (inblue) curves correspond to a quadratic expansion near B anda linear expansion near D , respectively. Middle row: when ψ B − ψ C > π/
2, the particle leaves the boundary at E where ψ E = ψ B − π/
2. It travels in a straight line through the bulkuntil F , then follows the boundary to D . Bottom row: when ψ B − ψ C > π/ D to end its jump in a non-neighboringconvex region at a point D (cid:48) with the same normal as B . and define τ α as the time it takes to cross region α . Re-gions (i) and (iii) correspond to the vicinity of B and D , respectively, while region (ii) contains the remainingmiddle part of the jump.Region (ii) is the fastest of the three, with a relativeangle φ ∼ s ∼ v . This leads toa crossing time τ ii ∼ (cid:96)/v where (cid:96) is the length of thejump, which is of the same order or smaller as the sizeof the box. Note that instances of the particle leavingthe boundary (see section A, and the middle and bottomrow of Fig. 12) do not modify the scaling of τ ii . Indeed,when flying through the box the particle travels at speed v for a distance at most of the order of the box size.Comparing with the reorientation time, we get τ ii /τ ∼ (cid:96)D r /v , which vanishes in the small noise limit. In region (iii), linearizing the gliding velocity around D leads to Eq. (4) with s = s D the arclength of point D . The dynamics is an exponential relaxation towards D , which takes a time τ iii ∼ R D /v where R D is theradius of curvature at point D [51]. Comparing with thereorientation time, we get τ iii /τ ∼ R D D r /v , which alsovanishes in the small noise limit.Region (i) is the slowest part of the jump. While B isa fixed point, it is also an inflexion point, and thus thelinear contribution in Eq. (4) vanishes. A second orderexpansion yields ˙ x = v dθ + v κx where x = s − s B , κ = ( d ψ/ds ), and v dθ is the initial velocity requiredto start the jump. Integrating leads to the displacement x ( t ) = (cid:114) dθκ tan (cid:16) v t √ κdθ (cid:17) (B1)However, evaluating τ i requires additional assumptionson dθ . In fact, noise and curvature both play a role insetting τ i . Although noise is initially the only contri-bution, its fluctuations soon get amplified by curvatureeffects, and understanding their interplay is necessary toevaluate τ i .For our purpose, however, it is sufficient to note thatthe curvature in region (i), however small, always makesthe particle progress faster than on a flat edge. The caseof a flat edge was discussed in Ref. [47] in the contextof polygonal boxes, and we briefly summarize it here. Inthe absence of curvature, the relative angle φ follows arandom walk. Its typical value then grows diffusively: φ ∼ ( D r t ) / . Integrating with respect to time givesthe typical distance travelled along the boundary after atime t : s ∼ v D / t / . Travelling a length (cid:96) then takesa typical time τ ∼ ( (cid:96)/v ) / D − / . The length of region(i) can be evaluated by noting it ends when ψ − ψ B ≈ κx is of order 1. This gives (cid:96) ≈ κ − / and τ i ∼ ( v κD r ) − / .Comparing with the reorientation time, we get τ i /τ ∼ ( v D r ) / ( κ ) − / , which goes to zero in the small noiselimit as well.In summary, in the small noise limit D r → d ψ/ds ) of the normal angle with respect toarclength. In other words, the flatter the region nearinflexion points, the slower the associated jumps. ACKNOWLEDGMENTS