The chaotic milling behaviors of interacting swarms after collision

We consider the problem of characterizing the dynamics of interacting swarms after they collide and form a stationary center of mass. Modeling efforts have shown that the collision of near head-on interacting swarms can produce a variety of post-collision dynamics including coherent milling, coherent flocking, and scattering behaviors. In particular, recent analysis of the transient dynamics of two colliding swarms has revealed the existence of a critical transition whereby the collision results in a combined milling state about a stationary center of mass. In the present work we show that the collision dynamics of two swarms that form a milling state transitions from periodic to chaotic motion as a function of the repulsive force strength and its length scale. We used two existing methods as well as one new technique: Karhunen-Loeve decomposition to show the effective modal dimension chaos lives in, the 0-1 test to identify chaos, and then Constrained Correlation Embedding to show how each swarm is embedded in the other when both swarms combine to form a single milling state after collision. We expect our analysis to impact new swarm experiments which examine the interaction of multiple swarms.


I. INTRODUCTION
Improved data gathering and analysis techniques have enabled researchers to collect data and analyze the motions of individual agents in biological flocks, and formulate more accurate, empirical models for collective motion strategies of flocking species such as birds, fish and insects, 1,2 .The derived mechanics has resulted in the translation of swarm theory to communicating robotic systems.Swarms of collaborating robots have been proposed for conducting tasks including search and rescue, density control, and mapping 3,4 .
From general models of swarm dynamics, one observes three basic swarming states or modes: flocking, in which a a) SK is a National Research Council Postdoctoral Research Associate.center of mass moves in a straight line with a constant velocity; milling, where the agents rotate coherently around a stationary center of mass; and rotating, where the center of mass itself rotates about a fixed point in space [5][6][7] .Much is known about the behaviors and stability of single swarms with physics-informed, nonlinear interactions [8][9][10][11] .They are able to converge to organized, coherent behaviors in spite of complicating factors such as communication delay, localized number of neighbors each agent is able to interact with, a consensus force, heterogeneity in agent dynamics, and environmental noise [12][13][14][15] .However, in many applications, beyond the dynamics of single swarms, multiple swarms interact and produce even more complex spatio-temporal behaviors.
New studies have begun to address swarm-on-swarm dynamics, and in particular the scattering of two large, colliding swarms with nonlinear interactions.Recent numerical studies have shown that when two flocking swarms collide, the resulting dynamics typically appear as a merging of the swarms into a single flock, milling as one uniform swarm, or scattering into separate composite flocks moving in different directions [16][17][18] .However, more detailed analytical understanding of how and when these behaviors occur is necessary, especially when designing robotic swarm experiments such as swarm herding and capture [19][20][21] and controlling their outcomes.
Building on early numerical insights, it has been shown that one can predict the critical swarm-on-swarm interaction coupling, below which two colliding swarms merely scatter, as a function of a physical swarm parameter 7 .In particular, analysis of the transient dynamics of two colliding swarms, seen in Fig. 1, has revealed the existence of a critical transition whereby the collision results in a combined milling state, where one swarm is embedded in the other.However, what is missing is a quantitative description of the post-collision asymptotic dynamics.
In the current work, our starting point is Fig. 2. We observe that after the colliding swarms have settled down to a milling behaviors, the dynamics change from periodic to quasi-periodic, and to chaotic motions for different values of repulsion strength.We examine the dimensionality of the combined mill by computing its Karhunen-Loeve dimension.It is a measure of the variance of the entire post-collision attractor, and computes how many expansion modes it takes to capture the dynamics.Next, we aim to numerically confirm that the two swarms near collisions displayed complex chaotic behavior with two methods.First, we determine whether the collision dynamics is chaotic, by using the existing binary test for chaos.Finally, we defined the constrained correlation embedding to quantify the dynamical complexity of the combined milling state.

II. SWARM MODEL
To make progress, we considered a class of well-known deterministic swarming models for both discrete 8,22,23 and continuous 9,23 systems consisting of mobile agents moving under the influence of self-propulsion, nonlinear damping, and pairwise interaction forces.In the absence of interactions, each swarmer tends to a fixed speed, which balances its selfpropulsion and damping, that is translationally invariant 24 .A simple model that captures the basic physics is where r i is the position-vector for the ith agent in two spatial dimensions, α is a self-propulsion constant, β is a nonlinear damping constant, and λ is a coupling constant [8][9][10][11] .The total number of swarming agents is N, and each agent has unit mass.Beyond providing a basis for theoretical insights, Eq.( 1) has been implemented in mixed-reality experiments with several robotics platforms including autonomous ground, surface, and aerial vehicles 14,25 .
An example interaction potential that we consider in detail is the Morse potential, ( This is a common model for implicit interactions with local repulsion and attraction length scales, scaled as l and 1, respectively 10,23 .In the following, we assume that two interacting swarms are subject to the same underlying physics, Eqs.(1-2), but with different initial conditions.Each swarm will be initially separated by a large distance compared to the interaction length scales, l and 1.Therefore, the interaction force an agent feels will be at early times effectively confined to their own swarm, given the exponential decay with distance implied by Eq.( 2).
For the swarm flocks that are nearly aligned upon collision, the relevant initial-condition parameter is the distance between the two flocks as they approach x = 0, regardless of the direction of their velocities.This distance is often called the impact parameter, ∆y, in classical mechanics 26 , and it represents the closest distance the two flocks would approach in the absence of interaction forces.Depending on the value of impact parameter and the coupling strength after collision, the two flocks typically translate as one, scatter or form a combined mill.
Based on the previous numerical experiments for different values of ∆y and λ 7 , we are particularly interested in collisions that result in a swarm milling state.Representative parameters for generating milling states upon collision are: l = 0.5, α = 1.0, β = 5.0, λ = 0.1, ∆y = 4.0, which are assumed throughout this work.See 7 for further details.Our aim is to vary the repulsive strength and observe the changes in the dynamic behaviors of milling after the collision.Having noticed that changing repulsion drastically changes the postcollision dynamics, we then turned our attention to quantify the complexity of the dynamics so that we can try to deduce the dynamic consequences as we sweep parameters.One useful technique is Karhunen-Loeve decomposition, which we address next.

III. DIMENSIONALITY ANALYSIS: KARHUNEN-LOEVE DECOMPOSITION OF SWARMS OF N-AGENTS
Karhunen-Loeve (KL) decomposition is a procedure of mode expansion for spatio-temporal processes that extracts the relevant degrees of freedom of the dynamics of a large data set.We measured the number of Karhunen-Loeve(KL) modes needed to capture most of the the dynamic variance.It is an effective way to learn model reduction while describing the space in which the asymptotic dynamics resides.Known as the method of snapshots 27 , KL analysis correlates dynamics in time while averaging in space.Note that KL mode decomposition essentially is Principle Component Analysis in the time domain 28 .In order to identify correlations, we compute the covariance matrix, where each entry is correlation of two snapshots.The modes are defined by the data and constitute a natural coordinate system that approximates from time-series data optimally in space with the L 2 norm.From Eq. 1, assume we have N agents, which are sampled at M time points.Since both positions and velocities are changing with time, we define the sequence of data snapshots to include both position and velocity by letting where j ∈ {1, 2, ..., M} 29 .We assume the field can be expanded in a split time-space sum of modes given by The matrix A j,i stores the elements of the amplitudes at time t j as rows.That is, each row of the array is one mode.The vectors ϕ i are of the same dimension as of z(t i ) filed at a given snapshot.We define the residual using M modes over time as Minimizing E generates an eigenvalue problem for the first temporal mode; i.e., we solve where C is an M × M covariance matrix.The i j th element is given by <> denotes an inner product.Since the covariance matrix is symmetric, the eigenvalues are real, non-negative.Then we have ordering such that The corresponding eigenvectors, α i , are the temporal modes.
To find the spatial modes, we make use of the fact that the modes from C are orthonormal; Finally, we have We may now approximate the original spatial-temporal field by examining the total variance, which may be computed from the eigenvalues.Let and normalize the eigenvalues by setting We examine the cumulative sum of λi so that the total variance exceeds a threshold, T.Here T is arbitrarily chosen to be 95 percent.After transients are removed, in order to capture 95 percent of the full data, for instance, Fig. 3 indicates that 4 KL modes are needed for C= 0.51, 6 modes for C = 0.70 and about 140 modes for C = 0.89.In other words, one can see that when the post collision dynamics is periodic, it only requires four modes to describe the effective limit cycle behavior.As the repulsion strength is increased, the dynamics lies on a torus, which may be thought of as two coupled oscillators.However, the torus then breaks up, as the repulsion is further increased and the resulting transition is to a chaotic behavior, corresponding to a substantial change in KL dimension, from 4 to 140.The rapid change in KL dimension signals a dramatic transition in dynamical complexity.As we increase the repulsion strength, the combined milling behavior becomes qualitatively more complex exponentially from collective oscillations.Since the dynamics becomes quite complex and seemingly chaotic, we next turn to formally extracting whether or not is is technically chaotic for large repulsion.

IV. BINARY TEST FOR CHAOS: 0 -1 TEST
For its ease of use and vast applicability, we employed the 0-1 test 30,31 to distinguish chaotic and non-chaotic parameter regions in our systems.We input one-dimensional time series φ (n), n ≤ N ∈ N where N denotes the number of data points used.and it produces two-dimensional Euclidean extension P(n) and Q(n).The main idea is to assess whether the trajectories of the systems are bounded.Let φ (n) be an observed time-series of the underlying system. 30.For an arbitrary a > 0 that is randomly chosen, usually from (0, π), we define P(n) and Q(n) for each a as We then compute the mean square displacement defined by requiring N ≫ n.The limit is guaranteed by calculating Eq.14 only for n ≤ n cut where N ≫ n cut , and n cut = N 10 31 .The test for chaos is based on the growth rate of the mean square displacement as a function of n.We calculated the asymptotic growth rate of Eq.14 defined as Eq.14 is bounded if φ (n) is regular, or it grows linearly in time if φ (n) is chaotic 32 .Therefore, the 0-1 test states that a value of K ≈ 0 indicates regular dynamics, and K ≈ 1 indicates chaotic dynamics.For continuous systems such as differential equations, one might encounter an inaccurate test result due to oversampling 31,33 .We address this issue by incorporating two routine checks suggested in 34 .If the underlying dynamical system is regular, then its PQ plot is bounded; it is chaotic if the plot is diffusive.Figure 4 identifies with the dynamics observation from the time-series in Fig. 2. Note that the 0 -1 test only distinguishes chaos from regular dynamics, and is unable to detect quasi-periodicity.In our implementation of the 0 -1 test, we used 100 data points, i.e., N a = 100, from the time-series of y component of a randomly chosen agent.It was collected from a post-collision milling state for various of repulsive strength C 7 , excluding the transients.We then calculated K a for each a randomly chosen from [ π 10 , 2π 10 ] in accordance with an observation described in 35 to avoid oversampling (Fig. 5).The indicator for chaotic behaviors, denoted as K, of our system for a given repulsive strength is then K = median(K a ).We computed such K as a function of varying repulsive strength C.
Figure 6 shows that the dynamics is regular for the repulsive strength C that are in the range of 0.50 to 0. 81-83.As the value of C increases, the value of K starts to tend to 1, indicating chaos from onward around C = 0.84.This roughly corresponds to the increasing trend of number of modes predicted by the KL decomposition depicted in Fig. 3.However, the 0-1 test does not provide any information about the specific properties or behaviors of the system beyond its chaotic or non-chaotic nature.It cannot tell us anything about the attractors, periodic orbits, bifurcations, or other features that may be present in the system such as quasi periodicity.For the current work, we present our 0 -1 test result to illustrate what we have observed, while other analytical techniques will be needed as a next step to fully understand and characterize the dynamics of the system under discussion.

V. THE LEVEL OF EMBEDDING OF ONE SWARM IN ANOTHER: CONSTRAINED CORRELATION EMBEDDING
In addition to detecting chaotic behaviors and categorizing dynamics state of milling after the collision, we are interested to see the level of complexity of interaction between the two swarmsby looking at how one swarm embeds itself in another in a similar way as shown in the appendix of previous work 36 , the idea of which originated from fractal dimension.Fractal dimension was developed as a way to quantify complexity of nature in relatively simple mathematical equations 37 .The estimation of the complexity becomes difficult in the case of the systems in which many factors of nonlinear interactions occur.We use the notion of dimension in a locally constrained manner to evaluate and quantify the level of interaction between two swarms, calling it"constrained correlation embedding." Without loss of generality, we count the number of agents of swarm A in a local neighborhood of an agent of swarm B. Let N be a total number of agents in each Swarm.For an agent a i ∈ A and an agent b j ∈ B, with the radius of a neighborhood ε, we define γ(ε) is the total count of red agents within an ε neighborhood of each agent of swarm B. We vary the radius of ε L times, and aim to see how the relative inclusion grows by plotting log of γ(ε)versus log of ε.As the value of ε grows, the number of agents from swarm A in the ε ball centered on an agent from swarm B grows as the power law γ(ε) ∝ ε d .The log-log plot behaves almost linearly until it saturates out for a large ε.The curve saturates at large ε because the ε balls will then engulf the whole swarm thus it cannot grow any further.At extremely small values of ε, the only inclusion in each ε is just the k-th agent from swarm B itself.For the range of ε, let maximum distance and minimum distance of a given ith agent from swarm B with respect to all agents in swarm A be Max and Min, respectively.Using the power law, we divided the difference of Max and Min equally spaced logarithmically into 32 38 , i.e., L = 32.We computationally approximate d, the slope of the line from the region that is most linear in the log-log plot, and we defined it to be the contained correlation embedding, much like the way correlation dimension is obtained 39 while it is done piece-wise 40 but locally.In order to evaluate the con-strained correlation embedding, we sampled every 10th time steps, totaling 1800 points of the time -series data of interacting swarms consisting of 50 agents each.
First, our intuition tells us that a constrained correlation embedding should be between 1 and 2, since we are on the plane.Figure 7(A) is one representation of log-log plot at C = 0.81 and the slope d reflects an embedding that is fractal.The implication by CCE is that the agents of one swarm are embedded in another in a complicated way.Furthermore, we see from Fig. 7(B), that the embedding grows linearly between C = 0.65 and 1.05,The steeper slopes indicate fast changes in the nature of the complexity.On the other hand, we see that the complexity of milling behaviors of a combined swarms saturates beyond the value C = 1.05 while the complexities up to C = 1.05 intensify for reach sampled value of C from C = 0.81.The same is true for the values of C below 0.65, which we interpreted as change of complexity that is more stable or minimal., consistent with our initial intuition.
We have seen from Fig. 6 that the milling state of two swarms exhibit chaotic motion as the value of C grows.But the constrained correlation embedding further illustrates increasingly complicated motions within the chaotic region.From Fig. 7(B), In references to the supplementary material videos submitted, we have detected and identified the milling behaviors that are consistent with the constrained correlation embedding tendency of Fig. 7(B): For C ∈ [0.81, 0.95], we say the combined swarm forms fragmented milling where the swarm mills but it still maintains the grouping of two.Then for C ∈ (0.95, 1.35], combined swarm forms (true) milling, where the combined swarm becomes more even and unified.Notice that the milling at C = 0.95 still maintains grouping of 2 while repulsive strength C = 1.05 and C = 1.35, the milling behaviors changes only slightly, reflecting the steepness of the slope.(See supplementary material for videos of the dynamics with varying repulsion.)

VI. CONCLUSION, DISCUSSION AND FUTURE WORKS
We have studied post collision dynamics of combined swarm milling states.We began our analysis by extracting the essential modal dynamics-information from swarm collision data through a Karhunen-Loeve decomposition, then performed the 0-1 test for detecting chaos; we further defined and used a constrained correlation method to see the level of local complexity as a result of two swarms colliding.We have observed that a motion of periodic orbits transitions to a motion on a torus before undergoing chaotic motion for a given range of repulsion strength.For a single swarm dynamics, we have not observed such transition to chaos using Eqs.1,2.Furthermore, our constrained correlation embedding method reveals that formation of the milling itself becomes more complex as a function of repulsive strength.CCE enables us to categorize the milling as fragmented or true milling based on the range of repulsive strength; where the simple binary test of chaos does not reveal the details of the complexity of the milling within chaotic region.
One of the implications of the current result is that there is most likely a limitation to the theoretical mean-field analysis of the particular interacting swarms systems we studied.The KL mode decomposition shows that the modes needed to capture 95 percent of the data variance increases drastically as the value of C increases, and reveals a sharp transition to high dimensions.It might be an indication that the bifurcation analysis at higher values of C with the usual mean-field theory not be very meaningful due to spatial averaging.We likely need to come up with a way that is more accommodating to the complexity and dynamical structural analysis of interacting swarms.Our future work includes the investigation of the sequence that leads to chaos, as we have observed and have shown a possible Ruelle-Takens-Newhouse scenario of chaos 41 due to torus breakup, which is accompanied by a large change in modal dimension as a function of C.Moreover, knowing the region of chaos in the swarm systems is a great advantage from a control theoretical point of view 42 .We may be able to stabilize unstable states of colliding swarm dynamics using various control techniques, such as those used in 43 , which is based on analyzing the dynamics of the main coherent structure in the data represented by the highest energy of KL mode.
Lastly, recall that from Eq. 4, KL mode decomposition relies on space -time splitting hypothesis.The concept of spacetime splitting refers to the separation of variables into spatial and temporal components.In some cases, there are dynamical systems in which the temporal behavior of the system cannot be separated from its spatial behavior.In other words, the evolution of the system over time is inherently tied to the way in which it moves through space, such as nonlinear waves 44 , turbulence flow 45 and chemical reactions 46 .Overall, many real-world dynamical systems do not exhibit the property of space-time splitting 47 .Nonetheless, in the current work, we present a case where the rapid change in KL dimension correlates with a dramatic transition in dynamical complexity.Combined with the constrained correlation embedding analysis, we have enhanced insights to the complexity analysis of chaos.

VII. SUPPLEMENTARY MATERIAL
The videos show the milling state of a combined swarms consisting of N=50 agents each.The values of Repulsive strength C for five videos are: 0.51, 0.70, 0.89, 0.95, 1.05 and 1.35 respectively.Other representative parameters for generating milling states after collision are fixed, as mentioned in Sec.II.

Figure 1 .
Figure 1.Four time-snapshots for different values of t, showing each swarm with different colors: red and blue circles.Velocities are drawn with arrows.Swarm parameters are: C = 1.15,N = 50.

Figure 3 .
Figure 3.The 0-1 chaos test result (red circles) plotted over KL mode decomposition (blue circles).The x axis represents values of repulsive strength, C, and y axis is the number of KL modes needed to capture 95 percent of the total variance.The figures outside of the graph axis, (A)-(C), are the multi-agent attractor time series projected onto the x-y plane corresponding to the time-series in Fig. 2. The figures inside the axes are the corresponding fixed time snapshots showing the individual agents.This is to show that in real time, combined milling of two swarms are fragmented while it is not clear from the multi-agent attractor time-series projected onto the x-y plane.

Figure 4 .
Figure 4. P vs Q plot obtained for (A): C = 0.51, (B): C = 0.70 and (C): C = 0.89.P and Q are obtained in the process of 0 -1 test.Underling system is regular if the PQ plot is bounded, and it is said to be chaotic if it is diffusive.

Figure 6 .
Figure 6.0-1 test results obtained from plotting K = median(K a ) for different values of the coupling strength C. We interpret our result as chaotic for K ≥ 0.9, indicated by green horizontal line.