Emergence of anomalous flocking in the fractional Cucker-Smale model

In this paper, we study the emergent behaviors of the Cucker-Smale (C-S) ensemble under the interplay of memory effect and flocking dynamics. As a mathematical model incorporating aforementioned interplay, we introduce the fractional C-S model which can be obtained by replacing the usual time derivative by the Caputo fractional time derivative. For the proposed fractional C-S model, we provide a sufficient framework which admits the emergence of anomalous flocking with the algebraic decay and an \begin{document}$\ell^2$\end{document} -stability estimate with respect to initial data. We also provide several numerical examples and compare them with our theoretical results.


1.
Introduction. Collective movements of a multi-agent(particle) system can be easily observed in our daily life. To name a few, herding of sheep, flocking of birds and swarming of fish, etc., [4,5,19,24]. These collective phenomena are often called flocking in broad sense, and it refers to the phenomena where self-propelled particles or agents adjust their movements to reach a velocity consensus using a limited environmental information and simple rules. In literature [5,19,24], several mechanical models were proposed for such a flocking modelling. Among others, we are mainly interested in the particle model proposed by Cucker and Smale [6]. Now, we begin with our discussion with the C-S model. Let (x i , v i ) ∈ R d × R d be the position and velocity of the i-th C-S particle, respectively. Then, the phasespace dynamics of the C-S particles is governed by the following Cauchy problem: where · is the 2 -norm in R d , and ψ : [0, ∞) → R is a communication weight function measuring the strength of coupling, and it is assumed to be bounded, nonnegative, Lipschitz continuous and monotonically decreasing so that global wellposedness is guaranteed by the standard Cauchy-Lipschitz theory. In recent years, there were extensive research on the emergent dynamics of (1) from various points of view (see a recent survey article [5] for details). Among others, if the monocluster flocking occurs for the C-S model, its decay mode should be exponential due to the linear velocity coupling in (1) 2 as noticed in [6,13,14]. This is the point where we address the following simple questions.
• (Q1): Can we modify the dynamics of (1) so that the resulting dynamics exhibits a slow anomalous mono-cluster flocking with an algebraic rate?
• (Q2): Can we incorporate a non-Markovian effect (memory) in (1)? In previous literature, all research on the C-S model use a Markovian approach by ignoring memory effect from the beginning. However, as can be seen in our daily life, one sometimes needs to make a choice(decision) based on one's own past memory, at least short-term memory. Thus, it is worthwhile to study such a non-Markovian effect in the context of flocking dynamics. For this, we introduce a fractional C-S model by replacing the temporal derivative with the Caputo fractional derivative: where D α * denotes the Caputo fractional derivative of order α ∈ (0, 1). So far, there are very few literature available for the C-S model with fractional derivatives, e.g., discretization of the fractional C-S model [10,11] and optimal control problem for C-S model with Riemann-Liouville derivative and constant communication weights [17]. Despite of recent active research on the fractional dynamical systems in science and engineering [7,9,15,23,25], as far as the authors know, emergent dynamics of the fractional flocking models has not been addressed in a full general setting.
In the previous works [10,11,17], authors discretized the fractional C-S model or assumed that communication weights are constant, which yields linear systems. So they could obtain the desired results from the analysis of the linear systems. In contrast to the aforementioned linear results, for our proposed nonlinear model, previous linear theories can not be applied to our case as it is. This is one of the main difficulties for the fractional Cucker-Smale model with nonconstant communication weights. To overcome this difficulty, we used the arguments in [21] together with the matrix analysis to derive our desired flocking estimates.
Moreover, due to Caputo's fractional derivative, we may expect that the resulting fractional C-S model might exhibit a slow relaxation toward the flocking state unlike the exponential flocking in the original C-S model (1). This algebraic relaxation is what we can additionally obtain by the introduction of fractional calculus in the flocking dynamics. Similar issues have been discussed in the fractional Kuramoto model [12], where the synchronization occurs algebraically fast. Now, we present a definition of mono-cluster flocking for the fractional C-S model as follows.
} be a C-S ensemble whose dynamics is governed by (2). Then, it exhibits a mono-cluster flocking if and only if the following two conditions hold.
Before we discuss our main results, we set the position and velocity configurations as follows: The main results of this paper are two-fold. First, we present a sufficient framework for the slow relaxation to the fractional C-S model. For a communication weight whose maximum and minimum are close or for an initial configuration where particles are initially located nearby with small velocity fluctuations, we derive a mono-cluster flocking estimates (see Theorem 3.3 for details): where E α (·) is the Mittag-Leffler function introduced in Definition 2.2, and λ is a positive constant to be defined later. Then, it follows from the property (3) of the Mittag-Leffler function that Second, we present an 2 -stability estimate with respect to initial data. Under the same condition as in the first result and min{ V 0 , Ṽ 0 } 1, we show that there exists a positive constant G independent of time t such that sup t≥0 The rest of this paper is organized as follows. In Section 2, we study several basic concepts on the fractional calculus, basic properties of the fractional C-S model and provide previous results for the Cucker-Smale model. In Section 3, we provide the results for the emergent behavior of the fractional C-S model. First, we consider allto-all case where the communication weights between agents are constant. Second, we provide our main results under general settings. In Section 4, we consider the 2 -stability of two solutions with respect to initial configurations. In Section 5, we provide numerical simulations for the fractional C-S model to support our analytic results. Finally, Section 6 is devoted to a brief summary of our main results and discussion on possible future works.

2.
Preliminaries. In this section, we provide basic results on the fractional calculus for self-containedness, and present basic properties of the fractional Cucker-Smale model (2).
2.1. Fractional calculus. In this subsection, we present several definitions and useful properties from fractional calculus to be used throughout the paper. First, we provide definition of Caputo fractional derivative used in our framework.

5468
SEUNG-YEAL HA, JINWOOK JUNG AND PETER KUCHLING Definition 2.1. [3,8,20] For a positive number α > 0, the Caputo derivative D α * f of order α is given by the following relation: as long as the R.H.S. are well defined, and Γ = Γ(z) is the Gamma function defined by the following integral: Next, we introduce the Mittag-Leffler function appearing as a solution of fractional differential equations with constant coefficients, and playing the same role of an exponential function for a linear ODE with constant coefficients. In the sequel, we list several basic properties of the Mittag-Leffler function and representation of the solutions to the fractional differential equations. We begin with definition of the Mittag-Leffler function.
Definition 2.2. [8,18,20] 1. For α, β ∈ C, we define the Mittag-Leffler function E α,β as follows: 2. A function f defined on (0, +∞) is completely monotone (CM), if it is a C ∞ -function and satisfies the following property: In the sequel, we provide some relations satisfied by the Mittag-Leffler function.
Then for an arbitrary integer p ≥ 1, the following assertions hold: 1. For |arg(z)| ≤ µ, Remark 2.1. As a direct application of the second statement in Proposition 2.2, we have where C is a positive constant. This can be seen easily as follows. Since |arg(−Ct α )| = π, we can use (2) of Proposition 2.2 to get the desired estimate (3).
The following proposition provides a situation where the Mittag-Leffler function satisfies the property so called "complete monotonicity (CM)" in Definition 2.2. Note that any CM function is nonnegative, monotonically decreasing and convex (see [18]). Throughout this paper, the monotonicity and nonnegativity of a CM function will be used repeatedly. Finally, we consider a general formulation of a solution for the following system of fractional differential equations: for T > 0, where Y : [0, T ) → R n , B : [0, T ) → R n and A is an n × n real constant matrix.

Proposition 2.4. [2]
Let Y be a solution to (4). Suppose that B belongs to a space is given by Then, the unique solution to (4) can be represented by where E α,β (A) is given by 2.2. The fractional Cucker-Smale model. In this subsection, we present the first moment estimate for the fractional C-S model.
Proposition 2.5. Let (X, V ) be a solution to (2). Then for any t ≥ 0, we have Proof. (i) We use the exchange symmetry i ↔ j to obtain the first relation: We use the representation formula (5) with to get the desired result.
(ii) We sum (2) 1 with respect to i to get Again, we apply the representation formula (5) with and use the result of (i) to obtain This yields the desired result.
Before we close this section, we recall the previous result for the mono-cluster flocking dynamics of (1). For a given configuration (X, V ), we set 12,14] Let (X, V ) be a solution to (1) with the initial data (X 0 , V 0 ) satisfying the following conditions: Then, there exists a positive constant x M > 0 such that 3. Mono-cluster flocking dynamics. In this section, we address the emergent behavior of the fractional Cucker-Smale model (2) with two types of communication weights, namely "time-independent weight" and "metric dependent weight". For the former case, we obtain exact decay rate of the velocity fluctuations toward the flocking velocity, whereas for the latter case, we provide such a decay estimate in a restricted manner due to the limitations in our analyses.
3.1. Time-independent weights. Consider a time-independent communication weight: For the moment, we assume that d = 1. Then, system (2) can be written as a vector form: where the coefficient matrix Ψ = (Ψ ij ) is given by By Proposition 2.4, the exact solution to (7) can be given by where E α (Ψt α ) is the Mittag-Leffler matrix-valued function defined in (6). We set From this representation and conservation of total momentum, we can easily find out that if 1 T V 0 = 0, then Thus, without loss of generality, we will assume that 1 T V 0 = 0. In the sequel, we have the following lemma.
Lemma 3.1. Suppose that the initial data and the communication weight satisfy and let (X, V ) be a solution to (7). Then, Proof. (i) We present upper and lower bound estimates for V (t) separately.
• Case A (upper bound estimate): We first look at the eigenvalues of Ψ in the space 1 ⊥ . Note that Ψ is symmetric, which makes Ψ orthogonally diagonalizable, and 1 is the eigenvector of Ψ corresponding to the eigenvalue 0 and by our assumption 1 T V 0 = 0, the velocity vector V (t) belongs to 1 ⊥ for all t > 0. Then, for any where ·, · denotes the standard inner product in R N and we used the zero row-sum property N j=1 Ψ ij = 0 (see (8)) and the symmetry of Ψ. From this, we can deduce that 0 is a simple eigenvalue of Ψ, and all eigenvalues except 0 are negative. Thus, the spectrum of Ψ is as follows: Hence, if we let b k be the unit eigenvector of Ψ corresponding to λ k , then we have From this decomposition, one has Therefore, we can obtain where we used that the matrix B is an orthogonal matrix satisfying the relation: This yields the desired result for the upper bound of (i).
• Case B (lower bound estimate): Note that for U = (u 1 , · · · , u N ) ∈ 1 ⊥ , Hence, we get This gives the desired result for the lower bound of (i).
(ii) It follows from the result (i) that we can obtain Here, it can be deduced from (i) and (ii) of Proposition 2.1 that Therefore, we use the above estimate to obtain This implies our desired result.
Moreover, we can also observe the emergence of flocking under a more general network, called a network with a spanning tree, as given in the following corollary: Suppose that the initial data and the communication weight satisfy and let (X, V ) be a solution to (7). Then, Proof. Since we can obtain the upper bound for X(t) and the lower bound of V (t) from the same argument in Lemma 3.1, we only prove the upper bound of Hence, we can find out that the spectrum of Ψ still consists of a simple eigenvalue, which is 0, and negative eigenvalues as follows: Since the rest of the proof can be obtained by following the proof of Lemma 3.1, here we omit the details and obtain the desired result. Now, we extend our argument to general d ≥ 1.

Theorem 3.2. Suppose that the initial data and the communication weight satisfy
and let (X, V ) be a solution to (7). Then, Then, we can rewrite the system (2) as If we set then it follows from the assumption that Now, we use Lemma 3.1 to each V k and X k to conclude the proof.
As a direct application of Theorem 3.2, we have the following corollary.
Corollary 3.2. Suppose that the following conditions hold: and let (x i , v i ) be a solution to (2). Then, ON THE FRACTIONAL CUCKER-SMALE MODEL 5475 3.2. Metric dependent weights. In this part, we consider a general case, i.e. the case where Ψ depends on the spatial variable X. As in the previous subsection, we first let d = 1, and then generalize the one-dimensional result to the multidimensional case. Note that the general representation of our system can be written as where Ψ(X(t)) is defined as Based on the proof of Theorem 2.6 in [21], we provide our desired result as below.
For the statement of the flocking result, let ψ = ψ(r) be a metric dependent communication weight, and for such ψ, we set We are now ready to provide our first main result on the emergence of monocluster flocking. Theorem 3.3. Suppose that the initial data and the communication weight satisfy the following relations: there exist ε ∈ (0, 1) such that where ψ * and X M are given by Let (X, V ) be a solution to (2) with d = 1. Then for any t ≥ 0, we have where λ := 2ψ * − ψ M is a positive constant. Before we provide the proof of Theorem 3.3, we set our framework for the proof. First, by the continuity of solution, there exists a δ > 0 such that where ψ * is the positive constant defined in Theorem 3.3. Now, we rewrite our system (9) as follows: where the matrixΨ is given by By Proposition 2.4, we can use the above relation to obtain the following system equivalent to (2): One objective for the proof of Theorem 3.3 is to show that the following set S ε is unbounded above: Note that the set S ε is nonempty by (10). To prove that S ε is unbounded above, we argue by a contradiction. So, we assume that T ε := sup S ε < ∞ and construct a sequence {V (j) } as follows: where U = V 0 . Here, we use the estimates in Lemma 3.1 to deduce thatΨ(X(t)) is positive semidefinite and the largest eigenvalue ofΨ(X(t)) can be bounded by For the estimates of {V (j) }, we provide a lemma.
With this setup, we now provide a proof of our main theorem.
Proof of Theorem 3.3. We use Lemma 3.4 to get where we used On the other hand, it follows from the assumptions on ψ and definition of ψ * that the ratio in the geometric sequence in (13) is strictly less than 1: where C is a positive constant. This implies that {V (j) } is a Cauchy sequence in the space of continuous functions on [0, T ε ]. Hence V (j) converges to a continuous function V * , and V * satisfies the second relation in (11). Thus, the uniqueness of solutions to fractional differential equations implies that V * is indeed a solution V to (9). Moreover, Since ψ M − 2ψ * = ψ M − (1 + ε)ψ 0 m < 0, the velocity alignment emerges at least in an algebraic rate. Here, we set λ := −(ψ M − 2ψ * ) > 0, and we can also obtain the boundedness of spatial variables as follows: On the other hand, it follows from the assumption on T ε that there exists a pair i = j such that However, we have This gives a contradiction, and hence, S ε is not bounded above. Thus, we can repeat all the process on [0, T ] for any T > 0 to yield the desired result.
The reason for this is as follows. Note that the above condition implies the existence of ε ∈ (0, 1) such that and for this ε, the other condition in Theorem 3.3 is automatically satisfied.

The result of Theorem 3.3 yields that
As a corollary of Theorem 3.3, we can also obtain a lower bound for the decay estimate of velocity fluctuation under more restricted conditions. Corollary 3.3. Suppose that the initial data and the communication weight satisfy the following relations: there exist ε ∈ 1 2 , 1 such that and let (X, V ) be a solution to (2) with d = 1. Then, we have Proof. Again, we construct a sequence {V (j) } as follows: for j ≥ 0. Then, it follows from the same argument in the proof of Theorem 3.3 that min This implies Ψ (X(t)) ≤ ψ M − ψ * , ∀ t ≥ 0.
Hence, we can follow the proof of Theorem 3.3 that for any T > 0, V (j) converges to V on [0, T ], which is a solution to (9), and the inequality on Lemma 3.4 holds on [0, T ]. Thus, one has where λ := 2ψ * − ψ M . On the other hand, it follows from Proposition 2.2-(2) that , as t → ∞.
Here, one has This implies 2ψ * −1 − λ −1 > 0 and we obtain our desired result. 4. Uniform 2 -stability estimates. In this section, we study uniform 2 -stability estimates for the solutions to system (2). Note that when the communication weights between particles are constant, it is easy to deduce the 2 -stability of solutions, since system (2) becomes a linear system in this case. Thus, we only cover the case with general communication weights. Here, we provide the detailed proof for the case when d = 1. Then, it is not difficult to deduce the results for the multi-dimensional case (d ≥ 2). Now, we state our second result on the 2 -stability.
Theorem 4.1. Suppose that the initial data (X 0 , V 0 ) and (X 0 ,Ṽ 0 ) satisfy the conditions in Theorem 3.3. Moreover, assume that there exists a constant η > 0 satisfying 2ψ Lip min{ V 0 , Ṽ 0 } < η(1 − η)λ 2 , where ψ Lip is a Lipschitz constant of ψ and λ is given in Theorem 3.3, and let (X, V ) and (X,Ṽ ) be two solutions to (2) with d = 1 corresponding to the initial data (X 0 , V 0 ) and (X 0 ,Ṽ 0 ), respectively. Then, there exists a positive constant G only depending on the initial configuration and ψ such that Proof. Without loss of generality, we may assume that V 0 ≤ Ṽ 0 , and set X := X −X and V := V −Ṽ .
Note that Φ(s)1 = 0 holds. First, we calculate Φ(s) for a later use. Note that the following estimate holds: We use the above relation to obtain that for U = (u 1 , · · · , u N ) ∈ 1 ⊥ , This yields Φ(s) ≤ 2ψ Lip X .
Hence, we use (14) to obtain Next, we construct a sequence {V (j) } of continuous functions iteratively: We may find that if the limit V := lim j→∞ V (j) exists, then V 0 = (1 + η) V 0 , and we can use a contradiction argument to yield

SEUNG-YEAL HA, JINWOOK JUNG AND PETER KUCHLING
Thus, it suffices to show that the limit exists, and it satisfies the desired inequality. Next, we claim: where C > 0 is a positive constant, and ρ is defined as Verification of (15): Similar to Lemma 3.4, we use induction on j to verify the relation (16).
• (Induction step): Suppose that (16) holds for j = n: Then, we use the above induction hypothesis to otain This verifies the claim for j = n + 1, and completes the induction procedure. Now, we use the same methodology in Theorem 3.3 to see that the sequence {V (j) } converges to the solution V on any finite interval [0, T ], and one has the following estimate: This implies the algebraic decay of V . Consequently, it yields the uniform boundedness of X as follows: Fianlly, we combine (17) with (18) to obtain

Now, we set
to get the desired estimate.
Remark 4.1. For a multi-dimensional case d ≥ 2, we again set V k andṼ k as in Theorem 3.2 and V k := V k −Ṽ k , and obtain Here, we use the convexity of 2 -norm to get the same inequality for V(t) . Then, we can follow the same argument in the proof of Theorem 4.1 to get the results for the multi-dimensional case.
5. Numercial simulations. In this section, we provide several numerical simulations to system (2). Note that one thing that distinguishes the fractional C-S model from the standard C-S model is the convergence mode toward velocity alignment. While the flocking emerges at the exponential rate in the standard case, it does at the algebraic rate in the fractional case. Thus, from numerical simulations, we would like to observe results not only about the velocity alignments, but also about the rate of velocity alignment. For numerical simulations, we used the scheme depicted in Appendix C in [8]. For readers interested in the scheme, we refer to Appendix A.  We present our results from the simulations in Figure 2-3. Here, we provide two figures for each case. The first one is a graph for V − v c over time, where v c is the mean velocity. This figure would show the asymptotic alignment of velocities of each particle. For the velocity alignment, we also consider the case α = 1 for comparison, which is actually the standard Cucker-Smale model. The second one is for over time, where d * y d * x denotes the following numerical differentiation:

Since our analytic result suggests
Thus, the graph of this quantity would show that the rate of velocity alignment is algebraic.   Here, we can observe that when ψ M does not exceed the twice of ψ m , the flocking emerges in an algebraic order (see Figure 2-(a) and 3-(a)). However, even for the case when ψ M exceeds the twice of ψ m , the flocking is still observed (see Figure  2-(b) and 3-(b)).

Nonnegative communication weights.
In this subsection, we present several results of numerical simulations, when the communication weight function ψ is just nonnegative. Here, we performed simulations under the following setup :  For the initial data, we take the following two types:  For the case (a), first we check whether the given data really satisfy our framework or not. One may obtain Thus, we can find out that ε = 0.8 satisfies our framework and hence our Theorem 3.3 guarantees that the flocking should emerge in this case. Now, we present our results for the simulations in Figure 5-6. Similar to the case where there is a positive lower bound, we present two figures for each case. In Figure  5-(a) and Figure 6-(a), we can observe that the flocking emerges asymptotically at the rate of t −α as expected. However, it can be also observed in Figure 5-(b) and Figure 6-(b) that the flocking emerges even though the flocking is not guaranteed.
However, if we choose an initial configuration where particles are far away from one another with relatively high velocity fluctuations, then it is possible to observe the case when flocking does not emerge. Results from numerical simulations for this situation is described in Figure 7. 6. Conclusion. In this paper, we have introduced the fractional Cucker-Smale model for the modeling of slow anomalous flocking, i.e., the velocity fluctuations around the asymptotic flocking velocity decays algebraically fast. For the proposed model, we provided a sufficient framework for the mono-cluster flocking and uniform 2 -stability estimate with respect to initial data. There are several interesting issues left for a future work. For example, the emergence of bi-cluster flocking, and meanfield limit of the fractional Cucker-Smale model as the number of particles tends to infinity. Furthermore, our numerical simulations show that, while the conditions in Theorems 3.3 are sufficient for flocking to emerge, they are not necessary. These issues will be treated in a future work.   (b) Although it is not guaranteed, the rate of velocity alignment is asymptotically t −α . Figure 6. Relaxation rate toward velocity alignment when ψ is just nonnegative.
Appendix A. Brief outline of the numerical scheme. In this appendix, we briefly explain the scheme that we used in Section 5 and for other schemes, we refer to [16]. Consider the following type of system of fractional differential equations: D α * y = f (t, y), y (j) (0) = y (j) 0 , j = 0, 1, · · · , [α].
For simplicity, we consider the uniform grid {t j = jh | j = 0, 1, · · · , n} with step size h on the time interval [0, T ]. Assume that we already have the values y j that approximate y(t j ) for j = 0, 1, · · · , k. Then in the scheme we will use, we would calculate the approximating value y k+1 from the following formula:  where a j,k+1 is given by if j = k + 1, and the corrector, y P k+1 , is given by