Collective behavior of self-propelled particles with multiple delays in cross-domain scenario

In cross-domain scenarios, the simultaneous presence of multiple sensing delays exerts a profound influence on collective behavior. Motivated by this, our paper presents a system based on self-propelled particles that consists of two swarms containing two intra-swarm sensing delays and an inter-swarm sensing delay. Three state emerges from the system, namely translating state, ring state, and rotating state. Employing mean field approximation and bifurcation analysis, we theoretically examine the parameter space’s boundaries that govern these three states. Our detailed findings emphasize that within the translating and rotating states, variations in the two intra-swarm delays can lead to the separation of the two swarms. Meanwhile, the extent of separation is influenced by the inter-swarm delay. Finally, we conduct numerical simulations to validate the accuracy of our theoretical derivations.


Introduction
Complex and dynamic collective behavior is a widespread phenomenon in nature.This phenomenon is observable in birds [1], fish [2], and even crowds [3].To investigate the mechanisms underlying collective behavior, researchers have developed various models based on self-propelled particles.Examples include the Reynolds model [4], the Vicsek model [5], and the Couzin model [6], among others [7][8][9][10].As a special case of self-propelled particles, active Brownian particles (ABPs) [11] and Active Ornstein-Uhlenbeck particles (AOUPs) [12] also provide useful theoretical tools for the study of active matter systems.A parental active model combining ABPs and AOUPs produced more complex collective dynamics, which may provide a more suitable model for a large class of biological and non-biological active matter system [13].These diverse approaches collectively converge on a similar conclusion: intricate collective behavior can emerge from simple local interactions among particles without the necessity of a central node.Building upon these pioneering works, researchers try to explain and reproduce different collective phenomena in nature through theory, particle simulation and even large-scale swarm robots [14][15][16].Concurrently, the application of these collective models for controlling swarm robots has gained significant attention as a contemporary research focus [17,18].
In both natural clusters and swarm robots, the presence of sensing delays among particles is inevitable.On the one hand, researchers have investigated the impact of sensing delays on collective behavior.Drawing from the Vicsek model, Mijalkov et al explored the effects of both positive and negative time delays on phototaxis swarm robots [19].Their findings revealed that a certain size of sensing delay can promote the aggregation of the swarm.Subsequent related works have verified and expanded this conclusion [20,21].Furthermore, utilizing the Vicsek model, Holubec et al demonstrated that the critical exponent of the Vicsek model can align with the critical exponent of biological swarms by introducing a certain degree of sensing delay [22].This result suggests the presence of potential sensing delays in natural biological swarms.On the other hand, researchers have also developed various algorithms to ensure the stability of swarms in diverse time-delay environments.For instance, Sun et al devised a cluster control methodology using an event-triggered control strategy, establishing sufficient conditions for achieving stable flocking motion in the presence of time-varying delays [23].Additionally, Yu et al introduced an innovative adaptive controller that enables swarms to attain consensus and reach specified areas despite time delays [24].Evidently, the presence of sensing delays among particles not only yields diverse collective behaviors but also poses challenges for the design of swarm controllers.Thus, further research in this area is warranted.
In practical applications, sensing delays can become more intricate.Cross-domain collaboration, e.g. two swarms performing the same task in two different domains in the aerial and underwater, is a research hot-spot in recent years.Due to the distinct environmental conditions of each domain, the intra-swarm sensing delay within these two swarms may differ.Moreover, for inter-group sensing between these swarms-the communication link spans two domains-the sensing delay could deviate from that encountered within each swarm.As far as we know, no work has been done to study the impact of sensing delays on collective behavior in such cross-domain scenario.Therefore, this paper intends to use a simple swarm model composed of only self-propelled force and symmetric attraction force to study the delayed collective behavior when there are multiple sensing delays exist in the system.We denote τ 1 and τ 2 as the sensing delays within the two swarms, and τ 3 as the sensing delay between swarms.The results reveal three distinctive states within our model: a ring state, a translating state, and a rotating state.Employing mean field approximation and bifurcation analysis, we theoretically scrutinize the boundary conditions associated with these three states.A comprehensive analysis of these states indicates that the separation occurs in the translating and rotating states when the intra-group sensing delays τ 1 and τ 2 differ.Finally, we validate our theoretical analysis using particle simulations experiments.

Model
In practical swarm scenarios, the introduction of sensing delays among particles can manifest in various ways.As illustrated in figure 1, consider a cooperative situation involving aerial and underwater components.Given the significant discrepancy in signal transmission rates between the air and underwater environments, the sensing delay within the two swarms differs.Let τ 1 and τ 2 denote the respective sensing delays within these two groups.Simultaneously, due to the dissimilarity between inter-group communication and intra-group communication, the sensing delay between the two groups contrasts with that within each group.We assume that the communication link lengths of the two swarms are kept constant and the sensing capabilities of the two swarms are the same, i.e. the inter-group sensing delay of both swarms is τ 3 Inspired by previous literature [25,26], we construct a system containing two swarms based on self-propelled particles with only symmetric attractive force.Then the motion equation of particles can be expressed as follows: ) Without losing generality, we only consider the motion of particles in two-dimensional plane.Therefore, r k i is a two-dimensional vector representing the position of particle i in swarm k.Our primary emphasis lies in investigating the impact of sensing delays on collective behavior.Consequently, we posit that both groups consist of an equal number, denoted as N, of particles, and the communication graph within each group is fully interconnected.The parameter a is introduced as a coefficient to modulate the interaction force among particles.Self-propulsion and velocity-dependent friction force is given by the first term (1 − ∥ṙ k i ∥ 2 )ṙ k i [27,28].If the particle is driven only by this term, the velocity of the particle will eventually stabilize at 1. Therefore, the particles can be maintained in a suitable velocity range to facilitate us to observe the collective behavior.The attraction force within the swarm is described by second term r k i (t) − r k j (t − τ k ), while the attraction force between swarms is described by the last term.

Mean field approximation
We first use the mean field approximation method to analyze the motion of the center mass of the swarm.Let denotes the center mass of swarm k.Thus the position of each particle can be rewritten as follows where δr k i represents the relative position of particle i from the swarm center k.And it is easy to get In order to facilitate the derivation process, let Substituting equation (2) into equation ( 1) can get the following equation Accumulating the above formula with respect to i = 1, . . ., N and dividing by N, we can obtain Suppose that N → ∞ and the perturbation term is ignored, the motion equation of the center mass of two swarms can be described as

Bifurcation analysis
Through a large number of particle simulations and the numerical simulation of equation ( 6) by Runge-Kutta method, we confirm that three states emerge from model ( 1), namely the translating state, the ring state, and the rotating state, as shown in figure 2. Although other states have appeared in particle simulations, they either only last a short time or appear in a harsh initial state.Therefore, these states are not studied in this paper.In this section, we first study the boundary of these three states in the parameter space by bifurcation analysis.To analyze the bifurcation of the equation ( 6), we must first identify the fixed point of the system.Evidently, the fixed point of the system is (R 1 , Ṙ1 , R 2 , Ṙ2 ) = (R 0 , 0, R 0 , 0), where R 0 represents a constant two-dimensional vector.This fixed point effectively characterizes a state wherein two swarms collectively move around the same central point.We confirm the existence of this state through numerical simulations and term it as the ring state.Subsequently, we explore the potential occurrence of Hopf bifurcation in the system.The presence of Hopf bifurcation signifies the emergence of a limit cycle, indicating that the system will exhibit periodic oscillations in addition to the ring state.By linearizing the equation ( 6) around the fixed point, we can obtain the characteristic equation, where J is the Jacobian matrix at the fixed point.In order to verify whether there is a characteristic root crossing the imaginary axis, λ = iω can be substituted into the characteristic equation.Then the characteristic equation can be written as follows, Separating the imaginary part and the real part of equation ( 8), we can get the following equation where . By adding the squares on both sides of equation ( 9), the following equation can be obtained, If the positive solution ω * > 0 exists in the above equation, then there is Hopf bifurcation in the system.By performing certain algebraic operations on equation ( 9), we can obtain the following conditions for the occurrence of Hopf bifurcation, Under the first curve of the Hopf bifurcation family, there is a fixed point in the motion equation of swarm center, that is, the ring state.The first Hopf bifurcation curve gives rise to the rotating state.The higher order Hopf bifurcation leads to faster angular velocity and instability of the rotating state.Equation ( 6) also has a translation solution, where R 10 and R 20 is the initial position of swarm 1 and 2, respectively.V 0 is the constant velocity of two swarms.This translation solution implies that two swarms move uniformly in the same direction at the same speed.Substituting equation ( 12) into equation ( 6), it can be obtained that the velocity V 0 satisfies The above equation should satisfy the condition: τ 1 < 4 a − τ 2 − 2τ 3 , that is, the translating state can only exist under this condition.While when τ 1 > 4 a − τ 2 − 2τ 3 , the solution V 0 is a pair of unphysical solutions with imaginary velocity.Therefore, the curve τ 1 = 4 a − τ 2 − 2τ 3 delineates a pitchfork bifurcation curve.In figure 3, we have plotted the Hopf bifurcation and the pitchfork bifurcation curves.The intersection of the Hopf bifurcation and the pitchfork bifurcation is referred to as the Bogdanov-Takens (BT) bifurcation point, occurring when ω = 0.The pitchfork bifurcation, in conjunction with the Hopf bifurcation and the BT point, delineates the boundaries of the three states observed in the system in the (a, τ 1 ) parameter space.
The bifurcation diagram of the system undergoes a highly intricate transformation as τ 2 and τ 3 vary, which results in the change in the feasible parameter space of the three states.At the same time, it can be seen that BT point also changes with the changes in parameters τ 2 and τ 3 .In figure 4, we have illustrated the values of a corresponding to the BT bifurcation points, denoted as a BT , under different τ 2 and τ 3 .The value of a BT decreases as τ 2 and τ 3 increase, a trend that is further confirmed by the data in figure 3. Additionally, when τ 2 and τ 3 are relatively large, no solutions exist for a BT , meaning that the Hopf bifurcation and the pitchfork bifurcation no longer intersect.It is also noticeable that the value of a BT is asymmetric with respect to τ 2 = τ 3 , which aligns with the distinct roles played by τ 2 and τ 3 in the system.
However, the bifurcation analysis on equation ( 6) solely on the motion of swarm center, disregards the influence stemming from the distribution of positions and velocities among particles.As a result, it does not provide an accurate portrayal of the boundary conditions for these three states.Consequently, we will delve into a more comprehensive examination of these three states, considering detailed conditions for their occurrence in conjunction with numerical simulations.

Collective motion
In the preceding section, we explored the conditions under which these three states transition through bifurcation analysis.In this subsequent section, our focus will shift towards a more intricate analysis of these three states, offering a deeper understanding of their characteristics and behaviors.

Translating state
When the initial velocities of the particles are similar and the initial positions are relatively concentrated, the system can evolve into a uniform translating state.Within this state, the two swarms move in the same direction at identical speeds.This phenomenon reveals that alignment interaction is not necessary to from flocking behavior, which is similar to the findings of several recent studies [29,30].In equation ( 13), we have derived the velocity of the swarms when the system stabilizes into the translating state.However, through particle simulation, we also observe that due to the distinct sensing delays of the two swarms, phase separation emerges within this state.Assuming the separation between the two swarms is denoted as d, by substituting equation ( 12) into equation ( 6), we can deduce that d adheres to the equation: As evident from equations ( 13) and ( 14), the main reason for the separation is the difference between τ 1 and τ 2 , and the degree of separation is also affected by τ 3 .To validate the theoretical outcomes stipulated by equations ( 13) and ( 14), we conduct numerical simulations.As illustrated in figure 5, the simulation findings broadly align with the theoretical predictions.

Ring state
When the system is stable in the ring state, all the particles will be evenly distributed on a circular orbit around a fixed point for periodic rotation.So we can get R 1 = R 1 = const.Substituting equation ( 5) into equation ( 4), we can get the following equation, The following equation can be obtained by simplifying equation ( 15) when assuming N → ∞, Transfer the above formula to the polar coordinate system δr i = (ρ i cos θ i , ρ i sin θ i ).Thus we have At the same time, the radius ρ 1 and ρ 2 are both constant with the same value in the stable ring state, that is, , and ρi = ρi = 0. Let θi = ω i , we can get the following equations: In a stable ring state, the rotational radius and angular velocity of the particles are exclusively dependent on the coefficient a and are unaffected by the sensing delay among particles.This outcome is rationalized by the dynamic motion of the ring state: when particles are uniformly distributed along the circular trajectory, the force exerted on the particles remains constant regardless of the magnitude of the time delay.For validation purposes, we employ numerical simulations, as illustrated in figure 6.The simulation findings corroborate the precision of the theoretical outcomes derived from equation (18).

Rotating state
The limit cycle generated by the Hopf bifurcation of the system indicates that the swarm center will undergo periodic circular motion.Different from the ring state where all particles are uniformly distributed on the ring, in the stable rotating state, due to the attraction force, all particles shrink to a point and move circularly.
Similarly, we can convert equation ( 6) into polar coordinate system (see appendix), ρk = In the stable rotating state, the pole diameter ρ 1 and ρ 2 are both constant but with the different value, and we set α = ρ2 ρ1 .The angular velocity of the two swarms is also a fixed value, that is, θ1 = θ2 = ω.Then we have θ k = θ k0 + ωt, where θ k0 is the initial angle of the center mass of swarm k.Through particle simulation, we find that there is a phase offset between the two swarms, let ∆θ = θ 1 − θ 2 represents the phase offset.Then for swarm 1, the equation ( 19) can be rewritten as ) . ( And the equation (19) for the swarm 2 can be rewritten as ) , ) . ( But by some algebraic operations of equations ( 20) and ( 21), we can get some transcendental equations to understand the properties of this state: sin ∆θ = where p(ω) = sin ω(τ1+τ2) 2 sin ω(τ1−τ2)

2
. By numerical methods we can solve the above equations and obtain the theoretical values of the rotating state properties.We have plotted the theoretical values and particle simulation results vary with τ 2 for different τ 3 in figure 7. It can be seen that when τ 3 = τ 1 , these four properties change more regularly with τ 2 , showing a trigonometric-like variation.At this point α is always less than 1, i.e. when τ 2 > τ 1 , the radius of swarm 2 is always smaller than swarm 1.In contrast, when τ 3 is larger than τ 1 , these four properties of the rotating state show a more complex change with τ 2 .

Conclusion and discussion
We have studied the impact of multiple sensing delays on collective behavior within a cross-domain scenario.We establish a system comprising two swarms with two intra-swarm sensing delay and an inter-swarm sensing delay.Within this framework, the system gives rise to three distinct states: a translating state, a ring state, and a rotating state.Theoretical analysis and numerical simulations are performed to study the properties of these three states.Results indicate that the two swarms separate from each other within both the translating and rotating states when the two intra-swarm delays are not the same.While the inter-swarm delay affects the degree of separation.The radius of the ring state solely correlates with the coefficient a, leading to both groups adopting identical circular trajectories within a stable ring state.Our results may guide the design of real flocking control methods and explain the behaviors in real cross-domain swarms.
Numerical simulations and theoretical derivations cross-validate that the interaction strength coefficients a together with the delay parameters lead to different states of the system, but their causes may be different.For the translating and rotating states, the relative positions of the two swarms remain constant, i.e. the interaction force between the two swarms is also fixed.Although the two swarms tend to converge at the same point under the attraction force, there is a difference in the time delay, so that the two swarms essentially reach agreement under the delayed information, which manifests itself as the occurrence of a separation of positions.For the ring state, on the other hand, since the particles of both swarms are uniformly distributed throughout the ring orbit, the force on the particles is the same at every moment.Therefore the difference in delay will not affect the nature of the ring state.Thus no separation of the two swarms occurs in the ring state.And the bifurcation, i.e. the boundary between the three states, is due to the trade-off between self-propelled and attraction force.Taking the translating state as an example, in this state the direction of the self-propelled force and the attraction force are opposite to each other, so when the attraction force is smaller than the self-propelled force, the system is able to maintain a stable translating state, while when the parameter increase causes the attraction force to be larger than the self-propelled force, the translating state will disappear, i.e. a pitchfork bifurcation occurs.
Our work introduces new insights into the understanding of collective behavior in scenarios involving multiple sensing delays, especially in cross-domain situations.However, there remains ample room for further investigation.Firstly, the model presented in this study only employs a simple symmetric attraction force.In real-world scenarios requiring complex tasks, swarm models are often more intricate.It remains uncertain whether the conclusions drawn in this paper can be extended to such complex models.Secondly, obstacle avoidance is a crucial capability for practical applications.Despite prior demonstrations that short-range repulsive forces do not fundamentally alter collective behavior [31,32], our multiple sets of swarm robots simulation experiments in Pybullet [33] suggest that obstacle avoidance might also be influenced by particle dynamics.Thirdly, in this study we assume that the sensing delays are all fixed, whereas time-varying delays are much more common in practice.From the bifurcation diagram (figure 3) and the properties obtained for three states, it can be seen that when the sensing delays change it has a large impact on the bifurcation as well as the nature of the system.And it is also unknown whether the time-varying sensing delays produce new states.Similarly, different sensing delays for inter-group are a more practical assumption but it is simplified in this study.How to carry out research on the above problems and how to carry out real swarm robots experiments will be the focus of our future research work.

Figure 1 .
Figure 1.Description of sensing delay in cross domain scenario.

Figure 2 .
Figure 2. Dynamic motion of the three states.(a) Translating state.(b) Ring state.(c) Rotating state.Under the action of attraction, the two swarms of translating state and rotating state converge to two points.

Figure 5 .
Figure 5. (a) The speed of the swarm center V0 and (b) the distance between two swarms d as a function of τ 2 for different τ 3. The lines are the theoretical results.Star and square markers are the numerical simulations.The parameter of the numerical simulations are set to N = 100, a = 0.1, τ1 = 0.5.The simulation experiments in this paper are carried out 5000 steps, and all the data are obtained from 30 independent experiments.

Figure 6 .
Figure 6.The comparison of the theoretical result (red line) and the particle simulation results (blue star).The parameters are set to τ1 = τ2 = τ3 = 1.

Figure 7 .
Figure 7. (a) The angular velocity ω (b) the phase offset of two swarms sin ∆θ, (c) the ratio of two radius α, and (d) the radius of swarm 1 ρ1 as a function of τ 2 corresponding to different τ 3. The lines are the theoretical results.The star, circle, and diamond markers are the simulation results.The parameters are set to a = 2.5, τ1 = 3.

Figure 8 .
Figure 8.Typical snapshot for the rotating state.The positions of the two swarms at the current time are denoted by solid circles, while the positions at the past are denoted by hollow circles.