EMERGENT DYNAMICS OF AN ORIENTATION FLOCKING MODEL FOR MULTI-AGENT SYSTEM

. We study the orientation ﬂocking for the deterministic counterpart of a stochastic agent-based model introduced by Degond, Frouvelle and Merino- Aceituno in 2017, where the orientation is deﬁned as a SO(3) matrix. Their proposed model can be reduced to the other collective dynamics models such as the Lohe matrix model and the Viscek-type model as special cases. In this work, we study the emergent dynamics of the orientation ﬂocking model in two frameworks. First, we present suﬃcient conditions leading to the orientation ﬂocking when the natural frequency matrices are identical. To be precise, we prove that all orientation matrices asymptotically converge to the common one, and the spatial position diameter remains uniformly bounded. Second, we show the emergence of orientation-locked states for non-identical natural frequency matrices, that is, the diﬀerence of any two orientation matrices tends to the deﬁnite constant matrix. On the other hand, we establish the ﬁnite-in-time stability with respect to initial data of the proposed orientation ﬂocking model. We also present the numerical results consistent with our rigorous analysis. Our work remains valid even for dimensions greater than three.

1. Introduction. Collective dynamics of self-propelled agents in biological complex systems can be easily found in our nature, for instance, swarming of fish, aggregation of bacteria, flashing of fireflies, etc. Indeed, there are many agentbased models describing the flocking and synchronization phenomena, such as the Kuramoto dynamics [30,31], the Cucker-Smale model [10], the Vicsek model [38], the swarm sphere model [33], etc, and these models have been extensively studied in literature [1,22,34]. In this paper, we are mainly concerned with an agent-based model for alignment of body attitudes which was proposed in [11] where the states of agents are described by the positions of their center of mass and body attitudes. For simplicity of modeling, other detailed internal structures are ignored at the level of modeling. To be specific, directions of agents are initially determined and move with the same speed, and agents try to adjust their body attitudes (or orientations) with their neighboring agents towards average orientation. Let x i ∈ R d be the position for its center of mass of i-th agent in d-dimensional Euclidean space, and A i ∈ SO(d) denote the orientation of the i-th agent. Then, the dynamics of (x i , A i ) is governed by the following coupled ODE system, called the orientation flocking model (OFM): Here, the term A i e 1 describes the direction of displacement for i-th agent where {e i } d i=1 is the canonical basis, and v 0 denotes a constant common speed of the agents. The skew-symmetric matrix H i denotes a generalized frequency-like matrix, κ measures the coupling strength between the agents, and the communication function ψ = ψ(r) describes the capacity of channel between agents determined by the relative spatial distances. In addition, we also assume that there exist two positive constants ψ m and ψ M such that 0 < ψ m ≤ ψ(r) ≤ ψ M , r ≥ 0. (2) Note that the condition (2) might appear to be a strong restriction from the modeling point of view, since any two particles can strongly interact even if their relative distance is large. In fact, one can consider the channel capacity in the Cucker-Smale model whose form is ψ(s) = (1 + s 2 ) −β for a long-ranged interaction β ∈ (0, 1/2) and a short-ranged interaction β ∈ [1/2, ∞).
In [11], the authors proposed two types of body attitude models on SO (3). For the first one, the dynamics of A i is given by the following gradient flow of the weighted average M i : where (ν, K) have a same role of (2κ, ψ) in (1), and P T A i denotes the projection operator on the tangent space T Ai so that A i always belongs to SO(d) for all time: for A ∈ SO(d) and M ∈ R d×d ,

EMERGENCE OF ORIENTATION FLOCKING 2039
The second one is also formulated by the gradient flow, however using the polar decomposition of M i : where PD(M ) denotes the orthogonal matrix which comes from the polar decomposition of the matrix M . Compared to the model (3), our flocking model (1) has an extra intrinsic component H i which plays the role of natural frequency matrix in the analogy with the Kuramoto model. We also refer the reader to [12] for an equivalent description of (5) using quaternions. For our model, we choose (1) instead of (5), since the polar decomposition process contains an implicit process. In order to provide explicit calculation, we only consider (1). See Proposition A.4 in [11] in which they showed that the trajectory can follow a geodesic in some special case. We also refer the reader to [35,37] for attitude synchronization models which contain the dynamics of SO(d) and [5] for a system of interacting orientable agents.
On the other hand, the authors in [24] introduced a new class of the matrix Lie groups, namely Lohe group, which can be generated by the continuous matrix-valued ODE system: For a general matrix Lie group G, H i belongs to the corresponding Lie algebra g, and we impose the following structure condition on G: so that a solution of (6) stays on the given Lie group G for all time. In particular, if we choose G = SO(d), then (6) becomes (1) 2 with ψ(r) ≡ 1 and A i = X i . See also [13] for the generalized model of (6). In this paper, we give a link between two aforementioned works [11] and [24]. The authors of [11] proposed the ODE systems (3) and (5) with a noise effect, formally derived their corresponding kinetic and hydrodynamic models, and studied the selforganized hydrodynamics for body attitude coordination. Hence, the analysis of the ODE model is left as a blank and we thus focus on the emergent dynamics of the agent-based model (1). Indeed, since system (1) is a matrix-valued ODE system, the techniques developed in [24,29] can be adopted throughout the paper. We classify our system (1) into two types: For the Lohe matrix model [32], H i plays the role of constant Hamiltonian associated with i-th quantum subsystem. In what follows, we will use the following notation: The main results of this paper are two-fold. First, we describe some of the emergent dynamics of system (1) with both identical and non-identical cases. For the identical case, we present a sufficient framework leading to the orientation flocking (see Definition 2.1): there exist positive numbers x ∞ > 0 and C > 0 such that where D(X(t)) and D(A(t)) are maximal diameters: Here, |·| and · are 2 -norm and Frobenius norm in R d and SO(d), respectively (see the end of this section). For the non-identical case, we show that the orientationlocked states (see Definition 2.1) emerge asymptotically in a large coupling regime.
In this case, we cannot guarantee that the diameter of position variable is uniformly bounded. Second, we present the finite-in-time stability estimate which is valid on any finite time interval. To be more specific, for a given T ∈ (0, ∞), we can find a positive constant G = G(T ) which does not depend on the number of particles N and initial data such that for any two solutions (x i , A i ) and (x i ,Ã i ), the following estimate holds for t ∈ [0, T ): As a direct application of this stability estimate, we can derive the corresponding kinetic model (39) of (1) and obtain a global measure-valued solution to (39) as in the Cucker-Smale model [17] using a particle-in-cell method. The rest of paper is organized as follows. In Section 2, we study basic properties of our system, such as an invariance property and present connections between our model and other collective models, e.g., the generalized Lohe matrix model and the Viscek-type flocking model. In Section 3, we present orientation flocking estimates with both identical and non-identical cases. In Section 4, we show that our system is stable with respect to the initial data in any finite time interval. In Section 5, we provide several numerical simulations consistent with the theoretical analysis obtained in Section 3. Finally, Section 6 is devoted to a brief summary of our main results and discussion of future works.
Notation: Throughout this paper, we use the 2 -norm and Frobenius norm for 2. Preliminaries. In this section, we briefly review basic properties of system (1) and study relations between our system (1) and other existing flocking models such as the Lohe matrix model and the the Viscek-type flocking model. First, we introduce several concepts related to orientation flocking as follows.
3. (Orientation locking): System (1) achieves the (asymptotic) orientation locking, if there exists a constant matrix F ∞ ij for each i, j ∈ {1, · · · , N } such that where e is any unit vector in R d .
, the following relations hold: Below, we recall the Grönwall-type lemma in [24]. For the convenience of the reader, we provide its proof.
Proof. We use Lemma 2.2 to see d dt Hence, if X > 0, then one has d dt X ≤ −γ X + Y , t > 0.
Next, we show the positive invariance of (1). More precisely, we prove that if A 0 i ∈ SO(d), then A i (t) ∈ SO(d) for all t ≥ 0. Proof. First, we recall the definition of G = SO(d) and its corresponding Lie algebra g = so(d): If the right-hand side of (1) 2 belongs to g, then we see thaṫ On the other hand, since H i is a skew-symmetric matrix and the following relation holds: we see that the right-hand side of (1) 2 indeed belongs to g. Hence, we conclude that A i ∈ G for all t ≥ 0.

Remark 2. (i)
In [11], the authors formulate model (1) in the following way: where P T A i defined in (4) denotes the projection operator onto the tangent space T Ai G. Hence, the solution A to (1) lies on SO(d) N for all time. Moreover in [25], they showed that (6) on the unitary group or the special orthogonal group also can be represented as a gradient flow.
(ii) For the global well-posedness of system (1), we just note that SO(d) is a compact manifold. Thus, global well-posedness directly follows from Remark 2.4 in [24].
Below, we study the solution splitting property of system (1) with identical H i 's: In this case, we consider two subsystems on SO(d) and R d × SO(d), respectively: and Then, we introduce two solution operators R H (t) and L(t) of two subsystems (9) and (10), respectively. For B i ∈ SO(d), R H (t) is a solution operator corresponding to (9): is a solution operator corresponding to (10): . Then, the solution A to (1) with (8) can be represented via a composition of two solution operators R H (t) and L(t) as follows.

Proposition 2. Suppose that the initial data and H i satisfy
for all i = 1, · · · , N , and let (X, A) be a global solution to (1). Then, we have Proof. Although the proof can be found in Proposition 4.1 in [24], we provide it for the convenience of readers. We denote Thus, the relation (11) holds.
Remark 3. Note that our solution splitting does not hold for the spatial variable Since the additional term Hy i appears in the relation above, our solution splitting property cannot be extended to the variables (x i , A i ).
Next, we introduce Lyapunov functionals for the flocking estimate and study their rates of changes over time. First, we set where · is the Frobenius norm defined in (7) and we use the relation Below, we derive differential inequalities for D(A) and D(A,Ã). Lemma 2.4. Let (X, A) and (X,Ã) be any two global solutions to (1)- (2). Then, the following estimates hold: Proof. Since the proofs are rather long, we present them in Appendix A.

Previous results.
We here provide previous results on the generalized Lohe matrix model (6) on G = SO(d) which can be regarded as a space homogeneous case of the OFM (1). Thus, it is worthwhile to mention the previous results of (6). For this, we denote the maximal quantities: Below, we state the main result without the proof and it can be found in Theorems 4.2 and 5.1 of [24]. Let X i be a solution to (6). Then, the following assertions hold. (i) Suppose that system parameters and initial data satisfy Then, the following estimate holds: (ii) Suppose that system parameters and nitial data satisfy Then, we verify lim

2.3.
Relations with other collective models. In this subsection, we briefly discuss the connections between our model and other synchronization models in literature.
For the case where the communication weight is a constant, system (1) 2 reduces to the generalized Lohe matrix model [24] with a network topology (a ik ): Now, we consider a two-dimensional case d = 2. In this case, we can parameterize H i and A i ∈ SO(2) as Then, we substitute the above ansatz into system (1) to recover the Viscek-type flocking model [19]: For ψ ≡ 1, (12) 2 is exactly the same as the Kuramoto model with zero natural frequencies which has been extensively studied in [2,3,4,8,15,16,14,18,21,22,23,28,36]. On the other hand, (12) can be obtained from Cucker-Smale flocking model with unit speed constraint in [7,26].
Before we end this section, we present the Grönwall-type estimates to be used later.
Lemma 2.6. Let y = y(t) be a nonnegative C 1 -function satisfying the following Riccati-type differential inequality: where p, q and r are positive constants.
(i) Suppose that r = 0. Then, one has (ii) Suppose that r > 0 and p 2 − 4qr > 0, and let y − and y + be two distinct positive roots of a quadratic equation qy 2 − py + r = 0: Furthermore, if we assume that y(0) < y + , then there exists a finite positive entrance time T * such that Proof. Although the proof can be found in Lemma 3.7 of [6], we provide it for the consistency of the paper.
(i) We set u = 1/y to derive the inequality for u and then directly integrate the resulting relation to obtain the desired estimate.
(ii) We split the proof into two cases.
• Case A: suppose that y(0) ≤ y − . For t = T such that y(T ) = y − , it follows from Hence, y(t) is restricted in the interval [0, y − ] for all time, i.e., y(t) ≤ y − for all t ≥ 0.
• Case B: suppose that y − < y(0) < y + . Then, it follows from (13) that y(t) starts to decrease strictly. It follows from Proposition 3.1 in [8] that there exists a finite entrance time T * such that Finally, we combine the cases A and B to obtain the desired result.
3. Emergence of orientation flocking. In this section, we present the orientation flocking estimate for (1). In what follows, we consider both identical and non-identical cases. Before we state the main result, we set where ψ m and ψ M are defined in (2) and ∆(ψ) is assumed to be positive from the condition (14) 1 .
3.1. Identical case. In this section, we consider the identical case, i.e., H i ≡ H. Then, we state our first main result on the emergence of orientation flocking of system (1).
Theorem 3.1. Suppose that system parameters and initial data satisfy Note that the condition (14) yields α > 0. Thus, we apply Lemma 2.6(i) to find the desired exponential decay of D(A): (ii) (Uniform boundedness of D(X)): For the group formation condition, we roughly estimate (16) as On the other hand, we have Now, we use (17) and (18) to obtain Thus, the group formation holds.

3.2.
Non-identical case. In this subsection, we study the non-identical case, i.e., D(H) > 0. Although we cannot guarantee the orientation alignment, we here show that orientation-locked states can emerge from some well-prepared initial data. Before we begin our discussion, we introduce the handy notations: Note that Λ(κ, ψ, H) can be positive for large κ > 0 and this can be guaranteed under the assumption (20). Below, we show that if the coupling strength κ is sufficiently large and initial diameter D(A 0 ) is small enough, then D(A) can be made sufficiently small.
Proposition 3. Suppose that system parameters and initial data satisfy and let (X, A) be a global solution to (1)- (2). Then, there exists a finite entrance time T * such that D(A(t)) < β − , t > T * .
On the other hand, the condition (20) implies, (20) and (21), we apply Lemma 2.6(ii) to conclude that there exists a finite entrance time T * such that y(t) < y − or D(A(t)) < β − , t > T * . (20) is necessary when we apply the Grönwall-type inequality Lemma 2.6 (ii). More precisely, the conditions (20) 1 and (20) 2 guarantee the existence of two positive roots of −py + qy 2 + r = 0 in (23) and the negativity of −p, respectively.

Remark 4. (i) The condition
(ii) For β − in (19), one has Hence, if the coupling strength is sufficiently large, then we can make β − sufficiently small and in turn, D(A(t)) can be also made sufficiently small.
We now arrive at our second main result which states that system (1) tends to the orientation-locked states under the constant communication weight and large coupling strength regime.
Theorem 3.2. Suppose that system parameters and initial data satisfy and let (X, A) and (X,Ã) be any two global solutions to (1)- (2). Then, there exists a constant limit matrix F ∞ ij ∈ SO(d) such that lim Proof. We split the proof into two steps.
• (Step A: convergence of D(A,Ã) towards zero): it follows from the condition (25) and Proposition 3 that there exists T * > 0 such that Then, we use Lemma 2.4(i): for t > T * , The coefficient −γ in (26) is negative due to (24). Hence, we find the relation which yields the zero convergence of D(A,Ã): • (Step B: convergence of A i A −1 j ): since we are only interested in the large-time behavior, without loss of generality, we may assume T * = 0. For any T ≥ 0, {A i (t + T )} is also a solution with the shifted initial data {A i (T )} due to system (1) being autonomous. Hence, we setÃ i (t) = A i (t + T ) to apply the result of Step A: To verify the convergence, we discretize (27) by setting t = n and T = 1. Then, (27) can be rewritten as and by the induction argument for m ∈ N, hence {A i A −1 j (n)} is a Cauchy sequence in the compact set D(A(t)) ≤ β − . Therefore, there exists a limit F ∞ ij ∈ SO(d) such that lim and its convergence rate is exponential due to (27).

Finite-in-time stability.
In this section, we study the finite-in-time stability of system (1) with respect to initial data. Before we proceed further, we present definition of stability of (1) with respect to the initial data. For , we set Z := (Z 1 , · · · , Z N ) and define an 2,∞ -like norm denoted by · 2,∞ : which does not depend on the number of oscillators and the initial data such that for any two global solutions Z andZ, if the following estimate hold, then we say that system (1) is finite-in-time stable with respect to the initial data.
Remark 5. If T can take on ∞, then we call the estimate (28) as the uniform-intime stability with respect to initial data. In fact, the uniform-in-time stability has been established in literature, for instance, [20] for the swarm sphere model and [27] for the Cucker-Smale model.
Before we begin a lengthy estimate, we briefly outline our strategy to derive the finite-time stability estimate (28) for any two solutions Z andZ of (1).
• (Estimate of |x i −x i |): We use (1) 1 and (29) to derive which yields the stability of spatial differences. Now, let us begin our stability estimate. Let Z andZ be two solutions to (1). Then, we use (40) to finḋ where ψ ik := ψ(|x i − x k |) andψ ik := ψ(|x i −x k |). For notational simplicity, we simply write Note thatP t i = P i and Q t ik = Q ki . Then, we calculate (30) 1 ×Ã −1 i + A i × (30) 2 to find the dynamics of P i : Below, we estimate the terms I 1k , k = 1, 2, 3, separately.
Lemma 4.2. Let Z andZ be two solutions to (1). Then, the terms I 1k , k = 1, 2, 3 can be rewritten as follows. (i) Proof. (i) By direct calculation, we see After tedious algebraic manipulation, we have (iii) By the same manipulation as in (ii), we have Lemma 4.3. Let Z andZ be two solutions to (1)- (2). Then, the following Grönwall inequalities hold: it can be directly integrated to yield that (ii) In (31), we use Lemma 4.2 to obtain We change the roles of A i andÃ i to see which also can be obtained by taking tilde notation into (32) if we allow to abuse the tilde notation. We add (32) and (33) to obtain On the other hand, sinceP i = P t i and the following relation holds: or equivalently, P iPi = −P i −P i , we take a trace on both sides of (34) to see Since ψ ik = ψ ki andψ ik =ψ ki hold, we write S ik := ψ ik Q ik . Then, (35) can be rewritten as Note that Hence, we have We obtain from (36) that Then, we give a rough estimate as Finally, we use Lemmas 4.2 and 4.3 to obtain the following theorem.
Proof. Note that there exists a uniform constant C > 0 such that where . We write P := max 1≤i≤N P i and use Lemma 4.3(ii) and (38) to obtain d P dt ≤ 4κ(1 + d)(C + ψ M ) P =: δ P , t > 0.
Then, Grönwall's lemma yields that for a given T ∈ (0, ∞), which can be rewritten as On the other hand, for spatial variables, we use Lemma 4.3 (i) to obtain Hence, we can find a constant G(T ) := max{1,

Remark 6.
The issue when one concerns with the (uniform-in-time or finite-intime) stability is whether the maximal time T can be ∞ or not. In addition, if two orientation matrices are not the same asymptotically, that is, if A i − A j does not converge to zero as time t goes to infinity, then |x i − x j | tends to infinity. Thus, when we deal with the stability with respect to the initial data, we always consider the case of the orientation alignment. In fact, Theorem 4.4 fails to give uniform-intime stability estimate due to the following technical reason. If we apply Lemma 2.3 in (32), then we obtain the following Grönwall-type inequality: This, however, does not guarantee the uniform stability not only for whole time but also for finite time. This failure comes from the zeroth order terms such as Q ik and Q ki in the right-hand side of (32). To overcome this, we split Q ik into higher order terms of P i as in (37).
Since the Frobenius norm of M d (R) corresponds to the 2 -norm of R d 2 , (1) can be regarded as a system defined on R d × R d 2 . Hence, if we adopt the Wasserstein framework in [27] to apply the finite-in-time stability, then we can identify the kinetic model of (1), which is valid on any finite time interval, as follows: for the one-particle distribution function f = f (x, A, H, t), where the projection matrix is defined as P A (A * ) = 1 2 (A * − AA t * A). We refer the reader to [11] for a detailed calculation of ∇ A and integration with respect to A.

5.
Numerical experiments. In this section, we provide several numerical simulations for the emergent dynamics of our system (1) for both identical case (H i ≡ H) and the non-identical case (D(H) > 0). We used the fourth-order Runge-Kutta method for numerical integrations.
Here, the initial data and the entries of skew-symmetric matrices H i are chosen randomly. Note that since we only consider a constant communication function ψ, the second equation becomes autonomous system in terms of the orientation A i . Thus, the initial data of position is only used to indicate relative positions of particles. Figure 1(a) describes an agent moving with velocity A i e 1 (red arrow) and the orientation determined by two vectors A i e 2 and A i e 3 . Note the positivity of orientation. In Figures 1(b)-(d), we can observe the orientation flocking of agents. In particular, Figure 1(d) exhibits the comparison of the decay of D(A(t)) obtained numerically and the analytical result in Theorem 3.1 which supports (16).
6. Conclusion. In this paper, we have studied flocking behaviors of the agentbased model [11] for the orientation flocking of a multi-agent system. We present sufficient frameworks leading to the emergence of orientation flocking and orientationlocked state. To be more precise, orientation flocking can emerge for the identical case H i ≡ H for all i = 1, · · · , N, that is, all orientations become eventually the same, as time goes to infinity. On the other hand, the orientation-locked state can occur for the non-identical case in a large coupling regime. In this case, the difference between any two orientation matrices converges to the definite constant nonzero matrix, whereas it converges to zero matrix for the case of orientation flocking. In addition, we show that our system is finite-in-time stable with respect to the initial data in any finite time interval, and as a corollary, we can identify its kinetic model in a Wasserstein framework. We also provided several numerical simulations supporting our analysis. Of course, there are still lots of interesting open questions such as the extension of stability estimate to the whole time interval, global-wellposedness, regularity and emergent dynamics of the kinetic model. These issues will be addressed in future works.
Appendix A. Proof of lemma 2.4. In this appendix, we provide the proof of Lemma 2.4.
(i) (Derivation of Grönwall'is differential inequality for D(A)): Note that We first observe that to obtain the dynamics for inverse matrix A −1 i : Let i, j be indices in {1, · · · , N }. Then, we have On the other hand, We apply Lemma 2.3 with X = Q ij to obtain For given t ≥ 0, we choose the index (i t , j t ) such that Then  This completes the derivation of the differential inequality (i) in Lemma 2.4.
(ii) (Derivation of Grönwall's inequality for D(A,Ã)): It follows from (41) that where M ij k is defined as Moreover, we rewrite M ij k in terms of differences between Q mn andQ mn : Since we assume that ψ ik is a constant, we have In (44), we use the relation (45) to obtain By a similar argument as in (i), we choose the indices (i t , j t ) such that