Clustering and flocking of repulsive chiral active particles with non-reciprocal couplings

Recently, non-reciprocal systems have become a focus of growing interest. Examples occur in soft and active matter, but also in engineered quantum materials and neural (brain) networks. Here, we investigate the impact of non-reciprocity on the collective behavior of a system of (dry) chiral active matter. Specifically, we consider a mixture of"circle swimmers"with steric interactions and non-reciprocal alignment couplings. Based on hydrodynamic equations which we derive from a set of Langevin equations, we explore the interplay of non-reciprocity, finite size, and chirality. We first consider, as a reference, one-species systems with reciprocal couplings. Based on a linear stability analysis and numerical simulations, we here observe three different types of collective behavior, that is, flocking, motility-induced phase separation, and a combination of both. Turning then to a non-reciprocal system, we find that non-reciprocity can turn otherwise stationary instabilities into oscillatory ones, affect the relative orientation of flocks, and, crucially, change the general type of instability. This illustrates the drastic impact of non-reciprocity and chirality on the emergent collective dynamics of active matter systems, with potentially far-reaching biological implications.

Given their non-equilibrium character and ubiquity in nature, various recent studies [1,2,5,[26][27][28] have addressed the collective dynamics of non-reciprocal (soft matter) systems based on a field-theoretical approach. While the detailed effect of non-reciprocity depends on the system considered, an overall finding is that nonreciprocity can drive time-dependent states. In fact, already for a two-component system of purely diffusive, conserved scalar fields, You et al. [2] have shown that nonreciprocity constitutes a generic route to traveling states. They found that a static demixed pattern can undergo a transition to a spatially inhomogeneous "run-andcatch" state, which breaks parity and time-reversal symmetry.
In the present work, we focus on active systems, whose constituents perform persistent motion due to an internal or external source of energy. Thus, even in the conventional case of reciprocal couplings, these systems are intrinsically out of equilibrium. It is now well established that active systems are capable of exhibiting a variety of non-equilibrium phase transitions and self-organization without external driving [29,30]. A growing number of studies is exploring the impact of non-reciprocal alignment between the active constituents. Examples include systems of (passive) dissenters in a flock of active particles [31] and generic phase transitions in nonreciprocal, active, two-species systems [26][27][28]. In particular, [26] demonstrated that non-reciprocity alone can destabilize the stationary (anti-)flocking state, characterized by (anti-)parallel motion of the particles of both species. The related phase transition is marked by an exceptional point in the space of field variables, resulting in a timedependent, so-called "chiral phase". Here, flocks of both species rotate at a constant speed with a fixed relative angle, although there is no intrinsic torque on the particle level.
Given these recent developments, the goal of our work is to examine the combined effect of two ubiquitous features of active matter systems, namely non-reciprocity and chirality of individual particles. Specifically, we consider a mixture of chiral active particles with non-reciprocal (anti-)alignment between particles of different species and mutual repulsion. In contrast to conventional "linear" swimmers that change their direction of motion only by diffusion or orientational (alignment) interactions (such as, e.g., active Brownian particles [32,33] or Vicsek particles [34,35]), chiral active particles (sometimes also named as "circle swimmers") additionally self-rotate with an intrinsic frequency [36,37]. The intrinsic rotation can be caused, e.g., by a chiral body shape as in anisotropic colloids [38], curved proteins [39] or artificial L-shaped particles [40,41]. In two dimensions, a chiral body shape indeed leads to a circular motion [36,37,42]. Other examples of chiral active particles are E. coli bacteria close to walls and interfaces [43][44][45][46], sperm cells [47,48], and particles actuated by rotating fields [49][50][51][52].
To examine the combined effect of non-reciprocity and chirality on the collective dynamics, we derive hydrodynamic equations for the density and polarization fields starting from the microscopic Langevin equations governing the motion of individual particles. We thereby employ a mean-field approximation and a truncation scheme to get rid of higher-order moments. While this strategy has been used before [26,31,53], including applications for circle swimmers [54][55][56], our approach additionally takes into account the impact of repulsive interactions.
The paper is organized as follows. In section 2, we start by introducing the microscopic model, followed by the derivation of the hydrodynamic equations. In section 3, we first consider a one-species chiral system with reciprocal alignment couplings, where we examine different scenarios of collective behavior and the impact of intrinsic rotation. We then turn to a two-species system with non-reciprocal couplings in section 4. Considering selected values of the intrinsic frequencies, we study in section 4, how non-reciprocity affects the collective behavior of the chiral mixture.
We close with a discussion of our results in section 5.

Particle-level description
We consider a two-dimensional system of chiral active particles comprising two species a = A, B. The N = a N a particles are located at positions r α (with α = i a = 1, ..., N a ) and move like active Brownian particles (ABP) with additional intrinsic torque. Thus, they rotate with a species-specific intrinsic frequency ω a , and self-propel with velocity v a along the instantaneous direction p α (t) = (cos θ α , sin θ α ) T , where θ α is the polar angle. The dynamics is then given by the overdamped Langevin equations (LE)ṙ where the sums over particles β = j b = 1, ..., N b couple the dynamics of particle α to the position and orientation of all other particles of both species b = A, B.
The translational LE (1a) involves the soft repulsive force within a cutoff-distance R r , where r αβ = |r αβ | = |r α − r β | and Θ(R r − r αβ ) = 1, if r αβ < R r , and zero otherwise. Note that, for simplicity, we assume that steric repulsion of strength k is the same for all particles (i.e. k and R r are speciesindependent). Further, we here assume a soft, piece-wise linear repulsive force [2,31], which allows for an analytical treatment in our continuum description (see section 2.2).
The rotational LE (1b) contains the torque, whose in-plane component is given by T ab al (r α , r β , θ α , θ β ) = K ab sin(θ β − θ α ) Θ(R θ − r αβ ), of strength K ab , which can be positive or negative. The sinusoidal dependence of the torque on the relative angle θ β − θ α is motivated by the expression derived from an orientation-dependent potential of the form p α · p β [53,62,66,74]. From Equation (3) it follows that particles of species a tend to orient parallel (align) or anti-parallel (anti-align) with neighboring particles (within radius R θ ) of species b when K ab > 0 or K ab < 0, respectively. For reciprocal couplings defined by the choice K ab = K ba , particles of species a align (or anti-align) with particles of species b in the same way as particles of species b with particles of species a. This is, in fact, the natural choice when the orientational coupling in Equation (3) is derived from a Hamiltonian, i.e. a many-body interaction potential. In the present work, we specifically allow for nonreciprocal orientational couplings, that is, K ab = K ba .
The mobilities are connected to thermal noise via µ r = β ξ and µ θ = β η, where β −1 = k B T is the thermal energy with Boltzmann's constant k B and temperature T .
Finally, we specify the propulsion velocity v a of individual particles appearing in LE (1a). In standard ABP-like models, one usually assumes that every individual particle self-propels with a constant speed v 0 . It turns out, however, that a constant self-propulsion speed combined with the mean-field approximation in our continuum approach does not reproduce the motility-induced phase separation (MIPS) into lowand high-density regions characteristic for ABPs [33,[75][76][77][78]. The underlying reason is that the mean-field approach disregards the structure of pair correlations, which become anisotropic due to activity. More specifically, previous numerical [58,63] and analytical [58][59][60]63] studies have shown that the force imbalance introduced by the self-propulsion of particles (there are more particles in front of the reference particle than behind it) causes an effective velocity reduction depending on the density of surrounding particles. In order to account for this effect, we replace the (speciesdependent) constant speed v a 0 with an effective density-dependent velocity of an active particle in an interacting repulsive system [53,[57][58][59][60][61]63]. Specifically, we assume that the particles self-propel with effective velocity [58][59][60][61]63] This choice expresses the fact that particles are slowed down in crowded situations, depending on the (species-independent) overall local particle density ρ = ρ(r) = a ρ a and velocity-reduction parameter ζ. Note that since we assume that particles of each species experience the same steric repulsion (see Equation (2)), each particle is slowed down with the same "frictional" parameter ζ, coupled to the overall density field ρ(r).

Coarse-grained description
Starting from the particle-level ("microscopic") Langevin equations (1a) and (1b), we derive a continuum model to study large-scale patterns and relevant mechanisms effecting the collecting behavior. To this end, we employ a coarse-graining strategy introduced in [79], yielding a mean-field Fokker-Planck equation for the one-particle probability density function (PDF) Relating the Fourier coefficients of f a (r, θ, t) to orientational moments, we then apply a closure relation neglecting higher-order moments to derive a hydrodynamic field description in terms of the particle density ρ a (r, t) = N a π −π f a (r, θ, t) dθ (6) and polarization density the latter measuring the overall orientation of particles at a certain position via w a /ρ a [80]. Our derivation closely follows the approaches of deriving a hydrodynamic version of the Vicsek model presented in [26] and of non-repulsive, reciprocally coupled chiral active particles in [54][55][56]. Details regarding the derivation of the hydrodynamic equations are given in Appendix A. Specifically, we find the continuity equation The flux involves the polarization density, which evolves according to In Equation (10), we have neglected the explicit terms of order O(∇w 2 ) and O(∇ρ∇w). The full equation is shown in Appendix B.
The density flux j a given in Equation (9) reflects that the motion of particles of species a in space arises from their self-propulsion in direction w a , whereby the particles are slowed down in crowded situations due to the density-dependent velocity.
Additionally, the flux comprises the drift of particles towards less crowded regions due to steric repulsion and translational diffusion. The change of the polarization density w a , described by Equation (10), originates from the competition between the tendency of particles to swim (with increasing speed) towards low-density regions (first term on r.h.s.), the rotation of the polarization with intrinsic frequency (second term), the decay of the polarization due to rotational diffusion (third term), and the orientational coupling of particles among all species (fourth term). The remaining (diffusional and non-linear) terms smear out low-and high-polarization regions.

parameter definition description
Pea v a 0 τ / Péclet number Ra 2 k µr R 4 r π ρ a 0 τ /(3 2 ) strength of the repulsive force z ζ ρ 0 τ / particle velocity-reduction due to the neighbors Dt ξ τ / 2 translational diffusion Ωa ωa τ intrinsic frequency g ab K ab µ θ R 2 θ π ρ b 0 τ /2 orientational coupling parameter equations are non-dimensionalized by choosing the Brownian time scale τ = 1/η as characteristic time scale and the particle radius as characteristic length scale. The particle and polarization densities of species a are scaled with the average particle density ρ a 0 . The remaining six dimensionless control parameters are the Péclet number Pe a = v a 0 τ / , R a = 2 k µ r R 4 r π ρ a 0 τ /(3 2 ) encoding the strength of the repulsive force, z = ζ ρ 0 τ / measuring the particle velocity-reduction due to the environment, the translational diffusion coefficient D t = ξ τ / 2 , Ω a = ω a τ for the intrinsic frequency, and g ab = K ab µ θ R 2 θ π ρ b 0 τ /2 as relative orientational coupling parameter. Thereby, g ab > 0 leads to an alignment and g ab < 0 to an anti-alignment of particles. The six control parameters are summarized in table 1.
To put the hydrodynamic equations derived here into the context of earlier continuum models for chiral active particles, we note that for one species of particles without steric repulsion (i.e. R a = z = 0), our hydrodynamic equations match those given by Liebchen and Levis in [55]. Further, when considering several different, non-mutually coupled species, our equations are in agreement with those derived by Fruchart et al. [26]. However, Fruchart et al. neither considered the presence of an intrinsic frequency (i.e. Ω a = 0) nor volume exclusion. In this work, we specifically allow for non-reciprocal couplings between different species of chiral particles, which do not only (anti-)align but also sterically interact due to the finite particle size.
In the present paper, we study the collective behavior based on the hydrodynamic Equations (8) -(10) following essentially two strategies. First, we employ a linear stability analysis starting from the uniform, rotationally isotropic state given by (ρ a , w a ) = (ρ a 0 , 0). This state corresponds, in fact, to the "trivial" solution of Equations (8) - (10). Second, we perform numerical simulations of the full hydrodynamic equations in two-dimensional periodic systems. For the numerical simulations, we use a pseudo-spectral code in combination with an operator splitting technique, allowing us to treat the linear operator exactly in the fourth-order Runge Kutta time integration. We choose the initial state to be a slightly perturbed disordered state of zero polarization w(r, 0) = 0 and constant density ρ(r, 0) = ρ 0 = 1.
The two-dimensional simulation box of size 100 × 100 is separated into 256 × 256 grid points.

One species
We start by investigating a system with one species (a = A) of repulsive chiral active particles, which allows us to focus on the effect and interplay of intrinsic frequency, steric repulsion, and reciprocal orientational couplings. On the basis of this section's results, we will then address the impact of non-reciprocity by studying a corresponding mixture in the subsequent section 4.
As expressed by Equation (11), the system is subject to perturbations involving all wave numbers k. We assume these perturbations to be of plane wave form with wave vector k, (complex) growth rate σ(k) and amplitudesρ(k) andŵ(k). Here, σ depends only on the wave number k = |k|, because we study the stability of the isotropic base state. We now insert the ansatz ρ(r, t) = ρ 0 + ρ (r, t), w(r, t) = w (r, t) into the evolution equations (8) and (10) (restricting them to one species) and assume ρ and w to be small. Linearization with respect to the perturbation then leads to a decoupling with respect to k. For each k, this yields an eigenvalue problem given by with where D = v 2 (ρ 0 )/(2 b) + D t , and effective velocity v(ρ 0 ) = Pe − z ρ 0 . From Equation (12), we can derive analytical expressions for the (complex) growth rates σ(k), which play the roles of eigenvalues. Note that the real parts of the three eigenvalues σ(k) determine the actual growth in time, whereas the imaginary parts are related to oscillatory behavior. We thus focus mainly on investigating the real parts, Re(σ).
In particular, the disordered state is linearly stable only if Re(σ(k)) < 0 for all k. In contrast, it is linearly unstable as soon as Re(σ(k)) > 0 for any k. In our investigation, we monitor all three functions Re(σ(k)). To interpret the behavior, we here assume that the largest value and corresponding eigenvector determine the type of emerging dynamics at short times. This assumption is later checked by means of numerical continuum simulations.
The eigenvalue equation (12) and, in particular, matrix (13) already allow us to deduce some general features of the linear dynamics. First, we note that the growth rates at k = 0 play a special role as they are generally related to the change of the spatially integrated value of the hydrodynamic quantities, corresponding to longwavelength perturbations. Importantly, such a perturbation can occur only in the polarization (signaling flocking) but not in the particle density. To see this, we recall that ρ(r, t) is conserved by means of the continuity equation (8), such that As seen from Equation (12), this implies that the growth rate of the corresponding mode at k = 0 must vanish. Hence, our system generally features at least one σ(k = 0) = 0. The two remaining eigenvalues at k = 0 are related to long-wavelength fluctuations of the polarizationŵ(k = 0) = w (r, t) dr, whereby Re(σ(k = 0)) > 0 or Re(σ(k = 0)) < 0 indicate that such a fluctuation grows or is suppressed.
Investigating the stability at finite k, on the other hand, is equivalent to taking the self-propulsion of swimmers (v(ρ 0 )) into account. This is because in matrix (13), all terms proportional to k or k 2 , which stem from gradient terms, are only non-zero for v(ρ 0 ) > 0, expect for the diffusional terms (∼ D t , R).
From matrix (13) we can further extract the parameters with the most important effects on the linear stability of the system. Specifically, it is seen that the translational diffusion coefficient D t and steric repulsion parameter R only appear in the diagonal entries and scale with k 2 , ensuring the stability at large wave numbers. Further, the average density ρ 0 only appears as a prefactor of R, z, and g AA , and thus has no direct effect itself. The Péclet number Pe only appears in v(ρ 0 ), such that only the difference between Pe and z ρ 0 has an impact on the dynamics, but not Pe itself.
Hence, all possible scenarios are governed by parameters z, g AA , and Ω A .

Results
As the analytical expressions for σ(k) at arbitrary frequencies Ω A are long, high-order polynomials, we here consider the exemplary case of Ω A = 0.    figure 2). Here, Re(σ(k)) ≤ 0 for all k, reflecting that perturbations of any type and at all wave numbers decrease in time. As soon as Re(σ) becomes positive, collective dynamics start to emerge.
real part of eigenvalues, Re(σ(k)) order parameters ∀ r, t after transient regime disordered We recall that, in an active fluid, such a polarization implies ordered motion of the particles. To find the corresponding conditions, we consider the complex growth rates at k = 0, The corresponding eigenvectors v(k = 0) = (ρ 0 ,ŵ x,0 ,ŵ y,0 ) T are v 1 (k = 0) = (1, 0, 0) T and v 2/3 (k = 0) = (0,ŵ x1,2 ,ŵ y1,2 ) T , whereŵ x1,2 ,ŵ y1,2 ∈ C. Hence, the eigenvalue σ 1 (k = 0) reflects the conservation of the particle density, whereas σ 2/3 (k = 0) indicate the change of the overall polarization. As seen from (15b), a flocking instability can only occur for strong alignment couplings characterized by g AA ρ 0 > 1. Furthermore, while the real parts of growth rates (15) are independent of the intrinsic frequency Ω A , the imaginary parts are not. Thus, Ω A does not affect the flocking instability itself, but makes the instability oscillatory. Indeed, as shown later in numerical simulations (section 3.2), the flocking phase is time-dependent in the sense that w(r, t) rotates in time. Microscopically, this means that the chiral active particles synchronize, such that they rotate coherently with a frequency, which is not necessarily the same as the intrinsic one. Within the flocking regime, the full real growth rates look like shown in figure 1(b) with two Re(σ(k = 0)) > 0, which decrease for increasing k. Since this type of instability is not affected by volume exclusion effects (modeled by parameters R and z), Liebchen and Levis [55] observed the same flocking instability (15b) for non-repulsive chiral systems.  Figure 2. Stability diagram obtained for the one-species system of repulsive chiral active particles with alignment interactions. We exemplarily show the result for intrinsic frequency Ω A = 0.1. Depending on the orientational coupling strength g AA and velocity-reduction parameter z, the system can exhibit different types of instabilities. The disordered state (blue) is stable for weak alignment, g AA ρ 0 < 1, or anti-alignment, g AA < 0. Increasing the alignment between particles (g AA ρ 0 > 1), the system undergoes a flocking transition (yellow). For specific values of z, the system exhibits MIPS -either pure for weak alignment (red) or in combination with flocking for strong alignment (orange). Analytically estimated stability conditions for z are plotted as white squares. The white crosses correspond to parameters used in the numerical continuum simulations in section 3.2. Other parameters are Pe = 1.5, R = 0.1, Dt = 0.01, and ρ 0 = 1.

MIPS
Taking volume exclusion into account, the second type of instability is related to (pure) MIPS, where the system separates into high-and low-density regions (see red regions in figure 2). The real growth rates typical for pure MIPS are shown in figure 1(c). The two eigenvalues with Re(σ(k = 0)) < 0 indicate that fluctuations of the overall polarization are suppressed. From (15b) it then follows that pure MIPS can only occur when the alignment coupling is weak (i.e. g AA ρ 0 < 1) or particles antialign (g AA < 0). We further see that the eigenvalue related to density fluctuations increases towards positive values at small k and exhibits a maximum at finite k > 0.
Such a maximum indicates a characteristic length scale of emerging patterns at small times. Importantly, we observe this type of instability only in a certain range of velocity-reduction parameters z. In fact, in case of non-chiral active particles, we can analytically compute a condition for z as outlined in Appendix C. For chiral active particles, we use the same approach (that is expanding the eigenvalues up to order k 2 ) to estimate values of z for which one eigenvalue becomes positive. This yields the white squares in figure 2. A more detailed discussion of the role of Ω A on MIPS is postponed to section 3.3. Here, we first address the effect of increasing alignment interactions g AA .
As mentioned before, pure MIPS can only occur when g AA ρ Re(σ(k)) at k = 0 (flocking) and a maximum at finite k > 0 (MIPS).
We note that a similar phenomenon, namely the emergence of clusters consisting of synchronized particles, has also been observed in the absence of steric repulsion at sufficiently large frequencies [55]. However, this so-called "microflock instability" is not related to the simultaneous emergence of flocking and MIPS reported in the present paper. Different to [55], we here consider steric repulsion between particles, leading to density-dependent velocity reduction and, consequently, MIPS. In contrast, the clustering in [55] appears as a "secondary" instability not of the isotropic state (considered here) but of the flocking state.

Numerical continuum simulations
While the linear stability analysis is a convenient tool to analyze fluctuations around the base state considered, we can neither use it to deduce the dynamics at long times nor the exact form of emerging patterns. In principle, these questions could be solved by performing particle-based simulations of the microscopic Langevin equations (1a) and (1b). However, extensive parameter scans and long simulation times would make a quick exploration of parameter spaces rather difficult. We therefore complement the previous analysis with numerical simulations of the full, non-linear hydrodynamic equations (8) - (10). Specifically, we consider the three parameter sets indicated by the white crosses in figure 2. This choice corresponds to parameters which lie well within the predicted stability regions of MIPS, flocking, and MIPS combined with flocking, respectively.
The emerging particle density and polarization density fields are time-dependent.
In the following, we focus on times after the initial transient regimes (t > t t ), whereby the time t t , after which the respective system has passed the transient regime, depends on the parameters. The subsequently shown snapshots represent instantaneous fields at t > t t . For an illustration of the actual time dependence, we provide videos of the different phases in the supplemental material.
We start by discussing the flocking phase. Representative snapshots of the spatially resolved particle density ρ(r, t) and the corresponding time-dependent  polarization density field w(r, t) are shown in figure 3(a). Figure 3(b) shows the probability distribution p(ρ) of the local density. As expected from our stability analysis, within the (pure) flocking state, ρ(r, t) is constant, as indicated by the single sharp peak in p(ρ). The spatially-averaged polarization density field rotates in time and can be described by Its amplitude is given by the (time-independent) absolute value measuring the ordering of particles and, hence, the "strength" of the flock formation.
Further, Ω flock = Ω(r, t) is the spatially-averaged rotation frequency of the flock.
In the flocking phase, | w | is non-zero (| w |/ρ 0 ≈ 0.65), reflecting ordered motion induced by sufficiently strong local alignment of particles (compare (15b)). These observations comply with the linear stability analysis (see figure 1(b)). The stability analysis further predicts that the emerging state is time-dependent since Im(σ(k = 0)) = 0 (see Equation (15b)), whereby the imaginary part is given by the intrinsic frequency Ω A . Indeed, we find in our numerical simulation that the direction of polarization rotates in time in counter clockwise direction with rotation frequency Ω flock = 0.075. Since Ω flock is the spatially-averaged rotation frequency with variance of O(10 −7 ) and stays constant after a transient regime (t > t t ), we deduce that the entire flock rotates with the same frequency Ω flock . This flock rotation frequency is somewhat smaller than the intrinsic frequency of each particle, Ω A = 0.1, which matches very well the prediction of linear stability analyses around a flocking state in reciprocal one-species chiral systems in [55]. Microscopically, the non-zero value of Ω flock implies that particles synchronize and thereby rotate in a coherent manner. For non-chiral active particles (Ω A = 0), numerical simulations show that the flock retains its direction of motion, as one would expect.
We now turn to (pure) MIPS. In this case, numerical simulations of the full hydrodynamic equations yield particle and polarization density fields as shown in the snapshots in figure 4(a). From ρ(r, t), we observe the emergence and growth of clusters with local density larger than ρ 0 . In fact, the cluster formation is a consequence of self-trapping mechanisms due to a reduction of self-propulsion velocity in crowded situations (e.g., [63,85]). As a result of the cluster formation, p(ρ) has two peaks: one at low density and the other one at larger density ( figure 4(b)). We also find that the  Equation (16), we can numerically extract the spatially-averaged rotation frequency Ω flock = 0.037 with variance var(Ω) = 0.004. The overall behavior is reminiscent of the "interrupted motility-induced phase separation" observed in experiments, which show that alignment of active Janus colloids interrupts the phase separation process, eventually leading to a fluctuating but non-increasing average cluster size [75].
Taken altogether, the results of our numerical simulations of the full, non-linear hydrodynamic equations are consistent with the predictions of the stability analysis regarding the type of emerging collective behavior. This holds even when both order parameters are involved, as in the case for the combined MIPS and flocking instability.

Effect of intrinsic frequency
We now come back to the question of how the intrinsic frequency Ω A of the chiral active particles affects MIPS. So far, this has been studied by particle-based simulations but not via a hydrodynamic theory [64]. Here, we investigate this question on the basis of a linear stability analysis. Results for the eigenvalues as functions of k are presented in figure 6, where we consider three values of alignment coupling.
In systems without alignment interactions (g AA = 0), the instability related to (pure) MIPS remains for small intrinsic frequencies Ω A = 0.25 ( figure 6(a)).
In this range, the (relevant) growth rate Re(σ(k)) is zero at k = 0 and has a maximum at finite wave number. However, when the particles rotate with a larger frequency (Ω A = 0.5, 0.75, 10), the maximum disappears and the disordered phase becomes stable. Thus, we conclude that the intrinsic frequency of chiral particles generally opposes MIPS. This prediction is consistent with results from particle-based simulations of non-aligning chiral active particles by Liao et al. [64], who showed that MIPS only occurs for small intrinsic frequencies.
Switching on the alignment coupling between the particles, but keeping it weak (e.g. g AA = 0.5), the linear stability analysis still predicts pure MIPS for slowly rotating active particles ( figure 6(b)). In fact, the range of Ω A where MIPS occurs is even extended relative to the case g AA = 0, as seen from the curve pertained to Ω A = 0.5 in figure 6(b).
Further increasing the alignment strength to g AA = 1.5, the system undergoes a flocking transition, characterized by Re(σ(k = 0)) > 0 for all Ω A considered. As

Two species
We now turn to a binary system (a = A, B) of chiral active particles, focusing on the impact of non-reciprocal alignment couplings. To this end, we employ again the previously established combination of linear stability analyses and numerical solutions of the full continuum equations (8) - (10). Clearly, the two-species system involves a large set of parameters, making a complete investigation of the full parameter space a rather overwhelming task. Here we therefore focus on some representative parameter combinations. We recall in this context that several parameters characterizing the two species individually have already set equal (this concerns the steric repulsion, the self-propulsion velocity, and diffusion). The main control parameters for the twospecies system are therefore the alignment coupling parameters g ab and the intrinsic frequencies Ω a . Using the results obtained in section 3 for the one-species system as a reference, our main goal is to explore the impact of non-reciprocity.

Linear stability analysis
In a binary system, the trivial solution to the hydrodynamic equations (8) - (10) is given (as in the one-species system) by a disordered state with zero polarization w A (r, t) = w B (r, t) = 0 and constant density ρ A (r, t) = ρ A 0 , ρ B (r, t) = ρ B 0 . Assuming again that the perturbations scale with ∼ e σ(k)t (compare Equation (11)), and linearizing, we obtain the eigenvalue equation In Equation (18), the eigenvector v = (ρ A ,ŵ A x ,ŵ A y ,ρ B ,ŵ B x ,ŵ B y ) T is now sixdimensional, containing the possible perturbations of the particle densities and the two components of the polarizations for each species. Further, the 6×6 matrix M 2 (k) can be written as a combination of two types of submatrices, Here, P a (k) involves the couplings within each species a = A, B with D = v 2 a (ρ 0 )/(2 b) + D t and effective velocity v a (ρ 0 ) = Pe a − z (ρ a 0 + ρ b 0 ). The other submatrix involves the couplings between the two species (ab = AB, BA) Note that, in case of non-reciprocal alignment couplings, g AB = g BA , and different intrinsic frequencies, Ω A = Ω B , the entire matrix M 2 (k) becomes non-symmetric.
We assume (as in the one-species case) that the six resulting eigenvalues indicate the instability in the two-species system. At k = 0, the growth rates are given by where we have chosen ρ A 0 = ρ B 0 = ρ 0 /2 for simplicity. The first two growth rates (22a) vanish due the conservation of the particle density. The real parts of the other four eigenvalues can become positive for strong orientational alignment coupling, see (22b).
As in the one-species case, this generally indicates a flocking instability. The flocking instability becomes oscillatory for non-zero imaginary parts of the eigenvalues (22b). This generally happens when Ω a = 0.
In addition to (pure) flocking, we also observe pure MIPS and combined flocking and MIPS. To identify these instabilities, we employ the same characteristics as outlined for the one-species system in section 3 and summarized in table 2.
We now discuss some specific features occurring in the two-species system. Maybe most intriguingly, the flocking instability can become oscillatory in the non-reciprocal binary system even in the absence of chirality, Ω a = 0. Specifically, this happens for "antagonistic" inter-species couplings (g AB g BA < 0) with In this case, the function C defined in Equation (23)  we assume that particles within both species weakly align (g AA = g BB = 0.5) and we chose a velocity-reduction parameter of z = 0.375 Pe a /ρ a 0 . In the one-species system with ρ A 0 = ρ 0 /2 and Ω A = 0.1, these parameters would result in a pure MIPS instability. Starting from this scenario, we now vary the inter-species coupling strengths g AB and g BA . Our results regarding the linear stability of the system are summarized in the diagram in figure 7. The exemplarily shown real parts of the growth rates in figure E1 in Appendix E share much similarities with the one-species case (figure 1).
When g AB = g BA (diagonal solid white line), the orientational coupling in the binary system is reciprocal. Therefore, when moving along the diagonal line, we only vary the strength of inter-species alignment (g ab > 0) or anti-alignment (g ab < 0).
The system becomes non-reciprocal if g AB = g BA . Note that the stability diagram is indeed symmetric under the exchange g AB ↔ g BA . At k = 0, this directly follows from the structure of Equations (22) and (23). For a more general discussion of this point, see Appendix H. analysis predicts either a stable disordered phase (blue) or MIPS (red, figure E1(b)).
In particular, the disordered phase is stable for weak inter-species anti-alignment, while the MIPS emerges for weak alignment.
To illustrate, in particular, the effect of non-reciprocity, let us consider the Clustering and flocking of repulsive chiral active particles with non-reciprocal couplings23 reciprocal situation generating MIPS combined with flocking (right upper corner in figure 7) as a starting point. Then, we gradually increase non-reciprocity by moving horizontally in the stability diagram towards smaller g AB while keeping the same g BA .
The linear stability analysis predicts that MIPS combined with a flocking instability first changes to pure MIPS until, eventually, the disordered phase is stabilized. This demonstrates the fact that non-reciprocity alone can dramatically change the character of the instability.
The results discussed so far pertain to a specific choice of intrinsic frequencies.
As we show in Appendix F, a different choice slightly shifts the instability regions, while the qualitative picture remains.

Numerical continuum simulations
Performing numerical simulations of the full, non-linear hydrodynamic equations (8) -(10) allows us to complement the linear stability analysis of the two-species system. Like in the one-species case, we choose parameter sets which lie well within the stability region of the respective phase transition in figure 7 (black crosses). In particular, we choose a weak intra-species alignment of g AA = g BB = 0.5. We now explore the effect of non-reciprocity by varying the inter-species coupling strengths g AB , g BA . Corresponding simulation videos can be found in the supplemental material.
Numerical simulation results of MIPS in a non-reciprocal two-species system.
(a) Snapshots of the time-dependent particle density ρ a (r, t) and polarization density field w a (r, t). is, as in the one-species system, reminiscent of "interrupted motility-induced phase separation" in aligning particle systems [75]. Within the clusters, both species form rotating flocks. Each individual flock can be described by Equation (16) Moving away from this reciprocal case by decreasing g AB relative to g BA , the instability changes to pure MIPS (upper middle cross in figure 7). Corresponding snapshots presented in figure 9(a) reveal that, indeed, both species form growing patterns with vanishing polarization. However, due to the non-reciprocal couplings between the species, the emerging patterns differ from each other. In particular, for the chosen parameters of antagonistic couplings with g AB = −0.5 and g BA = 1.5, particles of species A weakly anti-align with particles of species B, while the latter strongly align with particles of species A. The opposite goals combined with weak intra-species alignment prevents flocking of either species. Nevertheless, due to the stronger magnitude of |g BA | > |g AB |, one might assume that species B "wins" the competition, resulting in enhanced alignment of species B as compared to species A.
As already weak particle alignment promotes MIPS (see section 3.3), this could be the reason for the enhanced cluster formation of species B, while particles of species A accumulate at the edges of the clusters. This observation is also reflected in p(ρ a ) in figure 9(b). The peak at large densities of species B indicates the cluster formation, whereas species A is characterized by a more uniform density with a narrow peak at a slightly enhanced density. These numerical results match the linear stability analysis, which predicts the different cluster growth within the two species by means of the eigenvectors corresponding to the positive eigenvalues.
A further decrease of g AB relative to g BA (upper left cross in figure 7) eventually stabilizes the disordered phase. Here, the antagonistic couplings, g AB = −g BA = −1.5, are of equally strong magnitude, preventing not only flocking but also MIPS of either species.
Lastly, we turn to the anti-flocking regime (yellow region in figure 7). As a reference, we present in figure 10(a)representative snapshots of ρ a (r, t) and w a (r, t) (a) (b) Figure 11. Numerical simulation results of anti-flocking in a non-reciprocal two-species system after the initial transient regime, t > tt. (a) Snapshots of the time-dependent particle density ρ a (r, t) and polarization density field w a (r, t). White arrows indicate the instantaneous, local direction of w a . in the case of strong reciprocal anti-alignment of strength g AB = g BA = −1.5 (lower cross in yellow region in figure 7). Zoomed-in snapshots at two different times are shown in figure G2 in Appendix G. Both species form time-dependent patterns with non-vanishing polarization yet with different magnitude. Contrary to the previous cases of inter-species alignment (g AB , g BA > 0), where both species accumulated at the same places, particle clusters of different species seem to avoid each other in case of inter-species anti-alignment (g AB , g BA < 0). At the same time, particles of the slower rotating species A accumulate into denser clusters than particles of species B.
Even though the linear stability analysis does not predict MIPS, we observe cluster formation, which is, in accordance with our previous results (see section 3.3), weaker for larger intrinsic frequencies. The probability distributions p(ρ a ) in figure 10(b) further support these observations. Flocks of both species rotate with the same spatially-averaged frequency Ω A,B flocks = 0.19 > Ω A (variance var(Ω) = 0.02), which is time-independent (after the initial transient regime). Thereby, the anti-alignment makes the two flocks move in rather anti-parallel direction. The linear stability analysis predicts a relative angle of Θ lin = 172.3°, which agrees very well with the observed angle of Θ obs (t) = 172.8°∀ t > t t .
Introducing non-reciprocity while staying in the anti-flocking regime (upper cross in the yellow region in figure 7), we observe time-dependent patterns as shown in We close with two comments: The first one concerns the symmetry of the stability diagram in figure 7. Having seen the numerical simulation results, it is not surprising that although exchanging g AB ↔ g BA does not affect the linear stability of the system, the full, non-linear dynamics is not symmetric under the exchange (see also Appendix H).
The second point is related to the time-dependent phase termed "chiral phase" by Fruchart et al. [26]. This phase (not to be confused with the rotating flocking phases in our system) emerges for strong enough antagonistic couplings (see Equation (24)). To illustrate that our hydrodynamic model also captures these cases, numerical simulation results for a non-chiral system (Ω A = Ω B = 0) with stronger intra-species alignment are shown in Appendix I.

Conclusion
In this paper, we have studied the collective behavior of chiral active particles, interacting via volume exclusion and orientational couplings. In particular, we have allowed for non-reciprocal (anti-)alignment between active particles of different species.
We have started from a particle-level description of the swimmers in terms of Langevin equations, in which an effective, density-dependent propulsion velocity accounts for the trapping of particles in crowded situations due to steric repulsion.
We have then derived the corresponding coarse-grained description under the meanfield assumption and a scaling ansatz for higher-order orientational moments. The resulting hydrodynamic equations consist of a continuity equation for the (conserved) particle density field and a second time evolution equation for the (non-conserved) polarization density. These equations have been analyzed by linear stability analyses around the homogeneous, isotropic state, and by numerical solutions of the full, nonlinear equations.
We have first focused on the effect and interplay of the intrinsic frequency, steric repulsion, and reciprocal orientational couplings by studying the one-species system.
Here, we found three different types of non-trivial collective behavior: flocking, MIPS, and flocking combined with MIPS. Different to flocking, MIPS is strongly affected by the intrinsic frequency of the constituents. In particular, the intrinsic frequency generally opposes MIPS. However, this effect can to some extend be compensated by an increase of alignment strength, which promotes the emergence of MIPS. With this, our work presents the first hydrodynamic results for MIPS in systems of aligning chiral active particles. Despite the mean-field character of our hydrodynamic approach, the predictions qualitatively agree with particle-based simulations performed in earlier studies [64].
To explore the effect of non-reciprocal couplings between particles, we have then turned to a two-species system. Our results demonstrate that non-reciprocity has indeed a significant impact on the collective dynamics. Consistent with recent field-theoretical results for non-reciprocal systems (e.g., [2,26]), otherwise stationary instabilities can become oscillatory when couplings between particles are antagonistic.
The (repulsive) system considered here thus adds an important example of an active soft-matter system with non-trivial time-dependent states, here generated by the interplay of non-reciprocity and chirality. Additionally, non-reciprocal inter-species couplings affect the relative orientation of the formed flocks. However, the most severe effect of non-reciprocity is that it can even change the general type of instability.
In the present paper, we have mainly focused on a qualitative description of the emerging behavior in repulsive chiral active systems with non-reciprocal orientational couplings. As stated earlier in the paper (see section 3.2), it would be desirable to complement our hydrodynamic results with particle-based simulations of the underlying Langevin equations. In this way one could avoid the (mean-field like) approximations in the derivation of the continuum equations, which would allow to study the characteristics of the individual phases and transitions, such as the active self-assembly into flocks, in more detail. Of course, the price to pay are much larger computational costs needed to find relevant parameter sets and analyzing the results.
The present calculations could serve as a guideline for parameter sets and phenomena to be investigated. A further interesting question concerns the nature of the nonequilibrium transitions (or bifurcations in the continuum picture) separating different states. Work in these directions has just started (on the basis of simpler, non-reciprocal models [27,28]).
We also note that our findings are relevant for real chiral active mixtures such as anisotropic colloid systems [38], bacteria close to walls [43][44][45][46] or sperm cells [47,48]. Due to naturally occurring heterogeneity and, in particular, couplings mediated through non-equilibrium environments, non-reciprocity is indeed pervasive in active matter systems. Moreover, as pointed out in [1], non-reciprocity has farreaching biological implications as it might be crucial in order to understand directed information transmission in living systems.  [26,[54][55][56].
We then sum over all (N a ) particles α, take the ensemble average, and employ the mean-field approximation to approximate the two-particle PDF f ab (r, r , θ, θ , t). To treat the spatial integrals in Equation (A.2), we assume that the coupling ranges R r , R θ are small such that only near-by particles interact. Changing then to the particle distance d = |r − r| as new integration variable, we can perform Taylor expansions of the terms inside the integrals up to first order around d = 0. These expansions allow us to perform the spatial integration analytically. As a result we obtain the Fokker-Planck equation Note that the same Fokker-Planck equation can be obtained by determining drift and diffusion coefficients using the Kramers-Moyal expansion of the distribution function, as described for general settings in [81]. Importantly, one would also need the same mean-field assumption (A.3) as employed in this work. Following [80], the calculation of the remaining orientational integral in the , a lengthy but straightforward calculation yields the time evolution of the Fourier modesf a n (r, t) = π −π f a (r, θ, t) e inθ dθ, reading ∂ tf a n = − +i ω a nf a n − ξ ∂ z ∂ zf a n − η n 2f a n , where we omitted the (r, t)-dependence of the Fourier modesf a n . The identification of complex numbers with two-dimensional vectors allows us to relate the Fourier modes Clustering and flocking of repulsive chiral active particles with non-reciprocal couplings32 to moments of the one-particle PDF f a (r, θ, t). In particular, we can identify the particle density (related to mode n = 0) measuring the probability of finding a particle of species a at position r and time t, and the polarization density (related to mode n = 1) describing the average orientation of particles of species a via w a /ρ a .
The time evolution (A.6) of the Fourier modesf a n represents a hierarchy of equations, which requires a consistent closure scheme. Here we employ the scaling ansatz proposed in [80,82], which has already been used in numerous active matter systems, e.g. [26,31,83,84], including chiral active particle systems in [54][55][56]. The scaling ansatz assumes, first, that deviations from the isotropic state are so small that we can neglect moments of order n = 3 or higher (f a n = 0 for n ≥ 3). Second, the scaling ansatz assumes that the nematic order parameter, which is related to the second modef a 2 , changes adiabatically, i.e. ∂ tf a 2 = 0. These assumptions allow us to expressf a 2 solely in terms of the first Fourier modef a 1 . Specifically, we obtain from As a result of the closure relation, the full dynamics of the one-particle PDF is reduced to the dynamics of the particle density and polarization density. The concrete coupled equations (8)-(10) for the particle density ρ a and polarization density w a follow after some lengthy, but straightforward calculations and non-dimensionalization.
Appendix C. Motility-induced phase separation in weakly-aligning non-chiral active systems A useful reference for the behavior of one-species chiral active systems discussed in section 3.1.2 is the case of zero intrinsic frequency (Ω A = 0), which allows for analytically feasible expressions. The three eigenvalues of the non-chiral active system where and v(ρ 0 ) = Pe − z ρ 0 .
In the regime of weak alignment, i.e. g AA ρ 0 < 1, the eigenvalues σ 1 (k) and σ 3 (k) are negative at all k. In contrast, σ 2 (k) can be of typical MIPS form with σ 2 (k = 0) = 0 and a positive maximum at a finite wave number. To see this, we expand eigenvalues the (C.1) up to second order in k, yielding Thus, in the systems of non-chiral active particles with weak alignment interactions (g AA ρ 0 < 1), we can expect pure MIPS to occur for velocity-reduction parameters z 1 > z > z 2 with The expression on the right-hand side of Equation (C.5) resembles the one given by Sesé-Sansa et al. [63] for active particles without alignment interactions (i.e. g AA = 0).
Here, we non-dimensionalize our equations with the rotational diffusion strength η, such that in our instability condition (C.5), only the translational diffusion R ρ 0 + D t appears. Reference [63] uses a different re-scaling, such that z 1/2 depends on the product of rotational and effective translational diffusion coefficient (D r D).
The stability diagram in figure C1 shows in two-species systems can be obtained by looking at the respective eigenvectors. In particular, we here assume that the eigenvectors corresponding to the largest positive real part of the growth rates indicate the direction in the space of dynamical variables, in which perturbations grow the fastest. We further assume that the flock orientation of the individual species is given byŵ A 0 =ŵ A (k = 0) andŵ B 0 =ŵ B (k = 0). Then, the relative angle Θ betweenŵ A 0 andŵ B 0 suggests whether the flocks are oriented parallel, anti-parallel or enclose a certain angle with each other. Some general remarks regarding the relative orientation of flocks can be deduced when we consider two exemplary choices of the intrinsic frequencies, namely Ω A = Ω B and Ω A = −Ω B . For simplicity, we further set g AA = g BB and ρ A 0 = ρ B 0 = 1. When Ω A = Ω B , the eigenvalues at k = 0 are given by (see Equation (22b)) Thus, a flocking instability occurs for inter-species couplings fulfilling the condition ±Re √ g AB √ g BA > 1 − g AA . The corresponding eigenvectors are Since the total densities of both species are conserved quantities, the corresponding components are zero. To determine the relative angle between the flocks, we focus on the remaining components and calculate the angle between (Re(ŵ A x ), Re(ŵ A y )) T (flock A) and (Re(ŵ B x ), Re(ŵ B y )) T (flock B). The resulting relative angles of the growing eigenvectors are Hence, if Ω A = Ω B and g AB g BA > 0, the two flocks are either exactly parallel or exactly anti-parallel -independent of whether the inter-species couplings are reciprocal or non-reciprocal. Only when the intra-species couplings are large, specifically, g AA = g BB > 1, a flocking instability can occur for g AB g BA < 0 with relative angle Θ = π/2. These results are in accordance with findings by Fruchart et al. [26], who showed that for non-chiral active systems, inter-species couplings of the same sign always lead either to exactly parallel flocking or exactly anti-parallel anti-flocking.
For opposite chiralities, Ω A = −Ω B = Ω, the (degenerated) eigenvalues are given by (see Equation (22b)) such that a flocking instability occurs as soon as ±Re g AB g BA − Ω 2 A > 1 − g AA . As in the previous case of same chirality, we focus on the real part of the eigenvector components to calculate the relative orientation of flock A and flock B. For g AB g BA > 0, the relative angle is given by (D.6) Clustering and flocking of repulsive chiral active particles with non-reciprocal couplings37 Equation (D.6) reveals that in the case considered here (opposite chiralities), the angle between the flock directions can be different from 0 or π or π/2. The possibility of flocks, which move under a relative angle, has already been found in reciprocal chiral active mixtures [56,86]. Here, we show that such a "mutual flocking phase" can also occur in non-reciprocal systems with opposite chiralities, as long as g AB g BA > 0.
However, for strongly non-reciprocal couplings with g AB g BA < 0, the resulting relative angle of Θ = π/2 is again not affected by the coupling strengths.
In fact, comparing the predicted relative angles between flocks with those obtained from numerical continuum simulations, we find good agreement at short times.
Appendix E. Real parts of growth rates in two-species system 0 2 4 6 8 10 k  (a) (b) Figure G3. Numerical simulation results of anti-flocking in a non-reciprocal twospecies system after the initial transient regime. Left: Representative snapshots of the time-dependent particle density ρ a (r, t) and polarization density field (c) Figure H1. Numerical simulation results and growth rates in case of MIPS in a non-reciprocal two-species system. (a) Snapshots of the time-dependent particle density ρ a (r, t) and polarization density field w a (r, t). The particle clusters grow in time. White arrows indicate the instantaneous, local direction of w a . (b) Probability distribution p(ρ a ) and absolute value of spatially-averaged polarizations | w a | in regions of large densities (larger than right dashed vertical line) for species A and overall for species B. Inter-species coupling strengths are g AB = 1.5, g BA = −0.5. (c) Largest real parts of growth rates for g AB ↔ g BA . Other parameters are gaa = 0.5, z = 0.375 Pea/ρ a 0 , Ω A = 0.1, Ω B = 0.5, Pea = 1.5, Dt = 0.03, ρ a 0 = 1, and R = 0.1 with a = A, B.
As mentioned in section 4.1.1, the instabilities predicted by our linear analysis are symmetric under the exchange g AB ↔ g BA , even though the chiral active particles rotate with different frequencies Ω A = Ω B . This is explicitly shown in figure H1(c), where we plot an example for the dominant growth rate. Clearly, Re(σ(k)) is fully symmetric at all wave numbers.
However, when performing simulations of the full, non-linear hydrodynamic equations (8) -(10), we see differences in the snapshots of both systems in figure 9(a) and figure H1(a). While both systems form clusters with vanishing polarization, the probability distributions of particle densities differ (see figure 9(b) and figure H1(b)).
Appendix I. Numerical simulation of a time-dependent "chiral" phase in non-chiral, non-reciprocal two-species systems In numerical simulations of our full hydrodynamic equations (8) -(10), we can observe a time-dependent "chiral" phase already in the absence of intrinsic rotation Ω A = Ω B = 0 (i.e., in non-chiral active systems), when couplings between species are antagonistic (see Equation (24)). This intriguing phenomenon, predicted also in earlier field-theoretical studies [26], is illustrated for the present system in figure I1.
As shown in the two snapshots in figure I1(a) and (b), the resulting particle density and polarization density fields are time-dependent, whereby the direction of motion within the flocks continuously changes. In fact, after the initial transient regime, the flocks of both species move under a relative angle of Θ obs (t) = 89.9°∀ t, which is very well predicted by the linear stability analysis with Θ lin = 90°(see Appendix D).
(a) (b) Figure I1. Numerical simulation results of MIPS combined with flocking in a non-chiral, non-reciprocal two-species system. Snapshots of the time-dependent particle density and polarization density field at time (a) t = 490 τ . (b) t = 500 τ .