Stability of Translating States for Self-propelled Swarms with Quadratic Potential

Second order models of self-propelled agents with all-to-all gradient coupling have been extensively studied, notably by researchers in Prof. Bertozzi’s group ([6], for continuous models, and [7] for agent-based models). Let $K:[0,\infty)\to\mathbb{R}$ be a continuous scalar function, differentiable on $(0,\infty).$ Define $U:\mathbb{R}^{d}\to\mathbb{R}$ as $U(x)=K(|x|).$ With the possible exception of the origin, the potential $U$ is differentiable with $\nabla U(x)=K^{\prime}(|x|)\frac{x}{|x|}.$ Consider the system with identical self propelling terms and rotationally symmetric potential:

where $\alpha,\beta$ are positive constants.

The system (1) is in the class of (5) with $K(r)=\frac{1}{2}r^{2}$ and $\alpha=\beta=1.$ Other well known classes of radial potentials are Morse potentials where $K$ is the difference between two exponential decays (see [8], [9]), and power law potentials where $K(r)=c_{a}r^{a}-c_{b}r^{b},$ ( [18], Section 8; [19]; [4], 3.3-3.5). Given the symmetry of (5) it is natural to expect that all emergent patterns feature rotational symmetries; in that spirit, with the exception of [4] and [5], most theoretical results address particles that are uniformly distributed around their center of mass and their similarly restricted perturbations, or have assumptions about the degeneracies of the system that for (1) only hold for $n\leq 2$, [13]. We emphasize that we make no uniform distribution assumption in this paper, which explains the sharp distinction between the rates of convergence to limit configuration: exponential rates obtained by authors working exclusively in the subset with rotational symmetries, versus $\frac{1}{\sqrt{t}}$ in this paper (and in [5]). In fact (1) has very few configurations with symmetries:

It was shown in [5] that the only possible configurations of (1) with a stationary center of mass come from partitioning the agents $\{1,\dots n\}$ into subsets $\mathcal{F},\mathcal{L},\mathcal{R}$ where the agents in $\mathcal{F}$ have fixed positions coinciding with the center $R(0)$ at all times; each agent in $\mathcal{L}$ oscillates on a segment centered at $R(0)$, satisfying the scalar Lienard equation $\ddot{a}=(1-\dot{a}^{2})\dot{a}-a$ (which has an attracting cycle); each agent in $\mathcal{R}$ is asymptotic to a rotation $e^{\pm it}e^{i\theta_{\infty,k}}.$

Some authors, but not all, use the term mill for a rotating state; most authors define a mill as a rotating state whose angles $\theta_{0,k}$ are uniformly distributed in $[0,2\pi]$, a.k.a having rotational symmetry. All the rotating states of (1) are stable, per [5]. Configurations with rotational symmetries for this system and for other radial potentials have been extensively studied: see [16], [8], [9], [18] for some early theoretical and experimental results.

Note that for a flocking state $\ddot{r}_{k}=\dot{V}.$ If we subtract the equation for $\ddot{r}_{k_{1}}$ from that of $\ddot{r}_{k_{2}}$ we get that $r_{k_{1}}=r_{k_{2}}$ for all $k_{1},k_{2}$ meaning that all particles have collapsed onto the center of mass, moving together according to

Unless $V(t)=0$ for all $t$ (the trivial solution), the velocity of the flocking solutions maintains the initial direction of $V(0)=|V(0)|e^{i\Theta_{0}}$ with speed approaching unit speed at a rate of $e^{-2t}$ (since $|V|^{2}$ solves the logistic equation $\displaystyle\frac{d}{dt}(|V|^{2})=2(1-|V|^{2})|V|^{2}\;).$ A nontrivial flocking state has velocity approaching $e^{i\Theta_{0}}$ exponentially fast. Moreover, if $c_{\infty}$ denotes the absolutely convergent integral $\displaystyle c_{\infty}=\int_{0}^{\infty}\left[V(\tau)-e^{i\Theta_{0}}\right]d\tau$ then the flocking state’s center of mass has nearly rectilinear motion: $\lim_{t\to\infty}\left[R(t)-(c_{\infty}+te^{i\Theta_{0}})\right]=0.$ The limit holds since $R(t)-te^{i\Theta_{0}}-c_{\infty}=\int_{0}^{t}\left[V(\tau)-e^{i\Theta_{0}}\right]d\tau-\int_{0}^{\infty}\left[V(\tau)-e^{i\Theta_{0}}\right]d\tau=-\int_{t}^{\infty}\left[V(\tau)-e^{i\Theta_{0}}\right]d\tau.$

The main result of this paper is proving that the translating states are stable, in the sense that if the particles in (1) start close to each other, with velocities near some unit vector $(\cos\Theta_{0},\sin\Theta_{0})$ then the particles remain close to each other, with velocities that are asymptotic to a common unit vector, whose direction does not drift far from $\Theta_{0},$ as stated in Theorem 1.6. We also quantify the rate of convergence to the limiting configuration (Corollary 4.4.1), and we show that it is at a rate of $\displaystyle\frac{1}{\sqrt{t}},$ much slower than the conjectured exponential rate.

The proof of Theorem 1.6 is completed in three parts, corresponding to three dimension reductions of the system: form $4n$ to $4n-3$ in Section 2, to $2n$ in Section 3, and to $n$ in Section 4, with the stability of the translating states for (1) being equivalent to that of the origin for the reduced systems.

The dynamical system produced in Section 2 evolves in a $4n-3$ dimensional invariant subspace for a system taking the form $\dot{w}=Aw+G(w)$ where $A$ is the Jacobian matrix and $G$ is a nonlinear function whose Maclaurin expansion has terms or order 3 or higher. Counted with multiplicity, the eigenvalues of $A$ consist of $n$ pairs of purely imaginary roots $\pm i$, and $2n-1$ eigenvalues with negative real parts. The stability of the origin for systems whose Jacobian matrices have $q$ pairs of purely imaginary roots $(\pm i\omega_{1},\dots,\pm i\omega_{q})$ and $d-2q$ roots with negative real part has been extensively investigated; the problem is commonly known as the critical case with purely imaginary roots. Most authors apply elaborate iterative transformations from normal form theory to obtain systems for which stability criteria such as Molchanov and Lyapunov’s second method can be employed. These transformations require that the frequencies $\omega_{k}$ satisfy non-resonance conditions which are violated in our setting,([33], [34]), or alternatively, they assume the a priory existence of some quadratic Lyapunov functions ( [35]), rendering their methods inapplicable here despite concluding similar decay rate of $1/\sqrt{t}$ towards the origin. Our approach draws inspiration from a class of time-dependent systems with cubic nonlinearities instead:

where $x=(x_{1},\dots x_{n})^{T}$ is vector-valued from $\mathbb{R}$ to $\mathbb{R}^{n}$, the coefficient functions $b_{k},c_{k},d_{k}:\mathbb{R}\to\mathbb{R}$ are continuous and the remainder functions $R_{k}$ defined on $\mathbb{R}^{n}\times\mathbb{R}$ satisfy the $\textsl{o}_{3}$ condition: there exists a function $\epsilon:[0,\infty)\to\mathbb{R}$ with $\displaystyle\lim_{r->0}\frac{\epsilon(r)}{r^{3}}=0$ such that $|R_{k}(x,t)|\leq\epsilon(|x|)$ for all $(x,t)\in\mathbb{R}^{n}\times\mathbb{R}.$ In Section 3.3 we provide easily verifiable conditions on the coefficient functions that sufficiently ensure the asymptotic stability of the solution $x=\mathbb{0}_{n}$, with convergence rate of order $\frac{1}{\sqrt{t}}.$ We adapt the arguments from Section 3.3 to complete the proof of Theorem 1.6, in Section 4.1. We use the following to facilitate notations:

We conclude this introductory subsection by examining how other authors interpret stability for limit patterns and comparing their settings to our results, as stated in Theorem 1.6, Remarks 4.5, and 4.7, since there is no universally accepted definition.

For Cucker-Smale (CS) models with communication function $\psi$

the generally accepted definition of asymptotic flocking requires two conditions:

• •

  Velocity alignment: $\lim_{t\to\infty}\sum\limits_{j=1}^{n}|\dot{r}_{j}(t)-V(t)|^{2}=0$

• •

  Forming a finite-diameter group: $\sup_{0\leq t<\infty}\sum\limits_{j=1}^{n}|r_{j}(t)-R(t)|^{2}<\infty.$

Cucker-Smale models have constant mean velocity $V(t)=V(0)$ so the existence of a limiting value of $V$ is automatic; the problem reduces to an equivalent model with $V(t)=0.$ For the classical communications functions $\psi(s)=\frac{\alpha}{s^{\beta}}$ and $\psi(s)=\frac{\alpha}{(1+s^{2})^{\beta/2}}$ with $0\leq\beta\leq 1$ it is proven in [10] and [3] that there is unconditional flocking in the sense that all solutions exhibit asymptotic flocking with $|\dot{r}_{j}(t)-V(t)|$ being exponentially small.

For Kuramoto oscillators (KO)

the mean frequency $\omega=\frac{1}{n}\sum_{j=1}^{n}\omega_{j}$ is constant, ensuring that the mean phase $\frac{1}{n}\sum_{j=1}^{n}\theta_{j}$ is linear in time. For (KO) synchronization plays the role of (CS)’s velocity alignment, requiring that the angular frequencies $\dot{\theta}_{k}$ converge to the mean frequency. It was shown that under some mild conditions on the initial distribution of phase angles, even non-identical oscillators synchronize [31], and that the exponential rate of convergence proven in [30] for identical oscillators (i.e. when all $\omega_{k}=\omega$ ) extends to more general systems. Not all systems synchronize, and some do at much slower rates. Moreover, [30] shows that identical oscillators with initial phase angles confined in an semicircle have phase angles that converge exponentially to a common $\theta_{\infty}.$

For a swarm of identical planar agents coupled with a power law potential as in (5), some authors ([32]) define flocking configurations as being solutions where all agents have rectilinear motion, with a common velocity: $r_{k}(t)=r_{k,0}+v_{0}t$ for all $k.$ Necessarily $v_{0}=\sqrt{\frac{\alpha}{\beta}}e^{i\theta_{0}}$ for some $\theta_{0}\in[0,2\pi).$ Let $r_{\star}=(r_{k,0})_{k=1..n}$ in $\mathbb{R}^{2n}$ denote the initial locations of a flocking configuration. In [32] the authors define the stability of the flocking state evolving from $r_{\star}$ by referring to a manifold of configurations. They define the flock manifold $\mathcal{F}(r_{\star})\subset\mathbb{R}^{4n}$ as

and the ’local asymptotic stability’ of $\mathcal{F}(r_{\star})$ as the existence of a $\delta-$ neighborhood $U$ of $\mathcal{F}(r_{\star})$ such that for initial conditions $r_{0},\dot{r}_{0}$ in $U,$ the distance from the state variables $r(t),\dot{r}(t)$ to the manifold $\mathcal{F}(r_{\star})$ approaches zero as $t\to\infty.$ Note that this is a rather weak condition on the behavior of $V(t)$, as it only requires that $|V|\to\sqrt{\frac{\alpha}{\beta}}$ without imposing any conditions on $V$ having a steady direction.

One of the difficulties in working with translating and flocking states is the fact that they are neither fixed points nor periodic cycles. Moreover, for (1) the perturbations of the translating states have centers of mass that do not move with constant velocity (or in a predetermined direction) as in Cucker-Smale where excising a linear-in-$t$ term out of the motion transforms the flocking states into fixed points.

In [4] the authors prove that there exit constants $C_{r},C_{v},C_{a}$ that depend only on the number of particles, such that every solution to (1 ) satisfies $|r_{k}(t)-R(t)|\leq C_{r},\;|\dot{r}_{k}(t)-V(t)|\leq C_{v},$ and $|\ddot{r}_{k}(t)-\dot{V}(t)|\leq C_{a}$ for large enough $t.$ This suggests two possible workarounds for the unboundedness of flocking states: study (4) as a substitute for (1) or change the unknown functions of (1 ) from $r_{k}$ and $\dot{r}_{k}$ to $r_{k}(t)-R(t)$ and $v_{k}=\dot{r}_{k}.$ Both approaches lead to dynamics evolving in a compact phase space, with the translating states of (1 ) corresponding to fixed points with all velocities $v_{k}$ equal to the mean speed $V_{0},$ a unit vector: $|V_{0}|=1.$

This subsection briefly addresses some remaining obstacles in pursuing these model modifications. Although this subsection helps set the tone for the dimension reduction and the proof of Theorem 1.6 that begins in Section 2, the proof of the main results are independent of it.

Note that the fixed points of (4) with $|V_{0}|=1$ correspond to the translating solutions of (1). We want to emphasise that the perturbations that destabilize (4) don’t correspond to derivatives of solutions to (1). The main theorem (1.6) in this paper makes this point.

The second natural approach for transforming the unbounded solutions of (1) and its translating states into bounded trajectories and fixed points is to work with $\sigma_{k}=r_{k}-R=(\sigma_{x,k},\sigma_{y,k})$ and $\dot{r}_{k}=(v_{x,k},v_{y,k})$ as the unknown functions. Note that since $r_{n}-R=-\sum_{k=1}^{n-1}\sigma_{k}$ we can eliminate it from our unknowns, and work within a $4n-2$ dimensional phase space (if we were to keep the redundant $r_{n}-R$ as a state variable, we would have added two more zero eigenvalues to the system).

The system (1) becomes:

where $k=1..(n-1).$

The translating states of (1) with $V(t)=(\cos\theta_{0},\sin\theta_{0})$ correspond to the fixed points $P_{\theta_{0}}$ with $\sigma_{x}=\sigma_{y}=\mathbb{0}_{n-1},$ and $v_{x}=\cos\theta_{0}\;\mathbb{1}_{n},\;v_{y}=\sin\theta_{0}\;\mathbb{1}_{n}.$

Focus on the dynamics near $P_{0}=\mbox{col}(\mathbb{0}_{n-1},\mathbb{0}_{n-1},\mathbb{1}_{n-1},\mathbb{0}_{n-1},1,0).$ The other fixed points are obtained from $P_{0}$ via rotation of angle $\theta_{0}.$ The linearization of (13) at $P_{0}$ yields the Jacobian matrix

where $B_{n-1}$ is the square matrix $\textrm{I}_{n-1}-\frac{1}{n}\mathbb{1}_{n-1,n-1}$ and $\mathbb{O}=\mathbb{O}_{n-1,n-1}.$ We get that

The eigenvector for $\lambda=0$ is $\mbox{col}(\mathbb{0}_{n-1},\mathbb{0}_{n-1},\mathbb{0}_{n-1},\mathbb{1}_{n-1},0,1)$; it has the direction of the line tangent to the curve of fixed points $\{P_{\theta_{0}},\;\theta_{0}\in[0,2\pi]\}$ at $P_{0}.$ The eigenvalue $\lambda=0$ captures the rotation invariance of the model. Note that due to the $\pm i$ eigenvalues there are more degeneracies associated with (13) than the symmetries of the problem.

If in the system (13) all particles’ $y-$components are initially zero, they stay zero forever, with the remaining $2n-1$ dimensional system in $\sigma_{x},v_{x},v_{x,n-1}$ having the Jacobian $J_{x}$ at $\sigma_{x}=\mathbb{0}_{n-1},v_{x}=\mathbb{1}_{n-1},v_{x,n-1}=1$ equal to

and characteristic polynomial $(\lambda+2)(\lambda+1)^{2n-2}.$ In the absence of a $y-$ component, all trajectories approach the fixed point $P_{0}$ exponentially fast, meaning that the entire stable manifold for (13) at $P_{0}$ consists of the trajectories with no $y$ component. The flow on the $2n-1$ dimensional central manifold captures the drift in the direction of motion, and the perturbations/ oscillations that are orthogonal to the direction of motion.

In order to understand the dynamics, one needs a suitable approximation of a central manifold map $h$ whose graph is a central manifold, typically done using Taylor polynomials (as in [11], pages 3-5). Unfortunately, Taylor approximations are ill fitted to capture the correct scales of the flow in the presence of non-isolated fixed points, which is the case for (13). A presentation of that inadequacy is in [5] ¹¹1[5] provides an alternate construction for the central manifold for systems with non-isolated fixed points having no purely imaginary eigenvalues. Their method does not apply in the context of (13). The predicament of having non-isolated fixed points and imaginary eigenvalues of high multiplicity requires we further refine the model until all the translating states are identified with a single fixed point. That is done in Section 2.

For example, starting from the flocking solution $r_{k}(t)=(0,0,t)$ for all $t$ and $k$, given $\epsilon>0$ small, consider the perturbed initial conditions $r_{k,0}=(\epsilon\cos\theta_{k},\epsilon\sin\theta_{k},0)$ and $\dot{r}_{k,0}=(-\epsilon\sin\theta_{k},\epsilon\cos\theta_{k},\sqrt{1-\epsilon^{2}})$ where $\theta_{k}$ are angles such that $\displaystyle\sum_{k=1}^{n}e^{i\theta k}=0.$ Note that for these perturbations $R(0)=(0,0,0)$ and $V(0)=(0,0,\sqrt{1-\epsilon^{2}}).$ It is immediate to verify that the particles in the $\mathbb{R}^{3}$ version of (1) move with unit speed on helix-shaped trajectories according to

The perturbed trajectories have $R(t)=(0,0,t\sqrt{1-\epsilon^{2}})$ and therefore $|r_{k}(t)-R(t)|=\epsilon$ and $|\dot{r}_{k}(t)-V(t)|=\epsilon$ for all $t,$ failing the asymptotics (7) and (8).