R Natalini

M Ribot

M Twarogowska

A WELL-BALANCED NUMERICAL SCHEME FOR A ONE DIMENSIONAL QUASILINEAR HYPERBOLIC MODEL OF CHEMOTAXIS

Keywords: AMS Primary: 65M08; Secondary: 35L60, 92B05, 92C17

We introduce a numerical scheme to approximate a quasi-linear hyperbolic system which models the movement of cells under the influence of chemotaxis. Since we expect to find solutions which contain vacuum parts, we propose an upwinding scheme which handles properly the presence of vacuum and, besides, which gives a good approximation of the time asymptotic states of the system. For this scheme we prove some basic analytical properties and study its stability near some of the steady states of the system. Finally, we present some numerical simulations which show the dependence of the asymptotic behavior of the solutions upon the parameters of the system.

Introduction

The movement of bacteria, cells or other microorganisms under the effect of a chemical stimulus, represented by a chemoattractant, has been widely studied in mathematics in the last two decades, see [START_REF] Horstmann | From 1970 until present: the Keller-Segel model in chemotaxis and its consequences[END_REF][START_REF] Murray | Mathematical biology. I[END_REF][START_REF] Murray | Mathematical biology[END_REF][START_REF] Perthame | Transport equations in biology[END_REF], and numerous models involving partial differential equations have been proposed. The basic unknowns in these chemotactic models are the density of individuals and the concentrations of some chemical attractants. 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], where the evolution of the density of cells is described by a parabolic equation, and the concentration of a chemoattractant is generally given by a parabolic or elliptic equation, depending on the different regimes to be described and on authors' choices. The behavior of this system is quite well known now: 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] and references therein. However, a drawback of this model is that the diffusion leads alternatively to a fast dissipation or an explosive behavior, and prevents us to observe intermediate organized structures, like aggregation patterns.

keywords and phrases: hyperbolic system with source, chemotaxis, stationary solutions with vacuum, finite volume methods, well-balanced scheme. 1 Istituto per le Applicazioni del Calcolo "Mauro Picone", Consiglio Nazionale delle Ricerche, via dei Taurini 19, I-00185 Roma, Italy (roberto.natalini@cnr.it). 2 Laboratoire J. A. Dieudonné, UMR CNRS 6621, Université de Nice-Sophia Antipolis, Parc Valrose, F-06108 Nice Cedex 02, France & Project Team COFFEE, INRIA Sophia Antipolis, France (ribot@unice.fr). 3 INRIA Sophia Antipolis -Méditerranée, OPALE Project-Team, 2004, route des Lucioles -BP 93, 06902 Sophia Antipolis Cedex, France (monika.twarogowska@inria.fr).

In this paper, we consider a quasi-linear hyperbolic system 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, namely the formation of blood vessels networks during the embryonic development. The model forms a hyperbolic-parabolic system for the following unknowns: the density of endothelial cells ρ(x,t), their momentum ρu(x,t) and the concentration φ(x,t) of a chemoattractant. In one space dimension the system reads    ρ t + (ρu) x = 0, (ρu) t + ρu 2 + P (ρ) x = -αρu + χρφ x , φ t = Dφ xx + aρ -bφ.

(1.1) Loosely speaking, the movement of cells is directed by the gradient of the chemical mediator and is slowed down by the adhesion with the substratum. The positive constants χ and α measure respectively the strength of the cells response to the concentration of the chemical substance and the strength of the friction forces. Overcrowding of cells is prevented by a phenomenological density dependent pressure function given by the pressure law for isentropic gases

P (ρ) = ερ γ , γ > 1,ε > 0. (1.2)
Besides, the evolution of chemoattractant is given by a linear diffusion equation with a source term which depends on ρ: the chemoattractant is released by the cells, diffuses in the environment and is linearly degraded. The positive parameters D,a,b are respectively its diffusion coefficient, the production rate and the degradation rate and the production term is proportional to the cell density. This model was introduced to mimic 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 vascular, endothelial cells. Randomly seeded on the plain gel substratum, these cells migrate and interact together via chemotaxis signaling and, after a while, they aggregate to form a network of capillaries. To be more precise, different final stages are observed depending on the size of the initial density of cells. For low densities, only isolated, disconnected clusters are formed. Increasing the number of cells, a sharp percolative transition occurs and a capillary-like network appears with a characteristic chord length independent of the size of the initial density. Further increase of the number of cells leads to a continuum crossover characterized by the accumulation of additional cells on the network chords until the structure is no longer visible. Finally, for very high initial densities, a continuous carpet of cells with holes, the so called "Swiss cheese" configuration, is observed. From a mathematical point of view, emerging structured patterns, such as capillary-like networks, may be seen as the appearance of nonconstant asymptotic solutions with vacuum, namely solutions composed of regions where the density is strictly positive and of regions where the density (of cells) vanishes.

To reproduce the biological setting, we consider system (1.1) on a bounded domain [0,L] with no-flux boundary conditions, namely, for all time t > 0, ρ x (0,t) = ρ x (L,t) = 0 for the density, ρu(0,t) = ρu(L,t) = 0 for the momentum and φ x (0,t) = φ x (L,t) = 0 for the chemical concentration. From the analytical side, considering the Cauchy problem on the whole space, but in all space dimensions, it is possible to prove the global existence of smooth solutions, if the initial data are small perturbations of a strictly positive constant state, see [START_REF] Di | Analysis and Numerical Approximations of Hydrodynamical Models of Biological Movements[END_REF][START_REF] Cristiana | Existence and asymptotic behavior of solutions to a quasilinear hyperbolic-parabolic model of vasculogenesis[END_REF]. For the present case of the one-dimensional boundary value problem, when the differential part is linearized, it is possible to prove the global existence and the time asymptotic decay of the solutions, when the initial data consist of small perturbation of stable equilibrium constant states, see [START_REF] Francesca | Stability of constant states of qualitative behavior of solutions to a one dimensional hyperbolic model of chemotaxis[END_REF]. However, the analytical study in the present setting is still undone and is difficult in the presence of vacuum, since the hyperbolic part of the model degenerates as the eigenvalues coincide when the mass density vanishes, but see [START_REF] Jang | Well-posedness for compressible Euler equations with physical vacuum singularity[END_REF][START_REF] Gu | Well-posedness of 1-D compressible Euler-Poisson equations with physical vacuum[END_REF] for some rigorous results about the local existence of solutions for related models without chemotaxis.

Actually, preliminary numerical simulations show that, even if we start from strictly positive initial data, vacuum may appear in finite time. This occurrence is much more physically relevant with respect to the analogous situation in gas dynamics. If vacuum is not expected to appear really in gases, it is fully relevant when dealing with the density of cells, i.e.: there are admissible regions without cells, and in some sense it is the main goal of a biologically consistent model. This is a situation somewhat similar to what occurs when dealing with flows in rivers with shallow water type equations [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF]. Besides, also looking at the numerical approximation, dealing with vacuum needs for a special care, since we have to guarantee for the non negativity of the solutions [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF]. This can be understood at the level of the associated numerical flux. A numerical flux resolves the vacuum if for all values of the approximate solution, it is able to generate nonnegative solutions with a finite speed of propagation. In this paper it will be crucial to use schemes with this kind of property.

Another obstacle to a serious numerical exploration is given by the possible lack of proper resolution of nonconstant steady states. Actually, at least for not too large initial masses, solutions to system (1.1) are expected to stabilize for large times, around some global steady states of the system and, when the system reaches equilibrium, the flux is expected to vanish. However, most of the current schemes fail to reproduce this behavior. For instance, a classical centered discretization for the source is not precise enough near steady states. As well known, see for instance [START_REF] Natalini | An asymptotic high order mass-preserving scheme for a hyperbolic model of chemotaxis[END_REF][START_REF] Gosse | Maxwellian decay for well-balanced approximations of a supercharacteristic chemotaxis model[END_REF], this is a usual problem for schemes dealing with hyperbolic problems with source. This is why other approaches have been introduced to balance properly the fluxes and the sources, so giving a more accurate approximation at equilibria. In the case of a semilinear model, obtained from (1.1) by neglecting the drift term (ρu 2 ) x and taking a linear pressure, i.e. γ = 1, Natalini and Ribot [START_REF] Natalini | An asymptotic high order mass-preserving scheme for a hyperbolic model of chemotaxis[END_REF] have proposed an Asymptotically High Order method [START_REF] Aregba-Driollet | Asymptotic high-order schemes for 2 × 2 dissipative hyperbolic systems[END_REF] adapted to the case of a system with an external source term and set on a bounded interval. This technique increases the accuracy of the scheme for large times and yields a correct asymptotic stabilization near the nonconstant equilibria. In [START_REF] Gosse | Asymptotic-preserving and well-balanced schemes for the 1D Cattaneo model of chemotaxis movement in both hyperbolic and diffusive regimes[END_REF][START_REF] Gosse | Maxwellian decay for well-balanced approximations of a supercharacteristic chemotaxis model[END_REF], Gosse has studied the same problem with the well-balanced technique [START_REF] Greenberg | A well-balanced scheme for the numerical processing of source terms in hyperbolic equations[END_REF][START_REF] Gosse | Space localization and well-balanced schemes for discrete kinetic models in diffusive regimes[END_REF] in the framework of finite volume schemes, obtaining similar results. However, both techniques are limited to diagonalizable systems and are difficult to extend to the quasilinear model (1.1). Let us recall also that, in the case of the quasilinear system with linear pressure γ = 1 and without considering the presence of vacuum, Filbet and Shu [START_REF] Filbet | Approximation of hyperbolic models for chemosensitive movement[END_REF] have used the Upwinding Sources at Interface (USI) methodology [START_REF] Botchorishvili | Equilibrium schemes for scalar conservation laws with stiff sources[END_REF] to obtain a wellbalanced scheme. Their main ingredient was the hydrostatic reconstruction, which was originally applied to the Saint-Venant system by Audusse et al. [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF], implying the preservation of steady states with vanishing velocity. Actually, no boundary conditions were considered in [START_REF] Filbet | Approximation of hyperbolic models for chemosensitive movement[END_REF], and so the scheme was not tested against nonconstant steady states.

The main goal of this paper is quite different, and consists in introducing a new scheme which is able to deal both with vacuum problems and with problems near nonconstant steady states arising from the interaction with the source, which is composed by the friction and the forcing chemotactic terms. The strategy is to adapt the ideas in [START_REF] Bouchut | Upwinding of the source term at interfaces for Euler equations with high friction[END_REF], to design a scheme which is able to balance the numerical fluxes with the source term, considering new interface variables and vanishing velocities. Unlike [START_REF] Filbet | Approximation of hyperbolic models for chemosensitive movement[END_REF], we treat the damping term through the interface variables instead of using additional, fractional steps to integrate it in time. We prove that our scheme preserves the non negativity of the density and the stationary states with a vanishing velocity. Finally, the scheme we consider is able to treat vacuum states at the free boundary, which was not considered in any of the previous works.

Before we propose and study our scheme, we make a preliminary investigation on nonconstant stationary solutions with vacuum for system (1.1) and, in particular, on single bump solutions, which are positive only on one connected region. We restrict ourselves mostly to the case of a quadratic pressure γ = 2, which is the simplest one for finding explicit expressions. Under these two restrictions, we give a complete description of the stationary solutions with vacuum. We also numerically show that these solutions are stable and that they can be found as asymptotic states of the system (1.1), even in the case of strictly positive initial data. Other configurations with several regions of positive densities are also found numerically as asymptotic solutions of the system, but we cannot determine for the moment which are the parameters leading to one configuration or another. These results show that model (1.1) gives a realistic description of the vessels formations, which can be successfully tuned against experimental data. This paper is organized as follows: in Section 2, we describe the stationary solutions in the case of γ = 2 and of one region of positive density. Both cases of a lateral bump and of a centered bump are considered. Then, in Section 3, we explain the construction of the scheme we use to discretized system (1.1) and we prove that it preserves the nonnegativity of the density, and the stationary states with vanishing velocities. Finally, in Section 4, we present some numerical results, studying the accuracy of our scheme and showing that the stationary solutions of Section 2 are numerically stable. We also explore the dependence of asymptotic solutions on the parameters of the system, especially on the initial mass and on the value of the adiabatic coefficient γ.

Stationary solutions

In this section, we analyze the existence and the structure of stationary solutions to system (1.1) defined on a bounded interval [0,L] with no-flux boundary conditions given in one space dimension by

ρ x (0,•) = ρ x (L,•) = 0, ρu(0,•) = ρu(L,•) = 0, φ x (0,•) = φ x (L,•) = 0. (2.3)
Remark first that considering the evolution problem (1.1) with the previous boundary conditions (2.3) implies that the mass of the system is constant in time, namely that

M = [0,L] ρ(x,0)dx = [0,L]
ρ(x,t)dx, for all t ≥ 0.

(2.4)

Therefore, the mass M will be considered in what follows as a parameter which characterizes stationary solutions.

2.1. Preliminaries. The steady states of (1.1) are the solutions of the following stationary system (ρu) x = 0, (2.5a)

ρu 2 + P (ρ) x = -αρu + χρφ x , (2.5b)
-Dφ xx = aρ -bφ.

(2.5c) Using the boundary condition we have that ρu = 0, and equation (2.5b) becomes

P (ρ) x = χρφ x , (2.6
) which is automatically satisfied if ρ and φ are constant. More generally, for a pressure P (ρ) = ερ γ with γ > 1, we have two possible solutions for equation (2.6), which are ρ = 0 or

ρ γ-1 (x) = χ(γ -1) εγ φ(x) + K,
where K is an integration constant. Remark that functions ρ and φ are also constrained to respect relation (2.5c). As a consequence, nonconstant equilibria can be composed of intervals where the density is strictly positive, which we call bumps, and of intervals where the density vanishes. However, we are not able to determine a priori the number of such intervals as a function of the data. This is why we fix the number of intervals with positive density and we denote it by p ∈ N. In addition, we denote by x i the boundaries of these intervals, setting x 0 = 0 and x 2p-1 = L, so that the general form of nonconstant steady states is

ρ(x) = p k=1 ρ k (x)χ [x 2k-2 ,x 2k-1 ] (x),
where χ A is the indicator function of the interval A. On the intervals of the form [x 2k-2 ,x 2k-1 ], the density satisfies the relation

ρ k (x) = χ(γ -1) εγ φ k (x) + K k 1 γ-1
, where φ k is the solution to

-D (φ k ) xx =    aρ k -bφ k for x ∈ [x 2k-2 ,x 2k-1 ], k = 1,...,p, -bφ k for x ∈ [x 2k-1 ,x 2k ], k = 1,...,p -1.
Finding φ k from the above equations for a general adiabatic exponent γ > 1 is equivalent to solving a second order, nonlinear differential equation of the form

y xx = A 1 y + A 2 y 1 γ-1 + A 3
, where A 1 ,A 2 ,A 3 are constants. This is a non trivial task apart from the case γ = 2.

Moreover, a unique explicit solution can be obtained only in the case of one interval where the density is strictly positive. Otherwise there are more constants to determine than available equations. More precisely, we need to find 4p -2 constants coming from the integration, 2p -2 interface points x and p constants K. It gives 7p -4 parameters to determine. On the other hand, we have only 6p -3 equations coming from the boundary conditions, the continuity of φ and its first and second space derivatives and from the total mass conservation condition. These two quantities match only in the case of p = 1. For p > 1, we would have an infinite number of stationary conditions depending on p -1 parameters.

Here, we restrict ourselves to the case of only one interval where the density is positive. In this particular case, we can determine exactly the stationary solution as a function of the initial mass of the density M , the length of the domain L and the other parameters of the system α, χ, ε, γ, a, b and D.

2.2.

Quadratic pressure γ = 2: the case of a lateral bump. We give a precise characterization of the equilibria for P (ρ) = ερ 2 when we restrict our attention to solutions with only one bump, i.e.: p = 1. One of the possible forms is a lateral bump, that is: there exists x ∈ (0,L) such that the density satisfies

ρ(x) =    χ 2ε φ(x) + K, for x ∈ [0, x], 0, for x ∈ (x,L], (2.7)
and φ satisfies equation (2.5c). In the following proposition, we completely describe these steady states, as a function of the mass M of the density defined in equation (2.4).

Proposition 2.1. Consider system (2.5) on the interval [0,L] with boundary conditions (2.3) and with a quadratic pressure

P (ρ) = ερ 2 . If τ = 1 D aχ 2ε -b > 0 and L > π √ τ
, there exists a unique, positive solution with a density of the form (2.7) and of mass M . This stationary solution is defined by

ρ(x) =    χ 2ε φ(x) + K, for x ∈ [0, x], 0, for x ∈ (x,L],
(2.8a)

and φ(x) =              2εbK τ χD cos( √ τ x) cos( √ τ x) - aK τ D , for x ∈ [0, x], - 2εK χ cosh( b D (x -L)) cosh( b D (x -L)) , for x ∈ (x,L].
(2.8b)

The free boundary point x is given by the only value

x ∈ 1 √ τ (π/2,π) such that b τ D tan( √ τ x) = tanh( b D (x -L)) (2.8c)
and the constant K is equal to

K = D b M τ 3/2 tan( √ τ x) - √ τ x .
(2.8d)

If τ < 0, or τ > 0 but L < π √ τ
, then there is no one bump solution to the problem.

Proof. To prove the particular form of the equilibrium (2.8), let us first insert (2.8a) into (2.5c). Assuming τ = 1 D aχ 2ε -b > 0 and using boundary conditions φ x (0,•) = φ x (L,•) = 0, the concentration φ can be written as

φ(x) =      C 1 cos( √ τ x) - aK τ D , for x ∈ [0, x], C 2 cosh( b D (x -L)), for x ∈ (x,L].
Using the continuity of functions φ, φ x and φ xx at the interface point x, we obtain the values of the constants C 1 and C 2 , namely

C 1 = 2εbK τ χD cos( √ τ x) -1 , C 2 = - 2εK χ cosh( b D (x -L)) -1 .
which yield (2.8a) and (2.8b). The remaining equation gives the relation (2.8c) for the location of the interface x. The parameter K can be calculated from the total mass density

M = L 0 ρ(x)dx = x 0 χ 2ε φ(x)dx + K x giving (2.8d
). Now we have to determine whether there exists a solution to (2.8c)

or not. The function f (x) = b τ D tan( √ τ x) -tanh b D (x -L) is not defined on 1 √ τ π 2
+ πZ and has a positive derivative elsewhere. Its value

f (0) is positive and the sign of f ( 1 √ τ (π + kπ)), k ∈ Z, is the same as the sign of L - 1 √ τ (π + kπ). So, if there exists k ∈ N such that L ≥ kπ √ τ , f has exactly one zero in each interval of the form 1 √ τ ((π/2,π) + π{0,1,...,k -1}). All
these zeros are candidates to be a boundary x. However, only the smallest of these points, namely the one belonging to the interval 1 √ τ (π/2,π), guarantees the positivity of the function ρ defined by (2.8a). Remark that in the case when L < π √ τ , there is no stationary solution of the form (2.8).

Finally, in the case τ < 0, the same computation leads to write the concentration φ under the form

φ(x) =      C 1 cosh( √ -τ x) - aK τ D , for x ∈ [0, x], C 2 cosh( b D (x -L)), for x ∈ (x,L],
(2.9)

and the continuity of φ,φ x ,φ xx at x, gives the following relation for the interface point x

- b Dτ tanh( √ -τ x) = tanh b D (x -L) .
(2.10)

However, this equation has no solution for x ∈ (0,L), since the left-hand side is negative and the right-hand side positive.

2.3.

Quadratic pressure γ = 2 : case of a centered bump. Let us consider the case of a centered bump, namely the case of one interval where the density is strictly positive in the interior of the domain. Denoting by x, ȳ ∈ (0,L) the boundaries of this interval, with x < ȳ, the density of the nonconstant steady state can be written as

ρ(x) =        0, for x ∈ [0, x), χ 2ε φ(x) + K, for x ∈ [x, ȳ], 0, for x ∈ (ȳ,L].
(2.11)

In the following proposition, we describe these stationary solutions as a function of the mass M and of the length of the domain L.

Proposition 2.2. Let us consider the system (2.5) set on the interval [0,L] with boundary conditions (2.3) and with a quadratic pressure

P (ρ) = ερ 2 . If τ = 1 D aχ 2ε -b > 0 and L > 2π √ τ
then there exists a unique, positive solution with a density of mass M and of the form (2.11). This solution is given by the following expressions:

ρ(x) =        0, for x ∈ [0, x), χ 2ε φ(x) + K, for x ∈ [x, ȳ], 0, for x ∈ (ȳ,L],
(2.12a) and

φ(x) =                                      - 2εK χ cosh( b D x) cosh( b D x) , for x ∈ [0, x), 2εbK τ χD cos( √ τ (x -L 2 )) cos( √ τ (x -L 2 )) - aK τ D , for x ∈ [x, ȳ], - 2εK χ cosh( b D (x -L)) cosh( b D x) , for x ∈ (ȳ,L].
(2.12b)

The free boundary point ȳ is such that ȳ = L -x and so the solution is symmetric and the boundary x is given by the only value

x ∈ L 2 - π √ τ , L 2 - π 2 √ τ such that tanh( b D x) = b τ D tan( √ τ (x -L/2)). (2.12c)
The constant K is equal to

K = M τ (2x -L) b D -2 b D tanh( b D x) . (2.12d)
For other values of the parameters there are no solutions of this form.

Proof. We follow the same computations as in the case of a lateral bump of Subsection 2.2. Under the assumption τ > 0 and using boundary conditions (2.3), the concentration φ can be written as

φ(x) =              C 1 cosh( b D x), for x ∈ [0, x), C 2 cos( √ τ x) + C 3 sin( √ τ x) - aK τ D , for x ∈ [x, ȳ], C 4 cosh( b D (x -L)), for x ∈ (ȳ,L].
Solving the system of three equations given by the continuity of φ, φ x and φ xx at the interface point x (resp. ȳ) gives an expression for the three constants C 1 (resp. C 4 ), C 2 and C 3 as a function of x (resp. ȳ). Therefore, we obtain two expressions for C 2 and C 3 , which give us two non linear equations for x and ȳ, namely:

                   sin( √ τ x)tanh( b D x) -sin( √ τ ȳ)tanh( b D (ȳ -L)) = b Dτ cos( √ τ ȳ) -cos( √ τ x) , cos( √ τ x)tanh( b D x) -cos( √ τ ȳ)tanh( b D (ȳ -L)) = b Dτ sin( √ τ x) -sin( √ τ ȳ) .
This system can be easily rewritten as

       tanh( b D x) = tanh( b D (ȳ -L))cos( √ τ (ȳ -x)) - b Dτ sin( √ τ (ȳ -x)), b Dτ = tanh( b D (ȳ -L))sin( √ τ (ȳ -x)) + b Dτ cos( √ τ (ȳ -x)).
(2.13)

From (2.13) we get the relation

tanh( b D x) = b Dτ cos( √ τ (ȳ -x)) -1 sin( √ τ (ȳ -x)) = -tanh( b D (ȳ -L)),
which implies that ȳ = L -x with x ∈ (0,L/2

∈ L 2 - π √ τ , L 2 - π 2 √ τ .
The parameter K is computed thanks to the value of the mass M of the density and is given by equation (2.12d).

Remark 2.3. Computing solutions (x, ȳ) to systems (2.13) is, in general, a non trivial task. In the case of one lateral bump, we found the conditions that guarantee the existence of x giving nonnegative density everywhere. Nevertheless, due to the nonlinearity of equation (2.8c), we have to solve it numerically to find the value of x. In the case of one centered bump, we find out that the solution is symmetric, which simplifies the computations. The conditions for the existence of solutions can be found and the explicit solution can be calculated. However, the technique we use in that proof cannot be generalized to a higher number of bumps.

Remark 2.4. Notice that the solution for the lateral bump calculated with L/2 and M/2 and symmetrized to be defined on the whole interval [0,L] is equal to the one computed for the centered bump. Our computation was essentially aimed to exclude other cases.

Numerical approximation

Let us now explain how to construct a reliable scheme in order to perform numerical simulations of system (1.1). As a standard guess, we can expect that solutions stabilize on steady states. This scheme will also enable us to test the stability of the stationary solutions computed in the previous section.

System (1.1) couples equations of different natures, i.e. a quasi-linear system of conservation laws with sources, coupled with a linear parabolic equation for the evolution of the chemoattractant. The parabolic part can be approximated using, for example, the classical explicit-implicit Crank-Nicholson method. Now, denoting by U = (ρ,ρu) t the vector of the two unknowns, density and momentum, the hyperbolic part of system (1.1) can be written in the following form

U t + F (U ) x = S(U ), (3.14a)
where F is the flux function and S the source term, i.e.

F (U ) = F ρ F ρu = ρu ρu 2 + P (ρ) , S(U ) = 0 -αρu + χρφ x . (3.14b)
In this section we present a finite volume scheme for (3.14) defined on a bounded domain [0,L] with no-flux boundary conditions (2.3). The scheme needs to preserve the non negativity of density and all the steady states of the system.

3.1. Well-balanced scheme. According to the framework of finite volume schemes, we divide the interval [0,L] into

N cells C i = [x i-1/2 ,x i+1/2 ), cen- tered at nodes x i , 1 ≤ i ≤ N .
In the following, we will assume, for simplicity, that all the cells have the same length ∆x = x i+1/2 -x i-1/2 . We consider as a semi-discrete approximation of the solution U of system (3.14) on cell C i an approximation of the cell average of the solution at time t > 0, that is to say

U i (t) = 1 ∆x x i+1/2 x i-1/2 U (x,t)dx.
A general semi-discrete, finite volume scheme for (3.14) can be defined as

∆x d dt U i (t) + F i+1/2 (t) -F i-1/2 (t) = S i (t), (3.15)
where F i+1/2 (t) is an approximation of the flux F (U (x i+1/2 ,t)) at the interface point x i+1/2 at time t and S i (t) is an approximation of the source term S(U ) on the cell C i at time t.

A classical choice is to take F i+1/2 (t) = F(U i (t),U i+1 (t)), where F is any consistent C 1 numerical flux function for the homogeneous problem U t + F (U ) x = 0. The numerical source is given as S i (t) = S(U i (t)), where we discretized the derivative φ x with a space centered formula

φ x|C i ∼ φ i+1 -φ i-1 2∆x
, where φ i±1 is an approximation of function φ at points x i±1 . However, it is known that this kind of approximation produces large errors near nonconstant steady states.

Balancing the flux term and the source term increases significantly the accuracy near steady states, by imposing an exact discretization of the stationary solutions of the system. A possible approach is to calculate the flux terms F i±1/2 in (3.15) as a function of new interface variables U ± i+1/2 , i.e.

F i+1/2 = F(U - i+1/2 ,U + i+1/2
). These interface variables will be made precise later on and their computation will take into account the balance between the flux term and the source.

The technique is also to upwind the source term, defined as

S i = S - i+1/2 + S + i-1/2
, in the spirit of the USI method [START_REF] Botchorishvili | Equilibrium schemes for scalar conservation laws with stiff sources[END_REF][START_REF] Perthame | Convergence of the upwind interface source method for hyperbolic conservation laws[END_REF][START_REF] Th | Upwinding sources at interfaces in conservation laws[END_REF][START_REF] Bouchut | Upwinding of the source term at interfaces for Euler equations with high friction[END_REF]. We consider the following ansatz:

S - i+1/2 = 0 P ρ - i+1/2 -P (ρ i ) , S + i-1/2 = 0 P (ρ i ) -P ρ + i-1/2 . (3.16)
This ansatz is motivated by an exact discretization of the stationary part of system (3.14a), F (U ) x = S(U ), using that, in the case of a stationary solution, the momentum ρu vanishes thanks to boundary conditions.

Therefore, the final scheme can be written as :

∆x d dt U i + F(U - i+1/2 ,U + i+1/2 ) -F(U - i-1/2 ,U + i-1/2 ) = S - i+1/2 + S + i-1/2 .
(3.17) We will precise in the following subsection how to reconstruct the interface variables U ± i+1/2 , which contain information about the sources. This will be done according to the local equilibrium in order to make the scheme consistent with (3.14), preserving the non negativity of the density and preserving the steady states of (3.14) with a vanishing velocity.

3.2. Reconstruction. In order to complete the construction of the scheme, we define the interface variables U ± i+1/2 . To satisfy the well-balanced property and increase the accuracy of the approximation near nonconstant steady states, the reconstruction is obtained from the stationary system (ρu) x = 0, ρu 2 + P (ρ) x = -αρu + χρφ x .

(3.18)

The system (3.18) can be rewritten in terms of the internal energy function e(ρ) which, for a pressure law of isentropic gas dynamics (1.2), is defined by e ′ (ρ) = P (ρ) ρ 2 . We consider the function

Ψ(ρ) := e(ρ) + P (ρ) ρ = εγ γ -1 ρ γ-1 , for γ > 1,
and we divide the second equation of (3.18) by ρ, which leads to :

   (ρu) x = 0, u 2 2 + Ψ(ρ) -χφ x = -αu. (3.19)
We integrate now the previous system on [x i ,x i+1/2 ] (resp.[x i+1/2 ,x i+1 ]) to find the interface variables

U - i+1/2 = ρ - i+1/2 ρ - i+1/2 u - i+1/2 resp. U + i+1/2 = ρ + i+1/2 ρ + i+1/2 u + i+1/2
, and we obtain the two following equations for ρ - i+1/2 u - i+1/2 and ρ - i+1/2 :

ρ - i+1/2 u - i+1/2 = ρ i u i (3.20a) and εγ γ -1 ρ - i+1/2 γ+1 + C i (ρ - i+1/2 ) 2 + 1 2 ρ 2 i u 2 i = 0, (3.20b)
using (3.20a), with

C i = - εγ γ -1 ρ γ-1 i + χ(φ i -φ i+1/2 ) - 1 2 u 2 i + α x i+1/2 x i u(x)dx.
For an integer γ > 1, equation (3.20b) is a polynomial of order larger than two. The main difficulty in the reconstruction lies in finding its roots and checking that the form of U ± i+1/2 leads to a consistent scheme that preserves the non negativity of the density. Audusse et al. in [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF] introduced the hydrostatic reconstruction for shallow water equations, assuming that the velocity is zero at the steady states. This hypothesis simplifies equation (3.20b) such that an explicit solution can be found. This method generates a scheme, which is well-balanced at equilibria with vanishing velocity. Remark that, for a system of type (1.1) set on a bounded domain with no-flux boundary conditions (2.3), such equilibria are the only possible stationary solutions of the system.

However, in the quasilinear model of chemotaxis (1.1), the source in the momentum balance equation contains a damping term together with a chemotaxis term. The assumption u = 0 at a steady state cancels the friction term in the reconstruction and, to satisfy the consistency property, this term has to be added separately in the discretization of (3.14), namely :

∆x d dt U i + F i+1/2 -F i-1/2 = S - i+1/2 + S + i-1/2 + ∆x 0 -αρ i u i .
This approach was used by Filbet and Shu in [START_REF] Filbet | Approximation of hyperbolic models for chemosensitive movement[END_REF] in the case of a linear pressure function, i.e. γ = 1. In order to include the friction term into the reconstruction, we define the new interface variables taking into account stationary solutions with constant, instead of vanishing, velocity. This allows to deal with the damping term together with the chemotaxis term and the flux term, but also to simplify equation (3.20b). Therefore, considering a velocity constant in space and integrating equations (3.19), we obtain the two following relations for u - i+1/2 and ρ - i+1/2 :

u - i+1/2 = u i (3.21a)
and

Ψ(ρ - i+1/2 ) -Ψ(ρ i ) = -α x i+1/2 x i u(x)dx + χ(φ i+1/2 -φ i ). (3.21b)
Remark that for γ > 1 the function ρ → Ψ(ρ) is strictly increasing and continuous on [0,+∞) with a finite value at 0. So, there exists an inverse function Ψ -1 , which enables us to find a solution to this last equation.

It remains now to explain how to discretize the integral in (3.21b) and how to find the approximation φ i+1/2 . The integral in (3.21b) can be discretized by any consistent method, for example -α

x i+1/2 x i u(x)dx ≈ -α∆x(u i ) + ,
where X + = max(0,X), X -= min(0,X). The computation of φ i+1/2 is not completely plain. The values φ i , which approximate the function φ at points x i , are easily computed thanks to the parabolic equation for φ and we use these values to calculate φ i+1/2 . In order to preserve the non negativity of the density ρ, we take φ i+1/2 = min(φ i ,φ i+1 ). Other choices, as for instance taking the average between φ i and φ i+1 , do not guarantee this property.

In conclusion, the reconstruction of the densities becomes

   ρ - i+1/2 = Ψ -1 [Ψ(ρ i ) -α∆x(u i ) + + χ(min(φ i ,φ i+1 ) -φ i )] + , ρ + i+1/2 = Ψ -1 [Ψ(ρ i+1 ) + α∆x(u i+1 ) -+ χ(min(φ i ,φ i+1 ) -φ i+1 )] + , (3.22
) where the positivity-preserving truncations guarantee the non negativity of ρ.

3.3.

Properties of the semi-discrete scheme. In the following theorem, we prove some properties of the semi-discrete scheme defined by equations (3.17 (i) is consistent with (3.14) away from the vacuum, (ii) preserves the non negativity of ρ, (iii) preserves the steady states given by (3.19) with a vanishing velocity.

Proof. (i) To prove the consistency of the numerical scheme with system (3.14), we first need to show the consistency of the flux term, i.e. ∀i = 1,2,...,N lim

U i ,U i+1 →U,∆x→0 F U - i+1/2 ,U + i+1/2 = F (U ).
This is straightforward using Taylor expansions, since

U - i+1/2 = U i + O(∆x) and U + i+1/2 = U i+1 + O(∆x)
, and, therefore,

F U - i+1/2 ,U + i+1/2 = F (U i ,U i+1 ) + O(∆x).
The consistency finally comes from the consistency of the numerical flux F with the analytical flux F , i.e. F(U,U ) = F (U ).

We have now to prove the consistency of the discretization of the source term. To do so, we use the definition given by Perthame and Simeoni in [START_REF] Perthame | Convergence of the upwind interface source method for hyperbolic conservation laws[END_REF] and we show that ∀i = 1,2,...,N lim

U i ,U i+1 →U,∆x→0 1 ∆x S - i+1/2 + S + i+1/2 = S(U ).
Using equations (3.16), we find that

S - i+1/2 + S + i+1/2 = 0 P (ρ - i+1/2 ) -P (ρ i ) + P (ρ i+1 ) -P (ρ + i+1/2 )
.

We use now Taylor expansions of P (ρ ± i+1/2 ). We consider the density away from vacuum, so that, for ∆x small enough, the positivity-preserving truncations in the reconstruction (3.22) can be omitted and we obtain

P (ρ - i+1/2 ) -P (ρ + i+1/2 ) = P (ρ i ) + ρ i (-α(u i ) + + χ 2 (φ x,i -|φ x,i |))∆x -P (ρ i+1 ) -ρ i+1 (α(u i+1 ) -- χ 2 (φ x,i+1 + |φ x,i+1 |))∆x + O(∆x 2 )
which proves (i).

(ii) We write the first component of scheme (3.17) :

∆x d dt ρ i (t) + F ρ (U - i+1/2 ,U + i+1/2 ) -F ρ (U - i-1/2 ,U + i-1/2 ) = 0.
To prove the conservation of the positivity of the density for our scheme, let us show that, whenever ρ i (t) vanishes, the inequality

F ρ (U - i+1/2 ,U + i+1/2 ) -F ρ (U - i-1/2 ,U + i-1/2
) ≤ 0, holds true. Considering separately the two cases φ i > φ i+1 and φ i < φ i+1 , we can prove that, if ρ i = 0, then

Ψ(ρ - i+1/2 ) = [Ψ(ρ i ) -α(u i ) + ∆x + χ(min(φ i ,φ i+1 ) -φ i )] + = 0 which implies that ρ - i+1/2 = 0.
In the same way, we can prove that if ρ i = 0, then ρ + i-1/2 = 0. Since the numerical flux F preserves the non negativity of ρ for the homogeneous part of system (3.14), we have

F ρ (U i ,U i+1 ) -F ρ (U i-1 ,U i ) ≤ 0,
whenever ρ i (t) vanishes, which completes the proof of (ii).

(iii). We consider a discrete version of the stationary solutions defined by (3.19) with a vanishing velocity, satisfying therefore

u i = 0, Ψ(ρ i+1 ) -χφ i+1 = Ψ(ρ i ) -χφ i . ( 3
F(U - i+1/2 ,U + i+1/2 ) -F(U - i-1/2 ,U + i-1/2 ) -S - i+1/2 -S + i-1/2 =
0, which is clearly true using the consistency of F and equations (3.14b) and (3.16).

3.4.

Properties of the fully discrete scheme. Let us now consider a time discretization of system (3.14) with a given time step ∆t and discretization times t n = n∆t, n ∈ N. Using a standard approximation of the time derivative in scheme (3.17), the fully discrete scheme can be written as

U n+1 i = U n i - ∆t ∆x F(U n,- i+1/2 ,U n,+ i+1/2 ) -F(U n,- i-1/2 ,U n,+ i-1/2 ) + ∆t ∆x S n,- i+1/2 + S n,+ i-1/2 , (3.24)
where U n i is an approximation of the solution U of system (3.14) on cell C i at time t n , U n,± i+1/2 are the values of the interface variables at time t n , and

S n,- i+1/2 = 0 P ρ n,- i+1/2 -P (ρ n i ) , S n,+ i-1/2 = 0 P (ρ n i ) -P ρ n,+ i-1/2 . (3.25)
We can easily see that the time integration preserves two of the properties proved in Theorem 3.1, namely the consistency and the conservation of the stationary solutions. However, we have to find a suitable stability condition for the scheme, that is to say the relation between the space step and the time step to preserve the non negativity of the density. To establish it, we use the notion of invariant domain by interface given in [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF] and we follow the proof presented in [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF]. Indeed, proving directly that the scheme (3.24) preserves a convex domain, such as the positive half-plane, is a hard task since the stencil of scheme (3.24) is composed of three points. To simplify the computations, we consider the weaker notion of preservation of the domain by interface which consists in proving two inequalities involving two points each. However, it is proved in [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF] that these two notions are equivalent under a slightly more restrictive stability condition linking the time step, the space step and a numerical velocity. We first give the definition of a solver preserving the non negativity by interface for a fully discrete scheme, in the case of a homogeneous system. Definition 3.2. A solver F = (F ρ ,F ρu ) t for the homogeneous system U t + F (U ) x = 0 preserves the non negativity of ρ by interface with a numerical speed σ(U n i ,U n i+1 ) ≥ 0 if whenever the CFL stability condition σ(U n i ,U n i+1 )∆t ≤ ∆x holds, we have

     ρ n i - ∆t ∆x F ρ (U n i ,U n i+1 ) -ρ n i u n i ≥ 0, ρ n i+1 - ∆t ∆x ρ n i+1 u n i+1 -F ρ (U n i ,U n i+1 ) ≥ 0.
Let us assume that the numerical flux we use for the flux discretization preserves non negativity by interface. The following proposition gives the stability condition to conserve this property in the case of system (3.14) with source term. holds, we have

     ρ n i - ∆t ∆x F ρ (U n,- i+1/2 ,U n,+ i+1/2 ) -ρ n i u n i ≥ 0, ρ n i+1 - ∆t ∆x ρ n i+1 u n i+1 -F ρ (U n,- i+1/2 ,U n,+ i+1/2 ) ≥ 0. (3.27)
Proof. Assuming the CFL stability condition (3.26) holds, we have

     ρ n,- i+1/2 - ∆t ∆x F ρ (U n,- i+1/2 ,U n,+ i+1/2 ) -ρ n,- i+1/2 u n,- i+1/2 ≥ 0, ρ n,+ i+1/2 - ∆t ∆x ρ n,+ i+1/2 u n,+ i+1/2 -F ρ (U n,- i+1/2 ,U n,+ i+1/2 ) ≥ 0,
thanks to the Definition 3.2 of the preservation of the non negativity for the homogeneous system. It is equivalent to

       1 + ∆t ∆x u n i ρ n,- i+1/2 - ∆t ∆x F ρ (U n,- i+1/2 ,U n,+ i+1/2 ) ≥ 0, 1 - ∆t ∆x u n i+1 ρ n,+ i+1/2 + ∆t ∆x F ρ (U n,- i+1/2 ,U n,+ i+1/2 ) ≥ 0. ( 3

Numerical results

In this section, we analyze numerically the asymptotic behavior of system (1.1) using the scheme introduced in the previous section. First, we compare three different Riemann solvers, Roe, HLL and Suliciu and we show that Suliciu solver is the most adapted to treat vacuum. Then, in order to study the accuracy of our scheme, we compare it with a standard finite difference method with centered in space discretization of the source term. We will see that our scheme captures better the interface with vacuum, shows less diffusion than the standard finite difference method and gives a proper resolution of nonconstant steady states. Therefore, for the following numerical simulations, we use the Suliciu relaxation solver, described in [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF] for the system of isentropic gas dynamics, with an upwinding of the source term.

This scheme is accurate enough to study numerically the stability of the lateral bump and the dependence of the asymptotic solutions of system (1.1) on some of the parameters of the system. In particular, for γ = 2, we are interested in the asymptotic number of bumps for different lengths of the domain L and values of the chemotactic sensitivities χ. We also compare the behavior of the system for different values of the adiabatic exponent γ and, finally, we study the influence of the initial mass on the structure of asymptotic equilibria, comparing the results for γ = 2 and γ = 3.

4.1.

Comparison of different solvers for the homogeneous part. In order to obtain a numerical approximation using a finite volume scheme we have to calculate numerical fluxes between control cells. This fluxes evaluator is based on solving the Riemann problem at each facet in exact or approximate form. Hyperbolic, homogeneous part of the model (1.1) coincides with the isentropic gas dynamics system for which many Riemann solvers are available. However, due to the presence of the source term some of analytical properties of solutions are modified leading to numerical difficulties that cannot be handle by all known Riemann solvers. In particular, occurrence of vacuum and asymptotic, nonconstant states may cause instabilities and negative values of the density. This is the reason why we compare different approximate Riemann solvers, Roe's method, HLL and Suliciu solvers and explain our choice of Suliciu as the most adapted for numerical analysis in the following sections.

In the first simulation we consider system (1.1) with ε = D = a = b = L = 1, χ = 50, a quadratic pressure γ = 2 and the initial data defined as follows ρ 0 (x) = 1 + sin(4π|x -L/4|), φ 0 = 0, u 0 = 0. Figure 1 presents the L 2 (on the left) and L ∞ (on the right) numerical errors for the three previously mentioned solvers with well-balancing of the source term. Behavior at two different grid sizes ∆x = 0.05 and ∆x = 0.01 is studied. The reference solution is obtained using Suliciu solver and well-balancing reconstruction on the fine grid with mesh size ∆x = 10 -3 . The time step satisfies ∆t = 0.9∆x/λ. We observe similar behavior of the errors for all the solvers. They oscillate at the beginning of the evolution and stabilize with time. There is no significant difference between them, although the errors for Roe's method seem a little bigger than in the case of HLL and Suliciu. In the second test, we increase the adiabatic coefficient γ in the pressure function taking γ = 3. Figure 2 presents the density profiles at asymptotic states containing vacuum and steep gradients near the interfaces between regions where the density is strictly positive and regions where the density vanishes. We observe that Roe's method, which is based on a linearization of the system, fails and produces negative values of the density near vacuum. We mention here that the classical Stager-Warming flux splitting also oscillates at the interface with vacuum. The two other solvers, HLL and Suliciu, are stable, they preserve non negativity of the density and approximate the interface with high resolution. Suliciu's approach is to use a relaxation scheme in which the mass conservation equation is not relaxed. This makes it less diffusive than the HLL solver and allows to capture better contact discontinuities. Moreover, in general, the HLL solver needs very careful wave speed estimates. However, in the case of the isentropic gas equations, it is not clear a priori that there is a difference between these two solvers and, in our tests, we don't observe any significant distinction between these two methods. As our model contains source terms leading to steep gradients of the density and appearance of regions where the density vanishes, we decide to use Suliciu solver, which is adapted to treat vacuum.

Roe

HLL Suliciu Remark 4.1. Numerical simulations indicate that none of the three approximate Riemann solvers that we studied is stable for γ > 2 with large initial masses. Therefore, our numerical analysis of the system (1.1) is performed for initial masses small enough to assure the stability of the scheme. ) manages to capture a particular asymptotic behavior of system (1.1). More precisely, we show that the finite volume approach with an upwinding of the source term behaves better near nonconstant steady states than a classical finite difference centered discretization.

In this subsection, we consider the case of a lateral bump for γ = 2. Figure 3 displays asymptotic density profiles at time T = 100 obtained with various schemes -finite volume with upwinding, finite volume without upwinding and finite difference without upwinding -and the exact solution given by (2.8). On the left, we show the computations with a space step ∆x equal to 0.05 and on the right, equal to 0.01. The initial data are the same as in the previous subsection and the time step satisfies the same stability condition ∆t = 0.9∆x/λ. We see that the finite volume scheme with the wellbalanced property gives clearly the most accurate location of the interface. Indeed, studying the different approximations of the source term, we notice that the centered discretization produces a bigger error at the steady state. Then, comparing finite difference and finite volume, we observe that the finite difference method is characterized by a very high numerical diffusion, in contrast to the finite volume approach. However, decreasing the space step, all the schemes converge to the reference solution, even if the finite volume well-balanced scheme is still the most accurate to give the correct location of the free boundary.

Another important point in the numerical approximation of the solutions of system (1.1) defined on a bounded domain with no-flux boundary conditions is the accuracy of the velocity. In this case, the momentum should vanish at steady states. Figure 4 presents the asymptotic states of the density and the momentum obtained by two different methods : on top, the finite volume well-balanced scheme and on bottom, the finite difference scheme with a centered discretization of the source term. On the right, we ∆x = 0.05 ∆x = 0.01 can see the corresponding residues of the momentum. First, we observe that at time T = 100 the residues are less than 10 -8 and decreasing, which means that the steady state is reached. Then, we see that the approximation of the density is comparable for both schemes, which was also observed in the previous test. However, the finite difference approach produces an error in the momentum profile. Its L ∞ norm is close to one instead of being equal to zero. The two previous subsections justify the choice of a finite volume wellbalanced scheme with a Suliciu solver adapted to treat vacuum. This scheme is accurate enough to study the stability of the stationary solutions and the asymptotic behavior of the system and we use it in the following subsections. In what follows, the results are displayed at time T = 100 with a space step equal to 0.01 and a time step satisfying the previous stability condition. [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF].3. Stability of a lateral bump. The first issue we study is the stability of the stationary solutions computed exactly in subsection 2. We consider here system (1.1) with a quadratic pressure function P (ρ) = ερ with one lateral bump exists. The initial mass is taken equal to M 0 = 1 + 1 π . Assuming that the initial datum is not symmetric, the equilibrium has the form of the lateral bump (2.8) with the interface point x ≈ 0.3943. This value is obtained by solving numerically (2.8c).

We analyze the stability of this steady state under two different types of perturbations. First, we perturb the location of the interface point x. More precisely, we take initial data of the form (2.8), in which x is replaced by x * = x + δ. We also recalculate the parameter K, the density ρ and the concentration φ in order to have a perturbation with a zero mass. The parameter δ cannot be too large, since the definition (2.8) would lead to negative values of the density. This type of perturbation modifies the solution on its whole support. In the second test, we change only the density profile in a small region [x 1 ,x 2 ] satisfying 0 < x 1 < x 2 < x such that the initial density is defined as follows :

ρ 0 (x) =            ρ(x), for 0 < x ≤ x 1 , ρ(x 1 ), for x 1 < x ≤ x * , ρ(x 2 ), for x * < x ≤ x 2 , ρ(x), for x 2 < x ≤ x, 0, for x < x ≤ L, (4.29)
where ρ is the exact solution given by (2.8). The location of the jump x * ∈ (x 1 ,x 2 ) is chosen such that the mass of the perturbation is zero.

For the first type of perturbation with δ = 0.1, the results are presented at Figure 5. On the left, we can see the profiles of density (on top), concentration (in the middle) and momentum (on bottom) for the initial perturbed data in dotted blue lines (.), the asymptotic solution after perturbation in solid green lines (-) and the exact solutions in lines marked with red crosses (+). On the right, the corresponding residues are displayed. Initially, the residues of density have values of order 10 -2 , which suggests that the solution is still evolving. Then, at some point, they decrease rapidly to nearly 10 -14 and stabilize at this value, which means that the solution has reached asymptotically the expected steady state. We remark that the asymptotic profiles match perfectly the expected solutions, which confirms that the stationary solution with a lateral bump given by (2.8) is stable under this kind of perturbation. In Figure 6, we display the results obtained with the second type of perturbation (4.29) with [x 1 ,x 2 ] = [0.6x,0.8x], with the same curves as at Figure 5. We notice that the numerical results confirm the stability of the stationary solution composed of a lateral bump. The type of perturbation is different, but the mechanisms of convergence to equilibrium have similar features.

4.4. Dependence on the parameters χ and L of the asymptotic states. Now, in this section, we study which stationary solutions to system (1.1) are reached asymptotically as a function of some parameters of the system and, in particular, how many bumps the asymptotic profile contains. We showed analytically that, if τ ≤ 0, the equilibrium is given by a constant state. The same situation occurs if τ > 0 and L ≤ π/ √ τ . If we increase the size L of the domain, nonconstant stationary solutions with several bumps may appear. However, we are unable to determine analytically how many regions of positive density the asymptotic solution will contain and a numerical study is necessary here. More precisely, our aim is to analyze how the chemotactic sensitivity χ and the length of the domain L influence the number of bumps in the case of a quadratic pressure γ = 2. In Figure 7, we display asymptotic profiles of the density obtained for different lengths L of the domain. The initial data are taken equal to :

ρ 0 (x) = ξ (1 + sin(4π|x -L/4|)), φ 0 = 0, u 0 = 0,
where ξ is computed as a function of the length L of the domain in order to keep the initial mass constant and equal to M = 1.3183. We observe that if the length of the domain is not large enough, that is L ≤ π/ √ τ , the asymptotic state is constant, as expected. Increasing the size of the domain from L = 1 to L = 4 leads to the concentration of particles at the boundary and nonconstant solutions composed of one bump appear. A further increase, from L = 7 to L = 30, allows new bumps to appear and a two bumps solution is observed. We observe the same phenomenon when we increase the value of the chemosensitivity constant χ. We present at Figure 8 different asymptotic profiles obtained for different values of χ. Passing from χ = 3 to χ = 5 leads to a higher concentration of particles which, as a result, produces free space for new bumps. However, from χ = 5 to χ = 200, the global form of the solution, i.e. two lateral bumps on each side of the domain remains the same, even if the two bumps tend to become higher and narrower. We see that the free space available for the cells has an essential effect on the formation of nonconstant steady states. If the size of the domain is too small, only constant solutions are expected. A growth of the available space, by increasing directly the length of the domain or by concentrating particles in smaller regions, leads to the formation of new bumps. However, at some point, a further increase of the chemotactic force, even up to important values, does not produce any new bumps. So, we think that, on each domain, there is a maximum number of bumps once all the parameters are fixed. 4.5. Dependence on the adiabatic coefficient γ of the asymptotic states. The adiabatic exponent γ describes the response of cells to compression. For a high value of γ, the internal pressure repealing the cells is very strong. For example, for ideal gases, the value of γ is taken between 1 and 2, whereas in the Saint-Venant system, describing geophysical flows, it is equal to 2. In the case of cells, it should be much larger as they are less compressible than gas or water molecules.

In Section 2, we presented a general form of nonconstant steady states with several bumps for arbitrary γ. However, we obtained explicit solutions only in the case γ = 2. In this subsection, we study how the equilibria change for different values of γ and especially how the number of bumps varies.

We consider system (1.1) with ε = 1,D = 0.1,a = 20,b = 10,χ = 10 defined on the interval [0,3] and with initial data ρ 0 (x) = (1.5 + sin(4π|x -L/4|)), φ 0 = 0, u 0 = 0.

In Figure 9, we plot the asymptotic solutions for density (on the left) and chemoattractant concentration (on the right) for the following values of the adiabatic coefficient : γ ={2 (blue), 3 (green), 4 (red), 5 (cyan)}. We observe how the number of bumps changes, namely 4 bumps for γ = 2, 3 for γ = 3, 2 for γ = 4 (among which 2 lateral bumps in these 3 cases) and a constant profile for γ = 5. Large values of γ imply strong repealing forces at higher densities. It prevents the formation of high concentrations of cells. When the value of γ increases, the height of bumps decreases, while their support enlarges. Moreover, when the distance between the supports of two neighboring bumps becomes zero, they join together and, for γ high enough, the pressure forces are stronger than the chemotactic movement and constant steady states are observed. 4.6. Dependence on the initial mass of the asymptotic states. Experiments with endothelial cells performed by Serini et.al [START_REF] Gamba | Percolation, morphogenesis, and burgers dynamics in blood vessels formation[END_REF] showed that a vascular-like network develops only if the initial density of cells ranges from 100 to 400 cells/mm 2 . Below this interval, a disconnected structure is observed, while, above, a continuous carpet of cells with holes appears. System (1.1) was constructed to model the phenomenon of early formation of blood vessels and its solutions should also reflect this dependence on the initial mass. These experimental results suggest a particular behavior of steady states with several bumps of system (1.1). The regions where the density is strictly positive may correspond to the location of capillaries and would become thicker for large initial masses. Moreover, we expect to find a threshold value of the initial mass, above which constant equilibria are observed.

We performed some simulations of system (1.1) with γ = 2 and γ = 3, which are presented at Figure 10. On the left (resp. on the right), we can see the asymptotic profiles for density in the case γ = 2 (resp. γ = 3) for different values of the initial mass. The initial data are the same equal to ρ 0 (x) = ξ (1 + sin(4π|x -L/4|)), φ 0 = 0, u 0 = 0, and the initial density is only multiplied by different constants in order to change the value of the initial mass. In the first case γ = 2, we notice that the number of bumps is equal to 4 and remains exactly the same for the different masses. Moreover, the support of the bumps is also independent of the initial mass, whereas the height of the bumps increases with the initial mass. This can be seen theoretically, since the equations (2.8c) or (2.12c) determining the interface point x do not depend on the mass, while the other equations to determine the stationary density do. We notice that this behavior is not the one expected if we consider the experimental observations. However, in the second case γ = 3, the dependency on the mass is totally different and fits the experiments mentioned above. Indeed, when the mass increases, the supports of the bumps become larger and the bumps join together, until reaching the constant equilibrium for a mass large enough.

Theorem 3 . 1 .

 31 )-(3.16) with the reconstruction (3.21a)-(3.22). Let us consider the system (3.14) set on the interval [0,L] with boundary conditions(2.3). Let F = (F ρ ,F ρu ) t be a consistent, C 1 numerical flux preserving the non negativity of ρ for the homogeneous part of system (3.14), U t + F (U ) x = 0. The finite volume scheme (3.17)-(3.16) with the reconstruction (3.21a)-(3.22):

Proposition 3 . 3 .

 33 Let us assume that the homogeneous flux F preserves the non negativity of ρ by interface. Then the fully discrete scheme(3.24) with the reconstruction at interfaces given by (3.21a)-(3.22) preserves the non negativity of ρ by interface, which means than whenever the CFL stability condition σ(U n,- i+1/2 ,U n,+ i+1/2 )∆t ≤ ∆x(3.26)

Figure 1 .

 1 Figure 1. Time evolution of L 2 (on the left) and L ∞ (on the right) numerical errors of the density ρ for system (1.1) in the case P (ρ) = ερ 2 , with ε = D = a = b = L = 1, χ = 50, approximated using a finite volume, well-balanced scheme (3.24). Three different Riemann solvers: Roe method (in solid blue line), HLL solver (in dashed green line), Suliciu relaxation solver (in dotted red line) are compared for ∆x = 0.05 and ∆x = 0.01 .

Figure 2 .

 2 Figure 2. Density profiles for system (1.1) at time T = 5 in the case P (ρ) = ερ 3 , with ε = 1,D = 0.1, a = 20, b = 10, χ = 10, M = 1 approximated using a finite volume, well-balanced scheme (3.24). Three different approximate Riemann solvers: Roe method (on the left), HLL solver (in the middle), Suliciu relaxation solver (on the right) are compared.

4. 2 .

 2 Accuracy of the numerical approximation for the source term. Now, we analyze how the finite volume numerical scheme (3.24)-(3.25) with the reconstruction (3.21a)-(3.22

Figure 3 .

 3 Figure 3. Density profiles for system (1.1) with P (ρ) = ρ 2 , total mass M = 1 and L = 1. The other parameters of the system are D = 0.1, a = 20, b = 10, χ = 10. Comparison between different methods of approximation of a source term: (green) VW -finite volume method with well-balanced reconstruction; (pink) VC -finite volume method with centered in space discretization of the source term, (blue) DC -finite difference scheme with centered in space approximation of the source. Reference solution (red) is given by the exact solution (2.8)

Figure 4 .

 4 Figure 4. On the left: Asymptotic density (blue) and momentum (green) profiles for system (1.1) in the case P (ρ) = ερ 2 and ε = D = a = b = L = 1, χ = 50 obtained by: (on the top) a finite difference scheme with centered in space approximation of the source ; (on the bottom) a finite volume, well-balanced scheme. On the right: the corresponding residues of the momentum.

2

 2 and with the following parameters: ε = D = a = b = L = 1 and χ = 50. This choice guarantees that τ = aχ 2εD -b D > 0 and L > π/ √ τ so the nonconstant steady state

Figure 5 .

 5 Figure 5. On the left: Profiles for system (1.1) with γ = 2 and ε = D = a = b = L = 1, χ = 50 of density (on top), concentration (in the middle) and momentum (on bottom) : (+++ red) exact solutions (2.8); (... blue) initial data with a perturbation of the interface point x * = x + 0.1; (-green) asymptotic solutions corresponding to the previous initial data. On the right: corresponding residues.

Figure 6 .

 6 Figure 6. On the left: Profiles for system (1.1) with γ = 2 and ε = D = a = b = L = 1, χ = 50 of density (on top), concentration (in the middle) and momentum (on bottom) : (+++ red) exact solutions (2.8); (... blue) initial data perturbed in the domain [x 1 ,x 2 ] = [0.6x,0.8x];(-green) asymptotic solutions corresponding to the previous initial data. On the right: corresponding residues.

Figure 7 .Figure 8 .

 78 Figure 7. Density profiles of the asymptotic states of system (1.1) with P (ρ) = ερ 2 , χ = 3, ε = D = a = b = 1 and total mass M = 1.3183 for different values of the length of the domain L = {1,5,7,30}.

Figure 9 .

 9 Figure 9. Profiles of the density ρ (on the left) and of the chemoattractant concentration φ (on the right) at asymptotic states of system (1.1) with ε = 1,D = 0.1,a = 20,b = 10,χ = 10 on the interval [0,3] for different values of the adiabatic exponent γ ={2 (blue),3 (green),4 (red),5 (cyan)}.

Figure 10 .

 10 Figure 10. Density profiles in the case P (ρ) = ρ 2 (on the left) and P (ρ) = ρ 3 (on the right) for system (1.1) with ε = 1, D = 0.1, a = 20, b = 10, χ = 10 and an initial datum of the form ρ 0 (x) = ξ(1 + sin(4π|x -L/4|)). Comparison for different initial masses with ξ ={0.1(blue), 1 (green), 5 (red), 10 (cyan)}.

Acknowledgement. The authors thank François Bouchut for some useful suggestions. This work has been partially supported by the project PORAbruzzo and by the ANR project MONUMENTALG, ANR-10-JCJC 0103.