Emergent vortices and phase separation in systems of chiral active particles with dipolar interactions

Using Brownian dynamics (BD) simulations we investigate the self-organization of a monolayer of chiral active particles with dipolar interactions. Each particle is driven by both, translational and rotational self-propulsion, and carries a permanent point dipole moment at its center. The direction of the translational propulsion for each particle is chosen to be parallel to its dipole moment. Simulations are performed at high dipolar coupling strength and a density below that related to motility-induced phase separation in simple active Brownian particles. Despite this restriction, we observe a wealth of phenomena including formation of two types of vortices, phase separation, and flocking transitions. To understand the appearance and disappearance of vortices in the many-particle system, we further investigate the dynamics of simple ring structures under the impact of self-propulsion.


Introduction
Systems of self-propelled particles consist of a large number of motile constituents, each of which is capable of continuously converting energy from an internal source or the surroundings into mechanical motion. 1,2 Examples of biological self-propelled particles can be found over a wide range of length and time scales, from bird flocks, fish schools, mammalian herds, and pedestrian crowds in our daily life, to bacteria, sperm cells, and microtubules at the microscale. 3,4 Self-propelled particles can also be synthesized in the laboratory, famous examples including bimetallic nanorods, 5,6 spherical Janus particles, 7,8 and magnetic rollers. [9][10][11] It is now well established that already relatively simple systems of self-propelled particles can display complex collective behavior, giving rise to a still increasing scientific interest. 12 A prominent example of such complex behavior is motilityinduced phase separation (MIPS) which occurs in systems of disk-shaped, active Brownian particles above a critical density Φ crit . [13][14][15][16][17] Another "classical" example is the flocking transition in the Vicsek model, 18 a system of self-propelled point-like particles with ferromagnetic interactions.
The majority of theoretical and numerical studies on active particles assumes that each particle tends to "swim" along a straight line due to the translational self-propulsion. However, this assumption clearly becomes invalid once the chiral symmetry along the propelling direction is broken. Such broken symmetry often introduces a rotational self-propulsion, which, together with a translational self-propulsion, causes a single swimmer to move along a perfect circular path in the absence of thermal noise. 19 Hence, active particles simultaneously subject to the rotational and translational self-propulsion are generally referred to as "chiral active particles," 4 or "circle swimmers." 20 Famous examples include E. coli cells close to a surface, 21,22 FtsZ proteins, [23][24][25] and synthetic L-shaped particle. 26 Recent research has shown that chiral motion (i.e., circle swimming) can indeed significantly change the self-assembly dynamics in active systems. Already for the simplest models of chiral active particles, it has been shown analytically and numerically a Institut für Theoretische Physik, Technische Universität Berlin, Hardenbergstr. 36, D-10623 Berlin, Germany. E-mail: guo-jun.liao@campus.tu-berlin.de, klapp@physik.tuberlin.de that circle swimming generally suppresses MIPS, 27 induces clockwise vortices, 28 and yields hyperuniform structures with vanishing long-wavelength (but large local) density fluctuations. 29 Recent studies of circle swimmers with short-range anisotropic interactions have revealed even more intriguing features. For example, simulations have demonstrated that chiral active particles with an asymmetric 24,30,31 or elongated 32 shape display vortex structures. Chiral active particles with chemotactic alignment can form a global traveling wave. 33 Further, chiral active particles with polar (i.e., ferromagnetic Heisenberg-like) interactions exhibit rotating macrodroplets, 34,35 chimera states, 36-38 chiral selfrecognition, 39 dynamical frustration, 40 or a surface-population reversal effect in ring-shaped confinement. 41 Additional effects can arise in the presence of a rotating magnetic field, where diskshaped chiral active particles without anisotropic interactions can display chiral separation and polar ordering. 42 In the present study, we go one step further and consider chiral active particles with dipole-dipole interactions. These interactions differ from the previously considered anisotropic ones by the longer range and by a more complicated angular dependency. That is, the pair interaction depends not only on the angle between the two involved dipole moments, but also on their spatial configuration which eventually promotes chain formation. To explore the role of these factors for the resulting collective behavior is interesting not only from an academic point of view. Rather, many microswimmers have embedded (permanent or induced) magnetic or electric dipole moments giving rise to dipole-dipole interactions, examples being metallodielectric Janus spheres, 7,43 magnetotactic bacteria, 44,45 and magnetic rollers. [9][10][11] Therefore, the overall aim of this study is to enhance the understanding of the active matter systems involving circle swimming and dipoledipole interactions.
To this end, we consider a two-dimensional system of chiral active particles with dipolar interactions, each of which moves at a self-propulsion speed (characterized by the particle motility) and rotates at an "active" angular speed. The (permanent) dipole moment embedded in each particle is assumed to be directed along the direction of translational self-propulsion. Our model combines features of non-dipolar, disk-shaped circle swimmers 28 and dipolar active Brownian particles. 46 Therefore, we refer to our model as dipolar circle swimmers in the rest of this paper. Based on Brownian dynamics simulations, we explore a wide range of the particle motility and the active angular speed at a large dipolar coupling strength and a low density.
As a unique feature of dipolar circle swimmers, we observe two types of vortices, which we refer to as Type I and Type II vortices. Type I vortices consist of forward-moving particles whose dipole moments display counterclockwise patterns, while Type II vortices are backward-moving particles whose dipole moments show clockwise structures. As the motility increases from zero, Type II vortices vanish, and the system exhibits significant phase separation with the dense domain characterized by giant Type I vortices. We show that some aspects of the vortex formation can be explained by considering the dynamics of simple ring structures. At even higher motilities, the system displays flocking patterns, in which dipolar circle swimmers self-assemble into polar clusters and show a significant global orientational order.
The rest of this paper is organized as follows. In Sec. 2 we describe our methods of investigation, including the model, the simulation setup, and the target quantities. An overview of the system behavior (at a fixed density) is given in Sec. 3.1. In Sec. 3.2 − 3.5, we discuss specific aspects such as chain formation, clustering and phase separation, emergence of vortices, and orientational ordering. Finally, we summarize our findings and provide an outlook in Sec. 4.

Model and methods of investigation
Our model system of dipolar circle swimmers combines the previously investigated model of the dipolar active Brownian particles 46 and the model of spherical circle swimmers. 19,28,29,[47][48][49] Therefore, various methods of investigation can be transferred from previous works. In the following, we summarize the main points and refer the reader for details to ref. 28, 46.

Model
We consider N disk-shaped Brownian particles with diameter σ dispersed in a monolayer in the xy-plane. Each of the particles carries a permanent point dipole moment µ µ µ i located at its center. In addition, each particle is subject to a self-propulsion force F 0ê e e i (i = 1, . . . , N) and torque M 0ẑ z z. The particle orientationê e e i denotes the direction of self-propulsion force and is assumed to be directed along the unit dipole momentμ µ µ i at each instant of time, i.e.,ê e e i =μ µ µ i . The pairwise interaction potential for two swimmers i and j is given by u pair r r r i j , µ µ µ i , µ µ µ j = u sr r i j + u dd r r r i j , µ µ µ i , µ µ µ j , where u sr r i j with r i j = |r r r i j | = |r r r j − r r r i | stands for the shortrange (sr) steric repulsion, which prevents particles from overlapping. Specifically, we employ the Weeks-Chandler-Anderson potential 50 defined as The potential is truncated at a cut-off (c) distance r c = 2 1/6 σ , with σ being the particle diameter. The repulsive strength is described by ε * = β ε with β −1 = k B T representing the thermal energy (with k B denoting the Boltzmann constant and T the temperature).
The second term on the right-hand side of Eq. (1) represents the (long-range) dipole-dipole (dd) interaction, given by The strength of the dipole-dipole interaction is characterized by the parameter λ = β µ 2 σ −3 , where µ = |µ µ µ i | denotes the strength of each dipole moment.

Equations of motion
To investigate the dynamical self-assembly of dipolar circle swimmers, we perform extensive Brownian dynamics (BD) simulations involving N particles in a squared box with side length L. The motion of the ith particle is described by the coupled Langevin equations in the overdamped limit, 19 where r r r i denotes the particle position and ψ i represents the polar angle for the orientationê e e i = (cosψ i , sinψ i ). Focusing on ψ i as the relevant angle, we implicitly assume that the dipole moment µ µ µ i are oriented along in-plane directions (for justification, see Sec. 2.3). The dots above r r r i and ψ i on the left-hand side of Eq. (4) and (5) indicate time derivatives, and the potential energy for the ith particle U i is the sum over all the pairwise potentials between particle i and all the other particles j, i.e., with u pair being the pair potential [see Eq. (1)].
In Eq. (4), D t denotes the translational diffusion constant (we do not use the tensorial quantity due to the disk-like shape of the particles). Correspondingly, for the rotational motion, D r describes the rotational diffusion constant. To account for the Brownian motion, the random force ξ ξ ξ i (t) and torque Γ i (t) are Gaussian white noises, which have zero means and are delta correlated, ξ ξ ξ Here, the angle brackets · · · stand for ensemble average, and ⊗ denotes dyadic product. Finally, it is worth recalling the behavior of a single circle swimmer (without interacting with other swimmers). In the absence of thermal fluctuations and at M 0 > 0, the swimmer moves along a closed circular path of radius R = D t F 0 / (D r M 0 ) and rotates counterclockwise at a self-propulsion speed v 0 = β D t F 0 and angular speed ω 0 = β D r M 0 . 19

Parameter choice
Following ref. 46, we choose the strength of steric repulsion as ε * = 10 and D r = 2.57914D t /σ 2 . The dipolar coupling strength is set to λ = 10, such that the dipolar pair energy of two particles at contact is ten times the thermal energy. To model chain formation and orientational ordering in dipolar colloids, one typically considers dipolar coupling strengths between λ ≈ 1 and λ ≈ 10. [51][52][53] However, test simulations of the present model with λ varying from 0 to 10 show that vortices (which will be discussed in Sec. 3.4) only appear once λ ≈ 10. A further increase of λ > 10 requires even smaller time steps δt, which, together with the long-range feature of dipole-dipole interactions, makes the simulations even more computationally expensive. Therefore, in the present work the dipolar coupling strength is fixed to λ = 10. For passive monolayers of such strongly coupled dipolar particles, it is well known that there is a pronounced tendency to orient along in-plane directions; specifically, one observes self-assembly into chains and rings, [54][55][56][57][58] or dense planar ordered states. 54,59-61 Therefore, we here assume beforehand that the dipole moment µ µ µ i lies in the xy-plane, i.e., the z-component is neglected. Based on this assumption, a two-dimensional Ewald summation is employed to deal with the long-range feature of dipole-dipole interactions, as outlined in the Appendix of ref. 46.
Equations (4) and (5) are solved numerically via the Euler−Maruyama method 62 with a discrete time step δt = 2 × 10 −5 τ, where τ = σ 2 /D t denotes the Brownian diffusive time. We define the particle diameter σ , the thermal energy k B T , and the Brownian diffusive time τ as the units of length, energy, and time, respectively. All physical quantities in the system are then expressed in the units based on the dimensional combination of these three basic units. Consistent with our earlier work, 28 the impact of the translational and rotational self-propulsion with respect to the thermal noise is represented by the dimensionless motility v * 0 = v 0 σ /D t and the active angular speed ω * 0 = ω 0 /D r . Moreover, we set the particle number to N = 1156, unless otherwise stated. The sizes of the present simulations are limited to the order of 10 3 particles, due to the expensive computational cost resulting from the long-range feature of dipole-dipole interactions and the necessity of using small time step δt = 2 × 10 −5 τ for simulating active particles with strong steric repulsion (ε * = 10).
As the initial configuration for all simulations, we distribute particles on a square lattice and assign a random orientation to each of them. For reaching a steady state, we typically wait at least 5 × 10 5 steps. Then, we start to measure the statistical properties (see Sec. 2.4) every 500 steps and collect at least 1000 samples for each realization.
Our main focus of this study is on the pattern formation of dipolar circle swimmers in the regime of low densities. Specifically, the mean area fraction is chosen to be Φ = Nπσ 2 e f f / 4L 2 = 0.23, where the effective hard disk diameter is given by σ e f f ≈ 1.07851σ (see ref. 46 for details). We note that the density considered here is below the critical density of motility-induced phase separation (MIPS), Φ crit , in purely repulsive systems, which is in the range 0.28 Φ crit 0.4. 63,64 In other words, at Φ = 0.23 and in the absence of dipole-dipole interactions, MIPS does not occur even at extremely high motilities v * 0 .

Observables
In this section we discuss the relevant observables and their numerical analysis (more details can be found in ref. 28,46).

Percolation and polymerization
It is well established that passive dipolar spheres with strong dipole-dipole coupling have a tendency to aggregate into chains, rings, and percolated networks. 51,55,56,[65][66][67][68][69][70][71][72][73] Therefore, we expect similar patterns to emerge in the present system, at least at low motilities and angular speeds. A first measure is the percolation probability Π, defined as the probability of finding a cluster in a simulation snapshot that connects two opposite sides of the simulation box. 65 Following earlier studies of active colloidal systems, 14,74,75 a cluster is defined as follows: Two particles are regarded as being associated if their center-to-center distance is smaller than a "critical" distance r L . A cluster is then a set of particles that are mutually associated. Specifically, we set r L = 1.2σ .
As it turns out, the results do not change significantly if we choose a different value of r L , as long as it corresponds to a distance between the location of the first peak and the first valley of the pair correlation function for the corresponding "reference system" (v * 0 = 0, ω * 0 = 0, and λ = 0) at Φ = 0.23. Further, for the present system of dipolar circle swimmers, the location of the first peak and the first valley do not vary too much upon variation of v * 0 and ω * 0 . Therefore, we fix r L = 1.2σ throughout this work for simplicity.
The percolated dipolar networks may break, e.g., upon an increase of temperature. 68 In this case, the system may display a "string" fluid state with many "polymerized" chains, which are composed of dipolar spheres connecting with their neighbors in a head-to-tail fashion. In self-assembly studies of dipolar systems, these chain structures are commonly quantified by the degree of "polymerization," 46,64,68,76 In the above equation, N p denotes the number of particles in the "polymerized" chains. In the present study, a chain is defined as a cluster comprised of at least ten particles, which conforms with the large value of the dipolar coupling strength (λ = 10).
We note that the above definition of a chain cannot distinguish between an elongated chain and a compact disk-shaped cluster with local head-to-tail ordering. This problem cannot be solved by adding more complex criteria to the definition of clusters (as ref. 46 does), i.e., µ µ µ i · µ µ µ j > 0 and (µ µ µ i · r r r i j )(µ µ µ j · r r r i j ) > 0, since headto-tail ordering is present in all types of aggregates. The degree of polymerization φ p defined in Eq. (7) approaches unity if all particles self-assemble into string-like chains or disk-shaped clusters, and is zero if there is no chain structure. We consider the system to display a state with polymerized string fluids if φ p > 0.5, and no giant disk-shaped clusters (to be discussed later) are present. Chain Clustering, phase separation Orientational State formation and emergent vortices ordering φ e e e ≤ 0.5 Vortices with phase separation d e ≥ 3, double peaks in P (φ ) φ e e e ≤ 0.5 Micro-flocking -monotonic decay of nP (n) φ e e e > 0.5 Macro-flocking -broad shoulder or a peak φ e e e > 0.5 in nP (n) at large n Table 1 Characterization of the states of dipolar circle swimmers according to the target quantities described in Sec. 2.4.

Clustering and phase separation
As mentioned in Sec. 2.3, the density considered in this work is below the critical density of the MIPS occurring in non-dipolar active Brownian particles. Nevertheless, it has been shown that some anisotropic (e.g., polar) interactions can enhance the tendency for MIPS. 74,77 It thus seems worth checking for this phenomenon also in the present system. As a first step, we investigate the clustering behavior by measuring the fraction of the largest cluster 14,16,28,46,74,78 where n lc denotes the size of the largest cluster. The definition of a cluster is given in Sec. 2.4.1. The order parameter φ c reaches unity if the largest cluster is composed of all swimmers, and approaches zero if swimmers are homogeneously distributed. For the present dipolar system, we note that both, the fraction of the largest cluster φ c and the degree of polymerization φ p (see Eq. (7)), are close to unity if a single giant cluster is present. On the other hand, φ c 1 and φ p ≈ 1 indicate that the system forms many chains with intermediate size. We also note that a large value of φ c alone cannot distinguish between a state with percolated networks and a state with giant disk-shaped clusters. To solve this, we monitor, in addition, the percolation probability Π. Percolated networks are identified by φ c > 0.5 and Π > 0.5, while giant compact clusters are identified by φ c > 0.5 and Π ≤ 0.5.
Finally, to identify phase separation, we employ a Voronoi tessellation to obtain the probability distribution function P (φ ) of the local area fraction φ without a short-time average (see ref. 46 for details). The system is regarded to display phase separation if the density profile P (φ ) shows a double-peak structure. Subsequently, the coexisting densities are the density values corresponding to the peaks.

Emergent vortices
For active systems on a surface, a vortex structure is generally defined as a disk-shaped high-density region in which the particle velocities (or the coarse-grained velocity field) display circular patterns with a common center. Vortices emerge, for instance, in various systems composed of biological motile constituents, including microtubules 79 and bacterial suspensions. [80][81][82] Further, biological circle swimmers, such as FtsZ filaments, 23,83 and spermatozoa of sea urchins, 84 are found to be able to self-organize into an array of vortices. Vortex patterns can also be observed in synthetic colloidal systems of magnetic rollers. [9][10][11] Given these examples, it seems possible that also the present system of dipolar circle swimmers may self-organize into vortices, owing to the interplay between the active rotation and the dipole-dipole interactions. Inspired by ref. 80, we characterize the vortex structures by analyzing the orientational correlation function (recall that the direction of the particle's dipole moment coincides with its orientation,μ µ µ i =ê e e i ), If i = j, r i j = |r r r j − r r r i | = |r r r i − r r r i | = 0 andμ µ µ i ·μ µ µ j =μ µ µ i ·μ µ µ i = 1. Since particles do not overlap, we have r i j > 0 for i = j. Thus, we obtain C e e e (r = 0) = 1. If particle i and j are infinitely far apart (r i j → +∞), we assume that their average orientational correlation is very weak and is thus approximately zero at the given low density. Since C e e e (0) = 1 and lim r→+∞ C e e e (r) ≈ 0, a negative correlation at a finite distance r represents that particle orientations (i.e., propulsion directions) on average point in the opposite directions when separated by the distance r. Hence, the system is regarded to self-assemble into vortex structures if C e e e (r) shows negative values at finite distances. The distance corresponding to the minimum of C e e e (r) can be used to characterize the average diameter of vortices d e e e , defined by dC e e e (r) dr r=d e e e = 0.

Flocking behavior
Orientational ordering is commonly observed in a variety of active systems with anisotropic interactions. As a prominent example, the Vicsek model describes self-propelled particles whose velocities tend to align with those of their neighbors when perturbed by noise. 18 At high motilities, this model typically exhibits the socalled flocking behavior, in which particles self-organize into clusters with significant polar ordering and move collectively toward a certain direction. Moreover, the flocking states are found to persist in the presence of active rotation 34,35 and steric repulsion. 85 Although the dipole-dipole interactions in the present model are more complex than the Heisenberg-like interactions in the Vicsek model, we still expect to observe flocking patterns in the present system. This expectation is also motivated by our observation of flocking in dipolar active Brownian particles. 46 We characterize the emergence of flocking behavior by the global polar order parameter, 86 This order parameter represents the magnitude of the average orientation, which reaches unity if all particles move toward the same direction, and zero if the particle orientations are not correlated.
Depending on the length scale of polar clusters with respect to the particle number N, the flocking states can be further classified into micro-and macro-flocking states. To this end, we follow the treatment in ref. 46 and measure the distribution function of the cluster size P (n), which is the probability that a randomly selected cluster is composed of n particles. In the case of macro-flocking, the system is usually composed of several large polar clusters with many small polar clusters, such that P (n) is extremely small at large n. To solve this problem, we consider the weighted distribution function nP (n), whose value is proportional to the probability that a randomly selected particle belongs to a chain with a size n. We identify the swimmers are in a micro-flocking state if φ e e e > 0.5 and the characteristic cluster size in the weighted distribution function nP (n) does not scale with the particle number N. In contrast, the dipolar circle swimmers display a macro-flocking state if φ e e e > 0.5 and the size of the giant clusters corresponding to the peak structure in nP (n) scales with N.  Table 1. Overlapping symbols denote bistable states. The region surrounded by the dashed lines indicates a parameter regime where the system is percolated.

Results and discussion
In this section we present results from extensive BD simulations for a wide range of motilities (0 ≤ v * 0 ≤ 100) and angular speeds (0 ≤ ω * 0 ≤ 10) at a fixed area fraction Φ = 0.23 and dipolar coupling strength λ = 10. In the subsequent subsection 3.1, we first give an overview of the observed states, and then discuss in Sec. 3.3-3.5 various aspects in detail.

Non-equilibrium state diagram
Based on the targeted quantities described in Sec. 2.4, we have identified six states whose characteristics are summarized in Table 1. The resulting non-equilibrium state diagram spanned in the v * 0 , ω * 0 plane is shown in Fig. 1, and representative simulation snapshots are given in Fig. 2.
In the passive limit (v * 0 = ω * 0 = 0), the system displays percolated networks [see Fig. 2(g)], which is in agreement with simulation studies for a monolayer of passive dipolar spheres. 55,56 For a finite motility v * 0 and in the limit of linear swimmers (ω * 0 = 0), our model reduces to the model of dipolar active Brownian particles. 46 Upon increasing v * 0 this system undergoes a transition from a string fluid state characterized by chain-like structures [see Novel behavior emerges as the angular speed ω * 0 becomes nonzero. At vanishing motility (v * 0 = 0), a slight increase in ω * 0 from zero leads to the emergence of two types of vortex patterns, which will be discussed in detail in Sec. 3.4. As a "side effect" of the vortex formation, phase separation occurs at angular speeds ω * 0 = 2 − 3 [see Fig. 2(h)]. A further increase in ω * 0 renders the vortex patterns less pronounced, suppresses the phase separation, and breaks the percolated structures, as shown in Fig. 2(i).
Considering a finite motility of v * 0 = 40, the swimmers exhibit a string fluid state at small angular speeds (0 ≤ ω * Here we note that the parameter window of the angular speed ω * 0 for the vortex formation is narrower than the low motility case (v * 0 < 80). For fast rotation such as ω * 0 7, the system exhibits a macro-flocking state, where the sizes of flocking patterns are comparable to the simulation box [see Fig. 2(c)]. Interestingly, we also observe bistable states involving vortices and macro-flocking at v * 0 = 100 and ω * 0 = 5 − 6. This bistability will be further discussed in Sec. 3.3.

Chain formation
We start our detailed discussion of the system's behavior by focusing on the bottom left part of the state diagram in Fig. 1, where the key phenomenon is chain formation. We begin by investigating the percolation probability Π and the degree of polymerization φ p defined in Sec. 2.4.1. Figure 3(a) shows the percolation probability Π as a function of the motility v * 0 for various angular speeds ω * 0 . At v * 0 = 0 and ω * 0 = 0, Π approaches unity, which corresponds to the percolated networks well known for the passive case. 65 Upon an increase of v * 0 from zero at zero angular speed (ω * 0 = 0), Π quickly decays and reaches a value close to zero when  . Above this line, the system is percolated or polymerized, respectively (see Table 1). The solid lines are guides to the eye. v * 0 40. The corresponding snapshot is shown in Fig. 2(d). † For finite (fixed) values of ω * 0 , Π decreases more drastically with v * 0 , and this decrease becomes the more pronounced the larger ω * 0 is. This suggests that both, the motility and the active rotation, tend to suppress the percolation behavior (relative to the passive case).
Once Π has dropped to values below 0.5, the percolated networks have dissolved, yielding a fluid state with chains that do not span over the simulation box. To quantify these chain structures, we calculate the degree of "polymerization" φ p defined in Eq. (7). We recall that φ p alone cannot distinguish long string-like chains and compact disk-shaped clusters. To avoid possible confusion, in the following we only present the results for φ p when indeed chain structures are formed. Results plotted in Fig. 3(b) show that, indeed, the polymerization φ p is close to unity at zero motility and gradually decreases as v * 0 increases. Once v * 0 80, the chain structures essentially disappear. We also find that the degree of polymerization φ p is essentially unaffected by rotation for the small values of ω * 0 considered here (ω * 0 ≤ 2). The observed decrease of φ p at small ω * 0 (and Φ = 0.23) is, qualitatively, in a good agreement with earlier studies of dipolar active particles, where ω * 0 = 0. 46 Once ω * 0 > 2, the dipolar circle swimmers may self-assemble into patterns distinct from chain structures, such as † The arrows in Fig. 2(d) are chosen to be three times as large as the diameter of a disk-shaped swimmer for better visualization. However, this choice may create an unrealistic visual effect that some neighboring chains in Fig. 2(d) seem to be connected and span over the simulation cell. In fact, it is not the case, since the percolation probability Π approaches zero at v * 0 = 40 and ω * 0 = 0 [see Fig. 3(a)].
compact clusters, which we will discuss in the following sections.

Clustering and phase separation
It is well established that non-dipolar active Brownian particles can exhibit clustering and even phase separation at sufficiently large motilities (and appropriate densities). 14,16,28,74,78 In the absence of active rotation or alignment interactions, [14][15][16]78 the clustering behavior is purely induced by the interplay between the particle motility and the steric repulsion.
To quantify the clustering behavior in the present system, we first calculate the fraction of the largest cluster φ c , see Eq. (8) and Fig. 4. In the passive limit (v * 0 = 0 and ω * 0 = 0), one finds φ c ≈ 1. This is due to the large value of λ , which leads to strong headtail-alignment and thus, to the formation of percolated networks, with percolation probability Π ≈ 1 [see Fig. 2(g) and Fig. 3(a)]. Upon an increase of ω * 0 from zero at v * 0 = 0, the order parameter φ c decreases (as does Π), indicating that the percolated networks are suppressed by the active rotation. This seems plausible, since the active rotation tends to destroy the head-to-tail alignment. This leads to a breaking of percolated chains and, hence, a decrease of φ c .
At finite v * 0 the situation changes. Inspecting Fig. 3(a) again, we see that at v * 0 = 20, the system is only partially percolated (with Π 0.5) for ω * 0 2, and if v * 0 40, the system is not percolated at all (with Π ≈ 0) for all explored ω * 0 . We infer from the data that for v * 0 20, the values of φ c > 0.5 seen in Fig. 4 truly indicate the formation of giant compact clusters, rather than that of percolated networks. We also find from Fig. 4 that for v * 0 20, φ c varies non-monotonically with ω * 0 . We interpret this interesting observation such that giant clusters only appear at intermediate angular speeds (2 ω * 0 6), whereas the system is rather homogeneous (φ c ≤ 0.5) for slow and fast rotation. Further, with increasing motilities v * 0 , the window of ω * 0 values corresponding to giant clusters is shifted toward larger ω * 0 . We conclude that at finite ω * 0 , the motility tends to suppress (rather than enhance) the formation of giant clusters. A particularly complex situation occurs at very high motilities. For example, at v * 0 = 100 we observed that the system of dipolar circle swimmers can change within a single simulation from one state with φ c ≈ 1 (indicating giant clusters) to another state with φ c 0.5, or the other way around (snapshots not shown here). To check whether both of these states are steady states, we created multiple long-time realizations at v * 0 = 100 and various ω * 0 . Indeed, two different steady-state results were found from independent simulation realizations in the range ω * 0 = 5−6, indicating a bistability. Specifically, the realizations with φ c ≈ 1 are characterized by giant clusters, while those with φ c 0.5 correspond to a macro-flocking state, where macroscopic swarming is observed. The flocking behavior of dipolar circle swimmers will be discussed later in detail in Sec. 3.5.
Given the formation of giant clusters at suitable combinations of v * 0 and ω * 0 , it is interesting to explore whether this leads to phase separation. To this end, we plot in Fig. 5 the probability distribution function P(φ ) of the local area fraction φ for various angular speeds ω * 0 , taking the motility v * 0 = 40 as an example. Phase separation is indicated by a double-peak structure of P(φ ), where the coexisting densities are the density values corresponding to the peaks. From Fig. 5 we see that at zero and small angular speeds (ω * 0 1), the distribution function displays only a single peak located at φ ≈ 0.1, indicating that the system is essentially homogeneous. It is noted that at v * 0 = 40 and ω * 0 = 0 or 1, the system exhibits a state with string fluids (see Fig. 1), i.e., particles tend to form short linear chains. As a result, the distribution function is not symmetric and shows a pronounced tail at high densities. In contrast, at intermediate angular speeds (3 ω * 0 9), we observe that P (φ ) exhibits two well-defined peaks, showing that the dipolar circle swimmers phase-separate into dilute and dense domains. Finally, at a large angular speed ω * 0 = 10, the second maximum at high densities is only weakly pronounced, indicating that the dense domains almost disappear. Therefore, we expect that phase separation will eventually vanish upon further increase of ω * 0 . Finally, we plot in Fig. 6 the coexistence densities in the (v * 0 , φ ) plane for various angular speeds ω * 0 . At an intermediate angular speed ω * 0 = 3, the system displays phase separation for a broad range of the motility, v * 0 ≈ 0 − 80. Further, the density difference, ∆ φ = φ d − φ g , between the dense and the gas-like region first increases with v * 0 from zero to 20, and then decreases upon further increasing v * 0 from 20 to 80. This non-monotonic behavior of ∆ φ may be attributed to the fact that the parameter regime for ω * 0 = 3 and v * 0 = 20 − 80 is very close to the state boundary between string fluids and vortices with phase separation (see Fig. 1). As ω * 0 increases from 3 to 5, the motility range where phase separation occurs is shifted toward larger motilities (v * 0 = 20 − 100) with ∆ φ increasing monotonically with v * 0 . At large angular speeds such as ω * 0 7, the range of the motility shrinks to 30 v * 0 50. Moreover, the area enclosed by the curves of coexisting densities decreases with ω * 0 , suggesting that phase separation is in general suppressed by ω * 0 . The suppression of phase separation is also consistent with the vanishing peak of P(φ ) at high densities for ω 0 = 10 in Fig. 5. Indeed, at extremely large ω * 0 each of the particles tends to swim along a small circular path. Thus, it quickly alters its propulsion direction, causing the large dense domain to "melt" and break into small pieces. Similar observations regarding the impact of ω * 0 on phase separation have been made for systems of non-dipolar circle swimmers. 28

Emergent vortices
In this section we explore in detail an intriguing "byproduct" of the phase separation discussed in Sec. 3.3, namely, the formation vortices occupying the dense domain [for an illustration, see Fig. 2(e)].

Vortex types and sizes
We start our analysis of vortex structures by reconsidering the simulation snapshots in Fig. 2(h) and (e) with an alternative coloring scheme illustrating the motion of each particle, as shown in Fig. 7. Specifically, we characterize the motion of particle i relative to its orientationê e e i =μ µ µ i by a parameter g i (t), defined as where the (average) velocity v v v i,s (t) for particle i in a short time interval between t and t + ∆t s is given by The arrow corresponding to particle i in Fig. 7 is colored in red if g i (t) ≤ 0, meaning that the particle moves "backward" (against its orientation). In contrast, the arrow is colored in blue if the particle stops or moves "forward" (along its orientation). When calculating v v v i,s , we take into account the fact that an active rotating particle with zero motility needs more time to move a distance equal to its diameter than a highly motile circle swimmer. Therefore, we make the choice ∆t s = τ for the case of v * 0 = 0 and ∆t s = 0.01τ for v * 0 ≥ 40. Representative simulation snapshots at ω * 0 = 3 and three values of the motility. Each of the disk-shaped particles is represented by an arrow indicating the particle orientationê e e i . For better visual quality, all arrows are scaled by a factor of 3 and are formatted with red (blue) color if the particle moves against (along) its own orientation [see Eq. (12)]. The blue arrows display Type I vortices (a-c), and the red arrows exhibit Type II vortices (a).
All color-coded snapshots in Fig. 7 refer to ω * 0 = 3. At zero motility (v * 0 = 0), the particles have no tendency to move along their orientation. Therefore, the numbers of particles moving forward (blue) and backward (red) are on average equal, consistent with Fig. 7(a). A surprising result, however, is that forward-and backward-moving particles aggregate with particles displaying the same kind of motion, yielding vortices. In each instant of time, the chirality of the vortex structure formed by the forward-moving particles is counterclockwise, which we refer to as Type I vortices. In contrast, the chiral structure of the backward-moving particles is clockwise (Type II vortex). We note, however, that since the intrinsic active rotation drives each of the particles to rotate counterclockwise (see Sec. 2.2), the direction of overall rotation for both, Type I and Type II vortices, is counterclockwise. The detailed discussion about the physical origin of the emergent Type I and Type II vortices is postponed to Sec. 3.4.2.
Once the motility becomes non-zero, the particles tend to selfpropel along their orientation (instantaneous) orientation. We thus expect that values of g i > 0 become more and more relevant. Indeed, as seen from Fig. 7(b), we only observe forward-moving swimmers which are represented by blue arrows as v * 0 increases from zero to 40. Comparing Fig. 7(a) and (b), we further see that the vortex size significantly increases with v * 0 . At an even higher motility [v * 0 = 80, see Fig. 7(c)], the diskshaped vortex structures seen at lower motilities are slightly distorted and have a hole in the vortex center. To understand this, we recall that an isolated circle swimmer moves along a circular path with a radius R = v 0 /ω 0 . 19 Thus, R increases with v 0 . This ideal motion, however, is disturbed by the presence of other swimmers, which makes it increasingly difficult to stay in the center of a giant vortex when R becomes larger. This eventually leads to a hole in the vortex center. Further information on the vortex structure and size is provided by the orientational correlation function C e e e (r), whose detailed analysis is presented in Appendix A. Here, we focus on the vortex diameter d e e e , which is obtained from the minimum of C e e e (r) (see Sec. 2.4.3). Figure 8 shows d e e e as a function of the angular speed ω * 0 for various motilities v * 0 . Upon an increase of v * 0 , the range of ω * 0 where vortices appear is narrowed. For all explored motilities, d e e e generally decreases with the angular speed ω * 0 , and the curves are shifted toward a larger d e e e as v * 0 increases. Specifically, we observe a power law decay with an exponent ν ≈ −2 for v * 0 20. This exponent is different from that observed in studies of the rotating Vicsek model (and variants) where ν = −1. 34,35 Finally, at a fixed motility such as v * 0 = 40 and upon an increase of the angular speed from ω * 0 = 0 to ω * 0 = 3, the system undergoes a transition from a state with short linear chains (i.e., string fluids) to a state with giant vortices, as shown in Fig. 2(d) and 2(e).
Here we identify the transition by a drastic increase of φ c (see Fig. 4), and by the appearance of negative orientational correlations [see Fig. 12(b)]. Clearly, it would be very interesting to see whether this non-equilibrium transition is/resembles a first-order or continuous transition. However, such a study would require more extensive simulations and a full analysis of the order parameters, which is out of scope for the present paper. Fig. 9 (a, c) Sketches of two initial chain configurations composed of N t = 14 dipolar circle swimmers. The black arrows indicate the particle orientation at t/τ = 0. (b, d) The inverse of maximum gyration radius, R −1 g,M , between 30 ≤ t/τ ≤ 40 as a function of the motility v * 0 for initial configuration (a) (black circles) and (c) (red squares) at the angular speed ω * 0 = 3. Sketches show respective configurations characterized by R −1 g,M , where the red (blue) arrow indicates that the particle's velocity points against (along) its orientation [see Eq. (12) with ∆t s = 0.01τ] and the dashed arrows indicate the particle motion. The rectangles with curved angles relate data with corresponding sketches.

Ring argument
To understand the physical origin of the emergent vortices in the present system, we first recall that passive dipolar particles, at the large coupling strength considered here, self-assemble into chains, see Fig. 2(g). Moreover, particles in neighboring chains arranged side-by-side tend to point in anti-parallel directions. Inspired by this behavior we performed, at an intermediate angular speed ω * 0 = 3, test simulations of a simplified double-chain structure to investigate the impact of the motility v * 0 . As an initial configuration, N t dipolar particles are arranged into two anti-parallel chains lying side by side, see Fig. 9(a) and (c). Within each chain, the dipole moments are oriented head-to-tail. As illustrated in Fig. 9(a) and (c), we consider two different arrangements of these double-chain configurations, which are the mirror image of each other, meaning that they have opposite "chirality." The centerto-center distance between the two chains is chosen to be 2σ . Each of the chains is composed of N t /2 = 7 dipolar particles, such that the resulting chain length corresponds to the vortex diameter d e e e /σ ≈ 7 observed at v * 0 = 0 and ω * 0 = 3 (see Fig. 8). For simplicity, the thermal fluctuations are neglected.
Upon starting the simulations with ω * 0 = 3 and v * 0 = 0 we observe, for both initial configurations, that the swimmers rearrange themselves into a stable ring structure. The difference, however, is that the particles in the ring move forward (like in a Type I vortex) when starting from configuration (a), while the particles move backward when starting from the other initial configuration. This difference is due to the competition between the ring chirality favored by the initial configuration and the counterclockwise rotation supported by the active rotation. To further characterize the ring structure, we consider the gyration radius, given by The radius R g (t) is a constant once the chains form a stable ring, while R g (t) varies with time t if the chains remain separate and do not self-assemble into a ring even after a long time. With this in mind, we plot in Fig. 9(b) and (d) the inverse of the maximum gyration radius R −1 g,M /σ −1 within a time interval 30 ≤ t/τ ≤ 40 as a function of the motility v * 0 for the initial configuration sketched in Fig. 9(a) and (c), respectively. The choice of this time interval is based on the fact that no transient structure is observed at t/τ > 20. The quantity R −1 g,M /σ −1 is about 0.4 for a stable ring structure and reaches a value smaller than 0.1 if two chains do not form a ring. As v * 0 1, both Type I and Type II rings are stable, which conforms with the emergence of Type I and Type II vortices at v * 0 = 0 and ω * 0 = 3 [see Fig. 2(h) and Fig. 7(a)]. We also see that two chains no longer form a ring once v * 0 is greater than a "critical" motility v * 0,c . In this latter case, each chain swims along a circular path and does not collide with the other chain within the simulation time, see the sketches on the right side of Fig. 9(b) and (d). Importantly, v * 0,c for the Type II ring is much smaller (v * 0,c ≈ 3) than that for the Type I ring (v * 0,c ≈ 30). This explains why the system is dominated by the Type I vortices at intermediate motilities [see Fig. 2(e) and 7(b)]. Further increasing v * 0 leads to the breaking of both Type I and Type II rings [see Fig. 2(b)], suggesting that vortices will eventually vanish for large v * 0 . Indeed, for v * 0 80 the system can display flocking behavior, which we will further address in the following Sec. 3.5.
To check whether the vortex emergence is sensitive to the system size of the simulations, we performed various test simulations at v * 0 = 40 and ω * 0 = 4 with the particle number ranging from N = 100 to N = 2500. We observed that giant vortex structures already emerge at N ≈ 400 and persist up to N = 2500 (results not shown here).
The vortex formation in the present work is of fundamental difference from the emergence of the clockwise vortices in systems of simple (non-dipolar, disk-shaped) circle swimmers, 28 and the macroscopic droplets reported in studies of the rotating Vicsek model 34 and its variant. 35 For simple circle swimmers, the clockwise vortices do not appear at a density lower than the critical density of MIPS, Φ crit ≈ 0.3. 28 Further, their formation relies on the steric collision of particles in the dilute region with the boundary of giant clusters. In contrast, for dipolar circle swimmers, the vortex patterns appear already at densities lower than Φ crit . For circle swimmers with polar alignment, such as the rotating Vicsek model 34 and its variant, 35 macroscopic droplets emerge at small angular speed ω * 0 . Inside the droplets, particles align themselves with their neighbors and are directed along a single direction. This direction rotates in response to the active rotation exerted on each of the particles. These macroscopic droplets are significantly different from the giant vortices observed in dipolar circle swimmers, suggesting that the type of alignment interactions is crucial for the pattern formation of circle swimmers.

Orientational ordering
In this last section, we switch our focus onto motility-induced orientational ordering appearing at large v * 0 (see the top part of the state diagram in Fig. 1). We note that in the present model, the orientational order of dipole moments µ µ µ i implies coherent motion, i.e., flocking. The orientational order can be characterized by the global polar order parameter φ e e e , defined in Eq. (11).  Table 1]. Further, Fig. 10 reveals that the system reaches bistable states at v * 0 = 100 and ω * 0 = 5 − 6, consistent with earlier discussion of Fig. 4. Taking a closer look at Fig. 2(a)-(c), we find that the size of the ordered ("flocking") structures formed at large v * 0 significantly depends on ω * 0 . To further characterize this behavior we consider the (weighted) distribution of the cluster size, nP (n), and their dependence on the overall particle number N. Results for v * 0 = 100 are shown in Fig. 11. As can bee seen in Fig. 11(a), for zero rotation (ω * 0 = 0) the distributions decrease monotonically with the cluster size n and collapse onto one curve for N 900. This indicates that these are only small clusters whose size does not scale with the particle number N. According to Table 1, we classify such a situation as a "micro-flocking" state. In contrast, the weighted distribution function for fast rotation decays at small cluster sizes n, but exhibits a broad peak at large n [see Fig. 11(b)]. This peak corresponds to the large swarms shown in Fig. 2(c). On increas-ing the particle numbers N, the decay of nP (n) at small n becomes faster, while the peak at large n is shifted toward a larger n. In other words, the size of the formed structure scales with the particle number N, indicating the emergence of macroscopic swarming patterns. Based on Table 1, we identify this as a "macro-flocking" state. Fig. 11 Weighted distribution of cluster size at v * 0 = 100 and ω * 0 = 0 (a) and ω * 0 = 10 (b) for various particle numbers N.
The mechanism underlying the flocking behavior of dipolar active systems is quite complex already in the absence of circle swimming (that is, at ω * 0 = 0, see ref. 46): Short linear chains formed by dipolar active particles with head-tail orientation tend to align their velocities upon collisions (see Fig. 8 of ref. 46). In ref. 46, it is found that macro-flocking appears once the density increases from Φ = 0.23 to Φ = 0.58. This suggests that the transition from micro-to macro-flocking may be attributed to a densityinduced enhancement of particle collisions. For the present system of chiral dipolar active particles, a similar effect may take place: We suspect that the circle swimming of each particle leads again to an enhancement of collisions, similar to an increase of density. This might explain why the present, chiral system exhibits macro-flocking already at smaller densities than the corresponding non-chiral system (ω * 0 = 0). We note that our observation concerning the size of the flock-ing pattern contrasts the behavior previously observed in systems of chiral active particles with polar interaction. There, the size of flocking patterns is inversely proportional to the angular speed. 34,35 This suggests that different types of alignment interactions between active particles may significantly change the fundamental properties of flocking behavior.

Conclusions
In this work, we performed extensive Brownian dynamics (BD) simulations to investigate the pattern formation of dipolar circle swimmers dispersed on a monolayer. To this end, we explored a wide range of active angular speeds and motilities at a fixed density below that related to MIPS in simple systems and large dipolar coupling strength.
At a sufficiently large angular speed (ω * 0 1) and zero motility (v * 0 = 0), the system undergoes a transition from a state with percolated networks into a state with Type I and Type II vortices. Upon an increase of v * 0 from zero, the Type II vortices vanish, and the system displays phase separation with the dense domain characterized by giant Type I vortices. Based on test simulations of two anti-parallel chains composed of strongly coupled dipolar particles, we proposed a "ring" argument to unveil the underlying mechanism of the vortex formation and the disappearance of Type II vortices. In contrast with our model, the vortex structures are not observed in systems of circle swimmers with ferromagnetic Heisenberg-like interactions. 34,35 Instead, these systems display macroscopic polar droplets, in which swimmers move coherently along a certain direction and the direction rotates in accordance with the active rotation exerting on each circle swimmer. To further unravel the differences, we performed Brownian dynamics simulations of circle swimmers with ferromagnetic interactions decaying with a Yukawa potential (as a function of separation), 87,88 and, thus, decay with the particle distance. Our preliminary results (not shown here) do not show any vortex structures, either.
A further increase of the motility (v * 0 80) leads to two distinct flocking states, which can be distinguished by the cluster size distribution. Consistent with the behavior of dipolar active Brownian particles, 46 dipolar circle swimmers at zero and slow rotation display a micro-flocking state (ω * 0 3). In contrast, dipolar swimmers with fast rotation (ω * 0 7) exhibit a macro-flocking state. This is again distinctly different from circle swimmers with shortrange ferromagnetic interactions, in which the size of polar clusters decreases upon an increase of the active angular speed. 34,35 Hence, the type of anisotropic interactions that align particle velocities plays a vital role in determining the fundamental selforganization process for systems of circle swimmers.
For completeness, we also performed simulations for circle swimmers with truncated dipole-dipole interactions (results not shown here). Here, the Ewald summation was not employed, and the dipole-dipole interactions were truncated at r = 3σ . In this case, the flocking transition is shifted to a much larger motility (v * 0 ≈ 200). More importantly, we did not observe any vortex patterns throughout the explored parameter regime. This suggests that not only the angle dependency, but also the long-range character of dipole-dipole interactions is crucial for the collective behavior of dipolar active systems and, thus, should be not be neglected.
The model studied in the present paper does not account for hydrodynamic interactions between particles. However, hydrodynamic interactions can have a profound impact on the dynamical self-assembly of active colloidal systems, such as the suppression of MIPS [89][90][91][92] and the emergence of global polar ordering. [93][94][95][96][97] In systems of rotating particles, it is found that hydrodynamic interactions can induce cluster rotation, such that the overall cluster and the individual particles rotate in the same direction. 98 In the present work, we have seen that the giant Type I vortices and circle swimmers rotate in the same direction. Hence, we expect that hydrodynamics can further promote the formation of giant Type I vortices. Nevertheless, the detailed influence of hydrodynamics on dipolar circle swimmers remains to be unveiled by future works.
Furthermore, it is well established that mixtures of active and passive colloidal particles display fascinating collective behavior that is distinctly different from that of the corresponding onecomponent systems. 64,76,[99][100][101] Therefore, one future direction could be to consider mixtures of circle swimmers and passive dipolar particles, whose self-assembly process can be further controlled by the proportion of species.
Finally, the majority of the research works on systems of active particles focus on their collective behavior in two dimensions (2D), whereas the real-world suspensions of self-propelled colloids are often in three dimensions (3D). The dimensionality may play an important role in the motility-induced phenomena, such as the critical motility for MIPS 102 and the critical exponents for flocking transition. 86 In the present 2D model of dipolar circle swimmers, we assume that the particles are confined to a flat surface, and the rotation axis for each particle is restricted to be normal to the surface. It will be interesting to discover the similarities and differences in our model system when moving from 2D to 3D.

Conflicts of interest
There are no conflicts to declare.

A Orientational correlation function
To characterize the vortex structure, we calculate the orientational correlation function C e e e (r) defined in Eq. (9). Results are shown in Fig. 12. For passive dipolar particles (v * 0 = 0 and ω * 0 = 0), C e e e (r) decays with r and reaches zero at r/σ ≈ 10 [see Fig. 12(a)]. Upon an increase of ω * 0 from zero to 1, C e e e (r) decays faster than the case where ω * 0 = 0 and displays a negative correlation with a minimum C e e e (r) ≈ −0.15 at r/σ ≈ 10. As have been discussed in Sec. Comparing the results for v * 0 = 0 and v * 0 = 40, the vortex structures emerge at ω * 0 = 1 for v * 0 = 0, whereas for v * 0 = 40 the vortices appears at a larger angular speed ω * 0 = 3 [see Fig. 12(b)]. In particular, the vortex size d e e e becomes much larger and is close to the half of the side length of the simulation box (L/2 ≈ 33.69σ ), suggesting that there might be finite-size effect for d e e e . In other words, d e e e may scale with the total number of particles N, which requires further investigation. Nevertheless, since the long-range character of dipole-dipole interactions requires expensive computational resources, the present simulations are limited to the order of 10 3 particles. A further increase of ω * 0 from 3 to 10 causes the giant vortices to break into smaller vortices.
At high motilities such as v * 0 = 80, C e e e (r) drastically decays at short distances r/σ ≤ 3 and gradually reaches a value close to 0.3 at r = L/2 for zero, small and large angular speeds (ω * 0 = 0, 1 and 10). The non-vanishing, positive correlation function indicates the emergence of global orientational order, which is discussed in detail in Sec. 3.5. At the intermediate angular speeds (ω = 3 and 5), the minimum of C e e e (r) appear at r L/2.