Phase Transitions in a kinetic flocking model of Cucker-Smale type

We consider a collective behavior model in which individuals try to imitate each others' velocity and have a preferred speed. We show that a phase change phenomenon takes place as diffusion decreases, bringing the system from a"disordered"to an"ordered"state. This effect is related to recently noticed phenomena for the diffusive Vicsek model. We also carry out numerical simulations of the system and give further details on the phase transition.


Introduction
In many biological systems made of a large number of individuals such as cell populations [38], insect colonies [10] or vertebrate groups [29], agents are known to strongly interact with each other. Such social forces trigger the emergence of collective dynamics, where all individuals are moving coherently. Furthermore, it has been observed that such collective dynamics necessitate specific circumstances, such as e.g. a large density of individuals [10] while the same system exhibits disorganized dynamics when these circumstances are not met. Hence, when considering any sort of self-organized system, the question of how the model switches between disorganized and collective behavior becomes of tantamount importance.
Models of collective behavior abound in the literature, at the particle [1,2,16,22,25,36], kinetic [5,7,11,15,23] and hydrodynamic levels [6,21,28]. Among these models, the Vicsek model (VM) [8,21,36] is one of the simplest models exhibiting phase transitions between disordered to collective dynamics. In this model, particles moving with constant speed interact with their neighbors through local alignment and are subject to noise. Another extremely successful model is the Cucker-Smale model (CSM) [13,18,26,27]. In this model, particles can take all possible speeds and interact through local velocity consensus in a deterministic way, although noisy versions have also been considered [7,17].
In its original version, the CSM does not feature any self-propulsion. Particles move just because of inertia. In particular, if the initial total momentum of the particle system is equal to zero and in the absence of noise, particles become still in the large time limit. By constraining velocities to belong to a sphere, the VM exhibits a significantly different behavior, with phase transitions from disordered to collective states made possible [36]. This is why it is natural to equip the CSM with self-propulsion as proposed e.g. in [3,7] and to examine if such an augmented CSM exhibits similar phase transitions as the VM. Note that when the strength of the self-propulsion term is let to infinity, the noisy, self-propelled CSM converges to the VM as proven in [9]. Therefore, one may guess that this augmented CSM features phase transitions in at least some range of parameters. It is the goal of the present paper to prove this fact.
Phase transitions for the kinetic VM have been extensively studied in [19,20,24] in the spatially homogeneous case. Remember that the velocities are normalized, such that the velocity distribution is a probability distribution on the sphere. Supposing that the density is scaled to unity, then it is shown that there is a critical noise value D c of the noise intensity D such that for D > D c , only isotropic stationary distributions exist and are stable. By contrast, when D < D c , isotropic stationary distributions become unstable and a family of non-isotropic equilibria parametrized by a unit vector Ω on the sphere emerge and are stable. These equilibria are given by von Mises-Fisher distributions which are substitutes for the Gaussian distribution for probabilities defined on the sphere [37]. Note that [19,20,24] provide fully nonlinear (in)stability results which are obtained by taking advantage of the variational structure of the kinetic VM induced by a conveniently defined free energy functional. Now considering the noisy, self-propelled CSM, the first result towards the emergence of a phase transition is given in [3]. In this model, a subterfuge is used by considering a time-scale separation. For this purpose, the force acting on the particle is decomposed into the self-propulsion force on the one hand, and the force deduced from the combination of alignment with the neighbors and noise on the other hand. Each of these two forces is scaled with parameters ν and µ respectively. Now, considering again a spatially homogeneous situation for simplicity, [3] considers the limit µ → ∞ keeping ν fixed (i.e. when the alignment+noise contribution of the force is large compared to the self-propulsion force). In this limit, under a convenient time rescaling, the velocity distribution is shown to converge to a Maxwellian distribution whose mean velocity obeys an ordinary differential equation modelling the action of the self-propulsion force. Again, the number of equilibria that this ODE possesses depends on the noise intensity. Consistently with what was found for the VM, there again exists a threshold noise D c above which the only equilibrium mean-velocity is 0 and is stable, showing that no collective motion emerges. By contrast, for values of the noise intensity below D c , the zero mean-velocity equilibrium becomes unstable and a whole sphere of stable equilibria for the mean velocity emerges. The time-scale separation allows for an analytic determination of the equilibria of the interaction term (here alignment+noise), which turns out to be the Maxwellian. In particular, in the stationary states the mean velocity of this Maxwellian can be explicitly computed as a function of the noise. Without this scale separation hypothesis, the equilibria are given by a more complicated formula and their mean velocity is not explicitly known but is rather a solution of a nonlinear equation, making the analysis more difficult. The task tackled here is to provide a rigorous analysis of this equation.
In this work, we consider a noisy, self-propelled CSM. In this model, f is the distribution in both space x and velocity v at time t, and the model features a CSM term which aligns the velocity of individuals nearby in space, a term adding noise in the velocity, and a friction term which relaxes velocities back to norm one: Here K(x, y) is a suitably defined localization kernel and α and D are respectively the self-propulsion force and noise intensities. We have chosen scales such that the alignment force (modelled by the term (v − u f )f ) has intensity equal to 1. In this work, we focus on the spatially homogeneous case, where the model reduces to where and where f = f (t, v) is the velocity distribution at time t. Precisely, the goal of this work is to show that there is a phase transition between unpolarized and polarized motion as the noise intensity D is varied, for a specific range of the values of α. Therefore we achieve the goal of proving that the noisy self-propelled CSM behaves like the VM when the self-propulsion speed is large enough. Note that this problem has already been studied in dimension 1. Specifically, the existence of a transition between one and three stationary states has been proven in a one-dimensional setting by Tugaut in [32,33,34,35]. Here, we are interested in an arbitrary number of dimensions (practically 2 or 3) and we numerically demonstrate that there is indeed a phase transition by verifying that the stability of the isotropic equilibria changes as D crosses a threshold value. In addition, we analytically prove that for large noise D there is only one isotropic stationary solution, while for small D there is an additional infinite family of stationary states parametrized by a unit vector on the sphere, below referred to as the polarized equilibria. Therefore, although no analytic formula for the mean velocity of the polarized equilibria of the CSM exist, we know that at least this family of equilibria forms a manifold of the same dimension as the polarized equilibria of the VM. We remind that, in the VM case, the polarized equilibria were given by Von Mises-Fisher distributions about an arbitrary mean orientation Ω belonging to the unit sphere. Here, these polarized equilibria are still parametrized by a vector Ω of the unit sphere, but the precise expression of the mean velocity is not analytically known.
The paper is organized as follows. In Section 2, we introduce the homogeneous problem and discuss its steady states. We then focus on multiple dimensions and postulate that there are two regions of the parameter space, each with a different number of possible stationary solutions. We proceed to show that for small D, there is a manifold of equilibria parametrized by a vector on the unit sphere, while for large D, there can be only one. In section 3, we numerically validate the results of Section 2 and explore the influence of the parameter α and the stability of stationary states. We conclude and discuss future work in section 4. Some technical details are provided in Appendix A.

Phase Transition: Stationary States for the Homogeneous Case
In this section we focus on probability density solutions to the spatially homogeneous kinetic model associated to (1), that reads with t ≥ 0 and v ∈ R N . Here, the mean velocity u f is given by The second term on the right-hand side of (3) accounts for the tendency to align with the local velocity field, while the last term adds noise into the velocity component. The first term enforces the tendency to travel with unit speed. The kinetic equation satisfies the conservation of mass for all t ≥ 0. This equation lies in the general class of PDEs having a gradient flow structure, see [14], by writing the equation as Here, particles are thought to move under the effect of a confining potential given by an interaction potential of the form W (v) = |v| 2 2 , and with linear diffusion. The equation (3) is then seen as a continuity equation with a velocity field of the form −∇ v ξ, and thus there is a natural entropy for this equation given by the free energy of the system since ξ = δF δf . The second expression follows by expanding the square in the interaction potential. Actually, the dissipation of the free energy F[f ] in (5) along solutions is given by We will look for stationary solutions f (v) > 0 to (3) for u f =ū. Taking into account the dissipation property, they should satisfy ∇ v ξ = 0, or equivalently ξ = constant. Thus, stationary solutions are of the form with Z the normalization factor such that f has unit mass. Therefore, the set of stationary solutions of (3) can be parametrized by the set of mean velocitiesū ∈ R N such that Notice that by (7), H(ū, D) = 0 is equivalent toū = u fū given by (4). Let us point out thatū = 0 is always a solution corresponding to radially symmetric stationary states.  it is apparent that the sign of ∂H ∂u (0) shifts from negative to positive as D increases.
By choosing the axis, we may assume without loss of generality thatū points in the direction of the first axis or first vector e 1 of the canonical basis, and then let us denote the magnitude ofū by u ≥ 0. The full set of stationary solutions is obtained by composing fū with any rotation in R N , and thus yields an (N − 1)dimensional family of stationary solutions for eachū = ue 1 satisfying (7). Noticing that all components of H except for the first one vanish due to fū(v) being odd in v 2 , . . . , v N , we can restrict our attention to the first component of H. For the sake of simplicity we will denote by f the probability density given by (6) associated to the vectorū = ue 1 , and the real valued function whose roots have to be analyzed is the first component of H, given by with In figure 1, we plot H(u, D) in one dimension as a function of u for varying values of D. It is clear from the figure that for small values of D, H(u, D) has three roots, the zero root and two roots with identical speed; while for large values of D the only root is u = 0, and this can be deduced from the sign of ∂H ∂u (0, D). We will show analytically that this is the case in the next subsections and it will be further studied numerically in section 3.
Our goal is to show that, for given α > 0, as we vary the noise strength D, there is a region of the parameter space with only one possible solution, namely u = 0, and a region with at least two roots, u = 0 and u = u α,D > 0. In fact, we expect to have the unique homogeneous stationary state for large noise corresponding to a disordered state while for small noise we expect to have a nontrivial biased solution. The objective of the next two subsections is to show these facts for small and large noise.
The main theorem of this section can be summarized as follows:  Let us notice that the previous theorem does prove the appearance of a phase transition, although it does not give information about the critical value where it occurs. We will numerically show in section 3 that this phase transition is continuous and it happens at a sharp value of D as in [19] for the continuous Vicsek model. We will numerically demonstrate in subsection 3.3 the time asymptotic stability of the spatially homogeneous solution for large noise and the non homogeneous solution for small noise in one dimension using a Monte Carlo-like particle method.

D → ∞ Case: Unique Disordered State
Noticing that we may integrate by parts and obtaiñ Proof. Computing ∂H ∂u , we get The term 1 − |v| 2 obviously determines the sign of ∂H ∂u . We will compensate the positivity of this term on the unit ball with a piece of the integral outside. First, let us estimate the integrand inside the unit ball, We will also use that the other terms are close to 1 in any bounded region for D large enough. More precisely, for D ≥ D( ) large enough and |v| < 4.
Here, the notation η refers to the vector (3, 0, . . . , 0) and B 1 (η) denotes the Euclidean ball of radius 1 centered atη. We separate the integrand into three pieces corresponding to the sets A, B, and C. Since the integrand is negative in C, then we can estimate where I and II are the integrals restricted to A and B respectively. Taking into account (10) and (11), we control the first term as Similarly, the second term is bounded by (11) and since v 1 > 2 and |v| 2 ≤ 4 in B. In order to show that ∂H ∂u ≤ 0, we need only show that I + II ≤ 0: Here, we use that |A| = |B|.
The preceding claim proves that H(u, D) can have at most one nonnegative root, and since H(0, D) = 0, there can be no other positive root. Hence, in the case of D → ∞, there can only be the radially symmetric stationary solution associated to u = 0.

Laplace's Method
Laplace's method gives the asymptotics of integrals of the form for given functions f (x) and g(x) as D → 0. A usual statement of Laplace's Method that is commonly found in the asymptotic analysis literature (see [4,31] for instance) reads as follows: 3. Assume f is a smooth function that has a unique global minimum at x 0 , and that there exist We state and briefly prove a modified version of it including higher order terms which is well adapted to our arguments, avoiding general statements which become cumbersome in a multidimensional setting. We essentially follow the strategy in [4,Chapter 6]; see also [31] for multidimensional statements. For a multi-index β = (β 1 , . . . , β N ) we denote Of course, the constant M β is 0 whenever one of the β i is odd, and M β may also be expressed in terms of the Gamma function. The following calculation is at the basis of Laplace's method: ..,N be a multi-index and Let Q : R N → R be the quadratic function given by Proof. Carrying out the change of variables directly yields the given expression.
From the previous result one directly obtains the asymptotics of integrals with any polynomial g(x) instead of x β . By a linear change of coordinates it is then also simple to extend the result to include any positive definite quadratic form in the exponential, but we will not need that for our purposes. Extending the result to more general functions in the exponential requires a more careful argument that we give next. We begin by noticing that the only region asymptotically contributing to the integral is concentrated around x = 0: The constant implicit in the O notation depends continuously on the coefficients λ i .
Proof. Due to Lemma 2.4 we just need to estimate the difference to the integral over all of R N : where we have assumed D ≤ 1 since the statement concerns only the asymptotics as D → 0. This gives an explicit bound of the remainder term.
We now state the main result on Laplace's method that we use in this paper: Theorem 2.6. Let P : R N → R be a polynomial function given by where a 0 ∈ R, λ i > 0 for all i = 1, . . . , N and R contains only terms of degree three or higher. Assume that for some µ > 0, for all D ≤ 1. For n ∈ N we have the expansion where the numbers K i and the constant implicit in the O notation depend continuously only on µ and the coefficients of P and g, and M 0 := R N e −|x| 2 . In addition, K i depends only on the coefficients of g of degree at most 2i, and is equal to 0 if g has no such terms.
Remark 2.7. Condition (12) implies in particular that the global minimum of P is attained at x = 0. Moreover, it requires that the minimum be strict in a specific sense.
Proof. The constant term a 0 obviously gives the exponential factor e −a0/D , so me may assume that a 0 = 0. For D ≤ 1 we choose δ := D 1/3 , and break the integration into the region inside the ball B δ (0) and its complement. We first show that the integration outside this ball is very small as D → 0: using (12) and the inequality min(|x| 2 , 1) ≥ δ 2 for |x| ≥ δ = D 1/3 and D ≤ 1, then an argument very close to that in Lemma 2.5 shows that valid for all D ≤ 1. For the integral inside the ball of radius δ, denote where Q is the sum of all second-order terms of P and R is the sum of the remaining terms (of order greater than or equal to 3). We can expand e − R(x) D to obtain, for |x| ≤ δ, where the implicit constant depends only on a bound of R by |x| 3 (which can be chosen as a continuous function of its coefficients). We then have, using Lemma 2.5 in order to estimate the remainder term, We can use again Lemma 2.5 to estimate each term, since each of them is a quadratic exponential times a polynomial. Now, let us remember that M β = 0 in Lemma 2.5 whenever there is an index k such that β k is odd. This suggests rewriting the polynomial in the integrand as where the r i e (x) are even polynomials (all their monomials are even) and r i o (x) are odd polynomials. Therefore, Lemma 2.5 gives that Now, we can identify the full expansion in powers of D. The only term that contains terms of the lowest order is the first one: The rest of the terms have order strictly higher than this, and it is seen from Lemma 2.4 that their coefficients are continuous functions of the coefficients of P and g. We now observe that only integer powers of D will appear in the expansion due to the evenness of the polynomials in the remainder and Lemma 2.4, note that |β| is even for all the monomials in the expansion. One also sees that if g contains no terms of degree lower than or equal to i then every term in the expansion is of order at least D N/2 D i+1 , and hence the coefficient K i is equal to 0.

D → 0 Case: Multiple Solutions
We will now show the existence of a curve of nonsymmetric stationary states emanating from the stationary states for D = 0 (which are the measures δ u , for any u ∈ R N with |u| = 1). Since stationary states are determined by the roots of equation (8), we are interested in the behavior of H(u, D) as D → 0. The parameter 1/D appears inside the exponential, and the asymptotics of integrals of this form is given by Laplace's method, particularly by the statements given in the previous section. Intuitively, we expect the stationary distribution to approach a Dirac delta at the minimum of the polynomial P u (v) from eq. (9) as D → 0 (this will be rigorously justified by Lemma 2.6 as we will see below). Let us assume for the moment that for u > 0 there is a unique minimum of P u (v) that is achieved at v = v * (u) with v * k (u) = 0, k > 1 (we will also come back to this point next to check that this assumption is met). Therefore, if we want to find u > 0 such that H(u, 0) = 0, we can compute H(u, 0) formally at this point by substituting f by a Dirac delta at v * (u) in (8) In this section, we will rigorously justify this, and follow a perturbative argument to show that there is curve of solutions of H(u, D) = 0 that converges to (1, 0) as D → 0.
Global minima of P u . Let us first find the minima of the polynomial P u (v). Its gradient is The critical points of P u (v) are thus characterized by ∇ v P u (v)·e k = 0 for all k = 1, . . . , N (where {e 1 , . . . , e N } is the usual base of R N ). That is, which means that either v k = 0 for all k = 1, or |v| 2 = α−1 α = 1 − 1 α . For u > 0 we cannot find critical points for which |v| 2 = 1 − 1 α due to (15). Hence, in the case u > 0, critical points must satisfy v k = 0 for all k = 1 and for all α > 0. Therefore, in the case u > 0 all the critical points v = (v 1 , . . . , v N ) satisfy The case u = 1 can be explicitly solved since . It is simple to check that v 1 = 1 is the unique global minimum of P 1 (v). In general, for any u > 0 it can be seen that (16) has one positive root v 1 = v * 1 (u), and the remaining roots are either complex or negative, depending on the values of α and u (this can be checked by differentiating again in v 1 ). Since u > 0, it is easy to check that P u (−v 1 , v 2 , . . . , v N ) > P u (v 1 , v 2 , . . . , v N ) for all v 1 > 0, then the global minimum of P u (v) must be attained only at v = v * (u) = (v * 1 (u), 0, . . . , 0). We have then proved the following: Lemma 2.8 (Global minimum of P u ). For u > 0 and α > 0 the polynomial P u attains its global minimum only at is a continuous function of u > 0, is positive for u > 0, and v * 1 (1) = 1. Of course, v * also depends on α, but we omit this dependence in the notation since α is fixed in all arguments below. Notice that the continuity in u of v * (u) is a consequence of the continuity in u of the unique positive root of (16).
A decomposition of P u . For u > 0 we can write our polynomial P u (v) as with a 0 (u) = P u (v * (u)) and where Λ is the Hessian matrix of P u (v) at the global minimum v * (u). It can be calculated in terms of v * 1 as where the sum is over multiindices of the given order. In particular, at u = 1 we get Auxiliary functions. In order to make use of Theorem 2.6 here and in later sections, let us define the functions with k = 0, 1, 2. Applying Theorem 2.6 to F k (u, D) we conclude that and for u > 0 as D → 0, where c k (u), k = 0, 1, 2, are continuous functions of u, and the constants implicit in the O notation are uniform in a neighborhood of u = 1 (one can check that all conditions in Theorem 2.6 hold uniformly in a neighborhood of u = 1). The explicit expression of the first term in the expansion in Theorem 2.6 gives and thus by continuity we have c 0 (u) = 0 in a neighborhood of u = 1. For reference below, we take 0 > 0 such that (22)-(24) hold for |u − 1| < 0 . Analogously, we can use the expansion in Theorem 2.6 to obtain the first order term of the functions F 1 (u, D) and F 2 (u, D) at u = 1 to get and

(See Appendix A for the explicit calculations leading to this.)
Continuity of H as D → 0. The function H(u, D) is smooth with respect to u and D for all u > 0 and D > 0, as can be seen by standard arguments. Let us show that H(u, D) has a limit as D → 0 (which will enable us to define it by continuity at D = 0). It is easy to verify the following formulas that relate H to F 0 and F 1 : and We deduce from (27) taking into account (22) and (23) that for |u − 1| < 0 . As a consequence, by defining the function H(u, D) is continuous in (u − 0 , u + 0 ) × [0, +∞).

Differentiability in u.
It is straightforward to check that We proceed as before: applying (22)- (24) in (29) we obtain for |u − 1| < ε 0 . This shows the function H(u, D) (extended to D = 0 as in (28)) is differentiable with respect to u in a neighborhood of (1, 0). It is simple to check that This comes again from the explicit computation of the first term in the expansion of Theorem 2.6 applied to F 2 (u, D) which is given by the second moment centered at v * 1 (1) = 1 of Q 1 (v).

Differentiability in D.
In a similar way we can write Inserting (27) into (31), this is equivalently written as and with the decomposition (17) we get the expression Applying Theorem 2.6 to the two integrals in (33) to obtain as D → 0 where k 1 (u) and k 2 (u) are continuous functions for |u − 1| < ε 0 . In fact, using the expansion in Theorem 2.6 we obtain and whose expressions are derived in Appendix A. These expressions, together with the formulas (22) and (23), for |u − 1| < ε 0 . This shows that the function H(u, D) is differentiable from the right with respect to D in a neighborhood of (1, 0). Moreover, we find that for all α > 0, according to the values in Appendix A.
Existence of a curve of solutions close to (u, D) = (1, 0). Summarizing, we have proved that the function H(u, D) is a differentiable (C 1 ) function in a neighborhood of (1, 0), such that H(1, 0) = 0. Equation (30) implies that ∂H ∂u (1, 0) = 0 which allows us to apply the implicit function theorem, implying that there exists a curve u = u * (D) defined for D small enough such that H(u * (D), D) = 0. This shows the existence of a curve of non-symmetric stationary states emanating from the Dirac delta at v = e 1 for D = 0.

Numerical Results
In this section we numerically validate the results of the previous section, finding the bifurcation curves and numerically showing the stability of the stationary solutions. In addition to demonstrating the analytical results related to the parameter D, we explore the effect of the value of the parameter α on the critical noise threshold and the effect of both parameters on the stationary velocity profile. To emphasize the dependence of H(u, D) on the value of the parameter α we will use the notation H α (u, D) throughout this section. By examining the roots of H α (u, D), we numerically validate in both one and two dimensions the fact that there is more than one stationary state for small magnitudes of noise, and that there is only one stationary solution for large noise. Using both H α (u, D) and ∂Hα ∂u (0, D), we show where in α-D parameter space the transition from more than one to one stationary solutions occurs. We then numerically explore the α → ∞ limit.
To examine further properties of the system we use a Monte Carlo-like particle method to approximate the steady states and the transient behavior of the system, employing the Euler-Maruyama method to numerically solve the SDEs, see [30] for instance. With this framework, we are able to numerically validate the stability of the nonzero stationary solution when it exists, giving evidence that this bifurcation is indeed a phase transition. Using a large ensemble of independent runs, we also track the temporal evolution of both the average velocity and the free energy F defined in Section 2.
In order to efficiently compute the stationary states and the bifurcation diagram in two dimensions, we use radial coordinates. In fact, we can rewrite the function H α (u, D) in radial coordinates in any dimension as follows. Defining E D (r) = exp α D ( r 2 2 − r 4 4 ) − r 2 2D ) we can reexpress H α (u, D) in coordinates in terms of the angle with respect to the first axis as: Let us use the notation We further reduce to: Formula (36) in the two-dimensional case leads to an expression in terms of modified Bessel functions of the first kind and in the three-dimensional case they can be obtained explicitly in terms of hyperbolic sine and cosine functions. This is exploited in the two-dimensional numerics of subsections 3.1 and 3.2. Note that these expressions do not simplify the analytical discussion in Section 2 of the behavior of H α (u, D) for large and small noise D.

Bifurcations
We are able to numerically observe the bifurcations predicted by the analysis in the previous section. In Figures 2 and 3, we show in one and two dimensions, respectively, the root u(D) of H α (u, D) plotted as a function of D in solid lines. The curves were determined using a root-finding function on H α (u, D) for each fixed value of D. We also plot ∂Hα ∂u (u(D), D) as dotted lines for varying values of α. Here, the derivative is computed by substituting the root u(D) into the formula (29) of ∂Hα ∂u (u, D), given in subsection 2.3. Let us define the critical value of the noise D c as the noise at which u(D) attains zero for the first time. From Figures 2 and 3, it is clear that in both one dimension and two dimensions, ∂Hα ∂u (u(D), D) is equal to zero at D c . This tells us that the bifurcation branch intersects the zero branch perpendicularly, since This is unsurprising, given the form of H α (u, D) shown in figure 1: the slope of H α evaluated at the nonzero stationary solution should be negative, becoming increasingly shallow until it becomes zero. We also note that the bifurcation diagram is decreasing for D small as indicated by the formulas found in previous section, since ∂u

The role of α
According to the formulas (36), in two dimensions, I 2 0 (0) = π, I 2 1 (0) = 0, and I 2 2 (0) = π 2 , so F 0 (0) = π u for alpha=2 u for alpha=4 u for alpha=6 u for alpha=8 u for alpha=10 u for alpha=12 u for alpha=14 dH/du(u) for alpha=2 dH/du(u) for alpha=4 dH/du(u) for alpha=6 dH/du(u) for alpha=8 dH/du(u) for alpha=10 dH/du(u) for alpha=12 dH/du(u) for alpha=14 u for alpha=2 u for alpha=4 u for alpha=6 u for alpha=8 u for alpha=10 u for alpha=12 u for alpha=14 dH/du(u) for alpha=2 dH/du(u) for alpha=4 dH/du(u) for alpha=6 dH/du(u) for alpha=8 dH/du(u) for alpha=10 dH/du(u) for alpha=12 dH/du(u) for alpha=14  We plot points in parameter space where there exists a nonzero stationary state in red, and points where we find only the zero stationary state in blue. The line of demarcation between the two regions is created using a continuation method to find the root of ∂Hα ∂u (0, D). Taken together, these roots, which determine for which D the slope of H α at zero changes from positive to negative, define the line.
Thus, in two dimensions, This calculation highlights the dependence of the bifurcation curve on the parameter α. We demonstrate this numerically in both one and two dimensions in figures 2 and 3, respectively, where we can observe how the bifurcation curves change as we vary α.
In figure 4, we numerically determine where in α-D parameter space this bifurcation occurs in one dimension. This bifurcation diagram was found both analytically and numerically by Tugaut in [33,Subsection 4.1]. Here, we sample the parameter space, with α along the vertical axis and D along the horizontal one, and plot a blue dot when the point has more than one stationary solution and a red dot where it has only one. The black line was drawn using the continuation method to find the root of ∂Hα ∂u (0, D). It is clear from the figure that for D sufficiently large, we could consider α to be the bifurcation parameter, identifying a critical value α c at which a phase transition occurs. Unlike the case of increasing D, this bifurcation is from one stationary solution to more than one as α increases, and making this phase transition explicit remains to be explored. We finally mention that we observe that the critical value D c has a limit as α → 0, this fact was already studied analytically in [33,Subsection 4.1] showing that its limit value is 1/3 matching with our simulations.
It is interesting to note that the changes in the bifurcation curve lessen as α increases, see figure 5 in the two-dimensional case. This indicates that the bifurcation curves are approaching a limiting function as α → ∞. Letting α → ∞ means that the cruise speed term dominates the behavior of the system; intuitively speaking, as α → ∞, we recover a norm constraint in the velocity for particles. This intuition was in part rigorously proved by Bostan and Carrillo in [9], where they show that the kinetic equation (3) limits to the continuum Vicsek model in [19]. Here, we numerically conclude that a limiting phase transition curve does indeed exist in two dimensions and it qualitatively matches the one obtained in [24,19]. In fact, the critical noise value D c is converging towards the critical noise value 1/2 for the Vicsek model obtained in [24,Proposition 3.3].

Stability and Phase Transition
In order to numerically explore the stability of the stationary solutions, we approximate the solutions to (3) by a Monte Carlo-like method using the Euler-Maruyama numerical scheme. We use 10000 particles per run with a timestep of 0.01 and evaluate the average velocity at time 6000. This is enough for stabilization in time of the solutions except for noise values close to the critical noise parameter D c , which are well known to take longer to converge. In the particle simulations, we initialize the particles with velocities sampled from the Gaussian N (0.5, 1) in order to investigate which of the stationary solutions is stable. In figure 6, we plot this final average velocity over an ensemble of ten runs on top of three of the bifurcation curves studied in Figure 2. From the figure, it is clear that the particle simulations, initialized away from either stationary state, approach the nonzero stationary solution while such a solution exists. As expected, this indicates that, when it exists, the nonzero stationary solution is stable. This demonstrates that the zero stationary solution is unstable before the critical value of the noise and stable afterward, indicating the bifurcation we observe is indeed a phase transition.

Stationary Solutions
One validation of the efficacy of the numerical method with the particles is whether we are able to recover the stationary solutions for different values of D. In Figure 7, the dots show the final histograms at time 500 of an ensemble of 100 runs in one dimension with α = 2 and varying D. In solid lines, we plot the solution f u(D) given by equation (6), taking u(D) to be the values computed in subsection 3.1. We observe an impressive match between the histograms and the analytical distribution, further validating the Monte Carlo-like approach to the solutions.
With the formula for the steady-state and the stationary average velocities that we found numerically, we are able to consider also the shape of the stationary distributions in velocity space as α changes. The results shown in Figure 8 are for the one-dimensional case, though this could easily be done for the two-dimensional case, as well. We observe that the stationary distribution for D beyond the critical threshold D c can take forms which are not Gaussian around 0. In fact, they can be double-peaked, with the peaks at v = −1 and v = 1. We can see from Figure 8 that the preference for velocities of norm 1 is not strong enough to create a double-peaked distribution for the case of small α, but as α grows, there is a definitive preference for velocities with norm 1. Beyond the critical noise threshold, double peaks at velocities with norm one are clearly apparent even in cases in which u(D) = 0.

Free energy and average velocity in 1D
The particle simulations give us insight beyond indicating the stability of the stationary states. Using an ensemble of 100 runs, each with 10000 particles and a time step of 0.001, we are able to construct a histogram approximating the velocity profile f in one dimension. We used this to calculate the average velocity over time and to calculate the free energy given in equation (5) and plot it over the course of the simulation. In Figure 9, we show the average velocity over time of the particles in the ensemble of runs. We can see that Here, the average velocity over an ensemble of ten runs of the microscopic model is plotted (asterisks) over bifurcation curves (solid line) similar to those in figure 2. The particles were initialized with velocity sampled from N (0.5, 1). It is clear from this plot that the average velocity of the microscopic runs agrees with the nonzero stationary solution, indicating that this is the stable stationary state as long as it exists.
the average velocity initially dips for all of the values of D, presumably as the particles align, and then more slowly approaches the velocity that we know from figure 6 to be the stationary value of the average velocity.
In Figure 10, we plot the evolution of the free energy from the histograms of the particles. We again see an initial swift decline in the free energy, followed by a gradual decay of the free energy for small values of D, as is expected. For larger values of D, after this initial steep decay of the free energy we notice a gradual increase. Note that this particular method is not designed to preserve the decay of the free energy. In fact, we observe that our method initially undershoots the asymptotic value of the free energy for large noise since the limiting free energy value of the stationary state is accurately computed in view of the results in Figure  7. Free-energy decreasing deterministic methods [12] could be used to investigate this issue further.

Conclusion
In this work we have studied the stationary solutions of a noisy self-propelled Cucker-Smale model and proven than this model behaves similarly as the Vicsek model, i.e. it exhibits a phase transition when the noise intensity increases. For large noise intensity, we show the existence of a single equilibrium with zero meanvelocity, showing that no self-organized motion is possible. For small noise intensity, we show the existence of a family of polarized equilibria parametrized by a unit vector on the sphere. There are still a number of issues left. For instance, so far, we cannot rule out the existence of other families of equilibria for small noise, as our approach relies on a local approach using the implicit function theorem. Also, several phase transition points could exist with the emergence of other branches of equilibria. Note that these circumstances are not occurring in one dimension thanks to the results of Tugaut [32,33,34,35]. Another direction is to take advantage of the polarized equilibria to develop a hydrodynamic model in a similar spirit as [21] for the Vicsek model. Here again, we expect that the absence of analytical formulas for the mean velocity of the equilibria will generate additional difficulties. Numerically, it could be interesting, but challenging to compare the dynamics resulting from the Vicsek model and the self-propelled Cucker-Smale model with high precision.
The assumption of constant speed in the Vicsek model has often been disputed and this comparison would help determine what consequences this assumption has on the results. A similar comparison could be made at the level of the hydrodynamic models and would be equally useful.   In all of the figures, we plot the stationary distributions for several D values; from left to right, we fix α to be 0.001, 1, 2, 4, and 6, respectively. The stationary distributions are computed by substituting the average velocity which we found numerically into the formula for the stationary distribution (6).