Laurent Gosse
email: l.gosse@ba.iac.cnr.

MAXWELLIAN DECAY FOR WELL-BALANCED APPROXIMATIONS OF A SUPER-CHARACTERISTIC CHEMOTAXIS MODEL

Keywords: Chemotaxis modeling, discrete velocity kinetic model, stiffness, super-characteristic relaxation, stabilization in time, WB scheme, reflecting boundary conditions AMS subject classifications. 65M06, 35L60

 which relate the total initial mass of cells with the time-asymptotic values of the chemoattractant substance on each side of the domain. It is shown that the WB scheme stabilizes correctly onto zero-mass flow rate (hence Maxwellian) steady-states which agree with the aforementioned bifurcation diagrams. The evolution in time of residues is commented for every considered test-case.

Introduction.

1.1. Modeling of chemotaxis dynamics. This paper is concerned with the following semilinear model of chemotaxis movement, originally proposed by Greenberg and Alt [START_REF] Greenberg | Stability results for a diffusion equation with functional shift approximating a chemotaxis model[END_REF], where "particles" move at velocities ±λ, λ > 0:

   ∂ t ρ + ∂ x J = 0 ∂ t J + λ 2 ∂ x ρ = ρ∂ x ϕ -J ∂ t ϕ -D∂ xx ϕ = αρ -βϕ. (1.1)
One can rewrite it in a more mathematically tractable way by introducing its diagonal variables (its Riemann invariants, in hyperbolic terminology),

w = 1 2 ρ - J λ , z = 1 2 ρ + J λ ⇔ 0 ≤ ρ = w + z, J = λ(z -w), (1.2)
which satisfy a semilinear system of equations, similar to Cattaneo's model [START_REF] Hillen | Hyperbolic Models for Chemosensitive Movement[END_REF]:

   ∂ t w -λ∂ x w = -1 2λ G(∂ x ϕ; w, z) ∂ t z + λ∂ x z = 1 2λ G(∂ x ϕ; w, z) ∂ t ϕ -D∂ xx ϕ = αρ -βϕ.
(1.3) with non-negative values for D, α, β ≥ 0 and the conventional choice of turning rates,

G(∂ x ϕ; w, z) = (∂ x ϕ + λ)w + (∂ x ϕ -λ)z.
Besides the textbooks [START_REF] Bellouquid | Mathematical Modeling of Complex Biological Systems: A Kinetic Theory Approach[END_REF][START_REF] Perthame | Transport equations in biology[END_REF], general references about mathematical models for chemotaxis dynamics are for instance [START_REF] Chalub | Kinetic models for chemotaxis and their drift-diffusion limits[END_REF][START_REF] Donatelli | Singular convergence of nonlinear hyperbolic chemotaxis systems to Keller-Segel type models[END_REF][START_REF] Filbet | Derivation of hyperbolic models for chemosensitive movement[END_REF][START_REF] Greenberg | Stability results for a diffusion equation with functional shift approximating a chemotaxis model[END_REF][START_REF] Hillen | Hyperbolic Models for Chemosensitive Movement[END_REF][START_REF] Hillen | The one-dimensional chemotaxis model: global existence and asymptotic profile[END_REF][START_REF] Hillen | Existence of weak solutions for a hyperbolic model of chemosensitive movement[END_REF][START_REF] Hillen | Hyperbolic models for chemotaxis in 1-D, Nonlinear Analysis[END_REF][START_REF] Hyung | Global existence of classical solutions for a hyperbolic chemotaxis model and its parabolic limit[END_REF][START_REF] Keller | Initiation of slime mold aggregation viewed as an instability[END_REF][START_REF] Osaki | Finite dimensional attractor for one-dimensional Keller-Segel equations[END_REF][START_REF] Othmer | The diffusion limit of transport equations II: Chemotaxis equations[END_REF][START_REF] Schaff | Stationary solutions of chemotaxis systems[END_REF][START_REF] Segel | A theoretical study of receptor mechanisms in bacterial chemotaxis[END_REF] for the theoretical aspects and [START_REF] Chertock | A second-order positivity preserving central-upwind scheme for chemotaxis and haptotaxis models[END_REF][START_REF] Dolak | Cattaneo Models for Chemosensitive Movement: Numerical Solution and Pattern Formation[END_REF][START_REF] Filbet | Approximation of hyperbolic models for chemosensitive movement[END_REF][START_REF] Guarguaglini | Global stability of constant states and qualitative behavior of solutions to a one dimensional hyperbolic model of chemotaxis[END_REF][START_REF] Natalini | An asymptotic high order mass-preserving scheme for a hyperbolic model of chemotaxis[END_REF] for numerics. According to the analysis by Jin and Katsoulakis [START_REF] Jin | Hyperbolic systems with supercharacteristic relaxations and roll waves[END_REF], the system (1.3) ceases to be quasi-monotone for |∂ x ϕ| > λ and the relaxation process becomes super-characteristic (see also [START_REF] Greenberg | Congestion redux[END_REF]; in particular, useful estimates based on monotonicity as derived in [START_REF] Natalini | Weakly coupled systems of quasilinear hyperbolic equations[END_REF] fail). In such a regime, the decay estimates of the L p norms proved in §3 of [START_REF] Gosse | Asymptotic-Preserving and Well-Balanced scheme for the 1D Cattaneo model of chemotaxis movement in both hyperbolic and diffusive regimes[END_REF] do not hold and one expects an increase of the L ∞ norm of ρ and ϕ. In general, super-characteristic relaxation processes are unstable; however, for the particular model (1.1) inside a bounded domain and supplemented by Neumann boundary conditions on ϕ, it has been established, either analytically or numerically, in [START_REF] Guarguaglini | Global stability of constant states and qualitative behavior of solutions to a one dimensional hyperbolic model of chemotaxis[END_REF][START_REF] Natalini | An asymptotic high order mass-preserving scheme for a hyperbolic model of chemotaxis[END_REF] that certain asymptotic steady-states with zero velocity are stable (even if conventional upwind schemes can grossly miss them). It is therefore interesting to address the question whether or not the WB/AP discretization proposed in [START_REF] Gosse | Asymptotic-Preserving and Well-Balanced scheme for the 1D Cattaneo model of chemotaxis movement in both hyperbolic and diffusive regimes[END_REF] can be shown to be consistent with these peculiar large-time behaviors. This kind of numerical discretization for the model under consideration comes partly from the choice of the computational grid parameters. System (1.1) is concerned with the movements of cells which size is of the order of the micrometer µm (for instance, Escherichia coli is 2µm long). They move roughly at a velocity of the order of the µm per second [START_REF] Perthame | Transport equations in biology[END_REF] inside a box which characteristic length is around 5 centimeters. An aggregate of cells counts at least 100 individuals, meaning that its length should be around 10 -2 cm. Fixing the box as the interval x ∈ [0, 1], we get that in order to ensure that an aggregate is represented on a few grid points, we should have around 250 or 500 discretization points: this is the range of parameters used in the computations of sections 4 and 5. The usual CFL condition for hyperbolic systems relates the fineness of the space discretization with the maximal size of the time-step ∆t; hence, even assuming an unconditionally stable implicit integration of the diffusion equation on ϕ, the number of iterations necessary to reach steady-state will grow if finer grids are set up for (1.1). Using coarser grids fastens the convergence toward steady-states and WB/AP schemes ensure consistency with a correct long-time behavior despite the fact that fine layers may not be resolved.

1.2. Numerical stiffness and non-conservative products. Let us now briefly recall the ideas from [START_REF] Gosse | An asymptotic-preserving well-balanced scheme for the hyperbolic heat equations[END_REF][START_REF] Gosse | Space localization and well-balanced scheme for discrete kinetic models in diffusive regimes[END_REF][START_REF] Gosse | Asymptotic-Preserving and Well-Balanced scheme for the 1D Cattaneo model of chemotaxis movement in both hyperbolic and diffusive regimes[END_REF], which consist in passing from (1.3) to the non-linear and non-conservative (NC) system:    ∂ t w -λ∂ x w = -1 2λ j∈Z hG(∂ x ϕ; w, z)δ x -(j - 1 2 )h ∂ t z + λ∂ x z = 1 2λ j∈Z hG(∂ x ϕ; w, z)δ x -(j - 1 2 )h ∂ t ϕ -D∂ xx ϕ = αρβϕ, (1.4) where δ(.) stands for the Dirac mass in x = 0. This problem can be rewritten like in [START_REF] Isaacson | Convergence of the 2 × 2 Godunov method for a general resonant nonlinear balance law[END_REF], by introducing a steady discontinuous variable a(x) ∈ BV (R), replacing the righthand side of (1.3) by ∓ 1 2λ G(∂ x ϕ; w, z)∂ x a and adding the trivial equation ∂ t a = 0. For discontinuous w, z, this formulation can be unstable because of the products "Heaviside × Dirac"; however, it has been rigorously shown in [START_REF] Gosse | Space localization and well-balanced scheme for discrete kinetic models in diffusive regimes[END_REF] that these nonconservative products can be rigorously defined as weak limits in the framework of [START_REF] Ph | Representation of weak limits and definition of nonconservative products[END_REF] thanks to the uniform BV estimates which come from the linear convection in (1.3) (similar estimates for scalar balance laws are previously given in [START_REF] Gosse | Localization effects and measure source terms in numerical schemes for balance laws[END_REF]).

One may wonder what might be reason for substituting the homogeneous, but non-conservative system (1.4), at the place of the nonhomogeneous semilinear balance laws (1.3). Following [START_REF] Gosse | A Nonconservative Numerical Approach for Hyperbolic Systems with Source Terms: The Weil-Balanced Schemes[END_REF], this reformulation generally offers 2 main advantages:

• a good consistency with the long-time behavior as soon as the non-conservative products are defined in a correct manner. Thanks to the homogeneous character of (1.4), we know that in a bounded domain, the Godunov scheme has zero viscosity at steady-state because it contains only stationary shocks or contact discontinuities (if the domain is infinite, smooth traveling waves can also exist [START_REF] Härterich | Travelling waves in systems of hyperbolic balance laws[END_REF]). As the scheme uses an exact Riemann solver, jump discontinuities are resolved exactly based on Rankine-Hugoniot relations. Therefore, when the stationary discontinuities associated to the non-conservative products are resolved with integral curves of the steady-state equations, it is sure that these curves constitute the only way to pass from a computational cell to its neighbor when the numerical solution is stationary. This jump relation is precisely the one emanating from a limiting process when a(x) is regularized and the weak limit [START_REF] Ph | Representation of weak limits and definition of nonconservative products[END_REF] is studied for simple models: see [START_REF] Gosse | Localization effects and measure source terms in numerical schemes for balance laws[END_REF][START_REF] Gosse | Space localization and well-balanced scheme for discrete kinetic models in diffusive regimes[END_REF]. • a good robustness in stiff regimes (independently of ∆t) which is again a consequence of the fact that the source term is treated by means of a generalized Rankine-Hugoniot jump relation. Indeed, if one divides the source term by a small parameter, say 0 < ε ≪ 1, the only effect is to rescale the x variable (appearing in the steady-state equations) into x/ε, meaning that one goes much farer on its integral curve. However, since only convective waves play a role in the CFL condition, one deduces immediately that ∆t and h are restricted only by the eigenvalues of the Jacobian matrix of the fluxes. The non-conservative product is associated to the zero eigenvalue, thus has no influence on the time-step. This has been illustrated recently in [START_REF] Gosse | Transient radiative tranfer in the grey case: well-balanced and asymptotic-preserving schemes built on Case's elementary solutions[END_REF][START_REF] Gosse | Well-balanced schemes using elementary solutions for linear models of the Boltzmann equation in one space dimension[END_REF]. Stability and convergence properties have been shown for the Godunov scheme built onto the exact Riemann solver for (1.4) in [START_REF] Gosse | Asymptotic-Preserving and Well-Balanced scheme for the 1D Cattaneo model of chemotaxis movement in both hyperbolic and diffusive regimes[END_REF] in the sub-characteristic regime (since ∂ x ϕ is kept constant at each interface, see Remark 2 and [START_REF] Gosse | Asymptotic-Preserving and Well-Balanced scheme for the 1D Cattaneo model of chemotaxis movement in both hyperbolic and diffusive regimes[END_REF], only the hyperbolic part is rigorously exact), meaning that |∂ x ϕ| < λ. Therefore, in this paper, we concentrate on the opposite range of parameters, that is the super-characteristic region. The quantity |∂ x ϕ| can become stiff in such a regime: consequently, both the 2 aspects of WB schemes formerly recalled are likely to be needed.

Remark 1. Concerning the preservation of the stationary regimes, the superiority of schemes preserving all the steady-states has been recently studied in [START_REF] Ying | On the Advantage of Well-Balanced Schemes for Moving-Water Equilibria of the Shallow Water Equations[END_REF]; this property is also crucial in order to ensure that, within a parabolic rescaling, the WB scheme is asymptotic-preserving too. This fact has already been used in various contexts: see [START_REF] Gosse | An asymptotic-preserving well-balanced scheme for the hyperbolic heat equations[END_REF][START_REF] Gosse | Space localization and well-balanced scheme for discrete kinetic models in diffusive regimes[END_REF][START_REF] Gosse | Transient radiative tranfer in the grey case: well-balanced and asymptotic-preserving schemes built on Case's elementary solutions[END_REF]. One may want to use high-order schemes as the system (1.4) is just an usual hyperbolic system without source terms. The USI schemes proposed in [START_REF] Bouchut | Upwinding of the source term at interfaces for Euler equations with high friction[END_REF][START_REF] Filbet | Approximation of hyperbolic models for chemosensitive movement[END_REF] realize a trade-off between a treatment of the source term by means of a simpler jump relation (some flux terms are discarded in the steady-state equation) and overall consistency in the transient regime. They may be more vulnerable in stiff regimes because of the residual terms which are not treated through a jump relation.

Organization of the paper.

In §2, we mainly follow the results from [START_REF] Guarguaglini | Global stability of constant states and qualitative behavior of solutions to a one dimensional hyperbolic model of chemotaxis[END_REF] and recast the Godunov scheme inside a bounded domain with Neumann and Dirichlet boundary conditions (see (2.5) and (4.1)). The well-balanced treatment of boundary conditions is presented in §2.1 while the Appendix A recalls the most essential results from [START_REF] Gosse | Asymptotic-Preserving and Well-Balanced scheme for the 1D Cattaneo model of chemotaxis movement in both hyperbolic and diffusive regimes[END_REF]. A mass conservation property is proved in §2.2. The treatment of the diffusion equation on ϕ is presented in §2.3: implicit discretizations are selected in order to restrict the time step ∆t based on the hyperbolic CFL only. Remark 4 presents the Crank-Nicholson scheme, which is illustrated on a concrete example in Appendix B. Stability results for initial data satisfying smallness assumptions are stated in §3.2. Computational results with Dirichlet boundary conditions are presented in §4 in super-characteristic regime: residues decay rather smoothly and the mass flow rate is of the order of 10 -7 . In §5, we switch to the more delicate Neu-mann boundary conditions and aim at checking the consistency with the bifurcation diagrams reproduced on Fig. 5.1. The stability of a constant state with a big initial mass is studied in §5.2 for two different choices of α, β: one can see that many iterations are necessary in order to create a non-physical mass flux of the order of 10 -12 inside the computational domain. In §5.3, another test-case which is supposed to lead to a non-constant Maxwellian steady-state is presented: this is perhaps the most interesting test-case because the evolution of the residues actually displays 3 different regimes, see Fig. 5.5. The results of §5.4 are concerned with a similar experiment with an even bigger initial mass of cells: the scheme remains stable and mass flow rate inside the domain are shown to decay up to 10 -4 despite residues decay in an oscillatory way. We stress that in this case, |∂ x ϕ| is of the order of 3500 with λ ≃ 10, which is truly stiff. A last test-case is presented in §5.5 with α = β = 50 leads to similar results except that cells concentrate on the opposite side of the computational domain. We strongly insist on the fact that absolutely no modifications are made on the WB Godunov scheme studied in [START_REF] Gosse | Asymptotic-Preserving and Well-Balanced scheme for the 1D Cattaneo model of chemotaxis movement in both hyperbolic and diffusive regimes[END_REF] in order to produce these results except for the boundary conditions on each side of the computational domain.

2. Well-balanced boundary conditions for bounded domains. Let's denote w n j ≃ w(n∆t, jh), z n j ≃ w(n∆t, jh) for any j ∈ {1, 2, ..., 1 h } and n ∈ N; these numerical approximations induce piecewise constant functions w h (t, x) and z h (t, x):

w h (t, x) := w n j , z h (t, x) := z n j for t, x ∈ [n∆t, (n + 1)∆t[× (j - 1 2
)h, (j + 1 2 )h .

The Godunov scheme proceeds in defining a control cell ](j -1 2 )h, (j + 1 2 )h[×]n∆t, (n+ 1)∆t[ around each point x j = jh, solving a Riemann problem on both interfaces (j ± 1 2 )h and averaging: see Fig. 3.1 in [START_REF] Gosse | Asymptotic-Preserving and Well-Balanced scheme for the 1D Cattaneo model of chemotaxis movement in both hyperbolic and diffusive regimes[END_REF]. This notation is also used hereafter:

∀j, n ∈ 1, 2, ..., 1 h × N, (∂ x ϕ) n j+ 1 2 := 1 h ϕ n j+1 -ϕ n j . (2.1)
The piecewise-constant function ϕ h is defined similarly for t, x ∈ R + × [0, 1] whereas the "numerical space derivative" (∂ x ϕ) n j+ 1 2 can be computed on the staggered grid at the same locations than the Dirac masses on the right-hand side of (1.4).

2.1. Derivation of "stopping pressure" boundary conditions. In order to derive the correct expression for boundary conditions to be placed at each edge of the computational domain x ∈ (0, 1), we follow the approach of the "half Riemann problem" [START_REF] Dubois | Boundary conditions for nonlinear hyperbolic systems of conservation laws[END_REF]. Roughly speaking, the idea is to create 2 fictious computational cells containing exactly the right values able to forbid any mass flow in x = 0 and x = 1. We focus on the left side of the computational domain in order to find the correct values w 0 , z 0 which induce a stopping pressure because the reasoning is exactly similar for the right side. According to §3.1 and §3.3 in [START_REF] Gosse | Asymptotic-Preserving and Well-Balanced scheme for the 1D Cattaneo model of chemotaxis movement in both hyperbolic and diffusive regimes[END_REF] (see Appendix A where the main elements are recalled), and with the notation of Fig. 2.1:

∀n ∈ N, w * 1 2 -z n 0 = 0 = w n 1 -z * 1 2 , w * 1 2 + z n 0 = w n 1 + z * 1 2 . (2.2)
The first relation expresses the fact that the mass flow rate vanishes on both sides of the steady discontinuity emanating from the singular source term concentrated in x = h 2 ; the second, that the macroscopic density is constant. Therefore, one finds that z * = w n 1 . The well-balanced scheme for (1.3) reads for any n, j ∈ {1, ..., 1 h },

w n+1 j = w n j + λ∆t h w * j+ 1 2 -w n j , z n+1 j = z n j - λ∆t h z n j -z * j-1 2 , (2.3) x = h 3h 2 h 2 t = n∆t -λ λ w n 1 z n 1 δ x -3h 2 δ x -h 2 w n 1 z * 1 2 w * 3 2 z n 1 w n+1 1 z n+1 1 (only z jumps) (only w jumps) w n 2 z * 3 2
w * supplemented by the edge-values:

(∂ x ϕ) n 1 2 ≡ 0, (∂ x ϕ) n 1 h + 1 2
≡ 0 and,

z n+1 1 = z n 1 - λ∆t h (z n 1 -w n 1 ) , w n+1 1 h = w n 1 h + λ∆t h z n 1 h -w n 1 h . (2.4)
We recall that, according to the Lemma A.1, the interface values read:

       w * j+ 1 2 = 2 1+B n j+ 1 2 w n j+1 + 1-A n j+ 1 2 1+B n j+ 1 2 z n j , z * j+ 1 2 = - 1-B n j+ 1 2 1+B n j+ 1 2 w n j+1 + A n j+ 1 2 +B n j+ 1 2 1+B n j+ 1 2 z n j ,
with the parameters,

A n j+ 1 2 = 1 - h λ exp h λ 2 (∂ x ϕ) n j+ 1 2 , B n j+ 1 2 = 1 + h λ exp h λ 2 (∂ x ϕ) n j+ 1 2 .
Thus both equations (2.3)-(2.4) constitute a well-balanced approximation involving an exact NC Riemann solver for the system (1.1) posed in x ∈ (0, 1) and supplemented by Neumann boundary conditions (as advocated in [START_REF] Guarguaglini | Global stability of constant states and qualitative behavior of solutions to a one dimensional hyperbolic model of chemotaxis[END_REF][START_REF] Natalini | An asymptotic high order mass-preserving scheme for a hyperbolic model of chemotaxis[END_REF][START_REF] Schaff | Stationary solutions of chemotaxis systems[END_REF]): for any t ∈ R + ,

∂ x ϕ(t, x = 0) = ∂ x ϕ(t, x = 1) = 0 = J(t, x = 0) = J(t, x = 1). (2.5)
Observe that (2.2) means that these Neumann conditions are imposed on the staggered grid and not inside the border computational cells centered around x = h and x = 1 h .

Remark 2. In the expression of the coefficients

A n j+ 1 2 and B n j+ 1 2 , one consid- ers that (∂ x ϕ) n j+ 1 2
remains constant; hence comparing with the steady-state curves rigorously derived in [START_REF] Guarguaglini | Global stability of constant states and qualitative behavior of solutions to a one dimensional hyperbolic model of chemotaxis[END_REF], there is a discrepancy and only second order consistency actually holds. The numerical results displayed in §4-5 suggest that this compromise between well-balanced accuracy and overall feasibility is enough. These approximate jump relations were derived in [START_REF] Gosse | Asymptotic-Preserving and Well-Balanced scheme for the 1D Cattaneo model of chemotaxis movement in both hyperbolic and diffusive regimes[END_REF] where the usual sub-characteristic condition was assumed; here we test them in the context of high |∂ x ϕ| super-characteristic regime.

2.2. Mass-preservation (first moment of kinetic density) property. We begin by defining the approximate (piecewise-constant) function ρ h := w h + z h ; the motivation in this subsection is to prove that the scheme (2.3)-(2.4) generates a numerical process which is mass-preserving, even without any reference to the subcharacteristic condition which is crucial for the convergence results in [START_REF] Gosse | Asymptotic-Preserving and Well-Balanced scheme for the 1D Cattaneo model of chemotaxis movement in both hyperbolic and diffusive regimes[END_REF].

Theorem 1. For any h, ∆t ≥ 0, there holds for ρ h generated by (2.3)- (2.4):

∀t ∈ R + , 1 0 ρ h (t + ∆t, x)dx = 1 0 ρ h (t, x)dx. (2.6)
Proof. With obvious notation, we can deduce from (2.3)-(2.4) the numerical scheme acting on ρ n j := w n j + z n j :

ρ n+1 j = ρ n j - ∆t h λ (z n j -w * j+ 1 2 ) (z * j+ 1 2 -w n j+1 ) -λ(z * j-1 2 -w n j ) .
The flux conservation property on each side of the interface x = (j + 1 2 )h which reads

z n j -w * j+ 1 2 = z * j+ 1 2
w n j+1 has been shown in [START_REF] Gosse | Asymptotic-Preserving and Well-Balanced scheme for the 1D Cattaneo model of chemotaxis movement in both hyperbolic and diffusive regimes[END_REF]; it comes from the steady-state equations of (1.1) which are solved for computing the NC jump relations. The integrals in (2.6) are computed by summing the preceding expression on indexes j ∈ {1, 2, ..., 1 h } corresponding to the computational domain x ∈ (0, 1):

1 h j=1 hρ n+1 j = 1 h j=1 hρ n j -∆t λ(z * j+ 1 2 -w n j+1 ) -λ(z * j-1 2 -w n j ) = 1 h j=1 hρ n j + ∆tλ(z * 1 2 -w n 1 ) -∆tλ(z * 1 h + 1 2 -w n 1 h +1 ) = 1 h j=1 hρ n j + ∆tλ(z * 1 2 -w n 1 ) -∆tλ(z n 1 h + 1 2 -w * 1 h + 1 2 ) = 1 h j=1 hρ n j .
From the former subsection, we know that the edge values satisfy z *

1 2 = w n 1 and w * 1 h + 1 2 = z n 1 h + 1 2
thus, for any h, ∆t ≥ 0, we are done.

Since we didn't assume the sub-characteristic condition like in [START_REF] Gosse | Asymptotic-Preserving and Well-Balanced scheme for the 1D Cattaneo model of chemotaxis movement in both hyperbolic and diffusive regimes[END_REF], the system (1.3) isn't quasi-monotone in the sense of [START_REF] Natalini | Weakly coupled systems of quasilinear hyperbolic equations[END_REF], we can't hope to have the dissipation property in time for w h (t, .) L p + z h (t, .) L p for p > 1; in particular, we shall see in the numerical results that the L ∞ norm can grow strongly in certain cases. Concerning positivity-preserving properties, the situation is intricate: as the scheme (2.3)-(2.4) is not a convex combination, there's no hope to derive such a property for arbitrary values of ∂ x ϕ. Actually, we shall see in Remark 7 that negative values appear in transient regimes when the L 1 norm of ρ(t, .) is big enough so as to generate very steep gradients ∂ x ϕ and thus stiffen the problem up to a point where the approximation (2.1) breaks down. Positivity-preserving results, hypothesizing they may exist, will involve a careful analysis of the interaction between the well-balanced scheme on the hyperbolic part of (1.3) and, for instance, the "complete flux scheme" on the diffusion equation (which involves solving its steady-state problem at each interface) [START_REF] Ten Thije | The Finite Volume-Complete Flux Scheme for Advection-Diffusion-Reaction Equations[END_REF].

Remark 3. The proof does not exploit the expression of the coefficients A n j+ 1 2 and B n j+ 1 2 , meaning that (2.6) hopefully holds for general turning rates: see [START_REF] Hillen | Hyperbolic Models for Chemosensitive Movement[END_REF][START_REF] Hillen | Hyperbolic models for chemotaxis in 1-D, Nonlinear Analysis[END_REF][START_REF] Othmer | The diffusion limit of transport equations II: Chemotaxis equations[END_REF].

2.3. Implicit time discretization of the ϕ equation. Since we have in mind to study the asymptotics in time of the well-balanced scheme (2.3)-(2.4) on the kinetic variables w and z, it makes sense to substitute the explicit discretization on ϕ h (used in [START_REF] Gosse | Asymptotic-Preserving and Well-Balanced scheme for the 1D Cattaneo model of chemotaxis movement in both hyperbolic and diffusive regimes[END_REF] in order to be in position to take advantage of the Hoff-Smoller theory with an implicit one. Such a change allows to switch from a time-step ∆t submitted to the parabolic CFL restriction 2D∆t ≤ h 2 to a lighter hyperbolic one λ∆t ≤ h. Denoting Φ n := (ϕ n j ) j∈{1,..., 1 h } , w n := (w n j ) j∈{1,..., 1 h } and z n := (z n j ) j∈{1,..., 1 h } for any n ∈ N, the following equation completes the kinetic discretization (2.3)-(2.4):

Φ n+1 = (1 + β∆t)Id + D∆t h 2 L -1 (Φ n + α∆t w n + z n ) , (2.7)
where Id stands for the 1 h × 1 h identity matrix and L is a discrete approximation of the diffusion term acting on ϕ including convenient boundary conditions. For instance, a first-order Neumann boundary condition on ϕ can be implemented with the matrix,

L N 1 =          1 -1 0 • • • 0 -1 2 -1 . . . . . . 0 . . . . . . . . . 0 . . . . . . -1 2 -1 0 • • • 0 -1 1         
, whereas a more accurate second-order approximation leads to the following one [START_REF] Guarguaglini | Global stability of constant states and qualitative behavior of solutions to a one dimensional hyperbolic model of chemotaxis[END_REF]:

L N 2 =          2 3 -2 3 0 • • • 0 -1 2 -1 . . . . . . 0 . . . . . . . . . 0 . . . . . . -1 2 -1 0 • • • 0 -2 3 2 3          . (2.8)
Concerning the accuracy in time, there is an easy way to improve to second-order: it consists in advancing first by computing (2.7) for n = 0 but with ∆t/2 in place of ∆t. This way, one advances initially of half a time-step and the numerical space derivatives (∂ x ϕ) j+ 1 2 used in the parameters A j+ 1 2 , B j+ 1 2 are systematically computed on a staggered grid n + 1 2 corresponding to (n + 1 2 )∆t. Clearly, at the final step of the computation, one has to repeat the initial trick and advance in time of only half a timestep; this method is sometimes called the Strang time-splitting (or Strang-Marchuk splitting, see e.g. [START_REF] Faragó | Operator Splittings and Numerical Methods[END_REF] for more details).

Remark 4. A slightly more costly manner for treating the diffusion equation on ϕ is the unconditionally stable Crank-Nicholson scheme (sometimes called the "θscheme" with θ = 1 2 ) which consists in solving at each time step:

(1 + β 2 ∆t)Id + D∆t 2h 2 L Φ n+1 = (1 - β 2 ∆t)Id - D∆t 2h 2 L Φ n + α∆t w n + z n ),
with either L N 1 or L N 2 . However, we did not see qualitative differences between the outcome of this scheme and the one generated by (2.7) with a Strang splitting involved.

Residues can decay faster thanks to higher accuracy in time, but not systematically in a monotonic way. In particular, they still can display oscillations (like in the figures 4.2, 5.7 and 5.10): see an example in Appendix B.

3. Stability results for super-characteristic semi-linear relaxation.

3.1. General case of super-characteristic relaxation. In this subsection, we follow the analysis of [START_REF] Jin | Hyperbolic systems with supercharacteristic relaxations and roll waves[END_REF] in order to study the weakly nonlinear asymptotic limit for (1.1). First, for any 0 ≤ ε ≪ 1 and c ∈ R, we define new variables: ξ = 1 ε (xct) and τ = t. Let ρ > 0 be a fixed constant, we introduce the expansions:

ρ(t, x) = ρ + εA(τ, ξ), J(t, x) = ρ∂ x ϕ + εB(τ, ξ),
for some value of ∂ x ϕ. The equations on the traveling waves A and B read:

∂ τ A + 1 ε ∂ ξ (B -cA) = 0, ∂ τ B + 1 ε ∂ ξ (λ 2 A -cB) = A∂ x ϕ -B.
Let F = 1 ε (B -cA) be a modified flux term; this yields,

∂ τ A + ∂ ξ F = 0, ∂ τ F - 2c ε ∂ ξ F + 1 ε 2 ∂ ξ λ 2 ε (ρ + εA) -c 2 A = A ε (∂ x ϕ -c) -F.
We group the terms of the same order in the second equation:

∂ ξ ρ = 0, (λ 2 -c 2 )∂ ξ A -ε∂ ξ (2cF + A) = εA(∂ x ϕ -c) + O(ε 2 ).
Fixing c = ±λ leads to the relation ∂ ξ F = -∂ τ A and:

∂ τ (A + εF ) = (∂ x ϕ -c)A -εF 2c .
Thus, for -|c| < ∂ x ϕ < |c|, the subcharacteristic condition is satisfied and A decays exponentially in time. Oppositely, for |∂ x ϕ| > |c|, the leading terms in the ordinary differential equation on A + εF amplify and the original expansion for ρ and J will eventually break down. This is a super-characteristic regime but since the model (1.1) is semilinear, we do not obtain a Burgers-type equation as in [START_REF] Jin | Hyperbolic systems with supercharacteristic relaxations and roll waves[END_REF], Theorem 3.1.

Stability of constant

Maxwellian steady-states. Despite the preceding results, the computations in [START_REF] Guarguaglini | Global stability of constant states and qualitative behavior of solutions to a one dimensional hyperbolic model of chemotaxis[END_REF] indicate that the model (1.1) can stabilize also in the case of a super-characterisitic relaxation regime. Indeed, independently of the size of |∂ x ϕ|, it is proved that steady-states with zero velocity are stable provided their total mass is not too big and the perturbation oscillates so to have zero mean.

Theorem 2. (Guarguaglini et al. [START_REF] Guarguaglini | Global stability of constant states and qualitative behavior of solutions to a one dimensional hyperbolic model of chemotaxis[END_REF]) Let (ρ, 0, φ) ∈ R + × {0} × R + be a constant steady-state for (1.1) inside the interval x ∈ (0, 1) with Neumann boundary conditions (2.5) satisfying:

ρ < λ 2 α (β + Dπ 2 ). (3.1)
If the initial data (ρ 0 (x) + ρ, J 0 (x), ϕ 0 (x) + φ) ∈ H 1 (0, 1) is a perturbation of the aforementioned constant Maxwellian steady-state such that, for small ε ≥ 0,

1 0 ρ 0 (x)dx = 0, ρ 0 H 1 + J 0 H 1 + ϕ 0 H 1 ≤ ε,
then the corresponding solution (ρ, J, ϕ) of (1.1)-(2.5) decays exponentially toward the constant steady-state (ρ, 0, φ). There exist C, θ > 0 depending on ε such that:

∀t > 0, (ρ, J, ϕ)(t, .) -(ρ, 0, φ) H 1 ≤ C exp(-θt).
Moreover, in case the perturbation is chosen symmetric around x = 1 2 , the threshold (3.1) can be improved. The following result is an extension of Corollary 4.2 in [START_REF] Guarguaglini | Global stability of constant states and qualitative behavior of solutions to a one dimensional hyperbolic model of chemotaxis[END_REF]:

Corollary 1. For ε small enough and some 1 ≤ n ∈ N, assume that the initial perturbation (ρ 0 , J 0 , ϕ 0 ) ∈ C 1 (0, 1) 3 for (1.1)-(2.5) satisfies the symmetry properties:

∀m ∈ {0, 1, ..., n -1},    ρ 0 (x -x) = ρ 0 (x + x) J 0 (x -x) = -J 0 (x + x) ϕ 0 (x -x) = ϕ 0 (x + x)
for x = 2 m-n and x ∈ (0, x).

Then, if moreover 1 0 ρ 0 (x)dx = 0, the conclusion of Theorem 2 still holds with:

ρ < λ 2 α (β + 2 2n Dπ 2 ).
Proof. It proceeds by following the ideas of the proof in [START_REF] Guarguaglini | Global stability of constant states and qualitative behavior of solutions to a one dimensional hyperbolic model of chemotaxis[END_REF]; let us first denote W := (ρ-ρ, J, ϕ-φ) and study the case n = 1. For the choice of the right-hand side of (1.1), we have that a symmetric initial data W 0 (x) gives rise to a symmetric solution around x = 1 2 . Therefore, W (t, x)

x∈(0,x)
and W (t, x)

x∈(x,1)
solve the same Neumann IBVP in their respective intervals because the smoothness assumption implies that ∂ x ρ(t, x) = 0 = ∂ x ϕ(t, x) for all t > 0. Hence we apply Theorem 4.1 of [START_REF] Guarguaglini | Global stability of constant states and qualitative behavior of solutions to a one dimensional hyperbolic model of chemotaxis[END_REF] with L being the measure of the half computational interval (0, x). By induction, the general case n ∈ N is handled similarly by dividing into 2 the computational intervals and the space derivatives vanish thanks to the smoothness assumption. Therefore, these stabilization results do not contradict the instability of supercharacteristic relaxation which is addressed on the whole line. Indeed, from [START_REF] Guarguaglini | Global stability of constant states and qualitative behavior of solutions to a one dimensional hyperbolic model of chemotaxis[END_REF], one sees that the threshold ρ is a decreasing function of the size of the bounded interval inside which (1.1)-(2.5) is solved. However, there is a smallness assumption for ρ 0ρ in H 1 , meaning that beyond a critical value of n, it will become an obstruction for the increase of the threshold level ρ in Corollary 1 because ∂ x ρ 0 L 2 grows with n.

Remark 5. The symmetry conditions in Corollary 1 imply that (but are not equivalent to) W 0 (x) has 2 vanishing moments; indeed, for n = 1, we have

1 0 (x - 1 2 )ρ 0 (x)dx = 1 2 0 x ρ 0 ( 1 2 + x) -ρ 0 ( 1 2 -x) dx = 1 0 ρ 0 (x)dx = 0,
meaning that ρ 0 is orthogonal to the linear subspace of affine functions. For instance, the perturbation x → cos(8πx) is C ∞ , meets all the requirements for n = 2 and,

1 0 cos(8πx)dx = 0 = 1 0
x cos(8πx)dx, 1 0

x 2 cos(8πx)dx = 0.

The perturbation x → sin(4π|x -1 2 |) isn't smooth in x = 1 2 thus does not meet the requirements of the case n = 1 despite its symmetry properties (see Fig. 3

.1).

In case the problem (1.1) is posed on the whole real line x ∈ R, it is stated in §2 of [START_REF] Guarguaglini | Global stability of constant states and qualitative behavior of solutions to a one dimensional hyperbolic model of chemotaxis[END_REF] that stability for constant Maxwellian steady-states holds for ρ < β α λ 2 . 4. Numerical results with Dirichlet conditions on ϕ. For the sake of completeness, we also include the case where the equation on ϕ is posed in the interval x ∈ (-1, 1) and supplemented by Dirichlet boundary conditions:

∀t ∈ R + , ϕ(t, x = -1) = 0 = ϕ(t, x = 1). (4.1)
In this case, the corresponding L matrix at first order in space reads:

L D1 =          2 -1 0 • • • 0 -1 2 -1 . . . . . . 0 . . . . . . . . . 0 . . . . . . -1 2 -1 0 • • • 0 -1 2          ∈ R 2 h × 2 h ,
which is to be plugged into (2.3)-(2.4)-(2.7) in order to compute an approximate solution. Concerning the derivation of the "stopping pressure", one has to slightly modify (2.2): indeed, since the aforementioned matrix L D1 imposes that ϕ h (t, x = 0) = 0, there is no reason to suppose that (∂ x ϕ) n = 0. Thus, using that h(∂ x ϕ) n

1 2 = ϕ n 1
for Dirichlet conditions, (2.2) rewrites:

∀n ∈ N, w * 1 2 -z n 0 = 0 = w n 1 -z * 1 2 , exp(ϕ n 1 /λ 2 )(w * 1 2 + z n 0 ) = w n 1 + z * 1 2
.

In Fig. 4.1, we display the numerical approximations of the stationary solutions of (1.3) emanating from the simple set of discontinuous initial values with x ∈ [-1, 1]:

w h (t = 0, x) = 1 2 1135 + 15χ x∈[-0.65,0.45] , z h (t = 0, x) = 1 2 1135 -15χ x∈[-0.65,0.45] , ϕ h (t = 0, x) = 0, (4.2)
where χ A denotes the indicator function of the set A. The parameters in (1.1) are chosen as follows: D = α = β = 1 and λ = 10. The computational domain contains 256 grid points and the CFL number is fixed to 0.9, meaning that 1280∆t ≤ 1. Fig. 4.2 shows that one needs to iterate the whole numerical scheme up to t ≃ 10 in order to decrease the L 2 residues around the value of 10 -11 ; these residues are defined as:

∀t ∈ R + , R 2 (t) = w h (t + ∆t, .) -w h (t, .) L 2 + z h (t + ∆t, .) -z h (t, .) L 2 .
A key observation is that this corresponds to a macroscopic flux J h := λ(z hw h ) of the order of 10 can check on the top-left of Fig. 4.1 that w h = z h actually happens. The whole test-case is super-characteristic in the sense that the property λ ≥ |∂ x ϕ| does not generally hold, as can be seen by looking at the bottom-right of Fig. 4.1.

Remark 6. An interesting situation from a biological perspective is to impose a Dirichlet boundary condition on ρ and a Neumann condition on ϕ; this is usually referred to as the "killing boundary". This model can be recast in our well-balanced formulation too as we explain now: on the left side of the domain,

∀n ∈ N, w * 1 2 -z n 0 = w n 1 -z * 1 2 , exp(h(∂ x ϕ) n 1 2 /λ 2 )(w * 1 2 + z n 0 ) = w n 1 + z * 1 2
.

The Neumann condition on ϕ implies (∂ x ϕ) n

+ z n 0 = 0 = w n 1 + z * 1 2
, hence z * 1 2

= -w n 1 .

Numerical results with Neumann conditions on ϕ.

We switch now to the less stable case of the system (1.1) solved in x ∈ (0, 1) with the Neumann boundary condition on ϕ as prescribed in (2.5); these boundary conditions are implemented computationally through (2.2) and (2.8) inside the scheme (2.3)-(2.4)-(2.7).

Bifurcations of various types of equilibria for (1.1)

. By analogy with classical kinetic theory, we clearly distinguish between the 2 different types of timeasymptotic regimes which have been studied in e.g. [START_REF] Guarguaglini | Global stability of constant states and qualitative behavior of solutions to a one dimensional hyperbolic model of chemotaxis[END_REF]: Definition 1. For the system (1.1), or equivalently (1.3), inside a bounded domain, the following 2 types of equilibria are considered:

1. the "global (or simple) Maxwellians", which are asymptotic states satisfying J = 0 and not depending on the x variable. They read:

ρ ≡ C, ϕ ≡ C ′ , w = z ≡ C/2.

the "local

Maxwellians", which still satisfy J = 0, but do depend on the x variable. They are more general:

ρ = ρ(x), ϕ = φ(x), w = z = ρ(x)/2.
The figure 5.1 displays how the model bifurcates from the "global Maxwellians", which are the generic asymptotic profiles for small total mass (typically for small ρ h L 1 ) to more complex "local Maxwellians" when the mass increases. For small total mass, we are in the domain of application of the stability results shown in [START_REF] Guarguaglini | Global stability of constant states and qualitative behavior of solutions to a one dimensional hyperbolic model of chemotaxis[END_REF], and hopefully in the sub-characteristic regime because ∂ x ϕ is clearly related to the integral of ρ. The rigorous results of [START_REF] Gosse | Asymptotic-Preserving and Well-Balanced scheme for the 1D Cattaneo model of chemotaxis movement in both hyperbolic and diffusive regimes[END_REF] also apply to this class of asymptotic profiles. Beyond a critical value of the total mass, "local Maxwellians" also appear and they constitute the stable steady-states (the remaining global Maxwellians being stable only for symmetric perturbations). In general, these asymptotic states are reached after transient regimes which are super-characteristic. It has been shown in [START_REF] Guarguaglini | Global stability of constant states and qualitative behavior of solutions to a one dimensional hyperbolic model of chemotaxis[END_REF] that conventional upwind schemes fail to capture these local Maxwellians as they converge in time toward states endowed with a big mass flow, meaning that w = z. These local Maxwellians will be our main focus of interest in this section for testing our well-balanced method. More precisely, we have in mind to check numerically the bifurcation diagrams shown in [START_REF] Guarguaglini | Global stability of constant states and qualitative behavior of solutions to a one dimensional hyperbolic model of chemotaxis[END_REF] (and reproduced in Fig. 5.1) which consist in displaying the relation existing between the L 1 norm of 0 ≤ ρ h (t, .) (which is a constant with respect to the time t) and the value φ(x = 0) appearing at the steady-state of (1.1) once the parameters α, β and D are chosen. For instance, in the region corresponding to a small total mass of ρ h , the straight line on the left part of the diagram means that the asymptotic states are given by ρ h and ϕ h constant in x (the mass flux being zero). This is the domain of application of Theorem 2 and Corollary 1. When the total mass becomes higher, a bifurcation occurs and 3 states become admissible for the asymptotic values of ϕ h ; the value in the middle is unstable except when it is perturbed in a symmetric manner. We shall see in §5.3 that the other 2 admissible values for ϕ h correspond to the ones located at each side of the computational domain. If the total mass of ρ h is still augmented, the central branch bifurcates again into 2 other branches unstable except for symmetric perturbations, and a middle one which is completely unstable. In this region, numerical schemes will usually converge toward the extreme values of ϕ h (the smallest and the biggest ones, corresponding to the 2 preceding stable branches). All the numerical results (except in §5.1) are obtained with the following values:

1 h = 512, λ = 10, λ∆t = 0.9h.

Numerical stability of constant states.

In this first subsection, we aim at showing that the well-balanced scheme respects the stability of "global Maxwellian equilibria", that are the ones such that:

∀x ∈ [0, 1], w = z, ∂ x w = ∂ x z = 0, G(∂ x ϕ; w, z) = 0.
From a numerical viewpoint, these steady-states have the property that both the flux term and the source term vanish individually; this is clearly a particular subclass of the general steady-states of (1.3) which satisfy only to:

∀x ∈ [0, 1], w = z, ∂ x w = ∂ x z = 0, 2λ 2 ∂ x z = G(∂ x ϕ; w, z).
Global Maxwellian asymptotics are typically the only ones appearing for little initial mass, leading to non-stiff problems on which standard upwind schemes can give good results [START_REF] Guarguaglini | Global stability of constant states and qualitative behavior of solutions to a one dimensional hyperbolic model of chemotaxis[END_REF][START_REF] Natalini | An asymptotic high order mass-preserving scheme for a hyperbolic model of chemotaxis[END_REF]; for instance, one sees on Fig. 5.1 that in the case α = β = 1, prescribing 1 0 ρ h (t, x)dx ≤ 1100, or 1 0 ρ h (t, x)dx ≤ 100 for α = β = 50 leads to constant steady-states. The outcome of the well-balanced scheme with 1 0 ρ h (t, x)dx = 65 and α = β = 50 is displayed on Fig. 5.2 with a small number of grid points: 128h = 1 and t = 3. This corresponds to an illustration of Theorem 2. Another interesting test-case consists in investigating the stability of constant states with respect to symmetric perturbations with α = β = 1 while prescribing 1 0 ρ h (t, x)dx = 1135 and the same number of grid points. It turns out that the dynamic is similar, except that a very small destabilizing process is ignited after t ≃ 2 (around 2840 iterations with a CFL number of 0.9): see Fig. 5.3 for numerical results at much later time. 5.3. First choice: α = β = 1 and 1 0 ρ h (t, x)dx = 1135. This test-case appears in both [START_REF] Guarguaglini | Global stability of constant states and qualitative behavior of solutions to a one dimensional hyperbolic model of chemotaxis[END_REF][START_REF] Natalini | An asymptotic high order mass-preserving scheme for a hyperbolic model of chemotaxis[END_REF]; it is precisely with this choice of parameters that the spurious behavior of standard upwind numerical schemes is put at forefront. Indeed, in a way which is somewhat similar to the early experiments in [START_REF] Cargo | Un schéma équilibre adapté au modèle d'atmosphère avec termes de gravité[END_REF], it is found that conventional discretizations converge asymptotically toward solutions endowed with quite a big macroscopic flux J despite the fact it should vanish as t diverges. In particular, this means that these numerical steady-states are not Maxwellian because for a 2-velocities kinetic model, a Maxwellian solution is such that w = z and this implies trivially that J = λ(zw) = 0. On the contrary, one can see in Fig. 5.4 that the well-balanced scheme (2.3)-(2.4)-(2.7) endowed with second-order boundary conditions (2.8) stabilizes nicely onto a steady-state for which J h = λ(z hw h ) is of the order of 10 -5 at time t = 30: compare especially with Fig. 4 in [START_REF] Guarguaglini | Global stability of constant states and qualitative behavior of solutions to a one dimensional hyperbolic model of chemotaxis[END_REF] where a conventional scheme is shown to deliver a steady-state for which the macroscopic flux is of the order of 10 L. Gosse w h and z h are nearly Maxwellian as can be checked visually on the top left figure; moreover, the value of ϕ h (t = 30, x = 0) is compatible with the theoretical ones coming from the bifurcation diagram of [START_REF] Guarguaglini | Global stability of constant states and qualitative behavior of solutions to a one dimensional hyperbolic model of chemotaxis[END_REF], that is, around 1220. The residues' decay displayed in Fig. 5.5 is interesting because it neatly shows 3 distinct moments occurring during the stabilization process:

1. t ∈ (0, 2) roughly: the asymmetric perturbations are smoothed out progressively. The quantities ρ h and ϕ h remain very close to constant and macroscopic fluxes |J| are of the order of 50 around t = 2. Residues fall below 10 -5 but their evolution show high frequency oscillations (which seem to appear every time the perturbation touches one edge of the computational domain). 2. t ∈ (2, 13) roughly: The scheme didn't succeed in stabilizing the fluxes with ρ h and ϕ h nearly constants. This is consistent with Fig. 5.1 where one sees that 1 0 ρ h (t, x)dx = 1135 corresponds to a value located slightly beyond the first bifurcation pitchfork. Hence a new dynamic takes place with macroscopic fluxes J starting to increase together with residues, this time showing no high frequency oscillations. This is the moment during which all the cells agglutinate on the left side of the computational domain. Around t = 13, one has correct values, namely ρ h (t, x = 0) ≃ 2220 and ϕ h (t, x = 0) ≃ 1210. 3. t ∈ (13, 30): residues decay monotonically together with macroscopic fluxes J during the time the numerical scheme stabilizes globally. Values in x = 0 are only marginally adjusted and the kinetic densities w h , z h can be considered Maxwellian inside the whole computational domain.

5.4. Second choice: α = β = 1 and 1 0 ρ h (t, x)dx = 4100. This is a more difficult test-case as the mass initially prescribed goes beyond the second bifurcation point on Fig. 5.1. Numerical results generated by the well-balanced scheme are shown on Fig. 5.6 around time t = 35. The evolution is very different compared to the preceding test-case since high frequency oscillations appear for t ≥ 0.5 and show up until global stabilization. Numerical macroscopic fluxes J h stabilize around a value of the order of -5.10 -5 which can be considered as satisfying for such a problem involving a strong mass concentration around x = 0. Interestingly, one can observe that the values of ϕ h (t = 35, .) go roughly from 5260 down to 3550: this matches the two values of the red curve at the abscissa 4100 on the bifurcation diagram of Fig. 5.1. Hence this strongly suggests that beyond the first bifurcation point, the two branches of the pitchfork curve indicate the (stable) values ϕ(x = 0) and ϕ(x = 1). This sheds some light onto the meaning of both diagrams displayed in Fig. 5.1 for some given values of the parameters α, β and D. The curve of the residues shows rather high frequency oscillations: they seem to be stable and they survive as one passes from the implicit method (2.7) on ϕ h to the Crank-Nicholson scheme (see Remark 4). They also appear for different choices of the space step h, even if finer grids make them appear later and smaller.

Remark 7. We wish to propose in this remark an explanation for the appearance of high-frequency oscillations in the curve of the evolution of residues with respect to time based on several convenient numerical tests. They consist in simulating the well-balanced scheme for (1.3) coupled to the Crank-Nicholson discretization for the diffusion equation on ϕ with increasing total mass and D = α = β = 1. More precisely, we let this marching process go up to t = 2.5 and indicate in red the instants t n at which the numerical approximations w h (t , .) or z h (t , .) cease to be nonnegative. On Fig. 5 The curve deforms itself continuously as the first one is very similar to the one shown on Fig. 5.5 and the last one looks like Fig. 5.7. One sees that the increase of the total mass allows for more and more negative values for the kinetic densities w h , z h : this is a consequence of the concentration process in the vicinity of the borders of the computational domain whose strength increases with the initial total mass of ρ h . This accumulation process in ρ passes onto ϕ thanks to the nature of the diffusion equation thus some very big values of ∂ x ϕ appear. At this point, we touch 2 intrinsic limits of the present approach:

• the approximation of ∂ x ϕ as a constant in the non-conservative jump relations, already discussed in Remark 2. An exponential fitting (∂ x ϕ) j+ 1 2 = γϕ j exp(γ/2)/h, γ = log(ϕ j+1 /ϕ j ) gives only a marginal improvement.

• the dicretization of the diffusion equation on ϕ is not well-balanced, as already noted in [START_REF] Gosse | Asymptotic-Preserving and Well-Balanced scheme for the 1D Cattaneo model of chemotaxis movement in both hyperbolic and diffusive regimes[END_REF], and this may create issues in stiff regimes. A remedy may be contained in the recent paper [START_REF] Ten Thije | The Finite Volume-Complete Flux Scheme for Advection-Diffusion-Reaction Equations[END_REF], the "complete flux scheme". Let us also recall that existence and stability results given in [START_REF] Hillen | Hyperbolic models for chemotaxis in 1-D, Nonlinear Analysis[END_REF][START_REF] Hyung | Global existence of classical solutions for a hyperbolic chemotaxis model and its parabolic limit[END_REF] are for initial data close to equilibrium. We are in a very different situation here. 5.5. Third choice: α = β = 50 and 1 0 ρ h (t, x)dx = 125. This last test-case corresponds to a stabilization in time onto a non-constant steady-state after the first bifurcation point. The red curve indicating stable values for ϕ at each edge of the computational domain is a pitchfork with the inferior branch being very close to zero and the superior one which increases quite fast with respect to the initial mass. This hints that ρ is likely to concentrate strongly close to the borders. This is indeed what happens on Fig. 5.9 which display the results generated by the well-balanced scheme around time t = 30; as before, ϕ h (t, .) seems to connect both values of the pitchfork curve going this time from nearly zero up to 680. The macroscopic fluxes J h are of the order of 10 To establish strong convergence of w h and z h toward the unique solution of (1.3), we need a bound on the total variation of w h (t, .) and z h (t, .). Lemma A.3. Let β ≥ 0, ϕ 0 ∈ W 1,∞ (R) and w 0 , z 0 ∈ L 1 ∩ BV (R); under both the CFL condition λ∆t ≤ h and the subcharacteristic restriction λ ≥ sup j,n |(∂ x ϕ) n j+ 1 2 |, the following BV-bounds hold for any n ∈ N: T V w h (n∆t, .) + T V z h (n∆t, .) ≤ T V (w 0 ) + T V (z 0 ) + L λ n∆t w 0 L 1 (R) + z 0 L 1 (R) ,

(A.6)
where L is the Lipschitz constant of 1-A 1+B , 1-B 1+B depending on the values of (∂ x ϕ) n j+ 1 2 . The estimates given in (A.5) and (A.6) allow to establish convergence as h → 0 in the sub-characteristic regime by standard methods, see [START_REF] Gosse | Asymptotic-Preserving and Well-Balanced scheme for the 1D Cattaneo model of chemotaxis movement in both hyperbolic and diffusive regimes[END_REF] for details. Both these estimates fail when we enter the super-characteristic region λ < sup j,n |(∂ x ϕ) n j+ 1 2

|.

A.2. Diffusive scaling: Asymptotic-Preserving. Concerning the asymptoticpreserving property of the scheme (A.4), let us introduce a diffusive scaling related to a small parameter 0 < ε ≪ 1: x → x/ε, t → t/ε 2 . Based on ideas of [START_REF] Gosse | An asymptotic-preserving well-balanced scheme for the hyperbolic heat equations[END_REF][START_REF] Gosse | Space localization and well-balanced scheme for discrete kinetic models in diffusive regimes[END_REF], it has been proposed in [START_REF] Gosse | Asymptotic-Preserving and Well-Balanced scheme for the 1D Cattaneo model of chemotaxis movement in both hyperbolic and diffusive regimes[END_REF] to rewrite (A.4) as follows: w n jz n j-1 -

λ∆tC n ε,j-1 2 εh 1+B n ε,j-1 2
w n j .

(A.7)

The rescaled parameters read:

A ε = 1 -h ελ exp(h∂ x ϕ/ελ 2 ), B ε = 1 + h ελ exp(h∂ x ϕ/ελ 2 ).

It is possible (even if somewhat tedious) to establish L p and BV bounds of the type (A.5) and (A.6) uniform in ε for (A.7) under the usual parabolic CFL condition, ∆t = O(h 2 ). Moreover, a control on the deviation with respect to the Maxwellian distribution can be established: Lemma A.4. If w 0 -z 0 L 1 (R) = O(ε) (well-prepared initial data) and for ε small enough, there holds for any t > 0: w h (t, .)z h (t, .) L 1 (R) = O(ε). (A.8) This lemma is essential in order to prove rigorously the AP property of the scheme; indeed, by first letting ε → 0, one gets:

Theorem A.5. Under the assumptions of Lemma A.4, for any n ∈ N:

ρ n+1 j -ρ n j ∆t + (∂ x ϕ) n j+ 1 2 ρ n j+1 -(∂ x ϕ) n j-1 2 ρ n j h = λ 2 ρ n j+1 -2ρ n j + ρ n j-1 h 2 + O(ε).
In particular, the sequences w h , z h are relatively compact in L 1 loc (R + * × R) as ε → 0 with h > 0 fixed and the remaining term in O(ε) converges to zero in L

Fig. 2 . 1 .

 21 Fig. 2.1. Control cell for the well-balanced scheme with fictious values w 0 , z 0 in x = 0.

Fig. 3 . 1 .

 31 Fig. 3.1. Perturbation ρ 0 (x): sin(4π|x -1 2 |) (left) and cos(8πx) (right).

Fig. 4 . 2 .

 42 Fig. 4.2. Decay of L 2 residues as a function of time.

1 2 = 0 .

 120 The Dirichlet condition on ρ is expressed with w * 1 2

Fig. 5 . 1 .

 51 Fig. 5.1. Bifurcation diagrams (reproduced from [24]): α = β = 1(left) and α = β = 50 (right).

FluxFig. 5 . 2 .

 52 Fig. 5.2. Global Maxwellian for α = β = 50, t = 2: J h (left) and residues (right).

Flux

Fig. 5 . 3 .

 53 Fig. 5.3. Global Maxwellian for α = β = 1 around t = 15.

Fig. 5 . 5 .

 55 Fig. 5.5. Evolution of L 2 residues as a function of time.

Fig. 5 . 7 .Fig. 5 . 8 .

 5758 Fig. 5.7. Decay of L 2 residues as a function of time.

Fig. 5 . 9 .Fig. 5 . 10 .

 59510 Fig. 5.9. Cell density ρ h , chemoattractant ϕ h (left) and mass flow rate J h , (∂xϕ) h (right).

Fig. B. 2 .

 2 Fig. B.2. Evolution of L 2 residues as a function of time.

  -7 and this is a desirable feature as one expects the numerical scheme to stabilize in time onto a Maxwellian distribution which, in the present context of a discrete 2-velocities kinetic model, is characterized by the condition w h = z h . One

												2.0e-007
												1.5e-007
												1.0e-007
												5.0e-008
												0.0e+000
												-5.0e-008
												-1.0e-007
												-1.5e-007
												-2.0e-007
	-1.0	-0.8	-0.6	-0.4	-0.2	0.0	0.2	0.4	0.6	0.8	1.0	-1.0	-0.8	-0.6	-0.4	-0.2	0.0	0.2	0.4	0.6	0.8	1.0
		Density (w^h+z^h)				Initial density					Flux
		z^h moving right
		w^h moving left
												1000
												800
												600
												400
												200
												0
												-200
												-400
												-600
												-800
												-1000
	-1.0	-0.8	-0.6	-0.4	-0.2	0.0	0.2	0.4	0.6	0.8	1.0	-1.0	-0.8	-0.6	-0.4	-0.2	0.0	0.2	0.4	0.6	0.8	1.0
		Chemotractant										Space derivative
		Initial chemotractant									Initial space derivative

Fig. 4.1. Cell density ρ h , chemoattractant ϕ h (left) and mass flow rate J h , (∂xϕ) h (right).

  2 . It is therefore clear that our steady-state kinetic densities Cell density ρ h , chemoattractant ϕ h (left) and mass flow rate J h , (∂xϕ) h (right).

											0.0e+000
											-5.0e-006
											-1.0e-005
											-1.5e-005
											-2.0e-005
											-2.5e-005
											-3.0e-005
											-3.5e-005
											-4.0e-005
											-4.5e-005
											-5.0e-005
	0.0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1.0	0.0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1.0
		Density (w^h+z^h)				Initial density					Flux
		z^h moving right
		w^h moving left
											50
											0
											-50
											-100
											-150
											-200
											-250
	0.0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1.0	0.0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1.0
		Chemotractant										Space derivative
		Initial chemotractant									Initial space derivative
	Fig. 5.4.

  -7 and this is equivalent to say that the kinetic densities w h , z h are

	10												6e-007
													5e-007
	10
													4e-007
	10												3e-007
													2e-007
	10
													1e-007
	10	0.0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1.0	0e+000	0.0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1.0
			Density (w^h+z^h)				Initial density						Flux
			z^h moving right
			w^h moving left
													3000
													2500
													2000
													1500
													1000
													500
													0
													-500
		0.0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1.0		0.0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1.0
			Chemotractant											Space derivative
			Initial chemotractant										Initial space derivative

  1 .

											1.2e-009
											1.0e-009
											8.0e-010
											6.0e-010
											4.0e-010
											2.0e-010
											0.0e+000
											-2.0e-010
											-4.0e-010
											-6.0e-010
	0.0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1.0	0.0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1.0
		Density (w^h+z^h)				Initial density					Flux
		z^h moving right
		w^h moving left
											50
											0
											-50
											-100
											-150
											-200
											-250
											-300
	0.0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1.0	0.0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1.0
		Chemotractant										Space derivative
		Initial chemotractant									Initial space derivative

Fig. B.1. Cell density ρ h , chemoattractant ϕ h (left) and mass flow rate J h , (∂xϕ) h (right).

practically Maxwellian. The evolution of the residues on Fig. 5.10 is very reminiscent of the one shown in Fig. 5.7. The test-case is super-characteristic since |∂ x ϕ| ≫ λ.

Conclusion and outlook.

We addressed 2 main topics in this paper:

1. The numerical experiments of [START_REF] Guarguaglini | Global stability of constant states and qualitative behavior of solutions to a one dimensional hyperbolic model of chemotaxis[END_REF] show that conventional upwind schemes display a spurious behavior when it comes to approximate the long-time asymptotics of the system (1.1) in super-characteristic relaxation regime for which |∂ x ϕ| > λ. A new type of discretization (so-called Asymptotic High-Order schemes [START_REF] Aregba-Driollet | Asymptotic high-order schemes for 2 × 2 dissipative hyperbolic systems[END_REF][START_REF] Natalini | An asymptotic high order mass-preserving scheme for a hyperbolic model of chemotaxis[END_REF]) has been set up specifically to correct this deceptive feature. Here, we explained that this defect is mainly a well-balancing problem which is reminiscent of the spurious behavior of standard numerical schemes on shallow water equations with topography or some atmospheric modelings [START_REF] Cargo | Un schéma équilibre adapté au modèle d'atmosphère avec termes de gravité[END_REF], especially when it comes to decaying toward non-constant steady-states. A well-balanced and asymptotic-preserving scheme for (1.1) has been studied in [START_REF] Gosse | Asymptotic-Preserving and Well-Balanced scheme for the 1D Cattaneo model of chemotaxis movement in both hyperbolic and diffusive regimes[END_REF] in sub-characteristic regime; here we show numerically that it allows to simulate efficiently the correct long-time asymptotics prescribed in [START_REF] Guarguaglini | Global stability of constant states and qualitative behavior of solutions to a one dimensional hyperbolic model of chemotaxis[END_REF] with both constant and non-constant steady-states. 2. We also raised the fact that, since (1.1) is essentially a discrete velocity kinetic model, the spurious stationary states of conventional schemes with big mass flows inside the computational domain correspond to a strong deviation away from the Maxwellian distribution. Thus one can think that well-balanced discretizations may be necessary in the context of the numerical simulation of general kinetic problems for which a convergence onto the Maxwellian distribution is expected for large times. A first step in this direction has been taken for 1-D radiative transfer models in [START_REF] Gosse | Transient radiative tranfer in the grey case: well-balanced and asymptotic-preserving schemes built on Case's elementary solutions[END_REF] by taking advantage of the Case's eigen-decomposition in "elementary solutions" (sometimes referred to as Caseology [START_REF] Case | Elementary solutions of the transport equation and their applications[END_REF][START_REF] Pomraning | Weakly) three-dimensional Caseology[END_REF]). Let us quickly present how one can extend these results toward seemingly more involved chemotaxis models: following [START_REF] Othmer | The diffusion limit of transport equations II: Chemotaxis equations[END_REF], the Cattaneo system (1.3) is a 2-velocity restriction of more a general kinetic equation which reads,

This is a so-called "run and tumble" model: the parameter m > 0 measures the sensibility to ∂ x ϕ during the "tumbling phases" where the cell reorients itself in a direction of increasing chemoattractant concentration. Being able to build a well-balanced and asymptotic-preserving scheme for (6.1) means that one can compute easily the solutions of the forward-backward problem for its steady-state equations: for m = 0, this is done in [START_REF] Gosse | Transient radiative tranfer in the grey case: well-balanced and asymptotic-preserving schemes built on Case's elementary solutions[END_REF] by exploiting the seminal paper of K. Case [START_REF] Case | Elementary solutions of the transport equation and their applications[END_REF]. The handling of the flux term for m > 0 proceeds by observing that, assuming again that ∂ x ϕ is a constant at each interface, the integral 1 2λ λ -λ mξ ′ f (t, x, ξ ′ )dξ ′ equals a constant, J m . Thus its effect reduces to a translation of J m ξ∂ x ϕ with respect to the solutions computed in [START_REF] Gosse | Transient radiative tranfer in the grey case: well-balanced and asymptotic-preserving schemes built on Case's elementary solutions[END_REF]. This trick has been used in a different context in [START_REF] Gosse | Well-balanced schemes using elementary solutions for linear models of the Boltzmann equation in one space dimension[END_REF].

Appendix A. The sub-characteristic WB/AP scheme of [START_REF] Gosse | Asymptotic-Preserving and Well-Balanced scheme for the 1D Cattaneo model of chemotaxis movement in both hyperbolic and diffusive regimes[END_REF].

A.1. Hyperbolic scaling: Well-balanced. In order to solve the Riemann problem for (1.4), we must derive the jump relations ruling the non-conservative product. Let us denote by w, z the microscopic profiles which are shrunk inside the discontinuity: they satisfy a forward-backward problem for the stationary equations of (1.3) which reads, for x ∈ [0, h],

with inflow boundary data. Its solution reads, with obvious notation:

At this microscopic scale, the quantity ∂ x ϕ is a constant. Hence we have the result: Lemma A.1. For any h > 0, the stationary equations of (1.3) yield the following jump relations across the zero-waves of (1.4) located in (j - 1 2 )h, j ∈ Z:

with the notation:

In particular, the following important relation holds: (flux conservation)

It is interesting to observe that a simple linearization of exponentials gives:

These linearizations will be useful for establishing consistency as h → 0 since |∂ x ϕ| remains bounded for D, β ≥ 0 are big enough. In practice, a first order divided difference computed at each interface x = (j -1 2 )h and t = n∆t will appear in place of ∂ x ϕ; there are no discontinuities involved in the discretization of ϕ, obviously. Once the jump relations (A.2) are available, there is no obstacle in the derivation of the scheme (2.3). Clearly, it can be rewritten as a more standard upwind scheme:

)z n j-1 .

(A.4)

We stress that since A and B are computed at the borders of each control cell, the term ∂ x ϕ is very well defined because odd derivatives "live on the staggered grid". Lemma A.2. Let w 0 and z 0 belong to L p (R), 1 ≤ p ≤ ∞; under both the CFL condition λ∆t ≤ h and the subcharacteristic restriction λ ≥ |∂ x ϕ|, there holds:

Appendix B. Improvements coming from the Crank-Nicholson scheme. In Remark 4, a more sophisticated discretization for ϕ is briefly presented. Here, we want to illustrate the behavior of our well-balanced scheme (2.3)-(2.4) for (1.3) when it is completed by the Crank-Nicholson method. We selected initial data of the type (4.2), but with a stronger discontinuity:

ϕ h (t = 0, x) = 0.

The space and time steps correspond to a quite coarse grid: ): an interesting improvement of the Crank-Nicholson scheme is that residues decay more quickly and at time t = 30, they stall around 10 -12 . However, we still observe the decay decomposed into 3 different phases as explained in §5.2.