The Function of Sleep

1.   Introduction

Flocking is observed in large populations of social agents such as birds ^[2], fish ^[3] or insects ^[4]. It refers to the emergence of large scale spatio-temporal structures which are not directly encoded in the individual agents' behavior. Understanding how large scale structures appear from individual behavior has sparked a huge literature in the recent years concerned with both modelling ^{[5,6,7,8,9,10]} and experiments ^[11,12]. We refer the reader to the reviews ^[13,14,15]. Yet this phenomenon is still poorly understood. The modelling of individual behavior at the microscopic scale requires to resolve the motion of each individual in the course of time. This leads to so-called Individual-Based Models (IBM) (or particle models) consisting of a huge number of coupled ordinary or stochastic differential equations. By contrast, given the very large number of agents involved, the modelling of the macroscopic scale is best done through Continuum Models (CM). They describe the system as a fluid through average quantities such as the density or mean velocity of the agents. At the mesoscopic scale, kinetic models (KM) are intermediates ketween IBM and CM. They describe the agents dynamics through the statistical distribution of their positions and velocities. In the literature, IBM, KM or CM are often chosen according to the authors' preferences rather than following an explicit rationale. However, understanding flocking requires understanding how the individuals' microscale impacts on the system's macroscale and consequently demands the use of a consistent sequence of IBM, KM and CM models. This consistency can only be guaranteed if the passage from IBM to KM and from KM to CM can be systematically established. But this issue is seldom considered. One goal of this paper is to establish the consistency between a KM and a CM of flocking based on a variant of the celebrated Cucker-Smale model of consensus (Eqs (1.1) and (1.2) below). KM have been introduced in the last years for the mesoscopic description of collective behavior of agents/particles with applications in collective behavior of cell and animal populations, see ^[1,16,17,18] and the references therein for a general overview on this active field. These models usually include alignment, attraction and repulsion as basic bricks of interactions between individuals.

In this paper, we will discard the microscopic scale, i.e. the IBM. We will directly consider the mescoscale, i.e. the KM and will focus on the derivation of the CM from the KM. Of course, our KM has an underlying IBM. However the techniques involved in the passage from IBM to KM are quite different from those needed to pass from KM to CM. For this reason, considering the IBM would have brought us beyond the scope of the present paper. We refer to ^{[19,20,21,22,23,24,25,26,27]} and references therein for a derivation of KM from IBM. The derivation of the CM from the KM requires a spatio-temporal rescaling. Indeed, the KM still describes-although in a statistical way compared to the IBM---the microscopic dynamics of the particles. In particular, it is written in space and time units that are of the order of the particle interaction distances and times. To describe the macroscopic scale, one needs to introduce a change of variables by which the space and time units become of the order of the system scale, which is much wider that the particle-related scales. This rescaling introduces a small parameter, the ratio of the microscopic to the macroscopic space units. The derivation of the CM from the KM consists in finding the limit of the KM when this small parameter tends to zero. This is the goal of this paper. We will see that the CM (Eqs (1.5) and (1.6) below) corresponding to the considered modified Cucker-Smale model is a system consisting of equations for the density and mean orientation of the particles as functions of space and time. This system, referred to as the `Self-Organized Hydrodynamics (SOH)' bears analogies with the Euler equations of isothermal compressible gas dynamics, with the important difference that the average velocity is replaced by the average orientation, i.e. a vector of norm one. The SOH appears in a variety of contexts related to alignment interactions, such as repulsion ^[28], nematic alignment ^[29], suspensions ^[30], solid orientation ^[31], and can be seen as a basic CM of collective dynamics.

Here, we focus on the derivation of macroscopic equations (SOH) for the collective motion of self-propelled particles with alignment and noise when a cruise speed for individuals is imposed asymptotically for large times as in ^{[32,33,34,35,36,37,38]}. More precisely, in the presence of friction and self-propulsion and the absence of other interactions, individuals/particles accelerate or break to achieve a cruise speed exponentially fast in time. The alignment between particles is imposed via localized versions of the Cucker-Smale or Motsch-Tadmor reorientation procedure ^{[17,39,40,41,42,43]} leading to relaxation terms to the mean velocity modulated or not by the density of particles. By scaling the relaxation time towards the asymptotic cruise speed, or equivalently, penalizing the balance between friction and self-propulsion, this alignment interaction leads asymptotically to variations of the classical kinetic Vicsek-Fokker-Planck equation with velocities on the sphere, see ^{[1,10,44,45,46,47,48]}. It was shown in ^[37,38] that particular versions of the localized kinetic Cucker-Smale model can lead to phase transitions driven by noise. Moreover, these phase transitions are numerically stable in this asymptotic limit converging towards the phase transitions of the limiting versions of the corresponding kinetic Vicsek-Fokker-Planck equation.

In this work, we choose a localized and normalized version of the Cucker-Smale model not showing phase transition. More precisely, let us denote by $ f = f(t, x, v) \geq 0 $ the particle density in the phase space $ (x, v) \in \mathbb{R}^d \times \mathbb{R}^d $, with $ d \geq 2 $. The standard self-propulsion/friction mechanism leading to the cruise speed of the particles in the absence of alignment is given by the term $ {\rm div}_v \{f (\alpha - \beta |v|^2)v\} $ with $ \alpha, \beta > 0 $, and the relaxation toward the normalized mean velocity writes $ {\rm div}_v \{ f(v - \Omega[f])\} $. Here, for any particle density $ f(x, v) $, the notation $ \Omega[f] $ stands for the orientation of the mean velocity

Notice that we always have

Let us remark that the standard localized Cucker-Smale model would lead to $ \rho {\rm div}_v \{ f(v -u[f])\} $ while the localized Motsch-Tadmor model would lead to $ {\rm div}_v \{ f(v -u[f])\} $. Our relaxation term towards the normalized local velocity $ \Omega [f] $ does not give rise to phase transition in the homogeneous setting on the limiting Vicsek-Fokker-Planck-type model on the sphere according to ^[45] and it produces a competition to the cruise speed term comprising a tendency towards unit speed. Including random Brownian fluctuations in the velocity variable leads to the kinetic Fokker-Planck type equation

We include this equation in a more general family of equations written in a compact form as

where

for any density distribution $ f $ with $ V $ a general confining potential in the velocity variables and $ \eta > 0 $ (see Lemma 2.1 for more information on the type of potentials that we consider). In the particular example considered above we take $ V = V_{\alpha, \beta} (|v|): = \beta \frac{|v|^4}{4} - \alpha \frac{|v|^2}{2}\; $.

We investigate the large time and space scale regimes of the kinetic tranport Eq (1.1) with collision operator given by (1.2). Namely, we study the asymptotic behavior when $ \varepsilon \to 0 $ of

supplemented with the initial condition

The rescaling taken in the kinetic transport Eq (1.3) with confining potential $ V_{ \alpha, \beta} $ can be seen as an intermediate scaling between the ones proposed in ^[48] and ^[37]. The difference being that we have a relaxation towards the normalized mean velocity $ \Omega[f] $ rather than the mean velocity $ u[f] $ as in ^[37,48]. This difference is important since in the first case there is no phase transition in the homogeneous limiting setting on the sphere as we mentioned above, while in the second there is, see ^[38,45,48]. In fact, in ^[48] the scaling corresponds to $ \eta = 1 / \varepsilon $ in (1.3), that is the relaxation to the cruise speed is penalized with a term of the order of $ 1 / \varepsilon^2 $. Whereas in ^[37] the scaling correponds to $ \eta = \varepsilon $, that is the cruise speed is not penalized at all.

The methodology followed in ^[48] lies within the context of measure solutions by introducing a projection operator onto the set of measures supported in the sphere whose radius is the critical speed $ r = \sqrt{\alpha/\beta} $. These technicalities are needed because the zeroth order expansion of $ f_\varepsilon $ lives on the sphere. This construction followed closely the average method in gyro-kinetic theory ^[49,50,51].

However, in our present case we will show in contrast to ^[37,48] that there are no phase transitions which is in accordance with the results obtained in ^[46] for the kinetic Vicsek-Fokker-Planck equation with analogous alignment operator on the sphere. A modified version of (1.1) and (1.2) in which phase-transitions occur was studied in ^[38] whose analysis is postponed to a future work to focus here on the mathematical difficulties of the asymptotic analysis. Another difference in the present case is that the zeroth order expansion of $ f_\varepsilon $ will be parameterized by Von Mises-Fisher distributions in the whole velocity space, that is $ f(t, x, v) = \rho (t, x) M_{\Omega(t, x)}(v) $, with $ \rho $ and $ \Omega $ being, respectively, the density and the mean orientation of the particles. And where for any $ \Omega \in \mathbb{S}^{d-1} $ we define (see section 2)

and

The main result of this paper is the asymptotic analysis of the singularly perturbed kinetic transport equation of Cucker-Smale type (1.3). The particle density $ \rho $ and the orientation $ \Omega $ obey the hydrodynamic type equations given in the following result.

Theorem 1.1. Let $ f ^{{\rm in}} \geq 0 $ be a smooth initial particle density with nonvanishing orientation at any $ x \in \mathbb{R}^d $. For any $ \varepsilon > 0 $ we consider the problem

with initial condition

At any $ (t, x) \in \mathbb{R}_+ \times \mathbb{R}^d $ the leading order term in the Hilbert expansion $ f_\varepsilon = f + {f_1}_{\varepsilon} + \ldots $ is an equilibrium distribution of $ Q $, that is $ f(t, x, v) = \rho (t, x) M_{\Omega(t, x)}(v) $ with $ M_{\Omega(t, x)}(v) $ defined in (1.4), where the concentration $ \rho $ and the orientation $ \Omega $ satisfy

with initial conditions

The constants $ c_1, c_2 $ are given by

and the function $ \chi $ solves

where $ e(c, r) = \exp(rc/\sigma) \exp(- (r^2 + 1)/(2\sigma) - V(r)/\sigma) $.

Our article is organized as follows. First, in section 2 we state auxiliary results allowing us to discuss the kernel of the collision operator. Then in section 3 we concentrate on the characterization of the collision invariants. We prove that the generalized collision invariants introduced in ^[1] coincide with the kernel of a suitable linearised collision operator. We explicitly describe the collision invariants in section 4 and investigate their symmetries. Finally, the limit fluid model is determined in section 5 and we analyse its main properties.

2.   Preliminaries

Plugging into (1.3) the Hilbert expansion

we obtain at the leading order

whereas to the next order we get

where $ \, \text{d} Q_f $ denotes the first variation of $ Q $ with respect to $ f $. The constraint (2.1) leads immediately to the equilibrium

where

Indeed, by using the identity

we can recast the operator $ Q $ as

We denote by $ \mathbb{S}^{d-1} $ the set of unit vectors in $ \mathbb{R}^d $. For any $ \Omega \in \mathbb{S}^{d-1} $, we consider the weighted spaces

and

The nonlinear operator $ Q $ should be understood in the distributional sense, and is defined for any particle density $ f = f(v) $ in the domain

We introduce the usual scalar products

and we denote by $ |\cdot |_{M_{\Omega}}, \|\cdot \|_{M_{\Omega}} $ the associated norms. We make the following hypotheses on the potential $ V $. We assume that for any $ \Omega \in \mathbb{S}^{d-1} $ we have

Clearly (2.5) holds true for the potentials $ V_{\alpha, \beta} $. Notice that in that case $ 1 \in L^2_{M_{\Omega}} $ and $ |1|_{M_{\Omega}} = 1 $ for any $ \Omega \in \mathbb{S}^{d-1} $. Moreover, we need a Poincaré inequality, that is, for any $ \Omega \in \mathbb{S}^{d-1} $ there is $ \lambda_\Omega > 0 $ such that for all $ \chi \in H^1_{M_{\Omega}} $ we have

A sufficient condition for (2.6) to hold comes from the well-known equivalence between the Fokker-Planck and Schrödinger operators (see for instance ^[52]). Namely, for any $ \Omega \in \mathbb{S}^{d-1} $ we have

The operator $ \mathcal{H}_\Omega = - \sigma \Delta_v + \left [\frac{1}{4\sigma} |\nabla_v \Phi_\Omega |^2 - \frac{1}{2} \Delta_v \Phi_\Omega \right] $ is defined in the domain

Using classical results for Schrödinger operators (see for instance Theorem ⅩⅢ.67 in ^[53]), we have a spectral decomposition of the operator $ \mathcal{H}_\Omega $ under suitable confining assumptions.

Lemma 2.1. Assume that for $ \Phi_\Omega $ defined in (2.3) the function $ v \to \frac{1}{4\sigma} |\nabla_v \Phi_\Omega |^2 - \frac{1}{2}\Delta_v \Phi_\Omega $ satisfies the following:

a) it belongs to $ L^1_{{\rm loc}}(\mathbb{R}^d) $,

b) it is bounded from below,

c)

Then $ \mathcal{H}_\Omega ^{-1} $ is a self adjoint compact operator in $ L^2(\mathbb{R}^d) $ and $ \mathcal{H}_\Omega $ admits a spectral decomposition, that is a nondecreasing sequence of real numbers $ (\lambda_\Omega ^n)_{n\in \mathbb{N}} $, $ \lim_{n \to \infty} \lambda_\Omega ^n = \infty $, and a $ L^2(\mathbb{R}^d) $-orthonormal basis $ (\psi_\Omega ^n)_{n \in \mathbb{N}} $ such that $ \mathcal{H}_\Omega \psi_\Omega ^n = \lambda_\Omega ^n \psi_\Omega ^n, n \in \mathbb{N} $, $ \lambda_\Omega ^0 = 0 $, $ \lambda_\Omega ^1 > 0 $.

Let us note that the spectral gap of the Schrödinger operator $ \mathcal{H}_\Omega $ is the Poincaré constant in the Poincaré inequality (2.6). Notice also that the hypotheses in Lemma 2.1 are satisfied by the potentials $ V_{\alpha, \beta} $, and therefore (2.6) holds true in that case. It is easily seen that the set of equilibrium distributions of $ Q $ is parametrized by $ d $ parameters as stated in the following result.

Lemma 2.2. Let $ f = f(v)\geq 0 $ be a function in $ D(Q) $. Then $ f $ is an equilibrium for $ Q $ if and only if there are $ (\rho, \Omega) \in \mathbb{R}_+ \times \left (\mathbb{S}^{d-1} \cup \{0\} \right) $ such that $ f = \rho M_{\Omega} $. Moreover we have $ \rho = \rho[f]: = \int_{\mathbb{R} ^d} \!{f(v)} \; {\rm d}v $ and $ \Omega = \Omega[f] $.

Proof. If $ f $ is an equilibrium for $ Q $, we have

and therefore there is $ \rho \in \mathbb{R} $ such that $ f = \rho M_{\Omega}f $. Obviously $ \rho = \int_{\mathbb{R} ^d} \!{f(v)} \; {\rm d}v\geq 0 $ and $ \Omega [f] \in \mathbb{S}^{d-1} \cup \{0\} $. Conversely, we claim that for any $ (\rho, \Omega)\in \mathbb{R}_+ \times (\mathbb{S}^{d-1} \cup \{0\}) $, the particle density $ f = \rho M_{\Omega} $ is an equilibrium for $ Q $. Indeed, we have

We are done if we prove that $ \Omega [f] = \Omega $. If $ \Omega = 0 $, it is easily seen that

implying $ \Omega[f] = 0 = \Omega $. Assume now that $ \Omega \in \mathbb{S}^{d-1} $. For any $ \xi \in \mathbb{S}^{d-1} $, $ \xi \cdot \Omega = 0 $, we consider the orthogonal transformation $ \mathcal{O}_\xi = I_d - 2\xi \otimes \xi $. Thanks to the change of variable $ v = \mathcal{O}_\xi v^{\prime} $, we write

where we have used the radial symmetry of $ V $, $ \mathcal{O}_\xi \xi = - \xi $ and $ \mathcal{O}_\xi \Omega = \Omega $. We deduce that $ \int_{\mathbb{R} ^d} \!{f(v)v} \; {\rm d}v = \int_{\mathbb{R} ^d} \!{(v\cdot \Omega) f(v)} \; {\rm d}v\; \Omega $. We claim that $ \int_{\mathbb{R} ^d} \!{(v\cdot \Omega) f(v)} \; {\rm d}v > 0 $. Indeed we have

Obviously, we have for any $ v \in \mathbb{R}^d $ such that $ v \cdot \Omega > 0 $

implying that $ \int_{\mathbb{R} ^d} \!{(v\cdot \Omega) f(v)} \; {\rm d}v > 0 $ and

3.   Characterization of the collision invariants

In ^[1], the following notion of generalized collision invariant (GCI) has been introduced.

Definition 3.1. (GCI)

Let $ \Omega \in \mathbb{S}^{d-1} $ be a fixed orientation. A function $ \psi = \psi (v) $ is called a generalized collision invariant of $ Q $ associated to $ \Omega $, if and only if

for all $ f $ such that $ (I_d - \Omega \otimes \Omega) \int_{\mathbb{R} ^d} \!{f(v)v} \; {\rm d}v = 0 $, that is such that $ \int_{\mathbb{R} ^d} \!{f(v)v } \; {\rm d}v \in \mathbb{R} \Omega $.

In order to obtain the hydrodynamic limit of (1.3), for any fixed $ (t, x) \in \mathbb{R}_+ \times \mathbb{R}^d $, we multiply (2.2) by a function $ v \to \psi_{t, x} (v) $ and integrate with respect to $ v $ yielding

The above computation leads naturally to the following extension of the notion of collision invariant, see also ^[48].

Definition 3.2. Let $ f = f(v)\geq 0 $ be an equilibrium of $ Q $. A function $ \psi = \psi (v) $ is called a collision invariant for $ Q $ associated to the equilibrium $ f $, if and only if $ \mathcal{L}_f ^\star \psi = 0 $, that is

We are looking for a good characterization of the linearized collision operator $ \mathcal{L}_f $ and its adjoint with respect to the leading order particle density $ f $. Motivated by (2.1), we need to determine the structure of the equilibria of $ Q $ which are given by Lemma 2.2.

By Lemma 2.2, we know that for any $ (t, x) \in \mathbb{R}_+ \times \mathbb{R}^d $, there are $ (\rho (t, x), \Omega (t, x)) \in \mathbb{R}_+ \times (\mathbb{S}^{d-1} \cup \{0\}) $ such that $ f(t, x, \cdot) = \rho (t, x) M_{\Omega (t, x)} $, where

The evolution of the macroscopic quantities $ \rho $ and $ \Omega $ follows from (2.2) and (3.1), by appealing to the moment method ^{[54,55,56,57,58,59]}. Next, we explicitly determine the linearization of the collision operator $ Q $ around its equilibrium distributions. For any orientation $ \Omega \in \mathbb{S}^{d-1} \cup \{0\} $ we introduce the pressure tensor

and the quantity

We will check later, see Lemma 3.2, that the pressure tensor $ {\mathcal{M}}_{\Omega} $ is symmetric.

Proposition 3.1. Let $ f = f(v) \geq 0 $ be an equlibrium distribution of $ Q $ with nonvanishing orientation, that is

(1) The linearization $ \mathcal{L}_f = \, \mathit{\text{d}} Q_f $ is given by

where $ P_f : = I_d - \Omega[f] \otimes \Omega [f] $ is the orthogonal projection onto $ \{\xi \in \mathbb{R}^d : \; \xi \cdot \Omega[f] = 0\} $. In particular $ \mathcal{L}_{\rho M_{\Omega}} = \mathcal{L}_{M_{\Omega}} $.

(2) The formal adjoint of $ \mathcal{L}_f $ is given by

(3) We have the identity

Note that $ {\rm div}_v $ refers to the divergence operator acting on matrices defined as applying the divergence operator over rows.

Proof.

(1) By standard computations we have

and

Therefore we obtain

(2) We have

implying

(3) For any $ i \in \{1, ..., d\} $ we have

and therefore, since $ f = \rho M_\Omega $ satisfies $ \sigma \nabla_v f + f \nabla_v \Phi_\Omega = 0 $, we get

Notice that at any $ (t, x) \in \mathbb{R}_+ \times \mathbb{R}^d $, the function $ h = 1 $ is a collision invariant for $ Q $, associated to $ f(t, x, \cdot) $. Indeed, for any $ g = g(v) $ we have

implying that $ \int_{\mathbb{R} ^d} \!{\; (\mathcal{L}_{f(t, x, \cdot)}g)(v)} \; {\rm d}v = 0 $ and therefore $ \mathcal{L}_{f(t, x, \cdot)} ^\star 1 = 0 $, for all $ (t, x) \in \mathbb{R}_+ \times \mathbb{R}^d $. Once we have determined a collision invariant $ \psi = \psi (t, x, v) $ at any $ (t, x) \in \mathbb{R}_+ \times \mathbb{R}^d $, we deduce, thanks to (3.1), a balance for the macroscopic quantities $ \rho (t, x) = \rho[f(t, x, \cdot)] $ and $ \Omega (t, x) = \Omega[f(t, x, \cdot)] $, given by the relationship

When taking as collision invariant the function $ h(t, x, v) = 1 $, we obtain the local mass conservation equation

As usual, we are looking also for the conservation of the total momentum, however, the nonlinear operator $ Q $ does not preserve momentum. In other words, $ v $ is not a collision invariant. Indeed, if $ f = \rho M_{\Omega} $ is an equilibrium with nonvanishing orientation, we have

and therefore $ v $ is not a collision invariant.

We concentrate next on the resolution of (3.2). We will use the notation $ \partial_v \xi = \left(\tfrac{\partial \xi_i}{\partial v_j}\right) $ for the Jacobian matrix of a vector field $ \xi $ and $ {\rm div}_v $ for the divergence operator in $ v $ of both vectors and matrices with the convention of taking the divergence over the rows of the matrix. With this convention, we have

for all smooth functions $ g $, vector fields $ \xi, \eta $, and matrices $ A $. We now focus in finding a parameterization of the kernel of the operator $ \mathcal{L}_f ^\star $.

Lemma 3.1. Let $ f = \rho M_{\Omega} $ be an equilibrium of $ Q $ with nonvanishing orientation. The following two statements are equivalent:

(1) $ \psi = \psi (v) $ is a collision invariant for $ Q $ associated to $ f $.

(2) $ \psi $ satisfies

for some vector $ W \in \ker ({\mathcal{M}}_{\Omega} - \sigma c_1 I_d) $.

Moreover, the linear map $ W : \ker (\mathcal{L}_f ^\star) \to \ker ({\mathcal{M}}_{\Omega} - \sigma c_1 I_d) $, with $ W[\psi] : = \int_{\mathbb{R} ^d} \!{M_{\Omega} (v) \nabla_v \psi} \; {\rm d}v / c_1 $ induces an isomorphism between the vector spaces $ \ker (\mathcal{L}_f ^\star)/\ker W $ and $ \ker ({\mathcal{M}}_{\Omega} - \sigma c_1 I_d) $, where $ \ker W $ is the set of the constant functions.

Proof.

(1)$ \implies $(2) Since $ \psi $ is a collision invariant associated to $ f $, i.e. $ \mathcal{L}_f ^\star \psi = 0 $, and by the third statement in Proposition 3.1 we deduce (using also the first formula in (3.5) with $ f {\mathcal{M}}_{\Omega} /c_1 $ and $ \psi $)

Note that if $ \rho = 0 $ then $ f = 0 $ and $ \int v f \, \text{d} v = 0 $, implying that $ \Omega [f] = 0 $. Hence, since $ \Omega \neq 0 $, we have $ \rho > 0 $ and thus $ W[\psi] \in \ker ({\mathcal{M}}_{\Omega} - \sigma c_1 I_d) $, saying that (3.6) holds true with $ W = W[\psi] \in \ker({\mathcal{M}}_{\Omega} - \sigma c_1 I_d) $.

(2)$ \implies $(1) Let $ \psi $ be a function satisfying (3.6) for some vector $ W \in \ker ({\mathcal{M}}_{\Omega} - \sigma c_1 I_d) $. Multiplying (3.6) by $ f(v-\Omega) $ and integrating with respect to $ v $ yields (thanks to the second formula in (3.5))

which implies $ W[\psi] = W $ since $ {\mathcal{M}}_{\Omega} W = \sigma c_1 W $ by the assumption $ W \in \ker ({\mathcal{M}}_{\Omega} - \sigma c_1 I_d) $. Therefore $ \psi $ is a collision invariant for $ Q $, associated to $ f $

Remark 3.1. For any non negative measurable function $ \chi = \chi (c, r) :]-1, +1[\times] 0, \infty[\to \mathbb{R} $ and any $ \Omega \in \mathbb{S}^{d-1} $, for $ d \geq 2 $, we have

where $ |\mathbb{S}^{d-2}| $ is the surface of the unit sphere in $ \mathbb{R}^{d-1} $, for $ d \geq 3 $, and $ |\mathbb{S}^0| = 2 $ for $ d = 2 $. In particular we have the formula

where $ e(c, r) = \exp(rc/\sigma) \exp (- (r^2 + 1)/(2\sigma) - V(r)/\sigma) $.

Notice that thanks to (3.7) the coefficient $ c_1 $ does not depend upon $ \Omega \in \mathbb{S}^{d-1} $

In order to determine all the collision invariants, we focus on the spectral decomposition of the pressure tensor $ {\mathcal{M}}_{\Omega} $ for any $ \Omega \in \mathbb{S}^{d-1} $. In particular, the next lemma will imply the symmetry of the pressure tensor.

Lemma 3.2. (Spectral decomposition of $ {\mathcal{M}}_{\Omega} $) For any $ \Omega \in \mathbb{S}^{d-1} $ we have $ {\mathcal{M}}_{\Omega} = \sigma c_1 (I_d - \Omega \otimes \Omega) $. In particular we have $ \ker ({\mathcal{M}}_{\Omega} - \sigma c_1 I_d) = (\mathbb{R} \Omega)^\perp $ and thus $ {\rm dim} \left (\ker (\mathcal{L}_f ^\star) /\ker W \right) = {\rm dim} \ker ({\mathcal{M}}_{\Omega} - \sigma c_1 I_d) = d-1 $, cf. Lemma 3.1.

Proof. Let us consider $ \{E_1, \ldots, E_{d-1}\} $ an orthonormal basis of $ (\mathbb{R}\Omega) ^\perp $. By using the decomposition

one gets

We claim that the following equalities hold true

Formula (3.9) is obtained by using the change of variable $ v = (I_d - 2 E_j \otimes E_j) v^{\prime} $. It is easily seen that

and therefore we have

which implies (3.9). For the formulae (3.10) with $ i \neq j $, we appeal to the orthogonal transformation

Notice that $ \mathcal{O}_{ij} \xi = \xi $, for all $ \xi \in ({\rm span}\{E_i, E_j\}) ^\perp $, $ \mathcal{O}_{ij} E_i = - E_j $, $ \mathcal{O}_{ij} E_j = E_i $ and therefore

After this change of variable we deduce that

and also

Thanks to the equality $ \sum_{i = 1}^{d-1} (E_i \cdot v)^2 = |v|^2 - (v \cdot \Omega)^2 $, one gets

Coming back to (3.8) we obtain

We are done if we prove that

Notice that, using (2.4):

and therefore

By Lemmas 3.1 and 3.2 the computation of the collision invariants is reduced to the resolution of (3.6) for any $ W \in (\mathbb{R} \Omega) ^\perp $. Hence, if we denote by $ E_1, E_2, \ldots, E_{d-1} $ any orthonormal basis of $ (\mathbb{R} \Omega)^{\perp} $, we obtain a set of $ d-1 $ collision invariants $ \psi_{E_1}, \psi_{E_2}, \ldots, \psi_{E_{d-1}} $ for $ Q $ associated to the equilibrium distribution $ f $ such that $ E_i = W[\psi_{E_i}] $, $ i = 1, \dots, d-1 $. This set of collision invariants forms a basis for the $ \ker(\mathcal{L}_f ^\star) $. In the next section we will characterize this set of collision invariants and provide and easy manner to compute them (see Lemma 4.1).

We conclude this section by showing that in our case the set of all GCIs of the operator $ Q $ coincide with the kernel of the operator $ \mathcal{L}_f ^\star $.

Theorem 3.1. Let $ M_{\Omega} $ be an equilibrium of $ Q $ with nonvanishing orientation $ \Omega \in \mathbb{S}^{d-1} $. The set of collision invariants of $ Q $ associated to $ M_{\Omega} $ coincides with the set of the generalized collision invariants of $ Q $ associated to $ \Omega $.

Proof. Let $ \psi = \psi (v) $ be a generalized collision invariant of $ Q $ associated to $ \Omega $. We denote by $ \{e_1, \ldots, e_d\} $ the canonical basis of $ \mathbb{R}^d $. For any $ f = f(v) $ satisfying $ (I_d - \Omega \otimes \Omega) e_i \cdot \int_{\mathbb{R} ^d} \!{f(v) v } \; {\rm d}v = 0, 1\leq i \leq d $, we have

Therefore the linear form $ f \to \int_{\mathbb{R} ^d} \!{f(\sigma \Delta_v \psi - \nabla_v \Phi_\Omega \cdot \nabla_v \psi)} \; {\rm d}v $ is a linear combination of the linear forms $ f \to (I_d - \Omega \otimes \Omega) e_i \cdot \int_{\mathbb{R} ^d} \!{f(v) v } \; {\rm d}v $. We deduce that there is a vector $ \tilde{W} = (\tilde{W}_1, \ldots, \tilde{W}_d) \in \mathbb{R}^d $ such that

for any $ f $ and thus

implying that $ \psi $ satisfies (3.6) with the vector $ W = (I_d - \Omega \otimes \Omega)\tilde{W} \in (\mathbb{R} \Omega) ^\perp $, that is, $ \psi $ is a collision invariant of $ Q $ associated to $ M_{\Omega} $.

Conversely, let $ \psi = \psi (v) $ be a collision invariant of $ Q $ associated to $ M_{\Omega} $. By Lemmas 3.1 and 3.2 we know that there is $ W \in (\mathbb{R} \Omega)^\perp $ such that

Multiplying by any function $ f $ such that $ (I_d - \Omega \otimes \Omega) \int_{\mathbb{R} ^d} \!{f(v)v} \; {\rm d}v = 0 $ one gets

implying that $ \psi $ is a generalized collision invariant of $ Q $ associated to $ \Omega $.

4.   Identification of the collision invariants

In this section we investigate the structure of the collision invariants of $ Q $ associated to an equilibrium distribution $ f = \rho M_{\Omega} $. By Lemmas 3.1 and 3.2, we need to solve the elliptic problem

for any $ W \in (\mathbb{R}\Omega)^\perp $. We appeal to a variational formulation by considering the continuous bilinear symmetric form $ a_\Omega: H^1_{M_{\Omega}} \times H^1_{M_{\Omega}} \to \mathbb{R} $ defined as

and the linear form $ l : H^1_{M_{\Omega}} \to \mathbb{R} $, $ l (\theta) = \int_{\mathbb{R} ^d} \!{\theta(v) (v \cdot W) M_{\Omega} (v)} \; {\rm d}v, \theta \in H^1_{M_{\Omega}} $. Notice that $ l $ is well defined and bounded on $ H^1_{M_{\Omega}} $ provided that the additional hypothesis $ |v| \in L^2_{M_{\Omega}} $ holds true, that is

The above hypothesis is obviously satisfied by the potentials $ V_{\alpha, \beta} $. We say that $ \psi \in H^1_{M_{\Omega}} $ is a variational solution of (4.1) if and only if

Proposition 4.1. A necessary and sufficient condition for the existence and uniqueness of variational solution to (4.1) is

Proof. The necessary condition for the solvability of (4.1) is obtained by taking $ \theta = 1 $ (which belongs to $ H^1_{M_{\Omega}} $ thanks to (2.5)) in (4.2) leading to (4.3). This condition is satisfied for any $ W \in (\mathbb{R}\Omega)^\perp $ since we have

The condition (4.3) also guarantees the solvability of (4.1). Indeed, under the hypotheses (2.5) and (2.6), the bilinear form $ a_\Omega $ is coercive on the Hilbert space $ \tilde{H}^1_{M_{\Omega}} : = \{\chi \in H^1_{M_{\Omega}}\; :\; ((\chi, 1))_{M_{\Omega}} = 0\} $, i.e. we have:

Applying the Lax-Milgram lemma to (4.2) with $ \psi $, $ \theta \in \tilde{H}^1_{M_{\Omega}} $ yields a unique function $ \psi \in \tilde{H}^1_{M_{\Omega}} $ such that

Actually, the compatibility condition $ l(1) = 0 $ allows us to extend (4.4) to $ H^1_{M_{\Omega}} $. This follows by applying (4.4) with $ \tilde{\theta} = \theta - ((\theta, 1))_{M_{\Omega}} $, for $ \theta \in H^1_{M_{\Omega}} $. Moreover, the uniqueness of the solution for the problem on $ \tilde{H}^1_{M_{\Omega}} $ implies the uniqueness, up to a constant, of the solution for the problem on $ H^1_{M_{\Omega}} $.

As observed in (3.1), the fluid equations for $ \rho $ and $ \Omega $ will follow by appealing to the collision invariants associated to the orientation $ \Omega \in \mathbb{S}^{d-1} $, for any $ W \in (\mathbb{R} \Omega)^\perp $. When $ W = 0 $, the solutions of (4.1) are all the constants, and we obtain the particle number balance (3.4). Consider now $ W \in (\mathbb{R} \Omega)^\perp \setminus \{0\} $ and $ \psi_W $ a solution of (4.1). Obviously we have $ \psi_W = \tilde{\psi}_W + \int_{\mathbb{R} ^d} \!{\psi_W (v)M_{\Omega} (v)} \; {\rm d}v $, where $ \tilde{\psi}_W $ is the unique solution of (4.1) in $ \tilde{H}^1_{M_{\Omega}} $. It is easily seen, thanks to (3.4) and the linearity of (3.3), that the balances corresponding to $ \psi_W $ and $ \tilde{\psi}_W $ are equivalent. Therefore for any $ W \in (\mathbb{R} \Omega)^\perp $ it is enough to consider only the solution of (4.1) in $ \tilde{H}^1_{M_{\Omega}} $. From now on, for any $ W \in (\mathbb{R} \Omega)^\perp $, we denote by $ \psi_W $ the unique variational solution of (4.1) verifying $ \int_{\mathbb{R} ^d} \!{\psi_W (v) M_{\Omega} (v)} \; {\rm d}v = 0 $. The structure of the solutions $ \psi_W, W \in (\mathbb{R} \Omega)^\perp\setminus \{0\} $ comes by the symmetry of the equilibrium $ M_{\Omega} $. Analyzing the rotations leaving invariant the vector $ \Omega $, we prove as in ^[48] the following result.

Proposition 4.2. Consider $ W \in (\mathbb{R} \Omega)^\perp\setminus \{0\} $. For any orthogonal transformation $ \mathcal{O} $ of $ \mathbb{R}^d $ leaving $ \Omega $ invariant, that is $ \mathcal{O} \Omega = \Omega $, we have

where $ ^t \mathcal{O} $ denotes the transpose of the matrix $ \mathcal{O} $.

Proof. First of all notice that $ ^t \mathcal{O} W \in (\mathbb{R} \Omega) ^\perp \setminus \{0\} $. We know that $ \psi_W $ is the minimum point of the functional

It is easily seen that, for any orthogonal transformation $ \mathcal{O} $ of $ \mathbb{R}^d $ leaving the orientation $ \Omega $ invariant, and any function $ z \in \tilde{H}^1_{M_{\Omega}} $, we have, by defining $ z_{\mathcal{O}}: = z\circ \mathcal{O} \in \tilde{H}^1_{M_{\Omega}} $

Moreover, we obtain with the change of variables $ v' = \mathcal{O} v $ and using that $ M_\Omega (v) = M_\Omega (v') $:

Finally, we deduce that

for any $ z \in \tilde{H}^1_{M_{\Omega}} $, implying that $ \psi_W \circ \mathcal{O} = \psi_{^t \mathcal{O} W} $.

The computation of the collision invariants $ \{\psi_W\; :\; W \in (\mathbb{R} \Omega) ^\perp \setminus \{0\}\} $ can be reduced to the computation of one scalar function. For any orthonormal basis $ \{E_1, \ldots, E_{d-1}\} $ of $ (\mathbb{R} \Omega) ^\perp $ we define the vector field $ F = \sum_{i = 1}^{d-1} \psi_{E_i} E_i $. This vector field does not depend upon the basis $ \{E_1, \ldots, E_{d-1}\} $ and has the following properties, see ^[48].

Lemma 4.1. The vector field $ F $ does not depend on the orthonormal basis $ \{ E_1, \ldots, E_{d-1}\} $ of $ (\mathbb{R}\Omega)^\perp $ and for any orthogonal transformation $ \mathcal{O} $ of $ \mathbb{R}^d $, preserving $ \Omega $, we have $ F \circ \mathcal{O} = \mathcal{O} F $. There is a function $ \chi $ such that

and thus, for any $ i \in \{1, ..., d-1\} $, we have

Proof. Let $ \{F_1, \ldots, F_{d-1}\} $ be another orthonormal basis of $ (\mathbb{R} \Omega)^\perp $. The following identities hold

and therefore

For any orthogonal transformation of $ \mathbb{R}^d $ such that $ \mathcal{O} \Omega = \Omega $ we obtain thanks to Proposition 4.2

where, the last equality holds true since $ \{^t \mathcal{O} E_1, \ldots, \; ^t\mathcal{O} E_{d-1}\} $ is an orthonormal basis of $ (\mathbb{R} \Omega) ^\perp $. Let $ v \in \mathbb{R}^d \setminus (\mathbb{R} \Omega) $ and consider

Notice that $ E \cdot \Omega = 0, |E| = 1 $. When $ d = 2 $, since the vector $ F(v) $ is orthogonal to $ \Omega $, there exists a function $ \Lambda = \Lambda (v) $ such that

If $ d \geq 3 $, let us denote by $ ^\bot E $, any unitary vector orthogonal to $ E $ and $ \Omega $. Introducing the orthogonal matrix $ \mathcal{O} = I_d - 2 \; ^\bot E \otimes \; ^\bot E $ (which leaves $ \Omega $ invariant), we obtain $ F \circ \mathcal{O} = \mathcal{O} F $. Observe that

and thus

from which it follows that $ ^\bot E \cdot F(v) = 0 $, for any vector $ ^ \bot E $ orthogonal to $ E $ and $ \Omega $. Hence, there exists a function $ \Lambda (v) $ such that

We will show that $ \Lambda (v) $ depends only on $ v\cdot \Omega/|v| $ and $ |v| $. Indeed, for any $ d\geq 2 $, and any orthogonal transformation $ \mathcal{O} $ such that $ \mathcal{O} \Omega = \Omega $ we have $ F(\mathcal{O} v) = \mathcal{O} F(v) $, $ E (\mathcal{O} v) = \mathcal{O} E (v) $ because

implying that $ \Lambda (\mathcal{O} v) = \Lambda (v) $, for any $ v \in \mathbb{R}^d \setminus (\mathbb{R} \Omega) $. We are done if we prove that $ \Lambda (v) = \Lambda (v^{\prime}) $ for any $ v, v^{\prime} \in \mathbb{R}^d \setminus (\mathbb{R} \Omega) $ such that $ v\cdot \Omega /|v| = v^{\prime} \cdot \Omega /|v^{\prime}|, |v| = |v^{\prime}|, v \neq v^{\prime} $. It is enough to consider the rotation $ \mathcal{O} $ such that

The equality $ \mathcal{O} E = E ^{\; \prime} $ implies that $ \mathcal{O} v = v^{\prime} $ and therefore $ \Lambda (v^{\prime}) = \Lambda (\mathcal{O} v) = \Lambda (v) $, showing that there exists a function $ \chi $ such that $ \Lambda (v) = \chi (v \cdot \Omega/|v|, |v|), v \in \mathbb{R}^d \setminus \{0\} $.

In the last part of this section we concentrate on the elliptic problem satisfied by the function $ (c, r) \to \chi (c, r) $ introduced in Lemma 4.1. Even if $ \psi_{E_i} $ are eventually singular on $ \mathbb{R}\Omega $, it will be no difficulty to define a Hilbert space on which solving for the profile $ \chi $. We again proceed using the minimization of quadratic functionals.

Proposition 4.3. The function $ \chi $ constructed in Lemma 4.1 solves the problem

Proof. For any $ i \in \{1, \ldots, d-1\} $, let us consider $ \psi_{E_i, h} (v) = h (v \cdot \Omega /|v|, |v|) \frac{v \cdot E_i}{\sqrt{|v|^2 - (v \cdot \Omega)^2}} $ where $ v \to h (v \cdot \Omega /|v|, |v|) $ is a function such that $ \psi_{E_i, h} \in H^1_{M_{\Omega}} $. Observe that if $ h = \chi $, then $ \psi_{E_i, h} $ coincides with $ \psi_{E_i} $. Note that generally $ \psi_{E_i, h} $ are not collision invariants, but perturbations of them, corresponding to profiles $ h $. In this way, the minimization problem (4.7) will lead to a minimization problem on $ h $, whose solution will be $ \chi $. Notice that once that $ \psi_{E_i, h} \in H^1_{M_{\Omega}} $, then $ \int_{\mathbb{R} ^d} \!{\psi_{E_i, h}M_{\Omega} (v)} \; {\rm d}v = 0 $, saying that $ \psi_{E_i, h} \in \tilde{H}^1_{M_{\Omega}} $. We know that $ \psi_{E_i} $ is the minimum point of $ J_{E_i} $ on $ \tilde{H}^1_{M_{\Omega}} $ and therefore

A straightforward computation shows that

and

The condition $ \psi_{E_i, h} \in H^1_{M_{\Omega}} $ writes

which is equivalent, thanks to the Poincaré inequality (2.6) to

Based on formula (3.7), we have

and therefore we consider the Hilbert space $ H_d = \{h :\; ]-1, +1[\times] 0, \infty[\to \mathbb{R}, \|h \|_d ^2 < \infty \} $, endowed with the scalar product

for $ g $ and $ h $ in $ H_d $ and the norm given by

The expression $ J_{E_i}(\psi_{E_i, h}) $ writes as functional of $ h $

where

Coming back to (4.7) and using (4.5), we deduce that

Therefore, by the Lax-Milgram lemma, we deduce that $ \chi $ solves the problem (4.6).

5.   Hydrodynamic equations

After identifying the collision invariants, we determine the fluid equations satisfied by the macroscopic quantities entering the dominant particle density $ f(t, x, v) = \rho (t, x) M_{\Omega (t, x)} (v) $. As seen before the balance for the particle density follows thanks to the collision invariant $ \psi = 1 $. The other balances follow by appealing to the vector field $ F $ (cf. Lemma 4.1) and the details are given in section 5.1.

5.1.   Proof of Theorem 1.1

Applying (3.1) with $ \psi = 1 $ leads to (3.4). For any $ (t, x) \in \mathbb{R}_+ \times \mathbb{R}^d $ we consider the vector field

By the definition of $ F(t, x, \cdot) $, we know that

and therefore (3.1) implies

It remains to compute $ \int_{\mathbb{R} ^d} \!{\partial_t f F(t, x, v)} \; {\rm d}v $ and $ \int_{\mathbb{R} ^d} \!{v \cdot \nabla_x f F(t, x, v) } \; {\rm d}v $ in terms of $ \rho (t, x) $ and $ \Omega (t, x) $. By a direct computation we obtain

implying, thanks to the equalities $ \partial_t \Omega \cdot \Omega = 0 $ and $ v \cdot \partial_t \Omega = P_f v \cdot \partial_t \Omega $,

It is an easy exercise to show that the integral $ \int_{\mathbb{R} ^d} \!{\chi \left (\frac{v \cdot \Omega}{|v|}, |v|\right)M_{\Omega} (v) \frac{P_f v}{|P_f v |}} \; {\rm d}v $ vanishes and that the following relationship holds

Therefore, by taking into account that $ \partial_t \Omega \cdot \Omega = 0 $, the equality (5.2) becomes

Similarly we write (for any smooth vector field $ \xi (x) $, the notation $ \partial_x \xi $ stands for the Jacobian matrix of $ \xi $, i.e. $ (\partial_x \xi)_{ i, j} = \partial_{x_j} \xi_i $)

implying

where

Notice that in the above computations we have used $ (^t \partial_x \Omega) \Omega = 0 $ and

for any $ i, j, k \in \{1, \ldots, d-1\} $. Combining (5.1), (5.3) and (5.4) yields

or equivalently

where

5.2.   Properties of the hydrodynamic equations

Let us start by noticing that the system (1.5) and (1.6) is hyperbolic as a consequence of Theorem 4.1 in ^[60]. On the other hand, the orientation balance Eq (5.5) propagates the constraint $ |\Omega| = 1 $. Indeed, let $ \Omega = \Omega (t, x) $ be a smooth solution of (5.5), satisfying $ |\Omega (0, \cdot)| = 1 $. Multiplying by $ \Omega (t, x) $ we obtain

implying that $ |\Omega| $ is constant along the characteristics of the vector field $ c_2 \Omega (t, x) \cdot \nabla_x $ and thus $ |\Omega (t, \cdot)| = |\Omega (0, \cdot)| = 1 $, for all $ t \geq 0 $.

The rescaled Eq (1.3) can be considered as an intermediate model between the equations introduced in ^[37] and ^[48] when there are no phase transitions. In ^[48], the authors considered a strong relaxation towards the `terminal speed' (or cruise speed). Whereas in ^[37] the authors do not impose a penalization on the self-propelled/friction term. Our result could be applied to obtain the results in ^[48] without resorting to measures supported on the sphere by doing a double passage to the limit. First, passing to the limit $ \varepsilon \rightarrow 0 $ in (1.3), taking $ V = V_{\alpha, \beta} $, we obtain (1.5) and (1.6). Afterwards, we rescale $ V $ to $ \lambda \widetilde{V} $ in the system (1.5) and (1.6) and study the limit when $ \lambda \rightarrow \infty $. This amounts to study the behavior of the coefficients

and

when $ \lambda \rightarrow \infty $, where the function $ \chi_\lambda $ solves the elliptic problem

and $ e_\lambda (c, r) = \exp(rc/\sigma) \exp(- (r^2 + 1)/(2\sigma) - \lambda \widetilde{V} (r)/\sigma) $.

In order to analyse the asymptotic behavior of $ c_{1, \lambda} $ we introduce the following result.

Lemma 5.1. Let $ \varphi \in C^2 ((0, \infty); \mathbb{R}) $ and $ g \in C^0 ((0, \infty) \times (0, \pi); \mathbb{R}) $. Let us assume that

i) $ \int_{\mathbb{R}_+} \int_0^\pi \exp (\lambda \varphi (r)) | g (r, \theta)| \, \mathit{\text{d}} \theta \, \mathit{\text{d}} r < \infty $,

ii) The function $ \varphi $ has a unique global maximum at an interior point $ r_0 $,

iii) $ \int_0^\pi g (r_0, \theta) \, \mathit{\text{d}} \theta \neq 0 $.

Then the function $ G(\lambda) $ defined as

has the the following asymptotic behavior

as $ \lambda \rightarrow \infty $.

The proof of this result is a direct application of the Laplace method, see for instance ^[61]. As an immediate consequence of Lemma 5.1 we obtain

where $ r_0 $ is the minimum of the potential function $ V_{\alpha, \beta} (r) $. Let us note that the asymptotic study of the coefficient $ c_{1, \lambda} $ when $ \lambda \rightarrow \infty $ can also be performed using Lemma 5.1 for more general potentials than $ V_{\alpha, \beta} (r) $. In particular, we could also consider smooth potentials $ V (|v|) $ having a unique global minimum $ r_0 $ such that $ V' (r) < 0 $, for $ 0 < r < r_0 $, and $ V' (r) > 0 $, for $ r > r_0 $. On the other hand, the asymptotic study of $ c_{2, \lambda} $ could be addressed following similar techniques as in ^[60], however, we do not dwell upon this matter here and leave it for a future work.

6.   Conclusion

In this paper, we have considered a flocking model consisting of a modified Cucker-Smale involving noise, alignment and self-propulsion. We have investigated its macroscopic limit when the time and space variables are rescaled to the macroscopic scale. In this limit, the velocity distribution converges to an equilibrium which depends on the local density and local mean orientation of the particles. The density and mean orientation evolve in space and time according to a hyperbolic system of equations named the Self-Organized Hydrodynamics. This system is akin to the system of isothermal compressible gas dynamics with the important difference that the velocity is a vector of unit norm. It is structurally identical with that obtained from the Vicsek dynamics, with differences only in the values of the coefficients. An interesting open problem would be to investigate whether the coefficients of the two models coincide in the limit where the Cucker-Smale model converges to the Vicsek model. Future work will investigate different variants in the Cucker-Smale model, possibly involving phase transitions (which have been avoided in the present work by an appropriate expression of the alignment force).

Acknowledgments

PAS and PD acknowledge support by the Engineering and Physical Sciences Research Council (EPSRC) under grants no. EP/M006883/1 and EP/P013651/1. This work was accomplished during the visit of MB to Imperial College London with an ICL-CNRS Fellowship. JAC was partially supported by the EPSRC grant number EP/P031587/1. PD also acknowledges support by the Royal Society and the Wolfson Foundation through a Royal Society Wolfson Research Merit Award no. WM130048 and by the National Science Foundation (NSF) under grant no. RNMS11-07444 (KI-Net). PD is on leave from CNRS, Institut de Mathématiques de Toulouse, France. PAS and PD would like to thank the Isaac Newton Institute for Mathematical Sciences for support and hospitality during the programme ``Growth form and self-organisation'' when part of this work was carried over which was supported by: EPSRC grant numbers EP/K032208/1 and EP/R014604/1.

Data statement: No new data were collected in the course of this research.

Conflict of interest

The authors declare there is no conflict of interest.