F Berthelin

D Chiron

M Ribot

STATIONARY SOLUTIONS WITH VACUUM FOR A ONE-DIMENSIONAL CHEMOTAXIS MODEL WITH NON-LINEAR PRESSURE

In this article, we study a one-dimensional hyperbolic quasilinear model of chemotaxis with a non-linear pressure and we consider its stationary solutions, in particular with vacuum regions. We study both cases of the system set on the whole line R and on a bounded interval with no-flux boundary conditions. In the case of the whole line R, we find only one stationary solution, up to a translation, formed by a positive density region (called bump) surrounded by two regions of vacuum. However, in the case of a bounded interval, an infinite of stationary solutions exists, where the number of bumps is limited by the length of the interval. We are able to compare the value of an energy of the system for these stationary solutions. Finally, we study the stability of these stationary solutions through numerical simulations.

Introduction

The movement of cells, bacteria or other microorganisms following the gradient of a chemical concentration, known as chemoattractant, has been widely studied in mathematics in the last two decades [START_REF] Horstmann | From 1970 until present: the Keller-Segel model in chemotaxis and its consequences[END_REF][START_REF] Perthame | Transport equations in biology[END_REF]. Partial differential equations models have been proposed to describe the complex behavior of such system, whose two main unknowns are the density of cells or bacteria and the concentration of chemoattractant. One of the most considered models is the Patlak-Keller-Segel system [START_REF] Keller | Initiation of slime mold aggregation viewed as an instability[END_REF][START_REF] Patlak | Random walk with persistence and external bias[END_REF], where the evolution of the density of cells is described by a parabolic equation of driftdiffusion type, and the concentration of a chemoattractant is generally given by a parabolic or elliptic equation, depending on the different regimes to be described. The behavior of this system is quite well known now, at least for linear diffusions: in the one-dimensional case, the solution is always global in time [START_REF] Nagai | Large time behavior of bounded solutions to a parabolic system of chemotaxis in the whole space[END_REF], while in two and more dimensions the solutions exist globally in time or blow up according to the size of the initial data, see [START_REF] Calvez | The parabolic-parabolic Keller-Segel model in R 2[END_REF][START_REF] Calvez | Blow-up, concentration phenomenon and global existence for the Keller-Segel model in high dimension[END_REF]. However, a drawback of this model is that the diffusion leads alternatively to a fast dissipation or an explosive behavior, while in general we are interested in the creation of patterns, like in the vasculogenesis process.

keywords: chemotaxis, stationary solutions, vacuum, quasilinear hyperbolic problem, energy. 1 Laboratoire J. A. Dieudonné, UMR 7351 CNRS-Université de Nice-Sophia Antipolis, Parc Valrose, F-06108 Nice Cedex 02, France & Project Team COFFEE, INRIA Sophia Antipolis, France (bertheli@unice.fr, ribot@unice.fr). 2 Laboratoire J. A. Dieudonné, UMR 7351 CNRS-Université de Nice-Sophia Antipolis, Parc Valrose, F-06108 Nice Cedex 02, France (chiron@unice.fr).

In order to improve the modeling of the phenomena we deal with, two kinds of modifications of the Keller-Segel equation can be considered. The first one is to introduce a non linear pressure, as in [START_REF] Kowalczyk | Preventing blow-up in a chemotaxis model[END_REF][START_REF] Calvez | Volume effects in the Keller-Segel model: energy estimates preventing blow-up[END_REF]. In that case the Keller-Segel-like system with nonlinear diffusion reads as:

∂ t ρ = ∂ x (∂ x p(ρ) -χρ∂ x φ), δ∂ t φ = D∂ xx φ + aρ -bφ, (1.1)
where ρ is the density of cells, φ is the concentration of chemoattractant and χ, D, a and b are given positive parameters. If δ = 1, we consider a parabolicparabolic model and in the case where δ = 0, we deal with a parabolic-elliptic model. The coupling between the two equations can be described as follows: the drift term χ(ρφ x ) x stands for the movement of the cells towards the gradient of concentration of the chemoattractant, whereas the source term aρ indicates that the chemoattractant is produced by the cells themselves. The function p is a phenomenological, density dependent function, which is given by the pressure law for isentropic gases

p(ρ) = κρ γ , γ > 1, κ > 0, (1.2)
which takes into account the fact that cells do not interpenetrate since they have a non zero volume. Overcrowding may also be prevented by the modeling of a volume-filling effect, see [START_REF] Hillen | Global existence for a parabolic chemotaxis model with prevention of overcrowding[END_REF][START_REF] Painter | Volume-filling and quorum sensing in models for chemosensitive movement[END_REF]. Let us also mention some very recent works on a similar subject [START_REF] Burger | Stationary states and asymptotic behavior of aggregation models with nonlinear local repulsion[END_REF] where the authors study the steady states and their stability of a Keller-Segel type system with a more general potential. An adapted numerical scheme to this problem is developed in [START_REF] Carrillo | A finite-volume method for nonlinear nonlocal equations with a gradient flow structure[END_REF] using the gradient flow structure and interesting numerical simulations are presented.

The second type of modifications consists in studying hyperbolicparabolic or hyperbolic-elliptic models, which are more likely to show fine transient behaviors. We consider consequently a quasilinear hyperbolic model of chemotaxis introduced by Gamba et al. [START_REF] Gamba | Percolation, morphogenesis, and burgers dynamics in blood vessels formation[END_REF] to describe the early stages of the vasculogenesis process. In this model, three unknowns are present: the density of cells ρ(x,t), their momentum ρu(x,t) and the concentration φ(x,t) of the chemoattractant. The model reads as

   ∂ t ρ + ∂ x (ρu) = 0, ∂ t (ρu) + ∂ x (ρu 2 + p(ρ)) = χρ∂ x φ -αρu, δ∂ t φ = D∂ 2 xx φ + aρ -bφ, (1.3)
where the constants are D,χ,α > 0, a,b ≥ 0. The term αρu represents the damping force and the pressure p is still given by the pressure law for isentropic gases (1.2). This model of chemotaxis has been introduced to describe the results of in vitro experiments performed by Serini et al. [START_REF] Serini | Modeling the early stages of vascular network assembly[END_REF] using human endothelial cells which, randomly seeded on a matrigel, formed complex patterns with structures depending on the initial number of cells. However, as far as we know, only a few analytical results exist for system (1.3). In [START_REF] Russo | Analysis and numerical approximations of hydrodynamical models of biological movements[END_REF][START_REF] Russo | Existence and asymptotic behavior of solutions to a quasilinear hyperbolic-parabolic model of vasculogenesis[END_REF], the authors prove the global in time existence of solutions if the initial datum is a small perturbation of a small enough constant stationary solution. This proof is valid for the Cauchy problem set on the whole line R, when the vacuum is not reached. The main difficulty relies in the fact that some vacuum regions, i.e. where the density vanishes, appear in finite time and that the hyperbolic system (1.3) becomes degenerate in such a case. Although vacuum regions can be avoided in gas dynamics models, it is of great importance to have such a possibility in models for biology, where some regions may be clear of any cells, as in the vasculogenesis process. Some related results are given in [START_REF] Jang | Well and ill-posedness for compressible Euler equations with vacuum[END_REF] about the local existence of solutions for the compressible Euler equations with damping and vacuum, but without chemotaxis or in [START_REF] Gu | Well-posedness of 1-D compressible Euler-Poisson equations with physical vacuum[END_REF] where they study the 1D compressible Euler-Poisson equations with moving physical vacuum boundary condition in order to to describe the motion of a self-gravitating inviscid gaseous star. However, many articles deal with the study of stationary solutions with vacuum of the Euler-Poisson system in the context of polytropic gaseous stars [START_REF] Deng | Solutions of Euler-Poisson equations for gaseous stars[END_REF][START_REF] Rein | Non-linear stability of gaseous stars[END_REF].

As for the numerical point of view, some interesting numerical simulations in 2D and 3D can be found in [START_REF] Andreianov | Finite volume methods for degenerate chemotaxis model[END_REF] or [START_REF] Cavalli | 3D simulations of early blood vessel formation[END_REF]. In [START_REF] Natalini | A well-balanced numerical scheme for a one dimensional quasilinear hyperbolic model of chemotaxis[END_REF], an adapted numerical scheme based on unwinding technique was developed in order to compute efficiently the time behavior of the solutions of system (1.3) in 1D. This scheme was proved to preserve the positivity of solutions and a discrete version of stationary solutions of (1.3). Some numerical simulations for different values of γ show that we find asymptotically a wide range of stationary solutions and that most of them seem stable. This numerical scheme was improved in [START_REF] Natalini | A numerical comparison between degenerate parabolic and quasilinear hyperbolic models of cell movements under chemotaxis[END_REF] considering the large time-large damping limit of (1.3) towards the Keller-Segel system (1.1) with δ = 0, leading to an asymptoticpreserving scheme. Some numerical simulations enable to compare the asymptotic behavior of the hyperbolic system (1.3) and the asymptotic behavior of the parabolic system (1.1) and surprisingly, these behaviors may be different. However, the validation of the scheme is not so obvious and a first computation of some stationary solutions with vacuum in the case γ = 2 was made in [START_REF] Natalini | A well-balanced numerical scheme for a one dimensional quasilinear hyperbolic model of chemotaxis[END_REF] in order to validate the numerical simulations.

In this article, we push forward these computations in order to have a very precise idea of all the possible configurations of stationary solutions. We also compute an energy for all these solutions and we compare them. More precisely, we consider stationary solutions for (1.3) either on the whole line R or on a bounded domain Ω ⊂ R with no-flux boundary conditions

∂ x ρ = 0, ∂ x φ = 0, ρu = 0 on ∂Ω.
In both cases, the momentum ρu vanishes and the two density functions ρ and φ would satisfy the following system, defining β ≥ 0 as b = Dβ 2 :

∂ x p(ρ) = χρ∂ x φ, D∂ 2 xx φ + aρ -Dβ 2 φ = 0. (1.4)
We also impose that the densities remain positive, namely

ρ,φ ≥ 0 (1.5)
and that they have the following regularities

ρ ∈ C 0 (Ω), φ ∈ C 2 (Ω).
Some constant in space stationary solutions are obvious solutions to equation (1.4). Now, let us consider stationary solutions with vacuum.

In the following, we will call "bump" a region with a nonnegative density surrounded by two regions of vacuum and "half bump" will be a bump cut in its middle and stuck to an extremity of the interval.

Unless in the degenerate case γ = 1 [START_REF] Filbet | Approximation of hyperbolic models for chemosensitive movement[END_REF], where the solutions remain always strictly positive, the value of the parameter γ is not clear in the biological context and we have to keep all the values of γ strictly larger than 1 under consideration, since vacuum may appear for all these cases. Note that the previous mentioned numerical results are valid for different values of γ > 1. However, we mainly focus in this article on the pressure (1.2) for γ = 2, which is the only case with exact and tractable analytical expressions for the solutions of system (1.4). In particular, the first equation of (1.4) simplifies as ǫρ∂ x ρ = χρ∂ x φ which leads to ρ = 0 or ερ = χφ -K with K constant. These two relations give two types of intervals: the ones where ρ is null and the ones where ρ is strictly positive. We may then find some conditions to construct a solution with a transition between these two types of intervals. Among them, the following relation between the parameters of the system is necessary: aχbε > 0. Otherwise, when aχbε ≤ 0, only the already mentioned constant solutions exist. From these conditions, we may prove the following results.

First, on the whole line R, once the mass of the solution is fixed, there exists a unique stationary solution up to a translation, made of one bump. Then, on a bounded interval [0,L], things are much more complicated: first, once the length L of the interval is fixed, only a finite number of configurations are possible. More precisely, we can construct (at least) one k half bumps solution under the following condition linking the length of the domain and the parameters of the system: L ≥ kπ ω , where ω 2 = aχ Dε β 2 . Notice that, in this formula, a bump corresponds to two half bumps. Now, if the length of the domain is large enough and once the mass of the solution is fixed, we can construct two half bumps, one on the right of the interval and one on the left and we can also construct a unique bump. However, a striking fact is that there exists a continuum of two bumps solutions, parametrized by the length of vacuum in the interval. For the two bumps solution, we can also find a continuum of two half bumps solutions or of one bump and a half bump solutions. Some of these stationary solutions are displayed on Figure 3.

Finally, we may compare analytically the energy of all these stationary solutions and we show that the constant solution has the larger energy, whereas the half bump solution, when it exists, has the smallest energy and we expect, therefore, this solution to be the stable one. The numerical simulations presented in the last section exhibit some sets of parameters for which the half bump seems to be indeed the only stable solution. However, for some other sets of parameters, a wide range of stationary solutions look stable in our numerical experiments.

This article is organized as follows: in Section 2, we give two preliminaries for the following study, that is to say we define a suitable energy for system (1.3) and we also study the condition to have a transition between a region with a positive density and a region with a null density. In Section 3, we are then able to prove that, in the case of the whole line R, the only stationary solution is the one composed of one bump. The case of a bounded interval with no-flux boundary conditions is much more complicated and is detailed in Section 4: in this section, we prove the existence of different types of stationary solutions, we compute the energy of each of these solutions and we compare them. Finally in section 5, we present a numerical study of the stability of the previous mentioned stationary solutions.

Preliminaries

We begin this article with some preliminaries which will be useful in the following. In Subsection 2.1, we give an energy estimate and in Subsection 2.2, we study in details under which condition we may find a transition between a region with positive density and a region with a vanishing density.

2.1. Energy estimate. Let us define the energy density

j(ρ,u,φ)(x,t) = 1 2 ρu 2 + Γ(ρ) -χρφ + χDβ 2 2a φ 2 + χD 2a (∂ x φ) 2 , (2.1)
where ρΓ ′′ (ρ) = p ′ (ρ).

(2.2)

In particular for γ = 2, when p(ρ) = ε 2 ρ 2 , we have Γ(ρ) = ε 2 ρ 2 = p(ρ). Energy for a stationary solution. We begin by computing an expression of energy for the stationary solutions, which will enable us later on to compare the stationary states one to another. In the case of a stationary solution satisfying system (1.4) and u = 0, the energy is therefore defined as:

J = Ω Γ(ρ) -χρφ + χDβ 2 2a φ 2 + χD 2a (∂ x φ) 2 dx, (2.3)
where Ω is a subset of R. In the following, we will consider the case Ω = R or the case Ω = [0,L]. Notice that for a H 1 (Ω) solution of (1.4), the following equality holds

0 = Ω φ(D∂ 2 xx φ + aρ -Dβ 2 φ)dx = - R D(∂ x φ) 2 dx + Ω φ(aρ -Dβ 2 φ)dx that is to say Ω χDβ 2 2a φ 2 + χD 2a (∂ x φ) 2 dx = Ω χ 2 φρdx.
Thanks to this last equality, we may simplify (2.3) and we find the following expression for the energy of a stationary state:

J = Ω Γ(ρ) - χ 2 φρdx. (2.4)
Energy estimate. In the following proposition, we prove an energy estimate.

Proposition 2.1. Let us consider some smooth enough functions (ρ,u,φ) solutions of system (1.3) set on Ω ⊂ R. We have the following energy dissipation estimate for all (x,t)

∈ Ω × R + ∂ t j(ρ,u,φ) + ∂ x 1 2 ρu 3 + uΨ(ρ) -χρuφ - χD a ∂ t φ∂ x φ = -ρu 2 - χ a (∂ t φ) 2 ≤ 0, (2.5)
where j is defined at equation (2.1), Γ at equation (2.2) and

Ψ(ρ) = ρΓ ′ (ρ). (2.6)
Proof. From the first equation of system (1.3), we get

∂ t ( 1 2 ρu 2 ) + ∂ x ( 1 2 ρu 3 ) = u ∂ t (ρu) + ∂ x (ρu 2 )
and, thus, the second equation of (1.3) implies

∂ t ( 1 2 ρu 2 ) + ∂ x ( 1 2 ρu 3 ) = u(-∂ x p(ρ) + χρ∂ x φ -ρu). (2.7)
Noticing that

ρu∂ x φ = ∂ x (ρuφ) -φ∂ x (ρu) = ∂ x (ρuφ) + φ∂ t ρ = ∂ x (ρuφ) + ∂ t (ρφ) -ρ∂ t φ, equation (2.7
) can be rewritten as

∂ t ( 1 2 ρu 2 -χρφ) + ∂ x ( 1 2 ρu 3 -χρuφ) = -u∂ x p(ρ) -ρu 2 -χρ∂ t φ. (2.8)
Now, we can compute

∂ t Γ(ρ) + ∂ x (uΨ(ρ)) = (Ψ ′ (ρ) -Γ ′ (ρ))u∂ x ρ + (Ψ(ρ) -Γ ′ (ρ)ρ)∂ x u,
which becomes, thanks to (2.2) and (2.6),

∂ t Γ(ρ) + ∂ x (uΨ(ρ)) = p ′ (ρ)u∂ x ρ.
(2.9)

Combining equations (2.8) and (2.9), we have

∂ t ( 1 2 ρu 2 -χρφ + Γ(ρ)) + ∂ x ( 1 2 ρu 3 -χρuφ + uΨ(ρ)) = -χρ∂ t φ -ρu 2 . (2.10)
Considering now the third equation of system (1.3), we obtain

∂ t χb 2a φ 2 + χD 2a (∂ x φ) 2 + ∂ x - χD a ∂ t φ∂ x φ = χb a φ∂ t φ + χD a ∂ x φ∂ 2 tx φ - χD a ∂ 2 xt φ∂ x φ - χD a ∂ t φ∂ 2 xx φ = χ a ∂ t φ bφ -D∂ 2 xx φ = χ a ∂ t φ(aρ -∂ t φ). (2.11)
Finally, the addition of equations (2.10) and (2.11) gives the energy estimate(2.5), which concludes the proof.

Unfortunately, the variational structure obtained in the previous proposition is not enough to prove some stability results for the stationary solutions we deal with. Indeed, since ρ vanishes on some intervals, the condition ρ,φ ≥ 0 may be saturated and the natural set where we wish to define the energy functional is not an open subset of a linear function space.

Study of the transitions between vacuum and positivity region.

We consider in this subsection a stationary solution (ρ,φ) satisfying system (1.4) for γ = 2, i.e. with a pressure equal to

p(ρ) = ε 2 ρ 2 , (2.12)
that is to say

ερ∂ x ρ = χρ∂ x φ, D∂ 2 xx φ + aρ -Dβ 2 φ = 0, (S)
and we study under which condition a transition may occur between an interval where ρ is positive and an interval where ρ vanishes. The aim is to describe what happens in the case of a stationary "bump", that is to say a region with ρ positive, surrounded by two vanishing regions.

Proposition 2.2. Let M > 0. We consider a solution (ρ,φ) satisfying system (S) on an interval

Ω ⊂ R, with ρ ∈ C 0 (Ω), φ ∈ C 2 (Ω), ρ, φ ≥ 0 and such that Ω ρ(x)dx = M . We assume that ω 2 def = aχ Dε -β 2 > 0. (i) If ρ > 0 on an interval I =]x, x + l[⊂ R, then      φ(x) = -Acos(ω(x -x)) + B sin(ω(x -x)) + aK Dεω 2 , ρ(x) = - Aχ ε cos(ω(x -x)) + Bχ ε sin(ω(x -x)) + β 2 K εω 2 ;
(2.13)

with A,B,K some constants depending on I.

(ii) If ρ = 0 on an interval J =]x -d, x[ ⊂ R, then φ(x) = P + e β(x-x) + P -e -β(x-x) , (2.14)
with P + ,P -constants depending on J.

(iii) If we consider a transition at a point x between J =]xd, x[, where ρ vanishes, and I =]x, x + l[, where ρ is strictly positive, the constants of equations (2.13) and (2.14) satisfy the following inequalities:

A,B,P + ≥ 0, P -≤ P + .

(2.15)

(iv) The maximum length of the interval I =]x, x + l[ where ρ is positive is given by the following formula:

l = 2 ω π -arctan( B A
) .

(2.16)

(v) If we consider now the transition at a point x + l between I =]x, x + l[, where ρ is strictly positive and J ′ =]x + l, ȳ[, where ρ vanishes, and if we denote by (P ′ ) + ,(P ′ ) -the constants such that φ(x) = (P ′ ) + e β(x-ȳ) + (P ′ ) -e -β(x-ȳ) on J ′ , then (P ′ ) + = P -e β(ȳ-x-l) and (P ′ ) -= P + e -β(ȳ-x-l) .

(2.17)

(vi) Finally, the mass of a solution defined by expression (2.13) on an interval I =]x, x + l[ with l given by equation (2.16) is equal to

M = 2Aχ εω π -arctan( B A ) + B A . (2.18)
Proof. (i) Firstly, let us consider the case where ρ > 0 on an interval

I =]x, x + l[⊂ R.
From the first equation of (S), ρ ∈ C 1 (I) and, therefore, using expression (2.12) of the pressure we have on this interval that

p ′ (ρ)∂ x ρ = ερ∂ x ρ = χρ∂ x φ, D∂ 2 xx φ + aρ -Dβ 2 φ = 0. Thus ∂ x (ερ -χφ) = 0,
and

ρ = 1 ε (χφ -K), (2.19)
with K constant (depending on I).

Inserting this expression into the second equation of (S) gives

D∂ 2 xx φ + a ε (χφ -K) -Dβ 2 φ = 0,
that is to say

∂ 2 xx φ + aχ Dε -β 2 φ = aK Dε . (2.20)
It follows that the sign of the coefficient aχ Dε β 2 will be crucial, and we focus in this proposition on the case, denoted by Case (P):

aχ Dε -β 2 = ω 2 > 0.
The general solution of (2.20) is then given by (2.13).

(ii) Secondly, let us consider the case where ρ = 0 on an interval J = ]xd, x[⊂ R. The first equation of (S) is satisfied and the second one becomes ∂ 2 xx φβ 2 φ = 0, whose solution is given by (2.14).

(iii) We pursue the computations in Case (P) and we study the transition between J =]xd, x[ (where ρ vanishes) and I =]x, x + l[ (where ρ is strictly positive) at a point x. From the expressions given at equations (2.13) and (2.14) and using the continuity of ρ,φ,∂ x φ at x, we find the following equalities linking the different constants one to another:

           - Aχ ε + β 2 K εω 2 = 0, -A + aK Dεω 2 = P + + P -, Bω = β(P + -P -), that is to say              A = β 2 χω 2 K, A ω 2 β 2 = P + + P -, Bω = β(P + -P -). (2.21)
Now, let us find some sign conditions for the constants. On the right side of x, ρ is strictly positive, thus ρ ′ (0 + ) ≥ 0 which gives B ≥ 0 and thus P + -P -≥ 0. Moreover, constraints (1.5), namely ρ,φ ≥ 0, imply in particular φ(0) ≥ 0 , that is P + + P -≥ 0 and thus A ≥ 0. Now P + -P -≥ 0 and P + + P -≥ 0 imply P + ≥ 0. Summing up, the conditions are then A,B,P + ≥ 0, P -≤ P + .

(iv) Conversely, we are looking now on which interval ρ remains strictly positive and what is the right extremity of this interval. Using relations (2.21), we may rewrite the expression (2.13) for ρ as

ρ(x) = Aχ ε (1 -cos(ω(x -x))) + Bχ ε sin(ω(x -x)). (2.22)
We distinguish three cases:

Case (I): If A = 0 and B > 0, then ρ(x) > 0 if x ∈]x, x + π/ω[ and ρ(x + π/ω) = 0.
Case (II):

If A > 0 and B = 0, then ρ(x) > 0 if x ∈]x, x + 2π/ω[ and ρ(x + 2π/ω) = 0.
Case (III): If A > 0 and B > 0, simple computations lead to the following expression

ε χ ρ(x) = Asin(ω(x -x)) tan(ω(x -x)/2) + B A , (2.23
)

thus ρ(x) > 0 if x ∈]x, x + π/ω[, whereas (2.13) gives ρ(x + π/ω) > 0. Furthermore equation (2.23) also gives ρ(x) > 0 on ]x + π/ω,x * [ and ρ(x * ) = 0 with ω(x * -x)/2 = π -arctan( B A ).
Remark that, since φ = 1 χ (ερ + K) on I with K ≥ 0, we also have φ > 0 whenever ρ > 0.

Summing up cases (I), (II) and (III), we obtain from this study the length l of the interval I =]x, x + l[ where ρ is positive, namely:

l = x * -x = 2 ω π -arctan( B A ) for A ≥ 0,B ≥ 0.
Remark that from this expression, we recover cases (I) and (II), that is to say l = π/ω if A = 0,B > 0, l = 2π/ω if A > 0,B = 0. (v) Now, let us go further and study the transition at x + l from a region I =]x, x + l[ where ρ is positive to a region J ′ =]x + l, ȳ[ where ρ vanishes. From (2.14), the functions on the region J ′ are of the form ρ(x) = 0 and

φ(x) = (P ′ ) + e β(x-ȳ) + (P ′ ) -e -β(x-ȳ) .
Following the previous technique, we use the continuity of φ and ∂ x φ at point x + l and we obtain

   K χ = (P ′ ) + e β(l-ȳ+x) + (P ′ ) -e -β(l-ȳ+x) , Aω sin(ωl) + Bω cos(ωl) = β(P ′ ) + e β(l-ȳ+x) -β(P ′ ) -e -β(l-ȳ+x) .
(2.24)

Notice that the continuity of ρ gives no condition here. Using expression (2.16), it is straightforward to show that Aω sin(ωl) + Bω cos(ωl) = -Bω Note (see (2.21) and (2.24)) that both couples (X = P + ,Y = P -) and (X = (P ′ ) -e -β(l-ȳ+x) ,Y = P -= (P ′ ) + e β(l-ȳ+x) ) satisfy the linear system

   X + Y = K χ , β(X -Y ) = Bω.
(2.25)

By unicity (β = 0) of the solution to (2.25), we get formula (2.17).

(vi) The final step is to find a formula for the mass of the bump as a function of the constants A and B. From (2.22) and (2.16), a simple computation gives, in the case when A > 0,

M = x+l x ρ(x)dx = χ ε x+ 2 ω (π-arctan( B A )) x A(1 -cos(ω(x -x))) + B sin(ω(x -x))dx = 2Aχ εω π -arctan( B A ) + B A .
It is extended by continuity when A = 0 by:

M = 2Bχ εω . (2.26)
Remark 2.3. In the following, we call Case (P), the case of the previous proposition where aχ Dε β 2 = ω 2 > 0. We can follow the sketch of the previous proof in the other cases. However, the result differs and it is impossible to construct a bump in those cases. More precisely, Case (N ):

aχ Dε -β 2 = -ω 2 < 0. The general solution of (2.20) is then      φ(x) = Acosh(ω(x -x)) + B sinh(ω(x -x)) - aK Dεω 2 , ρ(x) = Aχ ε cosh(ω(x -x)) + Bχ ε sinh(ω(x -x)) - β 2 K εω 2 .
(2.27)

For the transition between J and I at a point x, relations (2.21) and inequalities (2.15) remain unchanged. However, the expression of ρ in

I =]x, x + l[ becomes ρ(x) = Aχ ε (cosh(ω(x -x)) -1) + Bχ ε sinh(ω(x -x)) (2.28)
and thus never vanishes for x > x.

Case (Z): aχ Dε β 2 = 0. The general solution of (2.20) is then

     φ(x) = A + B(x -x) + aK 2Dε (x -x) 2 , ρ(x) = Aχ -K ε + Bχ ε (x -x) + Kβ 2 2ε (x -x) 2 .
(2.29)

Relations (2.21) become      Aχ -K ε = 0, A = P + + P -, B = β(P + -P -),
that is to say

   Aχ = K, A = P + + P -, B = β(P + -P -).
(2.30)

Then, we shall have ρ ≥ 0 in [x, x + l] only if B ≥ 0, and then

ρ(x) = Bχ ε (x -x) + Kβ 2 2ε (x -x) 2
remains always positive.

Remark 2.4. In the case when γ = 2, equation (2.19) is replaced by

χφ = ε γ ρ γ + K
for some constant K. Reporting in the second equation of system (S) gives, with the condition φ ≥ K/χ,

∂ 2 xx φ + a D γ(χφ -K) ε 1 γ -β 2 φ = 0,
that we may recast under the form of a Newton type equation

2∂ 2 xx φ + d dφ 2aε 2 Dχ(γ + 1) γ(χφ -K) ε γ+1 γ -β 2 φ 2 = 0.
We know from numerical simulations [START_REF] Natalini | A well-balanced numerical scheme for a one dimensional quasilinear hyperbolic model of chemotaxis[END_REF][START_REF] Natalini | A numerical comparison between degenerate parabolic and quasilinear hyperbolic models of cell movements under chemotaxis[END_REF] that stationary solutions with vacuum also exist in this case. However, exhibiting them analytically is much more complicated than in the case γ = 2 since no explicit computations are available.

Stationary problem (S) on R

In this section, we consider the case where a solution of system (S) is defined on the whole line R. We first construct a single bump as a stationary solution of the system in the Case (P), when aχ Dε β 2 = ω 2 > 0, then we will show that we cannot find any solution with a higher number of bumps.

3.1.

A single bump on R in Case (P). Let us begin with Case (P), for which aχ Dε β 2 = ω 2 > 0 and let us prove that we can construct a "one-bump" solution in that case. Proposition 3.1. Let M > 0. We assume that aχ Dε β 2 = ω 2 > 0. Among the solutions (ρ,φ) satisfying system (S) on R,

with ρ ∈ C 0 (R), φ ∈ C 2 (R) ∩ H 1 (R), ρ, φ ≥ 0 and such that R ρ(x)dx = M ,
there exists a unique "one bump" stationary solution, up to a translation, defined on R by the following expressions:

ρ(x) =        0 for x < x, M I χ ε (1 -cos(ω(x -x))) + ω β sin(ω(x -x)) for x ≤ x ≤ x + l, 0 for x > x + l, and
φ(x) =                M I ω 2 β 2 e β(x -x) , for x < x, M I -cos(ω(x -x)) + ω β sin(ω(x -x)) + aχ εβ 2 , for x ≤ x ≤ x + l, M I ω 2 β 2 e -β(x -x -l) , for x > x + l, with I = 2χ εω π -arctan( ω β ) + ω β > 0 and l = 2 ω π -arctan( ω β ) .
The energy of such solution, as defined at equation (2.4), is equal to

J ♭ (M ) = - εω 3 M 2 4β 2 (π -arctan( ω β ) + ω β ) . (3.1)
Proof. Following the results of Proposition 2.2, we can use the following expressions for ρ and φ, defined on R,

ρ(x) =      0, for x < x, - Aχ ε cos(ω(x -x)) + Bχ ε sin(ω(x -x)) + β 2 K εω 2 , for x ≤ x ≤ x + l, 0, for x > x + l, and
φ(x) =      P + e β(x-x) + P -e -β(x-x) , for x < x, -Acos(ω(x -x)) + B sin(ω(x -x)) + aK Dεω 2 , for x ≤ x ≤ x + l, (P ′ ) + e β(x-ȳ) + (P ′ ) -e -β(x-ȳ) , for x > x + l.
Since φ ∈ H 1 (R), we have φ(x) → ±∞ 0, from which we deduce that P -= 0 and that (P ′ ) + = 0 with (2.17). We are then able to express all the constants given by (2.17) and (2.21) as a function of A as follows: A can be rewritten in term of the mass M using equation (2.18), namely

K = χω 2 β 2 A, P + = ω 2 β 2 A, B = ω β A, (P ′ ) -= ω 2 β 2 e -β(ȳ-x-l) A,
A = M 2χ εω π -arctan( ω β ) + ω β .
Using equation (2.4), we are now able to compute the energy of a "onebump" solution as a function of the mass M and of the parameters of the system:

J ♭ (M ) = 1 2 R (ερ 2 -χφρ)dx = 1 2 x+l x ρ(ερ -χφ)dx = - K 2 x+l x ρdx = - χω 2 2β 2 AM = - εω 3 M 2 4β 2 (π -arctan( ω β ) + ω β ) . Remark 3.2. The bump is symmetric with respect to x + l 2 = x + 1 ω π -arctan( ω β )
. This is clear outside the region [x, x + l] and for x inside this region, the symmetry derives from the formula

-cos(ω(x -x)) + ω β sin(ω(x -x)) = 1 1 + ω 2 β 2 cos(π -arctan( ω β )) × cos(ω(x -x)) + sin(π -arctan( ω β ))sin(ω(x -x)) = 1 1 + ω 2 β 2 cos ω(x -x) -π + arctan( ω β
) .

We exhibit on Figure 1 the graphs of the functions ρ and φ defined in Proposition 3.1, for two different values of the parameters ω and β.

3.2.

Several bumps on R in Case (P). In the following proposition, we prove that we cannot have a solution defined by a finite or a countable number of bumps on the whole real line R. The contradiction will be obtained using the expressions for the constants (2.17) and (2.21) and the signs (2.15), in the case of a finite number of bumps, or using the finite value of the total mass, in the case of a countable number of bumps.

Proposition 3.3. Let M > 0. We assume that aχ Dε -β 2 = ω 2 > 0. We consider a solution (ρ,φ) satisfying system (S) on R, with ρ ∈ C 0 (R), φ ∈ C 2 (R), ρ, φ ≥ 0 and such that R ρ(x)dx = M .
There exists no solution defined on R by a finite number N ≥ 2 or by a countable number of bumps as described in Proposition 2.2.

Proof. In the current proof, we will use the notations displayed on Figure 2.

(i) Assume that we have a finite number of bumps, numbered by 1,2,••• ,N with N ≥ 2. Since φ ∈ H 1 (R), we have φ(x) → ±∞ 0, which leads in expression (2.14) to P - 1 = 0 and to P + 2 = 0, P - 2 = P + 1 e -β(ȳ-x-l) , using equations (2.17). We obtain now, from equation (2.21), the following expression for the constant for the second bump:

x x n x n + l n x n+1 x n+1 + l n+1 n-1 x n-1 + l n-1 d d n n+1 , B n , K n n A A , B n-1 n-1 , K n-1 A , B
B 2 ω = β(P + 2 -P - 2 ) = -βP + 1 e -β(ȳ-x-l)
. We find therefore a contradiction between this expression and the signs (2.15) of the constants, namely P + 1 > 0 and B 2 > 0. (ii) Assume now that we have a countable infinite number of bumps, indexed by n, n ∈ Z (or n ∈ N). We define by M n the mass of the n-th bump. Equation (2.17) yields that the product P + n P - n is constant. This constant cannot be zero unless one of the terms P + n or P - n vanishes. In that case, either P - n = 0 or P - n+1 = 0, since P + n = 0 implies P - n+1 = 0 from equation (2.17). As before, the same contradiction for B n or B n+1 with the signs (2.15) -ω 2 ≤ 0, we can prove that there is no stationary solution. In this section, we will show that, unlike the previous results, it is possible to construct multi-bumps solutions of system (S):

ερ∂ x ρ = χρ∂ x φ, D∂ 2 xx φ + aρ -Dβ 2 φ = 0, (S)
∂ x ρ(x) = 0, ∂ x φ(x) = 0, on x = 0,L. (4.1)
The number of bumps of these solutions depends on the length of the interval and on the parameters of the system. We will study successively different types of solutions, namely constant solutions in Subsection 4.1, single bump in Subsection 4.2, a half-bump at the boundary in Subsection 4.3, two bumps in Subsection 4.4, two half bumps at the boundaries in Subsection 4.5, one bump and a half in Subsection 4.6. We display some of these solutions in Figure 3.

For each type of solution, we will also compute the corresponding energy and we will compare those energies one to another. In Figure 4, we present the above mentioned stationary solutions on a bifurcation diagram by plotting their energy as a function of the length L of the domain. It is quite easy to see that we may construct a continuous family of solutions with an arbitrary number of bumps by putting several bumps side by side as soon as the length L of the interval is sufficiently large. We plot these curves using these formulas which will be proved in the following:

(ρ,φ) = M L , aM LDβ 2 ,
whatever the sign of aχ Dε β 2 . The energy of this solution is given by

J Cte (M,L) = εM 2 2Lβ 2 β 2 - aχ εD . (4.2)
In particular, in the case when aχ Dε β 2 > 0, this energy is negative and increases with L up to 0 as L → +∞.

From now on, we concentrate on the case when aχ Dε β 2 > 0.

4.2.

A single bump inside the interval [0,L]. In this section, we study the possibility of having a single bump solution to (S) inside the interval [0,L]. For that purpose, we will use the results of Proposition 2.2 and the notations are explained in Figure 5. We split the interval [0,L] into three

d 0 + l 1 + d 1 = L. (4.3)
We introduce two functions, which will be useful in the following: In particular, we will show later on that constraint (4.3) may be written under the simple form H(d 0 ) = L.

H : R + ∋ d → 2d + 2 ω π -arctan ω β tanh(βd) (4.
In the following proposition, the construction of such a bump is given and it is shown that its energy is smaller than the energy of the constant solutions computed at the previous subsection. The energy of such solution is equal to

J 1♭ (M,L) = - εω 3 M 2 4β 2 g • H -1 (L) , (4.6)
where H and g are defined at equations (4.4) and (4.5), and satisfies the following properties:

(i) The function L → J 1♭ (M,L) defined on [2π/ω,+∞[ increases from - εω 3 M 2 4πβ 2 to J ♭ (M ), defined at equation (3.1). (ii) For every L > 2π ω , J 1♭ (M,L) < J Cte (M,L). (iii) J 1♭ (M, 2π ω ) = J Cte (M, 2π ω ) and 0 < ∂J 1♭ ∂L (M, 2π ω ) < ∂J Cte ∂L (M, 2π ω ).
Proof. On the interval ]0,d 0 [, which corresponds to the left hand-side of the interval [0,L], the function ρ vanishes and the equation (2.14) gives the following expression for the function φ: φ(x) = P + 1 e β(x-d 0 ) + P - 1 e -β(x-d 0 ) . The boundary condition ∂ x φ(0) = 0 implies that

P + 1 e -βd 0 = P - 1 e βd 0 ( 4
     φ(x) = -A 1 cos(ω(x -x)) + B 1 sin(ω(x -x)) + aK 1 Dεω 2 , ρ(x) = - A 1 χ ε cos(ω(x -x)) + B 1 χ ε sin(ω(x -x)) + β 2 K 1 εω 2 .
Let us study now how the transitions at points d 0 and d 0 + l 1 link the constants A 1 , B 1 , K 1 , P + 1 , P - 1 , P + 2 , P - 2 defined above. We use the relations (2.21), corresponding to the transition of ρ,φ,∂ x φ at d 0 , and we obtain the three following equations:

           A 1 = β 2 χω 2 K 1 , A 1 ω 2 β 2 = P + 1 + P - 1 , B 1 ω = β(P + 1 -P - 1 ).
(4.9)

The transition at point d 0 + l 1 and the relation (2.17) give P + 2 = P - 1 e β(L-d 0 -l 1 ) = P - 1 e βd 1 and P - 2 = P + 1 e -β(L-d 0 -l 1 ) = P + 1 e -βd 1 (4.10) using equation (4.3). Therefore, relations (4.7), (4.8) and (4.10) give

P + 1 (1 -e 2β(d 0 -d 1 ) ) = 0. ( 4

.11)

If P + 1 = 0, then P - 1 = P + 2 = P - 2 = 0 and we get the trivial solution A 1 = B 1 = K 1 = 0. Therefore, equation (4.11) implies that

d 0 = d 1 .
Now, we consider the length of the bump given by (2.16). It is equal to Notice that the function H depends only on the parameters ω, β of the system and not on the length L of the interval and is continuously increasing from 2π/ω to +∞. Indeed,

l 1 = 2 ω π -arctan( B 1 A 1 ) = 2 ω π -arctan( ω β tanh(βd 0 )) .
H ′ (d) = 2tanh 2 (βd) 1 + ω 2 β 2 tanh 2 (βd) ω 2 β 2 + 1 > 0 (4.13)
for d > 0. The graph of function H is given in Figure 6. Finally, using formulas (2.18), (4.9) and (4.7), the mass of the bump is equal to

M = 2χβ 2 εω 3 P + 1 (1 + e -2βd 0 ) π -arctan( ω β tanh(βd 0 )) + ω β tanh(βd 0 ) = 2χβ 2 εω 3 P + 1 (1 + e -2βd 0 )g(d 0 ). ( 4

.14)

We can now conclude on the existence of a one-bump solution. For L < 2π ω , we cannot find d 0 ≥ 0 such that H(d 0 ) = L, hence there is no solution to (4.12) and no bump solution. For L ≥ 2π ω , there is exactly one solution to (4.12). With this value of d 0 , relation (4.14) gives the value of P + 1 (since M is fixed). Then (4.7) gives P - 1 , (4.10) gives P - 2 and P + 2 and (4.9) gives A 1 , B 1 and K 1 . Thus we have find exactly one bump solution. Its energy (using (2.4)) is given by

J 1♭ (M,L) = 1 2 d 0 +l 1 d 0 ρ(ερ -χφ)dx = - K 1 2 x+l 1 x ρdx = - χω 2 2β 2 A 1 M = - εω 3 M 2 4β 2 g(d 0 )
, where H(d 0 ) = L.

At last, let us prove some properties on the energy J 1♭ (M,L). Since

g ′ (d) = ω 3 β 2 (1 -tanh 2 (βd))tanh 2 (βd) 1 + ω 2 β 2 tanh 2 (βd) = ω 3 β 2 tanh 2 (βd) cosh 2 (βd)(1 + ω 2 β 2 tanh 2 (βd)) > 0, (4.15)
and H is increasing, it follows that

J 1♭ (M,L) = - εω 3 M 2 4β 2 g • H -1 (L)
is increasing with L. Furthermore, as L → +∞, we recover the situation of section 3.1, and in particular:

J 1♭ (M,L) = - εω 3 M 2 4β 2 g • H -1 (L) → J ♭ (M ) = - εω 3 M 2 4β 2 (π -arctan( ω β ) + ω β )
.

We may also compute the following expansions when d → 0,

g(d) = π -arctan ω β tanh(βd) + ω β tanh(βd) = π + (ωd) 3 3 + O(d 5 ),
and

H(d) - 2π ω = 2d - 2 ω arctan ω β tanh(βd) = 2 3 (ω 2 + β 2 )d 3 + O(d 5 ).
Thus, as L → 2π ω ,

g • H -1 (L) = π + ω 3 ω 2 + β 2 L - 2π ω + O (L - 2π ω ) 5/3 , which implies J 1♭ (M,L) = - εω 3 M 2 4β 2 π + εω 6 M 2 4β 2 π 2 (ω 2 + β 2 ) L - 2π ω + O (L - 2π ω ) 5/3 .
Observe that the expansion of J Cte (M,L) given in (4.2) near L = 2π ω is equal to

J Cte (M,L) = - εω 2 M 2 2Lβ 2 = - εω 3 M 2 4πβ 2 + εω 4 M 2 8π 2 β 2 L - 2π ω + O (L - 2π ω ) 2 , which implies that ∂J 1♭ ∂L (M, 2π ω ) < ∂J Cte ∂L (M, 2π ω ).
Finally, we may compare J Cte (M,L) and J 1♭ (M,L) easily. Indeed, on half of the interval. We call this type of solution "half bump". We also prove in the following proposition that its energy is smaller than the energy of the constant solutions described in subsection 4.1 and than the energy of the bump of subsection 4.2 of the same mass on an interval of the same length.

J 1♭ (M,L) < J Cte (M,L) is equivalent to 2g(d) < ωH(d),
Proposition 4.2. Let L > 0 and M > 0.There exists a stationary solution to system (S)-(4.1) with half bump on the boundary of the interval [0,L] if and only if L ≥ π ω , in which case it is unique up to the symmetry with respect to L/2. Its energy is equal to

J 1 2 ♭ (M,L) = 1 2 J 1♭ (2M,2L) = 2J 1♭ (M,2L) = - εω 3 M 2 2β 2 g • H -1 (2L) (4.16)
and satisfies the following properties:

(i) The function L → J 1 2 ♭ (M,L) defined on [π/ω,+∞[ increases from - εω 3 M 2 2πβ 2 to 2J ♭ (M ) = - εω 3 M 2 2β 2 (π -arctan( ω β ) + ω β )
, where J ♭ (M ) is defined at

equation (3.1). (ii) For every L > π ω , the inequality J 1 2 ♭ (M,L) < J Cte (M,L) holds. (iii) J 1 2 ♭ (M, π ω ) = J Cte (M, π ω ) and 0 < ∂J 1 2 ♭ ∂L (M, π ω ) < ∂J Cte ∂L (M, π ω ).
(iv) For every L ≥ 2π ω , J 1 2 ♭ (M,L) < J ♭ (M,L). Proof. The expression of the energy follows by considering a bump of mass 2M in the interval of length 2L and dividing the result by two. Hence, points (i), (ii) and (iii) come from (i), (ii) and (iii) in Proposition 4.1 and form the fact that in view of (4.2), J Cte (2M,2L) = 2J Cte (M,L). It remains to show (iv), and this is a direct consequence of the following lemma. We apply this lemma with ℓ = L ≥ 2π/ω and we deduce from the monotonicity that

g(H -1 (L)) > g(H -1 (2L)) 2 ,
which is exactly the desired inequality for (iv), using formulas (4.6) and (4.16).

Proof of Lemma 4.3. We denote ψ(θ) = θ -1 g(H -1 (θℓ)) and

d = d(θ) = H -1 (θℓ). Then, for θ > θ 0 , ψ ′ (θ) = g ′ (d(θ))d ′ (θ) θ - g(d(θ)) θ 2 , hence θ 2 ψ ′ (θ) ′ = θ g ′ (d(θ))d ′ (θ) ′ .
Since d ′ (θ) = ℓ/H ′ (d(θ)), and using the derivatives of H and g expressed at equations (4.13) and (4.15), we deduce

g ′ (d(θ))d ′ (θ) = ℓω 3 2(β 2 + ω 2 ) (1 -tanh 2 (βd(θ))),
which is a decreasing function of θ (d increases). Therefore, θ → θ 2 ψ ′ (θ) is decreasing. For θ = θ 0 , d = 0 thus, by definition of θ 0 ,

θ 2 0 ψ ′ (θ 0 ) = θ 0 g ′ (d(θ 0 ))d ′ (θ 0 ) -g(d(θ 0 )) = ℓθ 0 ω 3 2(β 2 + ω 2 ) -π = πω 2 β 2 + ω 2 -π = - πβ 2 β 2 + ω 2 < 0.
As a consequence, θ 2 ψ ′ (θ) is negative for θ > θ 0 , and the result follows. 4.4. Two bumps inside the interval. We investigate in this section the existence of a solution consisting of two bumps inside the interval [0,L] as shown in Figure 8. We shall see that a two bumps solution is necessarily the concatenation of two one bump solutions as defined at Proposition 4.1, set on two smaller intervals. Moreover, there exists a one parameter family of two bumps solutions, parametrized by d 0 (see Figure 8). Then the parameter d 2 and the masses M 1 and M 2 of the two bumps are uniquely defined. To define the masses, the following function is needed: where g is defined at equation (4.5). Since g is increasing, from equation (4.15), it is clear that G is the product of two positive increasing functions, hence is a positive increasing function.

G(d) = cosh(βd) π -arctan( ω β tanh(βd)) + ω β tanh(βd) = cosh(βd)g(d),
In the following proposition, we will also show that among all the configurations with two bumps inside the interval, the one corresponding to d 0 = 0 is the one with the smallest energy, whereas the symmetric one, satisfying d 0 = d 2 is the one with the highest energy. In this case (see Figure 8), the solution may be parametrized by the parameter d 0 with 0 ≤ d 0 ≤ d ♯ ≡ H -1 (L/2), corresponding to the first bump. The second bump is defined by the parameter d 2 which satisfies the relation L = H(d 0 ) + H(d 2 ) and the masses of the two bumps are defined by

(M 1 ,M 2 ) = M G(d 0 ) G(d 0 ) + G(d 2 ) , G(d 2 ) G(d 0 ) + G(d 2 )
. (4.17)

The symmetric configuration with respect to L/2 is also a stationary solution.

Its energy is then given by

J 2♭ (M,L,d 0 ) =J 1♭ (M 1 ,L 1 ) + J 1♭ (M 2 ,L 2 ) = - εω 3 M 2 1 4β 2 g(d 0 ) - εω 3 M 2 2 4β 2 g(d 2 ) , (4.18)
and satisfies the following properties:

(i) The function d → J 2♭ (M,L,d) defined on [0,d ♯ ] is increasing.
(ii) For L > 4π ω , the following inequalities hold:

J Cte (M,L) > J 2♭ (M,L,d ♯ ) = 2J 1♭ (M/2,L/2) = 1 2 J 1♭ (M,L/2) > J 1♭ (M,L). (iii) J 2♭ (M, 4π ω ,0) = J Cte (M, 4π ω ) and 0 < ∂J 2♭ ∂L (M, 4π ω ,0) < ∂J Cte ∂L (M, 4π ω ).
(iv) As L → +∞, J 2♭ (M,L,0) → J ♭ (M ), defined at equation (3.1), and

J 2♭ (M,L,d ♯ ) → 1 2 J ♭ (M ) = 1 2 J 1♭ (M,∞).
Remark 4.5. In the case L = 4π ω , the continuum of solutions degenerates into a unique 2 bumps solution.

Proof. We use the notations of Figure 8. We begin to prove that a twobumps solution is necessarily the concatenation of two one bump solutions as defined at Proposition 4.1, set on two smaller intervals. The length of the total interval is

L = d 0 + l 1 + d 1 + l 2 + d 2 (4.19)
and we look for a relation linking parameters d 0 , d 1 and d 2 . Inspired by equations (4.7) and (4.8), the boundary conditions give

P - 1 = P + 1 e -2βd 0 , P - 3 = P + 3 . (4.20)
and the transition conditions (2.17) imply Let us determine all the constraints so that when we concatenate the two bumps, we obtain a solution in the interval [0,L]. First, note that since L j ≥ 2π/ω (j = 1, 2), this already imposes L = L 1 + L 2 ≥ 4π/ω. Second, the total mass M must be equal to

P + 2 = P - 1 e βd 1 ,P - 2 = P + 1 e -βd 1 ,P + 3 = P - 2 e βd 2 ,P - 3 = P + 2 e -βd 2 . ( 4
M 1 + M 2 . Once (L 1 ,L 2 ,M 1 ,M 2 ) are given such that L = L 1 + L 2 and M = M 1 + M 2 ,
the two bumps are completely determined, which means that, for the moment, we have two free parameters. The relation (4.12) holds for the two bumps, hence

H(d 0 ) = L 1 and H(d 2 ) = L 2 . (4.23)
Furthermore, since (ρ,φ) solves the system (S) in ]0,L 1 [ and in ]L 1 ,L[, (ρ,φ) is a solution near L 1 if and only if it is continuous. Indeed, their derivatives is zero by the Neumann condition (4.1). Since ρ is zero near L 1 on the right and on the left, we are reduced to verify the relation φ(0) = φ(L). Indeed, φ(0) (resp. φ(L)) is the value of φ at the left (resp. right) of L 1 due to the symmetry of the bumps, see Remark 3.2. Using equations (2.14) and (4.7), we have φ(0) = P + 1 e -βd 0 + P - 1 e +βd 0 = 2P + 1 e -βd 0 , and we express consequently the mass M 1 as a function of φ(0) thanks to equation (4.14):

M 1 = 2χβ 2 εω 3 P + 1 (1 + e -2βd 0 )g(d 0 ) = φ(0)
2χβ 2 εω 3 G(d 0 ). We have then the relations

M 1 = φ(0) 2χβ 2 εω 3 G(d 0 ) and M 2 = φ(L) 2χβ 2 εω 3 G(d 2 ), arguing symmetrically in [L 1 ,L]. Therefore, the condition φ(0) = φ(L) is equivalent to M 1 G(d 0 ) = M 2 G(d 2 ) (4.24)
and we have now only one parameter free.

Let us determine to which interval the parameters d 0 and d 2 (and the masses M 1 and M 2 ) belong. Since L > 4π/ω, we may define d > 0 as the solution to H( d) = L -2π/ω, which corresponds to the value of d 2 (resp. d 0 ) whenever d 0 = 0 (resp. d 2 = 0). Then, using equations (4.23), the relation

L = H(d 0 ) + H(d 2 ) ≥ 2π/ω + H(d 0 ) imposes d 0 ≤ d, and also d 2 ≤ d by symmetry. Notice that M 1 = M G(d 0 ) G(d 0 ) + G(d 2 ) = M 1 + G(d 2 )/G(d 0 ) is then a decreasing function of d 0 , which varies between M 1 + G(0)/G( d) = M 1 + π/G( d)
and

M 1 + G( d)/π
. Therefore, if we want to fix the masses M 1 and M 2 , with

M 1 + M 2 = M , they have to satisfy M 1 ,M 2 ∈ M 1 + G( d) π , M 1 + π G( d) , or, equivalently, M 1 M 2 ∈ π G( d) , G ( d) π .
We notice that the masses M 1 and M 2 are bounded away from 0 by a constant depending on L.

Then, recalling that we started from two one bump solution of mass M 1 (resp. M 2 ) set on an interval of length L 1 (resp. L 2 ), the two parameters

(d 0 ,d 2 ) have to solve M 1 G(d 0 ) = M 2 G(d 2 ) and L = L 1 + L 2 = H(d 0 ) + H(d 2 )
by (4.23) and (4.24). This system has a unique solution: if we let d 0 increase, then d 2 decreases by the second relation, since H is an increasing function; on the contrary, for M 1 and M 2 given, if we let d 0 increase, then d 2 increases by the first relation, since G is a decreasing function. Hence we find a unique solution d 2 as a function of d 0 .

As a consequence, we have obtained a continuous one parameter family of solutions with two bumps, that we may parametrize by

L 1 ∈ [2π/ω,L -2π/ω[ or, equivalently, by d 0 ∈ [0, d[. The parameter d 2 ∈ [0, d[ is given by the relation L = H(d 0 ) + H(d 2 )
. Moreover, the masses (M 1 ,M 2 ) must satisfy the relation M = M 1 + M 2 and the equation (4.24), which give exactly one solution computed at equation (4.17). Now, let us consider the energy of the two bumps solutions we constructed. One may wonder what are the variations of the energy when L 1 or, equivalently, d 0 varies and if a single bump solution has less energy than these solutions with two bumps. For the solution with two bumps, using the definition of the energy of a single bump (4.6), the energy is given by formula (4.18), that is to say Lemma 4.6. Let L ≥ 4π/ω and M > 0 be given. Then, the function

J 2♭ (M,L,d 0 ) =J 1♭ (M 1 ,L 1 ) + J 1♭ (M 2 ,L 2 ) = - εω 3 M 2 1 4β 2 g(d 0 ) - εω 3 M 2 2 4β 2 g(d 2 ) ,
J 2♭ (M,L,•) is increasing on [0,d ♯ ]. Proof. Since H ′ (d 0 ) + H ′ (d 2 )∂ d 0 d 2 = 0, we compute -H ′ (d 2 ) 4β 2 εω 3 dJ 2♭ dd 0 = H ′ (d 2 ) ∂ ∂d 0 M 2 1 g(d 0 ) + M 2 2 g(d 2 ) -H ′ (d 0 ) ∂ ∂d 2 M 2 1 g(d 0 ) + M 2 2 g(d 2 ) = H ′ (d 2 ) 2M 1 ∂ d 0 M 1 g(d 0 ) - M 2 1 g ′ (d 0 ) g 2 (d 0 ) + 2M 2 ∂ d 0 M 2 g(d 2 ) -H ′ (d 0 ) 2M 2 ∂ d 2 M 2 g(d 2 ) - M 2 2 g ′ (d 2 ) g 2 (d 2 ) + 2M 1 ∂ d 2 M 1 g(d 0
) .

We now report, from (4.17),

∂M 1 ∂d 0 = - ∂M 2 ∂d 0 = M G(d 2 )G ′ (d 0 ) (G(d 0 ) + G(d 2 )) 2 and ∂M 1 ∂d 2 = - ∂M 2 ∂d 2 = -M G ′ (d 2 )G(d 0 ) (G(d 0 ) + G(d 2 )) 2 to obtain - 4β 2 H ′ (d 2 ) εω 3 M 2 dJ 2♭ dd 0 = H ′ (d 2 ) (G(d 0 ) + G(d 2 )) 2 2G(d 2 )G ′ (d 0 ) G(d 0 ) + G(d 2 ) G(d 0 ) g(d 0 ) - G(d 2 ) g(d 2 ) -g ′ (d 0 ) G 2 (d 0 ) g 2 (d 0 ) - H ′ (d 0 ) (G(d 0 ) + G(d 2 )) 2 2G(d 0 )G ′ (d 2 ) G(d 0 ) + G(d 2 ) G(d 2 ) g(d 2 ) - G(d 0 ) g(d 0 ) -g ′ (d 2 ) G 2 (d 2 ) g 2 (d 2 )
.

We now use that G(d) = cosh(βd)g(d) to infer

-(G(d 0 ) + G(d 2 )) 2 4β 2 H ′ (d 2 ) εω 3 M 2 dJ 2♭ dd 0 = 2 cosh(βd 0 ) -cosh(βd 2 ) G(d 0 ) + G(d 2 ) G(d 2 )G ′ (d 0 )H ′ (d 2 ) + G(d 0 )G ′ (d 2 )H ′ (d 0 ) + H ′ (d 0 )g ′ (d 2 )cosh 2 (βd 2 ) -H ′ (d 2 )g ′ (d 0 )cosh 2 (βd 0 ). (4.25)
Using (4.15) and (4.13), we then deduce

H ′ (d 0 )g ′ (d 2 )cosh 2 (βd 2 ) = 2 ω 3 β 2 ω 2 β 2 + 1 tanh 2 (βd 0 )tanh 2 (βd 2 ) (1 + ω 2 β 2 tanh 2 (βd 0 ))(1 + ω 2 β 2 tanh 2 (βd 2 )) =H ′ (d 2 )g ′ (d 0 )cosh 2 (βd 0 ), so that the last line in (4.25) is zero. If d 0 < d 2 , that is d 0 < d ♯ ,
then the right-hand side of (4.25) is < 0 since G, G ′ and H ′ are positive, which is the desired result. Now, let us prove the various properties of the energy described in Proposition 4.4. For (ii), the first and last inequalities are a consequence of Lemma 4.3. (iv) follows directly from either the fact that we have symmetrized one bump with parameters (M/2,L/2), either the fact that (M 1 ,M 2 ) → (0,M ) and d 2 → +∞. We now compute the derivative for (iii). Observe first that when d 0 = 0, d 2 = d = H -1 (L -2π/ω), thus d → 0 as L → 4π/ω. We have already seen in section 4.2 that

g(d) = π + (ωd) 3 3 + O(d 5 ), thus G(d) = g(d)cosh(βd) = π + πβ 2 2 d 2 + (ωd) 3 3 + O(d 4 ).
Therefore,

4β 2 εω 3 M 2 J 2♭ (M,L,0) = - π + G 2 ( d) g( d) (π + G( d)) 2 = - 1 2π + ω 3 d3 12π 2 + O( d4 ). Since L - 4π ω = H( d) - 2π ω = 2 3 (ω 2 + β 2 ) d3 + O( d5 ), it follows that 4β 2 εω 3 M 2 J 2♭ (M,L,0) = - 1 2π + ω 3 8π 2 (ω 2 + β 2 ) L - 4π ω + O (L - 4π ω ) 4 .
As a consequence,

0 < ∂J 2♭ ∂L (M,4π/ω) = εω 6 M 2 32π 2 β 2 (ω 2 + β 2 ) < εω 4 M 2 32π 2 β 2 = ∂J Cte ∂L (M,4π/ω),
from the explicit expression of J Cte in (4.2). This completes the proof.

For (iv), notice that it is natural to recover the energy of a single bump when d 0 = 0 and L is large since the mass of the bump on the left tends 0 L 2L to zero. Despite our efforts, we have not been able to determine the variations of the least energy configuration L → J 2♭ (M,L,0), even though numerically, it seems convincing that it is increasing on a small interval and then decreasing. In this case (see Figure 8 for the notations), the solutions may be parametrized by the parameter 0 ≤ d 0 ≤ d ♯ ≡ H -1 (L) for the first half bump. The second half bump is defined by the parameter d 2 which satisfies the relation 2L = H(d 0 ) + H(d 2 ) and the masses of the two half bumps are given by equation (4.17). The symmetric configuration with respect to L/2 is also a solution.

Its energy is given by

J 2× 1 2 ♭ (M,L,d 0 ) = 1 2 J 2♭ (2M,2L,d 0 ) = 2J 2♭ (M,2L,d 0 )
and satisfies the following properties: (iii)

(i) The function J 2× 1 2 ♭ (M,L,•), defined on [0,d ♯ ], is increasing. (ii) J Cte (M,L) > J 2× 1 2 ♭ (M,L,d ♯ ) = J 1♭ (M,L).
J 2× 1 2 ♭ (M, 2π ω ,0) = J Cte (M, 2π ω ) and 0 < ∂J 2× 1 2 ♭ ∂L (M, 2π ω ) < ∂J Cte ∂L (M, 2π ω ).
(iv)

As L → +∞, J 2× 1 2 ♭ (M,L,d ♯ ) → J ♭ (M ) = J 1♭ (M,∞) and J 2× 1
2 ♭ (M,L,0) → 2J ♭ (M ). 4.6. One bump inside the interval and half bump on the boundary. Finally, we consider the solutions where we have one bump inside the interval and half bump on the boundary. Clearly, we may use a reflection across the boundary for the half bump to reduce to the case of two bumps, see Figure 10. We have consequently the following proposition: Proposition 4.8. Let L > 0 and M > 0. There exists a continuum of solutions to system (S)-(4.1) with one bump inside the interval [0,L] and one half bump on the boundary of [0,L] if and only if L ≥ 3π ω .

In this case, the solutions may be parametrized by 0 ≤ d 0 ≤ d ♯ ≡ H -1 (Lπ/ω) for the bump on the center of the interval or by 0 ≤ d 2 ≤ d ≡ H -1 (2L -4π/ω) for the half bump on the boundary, the other parameter being given by relation

L = H(d 0 ) + 1 2 H(d 2 )
. The mass M 1 of the bump and the mass M 2 of the half bump are defined by

(M 1 ,M 2 ) = M 2G(d 0 ) 2G(d 0 ) + G(d 2 ) , G(d 2 ) 2G(d 0 ) + G(d 2 ) . ( 4

.26)

The symmetric configuration with respect to L/2 is also a solution. Its energy J (1+ 1 2 )♭ (M,L,d 0 ) is given by

J (1+ 1 2 )♭ (M,L,d 0 ) = J 1♭ (M 1 ,L 1 ) + J 1 2 ♭ (M 2 ,L 2 ) = - εω 3 M 2 1 4β 2 g(d 0 ) - εω 3 M 2 2 2β 2 g(d 2 ) , (4.27)
and satisfies the following properties:

(i) The function J (1+ 1 2 )♭ (M,L,•) is increasing on [0,d ⋆ ] and decreasing on [d ⋆ ,d ♯ ], where d ⋆ ≡ H -1 (2L/3) is the value corresponding to d 0 = d 2 . (ii) J Cte (M,L) > J (1+ 1 2 )♭ (M,L,d ⋆ ) = 3J 1 2 ♭ ( M 3 , L 3 ) = 3 2 J 1♭ ( 2M 3 , 2L 3 ) = 2 3 J 1♭ (M, 2L 3 ) > J 1♭ (M,L). (iii)J (1+ 1 2 )♭ (M, 3π ω ,0) = J Cte (M, 3π ω ) and 0 < ∂J (1+ 1 2 )♭ ∂L (M, 3π ω ) < ∂J Cte ∂L (M, 3π ω ).
(iv)

As L → +∞, J (1+ 1 2 )♭ (M,L,0) → 2J ♭ (M ) = J 1 2 ♭ (M,∞), J (1+ 1 2 )♭ (M,L,d ⋆ ) → 2 3 J ♭ (M ) = 2 3 J 1♭ (M,∞) and J (1+ 1 2 )♭ (M,L,d ♯ ) → J ♭ (M ). Proof.
The argument is very close to the case of two bumps inside the interval explained in details in section 4.4, hence we only sketch the proof. The constraint on the length is

L = H(d 0 ) + 1 2 H(d 2 ),
and this gives the upper bounds d ♯ and d for d 0 and d 2 respectively. If M 1 (resp. M 2 ) denotes the mass of the bump (resp. half bump), then we must have

M = M 1 + M 2 and M 1 G(d 0 ) = 2M 2 G(d 2 ) ,
the second one being the matching condition for φ, thus relation (4.17) is replaced now by relation (4.26). Moreover, the energy of the one bump and one half bump solution is given by equation (4.27). To study the variations of J (1+ 1 2 )♭ (M,L,•), we perform some computations which are very similar to those of Lemma 4.6 and we obtain

β 2 H ′ (d 2 ) εω 3 M 2 dJ (1+ 1 2 )♭ dd 0 = -2 cosh(βd 0 ) -cosh(βd 2 ) (2G(d 0 ) + G(d 2 )) 3 G(d 2 )G ′ (d 0 )H ′ (d 2 ) + 2G(d 0 )G ′ (d 2 )H ′ (d 0 ) .
Since the bracket is a sum of positive terms, dJ (1+ 1

2 )♭ dd 0 has the same sign as As a consequence, except in the exceptional case ωL ∈ πN, we must have A = 0. If ωL ∈ πN, then there exists a continuous family of solutions

     φ(x) = -Acos(ωx) + aK Dεω 2 , ρ(x) = - Aχ ε cos(ωx) + β 2 K εω 2 .
We have represented some typical solutions ρ when (L,ω) = (π,1) and when (L,ω) = (2π,1) in Figure 11. Note that the mass of such a solution is given by

M = L 0 ρ(x) dx = β 2 KL εω 2 since ωL ∈ πN.
Taking into account the conditions ρ ≥ 0 and φ ≥ 0, we then obtain the fuzzy solutions in [0,L]

     φ(x) = -Acos(ωx) + εM χL , ρ(x) = - Aχ ε cos(ωx) + M L ,
with A < εM χL . All these solutions have the same energy

J Fuzzy (M,L) = - εω 2 M 2 2Lβ 2
, which is the same as the energy J Cte of the constant solution for this value of L. It is clear that if we slightly increase L, then all these solutions disappear, except the two extreme ones which vanish either at x = 0 either at x = π (the red and the cyan solutions), which are symmetric with respect to L/2 and will give the half bump solutions. The same phenomenon occurs when ωL = 2π, leading to the solutions with one bump inside and two half bumps on the boundary. Remark 4.9. In the Case (N ) and Case (Z), that is to say if aχ Dε β 2 = -ω 2 ≤ 0, we can prove that either ρ = 0 everywhere, either ρ > 0 everywhere. Therefore, the only nonzero solution is given by the constant state exhibited in Section 4.1.

Numerical simulations on a bounded interval [0,L]

Let us notice that finding analytically the stability of the previous exhibited stationary solutions is a hard task and out of scope for the moment. An important difficulty is that the stationary solutions lie on the boundary of the admissible set of functions, since the nonnegativity constraint on the density is saturated.

Since we prove in section 4 that the half bump solution has the smallest energy among the stationary solutions we consider, in the case when L ≥ π ω , we expect this solution to be a stable one. However, the stability of the other stationary solutions is not that clear and seems to depend deeply on the parameters of the system. In what follows, we present a numerical study of the stability of these solutions for three different sets of parameters :

•

Case 1: ε = 2, γ = 2, χ = 10, D = 0.1, a = 20, b = 10, α = 1 on a domain of length L = 3. • Case 2: ε = 2, γ = 2, χ = 50, D = 1, a = 1, b = 1, α = 1 on a domain of length L = 2. • Case 3: ε = 2, γ = 2, χ = 50, D = 1, a = 4, b = 1, α = 1 on a domain of length L = 2.
For one set of parameters (Case 2), the half bump seems to be indeed the only stable solution. However, for another set of parameters (Case 1), a wide range of stationary solutions look stable in our numerical experiments.

Generalizing the conditions found in Propositions 4.1 to 4.8, we may say that we can construct (at least) one k half bumps solution under the following condition linking the length of the domain and the parameters of the system: L ≥ kπ ω . In this formula, we assume that one bump accounts for two half-bumps. Hence, for the three previous sets of parameters and given length, we may construct the following types of stationary solutions (except from the constant solution) :

• Case 1: 1 ≤ k ≤ 28, that is to say we may construct up to 14 bumps solution (or equivalently 13 bumps inside the interval and 2 half bumps on the boundaries). • Case 2: 1 ≤ k ≤ 3, that is to say we may construct up to 1 bump inside the interval and 1 half bump on the boundary. • Case 3: 1 ≤ k ≤ 6, that is to say we may construct up to 3 bumps inside the interval (or equivalently 2 bumps inside the interval and 2 half bumps on the boundaries). In the following, we use an adapted numerical scheme based on a upwinding technique described in [START_REF] Natalini | A well-balanced numerical scheme for a one dimensional quasilinear hyperbolic model of chemotaxis[END_REF] and improved in [START_REF] Natalini | A numerical comparison between degenerate parabolic and quasilinear hyperbolic models of cell movements under chemotaxis[END_REF]. This scheme is constructed in order to deal with vacuum, to preserve the positivity of the solution and to be well-balanced, that is to say to preserve stationary solutions with constant velocities. The initial condition is a perturbation of the half-bump stationary solution and this stationary solution seems stable. 5.1. Case 1: stability of 1 half bump, 1 bump, 2 half bumps and 2 bumps. We begin with the first set of parameters (Case 1) and we study the stability of the half-bump solution on the boundary described at Proposition 4.2. In Figure 12, we show the results of a numerical simulation of system (1.3) with an initial condition equal to a perturbation of a half bump on a boundary. Note that the perturbation is chosen so that the initial condition is a nonnegative function. It is clear that the half bump stationary solution is stable under such a perturbation. Now, we consider the stability of various two half-bumps solutions. We compute here the asymptotic solutions of system (1.3) starting with three different two half bumps solution, namely the symmetric one with the maximal energy (on top on the left), the one with the minimal energy (on top on the right ) and a third one with an intermediate energy (on bottom). In Figure 13, we can see that all these stationary solutions are numerically stable, although we could expect only the one with minimal energy to be stable. See Proposition 4.7 for more details.

We continue with the study of the stability of the one bump stationary solution, which existence and unicity have been shown at Proposition 4.1. A first test, not presented here, shows that this solution remains stable under a slight perturbation. In Figure 14, we show two different tests: we take an initial condition equal to a translation of the one bump solution. If the translation is small enough (figure on the left), the asymptotic solution is still the one bump stationary solution, whereas if the translation is large enough (figure on the right), the asymptotic solution is one bump with a half bump on the boundary.

Finally, we study in Figure 15 the stability of the two-bumps solution. As demonstrated in Proposition 4.4, there is a continuum of such solutions. We first take as initial datum a perturbation of three different two-bumps solutions and the results presented in Figure 15 show that all these solutions are stable. More precisely, we have considered the two-bumps solution with a minimal energy (on top, on the left), symmetric with a maximal energy 3). On each subfigure, the initial condition is a perturbation of one two half-bumps stationary solution. All these stationary solutions seem stable. (on top, on the right) and with an intermediate energy (on bottom, on the left). Finally, if the initial datum of the simulation is a translation of the two bumps solution with an intermediate energy, we find asymptotically a solution with two bumps and one half-bump on the boundary.

5.2.

Case 2: stability of 1 bump. For this second set of parameters, only a few types of stationary solutions exist: the half-bump solution, the bump solution, the two half bumps solution and a solution made of one central bump and one lateral half-bump. In Figure 16, we display the 3). On top-left, top-right and bottomleft subfigures, the initial condition is a perturbation of one two bumps stationary solution. All these stationary solutions are stable. On the last subfigure (bottom-right), the initial condition is a translation of a two-bumps solution and we reach asymptotically a solution composed of 2 bumps and one lateral half-bump.

asymptotic solutions for two different initial conditions, that is to say the one bump stationary solution translated (on the left) or a two-bumps solution computed with another set of parameters and translated (on the right). For this set of parameters, the asymptotic solution we obtain is always a half bump on the boundary of the interval (on the left or on the right). From these simulations, it seems that the half-bump solution is the only stable stationary solution for the set of parameters of Case 2, unlike the results of Case 1.

5.3. Case 3: stability of 1 bump and 2 bumps. Finally, we present in Figure 17 the results of a stability study for stationary solutions in the case of an intermediate set of parameters where 2 bumps solutions exist. Beginning with a slight perturbation of the one bump solution (figure on the left), the asymptotic solution of system (1.3) is a half bump on the boundary of the interval. Therefore, the 1 bump solution seems to be unstable. However, the half-bump solution is not the only stable solution of the system as before, since the second test (figure on the right) shows that, starting from a two bumps solution, we find two half-bumps as an asymptotic solution of the system. The results for the other two-bumps stationary solutions are the same. numerically unstable (and a half-bump seems stable); on the right, a two bumps solution is numerically unstable (and two half-bumps seem stable).

  and, from equation (2.16), l = 2 ω πarctan( ω β ) . The remaining constant

Figure 1 .

 1 Figure 1. A one-bump solution defined by the function ρ (in red) and the function φ (in green) for two different values of the parameters (ω,β), namely (1,0.5) on the left and (1,10) on the right.

Figure 2 .

 2 Figure 2. Notations for the proof of Proposition 3.3.

4 .

 4 Stationary problem (S) on a bounded interval [0,L] ⊂ R : a bifurcation diagram

Figure 3 .

 3 Figure 3. Various stationary solutions: constant (in green), one lateral half bump (in red), one central bump (in blue), two bumps with minimal energy (in black), two symmetric bumps with maximal energy (in cyan). The parameters are the following: κ = 1, χ = 10, D = 0.1, a = 20, b = 10, L = 1 and M = 10.

•

  equation (4.2) for the energy of the constant solution, • equation (4.16) for the energy of the half bump solution, • equation (4.6) for the energy of the 1 bump solution, • expression of Proposition 4.7 for the energy of the 2 half bumps solution, • equation (4.18) for the energy of the 2 bumps solution.

Figure 4

 4 Figure 4 is therefore the summary of al the computed solutions in the next subsections.

Figure 4 . 4 . 1 .

 441 Figure 4. Bifurcation diagram of the energies of different types of stationary solutions (constant, half bump, 2 half bumps, 1 bump, 2 bumps...) as a function of the length L of the domain

Figure 5 .

 5 Figure 5. Notations used in Subsection 4.2 for the case of a single bump inside the interval [0,L]

Proposition 4 . 1 .

 41 Let L > 0 and M > 0. There exists a stationary solution of system (S) with one single bump inside the interval [0,L], satisfying (4.1)if and only if L ≥ 2π ω, in which case it is unique and symmetric with respect to L/2.

. 7 ) 2 ( 4

 724 and the second boundary condition ∂ x ρ(0) = 0 is satisfied since ρ = 0 on ]0,d 0 [. We now consider the interval ]d 0 + l 1 ,L[, which corresponds to the right hand-side of the interval [0,L]. We obtain the expression φ(x) = P + 2 e β(x-L) + P - 2 e -β(x-L) from equation (2.14). In the same way, the boundary condition ∂ x φ(L) = 0 leads to P + 2 = P - .8) and the second boundary condition∂ x ρ(L) = 0 is satisfied since ρ = 0 on ]d 0 + l 1 ,L[.Finally, on the center part ]d 0 ,d 0 + l 1 [ of the interval [0,L] we can write, using equation(2.13)

Figure 6 .

 6 Figure 6. Graph of function H

Figure 7 . 4 . 3 .

 743 Figure 7. Notations for the case of half bump on the boundary of the interval [0,L]

Lemma 4 . 3 .

 43 Let ℓ ≥ 2π/ω be given and let us define θ 0 ≡ 2π/(ℓω) ∈]0,1]. Then, the functionθ → g(H -1 (θℓ)) θ defined on ]θ 0 ,+∞[ is decreasing.

Figure 8 .

 8 Figure 8. Notations for the case of two bumps inside the interval [0,L]

Proposition 4 . 4 .

 44 Let L > 0 and M > 0. There exists a continuum of stationary solutions to system (S)-(4.1) with two bumps inside the interval [0,L] if and only if L ≥ 4π ω .

. 21 )

 21 Combining the relations (4.20) and (4.21) (and since we assume P 1 = 0), we getd 1 = d 2 + d 0 . (4.22) Let us now consider the point y ≡ 2d 0 + l 1 = L -(2d 2 + l 2 ) ∈ [0,L] (in view of (4.19) and (4.22)). The solution on [0,y] satisfies the Neumann boundary condition (4.1) on the boundary x = 0. The bump being symmetric, the derivatives of ρ and φ vanish at y. In particular, (ρ,φ) |[0,y] is actually a single bump solution in the interval [0,y]. Similarly, (ρ,φ) |[y,L] is also a single bump solution in the interval [y,L]. Conversely, let us decompose L = L 1 + L 2 , and assume that we have two single bumps: one in the interval [0,L 1 ], and the other one in the interval [L 1 ,L].

with d 2

 2 defined by L = H(d 0 ) + H(d 2 ) and M 1 ,M 2 by (4.17). The function J 2♭ is seen as a function of the single variable d 0 ∈ [0, d[. Moreover, since we may assume d 0 ≤ d 2 by symmetry, it is natural to define d ♯ to be the solution of L = 2H(d ♯ ) and to work only for d 0 ∈ [0,d ♯ ]. Observe that 2H(d ♯ ) = L = H( d) + H(0) ≤ 2H( d), thus, H being increasing, d ♯ ≤ d. Let us prove now the following lemma, which corresponds to point (i) of Proposition 4.4, that is to say the energy for a two bumps solution increases with d 0 .

Figure 9 .

 9 Figure 9. Technique to construct a solution with two half bumps on the boundaries of the interval [0,L]

4. 5 .

 5 Two half bumps on the two boundaries. Now we construct some solutions made of two half bumps on the two boundaries. For that purpose, we follow what we made in Section 4.3: to find a two half bumps solution of mass M on a domain of length L, we construct a two bumps solution of mass 2M on a domain of length 2L as explained in Section 4.4 and we keep half of the domain in order to obtain two half bumps on the boundaries. The technique is illustrated in Figure9. The results are straightforward following Proposition 4.4 and are summarized in the following proposition. Proposition 4.7. Let L > 0 and M > 0. There exists a continuum of stationary solutions to system (S)-(4.1) with two half bumps on the boundaries of the interval [0,L] if and only if L ≥ 2π ω .

2 Figure 10 .

 210 Figure 10. Notations for the case of one bump inside and half bump on the boundary of the interval [0,L]

d 2 - 2 = d 0 . 4 . 7 .

 2247 d 0 , and (i) follows since d ⋆ is the value of d 0 corresponding to the case d Fuzzy solutions. In this section, we consider some very particular solutions, such that ρ > 0 on the interval [0,L]. Then, from (2.13) with x = 0 and the Neumann conditions (4Asin(ωL) = 0.

Figure 11 .

 11 Figure 11. Graphs of some solutions ρ when (L,ω) = (π,1) (on the left) and when (L,ω) = (2π,1) (on the right).

Figure 12 .

 12 Figure 12. Case 1 -Initial condition and asymptotic solution for system(1.3).The initial condition is a perturbation of the half-bump stationary solution and this stationary solution seems stable.

Figure 13 .

 13 Figure[START_REF] Gu | Well-posedness of 1-D compressible Euler-Poisson equations with physical vacuum[END_REF]. Case 1 -Initial conditions and asymptotic solutions for system (1.3). On each subfigure, the initial condition is a perturbation of one two half-bumps stationary solution. All these stationary solutions seem stable.

Figure 14 .

 14 Figure 14. Case 1 -Initial conditions and asymptotic solutions for system (1.3). On each subfigure, the initial condition is a translation of the 1 bump-solution: a small translation (on the left) and a large translation (on the right).

Figure 15 .

 15 Figure[START_REF] Horstmann | From 1970 until present: the Keller-Segel model in chemotaxis and its consequences[END_REF]. Case 1 -Initial conditions and asymptotic solutions for system (1.3). On top-left, top-right and bottomleft subfigures, the initial condition is a perturbation of one two bumps stationary solution. All these stationary solutions are stable. On the last subfigure (bottom-right), the initial condition is a translation of a two-bumps solution and we reach asymptotically a solution composed of 2 bumps and one lateral half-bump.

Figure 16 .Figure 17 .

 1617 Figure[START_REF] Jang | Well and ill-posedness for compressible Euler equations with vacuum[END_REF]. Case 2 -On the left, the one bump stationary solution is unstable; on the right, the initial condition is a two bumps solution computed with the set of parameters of Case 1. In both cases, the asymptotic solution is a lateral half-bump, on the left or on the right of the interval.

Acknowledgement. This work has been partially supported by the ANR project MONUMENTALG, ANR-10-JCJC 0103.