Collective motion of active particles exhibiting non-reciprocal orientational interactions

We present a Brownian dynamics study of a 2D bath of active particles interacting among each other through usual steric interactions and, additionally, via non-reciprocal avoidant orientational interactions. We motivate them by the fact that the two flagella of the alga Chlamydomonas interact sterically with nearby surfaces such that a torque acts on the alga. As expected, in most cases such interactions disrupt the motility-induced particle clustering in active baths. Surprisingly, however, we find that the active particles can self-organize into collectively moving flocks if the range of non-reciprocal interactions is close to that of steric interactions. We observe that the flocking motion can manifest itself through a variety of structural forms, spanning from single dense bands to multiple moderately-dense stripes, which are highly dynamic. The flocking order parameter is found to be only weakly dependent on the underlying flock structure. Together with the variance of the local-density distribution, one can clearly group the flocking motion into the two separate band and dynamic-stripes states.

The interdisciplinary field of active matter [1][2][3][4][5] , situated at the crossroad of physics, chemistry, biology and engineering, examines individual and collective behavior of units capable of autonomous motion. The units can be diverse, including, for instance, bacteria that use body appendages to swim in viscous media, active colloids that exhibit propulsive motion by generating flow around their surfaces by means of self-phoresis, bird flocks and sheep herds. Many interesting properties of these systems can be captured already with simple models which do not account explicitly for the medium surrounding active units. Instead, on the coarse-grained level the medium serves as a source of fluctuations and drag forces acting on the units. Typically, active units are modeled as particles moving with a roughly constant speed along an orientation vector attached to them 6 , which can change its direction due to rotational diffusion and inter-particle interactions. On one side of the spectrum, there are models with pairwise interactions. Possibly the simplest example within this class of models involves isotropic active particles interacting with each other via steric interactions, which can lead to the phenomenon of motilityinduced phase separation 7-10 (MIPS) between dense particle clusters and dilute regions of freely moving particles. In a nutshell, MIPS takes place when active particles are sufficiently fast and numerous 11,12 , increasing chances that free particles join existing clusters compared to chances that particles within clusters escape them as a result of rotational diffusion. On the other side, there are models with social interactions 13,14 . The most notable example within this class is the Vicsek model 15 , incorporating alignment with average particle orientation of all neighbors situated within some fixed radial distance, which can result in coherent flocking motion [16][17][18][19] : particles spontaneously form traveling bands propagating along a common direction. Models fusing both type of interactions 20 are also found in literature. Furthermore, different types of orientational interactions have been introduced [21][22][23][24][25][26][27][28] . In this article we introduce a specific form of orientational interactions, which are non-reciprocal and avoidant. We demonstrate how they influence the collective dynamics of our model active particles. The orientational interactions are motivated by the fact that the two flagella of the alga Chlamydomonas interact sterically with nearby surfaces such that the alga experiences a torque, which rotates it away from the surface.
If pairwise interactions stem from a potential, they are reciprocal due to Newton's action-reaction principle. However, the reciprocity can be broken if the pairwise interaction originates in some non-equilibrium mechanism 29 . Numerous examples exist in literature, including carpets of microfluidic rotors 30 , droplets exhibiting predator-prey interactions driven by non-reciprocal oil exchange 31 , neural networks 32,33 and many others 34 . An eminent example of non-reciprocal interactions from the realm of active matter is found in binary mixtures of active colloids with phoretic interactions [35][36][37][38] . In essence, colloids with two distinct hemispheres, immersed in a solution containing an appropriate chemical solute, can produce and maintain a local chemical concentration www.nature.com/scientificreports/ gradient by means of, for instance, self-diffusiophoresis or self-thermophoresis. The chemical gradient implies a hydrodynamic flow of the surrounding solvent, which in turn propels the colloid forward so that the entire system remains force free. Two populations of different active colloids will exhibit a different response to self-generated chemical gradients, leading in general to non-reciprocal interactions, the strength of which can be controlled by varying colloidal mobilities and activities. Interestingly, it has been reported that interplay of hydrodynamics and boundaries can also lead to non-reciprocal interactions 30 .
In this article we study the collective motion of a 2D low Reynolds number suspension of self-propelled particles displaying non-reciprocal pairwise interactions. The particles perform active Brownian motion with a constant swim speed v and translational and rotational diffusion constants D and D R . Two types of pairwise interactions are present in the system as we describe in more detail in the "Methods" section. Firstly, the particles interact with each other via repulsive steric forces with a typical interaction range σ , which can be interpreted as an effective particle diameter. In the absence of other kind of interactions, the active suspension would be described by two dimensionless parameters: the area fraction of active particles and the Péclet number Pe = vσ/D . The Pe number is the ratio of characteristic times σ 2 /D and σ/v that the active particle needs to respectively diffuse and swim its own size σ . In this study we fix = 0.24 and Pe = 80 such that the system is in the regime where MIPS occurs.
In addition, the particles exhibit pairwise orientational interactions of short range R, which manifest themselves as torques acting on the particles. Motivated by the avoidant torque, the alga Chlamydomonas experiences close to a surface as explained below, we have formulated the orientational interaction of Eq. (5) in the "Methods" section. If the distance r ij between the particles i and j is smaller than a cutoff distance R, they are subjected to torques of the amplitude T 0 . In the general case, the orientational interaction is non-reciprocal, meaning that the torques T T T ij and T T T ji acting on the particles are not equal and opposite, T T T ij = −T T T ji . This stems from the form of Eq. (5) where the torque exerted on particle i depends on its orientation u i and the normalized separation vector r ij /|r ij | , but not on the orientation of particle j. The interactions are predominantly avoidant in nature. For example, interaction of such form can be motivated by looking at the dynamics of alga Chlamydomonas reinhardtii 39 . Chlamydomonas propagates itself through the fluid by breaststroke beating of a pair of frontal flagella 40 such that its far-field flow topology resembles that of a puller 41 : the flow is inward along the main body axis and outward in the perpendicular direction. Recent experiments have demonstrated that Chlamydomonas scatter off solid surfaces primarily due to contact interactions 42 between their flagella and the surface, while later investigations reveal that hydrodynamic interactions are needed for a complete quantitative description of how algae interact with cylindrical obstacles 43 . The contact interactions, on which we concentrate here, generate torques 44 that avert the algal body from a collision with the surface. To capture this type of interaction, we propose the model described by Eq. (5), in which case T T T ij ≡ T T T i would stand for the torque exerted on the particle and r ij ≡ r i for the distance vector of the particle i from the surface. It is now reasonable to suppose that an analog contact interaction through flagella exists among two algae situated close to each other. Since there are no useful experimental insights concerning such interactions between algae, we propose the torque of Eq. (5) as a first approach for describing these interactions, which certainly needs to be refined in future work for a quantitative description. Thus, as the flagellum of one Chlamydomonas touches the cell body of the other, and vice versa, the torques exerted on the algae will be non-reciprocal in the general case. Some typical examples of the encounter of two algae modeled by Eq. (5) are illustrated in Fig. 1, where the orientation of alga i is fixed and eight orientations of alga j are chosen. The directions of the torques acting on both algae and their strengths are indicated by curved red arrows. In particular, the avoidant torque on particle j is maximal if its orientation vector www.nature.com/scientificreports/ u j is perpendicular to the distance vector to particle i ( u j ⊥ r ij ) and, by symmetry, zero for u j r ij . The torque in Eq. (5) can be viewed as a simple and generic form of an avoidant non-reciprocal orientational interaction and, therefore, in the following we explore how it influences the collective motion of active Brownian particles. We describe active particle motion by a set of coupled overdamped stochastic equations, which are presented in the "Methods" section. Starting from an active suspension exhibiting MIPS as a reference state, we introduce non-reciprocal orientational interactions (5) and using the method of Brownian dynamics conduct a detailed parameter study of collective motion of the resulting active bath as a function of torque amplitude T 0 and interaction range R. As the orientational interactions are predominantly avoidant, one expects disruption of MIPS for strong enough torques T 0 , because active particles at the edge of a cluster turn away from neighboring particles and leave the cluster. Even though the disordered state is indeed the most common outcome for strong torques, as the state diagram in Fig. 2 shows, surprisingly, we find that there exists a narrow region in the state diagram, characterized by short interaction range R and moderate torques T 0 , where a transition to ordered collective motion occurs. In the ordered state, the particles move in flocks which can take a variety of structures spanning from a compact dense band to a single very dynamic porous band and multiple bands. This is the main result of our study.
The article is organized as follows. In the "Results" section we analyze different states appearing in our system and construct the corresponding state diagram in the parameter space of T 0 and R. We elaborate the implications of our findings, suggest possible future experiments, and offer our conclusions in the "Discussion" section. Lastly, the equations of motion of active particles and details of Brownian dynamics simulations are presented in the "Methods" section.

Results
The results of our Brownian dynamics simulations are summarized in a state diagram shown in Fig. 2. We observe four distinct states. First, for relatively weak torques, T 0 /k B T 0.5 , non-reciprocal interactions play only a minor role and the active bath exhibits MIPS irrespective of their range R. This regime is illustrated in the top left panel of Fig. 3 and in Video 1. Second, for sufficiently strong torques, T 0 /k B T 0.5 , and large interaction range R the non-reciprocal interactions prevent motility-induced clustering. A typical snapshot of the resulting dispersed state is shown in the top right panel of Fig. 3. Smaller unstable clusters are visible (see Video 2).
The most interesting behavior occurs when the range of non-reciprocal interactions is close to the range of steric interactions, in our case R/σ 1.3 . For moderately large torques, T 0 /k B T 5 , one observes a transition to coherent flocking motion of active particles. Depending on the exact parameter choice from this domain, the system can display an assortment of different flocking structures, which we roughly classify into two groups, which we will justify quantitatively below. In the region displayed in orange in Fig. 2 a dense band of particles develops, cruising collectively in a common direction against a dilute background. This state is depicted in the middle panels of Fig. 3. There, a coherently moving band is shown at two time instances, 3D −1 R apart from each www.nature.com/scientificreports/ other; see also Video 3. In addition to ordered bands, other structurally diverse states displaying coherent motion are possible in the system. Some examples are shown in the bottom panels of Fig. 3, which include a moderately dense stripe and multiple stripes. The dynamics of multiple stripes is shown in Video 4. Other irregular structures are also possible. They share the property that in contrast to the moving band, which always keeps its integrity, the stripes are very dynamic. As we demonstrate now, they can be grouped into one state, which we term flocking (dynamic) and which we show as blue region in the state diagram in Fig. 2. First, to quantify the degree of order in coherently moving states, we use the average normalized orientation of active particles as the polar order parameter where N is the number of particles in the system. If the individual particles point in random directions, the motion is disordered and the order parameter is approximately zero. On the other side, all particles pointing in roughly the same direction would correspond to ≈ 1 . In Fig. 4 we plot the polar order parameter as a function of torque amplitude T 0 for several interaction ranges R. The highest order is realized when the range R of non-reciprocal forces coincides with the range of steric forces 2 1/6 σ . Consistent with the findings of Fig. 2, the polar order is vanishingly small for T 0 /k B T < 5 in the MIPS and dispersed state. Around this torque value, abruptly increases over a narrow torque interval up to values 0.8 . Interestingly, a structural change in the state of coherent motion from a flocking band to other flocking patterns around T 0 /k B T = 20 is not www.nature.com/scientificreports/ accompanied with a change of , which stays approximately the same for a broad range of T 0 values, see Fig. 4. Increase in the range R of orientational interactions leads to a reduction of order in general. For sufficiently long interaction range (in our case R/σ 1.18 ) we observe that the order initially increases with T 0 towards some maximum value, after which decreases and eventually becomes vanishingly small with a further increase of torque amplitude in the reentrant dispersed state. Intriguingly, for all values of R in this domain, the maximum value of is roughly attained for torques from the interval 15 < T 0 /k B T < 20 . Finally, let us note that the range of torque values for which the system displays coherent motion (finite ) decreases with an increase of R. For R/σ ≥ 1.28 the order disappears for all T 0 . In the inset of Fig. 4 we plot the decrease of the maximum of the order parameter, max , with the increase of the interaction range R. Further work is needed to reveal if the decrease of max to zero in the vicinity of R/σ ≈ 1.28 is continuous or discontinuous. Second, for snapshots of the different states we perform a Voronoi tessellation of the particle configuration. From the local area A loc assigned to each particle, we calculate a local packing fraction or density φ = σ 2 π/(4A loc ) . For the whole time series of such snapshots, we then determine the density distribution p(φ) for a specific parameter set. Examples for the different states are plotted in Fig. 5 for R/σ = 1.18 and different interaction amplitudes T 0 . While the distribution for the MIPS state in plot (a) shows a low density peak and a very sharp peak at close-packing density (note, the close-packing peak is not at ca. 0.90 as expected for hard disks but rather at ca. 0.75 since we use the WCA potential of Eq. (4), with a potential value of ε = 100k B T at a particle distance of σ ; therefore, the effective particle diameter is larger than σ ), the dispersed state is characterized by a broad maximum of p(φ) at around φ = 0.15 (b). The density distribution of the flocking (band) state again has two peaks as shown in plot (c). However, compared to the MIPS state the maximum at low densities is sharper  www.nature.com/scientificreports/ and the second maximum, which belongs to the moving band, appears at a density smaller than the close-packing value and is broader. Finally, in the distribution for the flocking (dynamic) state in plot (d) the low-density peak is much less pronounced and the stripes are indicated by a broad shoulder or a very weak maximum (not shown).
We can now justify our classification of states quantitatively. In Fig. 6 we plot the polar order parameter versus the variance Var(φ) of the density distribution for constant interaction range R/σ = 1.18 and all the maximal torque values T 0 from Fig. 2. The points clearly group into the four states, already introduced. The flocking states have a large polar order parameter, while the compact bands and the dynamic stripes distinguish themselves by the variance of the density distribution, which obviously is larger for the bands. There is only one system at the transition from dispersed to flocking (band) ( T 0 = 7k B T ), which does not clearly fall into the latter state. The dispersed state and MIPS have nearly zero polar order, as already discussed, but MIPS has a larger variance due to the sharp peak at close-packing density.
Finally, we also checked the occurrence of the flocking states for two additional density values φ = 0.15 and 0.35 and we also wanted to know if the interaction range, where flocking occurs, expands. Since a complete scan of the parameter space is very time consuming, we chose the specific torque amplitude T 0 = 15k B T , where flocking occurs at the maximum interaction range of R/σ = 1.275 . For the higher density flocking sets in around R/σ = 1.3 and for the lower density below R/σ = 1.25 . So, there is not a very pronounced variation in the maximum interaction range.

Discussion
To summarize, we studied a system of active particles interacting with each other via conventional steric interactions; besides, we introduced the additional non-reciprocal avoidant orientational interactions. The latter were motivated by an example of flagellum-body interactions which are present in an encounter of two algae (cf. Fig. 1). We showed that motility-induced active particle clustering is disturbed already for moderately low non-reciprocal interaction amplitudes T 0 and that the active bath finds itself in most instances in a dispersed state, characterized by numerous and unstable smaller clusters. However, the dispersed state is not the only possible outcome. Remarkably, if the non-reciprocal interaction range R is close to the steric interaction range, we demonstrated that the active system displays coherent flocking motion for a finite range of amplitudes T 0 (cf. Fig. 4). The polar order is found to decrease with increasing R, and is largely independent on the underlying flock structure. Together with the variance of the local-density distribution, we can group the flocking motion into band and dynamic stripe states.
The occurrence of a band of flocking particles can at least qualitatively be understood as follows. When they all move in the same direction, the torques acting on each particle from its neighbors cancel. Now, whenever a particle orientation deviates from the direction of the moving band, the particle moves towards its direct neighbor and experiences the avoidant torque back to the moving direction. The same happens when the close-by neighbor on the other side is approached. Thus the parallel orientations of all the particles are stabilized by the avoidant orientational interactions with the neighbors; however, only when the torque amplitude exceeds k B T so that thermal orientational motion is irrelevant. If the interaction range R is too large, a particle experiences the avoidant torques from all the neighbors simultaneously. The torques cancel each other and the particle's orientation is not turned back to the common direction. Thus, the band disintegrates. Interestingly, it is also not stable for large torques and the reason for this instability under increased "orientational tension" requires further explanation. Thus, fully uncovering the microscopic mechanisms of the observed behavior will be one direction of our future efforts.
Collective behavior of Chlamydomonas algae has recently been the focus of experimental studies 45 . Importantly, the flocking motion observed in our simulation work has not been reported in the experiments so far. We hope that our work may stimulate future investigations aimed at establishing a connection between our www.nature.com/scientificreports/ parameters T 0 and R and interaction parameters of real systems. A possible experimental setup could be perhaps realized either through biological engineering (manipulating the length and mechanical properties of algal flagella) or by manufacturing artificial microswimmers following prescribed non-reciprocal interactions. The non-reciprocal avoidant torque of Eq. (5) has a very generic form, which one can modify. In particular, there is a symmetry in the torque values for a particle oriented towards or away from its neighbor. One can break this symmetry and make the avoidant torque larger when the particle is oriented towards its neighbor. Such a modification might, for example, be interesting to describe the social interaction of a prey turning away from a predator ( Supplementary Videos 1-4).

Methods
We consider a system consisting of N overdamped active particles propelling with a fixed speed v and variable orientations in two spatial dimensions. The time evolution of their positions r i = (x i , y i ) and orientations u i = (cos ϕ i , sin ϕ i ) is described by associated overdamped Langevin equations The particles possess translational and rotational mobilities µ and µ R and their translational and rotational diffusive motion is characterized by respective diffusion constants D = µk B T and D R = µ R k B T , with T being the ambient temperature. Here ξ i and η i denote two zero-mean and unit-variance Gaussian white noises, Greek letters denote Cartesian components. The particles interact with each other in two ways.
Firstly, the particles exhibit steric forces F ij = −∇ r i V ε,σ (r i − r j ) , originating from the Weeks-Chandler-Andersen potential V ε,σ (r) , which has strength ε at a characteristic distance σ: Secondly, the particles located within a distance R from each other experience pairwise orientational interactions. These interactions are modeled by subjecting the particle i, when it is close to particle j, to the torque with T 0 being the torque amplitude, R the interaction range, r ij = r i − r j and θ(r) the Heaviside step function. Note that the orientational interaction is non-reciprocal, T T T ij = −T T T ji , in the general case. Some representative examples of this interaction are sketched in Fig. 1.
We measure time in units of active particle reorientation time D −1 R , length in units of σ , and energy in units of k B T . Particle motion is simulated within a box of area A and is subjected to periodic boundary conditions. In the absence of orientational interactions, the state of the active bath is described by the area packing fraction of disks � = Nσ 2 π/(4A) and the Péclet number Pe = σ v/D.
We simulate N = 10 4 particles, and set ε/k B T = 100 and µ/µ R = σ 2 /3 , derived from the Stokes friction coefficients of spherical particles. The dimensionless parameters = 0.24 and Pe = 80 are chosen such that the bath exhibits motility-induced phase separation in the absence of orientational interactions, mathscrT 0 = 0 . The amplitude, mathscrT 0 , and range, R of the orientational interaction are then varied to explore the state diagram for fixed and Pe . The numerical integration of Eqs. (2)-(3) is carried out following an Euler scheme with a time step of 10 −5 D −1 R . Upon reaching the steady-state, the data are averaged over 5 independent simulation runs, each of duration 1000 D −1 R .