STABILITY OF CONSTANT STATES AND QUALITATIVE BEHAVIOR OF SOLUTIONS TO A ONE DIMENSIONAL HYPERBOLIC MODEL OF CHEMOTAXIS

. We consider a general model of chemotaxis with ﬁnite speed of propagation in one space dimension. For this model we establish a general result of stability of some constant states both for the Cauchy problem on the whole real line and for the Neumann problem on a bounded interval. These results are obtained using the linearized operators and the accurate analysis of their nonlinear perturbations. Numerical schemes are proposed to approximate these equations, and the expected qualitative behavior for large times is compared to several numerical tests.

1. Introduction. The term chemotaxis indicates the motion of a population driven by the presence of an external (chemical) stimulus, specifically in response to gradients of such substance, usually called chemoattractant. The basic feature of such phenomena is the presence of concentration effects (chemotactic collapse), possibly leading to non-uniform pattern formation. A classical introduction to these models is contained in [27]; a thorough and more recent overview of the subject can be found in [34] and in [37]- [38]. From a mathematical point of view, the basic aim is to find thresholds for structure formation and to describe such structures. Predictions given by mathematical models should of course be compared with reality, always taking into account that simplified chemotaxis models can describe only preliminary phases in aggregation phenomena. When a sufficiently large and complex structure is formed, new mechanisms appear and the model becomes no more adequate.
The basic unknowns in PDE models for chemotaxis models are the density of individuals of the population and the concentration of the chemoattractant. The first and most celebrated model of these phenomena is the Patlak-Keller-Segel model (see [33,25] or system (7) below). In this case, the basic assumption is that dynamic of individuals is described by a parabolic equation coupled with an additional equation for the chemoattractant, chosen to be elliptic or parabolic, depending on the different regimes to be described and on the taste of the authors. Chemotaxis appears as a cross-diffusion term in the equation for the population. A large amount of articles and studies analyzed the Patlak-Keller-Segel system. For a complete review on the mathematical results see [17,22,34].
Existence of stationary solutions for the parabolic Keller-Segel model (together with their stability and bifurcations) has been studied in [7] for the simplest model and later generalized in [35] for more general chemotactic sensitive functions. Notice that in two and three space dimensions it is well-known that solutions can blow-up in finite time, see again [22] and the more recent results in [5] and references therein. However in one space dimension the global existence of solutions for general initial data has been shown by Osaki and Yagi [29], by using energy methods. Moreover, they have shown the existence of a global attractor under certain assumptions on the chemotactic sensitivities and the initial data, see also [19] for some refined results. Different models, taking in account realistic effect preventing overcrowding, can provide global existence of solutions also in several space dimensions (see [18,32]).
One of the problems of the diffusion models is that they imply an unrealistic infinite speed of propagation of cells. This approximation can be accepted in some large time regimes, but it is usually just too rough to take into account the fine structure of the cell density for short times. To circumvent such flaw of the classic Keller-Segel model, models based on hyperbolic equations have been considered starting from the papers [36,12]. In such chemotaxis models, the population is divided in compartments depending on the velocity of propagation of individuals, giving raise to kinetic-type equations, either with continuous or discrete velocities.
The basic example of a discrete kinetic model for chemotaxis is where γ, D > 0 and f, g, h are smooth functions satisfying suitable assumptions to be introduced later on. System (1) has been first considered in [12] and generalized in [21]. The derivation of macroscopic models (i.e. Patlak-Keller-Segel system) from appropriate rescaled hyperbolic/kinetic equations has been considered in many works, see for instance [16,21,31,6,9,13,14,24]), showing that, heuristically, hyperbolic models can be interpreted as a description of chemotaxis phenomena at a mesoscopic scale. While the literature on the Patlak-Keller-Segel model refers mainly to multidimensional problems, many papers on hyperbolic models deal with the one-dimensional case, see [21,13,20,23]. Multidimensional hyperbolic systems for chemotaxis are considered in [8,14,10], mostly from a numerical point of view. More elaborated models have been introduced in [11,26] in the framework of vasculogenesis modeling, with the choice of an Euler-type term v 2 /u + P (u) in place of γ 2 u, arising from the assumption that closely packed cells show limited compressibility and generate a pressure-like term in the equations, see also [39,1] for more details on this type of derivation. The analytical study in several space dimensions for these hyperbolic models is far to be complete, and will be the object of further studies.
In this paper, we want to perform an analytical study of problem (1), under very general assumptions on the coefficients, with aim of improving on previous results in [21,20,23], by determining sharp stability results of constant states, and interpreting instability as the sign of a non-homogeneous pattern formation. In particular, we are able to remove the main assumptions in those papers, namely that the functions h and g and the constant γ in (1) satisfy the inequality γh(φ, ψ) ≥ |g(φ, ψ)|, for all φ, ψ ∈ R 2 . ( Although this assumption can be motivated for simple models by some microscopic probabilistic arguments, at the macroscopic level it is just an artificial restriction with respect to more realistic models, see for instance [11,26,8], and also with respect to the corresponding diffusion limit. To remove this assumption, we proceed as follows. After a basic linear analysis of the operator (1), performed in Section 2, in Section 3 we give the proof of the stability of small perturbation of the zero state in the case of the Cauchy problem on the whole real line. First, using the accurate description of the Green function for the damped wave equation given in [3], we obtain sharp decay estimates for the linearized operator. Therefore, by Duhamel's principle, we are able to prove the stability result with (power) time decay estimates.
In Section 4, we consider the same system (1) on a bounded interval with Neumann boundary conditions. In this case, we are able to prove the exponential time decay of the linearized operator, which, in turn, yields a stronger stability result, namely the global existence and decay of solutions corresponding to small perturbations of all stationary constant states which are set in the region of stability given by the linear analysis, see Section 2 below.
The analysis is complemented with numerical experiments with the final goal of understanding qualitative properties of the system, even in regimes where a rigorous analysis is still not available. First, in Section 5, we give some details about the whole set of stationary solutions of system (1) and what is known of their stability. Some numerical bifurcation diagrams of steady states are displayed. In Section 6, we introduce two explicit-implicit numerical schemes to approximate system (1), for the bounded interval, with Neumann boundary conditions. The second one is specially designed to preserve the asymptotic profiles, much in the spirit of [2]. Thanks to these approximations we are able to present in Section 7, accurate numerical results about the asymptotic behavior of solutions of (1). This behavior strongly depends on the mass of initial data u 0 and on the eventual symmetries of initial functions u 0 and φ 0 . We explore the stability of the first few branches of bifurcation of the stationary solutions, also making an accurate comparison with the behavior of the corresponding diffusion model.

2.
The one dimensional discrete kinetic model. Many species alternate periods of straight movements in a fixed direction with others during which a new direction is chosen, as a consequence of presence of some attracting or repelling chemical. A simple one-dimensional model describing such combination of behaviors is given by the following hyperbolic-parabolic system which is strongly inspired by the Goldstein-Kac model for correlated random walks. Such kind of models were originally considered in [36], and later reconsidered in [12]. More recently, some generalizations of these models have been considered in [21,20,15,23].
Functions u ± denote the densities of the right/left moving part of the total population and φ is the external chemotactic stimulus biasing the movement of the population itself. Here, parameters γ, τ, D are assumed to be strictly positive constants. They represent, respectively, characteristic speed of propagation of u ± , time-scale for the dynamics and diffusion coefficient for the chemoattractant. After rescaling of D and f , we can set τ = 1. The terms µ ± are called turning rates and they control the probability of transition from u + to u − and vice versa, i.e. the change of direction in the movement of a single individual. As a consequence, they are assumed to be functions of φ and ∂ x φ, i.e. µ ± = µ ± (φ, ∂ x φ). The function f describes the production/degradation of the chemoattractant. Details on the specific form of µ ± and f will be given later on.
Introducing the two moments where Assuming additionally h ≡ 1 and formally disregarding the term ∂ t v in the second equation of (4), the system reduces to the Keller-Segel parabolic system To deal with both the original Keller-Segel system (5) and hyperbolic system (4), let us rewrite this system as with δ ≥ 0.

Remark 1.
Assuming again h(φ, ψ) ≡ 1, the first two equations in (6) can be rewritten as a single second order nonlinear wave equation for u The Keller-Segel system (5) is again recovered for δ = 0.
Next we want to find some reasonable assumptions on functions g and h, recalling that µ ± describes the rate of transition from population u ± to u ∓ and that such transition is driven by the presence of a gradient of the chemoattractant φ.
Following [12], we can choose µ ± to be affine functions of ∂ x φ; for instance we can take where χ = χ(φ) > 0 is a smooth function. This choice corresponds to which in the parabolic limit δ = 0 of equation (5), just coincides with the classical Keller-Segel system. Coming back to our general framework, let us assume for the moment only These assumptions reflect the fact that, in the absence of gradients of chemoattractant, the rate of transition from state u ± to u ∓ is the same and that some transition is always occurring.

Remark 2.
It is meaningful to assume that functions g and h are such that the system (6) is invariant with respect to the transformation (x, t, u, v, φ) → (−x, t, u, −v, φ). Such invariance has to be interpreted as a "one-dimensional rotational invariance". This amounts in requiring In terms of the turning functions µ ± , (10) corresponds to the relation Hence, if µ ± are assumed to be non-negative, then condition (10) is equivalent to the symmetry condition µ ± (φ, ψ) = µ ∓ (φ, −ψ), considered in [21,20,23].
If assumption (9) holds, system (6) has constant solutions of the form We are interested in determining critical parameters for stability/instability of such states. In particular, transition from stability to instability is expected to be an indicator of the presence of non constant solutions, bifurcating from constant solutions. Linearizing system (6) at the constant steady state (U, 0, Φ), we obtain where h(Φ, 0) = β > 0 by (9), By rescaling variables, we could reduce the number of parameters. However, we prefer to keep all of them, in order to determine limiting behavior when one or more of such parameters are assumed to be negligible.
Remark 3. When considering system of the general form (1) under the assumption ∂ v F (U, 0) = 0 (satisfied by the Euler type model), the linearized equations have the form given in (11) with γ 2 = P (U ). Hence the analysis we perform in what follows for the linearized stability holds also in that case. Linear stability analysis for the Euler type model has been performed in [26] assuming the presence of a viscosity term in the equation for the momentum v.
Remark 4. The linearized system (11) has the invariance described in Remark 2, regardless of the specific choice of g and h. It possesses also additional invariances, due to the linear structure of the problem. As an example, the system remains unchanged when applying the transformation ( Such invariance has no specific counterpart at the level of the nonlinear problem (it would correspond to a transformation of the form ( The dispersion relation of (11) is given by Next, we want to determine stability/instability conditions on the state (U, 0, Φ). Such conditions can be derived by analyzing, for any k the values λ such p(λ, k) = 0. First of all, we compose the Routh-Hurwitz array for δ > 0: By the Routh-Hurwitz criterion (see for instance [27], vol. I, Appendix B.1), all of the eigenvalues λ = λ(k) have negative real parts if and only if the elements of the first column in the array (13) are positive. By (9), β > 0. Assuming b > 0, the stability condition is satisfied if a 0 > 0, i.e. if, for any k under consideration, there holds U < γ 2 (b + Dk 2 ) a χ .

A ONE DIMENSIONAL HYPERBOLIC MODEL OF CHEMOTAXIS 45
In the case of the Cauchy problem, k spans all over R \ {0} and the condition becomes To compare with the Keller-Segel system, it is interesting to consider the case δ = 0. Then a 2 = β and the Routh-Hurwitz array is and the stability condition is the one obtained in the case δ > 0 for b ≥ 0. When k takes values arbitrarily close to 0 (as in the case of the Cauchy problem), it is interesting to find expansions for k → 0. For k = 0, there holds As a consequence, if b > 0, there is only one branch λ = λ(k) such that λ(0) = 0 for both δ > 0 and δ = 0. Such branch has a second-order expansion that is independent of δ and it has a parabolic scaling For large frequencies, the behavior for δ = 0 differs from the one with δ > 0. Indeed, looking for the expansion for the three branches λ = λ(k) as k → ∞, we get In the case δ = 0, p degenerates to a polynomial of degree two. For γ 2 /β = D, asymptotic expansions for the corresponding branches are Such expansions show that, while for the hyperbolic model, it is expected the presence of (decaying) Dirac's delta terms in the fundamental solutions, in the Keller-Segel system solutions to linearized problem are smoothed by the presence of diffusion in both equations.
3. The Cauchy problem: Stability of the zero solution. The aim of this section is to investigate existence of global solutions to the Cauchy problem on the whole real line for system (1), in the case D > 0, and their decay to zero, when initial data are sufficiently small in some suitable norms. Let us remark immediately that we are not able to prove the stability of general constant states under the broader assumption (14). Additional comments are given in Remark 6 below.
A sharp result of local existence of solutions is crucial for our proof; for this reason, although the proof is rather standard, this result is presented below. By a simple rescaling, it is possible to symmetrize the transport matrix and to deal with system (1) with γ = 1. Moreover, we assume f, g, h ∈ C 1 (R 2 ) and we consider the initial conditions with the regularity assumptions Let us set w p = w L p (R×(0,t)) and w(t) p = w(t) L p (R) .
Remark 5. It is easy to show that, assuming regularity on initial data, we obtain regular solutions; in particular, [4]).
Next, we prove existence of global solutions by using a continuation principle.
is bounded, then, by arguing as in the proof of local existence, the same holds for u L ∞ (R×(0,T )) and v L ∞ (R×(0,T )) .
, v L ∞ (R×(0,T )) } and let > 0 be the maximal time of existence for solutions of the Cauchy problem (1), (16), (17) as initial data for a Cauchy problem with maximal time of existence T ≥ .
To prove the global existence of solutions for small initial data, we make some assumptions on the functions f, g, h we are led to consider the system This system is supplemented by the initial conditions As a consequence of (H g ), The proof of the global existence of solutions is built up on the precise decay estimates obtained in [3] for the Green function of the linearization of the hyperbolic part of (19). For this reason, we need some estimates, involving the decay in time of L ∞ and L 2 norms of solutions. Given α > 0, we introduce the following quantities We collect the estimates referred to function φ in the following proposition.
Moreover, if K is sufficiently small, then we have Proof. For the function φ the following expression holds For future references, we recall that the heat kernel Γ d and its derivatives satisfy As a consequence, by properties of function f and (26), we have Now, by standard inequalities (see for instance [3], Lemma 5.2) we have and (21) follows. By using (25), we also obtain and by arguing as above we obtain the inequality Inequalities (22)-(23) can be proved similarly by differentiating (25).
Next, let G : R 2 → R 2×2 be the Green kernel of the dissipative hyperbolic system In [3] it is proved that G can be decomposed as follows, denoting by χ I the characteristic function of the set I, and for some suitable values of the constants.
to be given, we can write the components (u, v) of the solution of (19) by means of the Duhamel formulas and take advantage of the previous estimate to prove the existence of global solutions to system (19), for suitably small initial data.
then there exists a unique global solution to the Cauchy problem (19)-(20) satisfying Moreover, all of the quantities M ∂xφ are uniformly bounded for t ∈ (0, +∞).
Proof. Fix K > 0 large enough and let T > 1. Consider a solution to system (19) By representing u thanks to Duhamel formula relative to the hyperbolic part of system (19), we obtain having used the expressions for G ij , the bounds for |R ij |, and (26). By Proposition 1, the first integral in (30), denoted below by I 1 , can be estimate by where M u and M v are calculated at time t. By Lemma 5.2 in [3] we obtain Next, considering the second integral in (30), denoted by I 2 , we obtain Finally, for t > δ > 0, there holds where The next step is the estimate of N u ; taking into account the expressions for G ij , the bounds for |R ij | and (26), we have By arguing as for the estimate of M 1 4 u , we obtain Similarly, we prove that inequalities (31) and (32) hold also for M where A ∞ is a positive constant depending on K, B d ∞ is a positive constant depending on K and on the data and C d is also a positive constant depending on the data.
Formula (33), for suitably small data, implies that M 1 4 On the other hand, when t > 1, this implies that u(t) ∞ , v(t) ∞ do not increase with t. Thanks to Proposition 1, the same holds for φ(t) ∞ , ∂ x φ(t) ∞ , then the claim follows from the continuation principle. In particular, it is easy to verify that a global estimate for φ 1,∞ implies global estimates for u 1 and v 1 .
Remark 6. The technique used in the above proof of stability for the zero constant state fails when one attempts to prove the stability of non-zero constant stationary states, even if small: in those cases, in the estimates of M 1 4 we have to consider also the linear term U ∂ x φ, which does not present enough polynomial decay. This is not surprising, since our argument is based on a precise representation of the solution to the linear system without the term U ∂ x φ. The linear stability analysis of the full operator, including this term, is absolutely not trivial in the Cauchy case and beyond the aims of the present paper. However, in section 4, this analysis will be done in the case of a bounded domain and Neumann boundary conditions.
In the last part of this section we compare the large time behavior of the solution u to system (19) with the solution of the parabolic Keller-Segel model where the functions f, g satisfy the assumptions (H f ), (H g ).
It is known that, for small initial data, the solutions to the above problem decay in time in L ∞ and L 2 norms in the same way as the solutions to problem (19), see for instance [28] and references therein. We will show that, under the assumption of small initial data, if Theorem 3.4. Let (u, v, φ) and (ũ,φ) be the global solutions respectively to system (19) and to system (34), under the assumptions (H f ), (H g ), (H h ) and (35). Then, there exist 0 , L > 0 such that, if The difference between u andũ can be expressed as follows It is easy to verify that, for suitable large c 1 , c 2 , for |x| > t and t > 1, we have Γ and R 11 verify the bound in (27) and (28). Also, it is easy to verify that, for suitable large c 1 , c 2 , for |x| > t and t > 1, where K h 12 and R 12 verify the bound in (27) and (29). Now, from (36), for t > 1, we have Proceeding as in the proof of Theorem 3.3, we are able to estimate u −ũ 2 for large t: Arguing as in Proposition 1 it is easy to show that where C K is a constant depending on K; then, by using the known decays of the L 2 and L ∞ norms of u,ũ, φ,φ, ∂ x φ, ∂ xφ , we obtain 4. The Neumann problem: Stability of constant states. In this Section, we consider the problem of asymptotic stability of constant states (U, 0, Φ), with f (U, Φ) = 0 under suitable smallness condition on U , as analyzed in the Section 2, in bounded intervals. If f ∈ C 1 and b := −f φ (U, Φ) > 0, by the Implicit Function Theorem, we deduce that there exists a branch of constant solutions Φ =Φ(Ũ ) for U close to U . Precisely, given L > 0, we turn to consider system with Neumann boundary conditions for u and φ, i.e.
Assuming the system to be satisfied until the boundary, the boundary values V 0 (t) := v(0, t) and V L (t) : In particular, if V i (0) = 0, then V i (t) = 0 for any t > 0.
Assume f, g ∈ C 2 , h ∈ C 1 . If (U + u, v, Φ + φ) denotes a perturbed solution, the perturbation (u, v, φ) satisfies where where χ, β, a, b have been previously introduced at equation (11). Here, the symbol F = O(|w| α ) means that for any compact set K there exists C = C(K) such that |F | ≤ C|w| α for any w ∈ K. In what follows, we assume functions g and h together with their derivatives to be globally Lipschitz with respect to ψ so that the constant C can be chosen depending only on the L ∞ norm of u and φ. System (38) is completed by initial and boundary conditions Since the mass of u is conserved, decay of the perturbation can be expected only if If L 0 u 0 (x) dx is small, we can consider as unperturbed constant stateŨ := U + 1 L L 0 u 0 (x) dx, so that condition (40) can be assumed without loss of generality, as soon as |u 0 | L 1 is sufficiently small. Setting w = (u, v, φ), system (38) can be rewritten as where operators L and F are defined by together with initial/boundary conditions given in (39). Under the stability condition formally introduced in Section 2, we can prove that constant states are asymptotically stable solutions (with perturbations decaying exponentially fast) of the Neumann initial-boundary value problem. Theorem 4.1. Assume f, g ∈ C 2 , h ∈ C 1 with g and h satisfying (9) and let (U, 0, Φ) be a constant steady state of (37) such that (43) Assume the stability condition Let w 0 = (u 0 , v 0 , φ 0 ) ∈ H 1 be such that (40) holds, and w the corresponding solution to (41) with initial/boundary conditions given in (39). Then there exists for some C, θ > 0.
The proof of the Theorem is based on a detailed analysis of the Green distribution of the linearized problem (see Theorem 4.2) and on the use of the Duhamel formula to deal with the nonlinear problem.
In the case of functions f and g satisfying the symmetry assumption described in Remark 2, we can relax, for symmetric perturbations of constant states, the stability condition (44) of Theorem 4.1 .
Corollary 1. Let f, g, h be such that (43) holds and assume, additionally, that condition (10) is satisfied. Let (U, 0, Φ) be a constant steady state such that Then the conclusion of Theorem 4.1 holds for initial data w 0 = (u 0 , v 0 , φ 0 ) which verify the assumptions of this theorem and such that Proof of Corollary 1. The symmetry assumption on the initial datum w 0 together with (10) guarantee that the same symmetry with respect to the middle point As a consequence w (0,L/2) and w (L/2,L) are solutions to the same Neumann initial-boundary value problem in the intervals (0, L/2) and (L/2, L), respectively. Since the solution w (0,L/2) satisfies the assumptions of Theorem 4.1 with respect to the interval (0, L/2). Hence the conclusion follows from (45) and from the symmetry of the solution.
Study of the linearized problem. First of all, let us consider the linearized problem ∂ t w + Lw = 0 (47) with initial/boundary conditions given in (39). Following the analysis performed in Section 2, we assume the stability condition (44). We aim to study the Green distribution of the problem, i.e. a distribution-valued map (x, t) → G = G(x, ·; t) such that the solution of (47)-(39) is given by w(x, t) = G(x, ·; t), w 0 .
Considering an even extension of u 0 , φ 0 to (−L, L) and an odd extension of v 0 to (−L, L) and then a 2L−periodic extension of the full initial condition w 0 , we can restrict our investigation to the case of 2L−periodic data. So, we look for solutions of the form w(x, t) = G per (x, ·; t), w 0 , where the Green distribution G per (x, ·; t) has to be determined. Now, we express the solution by means of Fourier series which yields a system of ODEs for the coefficients W k = W k (t): Here A(k) := A 0 + iA 1 k − A 2 k 2 and  If the convergence of the series is sufficiently strong to exchange the summation with the integral, the final form for the Green function G per for the periodic problem turns to be where the convergence of the series is intended in the distributional (bounded measure) sense, since in general we expect a singular part (Dirac mass) arising from the hyperbolic part of the differential operator. For simplicity, from now on, we denote K L with k. The stability condition (44) implies that all of the eigenvalues of A(k) have positive real part for k = 0. As a consequence, the contribution of such terms in the series for the Green function given by (48) is expected to be exponential decaying in time. Hence, it makes sense to decompose the Green function as  The final step consists in analyzing the behavior for k large. Let us determine the eigenvalues expansions as k → ∞. Using expression (12) of p(λ, k) found in Section 2, we look for eigenvalues expansion of the form where Λ(k) := diag(λ − (k), λ + (k), λ D (k)). Inserting in (12), we find the coefficients of the expansion where η := 1 D aχU + 1 2 β 2 . Since γ = 0, for large k the matrix A(k) is diagonalizable with diagonalizing matrix C = C(k). Hence, for large k, e −A(k) t = C(k) e Λ(k) t C(k) −1 .
Next, let us look for the expansion of the matrix C = C(k). First we expand the right eigenvalues of the matrix −A: r(k) = r 0 + k −1 r −1 + O(k −2 ). Plugging in the relation (A(k) + λ(k))r(k) = 0 with λ(k) = c 2 k 2 + c 1 k + c 0 + O(k −1 ), equating to zero the coefficient of powers of k, and using the explicit expressions for A 0 , A 1 , A 2 , we deduce that Similarly, for the branch λ D (k), we have Thus, for k → ±∞, we can choose eigenvectors such that Correspondingly, the expansion for C(k) −1 is Thus we recover the expansion for the matrix e −A(k) t where we used the notation A, B = AB − BA. Moreover, setting Λ ∞ (k) := Λ 2 k 2 + Λ 1 k + Λ 0 , we have Hence, we deduce the following asymptotic expansion as k → +∞ (50) With all of these expansions, we can prove the following results. Note that, in order to obtain estimate on G per itself, only the expansion to order O(k −1 ) is needed; the detailed form of the k −1 is needed to estimate the x−derivative of G per .
Theorem 4.2. The Green distribution G per of the linear problem (47) with periodic initial data and its x−derivative ∂ x G per can be decomposed as where η := 1 D aχU + 1 2 β 2 , and the functions R i = R i (x, y; t), i = 0, 1 satisfy R i (·, y; t) 2 ≤ C e −θ t for t > 0 and some constants C, θ > 0.
Proof. Thanks to (50), the decomposition in formula (49) can be refined as follows: for some θ > 0, where Fixed t, the term G 2 defines an L 2 function of the difference x − y with Fourier coefficients G k 2 with respect to the orthogonal basis {e i(x−y) k } k∈Z given by G 0 2 = 0 and G k 2 = O(k −1 ) e Λ∞(k) t for k = 0. Since Λ 1 is purely imaginary and Λ 2 is non positive, there hold |G 0 2 | 2 = 0 and |G k 2 | 2 = O(k −2 ) e −β t for k = 0. As a consequence, the term G 2 can be absorbed in the term R 0 of the statement. The term f D is smooth for every t > 0, thanks to the presence of the fast decaying term e −D k 2 t . Since we are interested in the values of k = 0, this term can be absorbed in the term R 0 given in (51).
Next, we prove (52). From (50), formula (48) can be written as: where The term G 1 has been discussed above. Differentiation with respect to x gives raise to the sum of a distributional derivative of the singular term G sing and a smooth term whose L 2 norm is exponentially decaying in time. The term G 3 can be easily dealt with. Indeed, after differentiation with respect to x, we Finally, let us consider the term G 2 . After differentiation with respect to x, we get ∂ x G 2 (x, y; t) = H 1 (x, y; t) + H 2 (x, y; t) where completing the analysis of the term H 1 . Concerning H 2 , there holds Collecting the terms, we obtain (52).
Let us denote by G sing = G sing (x, t), H sing = H sing (x, t) and R i = R i (x, t), i = 0, 1, the operators with kernels G sing , H sing , R i in the decomposition given by Theorem 4.2. We have the following decay results.

Corollary 2.
Under the assumptions of Theorem 4.2, we have for some C, θ > 0, Nonlinear stability. In order to prove nonlinear stability using the linear analysis just performed, we need some basic energy estimates.
Then there exists C, θ > 0 such that Proof of Lemma 4.3. Multiplying the equation for φ by e θ t ∂ xx φ and rearranging, we get . Hence, integrating with respect to (x, t) in I × [0, t], we get, for θ sufficiently small, Applying Young inequality, we get (56).
Next we consider the nonlinear stability problem.
Proof of Theorem 4.1. Let w = (u, v, φ) be the solution of with L and F defined in (42). Let T (t) w 0 be the solution to the linear problem ∂ t w = Lw with initial datum w 0 . The operator-valued function t → T (t) can be represented by means of the Green distribution G, decomposed in Theorem 4.2, as the sum of a projection term Π/2L, a singular part (composed by exponentially decaying Dirac's delta functions) and an exponentially decaying smooth part. The solution to (57) is implicitly given by the Duhamel formula.
By means of decomposition (51), the solution w satisfies since ΠF = 0 and (40) holds. Secondly, we have Taking the L 2 −norm on both sides and applying (55), we obtain e θ t w 2 (t) ≤ C w 0 2 + C t 0 e θ s F(w)(s) 2 ds.
Differentiating (58) with respect to x, and taking the L 2 norms, we get Summing up, we deduce By the structure of F, it follows F(w) To get rid of the term ∂ xx φ we use (56). For |(u, φ)| ∞ sufficiently small the term containing |∂ xx φ| 2 2 in the right-hand side of (56) can be absorbed in the left-hand side. Hence, the following holds Using the above relation, from (60), we deduce Let η(t) := sup s∈[0,t] e θ s w H 1 (s). Then from (61) it follows The estimate (62) implies that there exists some ε 0 > 0 such that, if η(0) ≤ ε 0 , then η(t) ≤ C η(0) for some C > 0 and for any t > 0.

Stationary solutions and bifurcation diagram.
The aim of this section is to describe the whole set of stationary solutions for system (4) considered on the interval [0, L] with no-flux boundary conditions v(0, ·) = v(L, ·) = 0 and ∂ x φ(0, ·) = ∂ x φ(L, ·) = 0 and with the special choice The first equation in (4) implies v to be constant. Assuming no-flux boundary conditions, it leads to v = 0 and the system satisfied by the steady states reduces to Such system is the same of the one satisfied by steady solutions of the Keller-Segel system (5) with Neumann boundary conditions ∂ x u(0, ·) = ∂ x u(L, ·) = 0 and ∂ x φ(0, ·) = ∂ x φ(L, ·) = 0. A precise description of bifurcation patterns for such solutions has been considered in [35]. From the first equation in (63), we deduce u = C exp γ −2 φ , C ∈ R, so that the component φ satisfies The function φ satisfies therefore the equation which enables us to draw the phase diagram (see Figure 1). It shows that according to the length of the curve, non constant steady states could exist. In fact, system (4) supports also non-homogeneous stationary solutions. In the following, all these stationary solutions (u, φ) would be displayed as points of a same graph with coordinates ( L 0 u(x) dx, φ(0)). On this graph, the constant stationary solutions would form the line ax − by = 0 (see figure 2).
More precisely, according to [35], for all k ∈ N * , there exists a bifurcation branch B k of non constant stationary solutions starting, in our case, at constant steady The derivatives of both solutions u and φ of branch B k get null at exactly k + 1 points of [0, L]. In [35], more precise results on the form of branches are also given; for example, two different branches do not cross.
Numerically, we obtain the following bifurcation diagrams by a shooting method: starting in the neighborhood of point (C k = U k exp γ −2 Φ k , Φ k ) of the constant stationary solutions line, we solve the problem (64) for all suitable C. We change the values φ 0 until obtaining a solution φ such that ∂ x φ(L, ·) = 0. Once the solution φ is found, we compute the corresponding function u = C exp γ −2 φ and its mass L 0 u(x)dx. We can therefore put the corresponding point on the diagram, as can be seen on figure 2, where we exhibit the two first branches B 1 and B 2 . Let us now discuss the stability of all these stationary solutions. Since the quantity L 0 u(x, t)dx is preserved in time, stability has to be intended as orbital stability, that is to say with perturbationsũ of function u such that L 0ũ (x)dx = 0. In [35], it is proved that, for the Keller-Segel equation (5), constant steady states are stable under the first bifurcation point, that is to say under condition (44). In the previous section, we have shown that the same result holds true for the hyperbolic system (4).
In [35], a first result on stability for non-constant stationary solutions is also given, that is to say that the first bifurcation branch B 1 is stable if and only if b < 2π 2 /L 2 . We were unable to give here an analytic proof of the analogous result for hyperbolic system (4).
However, we observe numerically that in the case where b ≥ 2π 2 /L 2 , the form of bifurcation diagram around the first bifurcation point changes as seen in figure  3. As explained in more details in section 7, the same stability holds for both parabolic and hyperbolic models; indeed, in the case b < 2π 2 /L 2 , B 1 is stable and the other branches are unstable as shown in figure 2. In the case b ≥ 2π 2 /L 2 , B 1 is unstable around the first bifurcation point and then stable and the other branches are unstable, as shown in figure 3.
Moreover, according to Remark 2, special attention has to be drawn on symmetric initial data. For example, in the case b < 2π 2 /L 2 , that is to say, the case of figure 2, if we begin with a symmetric initial datum u which mass is such that u(x)dx < U 2 , the solution will tend to constant steady state for long time, although the constant steady state is unstable. This is in full agreement with the conclusions of Corollary 1. Also, it should be noticed that this is not in contradiction with the stability of the non-constant steady state, since none of the admissible perturbations (i.e. small enough) added with the non-constant steady state would give a symmetric function. We can also understand this result in this way: since the symmetry is preserved with time and the stable stationary solution of first branch B 1 is not symmetric, the solution can not tend towards the stable steady state.
6. Numerical methods. The aim of these last sections is to give numerical evidences of results of Section 4, but also to explore numerically the behaviors of solutions in other cases where the asymptotic behavior has not yet been established by analytical tools. That is to say, we will confirm the stability of constant stationary solutions before the first bifurcation point and the instability after it. We also aim at finding the stability for the hyperbolic model of non constant stationary solutions mentioned in Section 5. For that purpose, let us introduce an appropriate hyperbolic numerical scheme. We consider the system (4) on the interval [0, L] with and with homogeneous Neumann boundary conditions for φ and homogeneous Dirichlet boundary conditions for v. These two conditions and the second equation of (4) lead to homogeneous Neumann conditions for u on the boundaries. In order to discretize it, we write it in more convenient form with diagonal variables w = 1 2 (u − v/γ) and z = 1 2 (u + v/γ), which correspond respectively to the functions u − and u + of Section 2. The system becomes Let us denote by δx the space step and by δt the time step. We consider the discretization points x j = j δx, 0 ≤ j ≤ N + 1, with x 0 = 0 and x N +1 = L.
The discretization times will be given by t n = n δt, n ∈ N. Let us denote by w n j (resp. z n j and φ n j ) the approximation of w(x j , t n ) (resp. z(x j , t n ) and φ(x j , t n )). The discretization vectors at time t n will therefore be W n = (w n 1 , . . . , w n N ) t , Z n = (z n 1 , . . . , z n N ) t and Φ n = (φ n 1 , . . . , φ n N ) t . We first discretize the two first equations of system (65) by using upwind explicit methods and the third one using Crank-Nicolson scheme in time and finite differences method of order two in space. The boundary conditions are treated as follows: v(0, t) = v(L, t) = 0 lead to z n 0 = w n 0 and w n N +1 = z n N +1 in discrete (w, z) variables. Since w n+1 0 and z n+1 N +1 are computed thanks to the upwind method, we simply compute w and z on boundaries as The function φ satisfies homogeneous Neumann boundary conditions and we use a finite difference discretization of order two of its derivative, which gives: This enables us to compute ∂ x Φ n whose j-th component is the centered approximation of order two of ∂ x φ(x j , t n ). This vector is needed for the discretization of the two first equations of system (65).
Consequently, let M be the classical finite differences N × N matrix, using discretizations (66) for the computation of the first and the last raws. We consider the following scheme (67) Since the spectrum of M is contained in the disk D(2, 2) = {λ ∈ C : |λ − 2| ≤ 2}, the matrix 1 + b δt 2 I + δt 2δx 2 DM is invertible without any condition on δt and δx.
Such scheme preserves some properties of the original system: constant stationary solutions, conservation of mass, symmetry with respect to the transformation (x, t, u, v, φ) → (L − x, t, u, −v, φ). It also gives some results for the functions u and φ, which are in good adequacy with the values obtained on bifurcation diagrams at Section 5.
However, the function v is not null at equilibrium as predicted by the continuous analysis. Indeed, rewriting scheme (67) in the initial variables u and v, we find the following equation for u n : We therefore see that when u reaches its equilibrium and when the gradient of u is large (which happens frequently for non constant steady state solutions), the space step δx needs to be very small in order to cancel function v.
Such kind of problems are not surprising when dealing with the numerical approximation of stationary solutions.
Even in the context of hyperbolic chemotaxis, some attempts were made to give an exact resolution of steady states, see [10] and references therein. However, those papers were more devoted to the study of the Cauchy problem and the ideas they proposed do not seem relevant to the present framework, where the treatment of the boundary conditions is more crucial. Here, we follow the ideas of [2] and we present a scheme which is more accurate for stationary solutions. In a following work [30], we will give more details on the numerical problems related to the simulation of this system. Now, let us denote by f the right hand side φ x u and in the following we shall omit the parabolic equation for φ which will be treated as before. We therefore consider the following system for x ∈ [0, L] with the boundary conditions v(0, ·) = v(L, ·) = f (0, ·) = f (L, ·) = 0. We reconsider the system in diagonalizing variables, that is to say Let us denote by ω = w z and we rewrite the system in the following form, see [2,30]: We denote by ω n i an approximation of ω(x i , t n ) and we consider schemes of the following form: In [30], we shall give conditions to ensure monotonicity, consistency of this scheme, third-order accuracy of the truncation error when restricted to stationary solutions. We also compute the matrices B and D in order to preserve mass, symmetry and constant stationary solutions. Here we just present the final scheme, which is given by q = γ and We also have the following restrictions on the time and space steps, which guarantee the monotonicity of the scheme [30]: h ≤ 24γ 11 In Figure 4, we can compare the computations for functions u and v at final time (i.e. when u has reached its equilibrium) with the two schemes (67) and (68)-(69) for the same parameters and the same initial data. On the left, the computed functions u are almost the same and in sharp agreement with the expected result; however, we can see on the right the difference between the computations of function v, which is expected to be equal to zero. While the error for the first scheme is quite large, it almost vanishes for the second scheme. 7.1. Asymptotic behavior of system (65). The first aim of this section is to present some numerical simulations in agreement with the previous theoretical results and with the stability results announced in Section 5. In this first part, stability will always mean stability for the system (4). In reality, we think that stability for system (4) and for system (5) are the same and we shall give a few numerical evidences of this fact in subsection 7.2. Let us begin with the analysis of the bifurcation diagram in Figure 2, which displays the line of constant steady states and the two first bifurcation branches B 1 and B 2 . As shown analytically in Section 4, the stationary solutions before the first bifurcation point are stable. Indeed, we check it on Figure 5 since any solution with a mass smaller than the first bifurcation point tends to the corresponding stable constant steady state. Now, let us see what happens for solutions with a mass between the first and the second bifurcation point. Analytically we only know that the constant stationary solution is unstable, but not for symmetric perturbations of a constant state, as stated by Corollary 1. The behaviors are shown in Figure 6: in both cases, the mass of the initial datum L 0 u 0 (x)dx = 1135 is between the values of the two first bifurcation points. The stable steady state is displayed as asymptotic profile in Figure 6(a); we notice that the boundary values obtained for φ, i.e. φ(0) and φ(L) are in good agreement with the values displayed on the bifurcation diagram of Figure 2.
Finally, for solutions with a mass after the second bifurcation point, the unstable constant stationary solution will never be reached asymptotically. In fact, numerical studies displayed in Figure 7 show that, generally speaking, the asymptotic solution will be the profile of branch B 1 , as in Figure 7      perturbations, the final profile will be the symmetric stationary solution of branch B 2 , see Figure 7(b). Now, let us see some numerical results in the case b ≥ 2π 2 /L 2 , corresponding to the second bifurcation diagram in Figure 3. Apart from the neighborhood of the first bifurcation point, the stability results are the same as in the case b < 2π 2 /L 2 . However, around this point, we can see on the enlargement that three stable solutions coexist for a same mass. This is displayed in Figure 8: in (a), the asymptotic profile is the constant stable steady state and in (b), the asymptotic profile is the non constant stable steady state of branch B 1 . However, Figure 9 shows that the two remaining stationary state are unstable, as predicted by [35] (it is the initial data of the evolution displayed in subfigure (a)). In that case, the solution tends asymptotically to the constant steady state. Moreover, a small perturbation of the same unstable stationary solution leads to the other non constant stable steady state. We conclude from these two numerical experiments that the basins of attraction of the two stable stationary solution are difficult to characterize. 7.2. Comparison between hyperbolic and parabolic models. Let us just briefly compare the difference of behaviors between solutions for the hyperbolic model (4) and the parabolic system (5). As displayed in Figures 10 and 11, the asymptotic behaviors of the systems seem to be exactly the same. However, as we can notice in Figure 11, the function u oscillates much more for the hyperbolic equation than for the parabolic equation. Therefore, the function u solution of the hyperbolic equation converges more slowly to the asymptotic state that the function u solution of the parabolic equation.         This is quantified at Figure 12. On this figure, we can see above the convergence speed of function u towards its asymptotic state and below the convergence speed of function φ. More precisely on the subfigure above, we display the quantity ln ||u − u as || ∞ ||u as || ∞ with respect to time and on the subfigure below, the quantity ln ||φ − φ as || ∞ ||φ as || ∞ with respect to time. In red, it is displayed the convergence speed obtained for equation (4) and in blue the one obtained for equation (5). Let We display in red the convergence speed obtained for equation (4) and in blue the one obtained for equation (5).
us remark that the speed of convergence of the hyperbolic model oscillates due to the corresponding oscillations of the solutions. Since we are in the case a = b = D = 1, γ = 10 and L 0 u 0 (x)dx = 10, the asymptotic states for both u and φ are the constant function equal to 10. We can remark clearly that the speed of convergence to equilibrium in model (4) is slower both for u and φ than for model (5). We also notice that after a given time, a residual error remains which is due to the numerical approximation and which therefore depends both on δt and δx. A linear regression enables us to find the coefficients a u,h , C u,h , a φ,h , C φ,h such that for hyperbolic model (4) and the coefficients a u,p , C u,p , a φ,p , C φ,p such that ||u − u as || ∞ ||u as || ∞ = C u,p exp(−a u,p t), ||φ − φ as || ∞ ||φ as || ∞ = C φ,p exp(−a φ,p t) for parabolic model (5). The values of the coefficients are given at Table 1. For the hyperbolic model, they have been computed using only the straight line of upper values in red of Figure 12. However, we notice that for short times, the evolution of solution u for the hyperbolic model (4) is just a power law. In fact, for short times, the solution has not been yet influenced by the boundaries and therefore, the convergence speed is the one of the Cauchy case, in agreement with  (70) and (71) Hyperbolic model (4) Parabolic model (5) Function u a u,h = 0.7408576 a u,p = 4.2311464 C u,h = 0.0004380 C u,p = 0.0052609 Function φ a φ,h = 0.7123929 a φ,p = 10.650477 C φ,h = 0.0000137 C φ,p = 0.0005583 Theorem 3.3, and even better since a linear regression gives that, for short times, ||u − u as || ∞ ||u as || ∞ = C u,h t −a u,h , with a u,h = 0.4899891 and C u,h = 0.0007756.  Finally, Figure 13 shows us the difference of behaviors for solutions u for short times: in case of the hyperbolic model (4), the solution oscillates a lot and the spurious bumps move until the boundaries. For parabolic model (5), the solution also oscillates, but the bumps remain localized in space and the solution dissipates up to the asymptotic state.