Swirling transition with social interactions: Analyzed by a sixth-order Landau-type model

Swirling groups of animals or bacteria are a common phenomenon in nature. It is thought that this collective organization occurs in the vicinity of a continuous transition between dynamic states to ensure robust group cohesion while allowing for high sensitivity to outside stimuli like predators. Here, we present Brownian dynamics simulations of active particles with social interactions which can form stable swirls. We observe a transition between swarming and swirling states and analyze these using a sixth-order Landau-type model. Our results suggest that the transition is weakly discontinuous. However, by lowering the rotational diffusion coefficient, it becomes continuous.

based on simple interactions between neighbouring entities within larger groups, named social interactions [31][32][33][34][35][36]. These interactions do not represent external forces, but rather describe the way individuals react to their respective environment. Agent-based models with social interaction rules have successfully been used to describe collective decision making in animal groups, for example, in schooling fish [37][38][39], in human pedestrians [40], and in human group coordination [41].
Swirling phenomena, in which individual entities create collectively ordered rotational motion, are a special case of this wider class of collective organization. They can be found at different scales in nature, from the behaviour of microscopic bacteria [42], up to ants [43], and larger animals like schooling fish [44]. It has been proposed that such systems exist in the vicinity of a critical point belonging to a phase transition, to ensure robust group cohesion while allowing for high sensitivity to outside stimuli like predators [31,45].
Active particles (APs) are a versatile tool to model and explore the collective dynamics of energy consuming, self-propelled entities with tailored interactions, both in experiment and theory [4]. In this letter we perform simulations of APs with social interactions that show a transition between a swarming and an ordered swirling state. We directly address the experimental work by Bäuerle 47003-p1 et al. [46], in which a set of social interaction rules was devised and applied to propelled silica particles, in order to replicate the collective swirling dynamics found in nature. The observed transition between the swarming and swirling states is a symmetry-breaking transition described by an order parameter O R , where the sign of O R corresponds to counter-or clockwise rotating states. According to the experiments [46,47] the transition displays critical dynamics, which can be accessed through a Landau-type description with an effective free energy f = aO 2 R + bO 4 R . Here, we perform simulations following ref. [46] and reanalyze the probability distribution for the fluctuating swirling order parameter which we write similar to a Boltzmann distribution and where Δ is the swirling control parameter to be defined below. Our analysis suggests that the effective free energy appearing in P (O R , Δ) has to be described by a Landautype expansion up to sixth order in O R , Note that odd terms in O R cannot appear since left-and right-turning swirls are equally probable. The expansion coefficients a, b, and c depend on Δ. Now, the effective free energy of eq. (2) allows a richer phenomenology of transitions between dynamic states compared to the fourth-order model, where only a continuous symmetry-breaking transition occurs at a = 0 from the swarming (a > 0) to the swirling state (a < 0) with b > 0. In contrast and as sketched in fig. 1, the following scenario occurs for the sixth-order Landau-type model with constant c > 0 [48]. For positive b one still observes the continuous transition at a = 0 into the swirling state. However, for a negative fourth-order coefficient b the line at a = 0 goes over to a line of discontinuous or first-order transitions that bends away from the negative b axis. Both lines meet at the tricritical point (a = 0, b = 0).
In the following we show results from Brownian dynamics simulations for the symmetry-breaking swirling transition. They seem to suggest that the experiments in refs. [46,47] might display a weak discontinuous transition, which can be described by a Landau-type model of sixth order. Furthermore, by lowering the rotational diffusion coefficient of the active particles, the system moves around the tricritical point in the upper half plane, when varying Δ, so that the transition becomes continuous (see fig. 1 for a schematic).
Methods. -Active particles (APs) are able to selfpropel and reorient themselves by absorbing energy from their surroundings [4,49]. Their motion is influenced by a combination of random fluctuations and active swimming called active Brownian motion. They represent a minimal model for active individuals and can thus, in combination with appropriate interaction rules, be used to investigate the principles underlying the collective behaviour observed in biological systems.
Zonal models with social interaction rules have repeatedly been used to describe the collective behaviour of living organisms [32,33,[37][38][39]. In our simulations, the two-dimensional model applied in the experiments conducted by Bäuerle et al. [46] is used. The rules defined therein are based on the avoidance and alignment between neighbouring entities. The motion of the APs is restricted to a plane and depends on the presence of other particles in circular concentric zones of different radii. The zones are depicted in fig. 2(a) and will be referred to as the zone of repulsion (ZOR) and the zones of detection of orientations (ZOO) and positions (ZOP). Note that this zonal model can be extended by introducing a field of vision for all particles, as opposed to using their entire surroundings as we do here [46]. The APs are characterized by their positions r = rr relative to the mean position of all particles, where r is the distance andr the unit directional vector. Furthermore, they have propulsion velocities v = vû with speed v and unit orientation vectorû. Now, we describe the social interaction rules that determine a target orientation towards which the APs reorient while experiencing a torque. If there are no APs within the ZOR of an AP i, its target orientation is determined by calculating the average relative position P i of all APs within the ZOP, as well as the average orientation û i of neighbouring APs within the ZOO, according to 47003-p2 Here, I ZOP,ZOO denote the respective sets of particles in the two zones around particle i and n ZOP,ZOO are their numbers ( fig. 2(b)). To determine the target orientation, P i is rotated by an angle Δ either to the left (+) or to the right (−), which generates two unit vectorsd ± i . The one closer to the average orientation û i is then taken as the target orientation for particle i. As will be shown later, the transition from swarming to swirling of the APs strongly depends on the magnitude of the deviation angle Δ. If other APs are present within the ZOR of an AP, another rule is applied. To avoid collisions, the target orientation is set along the reversed mean relative position vector of those APs within the ZOR as depicted in fig. 2(c).
The swirling state of a group with N APs is quantified by introducing the rotational or swirling order parameter where the unit vector e z is normal to the plane in which the particle motion takes place. It ranges from −1 to 1, where O R = 0 indicates swarming. We choose our parameters as in the experiments [46], where the APs have a diameter of σ = 6.3 μm. The radius of ZOP, in which the positions are detected, is set to R P = ∞ such that P i points to the centre of the entire group of APs, while the orientations are detected within the smaller ZOO with a radius of R O = 25 μm, approximately four times the AP diameter. For the ZOR the radius is set to R R = 8 μm. This results in a surface-tosurface distance of two particles placed at the ZOR center and rim, which amounts to approximately a quarter of the particle diameter. The APs self-propel with a velocity v = 0.5 μm s −1 and rotate by applying a maximal torque of Γ max = 25 k B T as measured in the experiments. They are subject to translational and rotational diffusion, with the diffusion constants D T = 1.4 · 10 −2 μm 2 s −1 and D R = 2.8 · 10 −3 s −1 as reported in the experiments [46]. Their maximal reorientation rate of ω max = Γmax kB T D R ≈ 4 • s −1 is determined using the rotational mobility from the Einstein relation, [46]. Due to the idealized particle interactions used in the simulations, it is necessary to enhance the rotational diffusion constant D R , which controls the systems rotational noise, in order to obtain a good match with the experiments, as discussed later in fig. 3. The best agreement is achieved , which we will use in our simulations. To keep the measured maximal reorientation rate ω max = Γmax kB T D R constant, the maximal torque is scaled down according toΓ max = Γ exp max D exp R /D 0 R . We perform Brownian dynamics simulations of the APs in two dimensions, where r(t) andû(t) = (cos θ, sin θ) denote their position and orientation vectors, respectively. To simulate the dynamics of an AP, one first determines the target orientation and with it the rotational velocity . Then the dynamics is described by two coupled Langevin equationṡ Here, the random velocities ξ T and ξ R represent thermal white noise with zero mean and normally distributed components with respective variances The random velocities of different particles are uncorrelated.
To solve eqs. (6) and (7), we apply a simple Euler integration, where W T i and W R are Gaussian random numbers with zero mean and unit variance. A time step size of Δt = 0.2 s corresponding to the image acquisition rate of 5 Hz in the experiments of ref. [46] is used, since it is small enough to ensure proper integration of the dynamic equations (10) and (11). The initial positions and orientations are randomly distributed. They are drawn from normal distributions with zero mean and unit variance for the orientations and a variance of 5 σ for the positions. Due to the finite 47003-p3 size of the simulated APs, collisions must be handled. This is done after each time step by translating overlapping particles symmetrically away from each other along the connecting line. Several iterations are used until no more overlaps are present. As in the experiments we simulate N = 50 APs.
Results. -We first discuss the phenomenology of our simulations. At large enough deviation angles Δ, we indeed observe the formation of stable swirls as shown in fig. 3(a). The trajectories of all 50 particles are depicted during a period of 100 s. The swirl does not arise from a perfect concentric rotation of the individual particles. Instead, due to translational/rotational diffusion and active motion of the interacting particles, their paths show frequent directional changes, which on larger time scales lead to loose adherence to the group rotation (see the supplementary video video1.mp4). If one follows the trajectories of neighbouring APs, frequent changes in their relative positions and even reversals in the swimming directions of single particles are visible, as shown in the supplementary video video1.mp4. These changes in the swimming direction of single APs can lead to changes in the rotational direction of the entire swirl, as the plots in fig. 3(b) depict. During these reversals, the group generally enters a metastable swarming or flocking state with O R ≈ 0 for a short period of time, after which rotational order establishes again. This occurs either in the opposite (top and the supplementary video video2.mp4) or same direction (bottom) as the previous stable swirl.
To quantify the strength of swirling motion, we use the swirling order parameter O R . However, averaging it over time gives zero, since O R and −O R are equally probable. For this reason we consider the mean of the absolute value, |O R | , which shows a strong dependence on the deviation angle Δ, as seen in fig. 3(c). At Δ = 0, the APs aim to swim in direction of the group centre, which leads to little rotational order but rather induces swarming. As Δ increases, the APs start to circle around their common centre, as indicated by the increasing |O R | . Utimately, it reaches a maximum at Δ max = 78 • ± 1 • , below the deviation angle Δ = 90 • , which would correspond to ideal circular motion of the APs. The reason is that the group of APs increases in size for larger deviation angles, as fig. 3(c) strikingly demonstrates, up to a point at which the cohesion of the swirl is lost since fewer and fewer APs are in the ZOO. This results in rapidly decreasing order above Δ max . To quantify the spatial extent of a swirl, we use the mean distance r of all particles from the group centre. Note that our simulations perfectly match the experimental results as fig. 3(c) demonstrates, which enables us to utilize them for further analysis of the swirling active particles.
To investigate the steady swirling motion in more detail, we plot the probability distributions P (O R , Δ) of the swirling order parameter O R for different deviation angles Δ in fig. 4. For small values of Δ, the distributions show  [46,47] suggest choosing an effective Landau-type free energy f of fourth order in O R for the symmetry-breaking transition. However, a careful inspection of the distributions in fig. 4 always shows a weak maximum at O R = 0, even after the two peaks at O R = 0 have developed. The weak maximum corresponds to the metastable swarming state, which the APs assume during reversals of the swirling direction. These maxima do not appear in the experimental results presented in [46]. We attribute this to the fact that, as mentioned by the authors, measurements in which the APs enter this metastable state had to be discarded due to the limited field of view in the experiment. In our simulations we can easily include them. However, the corresponding peaks cannot be described by a fourth-order model of f in O R . For this reason we fit our data with a sixth-order expansion of the effective free energy introduced in eq. (2) and obtain the probability distribution Here a, b and c are used as free fitting parameters in a least-square fit. They are plotted vs. Δ in fig. 4. The fits of P (O R , Δ) show good agreement with the data. In particular, the peaks around O R = 0, which also occur at large Δ and then correspond to the metastable swarming state, are taken into account. Note that the dependence of the coefficients a and b on Δ agrees with the simplest possible forms. Since the maximum of P (O R , Δ) at O R = 0 always requires a > 0 or a positive curvature of the effective free energy f , the coefficient a cannot be linear in Δ as in a conventional fourth-order Landau model. Accordingly, the next higherorder approximation is a parabolic shape, which is nicely demonstrated by the quadratic fit curve in fig. 4, right. Furthermore, the coefficient b has to change sign, which is necessary so that the two additional peaks at O R = 0 can develop. Here, the simplest approximation is a linear dependence in Δ, which again is in nice agreement with our findings as the fit in fig. 4, right shows. Finally, the sixth-order term in the effective free-energy expansion is necessary to maintain stability. Interestingly, the corresponding coefficient c is independent of Δ for Δ ≥ 20 • , where it is needed for stability.
We use the effective free energy f (O R , Δ) of eq. (2) in the fluctuating order parameter to further clarify the dynamic swirling transition. Following the explanations in the introduction, we plot the state diagram in the a-b plane using the constant c as determined above in fig. 6. For b > 0 the free energy predicts a continuous transition at a = 0, while for b < 0 the solid line at a c = b 2 /(4c) refers to a discontinuous transition. The blue dots belong to all the pairs of a, b values determined in fig. 4 for different Δ, while the dashed line uses the parabolic and linear fits for a, b, respectively. Starting from the the upper right at Δ = 0, the dots move to the left with increasing Δ, assume b < 0, and then cross the first-order transition line. This clearly shows that the swarming-swirling transition is discontinuous. The intersection between the system path and the transition line yields a transition point Δ * ≈ 27 • .
As a consequence of the effective free energy in eq. (2), the most probable order parameter bifurcates from O m R = 0 for Δ < Δ * into both coexisting stable swirling states with non-vanishing O m R for Δ > Δ * . The values for O m R in the ordered phase are determined by the two symmetric global minima of the free energy in eq. (2) and are located at In accordance with our previous discussion, the bifurcation is discontinuous, as illustrated in fig. 5. While the dots correspond to the values of a, b, and c determined for each Δ, the solid line uses the respective fit functions from fig. 4, right. The very good agreement shows the high quality of the fits. In a next step, we investigate whether varying the system parameters will influence the observed swirling transition. We repeat the above analysis for three reduced values of the rotational diffusion coefficient D R compared to our default value D 0 R , since we noticed that varying this parameter leads to crucial changes in the observed transition. The resulting system paths are also presented in fig. 6. One notices that when reducing the rotational diffusion coefficient, the system path in the a-b plane moves 47003-p5  to the left. In particular, at 0.7 D 0 R (red dots) the swirling transition is localized at the tricritical point, while for smaller coefficients the system paths intersect the secondorder branch of the transition line at a = 0, b > 0 and the dynamic state transition becomes continuous. Thus, in the region with a < 0 and b > 0 the effective free energy f has a local maximum at O R = 0 and the corresponding probability distributions P (O R ) a local minimum instead of a peak. This is the range in which a fourth-order Landautype model can successfully describe the system. However, as soon as b < 0, the sixth-order term in the free-energy expansion is necessary to maintain stability. For a < 0 the local minimum at O R = 0 persists until a becomes positive again. Then a peak of P (O R ) at O R = 0 occurs, corresponding to a metastable swarming state. Note, the respective transition points move to decreasing deviation angles when D R is lowered. For the smallest rotational diffusion coefficient 0.3 D 0 R it is located at Δ * < 15 • . Finally, we note that for decreasing D R the parameter c depends more and more on Δ and, as opposed to the default value D 0 R , is no longer constant.
Conclusion. -The simulations of active particles with social interaction rules presented in this letter can successfully reproduce the transition between unordered swarms and stable swirls observed in experiments for a group of 50 particles [46]. The two dynamic collective states are realized when varying the deviation angle Δ, which controls the swimming direction chosen by individual active particles according to the social-interaction rules governing their motion. For sufficiently large Δ, the active particle groups create ordered swirls, while no such behaviour can be observed when the deviation angle is too small. The symmetry-breaking transition between those states can be described with a Landau-type model of sixth order. The model depends on the swirling order parameter O R and the deviation angle Δ serves as a control parameter. This approach suggests a weak discontinuous transition for the swarming and swirling states occurring at Δ * ≈ 27 • . Furthermore, we realized that the rotational diffusion coefficient D R is the relevant parameter for influencing the observed swirling transition. Interestingly, when lowering D R , the transition becomes continuous.
In our simulations we closely follow the experimental work of ref. [46] and employ the reported parameters. Nevertheless, the analysis of ref.
[47] on the same system reports a critical behavior instead of the weakly discontinuous transition, which we observe. From the experimental side the difference might be due to the fact that the experimental setup does not fully take into account a possible metastable swarming state occuring at large Δ [46]. The existence of such a state at O R = 0 can only be described by a Landau-type model of sixth instead of fourth order. From the theoretical side, we note that our simulations nicely capture the experimentally observed dependence of |O R | on Δ, as fig. 3(c) shows. Still, simulations need to make simplifying assumptions so that experiments are not fully replicated. In our case this could be, for example, the implementation of particle collisions or the fact that hydrodynamic interactions between the colloids are completely omitted. To resolve the issue of a possibly weakly discontinuous transition, it might be worthwhile to further investigate the existence of a metastable swarming state in experiments.
There are three further directions we plan to pursue with this system. First, social interactions are not instantaneous. Therefore, we will implement a time delay with which the agent or active particle adjusts its new orientation in response to the surrounding particles. It will be interesting to explore how the observed dynamics change, depending on the magnitude of this time delay. Second, we will introduce variations in the deviation angle Δ for single agents, or some distribution in Δ, and monitor the influence on the swirling state. This approach will model the fact that not all agents, such as fish, behave completely identically or that some strongly deviate in their behavior from the group. Finally, variations in the detection radii and the viewing angle, which is introduced to define a field of vision for all particles [46], can be used to model the different perceptions that agents might have of their environment. * * * We thank Johan Dettmar for providing an initial version of the simulation code, on which this work is based, and Clemens Bechinger for helpful discussions. We also acknowledge funding from the Berlin University Alliance under grant No. 824 BUA-NUS 4.

Data availability statement:
The data that support the findings of this study are available upon reasonable request from the authors.