2 The Model

The model we study in this report has been introduced in [7]. Let us introduce the model and highlight its main features. Consider a system of \(N\) interacting particles evolving on the one-dimensional torus. For each \(i\in \{1, \dots, N\}\), denote by \(x_t^{i,N} \in \mathbb{T}\) and \(v_t^{i,N} \in \mathbb{R}\), respectively, the position and velocity of the \(i\)-th particle at time \(t\); the dynamics of \(x_t^{i,N}\) and \(v_t^{i,N}\) is described by the following system of SDEs. \[\begin{align} \mathrm{d} {x}_t^{i,N} & = v_t^{i,N} \, \mathrm{d} t,\tag{2.1} \\ \mathrm{d} {v}_t^{i,N} & = - v_t^{i,N} \mathrm{d} t + G\left(\frac{ \sum_{j=1}^N \varphi (x_t^{i,N} - x_t^{j,N}) v_t^{j,N}}{ \sum_{j=1}^N \varphi (x_t^{i,N} - x_t^{j,N})}\right) \mathrm{d} t + \sqrt{2\sigma} \mathrm{d} W_t^{i},\tag{2.2} \end{align}\] where the \(W^i\)’s are independent one-dimensional standard Wiener processes and \(\sigma>0\) is fixed. The first equation simply states that the change in position is given by the velocity multiplied by the change in time. The equation governing the velocity however, contains more complexities.

The interaction function \(\varphi: \mathbb{T} \to \mathbb{R}\) controls which other particles the i\(^{\text{th}}\) particle interacts with and how strongly. In [7], to avoid technical problems with the denominator of the argument of \(G\), the function\(\varphi\) is assumed to be bounded below. We further assume \(\varphi\) is symmetric so that whether a particle is “in front” or “behind” of another has no effect on its interaction strength. More formally, \[\begin{equation} \varphi(x) \geq \varepsilon > 0,\quad \varphi(x) = -\varphi(x). \tag{2.3} \end{equation}\]

When numerically simulating the system, we are able to drop the positivity assumption. If \(\varphi \equiv 1\), this corresponds to a mean-field interaction and all particles continually interact with all others equally, irrespective of their relative positions. This type of interaction has been extensively studied and indeed forms the basis for much of statistical physics. This type of “all to all” interaction is not necessarily meaningful for many systems in biology or in the social sciences—for example, people with very different opinions are unliekly to interact. More realistic is that the agent’s ability to interact varies with distance and direction, for example when birds move out of each other’s line of sight. A simple example could be an indicator function on some set \([-\gamma, \gamma]\). Any particles that lie more than a distance \(\gamma\) apart will no longer interact. For this reason, \(\gamma\) is often called the radius of interaction. This type of interaction function is common in opinion dynamics and leads to the formation of the 2R conjecture [8]. It could be argued that this form of interaction with an instant decline is implausible for many multi-agent systems and more relevant is a slow decay to no interaction. This prompts the introduction of smoothed interaction functions such as bump functions centred at zero. Figure 2.1 illustrates these standard choices.

Figure 2.1: Some typical choices of interaction functions for this and other particle models. The solid blue line a) shows a uniform interaction \(\varphi(x) \equiv 1\), dashed orange b) shows a hard cutoff indicator \(\varphi(x) = \mathbb{I}_{[-\pi\gamma,\pi\gamma]}(x), \gamma = 0.2\) and the dotted green line c) shows a smooth bump function interaction.

We thus see that the argument of the function \(G\) is a weighted local average of the particles’ velocities. Indeed if \(\varphi\equiv 1\), this is simply the mean velocity of the ensemble. The perhaps more interesting case is when the interaction is not uniform and there is some spatial heterogeneity.

Now we better understand the role of the argument of \(G\), what of \(G\) itself? The herding function, \(G:\mathbb{R}\to \mathbb{R}\), defines how an individual should respond to the weighted average velocity. That is, how does the individual react to the velocities of the particles surrounding it? First notice that when \(G\equiv0\), there is no interaction between particles and the system reduces to \(N\) uncoupled Ornstein-Uhlenbeck processes. If there is to be an interaction, it is natural that this function should at some point satisfy \(G(u)=u\), so that particles are less perturbed when they are moving with the crowd, and also that \(G\) is constructed so that an individual’s velocity moves towards that of the herd. This can be achieved under the following assumption: \[ \tag{2.4} G(u) = - G(-u)\,, \qquad \begin{cases} G(u) > u & \textrm{ if } 0\leq u<1\,, \\ G(u) < u & \textrm{ if } u>1\,. \end{cases} \] The choice of \(G(1)=1\) here is chosen arbitrarily, the model displays similar dynamics for any \(\xi\) such that \(G(\xi)=\xi\). We also restrict our attention to the case when there are only two solutions to this equation. This function thus herds the velocities of the particles towards a fixed point. Note also that \(G\) controls how quickly an individual moves towards the herd average, or the speed of reaction. Within the modelling literature, it is common to take a step function for \(G\): \[ G(u) = \frac{u+\beta \mathrm{sgn}(u)}{1+\beta}\, , \beta>0 \]
The parameter \(\beta\) controls the strength of the herding function, a typical value is \(\beta=2\). The advantage of this choice is the constant gradient for any \(u\neq0\), but the discontinuity at this point renders them less amenable to analysis—and potentially to numerics, as we will illustrate. For this reason, analysts prefer a smoothed version such as \[ G_{\alpha}(u) = \frac{\mathrm{atan}(\alpha u)}{\mathrm{atan}(\alpha)}\, , \alpha > 0 . \] A herding function such as this still satisfies the assumptions, whilst also being smooth. The parameter \(\alpha\) again controls the gradient, a typical value is \(\alpha = 1\). Figure 2.2 provides some illustrations and in Section 3 we will consider both of these types of herding functions.

Figure 2.2: Some typical choices of herding functions for this model. The solid blue line a) shows the step herding function, \(G(u) = \frac{u+\beta \mathrm{sgn}(u)}{1+\beta}, \beta=2\). The orange dashed line b) shows the smoothed herding function $G(u) = , =1 $ and the dotted green line shows \(G(u)=\frac{h+1}{5}u-\frac{h}{125}u^3, h=5\), the herding function used in the numerical simulations of [9]. The dotted grey line is \(G(u)=u\). Note the intersection between the grey line and others describes the value to which the average velocity will be herded. For a) and b), this is \(0, -1, +1\). For c), the values are \(0, \frac{\pm5}{\sqrt{3}}\).

Finally, we are left with two terms to discuss. The first term in (2.2) is a friction term, meaning that if there is no noise any individual particle will come to rest in the absence of neighbours and noise. If there are two particles in the system, there are two possibilities: either they do not interact sufficiently and both come to rest; or they coalesce and move together with common velocity. The latter term describes the noise in the system with \(\sigma\) representing the strength. A common characteristic of these types of models is a phase transition between order and disorder dependent on \(\sigma\). Interestingly, it was proved in [7] that this model displays unconditional flocking—there is no transition between ordered motion and disorder. Note that this noise is only present in the velocity equation (2.2), the position equation (2.1) is entirely deterministic. This means the system is not elliptic, and renders it immune to many analysis techniques.

When considering a particle model, one often looks to the empirical measure. This can be thought of as a histogram of the positions and velocities at each time \(t\). It is defined as \[ S^N_t(\mathrm{d} x, \mathrm{d} v) := \frac{1}{N}\sum_{i=1}^N \delta_{(x^{i,N}_t,v^{i,N}_t)}(\mathrm{d}x, \mathrm{d}v). \] Here \(\delta\) denotes a Dirac measure, that is on any measurable subset \(A\) of some set \(X\), \[ \delta_x(A) = \begin{cases} 0 \text{ if } x \not\in A,\\ 1 \text{ if } x \in A. \end{cases} \] In [7] it was shown that, under certain assumptions, \(\lim_{N\to \infty} S^N_t = f_t\) where \(f_t\) solves \[\begin{equation} \tag{2.5} \partial_t f_t(x,v)=- v\,\partial_x f_t(x,v) - \partial_v\big\{\big[G(M_{f_t}(x)) - v\big] f_t(x,v) \big\} + \sigma\,\partial_{vv}f_t(x,v)\,, \end{equation}\] where \[\begin{equation} \tag{2.6} M_f(x) := \frac{\int_{\mathbb{T}}\!\mathrm{d} y \int_{\mathbb{R}}\!\mathrm{d} w\, f(y,w)\,\varphi(x-y)\, w}{\int_{\mathbb{T}}\!\mathrm{d} y \int_{\mathbb{R}}\!\mathrm{d} w\, f(y,w)\,\varphi(x-y)}\,. \end{equation}\] This is a nonlinear kinetic equation, the analysis of which will help inform our approach to the particle system. We will now summarise what is known about this PDE. The well-posedness of (2.5) was proved in [7] and that a unique solution exists in a weighted \(L^1\) space \[ L^1\left(\sqrt{1+v^2}\mathrm{d}x \mathrm{d}v\right) := \bigg\lbrace f:\mathbb{T}\times\mathbb{R} \to \mathbb{R} : \int_{\mathbb{T}}\int_{\mathbb{R}} |f| \sqrt{1+v^2} \mathrm{d}v\mathrm{d}x< \infty \bigg\rbrace. \] A complete analytic description of the long-time behaviour of (2.5) seems to be out of reach (we explain below why this is the case); however, the asymptotic behaviour of this dynamics has been characterized in some simplified cases, namely 1) when the density \(f_t\) is space-homogeneous, i.e. when it does not depend on the space-variable \(x\); 2) when the function \(\varphi\) is a perturbation of the mean field case and 3) when the non- linearity (2.6) is simplified. Let us comment on these three cases in turn.

1. If we consider only space- homogeneous densities, i.e. solutions of the form \(f_t(v)\), then the evolution (2.5)-(2.6) reduces to the (much) simpler dynamics \[\begin{equation} \tag{2.7} \partial_t f_t(v)=- \partial_v\big\{\big[G(M(t)) - v\big] f_t(v) \big\} + \sigma\,\partial_{vv}f_t(v)\,, \end{equation}\] where \(M(t):= \int_{\mathbb{R}}v\, f_t(v) \mathrm{d}v\) solves the one-dimensional ODE \[\begin{equation} \tag{2.8} \dot M (t)= G(M(t))-M(t) . \end{equation}\] In other words, the nonlinearity \(M_{f_t}(x)\) reduces to being the average velocity \(M(t)\); because one can close an equation on \(M(t)\), one can now regard \(M(t)\) as being a time-dependent coefficient; hence, in this space-homogeneous case, the nonlinear dynamics (2.5) reduces to a linear, non-autonomous PDE. In this case one can easily prove that there exist three stationary solutions of (2.7), namely three Gaussians \[\begin{equation} \tag{2.9} \mu_{\pm}= \frac{1}{Z_{\pm}} e^{-\frac{(v\mp 1)^2}{2 \sigma}} \quad \mbox{ and } \quad \mu_{0}= \frac{1}{Z_{0}} e^{-\frac{v^2}{2 \sigma}} \, , \end{equation}\] where \(Z_{\pm}\) and \(Z_0\) are appropriate normalization constants. Note that these are the only three stationary solutions of the dynamics (2.7)-(2.8), irrespective of the value of \(\sigma\). This is the unconditional flocking mentioned earlier. Moreover, if the initial profile \(f_0\) has positive mean, i.e. if \(M(0)>0\) (\(M(0)<0\), respectively) then the dynamics converges exponentially fast to the stationary measure with positive mean, \(\mu_+\) (\(\mu_-\), respectively). See Appendix A for an explanation.

2. As we have mentioned, the study of the long-time behaviour of the full model (2.5)-(2.6) is quite challenging. In particular one can easily see that the three measures (2.9) are still stationary solutions of (2.5)-(2.6); however the authors of [7, 9] were not able to prove that these are indeed the only three stationary measures; this is a conjecture that we aim to confirm numerically. The fact that the measures (2.9) are the only stationary solutions of (2.5)-(2.6) has been proven analytically only in the case in which the function \(\varphi\) is chosen to be a small oscillation around the constant function \(\varphi \equiv 1\); to be more precise, this has been proven (under some assumptions on the regularity of the solution), only in the case in which \[ \varphi(x)= 1+ \lambda \psi(x), \quad \mbox{where } 0<\lambda \ll 1, \int_{\mathbb{T}} \psi(y) \, dy =0 \,. \] The difficulty in solving the stationary problem associated with (2.5)-(2.6) comes from the fact that such a model is not in gradient form.

3. In [9], the authors introduced a model similar to (2.5)-(2.6), namely they study the following non-linear PDE \[\begin{equation} \tag{2.10} \partial_t f_t(x,v)=- v\,\partial_x f_t(x,v) - \partial_v\big\{\big[G(\tilde{M}_{f_t}(x)) - v\big] f_t(x,v) \big\} + \sigma\,\partial_{vv}f_t(x,v)\,, \end{equation}\] where this time nonlinearity \(\tilde{M}_f\) is given by \[\begin{equation} \tag{2.11} \tilde{M}_f(x) := {\int_{\mathbb{T}}\!\mathrm{d} y \int_{\mathbb{R}}\!\mathrm{d} w\, f(y,w)\,\varphi(x-y)\, w}\,. \end{equation}\] This corresponds to the following particle system \[\begin{align} \mathrm{d} {x}_t^{i,N} & = v_t^{i,N} \, \mathrm{d} t, \tag{2.12} \\ \mathrm{d} {v}_t^{i,N} & = - v_t^{i,N} \mathrm{d} t + G\left(\frac{ \sum_{j=1}^N \varphi (x_t^{i,N} - x_t^{j,N}) v_t^{j,N}}{N}\right) \mathrm{d} t + \sqrt{2\sigma} \mathrm{d} W_t^{i}.\tag{2.13} \end{align}\] Comparing this with (2.1)-(2.2), notice that the only difference is in the denominator of the argument of the herding function \(G\). This choice of normalisation is akin to scaling the interaction by the total number of particles in the system. For this reason, we call this global scaling. In contrast, we call the interaction term in (2.1)-(2.2) locally scaled as there the interaction is normalised by the amount of particles that are interacting with an individual. Scaling globally may seem counterintuitive: the dynamics of a flock of birds in one continent are surely not affected by a flock in another continent. We will show that in some regimes (for example with a large particle count \(N\) and large support of \(\varphi\)) these models are clearly very similar. We only expect a difference in the dynamics for smaller particle numbers, or where there is a large distance between agents. Global scaling is a common modelling assumption in the analysis literature as it reduces the convergence of the measure to a well-studied case—that of the mean field interaction. In Section 3 we will investigate numerically the difference in these two models. For the globally scaled model (2.10), it has been shown that the only stationary states are (2.9) [9].

Less is known about the particle model (2.1)-(2.2). When noise is present in the system, the stationary states of (2.5) correspond to metastable states in the particle system. It can be shown that when \(\varphi\equiv1\) this system has an invariant measure with mean zero, see Appendix B. We conjecture that the same is true when \(\varphi \not\equiv 1\), and aim to prove this in a future work. If there is no noise in the system and a uniform interaction is assumed, then the velocity of all particles converges to the point \(\xi\), where \(\xi\) is one of the points such that \(G(\xi)=\xi\).

In summary, the two models presented in [7] & [9] are very similar, however this discrepancy in the argument of the herding function will in fact have large effects on the dynamics. Furthermore, the analytic results for the locally scaled particle system (2.1)-(2.2) and its corresponding kinetic PDE (2.5) do not give a complete picture of the dynamics. We thus turn to a numerical approach to inform the analysis in the case when the interaction is not uniform. Our first question is whether a change in \(\varphi\) affects the rate of convergence to stationarity. Second, are there any spatially inhomogeneous solutions to the kinetic equation? Do particle distributions always tend towards being uniform in space? In addition, when \(\varphi\not \equiv1\), the locally scaled (2.1)-(2.2) and globally scaled (2.12)-(2.13) models are no longer equivalent. What differences in behaviour do they exhibit?

References

[7] Buttà, P., Flandoli, F., Ottobre, M. and Zegarlinski, B. (2019). A non-linear kinetic model of self-propelled particles with multiple equilibria. Kinetic & Related Models 12(4) 791–827.

[8] Blondel, V. D., Hendrickx, J. M. and Tsitsiklis, J. N. (2007). On the 2R conjecture for multi-agent systems. In 2007 european control conference (ecc) pp 874–81.

[9] Garnier, J., Papanicolaou, G. and Yang, T.-W. (2019). Mean field model for collective motion bistability. Discrete & Continuous Dynamical Systems - B 24 851.