Nonlinear stability of flock solutions in second-order swarming models

In this paper we consider interacting particle systems which are frequently used to model collective behavior in animal swarms and other applications. We study the stability of orientationally aligned formations called flock solutions, one of the typical patterns emerging from such dynamics. We provide an analysis showing that the nonlinear stability of flocks in second-order models entirely depends on the linear stability of the first-order aggregation equation. Flocks are shown to be nonlinearly stable as a family of states under reasonable assumptions on the interaction potential. Furthermore, we numerically verify that commonly used potentials satisfy these hypotheses and investigate the nonlinear stability of flocks by an extensive case-study of uniform perturbations.


Introduction
Self-organization and pattern formation are ubiquitous in nature and science, ranging from animal aggregation [1][2][3][4][5] and biological systems [6,7] to self-assembly of nano-particles [8,9]. The intense studies during the last two decades are on both individual-based systems and continuum equations [10][11][12][13]. One of the essential features in these models is the non-zero characteristic speed of the individual agents, which has been modelled in different ways. The speed can be assumed to be constant with a direction based on the averages of the neighbours [14] or to be driven by random noise [15]. On the other hand, a large class of models consist of self-propelled particles powered by biological or chemical mechanisms with friction forces, resulting in a preferred characteristic speed.
In general, the particles with non-zero equilibrium speed do not form any recognizable patterns [4,16], and interactions within the group have to be included to generate interesting spatial configurations. Most of these interaction forces have been taken into account in the combination of three effects: alignment, repulsion, and attraction; also called the ''first principles of swarming''. The basic mechanisms account for collisional avoidance and comfort regions (repulsion), grouping and socialization (attraction), and mimetic synchronization (alignment). The combination of these three effects goes back to fish school modelling [17,18] and these basic ideas have been improved and applied to several animal species including different mechanisms and interactions more adapted to particular living organisms, see for instance [19][20][21][22][23].
In this paper, we focus on the cases with velocity-independent interactions. More precisely, by introducing a pairwise symmetric interaction potential W (x) = U(|x|), we consider the two-dimensional model where x i , v i ∈ R 2 , i = 1, . . . , N are the positions and velocities of the individual particles and α, β are effective values for self-propulsion and friction forces, see [10][11][12] for more discussions. For this relatively simple system, a variety of patterns are observed, for instance coherent moving flocks and rotating mills [10]. Other patterns like rings and clumps were discovered in [11,12], with the introduction of the concept of Hstability from statistical physics to better characterize the patterns in the parameter space. Delta rings, with mass uniformly distributed on a circle, have been further studied thoroughly for power-law like potentials, we refer to [24][25][26][27][28][29][30] for details. While the stability and bifurcation of the ring solutions can be investigated in a straightforward way because of their explicit particle representations on a circle, there are few studies on the more prevalent compact steady solutions, like flocks or mills. The reason lies in the difficulties to solve some complicated integro-differential equations for most of the popular choices of the potential W . Only for certain particular potentials like the quasi-Morse type proposed in [31], the existence and uniqueness of the coherent moving flocks can be established rigorously [32].
In this paper, we focus on the stability of flock solutions of (1) for general potentials defined below.

Definition 1.
A flock solution of the particle model (1) is a spatial configurationx with zero net interaction force on every particle, that translates at a uniform velocity m 0 ∈ R 2 with |m 0 | =  α β , hence (x i (t), v i (t)) = (x i − tm 0 , m 0 ).
We note, that the spatial configurationx is a stationary state to the first-order interacting particle system which is analyzed for example in [8,9]. In this work, we focus on general flock solutions whose spatial configurations tend to a compactly supported particle density in the continuum limit. One conclusion from our analysis is the somehow unexpected deep relation between the linear stability of the flock spatial configuration as steady state for the first-order swarming model (2) and the nonlinear stability of the family of associated flock solutions for the second-order swarming model (1). This relation was already found in the linear stability analysis around flock solutions in [30]. There the authors showed that the linearization of (2) around the equilibrium statê x has a positive eigenvalue if and only if the linearization of (1) around the steady flock solution in the co-moving frame has an eigenvalue with positive real part. Moreover, they show that if the equilibrium statex is linearly stable for (2), then the associated flock solution is always linearly unstable due to the presence of a generalized eigenvector associated to the zero eigenvalue of the linearization of the flock solution in the co-moving frame due to symmetries.
However, it is in this work where we clarify completely this deep imbrication by showing the nonlinear stability of the flock solutions under mild hypotheses on the linearization of (2) about the spatial configurationx. Actually, the basic hypotheses, apart from few technical assumptions, are that the steady statex of (2) is linearly asymptotically stable except the obvious symmetries: translations and rotations. Therefore, we can apply our theorem to verify the nonlinear stability of flock solutions to more biologically adapted potentials such as the Morse and the Quasi-Morse potentials [11].
The main result of this paper asserts, except few technicalities, that the family of flock solutions associated to a linearly asymptotically stable steady state of (2) is asymptotically stable for the dynamics of (1). Here, the asymptotic stability of the family of flock solutions means that any small enough perturbation in (x, v)-space at any time t 0 will, under the dynamics of the system (1), relax towards (likely) another flock solution in the family at an exponential rate as t → ∞. Let us finally emphasize that the most rigorous way of stating our main theorem uses advanced concepts of dynamical systems [33]. Our main theorem can be rephrased as follows: the family of flock solutions to (1) associated to a linearly (except symmetries) asymptotically stable steady state of (2) forms a normally hyperbolic invariant manifold for the system (1) with an empty unstable manifold.
This result has important applications, especially in the study of flock patterns of (1) using particle simulations. In general, the desired patterns are observed only after certain transit dynamics, and with carefully prepared initial data. With the spatial configurations from the first order system (2) (if they exist), the computational load is reduced significantly, making possible to study the existence of these patterns in great details.
After some preliminary notations and assumptions, the main result of this paper is stated in Section 2. The formulation of the linearized dynamics when the mean velocity is transformed to zero is given in Section 3, where the non-trivial relation between the Jacobian matrices of the two systems (1) and (2) is clarified. A comparison of this work to the techniques used in [30] also is given. In Section 4, we use the eigenspace structures of these Jacobian matrices to reveal the connection between the nonlinear stability of the flock solutions for (1) and the linear stability of their spatial configuration for (2), leading to the proof of the main theorem. Section 5 is devoted to numerical experiments to check the validity of the assumptions on the linear stability for the steady state of the system (2) for different biologically reasonable potentials.

Main result
In this section, we formulate the main theorem of this work. We will start with some assumptions on the spatial configuration of flock solutions, which is a stationary state of the first-order swarming model (2). The hypotheses assume this spatial configuration to be stable under the dynamics of (2). Given a steady spatial particle configuration to (2), i.e, a set of particle positions ( we denote byx the 2N-dimensional spatial configuration vector (x 1 , . . . ,x N ). We can now write the linearization of the system (2) at the steady configurationx as where G(x) denotes the Jacobian matrix associated to (2) and h ∈ R 2N . The 2N × 2N-Jacobian matrix is explicitly given by G(x) = (G ij ) with G ij being the 2 × 2-blocks defined as with Hess W denoting the Hessian matrix of W .
We will make extensive use of the following standard notation: given two matrices A ∈ R n×m and B ∈ R p×q , the Kronecker product A ⊗ B is defined as the matrix A ⊗ B = (a ij B) ij ∈ R np×mq . Let us also denote by 0 n , 1 n the column vectors of length n with all entries equal to 0 or 1 respectively, and I n the identity matrix of size n. Let us now specify the precise assumptions onx, and thus indirectly on the potential W .

Assumption 1.
Given an interaction potential W , we assume the existence of a linearly asymptotically stable (except symmetries) stationary state not lying on a line, i.e., a spatial configurationx such that the following assumptions hold: • (H1)x is a stationary state of (2), i.e.x satisfies (3). • (H2) The eigenspace for the zero eigenvalue of G(x) is spanned by the eigenvectors representing invariance of G(x) with respect to translations and rotations in R 2 , and therefore is three dimensional, dim(Eig(G(x), 0)) = 3. • (H3) All other eigenvalues of G(x) have negative real parts, i.e. Re(λ i ) < 0 ⇔ λ i ̸ = 0 for all i ∈ σ (G(x)). • (H4) Not all points ofx lie on a straight line in R 2 .
Hypotheses (H1)-(H3) imply thatx is linearly asymptotically stable for (4) except for the obvious symmetries, or equivalently that the Jacobian matrix G(x) is negative definite except for the null eigenspace in (6). As a consequence,x is nonlinearly asymptotically stable for the dynamics of (2) except for translations and rotations. More rigorously speaking and making use of more advanced dynamical systems theory [33], properties (H1)-(H3) imply that the family of stationary states of (2) obtained by rotations and translations fromx given by forms a normally hyperbolic invariant manifold for the system of ODEs in (2) for which the unstable manifold is empty. Here, is the 2 × 2 rotation matrix of angle φ. This stability concept will be further discussed and used later on for the second order model (1). Let us remark, that assumptions (H1)-(H4) are very mild in the sense that they are natural properties of stationary states of (2) found through particle simulations from generic initial configurations, e.g. compactly supported flocks. Particular special cases satisfying hypotheses (H1)-(H4) are the so-called Delta rings: particles arranged equidistantly on a circle. We will discuss in the last section more biologically relevant potentials, such as Morse-like potentials [8][9][10][11]31,32], for which we can numerically find stationary statesx satisfying hypotheses (H1)-(H4). Now, let us make precise the concept of the family of flock solutions of (1) associated to the stationary statex of (2).
Let us denote by V s the set of possible steady velocity vectors for a flock solution of speed s =  α β to (1) defined as Then the family of flock solutions of (1) associated to the steady statex of (2) is given by It is straightforward to check that all z ∈ Z F are solutions to (1) in vectorial form. The next section will be devoted to analyse the system (1) in suitable coordinates to study their stability. A change of variables to (1) that eliminates the translation of one particular member of the flock family is given by the co-moving frameṽ i := v i − m 0 transforming the flock (x * + v * t, v * ) onto a stationary state. However, it is very natural to think of small perturbations on (x * + tv * , v * ) that, under the dynamics of (1), tend to possibly another spatial configuration in RT (x) translating at a different velocity generated bym 0 , i.e. to another member of the family of flocks Z F . Hence, we cannot expect linear stability for general perturbations for such a change of variables. . . . Instead, we propose a new change of variables measuring the difference to mean velocity in velocity space, and aim to establish nonlinear stability in the entire family Z F .

Definition 2.
Defining the mean velocity of N particles as we introduce the change of variables with respect to the mean velocity m(t) as: We will work with the dynamical system (1) and its linearization in vector form, and thus we need to introduce some matrix notation used in the following.

Definition 3 (Notation).
1. 0 n×p is a block of zeros of size n × p.
2. Writing δ i n ∈ R n , we refer to a vector of length n with value 1 as its i-th entry and zeros elsewhere.
3. For the sake of simplicity, we occasionally write column vectors as row vectors neglecting the transpose operation on its components, e.g.
4. For a square matrix A ∈ R n×p or a column vector b ∈ R n , we denote by the matrix obtained from A by eliminating the last two rows and analogously for vectors. The k-th row or column of A is labelled as a k,· or a ·,k respectively.

Lemma 1.
Using the change of coordinates (8), the microscopic model (1) reads For any mean velocity m 0 ∈ V s , the vectorQ = (x, 0 2N , m 0 ) T is a stationary state of (9). Moreover, the Jacobian matrix F obtained by linearizing (9) around a stationary stateQ with an arbitrary orientation m 0 is given by Proof. By definition using that the net total force in the system is zero by symmetry, i.e. Newton's action-reaction principle. From this, the equation for is a stationary state of (9) since G(x) = 0 and all self-propulsion terms vanish. Let us drop the tildes henceforth for simplicity and let us compute the Jacobian matrix of the system (9). The first line of blocks in F is trivial. In the second line of blocks, we first have G(x) since the dependency on x in (9) is identical to (1), hence one obtains the Jacobian of the first-order particle system (2) evaluated at the stationary spatial configuration. To compute the second block, consider an arbitrary particle i and its velocity vector Then, identifying the right-hand sides of (9) with the time derivative symbols, we get: for j ̸ = i and analogous expressions for dt . Evaluating atQ the last expression in all terms vanishes since |m 0 | 2 = α β , which in total gives the structure of the block using the Kronecker product notation. In the third line of F , the second block follows analogously. For the last block, we compute which reduces to ∂ m 1 dm 1 dt = −2βm 2 1 atQ , and analogously for the remaining partial derivatives. This completes the computation of the structure of F .
Let us point out that the family of flock solutions to (1) is translated through the change of variables (8) to the set of stationary solutions of (9) given bỹ Now, we have almost all the ingredients to write the main result of this paper that holds under an additional technical assumption on the interplay between the linearization and the orientation of the flock. Assumption 2 (H5). All eigenvectors to non-zero eigenvalues of G(x * ) are not particle-wise orthogonal to the fixed meanvelocity m 0 . This means, that for pairs of entries 2i − 1, 2i of any eigenvector w of G(x * ), we have We will understand the role of this additional assumption below, but let us point out that this assumption is generically met by the members of the flock familyZ F . The main result of this work reads as: Theorem 1. Assume that the symmetric pairwise interaction potential W allows for a stationary spatial configurationx of (2) satisfying assumptions (H1)-(H4). Consider the second-order swarming model (1) with α, β > 0, then the family of flock solutions Z F of (1) defined in (7) is locally asymptotically stable for all configurations x * satisfying (H5) in the following sense: Any small enough perturbation in (x, v)-space at any time t 0 of the flock solution z * ∈ Z F associated to x * will, under the dynamics of the system (1), relax towards another flock solutionz ∈ Z F at an exponential rate as t → ∞.
Let us remind that assumptions (H1)-(H4) imply in particular that the family of stationary states RT (x) forms a normally hyperbolic invariant manifold for (2) with empty unstable manifold. Therefore, our main theorem can be rephrased as stating that the family of flock solutions associated tox of the second order model (1) has essentially the same property. More precisely, we have the following result.

Corollary 1.
Assume that the symmetric pairwise interaction potential W allows for a stationary spatial configurationx of (2) satisfying assumptions (H1)-(H4). Then the subset of stationary states ofZ F in (11) satisfying (H5) is a normally hyperbolic invariant manifold of (9) with empty unstable manifold.
Next two sections are devoted to the proof of these results. We start in next section by restricting the dynamics to a physically consistent subspace.

Physically consistent dynamics and its linearization
Note that (9) is a 4N + 2 dimensional system and thus F ∈ R (4N+2)×(4N+2) by having added the mean velocity as an additional variable. However, a stability analysis of (10) would include unphysical perturbations where particle velocities do not match mean velocity. The dynamics we are interested in are 4N-dimensional, and thus we need to reduce (10) to physically consistent dynamics. Clearly, mean velocity consistent states are invariant under the dynamics of (9), since We hence define a base for the 4N-dimensional subspace of mean velocity consistent states by removing the redundant equation of the N-th particle.
Any physically relevant perturbation of (1) is projected onto span(B) via P : Clearly, B is isomorphic to a base of R 4N and the relevant invariant subspace for the dynamics of (9). As mean velocity consistency is a linear algebraic constraint to (9), it turns out that span(B) is also an invariant subspace of F , which means we can restrict the stability analysis of (10) to physically relevant perturbations by computing the matrix representation of where ⌈·⌉ is used according to Definition 3.4.

Proof.
We compute the image of base vectors b i under F . First, consider b i with a non-zero entry in a spatial coordinate (i = 1, . . . , 2N). Then Fb i = (0 2N , g ·,i , 0 2 ) T consists of the i-th column of G(x) = (g i,j ) ij . Since the sum of each column of (5) is zero, we can rewrite this in terms of B as which gives the first block column of (12). Next, consider b i with a non-zero entry in a velocity component k and assume i to be odd (i = 2N + (2k − 1), k = 1, . . . , N − 1). Then we can compute Fb i to get and zero in all other rows, which is then multiplied with −2β(m 0 ⊗ m T 0 ) to select its first column twice (as i is odd), indicated by the subscript () ·,1 . This rewrites in terms of B as We now derive two corollaries on the linear structure of F B B , which will be essential to establish our stability result. To do so, we first rewrite (12) as since the two blocks separate. We first investigate the eigenspace associated to the zero eigenvalue of F B B . . . .

Lemma 3. Suppose that G(x) satisfies (H1)-(H4), then F B B does not possess a generalized eigenvector for eigenvalue zero.
Proof. Since the lower 2 × 2 block of F B B in (13) can be diagonalized with eigenvalues 0 and −2α, then we restrict our analysis to H. To simplify the notation, we will skip the dependency of the matrix G(x) on the stationary state in this proof.
If z = (x, v) T ∈ Eig(H, 0) is a zero eigenvector, then v = 0 2N−2 and x ∈ Eig(G, 0) by (12). Suppose (u, w) T is a generalized eigenvector, then which is equivalent to Due to the definition of G in (5), we get Adding two rows using (15), (14a), and (14b), (14c) is equivalent to We left-multiply with x and get since G symmetric and x ∈ Eig(G, 0). This implies x = [a 1 . . . a N ] T ⊗m ⊥ for some constants a i since the right-sided expression is quadratic and hence x has to be particle-wise orthogonal to m 0 for the equation to hold (m ⊥ 0 := (−m 0,2 , m 0,1 ) T ). On the other hand, (14b) implies that x ∈ span(w 3 ), since w 3 , as given in (6), is the only zero eigenvector of G whose component sums can be zero. Hence, all position vectors are tangential to a fixed mean velocity m 0 , which is a contradiction to (H4).
. . . Next, we investigate the eigenvalues of F B B under another hypothesis linking the eigenvectors of G(x) with the orientationally arbitrary but fixed mean velocity m 0 .
Eq. (18c) rewrites as By the same argument used in Lemma 3,(19) is equivalent to due to (18b) and the properties of G. We left-multiply withx to obtain assuming without loss of generality the normalization ⟨x, x⟩ =x T x = 1 and denoting the complex scalar product as ⟨·, ·⟩.
Therefore, we deduce Note that A ≥ 0, B ≥ 0 sincex T Gx ≤ 0 for any x due to (H3). We study three cases: its square root is purely imaginary and thus Re(µ i ) < 0. Case 2: Suppose A = 0, B > 0. Then µ = ± √x T Gx ∈ iR. Since A is a sum over squared moduli of N complex numbers, Since m is fixed, this holds if and only if x is particle-wise orthogonal to m 0 . Now, inserting A = 0 back into (20) gives Gx = µ 2 x, µ 2 < 0 which implies that such an x is also an eigenvector for G. This is a contradiction to (H5) and no such case exists.
Together, this shows that To complete the study of Eig(F B B , 0) suppose now that µ = 0. Then A 2 −B ≥ 0 and 0 = −A+ , which takes us back (22b). Therefore, no further eigenvector to eigenvalue zero exists, which completes the proof.
Let us remark that the 4 basis eigenvectors associated to the zero eigenvalue of F B B represent linearized flow within the set of stationary flock solutionsZ F in (11). From the change of variables (8) and Lemma 1, we know that the family of flock is invariant under translation and spatial rotation, and thus z 1 , z 2 , and z 3 are the linearized flow along these invariances. The fourth eigenvector z 4 represents linearized rotation of the mean velocity m 0 .. . . [30]). The decisive technical difference used in the analysis presented here compared to [30] is the change of variables with respect to mean velocity. Here, we are able to show the absence of a generalized eigenvector in Lemmas 3 and 4. In the co-moving frame used in [30], this is not possible as a particular mean velocity is fixed, therefore the instability of a flock solution in the co-moving frame is shown. Instead, the instability criterion are in fact established proving the existence of a generalized eigenvectors associated to the zero eigenvalue (see [30,Remark 3.6]). The difference of the two approaches can also be seen in the proof of Lemma 4, which bears some technical similarity with [30, Lemma 3.1], but additionally uses the mean-velocity consistency in (18b) and (19), (20).

Proof of the main result
Changing coordinates to (8) eliminates the linear translation at constant speed and flock solutions become constant solutions which we identify with their profile and velocity at t = 0. The image of Z F under (8) isZ F in (11).Z F is a manifold of stationary states of the flow of system (9). We aim to show thatZ F is a normally hyperbolic invariant manifold according to [33] for the time flow map of (9) restricted on span(B), which we denote by F t B . First,Z F is trivially invariant, F t B (Z F ) =Z F for all t ∈ R, as it is entirely stationary by Lemma 1. Choose Q * ∈Z F such that (H5) is fulfilled. Since eigenvectors are distinct, (H5) holds true for a neighbourhood of Q * . There is no tangential flow alongZ F , and in particular Lemmas 3 and 4 tell us, that F B B decomposes into the direct sum of Eig(F B B , 0), which equals the tangential space ofZ F at Q * , and a subspace of eigenvectors whose eigenvalues have negative real part and span the non-tangential subspace ofZ F at Q * . It is important to stress, that this splitting is preserved under the flow F t B , since Q * is stationary. Adapting from [33], all prerequisites forZ F to be a normally hyperbolic invariant manifold are given. In the language of dynamical systems theory, the tangent bundle of R 4N at Q * has an F B -invariant splitting into tangential and stable sub-bundles. A necessary inequality between the tangential flow and the stable decay rate is trivially given in our case, since flow alongZ F is zero and as a consequence, we have a spectral gap. Note, that thoughZ F is not compact we can easily select a compact sub-manifold around Q * or perform a symmetry reduction and eliminate the translational invariance by further reducing span(B).
We can hence apply the theorem of Hadamard and Perron as in [33,Theorem 4.1], which, adapted to our situation, gives the following: There exists locally an F B -invariant stable smooth manifold in a neighbourhood ofZ F , where the splitting into stable and stationary of the linearized flow is preserved. The space between both manifolds is laminated by invariant fibres tangent to the stable subspace at Q * , which are curves of the nonlinear flow towardsZ F . Hence, for a small perturbation of Q * within that laminated space, we are sure to be sitting on some fibre, which might have Q * as its base point onZ F or a point in a neighbourhood of Q * . In either case, there is transport along the fibre towardsZ F , for which an exponential decay estimate holds. This estimate is obtained by closing (23) on the stable sub-bundle. In other words, there exists locally a decomposition into fast and slow (stationary, in our case) coordinates of the nonlinear flow, such that the dynamics is characterized by the linearized flow. Going back to original coordinates (1), we obtain the statements in Theorem 1 and Corollary 1.

Remark 2.
The precise meaning of ''relaxation in time'' towards a different flock solution in the statement of Theorem 1 is given by the notion of stable hyperbolicity of the family of flock solutions in coordinates (8), as used in the proof. In the original coordinates, there is no notion of ''closeness'' between elements of (7), since e.g. a change of orientation will cause two flocks to arbitrarily diverge in euclidean norm as t → ±∞. Relaxation thus means, that the trajectory of a perturbed flock will approach a trajectory of another flock solution with typically different orientation and similar spatial configuration. Theorem 1 also tells us that the manifold is stable towards small perturbations of the flow F B , though that case is not considered here. Note, that this allows for flock stability under small perturbations of the potential. Table 1 Numerical results for the largest non-zero eigenvalue µ 4 of the Jacobian G({x}) for varying potentials W and number of particles N. The numerical value of the zero eigenvalue µ 3 is shown as a measure of the exactness of {x}. The free scaling parameter D is set to normalize the smallest eigenvalue to min(µ i ) = −1. As µ 4 is negative, hypotheses (H2)-(H3) of Theorem 1 are numerically validated for the cases shown.

Numerical experiments
In this section, we numerically investigate the hypothesis of Theorem 1 for a number of commonly used interaction potentials and study the nonlinear stability of flocks in terms of their polarization.
The key advantage of Theorem 1 is that essentially linear (asymptotic except symmetries) stability of the spatial configuration of a flock in the first-order aggregation model (2) is inherited nonlinearly by the family of flock solutions to the second-order model (1). In order to demonstrate that the stability theorem applies to a given potential, (H1)-(H4) have to hold forx being a stationary state of (2). We compute a numerical approximation {x} ofx by a standard fourth-order Runge-Kutta time integration of (2) with random initial data and an interaction potential scaled with a factor D to relate the amplitude of the right-hand side of (2) with the accuracy of the solver: Note, thatx is independent of D. We hence obtain a spatial configuration for which (H1) holds numerically. Let us mention that (H4) is clearly satisfied for all potentials considered here. Hypotheses (H2)-(H3) require the Jacobian G({x}) to possess negative eigenvalues except for three zero eigenvalues associated to translation and rotation. As G({x}) ∈ R 2N×2N is symmetric it is sufficient to apply standard methods to compute its spectrum, at least for the number of particles we consider.
Let (µ 1 , µ 2 ; µ 3 ) denote the eigenvalues of G({x}) associated to translation and rotation and µ i , i ≥ 4 all other eigenvalues sorted in descending order. In Table 1 we show µ 4 for different potentials and numbers of particles, together with |µ 3 | to illustrate the numerical accuracy of {x} and the scaling parameter D. For the comparability of the potentials, D is set such that µ 2N = 1 for all cases shown. The potentials used in Table 1 are: where K 0 is the modified Bessel function of second kind (for detailed discussions of the respective potentials see [32,31,12,34]).
It can be seen that for most cases considered, eigenvalues other than µ 1 , µ 2 , µ 3 are negative and hence Theorem 1 applies. The difference between the numerical eigenvalue µ 3 , which is equal to zero in theory, and µ 4 is of order four or higher, which shows that the negativity of µ 4 is not a numerical artefact. However, in some rare particular cases this conclusion cannot be drawn as µ 3 and µ 4 are of the same order (see the case N = 40, Log-Newtonian). We remark, that the eigenvalue computations of Table 1 require a highly accurate computation of {x}, which is here obtained with a scaling of D = 500 and a time interval T = [0, 500].
Next, we aim to study the stability of a flock in the nonlinear dynamics of (1) for the Morse potential under uniformly distributed perturbations in space and velocity. Theorem 1 provides a theoretical bound on the decay rate into equilibrium, which is very small as seen in Table 1. For random perturbations it is natural to assume that the flock is ''more stable'' than the theorem asserts. To study whether a flock remains a flock after an initial perturbation and the subsequent dynamics of (a) Statistical comparison of initial polarization vs. minimal polarization in time evolution.
(1), we use the notion of polarization of a group of individuals, which is defined as for N agents with phase-space coordinates at (x i , v i ) ∈ R 4 , see e.g. [2] for instance. For an aligned flock, the polarization equals one whereas e.g. a rotating mill possesses zero polarization. We perform a sequence of simulations of (1) with randomly perturbed initial data. The stationary flock solution and varying perturbation strength a. Stability of the flock in the time evolution [0, T ] can estimated as follows: If the support of the solution remains bounded and the polarization is equal or larger than the change in polarization induced by the initial perturbation, is it reasonable to assume that the ensemble will again tend to form a flock solution. If however, the polarization drops to significantly low or near zero values, the particles have evolved into an unpredictable state whose convergence behaviour as T → ∞ is unknown. In Fig. 1(a), we plot the initial polarization induced by the perturbation against the minimal perturbation that occurred during a time integration of (1) with T = 100. The figure is based on 25,000 simulations with uniformly distributed perturbations strength a and shows the statistical mean and the 5% tails of the statistical distribution. In Fig. 1(b), we show the correlation between the polarization induced by the initial perturbation and the averaged l 2 -norm of the perturbation per particle 1 N  i ∥(∆x i , ∆v i )∥ 2 , using the same statistical measures. The potential parameters are identical to Table 1 and the number of particles is N = 100. We observe from Fig. 1, that the probability of loosing additional polarization in the flock nonlinearly increases with the strength of the initial perturbation. Below a critical threshold of approximately 0.6, the flock will most likely drop polarization with an increasing probability of near zero values (e.g. unpredictable state). From Fig. 1(b) we observe that a perturbation-induced polarization of 0.6 is on average generated by a strong average perturbation in (x, v) of approximately 1.34 per particle, which on average is a perturbation of approximately 0.95 in position per particle. However, the mean nearest neighbour distance of the stationary state considered here is approximately 0.252. Hence the nonlinear stability observed in this experiment well exceeds the validity regime of Theorem 1, which is expected to hold for spatial perturbations that are small compared to the nearest neighbour distances.

Conclusions & perspectives
In this paper, we have shown that flock solutions are fully nonlinearly stable under small enough position-velocity perturbations, in the sense that a perturbed flock configuration will relax to another flock configuration inZ F . Stability is inherited from the first-order particle system to the second-order model. Hence the flock might slightly rotate its position, translate its centre of mass, change the orientation of its mean velocity, or any combination of the three, but it will remain a flock solution. This result has been established for general interaction potentials W satisfying assumptions (H1)-(H4). We numerically investigated the spectrum of the Jacobian G(x) evaluated at a stationary statex for some commonly used potentials, and found them to possess the desired eigenvalue structure. Flocks for Morse-type potentials are linearly unstable in the co-moving frame second-order model [30], but nonlinearly asymptotically stable in the sense of our main result in Theorem 1. Finally, we investigated the nonlinear stability of flocks under uniformly distributed perturbations, and numerically validated the nonlinear stability of flocks in the sense of polarization.
Our findings are completely consistent with one's expectations from numerical experiments and, to the best of our knowledge, is the first of its kind for general flocks in second-order models. The task of theoretically proving nonlinear stability properties such as (H1)-(H4) of potentials and their steady solutions for the first-order aggregation equation remains largely open, and is inevitably complex given the extensive variety of stationary states (see e.g. [28]). An extension of the obtained result to three or higher dimensional space seems achievable, however with greater technical effort due to the increasing number of degrees of symmetry for the flock solution and the weakening of the interaction potential observed in [35]. The investigation of stability of other coherent patterns of swarming models such as rotating mills is another highly interesting open problem where the techniques used in this work clearly cannot be applied straightforwardly.