Exponential asymptotic flocking in the Cucker-Smale model with distributed reaction delays

We study a variant of the Cucker-Smale system with distributed reaction delays. Using backward-forward and stability estimates on the quadratic velocity fluctuations we derive sufficient conditions for asymptotic flocking of the solutions. The conditions are formulated in terms of moments of the delay distribution and they guarantee exponential decay of velocity fluctuations towards zero for large times. We demonstrate the applicability of our theory to particular delay distributions - exponential, uniform and linear. For the exponential distribution, the flocking condition can be resolved analytically, leading to an explicit formula. For the other two distributions, the satisfiability of the assumptions is investigated numerically.


Introduction
Individual-based models of collective behavior attracted the interest of researchers in several scientific disciplines. A particularly interesting aspect of the dynamics of multi-agent systems is the emergence of global self-organizing patterns, while individual agents typically interact only locally. This is observed in various types of systems -physical (e.g., spontaneous magnetization and crystal growth in classical physics), biological (e.g., flocking and swarming, [2,30]) or socio-economical [21,25]. The field of collective (swarm) intelligence also found many applications in engineering and robotics [16,19]. The newest developments in the mathematical approaches to the field are captured in, e.g., [1,26,4,24,27,28,29,13,20,10,23].
The Cucker-Smale model is a prototypical model of consensus seeking, or, in physical context, velocity alignment. Introduced in [8,9], it has been extensively studied in many variants, where the main point of interest is the asymptotic convergence of the (generalized) velocities towards a consensus value. In this paper we focus on a variant of the Cucker-Smale model with distributed delay. We consider N ∈ N autonomous agents described by their phase-space coordinates (x i (t), v i (t)) ∈ R 2d , i = 1, 2, · · · , N , t ≥ 0, where x i (t) ∈ R d , resp. v i (t) ∈ R d , are time-dependent position, resp. velocity, vectors of the i-th agent, and d ≥ 1 is the physical space dimension. The agents are subject to the following dynamicṡ for i = 1, 2, · · · , N , where | · | denotes the Euclidean distance in R d . The parameter λ > 0 is fixed and P is a probability measure on [0, ∞). For simplicity we consider constant initial datum on (−∞, 0] for the position and velocity trajectories, with (x 0 i , v 0 i ) ∈ R d × R d for i = 1, 2, · · · , N . The function ψ : [0, ∞) → (0, ∞) is a positive nonincreasing differentiable function that models the communication rate between two agents i, j, in dependence of their metric distance. For notational convenience, we shall denote ψ ij (t) := ψ(|x i (t) − x j (t)|).
In our paper we shall introduce the following three assumptions on ψ = ψ(r), namely, that ψ(r) ≤ 1 for all r ≥ 0, which clearly does not restrict the generality due to the freedom to choose the value of the parameter λ > 0. Moreover, we assume that there exist some γ < 1 and c, R > 0 such that and that there exists α > 0 such that ψ ′ (r) ≥ −αψ(r) for all r > 0. (6) The prototype rate considered by Cucker and Smale in [8,9] and many subsequent papers is of the form with the exponent β ≥ 0. The assumption (5) is verified for (7) if β < 1/2, while assumption (6) is satisfied for all β ≥ 0 by choosing α := 2β. Let us point out that the results of our paper are not restricted to the particular form (7) of the communication rate.
In real systems of interacting agents -animals, humans or robots, the agents typically react to the information perceived from their surroundings with positive processing or reaction delay, which might have a significant effect on their collective behavior. The system (1)-(2) represents a model for flocking or consensus dynamics where the reaction (or information processing) delay is distributed in time according to the probability distribution P . The main objective in the study of Cucker-Smale type models is their asymptotic behavior, in particular, the concept of conditional or unconditional flocking. In agreement with [8,9] and many subsequent papers, we say that the system exhibits flocking behavior if there is asymptotic alignment of velocities and the particle group stays uniformly bounded in time.
The term unconditional flocking refers to the case when flocking behavior takes place for all initial conditions, independently of the value of the parameters λ > 0 and N ∈ N. The celebrated result of Cucker and Smale [8,9] states that the system (1)-(2) without delay (this corresponds to the formal choice dP (s) := δ(s)ds, with δ the Dirac delta measure) with the communication rate (7) exhibits unconditional flocking if and only if β < 1/2. For β ≥ 1/2 the asymptotic behavior depends on the initial configuration and the particular value of the parameters λ > 0 and N ∈ N. In this case we speak about conditional flocking. The proof of Cucker and Smale (and its subsequent variants, see [15,14,3]) is based on a bootstrapping argument, estimating, in turn, the quadratic fluctuations of positions and velocities, and showing that the velocity fluctuations decay monotonically to zero as t → ∞.
The presence of delays in (1)-(2) introduces a major analytical difficulty. In contrast to the classical Cucker-Smale system (without delay), the quadratic velocity fluctuations are, in general, not decaying in time, and oscillations may appear. In fact, oscillations are a very typical phenomenon exhibited by solutions of differential equations or systems with delay, see, e.g., [12]. In [18] we developed an analytical approach for the Cucker-Smale model with lumped delay (corresponds to the formal choice dP (s) := δ(s − τ )ds in (1)-(2), with a fixed τ > 0). It is based on the following two-step procedure: first, construction of a Lyapunov functional, which provides global boundedness of the quadratic velocity fluctuations. Second, forward-backward estimates on appropriate quantities that give exponential decay of the velocity fluctuations. The main goal of this paper is to generalize the approach to the case of distributed delays with a general probability measure P . A demonstration of the approach to the scalar negative feedback equation with distributed delay, which can be seen as a special case of (1)-(2) with ψ ≡ 1 and N = 2, was recently given in [17].
Flocking in Cucker-Smale type models with fixed lumped delay and renormalized communication weights was recently studied in [22,5]. Both these papers consider the case where the delay in the velocity equation for the i-th agent is present only in the v j -terms for j = i. This allows for using convexity arguments to conclude a-priori uniform boundedness of the velocities. Such convexity arguments are not available for our system (1)- (2). In [6] the method is extended to the mean-field limit (N → ∞) of the model. In [7] the authors consider heterogeneous delays both in the x j and v j terms and they prove asymptotic flocking for small delays and the communication rate (7). A system with time-varying delays was studied in [27], under the a-priori assumption that the Fiedler number (smallest positive eigenvalue) of the communication matrix (ψ ij ) N i,j=1 is uniformly bounded away from zero. The same assumption is made in [11] for a Cucker-Smale type system with delay and multiplicative noise. The validity of this relatively strong assumption would typically be guaranteed by making the communication rates ψ ij a-priori bounded away from zero, which excludes the generic choice (7) for ψ. Our approach does not require such a-priori boundedness.
Cucker-Smale systems with distributed delays were studied in [23] and [29]. In both works, the delay is present in the expression for v j only, while v i in (2) is evaluated at the present time v i (t). The L ∞ analysis in [29] is based on a system of dissipative differential inequalities for the position and velocity diameters, leading to a nonexplicit "threshold on the time delay". The work [29] introduces hierarchical leadership to the distributed delay system. For the case of free will ultimate leader (i.e., it can change its velocity freely), a flocking result is given under a smallness condition on the leader's acceleration. To our best knowledge, the Cucker-Smale system of the form (1)- (2), where distributed delay is present in both the v j and v i terms on the right-hand side (2), has not been studied before. This paper is organized as follows. In Section 2 we formulate our assumptions and the main flocking result. In Section 3 we provide its proof divided into three steps -uniform bound on the velocities by a Lyapunov functional, forward-backward estimates, and exponential decay of the velocity fluctuation. Finally, in Section 4 we demonstrate the applicability of our theory to particular delay distributions -exponential, uniform and linear. For the exponential distribution, the flocking conditions can be resolved analytically, leading to an explicit formula. For the other two distributions, the satisfiability of the assumptions is tested numerically.

Main result
Let us first introduce several relevant quantities. For t ∈ R we define the quadratic fluctuation of the velocities, and the quantity Moreover, we introduce the moments of the probability measure P . The k-th order moment for k ∈ N shall be denoted M k , Finally, we shall need the moment K[κ], defined as We also introduce the quantity L (0), which depends on the (constant) initial datum and the second-and third-order moments of P , Our main result is the following: Theorem 1. Let the communication rate ψ = ψ(r) verify the assumptions (4)- (6). Let the parameter λ > 0 and the probability measure P be such that If there exists κ > 0 such that the conditions and are mutually satisfied, with α > 0 given by (6), then the solution of the system (1)-(2) subject to constant initial datum (3) exhibits flocking behavior in the sense of Definition 1. Moreover, the quadratic velocity fluctuation V = V (t) decays monotonically and exponentially to zero as t → ∞.
The above theorem deserves several comments. First, the conditions (12)- (14) relate the value of the parameter λ > 0, the moments of the probability measure P and the fluctuation of the initial datum through L (0). As we shall demonstrate in Section 4, for particular choices of the distribution P they lead to systems of nonlinear inequalities in terms of the distribution parameters and the fluctuation of the initial datum. These can be sometimes resolved analytically, leading to explicit flocking conditions. This is the case for the exponential distribution, as we shall demonstrate in Section 4.1. However, even if the nonlinear inequalities turn out to be prohibitively complex to be treated analytically, they are well approachable numerically. We show this for the uniform and linear distributions in Sections 4.2 and 4.3.
For fixed λ > 0 and P , the expression (14) can be interpreted as a smallness condition on the fluctuation of the initial datum mediated through the value of L (0). In fact, at the price of the flocking condition getting slightly more restrictive, the term L (0) in (11) can be replaced by a more intuitive expression. Indeed, with the bound ψ ≤ 1 given by (4) and the constantness of the initial datum, we have Consequently, we have and (14) is satisfied whenever We find this expression more appealing since it illustrates the necessity of smallness of the initial velocity fluctuation V (0) as a condition for asymptotic flocking. We admit that the assumption about the constantness of the initial datum on (−∞, 0], or on the interval corresponding to the support of the measure P , can be perceived as too restrictive. In fact, the methods we present in this paper can be generalized to the case of nonconstant initial data, as we demonstrated in [18]. However, since this would significantly increase the technicality of our exposition, we elect to focus on the essence of the method and thus restrict ourselves to the constant initial datum. We note that by the rescaling of time t → λt, of the velocities v i → λ −1 v i and of the probability measure P , the parameter λ is eliminated from the system (1)- (2). Nonetheless, for the purpose of compatibility with previous literature, we shall carry out our analysis for the original form (1)-(2). The scaling invariance shall become evident in Section 4, where we shall formulate the flocking conditions in terms of properly rescaled parameters of the probability distribution P and in terms of V (0)/λ 2 .
Finally, we note that the symmetry of the particle interactions ψ ij = ψ ji implies that the total momentum is conserved along the solutions of (2), i.e., Consequently, if the solution converges to an asymptotic velocity consensus, then its value is determined by the mean velocity of the initial datum.

Asymptotic flocking
The proof of asymptotic flocking for the system (1)-(2) will be carried out in three steps: First, in Section 3.1 we shall derive an uniform bound on the quadratic velocity fluctuation V = V (t) by constructing a suitable Lyapunov functional. Then, in Section 3.2 we prove a forward-backward estimate on the quantity D = D(t) defined in (9), which states that D = D(t) changes at most exponentially locally in time. Finally, in Section 3.3 we prove the asymptotic decay of the quadratic velocity fluctuation and boundedness of the spatial fluctuation and so conclude the proof of Theorem 1.

Lyapunov functional and uniform bound on velocity fluctuations
We first derive an estimate on the dissipation of the quadratic velocity fluctuation in terms of the quantity D = D(t) defined in (9).
Lemma 1. For any δ > 0 we have, along the solutions of (1)-(2), Proof. We have where the last equality is due to the conservation of momentum (16). With (2) we have For the first term of the right-hand side we apply the standard symmetrization trick (exchange of summation indices i ↔ j, noting the symmetry of ψ ij = ψ ji ), Therefore, we arrive at For the last term we use the Young inequality with δ > 0 and the bound ψ ≤ 1 by assumption (4), Next, we use (2) where [t − s] + := max{0, t − s} and we used the fact that the initial datum for the velocity trajectories is constant. Taking the square in (19) and summing over i we have The first inequality in (20) is Cauchy-Schwartz for the sum term, i.e.
We now define for t > 0 the functional where V = V (t) is the quadratic velocity fluctuation (8) and D = D(t) is defined in (9). Note that for t = 0 the above formula reduces to the formula (11) for L (0) (recall the constantness of the initial datum).  (17) we eliminate the integral term and obtain d dt We observe that the right-hand side is nonpositive if (22) is verified and conclude.
Consequently, if (22) holds, then the velocity fluctuation V = V (t) ≤ L (t) is uniformly bounded from above by L (0) for all t ≥ 0.

Remark 1.
Having established the decay estimate (23), one might attempt to apply Barbalat's lemma to prove the desired asymptotic consensus result, assuming merely the validity of (22). Indeed, with the uniform bound on velocities provided by Lemma 2 and the properties of the interaction rate ψ, one can prove that the second-order derivative d 2 dt 2 L (t) is uniformly bounded in time, which implies that d dt L (t) → 0 as t → ∞. This in turn gives D(t) → 0 as t → ∞. However, since ψ is not a priori bounded from below (and the uniform velocity bound allows for a linear in time expansion of the group in space), this does not imply that the velocity fluctuation V (t) decays asymptotically to zero.
Proof. For better legibility of the proof, let us adopt the notational convention that all quantities marked with a tilde are evaluated at time point t − s, i.e., v i : With this notation, we have and differentiation in time and triangle inequality gives, for t > 0, By assumption (6), |ψ ′ (r)| ≤ αψ(r) for r ≥ 0, we have for the first term of the right-hand side where in the second inequality we used the bound provided for t − s > 0 by Lemma 2. For t − s ≤ 0 it holds trivially due to the constantness of the initial datum.
For the second term of the right-hand side of (25) we apply the symmetrization trick, and estimate using the Cauchy-Schwartz inequality with some ε > 0 and the bound ψ ≤ 1 imposed by assumption (4), The first term of the right-hand side is equal to 2εD(t), while for the second term, for t − s > 0, we have with (2), the Jensen inequality and the bound ψ ≤ 1, For t − s < 0 we have d vi dt ≡ 0 due to the constant initial datum. Combining the above estimates in (25), we finally arrive at which immediately gives (24).
The following lemma constitutes the core of the forward-backward estimate method and was proved in [18,Lemma 3.5]. We reprint it here for the sake of the reader.
If there exists some κ > 0 such that then the following forward-backward estimate holds for all t > 0 and s > 0 e −κs y(t) < y(t − s) < e κs y(t).
Proof. Due to the assumed continuity of y(t) andẏ(t) on (0, ∞), (27) implies that there exists T > 0 such that We claim that (29) holds for all t ∈ R, i.e., T = ∞. For contradiction, assume that T < ∞, then again by continuity we have Integrating (29) on the time interval (T − s, T ) with s > 0 yields e −κs y(T ) < y(T − s) < e κs y(T ).
We now apply the result of Lemma 4 to derive a backward-forward estimate on the quantity D = D(t) defined in (9).

Decay of the velocity fluctuations and flocking
In order to bound D = D(t) from below by the quadratic velocity fluctuation V = V (t), we introduce the minimum interparticle interaction ϕ = ϕ(t), and the position diameter We then have the following estimate: Lemma 6. Let the parameter λ > 0 satisfy Then along the solutions of (1)-(2) we have Proof. Since, by assumption, ψ = ψ(r) is a nonincreasing function, we have with d X = d X (t) defined in (35). Moreover, we have for all i, j = 1, · · · , N , and Lemma 2 gives and integrating in time and taking the maximum over all i, j = 1, · · · , N yields which combined with (37) directly implies (36).
We are now in position to provide a proof of Theorem 1.
Proof. Let us recall the estimate (17) of Lemma 1, We apply (32) to the integral term Consequently, we have Optimizing in δ > 0 gives δ := λ K[κ], where we used the definition (10) of K[κ] := ∞ 0 s e κs −1 κ dP (s), so that By the definition (34) of the minimal interaction ϕ = ϕ(t) we have the estimate where for the last inequality we used (36) and the monotonicity of ψ, for all s > 0, and, consequently, Inserting into (38) yields Denoting ω := −2λ 2λ K[κ] − 1 > 0 and integrating in time, we arrive at Consequently, if ∞ ψ(s)ds = ∞, we have the asymptotic convergence of the velocity fluctuation to zero, lim t→∞ V (t) = 0. By assumption 5, namely that ψ(r) ≥ Cr −1+γ for all r > R, we have, asymptotically for large t > 0, Consequently, from (40), A slight modification of the proof of Lemma 6 gives The integral on the right-hand side is uniformly bounded, implying the uniform boundedness of the position diameter d X (t) ≤d X < +∞ for all t > 0. This in turn implies ϕ(t) ≥ ψ(d X (t)) ≥ ψ(d X ), so that (39) is replaced by the sharper estimate Thus we finally have, for all t > 0, and conclude the exponential decay of the velocity fluctuations.

Examples of delay distributions
In this section we demonstrate how the flocking conditions (12)- (14) of Theorem 1 are resolved for particular delay distributions -exponential, uniform on a compact interval and linear. The conditions (12)- (14) lead to systems of nonlinear inequalities in terms of the distribution parameters. For the exponential distribution they can be resolved analytically, leading to an explicit flocking condition. For the uniform and linear distributions they can be recast as nonlinear minimization problems and easily resolved numerically, using standard matlab procedures.

Exponential distribution
We first consider the exponential distribution dP (s) = µ −1 e −s/µ ds with mean µ > 0. We have for κ < µ −1 , Condition (12) reads 2 √ 2λµ ≤ 1, and (13) and (14) are satisfied if there exists κ > 0 such that As our goal is to visualize the dependence of the flocking condition on the initial velocity fluctuation V (0), we shall work with the more restrictive version of (14) given by (15), which for the exponential distribution reads Moreover, due to scaling properties, it is more convenient to investigate the flocking conditions in terms of the product λµ and rescale V (0) by λ 2 . In this form the flocking conditions read The second condition in (41) is easily resolved for κµ, κµ ≤ 1 − 2(λµ) 2 − 2λµ (λµ) 2 + 1, and the third condition is readily shown to be equivalent to 3 32(λµ) 2 + αλµ 2 1 + 12 Consequently, κµ can be eliminated from (41) and we arrive at Condition (15) reads Deciding satisfiability (in terms ofκ > 0) of the above conditions seems to be prohibitively complex for the analytical approach. However, the problem is well approachable numerically. For each pair (a, b) satisfying (43), the conditions (44)-(45) can be recast as a minimization problem inκ, and deciding satisfiability accounts to checking if the minimum is negative. The minimization problem can be efficiently solved using the matlab procedure fminbnd if we provide lower and upper bounds onκ. These can be obtained analytically. Indeed, carrying our Taylor expansion of the exponentials in (44) we see that Combining this estimate with (44) gives a necessary condition for its satisfiability in terms of explicit (in a and b) lower and upper bounds onκ, which are roots of the corresponding quadratic polynomial. We do not print the rather lengthy algebraic expressions here; let us just mention that an immediate rough lower bound isκ ≥ 2(a + b). We carried out two numerical studies. First, we fixed the values of α := 1 and V (0) λ 2 := 1 and plotted the critical value of the interval length (b − a) in dependence of the value of a > 0, see Fig. 2. We see that the flocking conditions (44)-(45) are satisfiable for a at most approx. 0.16, while for a approaching zero, the interval length can go up to approx. 0.26. In the second study, we fixed a := 0 and plotted critical value of the initial fluctuation V (0)/λ 2 in dependence on the interval length b > 0, Fig. 3.  Due to the scaling relations, it is again convenient to express the flocking conditions in terms of a := λA, κ := κλ −1 and V (0)/λ 2 . Condition (12) reads then a ≤ 3/2, while (13) and (15)  We are interested in the dependence of the critical value of the rescaled initial fluctuation V (0)/λ 2 on the parameter value a. We approach the above satisfiability problem numerically, in two steps. First, we observe that for any fixed a ∈ (0, 3/2), the function f a (κ) := 4 κ 2(e aκ + 1)

Linear distribution
is an increasing function of κ > 0; this is easily seen carrying out the Taylor expansion of the exponentials. Moreover, limκ →0+ f a (κ) = 2a 2 /3 < 1. Consequently, there existsκ a > 0 such that the first condition of (46) is equivalent toκ ∈ (0,κ a ). The value ofκ a is conveniently calculable using the matlab procedure fminsearch, profiting from the monotonicity of the function f a . In the second step, we numerically solve the maximization problem employing the matlab procedure fminbnd. Then the second condition of (46) is easily resolved for the critical value of V (0)/λ 2 . The outcome of this procedure for α := 1 is plotted in Fig. (4).