Kinetic equations and self-organized band formations

TL;DR

This paper investigates kinetic equations related to self-organized band formations in collective motion, presenting a numerical scheme that preserves physical properties and analyzing the emergence of stable band patterns.

Contribution

It introduces a numerical scheme for the kinetic Vicsek model that maintains key physical properties and provides insights into the formation of stable band patterns.

Findings

Numerical scheme preserves energy and positivity.

Pattern of bands emerges stably in simulations.

Insights into band formation mechanisms.

Abstract

Self-organization is an ubiquitous phenomenon in nature which can be observed in a variety of different contexts and scales, with examples ranging from fish schools, swarms of birds or locusts, to flocks of bacteria. The observation of such global patterns can often be reproduced in models based on simple interactions between neighboring particles. In this paper we focus on two particular interaction dynamics closely related to the one described in the seminal paper of Vicsek and collaborators. After reviewing the current state of the art in the subject, we study a numerical scheme for the kinetic equation associated with the Vicsek models which has the specificity of reproducing many physical properties of the continuous models, like the preservation of energy and positivity and the diminution of an entropy functional. We describe a stable pattern of bands emerging in the dynamics proposed by Degond-Frouvelle-Liu dynamics and give some insights about their formation.

Figures

Tables

Table 1. Table 1 : Parameters used in the simulations for Figures 1 - 2

Description	notation	value
Number particles	$N$	$30,000$
Strength alignment	$\mu$	$100$
Noise intensity	$\sigma$	$20$
Radius interaction	$R$	$.02$
Length domain	$L$	$4$
Time step	$\Delta$t	${10}^{-2}$

Equations

$$\begin{array}[]{rcl}{\bf x}_{i}^{\prime}&=&c\,\omega_{i}\\ \mathrm{d}\omega_{i}&=&P_{\omega_{i}^{\perp}}(\mu\,\Omega_{i}\mathrm{d}t+\sqrt{2\sigma}\circ\mathrm{d}B_{i}^{t}),\end{array}$$

$$P_{{\omega }_i^{\perp }}=\text{mbox}Id-{\omega }_i\otimes {\omega }_i,$$

$${\Omega }_i=\frac{{\mathbf{j}}_i}{|{\mathbf{j}}_i|},{\mathbf{j}}_i=\sum _{j,|{\mathbf{x}}_j-{\mathbf{x}}_i|\le R}{\omega }_j,$$

$$\begin{array}[]{rcl}{\bf x}_{i}^{\prime}&=&c\omega_{i}\\ \mathrm{d}\omega_{i}&=&P_{\omega_{i}^{\perp}}(\mu\,{\bf j}_{i}\mathrm{d}t+\sqrt{2\sigma}\circ\mathrm{d}B_{i}^{t}).\end{array}$$

$$\overline{\rho }(x,t)\Delta x$$

$$\overline{\rho }(x,t)\overline{u}(x,t)$$

$$\text{Average number of neighbors}\approx \frac{|B(0,R)|}{L^2}\times N=2.36,$$

$${\partial }_{t}f+c\omega \cdot {\nabla }_{\mathbf{x}}f=-\mu {\nabla }_{\omega }\cdot (P_{{\omega }^{\perp }}({\Omega }_f)f)+\sigma {\Delta }_{\omega }f,$$

$${\Omega }_f(\mathbf{x},t)=\frac{{\mathbf{j}}_f(\mathbf{x},t)}{|{\mathbf{j}}_f(\mathbf{x},t)|}\text{with}{\mathbf{j}}_f(\mathbf{x},t)={\int }_{y\in B(\mathbf{x},R),\omega \in {\mathbb{S}}^{d-1}}\omega f(\mathbf{y},\omega ,t)\mathrm{d}\mathbf{y}\mathrm{d}\omega ,$$

$${\partial }_{t}f+c\omega \cdot {\nabla }_{\mathbf{x}}f=-\mu {\nabla }_{\omega }\cdot (P_{{\omega }^{\perp }}({\mathbf{j}}_f)f)+\sigma {\Delta }_{\omega }f.$$

$${\partial }_{t}f=Q(f),$$

$$Q(f)$$

$$\mu_{f}=\left\{\begin{array}[]{ll}\mu&\text{Vicsek dynamics}\\ &\\ \mu|{\bf j}|&\text{DFL dynamics}.\end{array}\right.$$

$$\phi(\omega)=\left\{\begin{array}[]{ll}\langle{\bf j}_{f},\omega\rangle&\text{Vicsek dynamics}\\ &\\ \langle\Omega_{f},\omega\rangle&\text{DFL dynamics},\end{array}\right.$$

$$Q(f)=\sigma {\nabla }_{\omega }\cdot \left(M_f{\nabla }_{\omega }\cdot \left(\frac{f}{M_f}\right)\right),\text{with}{M}_f(\omega )={\mathrm{e}}^{\frac{\mu }{\sigma }\varphi (\omega )},$$

$${\int }_{\omega }{\partial }_{t}f\frac{f}{M_f}\mathrm{d}\omega =-\sigma {\int }_{\omega }{M}_f{\left|{\nabla }_{\omega }\left(\frac{f}{M_f}\right)\right|}^2\mathrm{d}\omega \le 0.$$

$$Q(f)$$

$${\int }_{\omega }{\partial }_{t}f\ln \left(\frac{f}{M_f}\right)\mathrm{d}\omega =-\sigma {\int }_{\omega }f{\left|{\nabla }_{\omega }\ln \left(\frac{f}{M_f}\right)\right|}^2\mathrm{d}\omega \le 0.$$

$$\mathcal{F}[f]={\int }_{\omega }f\ln f\mathrm{d}\omega -\frac{\mu }{\sigma }{\Phi }_f,$$

$$\Phi_{f}=\left\{\begin{array}[]{ll}|{\bf j}_{f}|&\text{Vicsek dynamics}\\ &\\ \frac{1}{2}|{\bf j}_{f}|^{2}&\text{DFL dynamics}.\end{array}\right.$$

$$\frac{d}{dt}\mathcal{F}=-\sigma {\int }_{\omega }f{\left|{\nabla }_{\omega }\ln \left(\frac{f}{M_f}\right)\right|}^2\mathrm{d}\omega \le 0.$$

$${\int }_{\omega }{\partial }_{t}f\ln \left(\frac{f}{M_f}\right)\mathrm{d}\omega$$

$$|{\mathbf{j}}_{f+h}{|}^2$$

$${\int }_{\omega }({\partial }_{t}f)\varphi (\omega )\mathrm{d}\omega ={\int }_{\omega }({\partial }_{t}f)\frac{\delta \Phi }{\delta f}(\omega )\mathrm{d}\omega =\frac{\text{d}}{\text{d}t}\Phi (f(t)).$$

$${\int }_{\omega }{\partial }_{t}f\ln \left(\frac{f}{M_f}\right)\mathrm{d}\omega$$

$$\mathcal{A}=\left\{f\in L^2({\mathbb{S}}^{n-1})|{\int }_{{\mathbb{S}}^{n-1}}\omega f(\omega )\mathrm{d}\omega =\mathbf{j}\text{and}{\int }_{{\mathbb{S}}^{n-1}}f(\omega )\mathrm{d}\omega =1\right\}.$$

$$\underset{\mathcal{A}}{\inf }\left\{{\int }_{\omega }f\ln f\mathrm{d}\omega \right\}={\int }_{\omega }{M}_{\ast }\ln M_{\ast }\mathrm{d}\omega ,$$

$$H[f]=\int f\ln f,\alpha [f]={\left|{\int }_{{\mathbb{S}}^{n-1}}\omega f-\mathbf{j}\right|}^2-1,\beta [f]={\left({\int }_{{\mathbb{S}}^{n-1}}f-1\right)}^2.$$

$${\frac{\delta H}{\delta f}|}_{f_{\ast }}$$

$$\Rightarrow \ln f_{\ast }+1$$

$$f_{\ast }(\omega )=C{\mathrm{e}}^{-2{\lambda }_1\omega \cdot \mathbf{j}},$$

Full text

1 Introduction

Swarming dynamics have attracted a lot attention in recent years raising the question how simple interaction rules could lead to complex pattern formation vicsek_collective_2012 . One of the main difficulty is to link the individual behaviors of agents and the pattern formations observed at a larger scale. Fortunately the framework of kinetic equations allows such transition between microscopic and macroscopic dynamics. Among the many existing swarming models camazine_self-organization_2001 , the Vicsek model vicsek_novel_1995 is one of the most popular since it describes a rather simple dynamics (there is only alignment) with few parameters but it is however able to generate complex pattern which are challenging to predict analytically. The Vicsek model has been well studied both numerically nagy_new_2007 ; chate_modeling_2008 ; gregoire_onset_2004 or analytically gamba_global_2016 ; bolley_mean-field_2012 ; figalli_global_2018 and the derivation of its kinetic and macroscopic equation is well-understood degond_hydrodynamic_2013 ; degond_continuum_2008 . However, as noted first by Chaté and Grégoiregregoire_onset_2004 , there exists a certain regime where the Vicsek model leads to the formation of traveling bands. Many numerical studies have been conducted to better analyze the formation of these bands at the particle level but none have been proposed so far to study the bands using a kinetic or macroscopic framework. This manuscript aims at proposing a first study on such band formation from the angle of kinetic equations.

After the discovery of band formation in the Vicsek dynamics by Chaté and Grégoire gregoire_onset_2004 , there has been a debate aldana_phase_2003 ; chate_comment_2007 about the order of the phase transition in the Vicsek model (continuous or discontinuous). As there was no analytic framework available, the conjecture could be only based on (particle) numerical simulations. However, the derivation of kinetic and macroscopic equation for the Vicsek model degond_continuum_2008 ; degond_hydrodynamic_2013 indicated that in a dense regime of particles (the so-called moderately interacting particle oelschlager_law_1985 ), the Vicsek model has a continuous transition from order to disorder. In this regime of high density, no phase transition or band formation could be observed. A major discovery was then provided by Degond, Frouvelle and Liu frouvelle_dynamics_2012 ; degond_macroscopic_2013 ; degond_phase_2015 where a modification of the (continuous) Vicsek was considered: alignment is proportional to the local density rather than to the mean direction. In their dynamics, a phase transition occurs: at low density, the velocity distribution becomes uniform, whereas at large density, the dynamics converge to a so-called von Mises distribution. This analytic result was only proven in an homogeneous setting (no spatial variable). Thus, it is still unknown what effect a transport term would have on the dynamics. This is however a very challenging question as the transport term breaks the entropy dissipation. In this manuscript, we propose to investigate numerically the Degond-Frouvelle-Liu (DFL) dynamics in a non-homogeneous setting.

Starting from the kinetic equation associated with the Vicsek model, we first review some properties of the collisional operator (entropy dissipation) that will be central for the building of our numerical scheme. Most of these estimates are built on the Fokker-Planck structure of the operator. We do take advantage of this formulation in the design of our numerical scheme. The key properties of the collision operator (positivity preserving, entropy dissipation) are also satisfied by the discrete operator. Since we aim at analyzing the long-time behavior of the solution, it is essential to preserve these properties. For instance, several papers have already proposed to solve numerically the kinetic equation associated with the Vicsek model using other methods (spectral method gamba_spectral_2015 , particle method dimarco_self-alignment_2016 , discontinuous Galerkin filbet_discontinuous-galerkin_2017 ). But we would rather have lower accuracy and a preserving numerical scheme to study the long-time behavior of the solution (even though our scheme is still second order accurate in the velocity-variable). We then explore the dynamics of the kinetic equation in various regimes. In the original Vicsek model, no band formations are observed, the spatial density becomes homogeneous while the velocity distribution becomes distributed according to a (global) von Mises distribution. In the Degond-Frouvelle-Liu (DFL) dynamics, however, when the density is above a threshold, band formation occurs starting from random initial configuration. As far as the authors know, this is the first time such band formations are observed at the kinetic level. Band were also observed in filbet_discontinuous-galerkin_2017 but there were only ’transient’, the density profile would become flat after a long time. Here, the density profile is not flattening out, but instead is becoming more and more concentrated. Numerically, we have to introduce an adaptive time step to deal with a demanding CFL condition because of this very phenomenon.

Although our numerical investigation suggests that band formation emerges from the DFL dynamics, it would be crucial to also develop an analytic framework to further understand this phenomenon. Our results indicate that the transport operator could further amplify the concentration originating from the alignment operator. From these observations, it seems unlikely that there exists an analytic profile for these bands. But the question remains open. Similarly, we could perform simulation in dimension 3, but the discretization of the unit sphere $\mathbb{S}^{2}$ is more delicate than $\mathbb{S}^{1}$ (there is no ’uniform grid’ on $\mathbb{S}^{2}$) and thus having discrete entropy dissipation or symmetry preserving would be more challenging. Finally, higher order accuracy in time should also be investigated using for instance gottlieb_strong_2001 ; carrillo_finite-volume_2015 .

2 Microscopic description

2.1 Vicsek model

The Vicsek model vicsek_novel_1995 ; degond_continuum_2008 at the particle level describes the motion of $N$ particles with position ${\bf x}_{i}\in\mathbb{R}^{d}$ (with $d=2,3$) and a direction $\omega_{i}\in\mathbb{S}^{d-1}$ (i.e. $|\omega_{i}|=1$). The evolution of the particles is given by the following system:

[TABLE]

where $c>0$ is the speed of the particle, $\mu$ is the strength of the alignment interaction, $\sigma$ is the intensity of the noise, $\mathrm{d}B_{i}^{t}$ are independent white noises, $P_{\omega_{i}^{\perp}}$ is the orthogonal projection on the orthogonal of $\omega_{i}$ defined as:

[TABLE]

which ensures that $|\omega_{i}(t)|=1$ over time, and $\Omega_{i}$ is the average direction of the particle $i$:

[TABLE]

where $R$ is the radius of interaction.

2.2 Degond-Frouvelle-Liu (DFL) dynamics

Degond, Frouvelle and Liu frouvelle_dynamics_2012 ; degond_phase_2015 proposed a modification of the dynamics where the alignment interaction $\mu$ is proportional to the norm of the flux ${\bf j}_{i}$:

[TABLE]

This modification has several consequences: i) unlike the Vicsek dynamics (1) which is ill-defined when the flux ${\bf j}_{i}$ equal zero ($\Omega_{i}$ not defined), the DFL dynamics does not have any singularity, ii) there is a phase transition in the dynamics (3) as the number of particles increases (or similarly as $\mu$ increases). The kinetic description of this dynamics will allow to better explain this phase transition (see section 3.2).

2.3 Band formation

Band formations have been first analyzed by Grégoire and Chaté gregoire_onset_2004 in the case of the original discrete Vicsek model and several numerical studies have been conducted since chate_modeling_2008 ; nagy_new_2007 ; aldana_phase_2003 . To motivate our study, we present numerically an example of such band formation in the context of the continuous dynamics (1).

The numerical simulation presented in this subsection is performed with $N=30,000$ particles on a square domain with length $L=4$ and periodic boundary condition. Initially, particles are distributed at random in space and velocity. Table 1 gives the list of values for the parameters. We observe in Figure 1 the formation of a traveling wave moving in the $x$-direction. To further quantify this formation, we estimate the average density $\rho$ and velocity $u$ in the $x$-direction:

[TABLE]

where $x_{i}$ and $\cos\theta_{i}$ are respectively the $x$-component of the position vector ${\bf x}_{i}$ and velocity $\omega_{i}$. We give an example of such $\rho$ and $u$ in Figure 2.

We notice that the regime in which the band formation occurs is far from being dense. Indeed, in a homogeneous setting, the average number of neighbors is given by:

[TABLE]

therefore we are far from being in kinetic regime (let alone macroscopic region). Thus, the validity of the kinetic equation associated with the dynamics (described in the next section) is questionable in this regime. Particles are not necessarily ’moderately interacting’ oelschlager_law_1985 .

3 Kinetic description

3.1 Introduction

The kinetic equation associated with the Vicsek dynamics (1) is described through the density distribution $f({\bf x},\omega,t)$. As the number of particles $N$ tends to infinity, the particle dynamics converge to the solution of a deterministic equation given by bolley_mean-field_2012

[TABLE]

where $c>0$ is the speed of the particles, $\mu>0$ is the intensity of the relaxation toward the mean velocity, $\sigma>0$ is the diffusion coefficient, $P_{\omega^{\perp}}$ the projection operator (2), $\Omega$ is the mean velocity at the point ${\bf x}$

[TABLE]

$R>0$ being the radius of interaction.

The DFL dynamics lead to a similar kinetic equation except for the transport term in $\omega$:

[TABLE]

in other words the strength $\mu$ of alignment is now proportional to $|{\bf j}_{f}|$.

3.2 Homogeneous case

To investigate kinetic equations, we study the homogeneous case, assuming that $f$ is independent of ${\bf x}$. Thus, the kinetic equations (4) and (5) become:

[TABLE]

with:

[TABLE]

and $\Omega_{f}=\frac{{\bf j}_{f}}{|{\bf j}_{f}|}$ with ${\bf j}_{f}=\int_{\omega\in\mathbb{S}^{d-1}}\omega f(\omega)\,\mathrm{d}\omega$ and

[TABLE]

The operator $Q$ defined by (7) can be written as a Fokker-Planck type operator. Introducing:

[TABLE]

with $\langle,\rangle$ the usual scalar product in $\mathbb{R}^{n}$, we find:

[TABLE]

using the identity $\nabla_{\omega}\langle{\bf u},\omega\rangle=P_{\omega^{\perp}}({\bf u})$. We deduce a first identity:

[TABLE]

Unfortunately, the left-hand side of (10) cannot be written as a total time derivative and thus we cannot deduce any entropy decay. The trick is to notice the following:

[TABLE]

Therefore,

[TABLE]

Thanks to the property of the logarithm, the left-hand side can now be written as a total time derivative and we deduce the following proposition.

Proposition 1

Suppose $f$ is a solution to the homogeneous kinetic equation (6) and consider the free energy:

[TABLE]

with:

[TABLE]

It satisfies:

[TABLE]

Proof

We remark that left-hand side of (11) is a total derivative:

[TABLE]

using the conservation of mass $\int_{\omega}\partial_{t}f\,\mathrm{d}\omega=0$. Notice moreover that $\phi$ (8) can be expressed as gradient (making the dynamics (7) a gradient flow as noted in figalli_global_2018 ). Indeed, taking $f+{h}$ a small perturbation of $f$, we have:

[TABLE]

thus $\frac{\delta|{\bf j}_{f}|^{2}}{\delta f}(\omega)=2\langle{\bf j}_{f},\omega\rangle$. In particular $\phi=\frac{\delta\Phi}{\delta f}$ with $\Phi$ given by (13). We deduce:

[TABLE]

Therefore, we obtain:

[TABLE]

with $\mathcal{F}$ given by (12).

3.3 Phase transition

Since the dynamics (7) have entropy, one can study the long-time behavior and deduce the convergence toward an equilibrium given as the minimizer of the free energy $\mathcal{F}$ (12). But first we need to identify the minimizers of $\mathcal{F}$. To do so, we notice that once we fix the flux ${\bf j}_{f}$, the minimizer would be given by von Mises.

Lemma 1

Fix ${\bf j}$ with $0<|{\bf j}|<1$ and consider the affine space:

[TABLE]

Then:

[TABLE]

where $M_{*}$ is the von Mises (9) satisfying $\int_{\omega}\omega M_{*}(\omega)\,\mathrm{d}\omega={\bf j}$.

Proof

Assume there exists a minimizer $f_{*}$ and rewrite the constraint as

[TABLE]

Denote $\lambda_{1}$ and $\lambda_{2}$ the Lagrange multiplier associated with $f_{*}$:

[TABLE]

since $\int f_{*}=1$. Thus, taking the exponential leads to:

[TABLE]

and therefore $f_{*}$ is a von Mises distribution.

As a consequence of Lemma 1, we can restrict the search of minimizers of the free energy $\mathcal{F}$ on von Mises distributions. In Figure 3, we estimate numerically the entropy $\int_{\omega}M\ln M$ of von Mises distributions depending on their average velocity $|{\bf j}|=|\int_{\omega}\omega M|$ along with its approximation near ${\bf j}\approx 0$:

[TABLE]

We deduce that the free energy $\mathcal{F}$ for the Vicsek model will never have a minimum at ${\bf j}=0$ meaning that the uniform distribution is never stable. However, for the DFL dynamics, when the diffusion $\sigma$ is large, the free entropy $\mathcal{F}$ will be minimum at ${\bf j}=0$ and therefore the uniform distribution will become the stable equilibrium. These two situations are depicted in Figure 4.

4 Numerical scheme

Several schemes have already been proposed to study the kinetic equation (4) using spectral method gamba_spectral_2015 , discontinuous Galerkin filbet_discontinuous-galerkin_2017 or semi-Lagrangian dimarco_self-alignment_2016 . However, we are now interested in the long time behavior of the solution, thus we would like to design a numerical scheme with several properties:

• •

conservative, preserve positivity (under some CFL condition)

• •

satisfy a discrete version of inequalities (10) and (11)

In the following, we study the 2D scenario, taking advantage of the fact that the velocity space $\omega\in\mathbb{S}^{1}$ can be parametrized using polar coordinates by $\theta\in\mathbb{R}/2\pi\mathbb{Z}$ with $\omega=\left(\begin{matrix}\cos\theta\\ \sin\theta\end{matrix}\right)$. Thus, the kinetic equation (4) becomes:

[TABLE]

where $\bar{\theta}$ is such that:

[TABLE]

and $\mu_{f}$ is either a constant (Vicsek model) or proportional to $|{\bf j}|$ (DFL dynamics). Our numerical scheme is then based on a splitting method solving separately:

• •

the transport part

[TABLE]

• •

the collision part :

[TABLE]

where $Q(f)=-\mu_{f}\partial_{\theta}(\sin(\bar{\theta}-\theta)f)+\sigma\partial_{\theta}^{2}f$.

4.1 Collision operator

In this section we focus on numerically solving equation (15). We write:

[TABLE]

where $M_{f}$ is the Von Mises distribution :

[TABLE]

and $C_{0}$ is a normalization constant. Here, we simply take $C_{0}=1$.

4.1.1 Discretization in $\theta$

Fix $N>0$ and consider a uniform discretization of the interval $[0,2\pi)$ with $\theta_{k}=i\Delta\theta$ and $\Delta\theta=\frac{2\pi}{N}$. Denote: $f_{k}=f(\theta_{k})$ and $f_{k+\frac{1}{2}}=f(\theta_{k}+\Delta\theta/2)$ and similarly $M_{k}=M_{f}(\theta_{k})$. To approximate $Q$ (which is a differential operator), we use the second order approximation:

[TABLE]

which gives

[TABLE]

with

[TABLE]

The discrete operator $Q_{N}$ can be identified with a square $N\times N$ matrix:

[TABLE]

and

[TABLE]

with a slight abuse of notation such as $M_{-\frac{1}{2}}=M_{N-\frac{1}{2}}$ (by the periodicity of $M_{f}$).

The discrete operator $Q_{N}$ reproduces many features of the differential operator $Q$. Define the scalar product:

[TABLE]

the operator $Q_{N}$ then satisfies the discrete equivalent to (10) and (11).

Proposition 2

The operator $Q_{N}$ (4.1.1) is symmetric with respect to this scalar product:

[TABLE]

and satisfies:

[TABLE]

[TABLE]

Proof

Take any vectors $u$ and $v$, then the Abel formula (discrete integration by parts) gives:

[TABLE]

From (4.1.1), we also deduce (17).

Moreover, using once again the Abel formula:

[TABLE]

Thus, denoting $x=\frac{u_{k+1}}{M_{k+1}}$ and $y=\frac{u_{k}}{M_{k}}$, we have an expression of the form:

[TABLE]

for any $x,y>0$. We deduce (18).

4.1.2 Explicit Euler

The Euler method can be used to discretize in time the collisional part of the kinetic equation (15):

[TABLE]

A sufficient condition for the $L^{\infty}$ stability of the scheme is to have the matrix $\text{Id}+\Delta tQ_{N}$ positive (i.e. all coefficients positive). This sufficient condition leads to the following CFL condition:

[TABLE]

which is usual for diffusion type operator. Moreover, if the CFL condition is met, then positivity and mass are preserved.

Remark 1

We can find an explicit sufficient condition to guarantee the CFL condition (19). Indeed, writing:

[TABLE]

where we have used the identity $\cos\alpha-\cos\beta=-2\sin\frac{\alpha+\beta}{2}\,\sin\frac{\alpha-\beta}{2}$. We deduce

[TABLE]

and find:

[TABLE]

This leads to the tractable (sufficient) CFL condition:

[TABLE]

4.1.3 Adaptative time step for the collision

One of the difficulties in computing an approximate solution to (5) is coping with the associated CFL condition (19). Indeed, the existence of a locally high $|j(f)|$ greatly decreases the right-hand side of (19), which penalizes the whole algorithm. Hence using a global CFL condition for the transport part and the collision part can lead to extremely long computation time.

We propose to decouple the time steps for the transport part and the collision equation at each time step, by using an adaptive method for the latter. Technically, we use the maximal time step associated with the CFL condition to solve the transport part (14), $\Delta t=\Delta x$, which incidentally has the advantage of minimizing numerical diffusion. Then, for each $(x_{i},y_{j})$, we consider (15) as a differential equation with final time $\Delta t$, which we solve by using the method described in Section 4.1 with a variable time step $\delta t$ that needs to be recomputed at each time $0\leq s\leq\Delta t$:

[TABLE]

This method also works for a constant relaxation $\mu(f)=\mu_{0}$, and can be preferred because it minimizes the numerical diffusion in the transport equation. Particularly when the constant $\sigma=\frac{\mu}{D}$ is large, in which case the collision CFL (20) is much smaller than the transport CFL (21).

A comparison between the errors done by the standard and adaptive schemes, respectively, is shown in Figure 5.

4.2 Numerical scheme for the transport operator

We use an upwind finite-difference method to solve the transport equation (14). More accurate schemes are available (such as WENO scheme shu_high_2009 ), however stability and positivity preserving are more essential in our study since we investigate the long-time behavior of the dynamics. We fix $M>0$ and consider a uniform discretization of the interval $[0,L)$ in $M$ points with $x_{i}=i\Delta x$, $y_{j}=j\Delta y$, and $\Delta x=\Delta y=\frac{L}{M}$. To discretize the kinetic equation, we use:

[TABLE]

and similarly for $\sin(\theta)\partial_{y}f$. Using this discretization, the standard Euler scheme gives as CFL condition:

[TABLE]

4.3 Summary full scheme

The full algorithm is finally a splitting between the transport and collision part. Notice that the time step $\Delta t$ should satisfy both CFL conditions (19) and (21). In general, the collisional CFL (19) is more restrictive. Therefore, the transport equation will be solved with a small CFL corresponding to large numerical viscosity. Since we aim at studying the large time behavior of the dynamics, this numerical viscosity might drastically change the outcome. Thus, we propose to use an adaptive method for the collisional operator. The idea is to simply iterate $K$ ‘small’ steps $\delta t=\Delta t/K$ to update the collision part choosing $K$ such that $\delta t$ satisfies the CFL condition (20).

5 Numerical experiments

5.1 Homogeneous case

To first investigate our numerical scheme, we study the homogeneous equation, thus solving only the collision operator (15). We present our numerical experiments with the Vicsek model (4), however results are similar with the DFL dynamics except that the time step $\Delta t$ may have to be adapted (since the CFL condition depends on $|{\bf j}|$ which varies over time).

As a first sanity check, we estimate the accuracy of the scheme. With this aim, we fix a final time $T=1$ and time step $\Delta t=.001$. Then, we vary the meshgrid in $\theta$, taking $\Delta\theta\in\{\frac{2\pi}{8},\frac{2\pi}{16},\dots,\frac{2\pi}{128}\}$ and estimate the $L^{2}$ error with the reference solution $f_{ref}$ computing with $\Delta\theta=\frac{2\pi}{256}$. For the initial condition, we use a smooth initial condition:

[TABLE]

We use a rather complicated expression to make sure that $f_{0}$ is non-symmetric. When $f_{0}$ is symmetric, the mean direction ${\bar{\theta}}$ is preserved over time, thus the Vicsek dynamics (6) becomes a linear evolution equation. Twisting the initial condition $f_{0}$ guarantees to have a fully non-linear equation.

In Figure 6-left, we plot the initial condition $f_{0}$ along with the reference solution $f_{ref}$ at $t=1$. The $L^{2}$ error for various discretizations is given in log scale in Figure 6-right. We observe that the error is decaying quadratically as expected.

Moreover, we also investigate the large-time behavior of the solution. First, we measure the evolution of the free entropy $\mathcal{F}$ over time and we observe that it is strictly decreasing (fig. 7-left). Second, we estimate the rate of convergence of $f(t)$ toward an equilibrium distribution. Using semi-log scale in fig. 7-right, we observe a linear decay indicating that the convergence is exponential.

5.2 Band formation

In the Vicsek model (4), we did not observe the formation of any bands. Rather, the dynamics always converge to a robust global alignment dynamics, where the spatial distribution (first moment of $f$) converges to a constant. The typical long-time behavior is represented in Figure 8. We postulate that the long-time behavior of this equation is just to converge to a uniform distribution of Von Mises equilibria to the homogeneous equation.

On the other hand, the DFL model (5) of interaction leads to the observation of bands. Typically, for a fixed set of parameters, we observed two different scenarios regarding the behavior of the local density $\rho$ and mean value ${\rho}$

[TABLE]

For a fixed mean value ${\rho}$, when the strength of interaction $\mu$ is small compared to the diffusion parameter $\sigma$ (i.e. $\mu\gg\sigma$), we observe that the solution converges to a uniform steady state. However, when $\mu\ll\sigma$, we observe the formation of bands as shown in Figure 9 (in which the x-axis has been reversed to provide better aesthetics). Thus, we retrieve an equivalent of the phase transition dynamics noticed by Frouvelle and Liu in frouvelle_dynamics_2012 . Those bands were first noticed starting from a random initial condition (Figure 9). Even though they literally emerge from chaos, they appear to be only meta-stable, as their small inhomogeneity in the direction perpendicular to the propagation amplifies slowly by attracting the neighbor particles and finally lead to high and localized concentrations as shown in Figure 9. At this point the computation is difficult to continue due to the extremely high computation times required by the CFL condition (20). Starting from an initial condition which is homogeneous in one direction (e.g. in $y$), however, the observed bands are very stable in time and can be kept alive for apparently an arbitrarily long time (the homogeneity being preserved by our scheme). Such a band is represented in Figure 9. The initial condition we used is the following:

[TABLE]

Bands were also observed in a modified one-dimensional model which we encoded to take advantage of the preservation of homogeneity in one direction. A resulting band is presented in Figure 9 and 9.

Numerical evidences show that the bands cannot be understood as traveling wave solutions to the kinetic equation (4), as one may believe at first sight. Indeed, there remains an inner motion inside the bands, that we can reveal by monitoring the maximal value of $\rho$ through time (see Figure 10). This reveals an asymptotically periodic behavior that strongly resembles the notion of pulsating fronts, which has been extensively studied in the context of reaction-diffusion equations xin_front_2000 . A deeper analytical understanding of this phenomenon is left for future work.

Finally, to strengthen the link between the phase transition and the formation of bands, we show in Figure 11 two kinds of entropy computed for a range of values of the diffusion coefficient $d$ and the mean value of the initial condition ${\bar{\rho}}$. Figure 11 represents the entropy of $f$ computed against the uniform distribution of the same mass:

[TABLE]

Figure 11 represents the generalized entropy of $f$ computed against the corresponding Von Mises distribution:

[TABLE]

where $M[{\bar{\rho}}](\theta)=2\pi{\bar{\rho}}\frac{\exp\left(\frac{\mu\kappa}{d}\cos(\theta)\right)}{\int_{0}^{2\pi}\exp\left(\frac{\mu\kappa}{d}\cos(\theta)\right)\mathrm{d}\theta}$ is the only candidate as stationary Von Mises distribution, $\kappa$ satisfying the compatibility condition

[TABLE]

Let us recall that the latter has only one solution $\kappa=0$ when $\sigma\geq\pi\mu{\bar{\rho}}$, and has exactly one positive solution when $\sigma<\pi\mu{\bar{\rho}}$ frouvelle_dynamics_2012 .

The match between the two plots suggests that the latter stationary state candidate is never stable except when $\kappa=0$. Indeed for $\sigma\geq\pi\mu{\bar{\rho}}$, the uniform distribution is stable for the homogeneous problem and $\kappa=0$; this corresponds to the top-left part of Figure 11, which suggests that this stability is transferred to the inhomogeneous problem (5). For $\sigma<\pi\mu{\bar{\rho}}$, however, we have $\kappa>0$ and the stable state for the homogeneous problem is described by the corresponding Von Mises distribution $M[{\bar{\rho}}](\theta)$; Figure 11 suggests that the inhomogeneous problem behaves otherwise, neither the uniform nor the homogeneous Von Mises distribution corresponding to the long-time behavior of the equation, except possibly in a very small area near $\sigma\approx\pi\mu{\bar{\rho}}$ (which appears more clearly in the log-plots 11 and 11). Instead, the instability of both homogeneous stationary states could be at the origin of the formation of bands.

6 Conclusion

In this manuscript, we have introduced a numerical scheme to solve kinetic equations related to the Vicsek model. Despite the additional difficulty to treat non-conservative dynamics, the numerical scheme is able to preserve positivity and entropy which allow to study the long-time behavior of the dynamics. Of particular interest in the formation of bands that were first observed at the microscopic model (i.e. particle simulations). The scheme is able to capture such band formation at the kinetic level.

It remains many open questions related to the kinetic equation of the Vicsek model. First, we still do not know if asymptotically there is an analytical description of the band. Having an Ansatz will help to study the convergence of the dynamics toward a band formation. However, the lack of entropy in the non-homogeneous case makes this task complicated. Another open problem is to investigate how the Vicsek model in a certain regime (i.e. low density, large coefficients) is able to create bands. We do not observe such band formation for the original Vicsek model at the kinetic level, bands only form with the DFL dynamics. This might indicate that the propagation of chaos, which links particle and kinetic formulations, is not valid in this regime.