A SURVEY ON WELL-BALANCED AND ASYMPTOTIC PRESERVING SCHEMES FOR HYPERBOLIC MODELS FOR CHEMOTAXIS

Chemotaxis is a biological phenomenon widely studied these last years, at the biological level but also at the mathematical level. Various models have been proposed, from microscopic to macroscopic scales. In this article, we consider in particular two hyperbolic models for the density of organisms, a semi-linear system based on the hyperbolic heat equation (or dissipative waves equation) and a quasi-linear system based on incompressible Euler equation. These models possess relatively stiff solutions and well-balanced and asymptotic-preserving schemes are necessary to approximate them accurately. The aim of this article is to present various techniques of wellbalanced and asymptotic-preserving schemes for the two hyperbolic models for chemotaxis.

A first example of organism moving under the effect of a chemoattractant, is a bacteria, Escherichia coli, endowed with a flagella, which follows a run-and-tumble process : the bacteria alternates periods of motion as a straight line and periods of rotations. In an homogeneous medium, the movement of these bacteria may be modeled by a diffusion equation, linked to random walk. When these bacteria feel some chemical fluctuations in the medium, the random walk is biased and the run phases may be extended or shortened, see [25,95] for examples of models for this phenomenon. Another bacteria, Bacillus subtilis, shows some complex patterns of collective behavior, like dendritic ramifications, under the influence of nutrients, see [91] for a model and numerical simulations of this phenomenon.
A second kind of organism using chemotaxis is amoebae, like Dictyostelium discoideum. When a threshold of concentration is reached, these amoebae change their frequency of emission of a chemoattractant, called AMPc, from stochastic to oscillatory, leading to the formation of a whole fruiting body.
We may also mention the vasculogenesis process, that is to say the formation of blood vessels by epithelial cells through the emission of a chemoattractant, called VEGF-A. In vitro experiments have been performed [99] and have shown that vascularized networks develop only if the initial density of cells lies between two threshold values. A review on mathematical models for the formation of vascular networks can be found in [97] and 3D numerical simulations in [30]. 0.2. Different types of mathematical models for chemotaxis. Many mathematical models have been proposed and studied to describe those chemotaxis phenomena. One may consult the following reviews and books on the subject : [68,69,89,80,65].
The most famous model is the following Patlak-Keller-Segel system [75,88], where the evolution of the cells density ρ is described by a parabolic equation and the chemoattractant concentration φ 1 is given by a parabolic or an elliptic equation : (1) ∂ t ρ = ∇ · (D ρ ∇ρ − χρ∇φ), where D ρ is the diffusion coefficient of the cells, χ is the chemosensitivity coefficient of cells, D is the diffusion coefficient of the chemoattractant, a is the production rate and b is the degradation coefficient of chemoattractant. Note that the quantities ρ and φ are coupled through the term χ∇ · (ρ∇φ) in the first equation, which describes the movement of cells according to the gradient of chemoattractant and are also coupled through the term aρ in the second equation, which expresses the production of chemoattractant by the cells. Note that the parameter δ may be equal to 0 (resp. 1) if we consider an elliptic (resp. parabolic) equation for φ. The behavior of system (1) with δ = 0 is well-known : in dimension 1, the solution is always global in time [81,85], whereas in higher dimensions, the solution may exist globally in time when the norm of the initial datum is small or may blow up when the norm of the initial datum is large [14,72]. This dichotomy may be explained by the competition between the two terms in the equation for ρ : the advection term leads to a concentration of the density of the cells, on the contrary of the diffusion term.
When δ = 1, the study of the system is not so plain [24] and shows some unexpected behaviors : blow-up in 1D for an appropriate non-linear diffusion coefficient [33], global in time 2D solutions for some large initial data [13], blow-up of 2D solutions for some small initial data [101]...
From a biological point of view, blow-up of cells density is not pertinent and we aim at limiting high concentrations and observing the formation of patterns or permanent structures. Modifications of Keller-Segel model (1) have therefore been proposed, see a review in [65], for example by introducing a "volume-filling" effect, which takes the form of a non-contant chemosensitivity coefficient in (1). One may prove the global-in-time existence of solutions, even for dimensions higher than 1 [27,10]. Another choice is to introduce a non-linear pressure term in the diffusion, which takes into account the volume of cells which cannot be considered as simple material points [76,23]. We obtain the following modified Keller-Segel model : where p is a phenomenological function depending on density. The choice of a pertinent pressure term in a biological context is not obvious and a possible choice for function p as a first approximation is the pressure law for isentropic gaz : Another possibility of modifications in (1) is to change the nature of the equation for ρ. Some hyperbolic equations for ρ, see [90,34] or some hyperbolic systems, linked with system (1), have been proposed. The idea is that hyperbolic equations or systems will produce less diffusive solutions, that may present a finer description of the behavior of cells, at least in the transient regime. The aim of this article is to explain how to deal accurately with these systems at the numerical level.
Let us consider first the semi-linear hyperbolic system for ρ obtained in [58] from the following 1D-kinetic discrete models with two velocities ±λ : The two first equations of the system stand for the evolution of the density of two populations of cells : the one evolving with a positive velocity λ and the other with the opposite negative velocity −λ. The source terms express the change of directions of the two populations, either due to natural change or due to the effect of the chemoattractant. With the change of variables ρ = w + z and J = λ(z − w), we obtain the following equivalent 1D semi-linear hyperbolic system [98,58,87,64] : where ρ is the total density of the population of cells and J is the flux of cells. This system is very close to the dissipative wave equation (also called hyperbolic heat equation, Cattaneo equation or P 1 system), except for the chemotactic term. The supplementary unknown of the system with respect to (1) is the mean flux J. This system has been studied at the analytical level in [60,62,66,67,70,71] and numerically in [82] or [50,51] in the 1D case. The multi-dimensional case has also been considered in [38,46,40,63], essentially at a numerical level, and analytically in [37,38]. However, the analytical results are much more difficult to obtain than for equation (1) based on the Euler equation for compressible gaz with a pressure of type (3). Here, the unknown ρu is the cells momentum and α is a friction coefficient with the medium. System (5) may be seen as a linearization near u = 0 of system (6) with J = ρu. This system is difficult to study analytically and numerically since vacuum regions may appear in the case γ > 1. Some analytical results may be found in [37,39] and 1D numerical studies in [84,83] or in [12] and in [47] for the particular case γ = 1. Note that multi-dimensional numerical simulations have been displayed in [6,30].
Finally, we may mention that kinetic systems have been also considered in order, for example, to model efficiently the run-and-tumble process, see [1,86] or [25], where interesting travelling waves solutions are described. 0.3. Links between the mathematical models presented previously. However these models, although of different natures, are tightly linked one to another and share similar features.
First, we may see that stationary solutions are common to several introduced models. For example, on a bounded interval with appropriate boundary conditions or on the whole line with appropriate vanishing functions at infinity, hyperbolic semi-linear system (5) and hyperbolic quasilinear system (6) with a linear pressure (i.e. γ = 1 in eq. (3)) have the same stationary solutions than parabolic Keller-Segel system (1), but their stability is not necessary the same, see [96] for mode details on that subject. Moreover, in the case of a non-linear pressure, stationary solutions of quasi-linear system (6) coincide with stationary solutions of modified Keller-Segel equation (2) : these solutions are described in details in [11] in the case of a bounded interval and have been studied on R in [9,22].
Then, let us mention that some of these models may be seen as the limit of some other models in the appropriate scaling. We consider equation (4) in the parabolic scaling, that is to say t → ε 2 t, x → εx and a ε = ε 2 a, b ε = ε 2 b : The solutions of rescaled system (7) converge towards the solutions of Keller-Segel equation (1), when ε tends to 0. In the same way, we can consider the large time-large friction limit of system (6) : which solutions converge when ε goes to 0 towards the solutions of system (2) with δ = 0, see [36]. Note that all the models considered here have also been derived from kinetic models [98,58,46,41,31,89,87], see also [46] for a derivation of the two hyperbolic systems (5) and (6) from the same kinetic equation, using two different scalings. 0.4. Numerical methods for chemotaxis models. Since all these models present complex behaviors and since they may be difficult to handle asymptotically, numerical studies are particularly important in order to study their behavior.
From a numerical point of view, Keller-Segel system (1) has been massively studied, using different kinds of methods : in [94], a simple finite differences method is proposed, with a upwind technique on the drift term and with an implicit treatment of the diffusion term in 1D and in 2D. Some mixed finite elements methods in 2D have also been used in [79] and a general framework for finite volumes methods for Keller-Segel equation (1) is described and analyzed in [45]. In [6] a finite volume scheme is also derived and analyzed for a degenerate parabolic system in dimension 2 or 3; the corresponding numerical simulations show interesting patterns formation. We may also mention some variations around upwind scheme [44,32] or discontinuous Galerkin methods [43].
However, numerical methods taking into account some particular features of the models are expected to be more accurate than the previous ones. A first class of methods take advantage of the flow gradient structure of the system. For example, in [15] or in [26], schemes based on a gradient-flow reformulation of some generalized Keller-Segel model with general kernels and a non linear diffusion are constructed. In [28], Filippov flux theory, which gives equivalent solutions to the ones built with gradient-flow theory, is applied to derive schemes for multi-dimensional aggregation equations with a general potential 0.5. Well-balanced and asymptotic-preserving schemes. Another idea is to use the knowledge on stationary solutions to construct accurate schemes. The aim of these well-balanced schemes is to preserve exactly the steady states or to increase the accuracy of computation around steady states, by making an efficient balance between the flux term and the source term; see [52] for an interesting review on well-balancing techniques and accurate preservation of stationary solutions. To do so, many different techniques have been introduced : well-balancing techniques [59,52], asymptotic high order techniques [7], upwinding sources at interfaces techniques [18,92,20], etc... All these approaches present different ways to ensure the preservation of stationary solutions, but share the common philosophy of balancing the flux term and the source term. The advantages of these techniques are numerous : by approximating accurately stationary solutions, it enables to handle stiff source terms. Note also that it improves the general convergence of the scheme, see [2,5,3,4] for precise estimates on that point.
Finally, we explore wether the constructed well-balanced schemes are able to reproduce the limits detailed in subsection 0.3, that is to say wether the schemes are asymptotic-preserving, see [74] for a review on asymptotic-preserving methods. Since stationary solutions are expected to be the asymptotic limits of the solutions of the system, we may hope that the scheme are asymptoticpreserving and we check in most of these cases that it is indeed the case. The only remaining concern is to check if the stability condition needed to assess the parabolic limit does not depend too drasticly on the parameter ε.
Note that, although most of the models have been set in the introduction in the multi-dimensional setting, the numerical schemes presented in the article are all designed in the one-dimensional case. The two systems (5) and (6) are studied on a bounded domain [0, L] complemented with no-flux boundary conditions, namely, for system (5): and for system (6), We adopt in the sequel of the article the following notations : we consider systems (5) and (6) on interval [0, L]. We denote by ∆t the time step and by ∆x the space step used in the discretization of the systems. The discrete times will be defined by t n = n∆t, n ∈ N and the discrete points by Eventually, in the finite volume framework, we will consider some cells centered in The article is organized as follows : in section 1, we present two schemes for the semi-linear hyperbolic system (5), an asymptotic high order scheme [82] and a well-balanced scheme [50] and we compare them. In section 2, we also present two schemes for the quasi-linear hyperbolic system (6), an upwinding sources at interfaces (USI) scheme [84,83] and an approximated Riemann solver [12]. We complete this article by a brief review of well-balancing techniques for kinetic equations and by a conclusion in section 3.

1.
Well-balanced and asymptotic-preserving schemes for semi-linear system (5) with a linear pressure term As explained in [82], a standard upwind scheme applied to the equations of system (4) fails in computing the solutions of the system. Indeed, stationary solutions admit high gradient values, which are difficult to compute in the presence of viscosity and which lead to non vanishing values for the flux J. We therefore need to add in the scheme some techniques that ensure a better accuracy when the solution is close to stationary solutions.
A complete study on stationary solutions for Keller-Segel system (1), which are the same as for semi-linear system (5), may be found in [96]. Some examples of stationary solutions are displayed in Figure 1. Two different techniques are explained in details in this section. They both use system (4), which is equivalent to system (5) in 1D expressed in diagonal variables. In both cases, we assume that function φ is an auxiliary function satisfying equation : Let us denote by φ n i an approximation of φ(x i , t n ) and by Φ n the vector composed of all the components φ n i . Function φ is computed by some classical finite differences numerical scheme, for example with a time implicit-explicit scheme : where I is the identity matrix and M is the matrix for the discretization of the 1D Laplacian with Neumann homogeneous boundary conditions. We will present the two schemes with α = 1 for simplicity reasons.
1.1. Asymptotic High Order (AHO) schemes for the semi-linear system (5). The asymptotic high-order scheme presented here and explained in [82] is derived in the framework of finite differences technique. The scheme is not well-balanced, in the sense that stationary solutions are not exactly preserved by the scheme, but the order of the scheme improves when the solution of the system becomes close to a stationary solution. The high order around stationary solutions is achieved using two ingredients : the approximation of the source on a larger stencil and the computation of the truncation error for stationary solutions.
We rewrite system (4) as follows : where f = χρ ∂ x φ is considered as a source term, that is computed in a foregoing step.
Let us denote by ω = w z and consider system (10) in the following form According to the finite differences formalism, ω n i and F n i stand for approximations of ω(x i , t n ) and F (x i , t n ), and we consider schemes to approximate system (11) under the general viscous form [49]: where B and D are 2 × 2 matrices acting on the vectors ω n i+ and F n i+ respectively at the points x i+ with = −1, 0, 1. The transport part of the scheme is discretized by an upwind strategy with artificial viscosity q ∈ R.
Remark the discretization of the source term of scheme (12) as This approximation uses a three-points stencil corresponding to the stencil related to the upwinding of the transport parts. It also introduces new degrees of freedom through the 6 matrices B and D , that we will use in order to increase the order of the scheme around stationary solutions. These matrices should satisfy natural conditions, like : in order to have the consistency of scheme (12) with system (11), and (14) in order to have the monotonicity of the scheme (see [7]).
Second-and third-order schemes. The main point of scheme (12) is now to compute the local truncation error for a stationary solution denoted byω. We setω i =ω(x i ). Using a Taylor expansion and the fact that (Λ −1 B) 2 = 0, we obtain : We therefore add to conditions (13) the two following conditions such that the ∆x term cancels and such that the second-order accuracy holds on every stationary solution: Due to a symmetry invariance of system (5) that we want to keep at the discrete level, other conditions on the coefficients of the matrices are also required, see [82] for more details.
Finally, one possible scheme (with q = λ) can be rewritten as: which for this particular example turns out to be exactly the scheme first proposed by Roe [93].
For monotonicity reasons, we should have the following conditions on the time step and the space step : ∆x ≤ 4λ and ∆t ≤ ∆x λ + ∆x/4 . Now, if we compute the truncation error one order further, we get : Remark that it is not possible to cancel straightforwardly the constant term, the ∆x term and the ∆x 2 term, but we can generalize the consistency condition (13) and the order conditions (15) to have a scheme of third order, as follows : for some suitable matrices C and E and Therefore, once E is chosen, we should take In order to satisfy the consistency condition (17), the monotonicity conditions (14) and the order condition with q = λ, we should take E = 1 12λ 1 1 1 1 with the following conditions on the time and space steps: ∆x ≤ 6λ, ∆t ≤ ∆x λ + ∆x/3 .
We obtain consequently the following scheme : which is of order 3 around stationary solutions.
Boundary conditions. Finally, the last concern of this work is to introduce adapted numerical boundary conditions for the scheme we derived. Note that an important feature of the continuous model (5) is the mass conservation, namely : In [82], we have computed boundary conditions in order to preserve exactly a discrete version of the mass. Namely, denoting by I n 1 a discrete version of I(t n ) = ρ(x, t n ) dx, we propose boundary conditions that ensure exactly I n+1 1 − I n 1 = 0. In the case of Roe's scheme (16), we obtain the following conditions : and in the case of the third-order AHO scheme (18), We can prove [82] that these boundary conditions are consistent with continuous boundary conditions (8).
Properties. To conclude, we may prove, see [82], that the resulting schemes (16) and (18) complemented with the previously derived boundary conditions are L 1 and L ∞ stable and consistent, and therefore convergent. We present in Figures 2 and 3 some numerical results obtained with the Asymptotic High Order schemes presented above.  (red stars) and AHO scheme of order 3 (blue diamonds) and on the right, the same for J.

1.2.
Well-balanced scheme for the semi-linear system (5). Now, we present another numerical scheme, in the finite volumes setting, derived to approximate the solutions of system (5) in [50]. Following the Greenberg-Leroux strategy [59], [55], we approximate the system (4) written as : , by a non-linear and non-conservative system which solutions can be defined, see [50] : where δ is the Dirac function. The source terms present in system (19) lead to the presence of zerowaves, that is to say stationary solutions, in the solutions of the Riemann problem. This allows to define different values on the right and on the left sides of the interface through the computation of the stationary solutions on [0, ∆x], as explained in Figure 4. Note that this 2-speeds scheme can be seen as a restriction to 2 points discretization of scheme based on Case solutions [52]. Another way to interpret the well-balanced strategy is to view it as a splitting scheme in space. The first step is devoted to the resolution of the differential system of stationary solutions on a cell of the discretization. This step leads to the definition of intermediate values at interfaces. The second step computes the exact resolution of the homogeneous transport part using the previous computed values at interface.
In this scheme, as before, φ would still be considered as an auxiliary function, which will be computed by standard methods for parabolic reaction-diffusion equations. Therefore, we consider that ∂ x φ n j+1/2 has been previously computed. Now, as explained before, we solve the following stationary system, knowing the values of boundary conditions z(0) and w(∆x) : With the change of variables ρ = w + z and J = λ(z − w), we obtain the equivalent following system : that is to say the stationary solutions of system (5) which can be solved as : Since we evaluate the expressions at point ∆x which is considered to be small, the second expression can be approximated by a simple expansion and we obtain the following approximation : Rewriting these approximations in the w − z variables and solving the linear system to obtain approximations of w(0) and z(∆x) as functions of z(0) and w(∆x), we obtain From these values, we can deduce some intermediate values at the interface : which are used in a second step :

Properties. The resulting scheme (20)-(21) has the following properties [50] :
• steady states on w and z are exactly preserved, • the subcharacteristic condition |∂ x φ| ≤ λ, necessary to guarantee the existence of solutions to (5), holds at the discrete level, under some suitable conditions, • under CFL and subcharacteristic conditions, L p estimates, for 1 ≤ p ≤ ∞, and BV bounds of the numerical solution w ∆x (t, ·) and z ∆x (t, ·) are proved, • solutions of the scheme converge towards solutions of system (5). The scheme is also asymptotic preserving. Namely, under some parabolic CFL condition ∆t = O(∆x 2 ) and assuming that ε ≤ λ ∆t ∆x , the limit ε → 0 in the implicit-explicit scheme which turns out to be an approximation of system (7) in the parabolic scaling, leads to a discretization of the Keller-Segel equation (1). Although the properties of the scheme have been shown under the subcharacteristic condition, a numerial study of the previous scheme is also performed in [51] in the super-characteristic regime. Note that in [2,3,4,5], some precise and very useful error estimates on well-balanced schemes have been established for scalar balance laws and relaxation systems.
Numerical results for the solutions of system (5) obtained with this scheme in [50] are similar to the ones displayed in Subsection 1.1.

Conclusion.
We have presented in this section two numerical schemes for an hyperbolic semilinear system of chemotaxis (5) : an asymptotic high order scheme and a well-balanced scheme. The two schemes are based on the same idea, that is to say a balance between the flux term and the source term in order to preserve stationary states, which results in stabilizing the numerical flux.
However, we may point out some differences between the two schemes : concretely, the AHO schemes result in a change in the discretization of the source term, whereas the well-balanced technique relies on a change of discretization of the flux. The well-balanced strategy needs the exact computation of the stationary solutions, whereas the AHO scheme requires only a small consistency error for a stationary solution, described by its differential equation. AHO schemes may be of order 2 and 3, when well-balanced are restricted to order 2. However, the mass conservation is straightforward for the well-balanced setting and it should be imposed through appropriate boundary conditions for the AHO scheme.

2.
Well-balanced and asymptotic-preserving schemes for quasi-linear system (6) with a non-linear pressure term In this section, we consider now the quasi-linear system (6) with a non-linear pressure term. Note that, in the previous section, in the case of system (5), the linearity of the transport part was very useful to perform some explicit computations in the derivation of both schemes. Unfortunately, this will not be possible any more in the case of the quasi-linear system (6) and other schemes need to be thought of. We will present two schemes here : an upwinding sources at interface scheme and an approximated Riemann solver.
We consider system (6) in the one-dimensional case, which writes as : set on the interval [0, L] with boundary conditions (9). With these boundary conditions, the mass, denoted by M is also preserved for this system, namely : ρ(x, t) dx, for all t ≥ 0.

Description of the stationary solutions for models with a non-linear pressure term.
To begin with, let us describe the stationary solutions of (22) with a non-linear pressure of the form (3) and with boundary conditions (9). They satisfy the following equations : which leads to the following expressions thanks to the boundary conditions : Equivalently, defining the energy function e such that p (x) = xe (x), that is to say we may rewrite system (23) as : and equations (24) as : In the particular case when γ = 2 the computations become simple, since e is linear and we can therefore replace the expression of ρ as a function of φ in the equation for φ, which is in such a case a linear second-order differential equation, that can be solved explicitly. In the case when τ = 1 D aχ 2κ − b > 0, we obtain : From this expression, we can compute stationary solutions of the system (6) by using these expressions on different neighboring intervals. Hence, we can prove that, in the case when L > π  Figure 5.
To conclude, we see that there are an infinite number of nonconstant stationary solutions and that the set of stationary solutions is complex. Therefore, in order to compute accurately the asymptotic behavior of solutions of system (6) towards one of these stationary solutions, we aim at deriving well-balanced schemes.

2.2.
Upwinding Sources at Interfaces (USI) scheme for the quasi-linear system (6). In [84,83], we proposed and studied a Upwinding Sources at Interfaces (USI) scheme, constructed following the ideas of [18,20,8,47]. The idea is to add some new values of the unknowns at the interface and to compute them in order to preserve exactly the stationary solutions.
Omitting the equation for φ and denoting by U = (ρ, ρu) t the vector of the two unknowns, density and momentum, the hyperbolic part of system (6) can be written in the following form where F is the flux function and S the source term, i.e. (27b) The steady states are given by : We are now considering a finite volume framework : we denote by ρ i (t) and u i (t) approximations of the mean value of functions ρ and u on the cell C i at time t and by ρ n i and u n i the approximations of the mean value of functions ρ and u at time t n on the cell C i . We also define a numerical flux F, consistent with flux F , and we denote by S i (t) an approximation of the mean value of the source term on cell C i and time t. The semi-discrete version of the scheme can be written as : where F i+1/2 (t) is an approximation of the flux F (U (x i+1/2 , t)) at the interface x i+1/2 at time t. A possible choice is to set F i+1/2 (t) = F(U i (t), U i+1 (t)) and S i (t) = S(U i (t)). However, this strategy produces large errors when the solutions come close to the stationary solutions.
Following [18,20,8,47], we compute the fluxes at interfaces F i±1/2 (t) under the form F(U − i+1/2 , U + i+1/2 ), where U ± i+1/2 are new interface variables at x i+1/2 , to be computed thanks to the equations of the stationary solutions. At the same time, we also need to discretize the source term using again the equation of the stationary solutions : we split the source term into two terms S i = S + i−1/2 + S − i+1/2 , with , these expressions being obtained by an integration of the equation (28) or (29) for stationary solutions, respectively on the two intervals Therefore, the semi-discretized scheme can be written as : As explained before, these values are obtained through the integration of equation (28) or equation (29) on [x i , x i+1/2 ] (resp.[x i+1/2 , x i+1 ]). In order to include the damping term −αρu inside the reconstruction, we consider non-constant stationary solutions with constant velocity, instead of vanishing velocity, and we obtain (see [84] for more details) : and, integrating eq.(29), where ρ i+1/2 = 1 2 (ρ i + ρ i+1 ). Note that the positive part and the discretization of ∆x∂ x φ as The fully discretized system can be derived from (31) as : where U n i is an approximation of the mean value of U 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 are the values of the source term approximations (30) at time t n .
Another approach can be considered, see [83], inspired by [47] in the case of a linear pressure function, i.e. γ = 1. The idea is to include the damping term inside the discretization rather than in the reconstruction, as follows : where the new interface variablesÛ ± i+1/2 are computed without the friction term. The interface values are defined in this case by : , using an integration of eq.(29) with null momentum, and by using an integration of eq.(28) with null momentum. The corresponding fully discretized scheme can be obviously deduced from (36). In [83], this approach was studied in order to treat implicitly in time the damping term −αρu, which enables to compute accurately solutions with a stiff momentum.
Properties. Under appropriate conditions, see [84,83], the semi-discretized finite volume scheme (31)- (30) with the reconstruction (32)- (33) or the reconstruction (32)-(34): (1) is consistent with (27) away from the vacuum, (2) preserves the non negativity of ρ, preserves the steady states given by (28) with a vanishing velocity, (4) is asymptotic-preserving in the parabolic scaling with (2). Under suitable conditions on the numerical flux F and under a stability condition of type σ(U n,− i+1/2 , U n,+ i+1/2 )∆t ≤ ∆x the fully discrete scheme (35) with the reconstruction at interfaces given by (32)- (33) or (32)- (34) preserves also the non negativity of ρ. Some comparisons of numerical solutions obtained with various schemes are represented in Figure  6. ∆x = 0.05 ∆x = 0.01 Figure 6. Comparison of the computation of the stationary function ρ for system (6) for different schemes (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. 3. Approximated Riemann solvers for the quasi-linear system (6). Now, let us detail the scheme proposed by Berthon, Crestetto and Foucher in [12], based on an approximated Riemann solver strategy. They construct a well-balanced scheme based on Riemann approximate solver for the quasi-linear system (6) for a non-linear pressure with γ = 2. For that purpose, they use some explicit formulas coming from the computation of the stationary solutions, that can be found in [84] or in [11]. Let us define W = (ρ, ρu, φ). The scheme in [12] is based on an approximate Riemann solver R(x/t, w L , w R ), which is an approximation of the exact Riemann solver W (x, t; w L , w R ) defined as the exact solution at time t of equation (6) with an initial datum under the form The approximate Riemann solver is chosen under the following form : Note that a discontinuity with null velocity has been added and that expression (37) contains some quantities to be defined : the two velocities λ L < 0 < λ R , and the intermediate states w * L and w * R , see Figure 7. Since these two intermediate states are made of three components each, we need to compute 6 quantities. The two velocities will be fixed later on in order to enforce some positivity conditions.
The integral consistency condition 1 ∆x gives 3 equations, providing an approximated computation of the right-hand side of the equality. This computation is usually straightforward and comes from the hyperbolic system under consideration, but it is more complicated in the present case because of the presence of the source term and of the coupling with the equation for φ. Consequently, after an approximation of the source terms involving the quantity φ and some simple computations, the following equations are obtained : Note that the quantities {ρ∂ x φ}, (∂ x φ) R and (∂ x φ) L are not defined yet. Then, the quantities ρu and e(ρ) − χφ are natural invariants across the stationary wave [19,59], as well as the variations of φ, according to the diffusion equation. This leads to the following relations after a linearization of energy function e : Denoting by we may find a solution to equations (38) and (39) under the form : Now, in order to have a well-balanced property, that is to say to have w * R = w R and w * L = w L whenever u L = u R = 0 and e(ρ L ) − χφ L = e(ρ R ) − χφ R , we can define : Then, some considerations on positiveness of densities have to be taken into account : we denote and we can choose λ L < 0 < λ R such that ρ HLL ≥ 0. We hence modify ρ * R and ρ * L ad follows : which ensures the positivity of ρ * * L and ρ * * R and guarantees the first relation of eq. (38). Similarly, we defineφ and we may choose ∆t, λ L and λ R such thatφ ≥ 0. We set Note that, for the moment, the two quantities (∂ x φ) R and (∂ x φ) L of eq. (38) are not defined yet, but the following condition has been used in what preceeds : The approximate Riemann solver we defined satisfies the following properties : solutions remain positive and density and momentum stay at rest, starting from a stationary solution. Finally, the scheme may write as follows, after a L 2 -projection over piecewise constant functions : where W * ,L,R i+1/2 denotes the left and right values of the approximated Riemann solver R(x/∆t, W n i , W n i+1 ). Finally, we define (1 + ε(∆x)) , where ε(∆x) tends to 0 when ∆x → 0 to enhance the consistency condition. In details, ε(∆x) is chosen in order to satisfy the well-balanced property for φ, differentiating cases given by the exact computation (26) of stationary solutions for γ = 2 : Properties. Under some conditions, see [12], among which the CFL-like condition : the resulting scheme has the following properties : • the scheme preserves the non-negativity of ρ and φ, • density and momentum of stationary solution (25) stay at rest, • starting from adapted initial data, the scheme is well-balanced for γ = 2 : W i+1 n = W i n for all i and for all n, • assuming that λ R = −λ L = λ, and under the CFL condition ∆t ε∆x max i∈Z |λ i+1/2 | ≤ 1 2 , the scheme is asymptotic preserving in the parabolic scaling with equation (2). Note that numerical results obtained with this scheme in [12] are similar to the ones presented in Subsection 2.2.
2.4. Conclusion. As before, we remark that the two schemes presented in this section have the same goal, that is to say the accurate computation of the solutions of the scheme by preserving the non-constant steady states. The techniques are slightly different, even if they relie on the computation of intermediate values that take into account the stationary solutions. Let us point out two main differences between the two schemes. The USI scheme of section 2.2 can be computed for any values of γ and only needs the equation of the stationary solution, whereas the approximated Riemann solver of section 2.3 requires the exact computation of the steady states, which can be done only when γ = 2. However, the approximated Riemann solver gives a well-balanced scheme for the full system on (ρ, ρu, φ), whereas the USI scheme is derived only for the hyperbolic part of the system on (ρ, ρu).

Further works and conclusion
Note that this paper concentrates on well-balancing techniques and asymptotic preserving schemes for hyperbolic systems for chemotaxis. Recently, various articles have studied numerical techniques for kinetic systems for chemotaxis. After some first articles [48,100], using for example semilagrangian methods, well-balanced and asymptotic preserving schemes have been derived. For example, an even and odd formalism may be used to construct an asymptotic preserving scheme towards the parabolic modified Keller-Segel system [29] or Volpert calculus enables to derive an asymptotic preserving scheme in the hydrodynamic limit [73]. Some well-balanced schemes may be found in [53] or, more recently in [56], where the authors introduce a scheme, which is at the same time well-balanced, following the approach of [53], and asymptotic-preserving using Volpert calculus, as in [73]. Another technique consists in constructing a two-steps scheme , combining a well-balanced step and an asymptotic-preserving scheme [42].
All the models presented and mentioned in the paper have already had many applications in various domains : modeling pursuit-evasion phenomena [57,78] by considering the interaction of two populations, each of them under the influence of chemotaxis. Recently, some articles also study the chemotaxis models on graphs, in order, for example, to represent the evolution of fibroblasts on the fibers during the wound healing process, see for example [21,61] for the semi-linear system (5) or [16] for equation (1) or [17] for kinetic models.
To conclude, the challenge in the future would be the construction of similar well-balanced schemes in the multi-dimensional case, and to begin with, in the 2D case. The 2D and 3D cases are relevant settings for real models for biology, but the step from 1D well-balanced schemes to 2D well-balanced schemes is stiff. Some first definitions of Godunov schemes in the 2D case have been proposed in [35,54] and are likely to produce new ideas to derive 2D well-balanced and asymptotic-preserving schemes.