Bruno Després
email: despres@ann.jussieu.fr

Christophe Buet

The structure of well-balanced schemes

Keywords: well balanced schemes, Friedrichs systems, conservative formulation, finite volume schemes. 2000 MSC: 65J10, 65N06, 65N99

Introduction

Our interest in this work is the mathematical structure of Friedrichs systems with linear relaxation, having in mind the definition of well-balanced schemes which is a hot topic nowadays [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF][START_REF] Gosse | An asymptotic-preserving well-balanced scheme for the hyperbolic heat equations[END_REF][START_REF] Jin | Asymptotic preserving (AP) schemes for multiscale kinetic and hyperbolic equations: a review[END_REF][START_REF] Buet | Design of asymptotic preserving schemes for the hyperbolic heat equation on unstructured meshes Numerish Mathematik[END_REF][START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws, Exponential-Fit, Well-Balanced and Asymptotic-Preserving[END_REF]. The particular case of Friedrichs systems of large size is challenging for numerical methods, and quite interesting since such large size systems are commonly encountered in the approximation of kinetic equations by moment methods [START_REF] Gosse | Well-balanced schemes using elementary solutions for linear models of the Boltzmann equation in one space dimension[END_REF][START_REF] Buet | Asymptotic Preserving Schemes on Distorted Meshes for Friedrichs Systems with Stiff Relaxation: Application to Angular Models in Linear Transport[END_REF]. In this work we will consider the S n model for radiation or neutrons propagation, which is a paradigm for approximations of kinetic equations [START_REF] Chandraseckhar | Radiative transfer[END_REF]. Non linear extensions such as well-balanced schemes for shallow water equations or Euler equations are discussed in [START_REF] Greenberg | A well-balanced scheme for the numerical processing of source terms in hyperbolic equations[END_REF][START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF][START_REF] Muller | Well-balanced simulation of geophysical flows via the shallow water equations with bottom topography: consistency and numerical computations[END_REF][START_REF] Ortleb | A well-balanced DG scheme with unconditionally positive implicit time integration[END_REF][START_REF] Coquel | Asymptotic preserving scheme for Euler system with large friction[END_REF].

Our generic model problem is a linear system with relaxation in two dimensions

∂ t U + ∂ x (A(x)U) + ∂ y (B(x)U) = -R(x)U, U(t, x) ∈ R n , x = (x, y) ∈ R 2 , (1)
where the unknown is the function U(t, x). The matrices A(x), B(x), R(x) ∈ R n×n may be functions of the space variable, even if they will take as constant in most of this work. In all applications we have in mind the symmetric part of the relaxation matrix on the right hand side is non negative in the sense that (V, RV) ≥ 0 for all V ∈ R n , that is

R + R t ≥ 0. ( 2
)
These matrices are symmetric for a Friedrichs system, that is A = A t , B = B t . The size n can be arbitrarily large. A more specific example which serves as a guideline in this text is the hyperbolic heat equation with n = 2

∂ t p + ∂ x u = 0, ∂ t u + ∂ x p = -σu. (3)
Here U = (p, u) t , A = 0 1 1 0 and R = 0 0 0 σ . Such linear systems are representative of the linearization of non linear systems with relaxation. Indeed consider the p-system with friction σ ≥ 0

∂ t τ -∂ x u = 0, ∂ t u + ∂ x p(τ) = -σu.
Linearization τ = τ 0 + τ 1 + . . . and u = u 0 + u 1 + . . . around an equilibrium ∂ x τ 0 = ∂ t τ 0 = u 0 = 0 yields

∂ t τ 1 -∂ x u 1 = 0, ∂ t u 1 -c 2 0 ∂ x τ 1 = -σu 1 ,
where c 0 > 0 is the speed of sound. This linear system can be rewritten under the form of the Friedrichs system with linear relaxation using the symetrized variables U = (c 0 τ 1 , u 1 ) t and the matrices A = -c 0 0 1 1 0 and R = σ 0 0 0 1 .

Remark 1. Adding gravity as in ∂ t τ -∂ x u = 0 with ∂ t u + ∂ x p(τ) = g -σu changes the structure of the linearized equations

∂ t τ 1 -∂ x u 1 = 0, ∂ t u 1 -∂ x c 0 (x) 2 τ 1 = -σu 1
where the speed sound is no more constant c 0 = c 0 (x) = -p (p -1 (xg)). For simplicity we will not consider this case hereafter, but we mention in our conclusion how the ideas developed in this work can be adapted without major difficulties.

A standard method for the discretization of problems like (1) relies on the splitting method: first, use a Finite Volume technique for the approximation of the homogeneous equation ∂ t U + ∂ x (A(x)U) + ∂ y (B(x)U) = 0, without right hand side; second, add the right hand side solving ∂ t U = -R(x)U. Even if efficient for most cases, a splitting method does not respect by construction the stationary states defined by

U = x → U(x); ∂ x (A(x)U) + ∂ y (B(x)U) = -R(x)U .
It has been observed in many cases that splitting methods may generate large numerical errors and unphysical oscillations, which therefore must be controlled in a way or another: a recent theoretical contribution on this topic is in [START_REF] Amadori | Error estimates for well-balanced and time-split schemes on a damped semilinear wave equation[END_REF]. In this direction the so called well-balanced techniques aim at combining Finite Volume techniques with the knowledge of the right-hand side so as to obtain new schemes which are exact for initial data in U. We refer to [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF][START_REF] Gosse | An asymptotic-preserving well-balanced scheme for the hyperbolic heat equations[END_REF][START_REF] Jin | Asymptotic preserving (AP) schemes for multiscale kinetic and hyperbolic equations: a review[END_REF][START_REF] Buet | Design of asymptotic preserving schemes for the hyperbolic heat equation on unstructured meshes Numerish Mathematik[END_REF] for examples. A recent and comprehensive state-of-the-art is to be found in [START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws, Exponential-Fit, Well-Balanced and Asymptotic-Preserving[END_REF]. Moreover a efficient well-balanced scheme is often the starting point of a rigorous asymptotic preserving scheme, a topic that will not be developed in this work but for which we refer the reader to [START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws, Exponential-Fit, Well-Balanced and Asymptotic-Preserving[END_REF][START_REF] Jin | Asymptotic preserving (AP) schemes for multiscale kinetic and hyperbolic equations: a review[END_REF][START_REF] Buet | Design of asymptotic preserving schemes for the hyperbolic heat equation on unstructured meshes Numerish Mathematik[END_REF]. All these new methods are quite complex to construct and to analyze, and therefore are difficult to understand. This is why more mathematical analysis is needed to explore the fundamental structures of these techniques.

The basis of the method used hereafter tries to explore such a structure: it is reminiscent in some sense of the seminal work [START_REF] Cargo | Un schéma équilibre adapté au modèle d'atmosphère avec termes de gravité. (in French) [A well-balanced scheme for a model of an atmosphere with gravity[END_REF], since we begin to modify the equation in a conservative way so that usual Finite Volume schemes can be directly used to obtain methods which are exact for initial data in U. For this purpose we will use the dual equation

∂ t V + A t (x)∂ x V + B t (x)∂ y V = R t (x)V, V(t, x) ∈ R n . (4)
A fundamental property is

∂ t (U, V) + ∂ x (A(x)U, V) + ∂ y (B(x)U, V) = 0 (5)
for all U solutions of the primal equation ( 1) and V solutions of the dual equation [START_REF] Bouchut | Nonlinear stability of finite volume methods for hyperbolic conservation laws, and well-balanced schemes for sources[END_REF]. The remarkable fact is that identity [START_REF] Buet | Design of asymptotic preserving schemes for the hyperbolic heat equation on unstructured meshes Numerish Mathematik[END_REF] is in conservation form. It means that if one has enough knowledge of the solutions of the dual equation (4)

V = x → V(x); A(x) t ∂ x V + B(x) t ∂ y V = R t (x)V ,
then it is possible to replace the non conservative primal equation (1) by the conservative identity [START_REF] Buet | Design of asymptotic preserving schemes for the hyperbolic heat equation on unstructured meshes Numerish Mathematik[END_REF]. That is instead of analyzing the primal set U, we put the emphasis on the dual set V of stationnary states of (4) which are now test functions. This is the basis of this work. Concentrating of Finite Volume techniques and in view of formula [START_REF] Buet | Design of asymptotic preserving schemes for the hyperbolic heat equation on unstructured meshes Numerish Mathematik[END_REF], it is possible to conjecture that any well balanced Finite Volume solver for the primal equation (1) can be recast as a standard Finite Volume solver for the conservative formulation (5) (more precisely written as equation [START_REF] Chandraseckhar | Radiative transfer[END_REF] in the core of the paper). In what follows we more modestly discretize directly [START_REF] Buet | Design of asymptotic preserving schemes for the hyperbolic heat equation on unstructured meshes Numerish Mathematik[END_REF] with usual conservative finite volume solvers, and deduce well-balanced solvers for the non conservative primal formulation [START_REF] Amadori | Error estimates for well-balanced and time-split schemes on a damped semilinear wave equation[END_REF]. We study two families of solvers, which are natural in our context. For the second family called the two-states solver, we show that it corresponds to a well balanced finite volume scheme based on space localization of the source term at the interfaces, see [START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws, Exponential-Fit, Well-Balanced and Asymptotic-Preserving[END_REF]. But in our approach there is absolutely no need of the localization method. Another asset of this method is the possibility to treat general meshes in higher dimensions. Notice that this approach is not restricted to Finite Volume techniques since the starting point is a modification of the equation.

This work is organized as follows. In section 2 we detail in one dimension the structure [START_REF] Buet | Design of asymptotic preserving schemes for the hyperbolic heat equation on unstructured meshes Numerish Mathematik[END_REF] and propose a new conservative formulation of the initial non conservative problem [START_REF] Amadori | Error estimates for well-balanced and time-split schemes on a damped semilinear wave equation[END_REF]. Section 3 is dedicated to the discretization (still in dimension one) by means on standard Finite Volume techniques applied to the conservative system. We detail two different solvers. The main theorem about the well-posedness of the second family of solvers is given in this section. The example of the hyperbolic heat equation shows the second family is the same as the Gosse-Toscani scheme. The application of these ideas is performed in detail for the S n model in section 4. The singular case det(A) = 0 is briefly detailed in section 5 and a possible multiD extension in section 6. Some conclusions are drawn in section 7.

For the simplicity of the presentation (and only for that reason), we will consider that A and B are now constant in space matrices. Other matrices will be constructed which are non constant in space.

The interested reader can find many numerical illustrations which ilustrate the theory of this paper in [START_REF] Buet | Design of asymptotic preserving schemes for the hyperbolic heat equation on unstructured meshes Numerish Mathematik[END_REF][START_REF] Coquel | Asymptotic preserving scheme for Euler system with large friction[END_REF][START_REF] Gosse | Well-balanced schemes using elementary solutions for linear models of the Boltzmann equation in one space dimension[END_REF][START_REF] Gosse | Transient radiative transfer in the grey case: Well-balanced and asymptotic-preserving schemes built on Case's elementary solutions[END_REF][START_REF] Greenberg | A well-balanced scheme for the numerical processing of source terms in hyperbolic equations[END_REF][START_REF] Jin | Asymptotic preserving (AP) schemes for multiscale kinetic and hyperbolic equations: a review[END_REF] and references therein. We particularily quote [START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws, Exponential-Fit, Well-Balanced and Asymptotic-Preserving[END_REF] for the variety of the models and of the schemes, and [START_REF] Buet | Asymptotic Preserving Schemes on Distorted Meshes for Friedrichs Systems with Stiff Relaxation: Application to Angular Models in Linear Transport[END_REF] where numerical tests are performed for Friedrichs system (including the Gosse-Toscani scheme which is a special case of our two-states solver below).

Conservative formulation in one dimension

In this section the matrices may be non symmetric. The set of stationary solutions of the dual equation is

V = x → V(x); A t ∂ x V = R t V . (6)
This is a vectorial space of dimension p with 0 ≤ p ≤ n. Let us denote a basis as V 1 (x), . . . , V p (x), so that V = Span V 1 (x), . . . , V p (x) .

Proposition 1 (Evident). Solutions of the primal equation (1) satisfy p linearly independent conservation laws

∂ t (U, V i ) + ∂ x (AU, V i ) = 0, i ≤ p. (7)
One can define

α i = (U, V i ), the vector α =           α 1 . . . α p           ∈ R p and the matrix P(x) =           V 1 (x) t . . . V p (x) t          
∈ R p×n so that (7) can be rewritten in a more compact form

∂ t α + ∂ x (P(x)AU) = 0. (8)
An even more compact formulation is possible if A is non singular.

Proposition 2. Assume A is non singular, that is det(A) 0. Then p = n and the matrix P can be represented with the matrix exponential

P(x) = e RA -1 x . ( 9
)
Moreover the system (8) can rewritten as

∂ t α + ∂ x (Q(x)α) = 0, ( 10
)
where the change of unknown is α = P(x)U ⇐⇒ U = P -1 (x)α and the matrix Q

(x) = P(x)AP -1 (x) is similar to A. Proof. Any V ∈ V is solution of ∂ x V = A -t R t V, so can be represented as V(x) = e A -t R t x W where W ∈ R n is arbitrary. So p = n. Moreover (5) yields ∂ t U, e A -t R t x W + ∂ x AU, e A -t R t x W = 0. Since W is arbitrary, it shows that ∂ t e RA -1 x U + ∂ x e RA -1 x AU = 0.
The last part of the claim is immediate. The proof is ended.

Remark 2. Another proof of the identity [START_REF] Chandraseckhar | Radiative transfer[END_REF] uses an idea commonly used in the numerical theory of Fokker-Planck equations. It is based on the identity

A∂ x U + RU = Ae -A -1 Rx ∂ x e A -1 Rx U = e -RA -1 x A∂ x e A -1 Rx U .
It yields the conservative equation ∂ t e RA -1 x U + A∂ x e A -1 Rx U = 0 which is exactly equal to [START_REF] Chandraseckhar | Radiative transfer[END_REF]. This is restricted of course to non singular matrices A. We will not develop this angle of attack since it is less clear how to generalize it to the singular case where det(A) = 0 and to the multidimensional case, as it will be done at the end of this work. Moreover this approach completely misses a central idea in this work which is that the matrix P(x) of the change of basis of the linear transformation is a fundamental object that needs specific examination.

Remark 3. The hypothesis that A is non singular appears at many places in this work. If the linear model derives from an underlying non linear model, it is equivalent to state that this underlying non-linear model is non resonant. The more difficult resonant case is briefly studied at the end of this work in section 5.

Considering our main example which is the hyperbolic heat equation ( 3), one has that RA -1 = 0 0 σ 0 . This matrix is nilpotent so

P(x) = I + xRA -1 = 1 0 σx 1 . ( 11
)
Therefore the unknown of the new formulation is α = (p, xσp + u) t and the matrix of equation ( 10) is

Q(x) = P(x)AP -1 (x) = P(x)AP(-x) = -xσ 1 1 -x 2 σ 2 xσ .
Proposition 3. Stationary solutions of the fully conservative formulation (10) correspond to stationary solutions (i.e. well balanced solutions) of the primal non conservative equation [START_REF] Amadori | Error estimates for well-balanced and time-split schemes on a damped semilinear wave equation[END_REF].

Proof. This is already a consequence of the previous proposition. A direct verification is as follows. Stationary solutions of the conservative formulation (10) satisfy Q(x)α(x) = W where W ∈ R n is arbitrary. In terms of the primal variable it writes P(x)AU(x) = W, that is e RA -1 x AU(x) = W. The differentiation with respect to x yields the identity e RA -1 x (RU(x) + A∂ x U(x)) = 0 which shows that U ∈ U. The proof is ended.

Finite volumes and Riemann solvers in one dimension

In what follows we start from [START_REF] Chandraseckhar | Radiative transfer[END_REF] and shall detail the structure of some Finite Volume techniques. Additionally to the fact that A is constant, we make the assumption that A is symmetric and non singular

A = A t ∈ R n×n and det(A) 0. ( 12
)
For the simplicity of notations the right hand side matrix will also be considered constant, R = R t . However it must be noticed that the fundamental assumption is more the dissipativity of the symmetric part of R, as in [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF]. It is fundamental to prove the main result in theorem 1.

Let us define a new variable

β = Q(x)α = P(x)AU (13)
so that (10) recasts as ∂ t α + ∂ x β = 0. The standard explicit Finite Volume discretization on a grid with varying mesh size

∆x j = x j+ 1 2 -x j-1 2 is α n+1 j -α n j ∆t + β n j+ 1 2 -β n j-1 2 ∆x j = 0 ( 14
)
where β n j+ 1 2 is the flux at time step t n = n∆t. The definition of the scheme relies on the definition of the flux. This will be performed in the context of simplified Riemann solvers as sketched in figure 1. The construction of the Riemann solver is based on the spectral decomposition of the matrix Q(x). Let us denote x * = x j+ 1 2 the point between cell L = j and cell R = j + 1. Since A is symmetric, it admits a basis of real orthonormal eigenvectors, that is Au p = λ p u p where λ p 0 since A is not singular. The point is that Q is similar to A, but is not symmetric. The (right) eigenvectors of Q are

Intermediate state R L

Q(x * )r * p = λ p r * p , r * p = P(x * )u p . (15)
The (left) eigenvectors are

Q t (x * )s * p = λ p s * p , s * p = P(x * ) -t u p . (16)
Based on this decomposition, we construct hereafter two different solvers which are referred to as the one-state solver and the two-states solvers. This distinction originates in a previous work [START_REF] Després | Lois de conservation eulériennes, lagrangiennes et méthodes numériques[END_REF].

A one-state solver

Since Q does not depend on the time variable, we also note that β is solution of the autonomous equation

∂ t β + Q(x)∂ x β = 0. ( 17
)
Let us assume that the local variations of matrix x → Q(x) are smooth enough, so that it is reasonable to freeze the matrix at the interface x j+ 1 2 still maintaining a correct accuracy. This idea is to be compared with similar ones in [START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws, Exponential-Fit, Well-Balanced and Asymptotic-Preserving[END_REF][pages 64-66] or [START_REF] Huang | A conservative, piecewise-steady difference scheme for transonic nozzle flow[END_REF][START_REF] Vanleer | On the relation between the upwind differencing schemes of Engquist-Osher, Godunov and Roe[END_REF]. In our case we modify locally equation [START_REF] Gosse | An asymptotic-preserving well-balanced scheme for the hyperbolic heat equations[END_REF] by

∂ t β + Q(x * )∂ x β = 0, x * = x j+ 1 2 .
Here locally means at the interface between cell j and cell j + 1. Solving this equation in the interval (x * -, x * + ) with Riemann data β(t = 0) = β L for x < x * and β(t = 0) = β R for x * < x, one can compute the Riemann invariants s * p , β using the equation ∂ t s * p , β + λ p ∂ x s * p , β = 0. The solution β * = β(t, x * ), t > 0 is provided by the solution of the linear system

       s * p , β * -β L = 0, λ p > 0, s * p , β * -β R = 0, λ p < 0. ( 18
)
This linear system is non singular since the eigenvectors (s * p ) are linearly independent. It defines a function ϕ :

R n × R n × R → R n such that ϕ(β L , β R , x * ) = β * . ( 19
)
The usual consistency property of Riemann solver writes

ϕ(β, β, x) = β ∀(x, β). (20)
One may call this a one-state Riemann solver because the eigenvectors are common to the same matrix Q(x * ). This will not be the case in the next section. The scheme writes

α n+1 j -α n j ∆t + ϕ(β n j , β n j+1 , x j+ 1 2 ) -ϕ(β n j-1 , β n j , x j-1 2 ) ∆x j = 0. ( 21
)
In terms of the original variable U j = P(x j ) -1 α j , the scheme writes

U n+1 j -U n j ∆t + P(x j ) -1 ϕ(β n j , β n j+1 , x j+ 1 2 ) -ϕ(β n j-1 , β n j , x j-1 2 ) ∆x j = 0. ( 22
)
Proposition 4. The scheme ( 22) is well-balanced.

Proof. Being well-balanced means being exact for stationary solutions. So let us assume that U 0 j = P(x j ) -1 α 0 j with α 0 j = Q(x j ) -1 γ where γ ∈ R n is a given vector. By construction

β 0 j = Q(x j )α 0 j = γ is constant. Therefore ϕ(β 0 j , β 0 j+1 , x j+ 1 2 ) = γ, ∀ j.
In view of ( 22) the scheme is stationary U 1 j = U 0 j for all j. The proof is ended. The well-balanced scheme ( 22) is written as a multiplicative modification of a standard finite volume solver, where multiplicative means that the numerical finite volume discretization of ∂ x (AU) is premultiplied on the left by the matrix P(x j ) -1 . With this respect, the initial additive modification, that is +RU, has been changed into a multiplicative one. A complementary understanding of the scheme ( 22) is provided after rewriting it in a more standard additive formulation. To do so we rewrite [START_REF] Jin | Numerical schemes for hyperbolic conservation laws with stiff relaxation terms[END_REF] as

U n+1 j -U n j ∆t + P(x j+ 1 2 ) -1 β * j+ 1 2 -P(x j-1 2 ) -1 β * j-1 2 ∆x j + P(x j ) -1 -P(x j+ 1 2 ) -1 ∆x j β * j+ 1 2 + P(x j-1 2 ) -1 -P(x j ) -1 ∆x j β * j-1 2 = 0 (23)
where

β * j+ 1 2 = ϕ(β n j , β n j+1 , x j+ 1 2 )
. We define

U * j+ 1 2 = A -1 P(x j+ 1 2 ) -1 β * j+ 1 2 , (24)
and ∆x ± j = x j± 1 2x j so that ( 23) can be rewritten in a slightly more conventional form

U n+1 j -U n j ∆t + A U * j+ 1 2 -U * j-1 2 ∆x j + P(x j ) -1 P(x j+ 1 2 ) -I ∆x j AU * j+ 1 2 + I -P(x j ) -1 P(x j-1 2 ) ∆x j AU * j-1 2 = 0. ( 25
)
Here ,

A U * j+ 1 2 -U * j-1 2 ∆x j
is the discretization of the divergent term A∂ x U, and

P(x j ) -1 P(x j+ 1 2 )-I ∆x j AU * j+ 1 2 + I-P(x j ) -1 P(x j- 1 2
)

∆x j AU * j-1 2
is the discretization of the source RU.

Proposition 5. The quantity U * j+ 1 2 is solution of the linear system            u p , U * j+ 1 2 -e -A -1 R∆x + j U j = 0, λ p > 0, u p , U * j+ 1 2 -e -A -1 R∆x - j+1 U j+1 = 0, λ p < 0. (26)
Proof. We make the elimination β *

j+ 1 2 = P(x j+ 1 2 )AU * j+ 1 2
and β j = P(x j )AU j in [START_REF] Greenberg | A well-balanced scheme for the numerical processing of source terms in hyperbolic equations[END_REF]. Therefore the first line of [START_REF] Greenberg | A well-balanced scheme for the numerical processing of source terms in hyperbolic equations[END_REF] for λ p > 0 rewrites as P(x j+ 1 2 ) -t u p , P(x j+ 1 2 )AU *

j+ 1 2 -P(x j )AU j = 0, or u p , AU * j+ 1 2 -P(x j+ 1 2 ) -1 P(x j )AU j = 0, or Au p , U * j+ 1 2 -A -1 P(x j+ 1
2 ) -1 P(x j )AU j = 0. Since Au p = λ p u p and

A -1 P(x j+ 1 2 ) -1 P(x j )A = A -1 e -RA -1 x j+ 1 2 -x j A = e -A -1 R∆x + j ( 27
)
it yields the first line of the claim. The second part of claim is proved in a similar fashion.

Remark 4 (Elimination of the matrix exponentials). We notice also that the formula ( 27) can recast as

P(x j+ 1 2 ) -1 P(x j ) = P(x j -x j+ 1 2 ) = e -RA -1 ∆x + j ( 28
)
so that the linear system ( 26) is also rewritten as

           u p , U * j+ 1 2 -P(x j ) t P(x j+ 1 2 ) -t U j = 0, λ p > 0, u p , U * j+ 1 2 -P(x j+1 ) t P(x j+ 1 2 ) -t U j+1 = 0, λ p < 0. ( 29
)
It means that the matrix exponentials can be eliminated in the fluxes. We notice also that the matrix exponentials do not show up in the scheme itself [START_REF] Tavelli | A high order semi-implicit discontinuous Galerkin method for the two dimensional shallow water equations on staggered unstructured meshes[END_REF]. This remark will be fundamental in the multiD example at the end of this work: we will see that the matrix exponentials are not well defined, while the matrix P has a natural multiD definition.

Continuing to use the matrix exponentials in one dimension since there are well adapted to our purposes, the formulas ( 23) can now be interpreted.

• If R = 0 which means no relaxation, the matrix exponential terms like e -A -1 R∆x + j are equal to the identity matrix. In this case U *

j+ 1 2 = U Riemann j+ 1 2
is the solution of the standard Riemann solver

           u p , U Riemann j+ 1 2 -U j = 0, λ p > 0, u p , U Riemann j+ 1 2 -U j+1 = 0, λ p < 0. ( 30
)
This construction corresponds to figure 1.

• In all cases the matrix exponentials are related to the solution of the stationary direct problem A∂ x W = -RW. The solution of this auxiliary equation can be written as W(x j+ 1 2 ) = Π(-∆x + j )W(x j ) where Π(-∆x + j ) = e -A -1 R∆x + j stands for the propagator over the distance ∆x + j . So one can rewrite [START_REF] Vanleer | On the relation between the upwind differencing schemes of Engquist-Osher, Godunov and Roe[END_REF] as the solution of the standard Riemann solver

           u p , U Riemann j+ 1 2 -U modified j = 0, λ p > 0, u p , U Riemann j+ 1 2 -U modified j+1 = 0, λ p < 0, (31)
where the modified states are designed by propagating the cells centered states. This is illustrated in figure 2. • The sink term in ( 23) is S = S 1 + S 2 . The first term is

R

S 1 = P(x j ) -1 P(x j+ 1 2 ) -I ∆x j AU * j+ 1 2 = e RA -1 ∆x + j -I ∆x j AU * j+ 1 2 .
A Taylor expansion with respect to mesh size yields

S 1 = ∆x + j ∆x j RU * j+ 1 2 + O(h) where h = sup j ∆x j . A similar analysis yields S 2 = -∆x - j ∆x j RU * j-1 2 + O(h). Considering that U * j+ 1 2 and U * j-1 2
are homogeneous to U and that ∆x + j -∆x - j = ∆x j , one sees that S = S 1 + S 2 ≈ RU is consistent with (the opposite of) the relaxation term.

Let us now particularize these quantities for the example of the hyperbolic heat equation [START_REF] Birkhoff | Harmonic solutions of transport equations[END_REF]. The eigenvectors of A are u 1 = (1, 1) associated to λ 1 > 0 and u 2 = (1, -1) associated to λ 2 = -1. The standard Riemann solver (30) can be written under the form

           p Riemann j+ 1 2 -p j + u Riemann j+ 1 2 -u j = 0, p Riemann j+ 1 2 -p j+1 -u Riemann j+ 1 2 -u j+1 = 0, that is p Riemann j+ 1 2 = p j + p j+1 2 + u j -u j+1 2 and u Riemann j+ 1 2 = u j + u j+1 2 + p j -p j+1 2 .
The modified solver (31) is function of

U modified j = e -A -1 R∆x + j U j where e A -1 Ry = I + yA -1 R = 1 σy 0 1 . ( 32
)
So U modified j = p j -σ∆x + j u j , u j and U modified j+1 = p j+1 -σ∆x - j+1 u j+1 , u j+1 . The new solver (31) writes

                   p * j+ 1 2 = p j + p j+1 2 + (1 -σ∆x + j )u j -(1 + σ∆x - j+1 )u j+1 2 , u * j+ 1 2 = p j -p j+1 2 + (1 -σ∆x + j )u j + (1 + σ∆x - j+1 )u j+1 2 . ( 33
)
This solver is identical to the "Piecewise Steady Approximation", an idea that goes back to [START_REF] Huang | A conservative, piecewise-steady difference scheme for transonic nozzle flow[END_REF], see also [START_REF] Jin | Numerical schemes for hyperbolic conservation laws with stiff relaxation terms[END_REF].

A two-states solver more adapted to discontinuous coefficients

The previous hypothesis that x → Q(x) is continuous with small variation is not always reasonable, in particular for more general problems with discontinuous coefficients as illustrated in figure 3. To illustrate this situation we consider the simple model equation

∂ t U + A∂ x U = -R(x)U (34)
with R(x) = R 1 for x < 0 and R = R 2 for 0 < x. In this simple example one has that A = A 1 = A 2 . Still assuming that A is non singular, one gets that P(x) = e R 1 A -1 x for x < 0 and P(x) = e R 2 A -1 x for 0 < x. If R 1 is very different from R 2 , the matrix P does not have a continuous derivative at x = 0. In this case it is possible to imagine that Q(x) might have strong local variation, so one may question the accuracy of the one-state solver. We now develop a method which can be used in the situation described previously but also extend to more general cases such as

A 1 A 2 . x>0 Zone 1 Zone 2 A =A R = R 2 2 R = R 1 1 A =A x x<0 x=0
Figure 3: Discontinuity of the coefficients. The analysis is performed for A = A 1 = A 2 . On physical grounds, such problems are very similar to scattering problems that can be analyzed with transmissions and reflections operators. This correspondance will be evidenced in section 4.

Our approach is here to upwind the eigenvectors in [START_REF] Greenberg | A well-balanced scheme for the numerical processing of source terms in hyperbolic equations[END_REF], and to consider instead

       s L p , β * * -β L = 0, λ p > 0, s R p , β * * -β R = 0, λ p < 0, (35)
where s L p is a generic eigenvector of the matrix Q(x L ) and s R p is a generic eigenvector of the matrix Q(x R ). We call this a two-states Riemann solver because the eigenvectors are different. The linear system is invertible if and only if the vectors s L P (for λ p > 0) and s R P (for λ p < 0) are linearly independent.

Theorem 1. Assume R + R t ≥ 0. Then the family s L p λ p >0 ∪ s R p λ p <0
is linearly independent, and so the two-states solver (35) is well defined.

Proof. It is clear that if R vanishes, then these eigenvectors are equal to the ones of the matrix A, and so are linearly independent. The key of the proof is the use of the dissipativity hypothesis R + R t ≥ 0 in an appropriate manner.

Let s i (x) = P(x) -t u i denotes a generic eigenvector of Q(x) t , with the convention that s L i is an eigenvector of the matrix

Q(x L ) for 1 ≤ i ≤ k (that is λ L i > 0), and s R i is an eigenvector of the matrix Q(x R ) for k + 1 ≤ i ≤ n (that is now λ R i < 0)
. One has to take care that the convention is here the opposite of the usual one since eigenvalues with low indices are positive and eigenvalues with higher indices are negative. We will show that the only real solution of the

equation k i=1 α i s L i + n i=k+1 α i s R i = 0 is α i = 0 for 1 ≤ i ≤ n, which proves the claim. Let us set z = k i=1 α i s L i = - n i=k+1 α i s R i . (36)
We denote x * a point between the left cell and the right cell:

x L < x * < x R .
• Let us study the function

f (x) =         Ae A -1 Rx *         k i=1 α i s i (x)         , e A -1 Rx *         k i=1 α i s i (x)                 =         Ae A -1 R(x * -x)         k i=1 α i u i         , e A -1 R(x * -x)         k i=1 α i u i                 .
One has that (assuming u i = 1)

f (x * ) =         A         k i=1 α i u i         ,         k i=1 α i u i                 = k i=1 λ i |α i | 2 ≥ 0. ( 37
)
On the other hand a direct calculation shows that

f (x) = -2         Re A -1 R(x * -x)         k i=1 α i u i         , e A -1 R(x * -x)         k i=1 α i u i                 = -         R + R t e A -1 R(x * -x)         k i=1 α i u i         , e A -1 R(x * -x)         k i=1 α i u i                 Here R + R t ≥ 0. So f (x) ≤ 0 and therefore f (x L ) ≥ f (x * ) ≥ 0. ( 38
)
• Similarily we define

g(x) =        Ae A -1 Rx *        n i=k+1 α i s i (x)        , e A -1 Rx *        n i=k+1 α i s i (x)               =        Ae A -1 R(x * -x)        n i=k+1 α i u i        , e A -1 R(x * -x)        n i=k+1 α i u i               .
One has that

g(x * ) =        A        n i=k+1 α i u i        ,        n i=k+1 α i u i               = n i=k+1 λ i |α i | 2 ≤ 0. ( 39
)
We also have that g (x) ≤ 0. Therefore 0

≥ g(x * ) ≥ g(x R ). ( 40
)
• But due to (36), one has that f (x L ) = g(x R ). By comparison with (38) and (40) it yields f (x L ) = g(x R ) = 0. Therefore (38) implies that f (x * ) = 0 which shows (for example with the help of (37)) that α i = 0 for i ≤ k. Similarly one has that g(x * ) = 0 which shows that α i = 0 for k + 1 ≤ i. The proof is ended.

The solution of the linear system (35) can be written with a function ψ :

R n × R n × R × R → R n such that ψ(β L , β R , x L , x R ) = β * * . ( 41
)
We call it the two-states solver. The usual consistency property of Riemann solvers writes

ψ(β, β, x, y) = β ∀(x, y, β). (42)
It is a direct consequence that if one plugs

β L = β R = β in (35), then β * * is a solution, independently of x L < x R .
Since this solution is unique by the theorem, it shows (42). With this notation for the Riemann solver, the scheme (35) writes

α n+1 j -α n j ∆t + ψ(β n j , β n j+1 , x j-1 , x j ) -ψ(β n j-1 , β n j , x j , x j+1 ) ∆x j = 0. ( 43
)
Proposition 6. The scheme (43) is well-balanced (same proof as the one of proposition 4).

As before it is possible to get a more classical understanding of (43) by rewriting it as

U n+1 j -U n j ∆t + P(x j+ 1 2 ) -1 β * * j+ 1 2 -P(x j-1 2 ) -1 β * * j-1 2 ∆x j + P(x j ) -1 -P(x j+ 1 2 ) -1 ∆x j β * * j+ 1 2 + P(x j-1 2 ) -1 -P(x j ) -1 ∆x j β * * j-1 2 = 0 (44)
where

β * * j+ 1 2 = ψ(β n j , β n j+1 , x j-1 , x j ).
The generic flux is identified as

P(x j+ 1 2 ) -1 β * * j+ 1 2 .
To make this statement more explicit we define 45) is solution of the linear system

U * * j+ 1 2 = A -1 P(x j+ 1 2 ) -1 β * * j+ 1 2 . ( 45
) Proposition 7. The quantity U * * j+ 1 2 in (
           u p , e A -1 R∆x + j U * * j+ 1 2 -U j = 0, λ p > 0, u p , e A -1 R∆x - j+1 U * * j+ 1 2 -U j+1 = 0, λ p < 0, ( 46
)
where the matrix exponentials are the opposite of the ones in the one-state solver [START_REF] Vanleer | On the relation between the upwind differencing schemes of Engquist-Osher, Godunov and Roe[END_REF].

Remark 5 (Elimination of the matrix exponentials). As it was done in remark 4 for the one-state solver, it is possible to rewrite (46) using only the matrix P. It yields

           u p , P(x j ) -t P(x j+ 1 2 ) t U * * j+ 1 2 -U j = 0, λ p > 0, u p , P(x j+1 ) -t P(x j+ 1 2 ) t U * * j+ 1 2 -U j+1 = 0, λ p < 0. (47)
Proof. Indeed (35) recasts as

           P(x j ) -t u p , P(x j+ 1 2 )AU * * j+ 1 2 -P(x j )AU j = 0, λ p > 0, P(x j+1 ) -t u p , P(x j+ 1 2 )AU * * j+ 1 2 -P(x j+1 )AU j+1 = 0, λ p < 0, that is            u p , P(x j ) -1 P(x j+ 1 2 )AU * * j+ 1 2 -AU j = 0, λ p > 0, u p , P(x j+1 ) -1 P(x j+ 1 2 )AU * * j+ 1 2 -AU j+1 = 0, λ p < 0, or also            Au p , A -1 P(x j ) -1 P(x j+ 1 2 )AU * * j+ 1 2 -U j = 0, λ p > 0, Au p , A -1 P(x j+1 ) -1 P(x j+ 1 2 )AU * * j+ 1 2 -U j+1 = 0, λ p < 0.
Since Au p = λ p u p it can written as

           u p , A -1 P(x j ) -1 P(x j+ 1 2 )AU * * j+ 1 2 -U j = 0, λ p > 0, u p , A -1 P(x j+1 ) -1 P(x j+ 1 2 )AU * * j+ 1 2 -U j+1 = 0, λ p < 0.
One has the identity A -1 P(x j ) -1 P(x j+ 1 2 )A = A -1 e RA -1 ∆x + j A = e A -1 R∆x + j from which one gets the first line of the claim. The second line is obtained similarly. The proof is ended.

The interpretation of two-states solver, that is (44) with (46), is quite close to the interpretation of the one state solver. We nevertheless observe that the analogue of ( 30) is now

           u p , U Riemann and modified j+ 1 2 -U j = 0, λ p > 0,
u p , U Riemann and modified

j+ 1 2 -U j+1 = 0, λ p < 0, (48)
where the modified states are the propagation of the intermediate state. The modification is the backward propagator Π(∆x + j ) which is the inverse of the direct propagator Π(-∆x + j ) that corresponds to the one-state solver. This is illustrated in figure 4. Remark 6 (Interpretation of figure 4). In the language of the localization method of L. Gosse, the explanation of the figure is that a stationary state is inserted into the Riemann solver, between the left and right states: this construction can somewhat be considered as arbitrary. Other interpretations of the figure are probably possible, but we will not pursue them in this work since they are not needed. With the proposed approach, the internal stationary state does not really exists: it is just the consequence of a certain linear algebra with standard Riemann solvers and standard conservative methods.

Going back once again to the hyperbolic heat equation and using (32), the system (46) rewrites as

           p * * j+ 1 2 + σ∆x + j u * * j+ 1 2 -p j + u * * j+ 1 2 -u j = 0, p * * j+ 1 2 + σ∆x - j+1 u * * j+ 1 2 -p j+1 -u * * j+ 1 2 -u j+1 = 0. The solution writes          p * * j+ 1 2 = 1-σ∆x - j+1 2+σ∆x + j -σ∆x - j+1 (p j -u j ) + 1+σ∆x + j 2+σ∆x + j -σ∆x - j+1 (p j+1 -u j+1 ), u * * j+ 1 2 = 1 2+σ∆x + j -σ∆x - j+1 (u j + u j+1 + p j -p j+1 ).
This solver coincides exactly with the Gosse-Toscani solver in [START_REF] Gosse | An asymptotic-preserving well-balanced scheme for the hyperbolic heat equations[END_REF].

Application to the S n model

The S n model is a natural way to approach the transfer equation

∂ t I + µ∂ x I = σ (< I > -I) , -1 ≤ µ ≤ 1, (49)
where the mean value is < I >= 1 2 1 -1 I(µ)dµ and the velocity µ = cos θ is representative of a direction. This is an integro-differential equation which has been widely studied in relation with the theory of propagation of light, neutrons and other types of particles. We refer to the seminal contributions of [START_REF] Chandraseckhar | Radiative transfer[END_REF][START_REF] Case | Linear transport theory[END_REF][START_REF] Dautray | Mathematical Analysis and Numerical Methods for Science and Technology[END_REF] and to the more dedicated work [START_REF] Waterman | Matrix exponential description of radiative transfer[END_REF]. Most of the algebra we show come from [START_REF] Chandraseckhar | Radiative transfer[END_REF] and is given to get complete interpretation of some recent schemes proposed in [START_REF] Gosse | Transient radiative transfer in the grey case: Well-balanced and asymptotic-preserving schemes built on Case's elementary solutions[END_REF].

One usually choose a set of velocities 0 < µ 1 < µ 2 < . . . < µ n ≤ 1 and positive weights w i > 0 such that 2 i=1 w i = 1 and 2 i=1 w i µ 2 i = 1 3 . The velocities are actually the cosine of some angles and this is why there are less than one. The intensity is represented by

I(x, t, µ) = i=1 w i f i (x, t)δ(µ -µ i ) + i=1 w i g i (x, t)δ(µ + µ i ), (50)
where δ indicates a Dirac mass. The f i and g i are solutions to

       ∂ t f i + µ i ∂ x f i = σ j=1 w j ( f j + g j ) -f i , ∂ t g i -µ i ∂ x g i = σ j=1 w j ( f j + g j ) -g i . (51)
Let us define the vector

U = f , g t ∈ R 2n where f = ( f 1 , . . . , f n ) t ∈ R n and g = ( g 1 , . . . , g n ) t ∈ R n . Equations (51) are recast as a system ∂ t U + A∂ x U = -R U (52)
where

A = D 0 0 -D ∈ R 2n×2n (with D = diag(µ i ) ∈ R n×n
= ( f, g) t = ( √ E f , √ E g) t ∈ R 2n which is solution of ∂ t U + A∂ x U = -RU, R = σ (I d -w ⊗ w) , (53)
where the relaxation matrix is symmetric. Since the coefficient is non negative, that is σ ≥ 0, and the new vector w = ( √ w 1 , . . . , √ w n , √ w 1 , . . . , √ w n ) ∈ R 2n is such that w = 1, the relaxation matrix is now non negative, that is R ≥ 0. Under this symmetrized form, one may apply directly the previous theory since A is non singular. It remains to characterize the one-state solver and the two-states solver. This will be performed using some basic symmetry principles in combination with the transmission and reflexion operators which are fundamental in the theory of light propagation. Some parts of the following analysis can be adapted to more general systems as in [START_REF] Arnold | Large-time behavior of discrete equations with nonsymmetric interactions[END_REF][page 1558].

The transmission and reflexion operators

Denoting Q = A -1 R, the matrix exponential e -Qx = e -A -1 Rx describes the result of the propagation of rays of light (50) through a slab [0, x] of material with absorption σ. The orientation of the slab is positive, that is x ≥ 0. One has the relation

f (x) g(x) = e -A -1 Rx f (0) g(0) . ( 54
)
The theory of light shows it is always possible to define a transmission operator (matrix) T (x) ∈ R n×n and a reflection operator (matrix) S(x) ∈ R n×n so that one can compute the pair ( f (x), g(0)) in function of ( f (0), g(x)). These operators can be defined as a consequence of an energy identity which is obtained by integration by part for stationary solutions of (53)

(D f (0), f (0)) + (Dg(x), g(x)) = (D f (x), f (x)) + (Dg(0), g(0)) + x 0 (RU(s), U(s))ds. ( 55
)
Since R ≥ 0 and D > 0, it shows that all quantities vanish if f (0) = g(x) = 0. So the linear transformation

f (0) g(x) = L f (0) g(0)
is non singular (that is invertible), and all unknowns can be expressed in function of the vector

f (0) g(x)
. It means one can define operators (matrices)

T 1 (x), T 2 (x), S 1 (x), S 2 (x) ∈ R n×n such f (x) g(0) = T 1 (x) S 1 (x) S 2 (x) T 2 (x) f (0) g(x)
.

Based on symmetry considerations for the physical problem, one has that T 2 (x) = T 1 (x) and S 2 (x) = S 1 (x). Therefore one can write f (x)

g(0) = T (x) S(x) S(x) T (x) f (0) g(x) (56)
after dropping the indices. The transmission or transfer operator is T (x), the reflexion or scattering operator is S(x). The matrix in (56) is the scattering matrix. The input is the pair ( f (0), g(x)) and the output is ( f (x), g(0)).

Of course T (0) = I d and S(0) = 0. Some basic properties can be proved, which are used to justify the fundamental composition relations (64).

Proposition 8. The transmission operator is non singular and the reflexion operator is strictly bounded, that is det (T (x)) 0 and D

1 2 S(x)D -1 2 < 1, ∀x ≥ 0. ( 57
)
Proof. The proof is by contradiction.

• Assume T to be singular. There would exist a non zero vector W ∈ R n such that

T (x) S(x) S(x) T (x) W 0 = 0 . . . , W 0.
It yields a non zero solution of (52) with vanishing Cauchy data at the end point x, that is f (0) = W 0 and f (x) = g(x) = 0. But this is impossible (54). This is why the transmission operator is non singular.

• Assume now that D 1 2 S(x)D -1 2 ≥ 1. There would exist a non zero vector Z ∈ R n such that D

1 2 S(x)D -1 2 Z ≥ Z 0. We set W = D -1 2 Z so that D 1 2 S(x)W ≥ D 1 2 W with W 0. ( 58
)
On the other hand the identity

T (x) S(x) S(x) T (x) 0 W = S(x)W T (x)W
can be interpreted with the energy identity (55), where the input is ( f (0), g(x)) = (0, W) and the ouput is ( f (x), g(0)) = (S(x)W, T (x)W). The energy identity yields the inequality D

1 2 W 2 ≥ D 1 2 S(x)W 2 + D 1 2 T (x)W 2 . Using (58), one obtains D 1 2 T (x)W = 0 and W = 0 since D 1 2 T (x)
is invertible. This is a contradiction. The proof is ended.

Proposition 9. All coefficients of T (x) and S(x) are non negative, and one has the identity

(T (x) + S(x)) t z = z ( 59
)
where z = µ 1 √ w 1 , . . . , µ n √ w n ∈ R n .
Proof. A proof of (59) can be performed using the identity Rw = 0. But it is better for physical intuition to start from the stationary equations

       µ i ∂ x f i = σ j=1 w j ( f j + g j ) -f i , -µ i ∂ x g i = σ j=1 w j ( f j + g j ) -g i . ( 60
)
It yields after integration i=1 µ i f i (0)

+ g i (x) = i=1 µ i f i (x) + g i (0) . The change of unknowns f i = w -1 2 i f i and g i = w -1 2 i g i yields i=1 z i ( f i (0) + g i (x)) = i=1 z i ( f i (x) + g i (0)), that is in vectorial form f (x) g(0) , z z = f (0) g(x) , z z .
The representation (56) yields

T (x) S(x) S(x) T (x) f (0) g(x) , z z = f (0) g(x) , z z
which can be simplified into

f (0) + g(x), (T (x) + S(x)) t z = ( f (0) + g(x), z) .
Since this is true for all f (0) + g(x) ∈ R n , it shows the identity (59).

Next we define h i (xy) = g i (y) for y ∈ [0, x]. The system (60) recasts as (0

≤ y ≤ x)        µ i ∂ x f i (y) + σ f i (y) = σ j=1 w j f j (y) + h j (x -y) , µ i ∂ x h i (y) + σ h i (y) = σ j=1 w j f j (x -y) + h j (y) .
For such a system a standard property is that if f i (0) ≥ 0 and h i (0) ≥ 0 for all i, then f i (x) ≥ 0 and h i (x) ≥ 0 for all i.

It shows the first part of the claim after a change of unknowns. The proof is ended.

Proposition 10. Assume that x ≥ 0. One has the representation formula

e -A -1 Rx = T -ST -1 S ST -1 -T -1 S T -1 (x). ( 61
)
The inverse formula is

e A -1 Rx = T -1 -T -1 S ST -1 T -ST -1 S (x). ( 62
)
This formula shows that the matrix on the left hand side which belongs to R 2n×2n can be represented with two matrices in R n×n , thus resulting in a reduction of the size of the basic objects used in the representation formulas.

Proof. Indeed (56) shows that g(x) = T -1 (x)(-S(x) f (0) + g(0)) and

f (x) = T (x) -S(x)T -1 (x)S(x) f (0) + T -1 (x)S(x)g(0).
It ends the proof of (61). That (62) is inverse formula of (61) is evident.

Proposition 11. The transposed operators are T t = DT D -1 and S t = DSD -1 .

Proof. Since e -A -1 Rx t = e -RA -1 x = Ae -A -1 Rx A -1 and A is block diagonal with D and -D on the diagonal, one has from (61) T -t = (-D)T -1 (-D) -1 which yields the first relation. The second relation can be deduced from -T -1 S t = D ST -1 (-D) -1 which yields the second relation after elimination of the transmission operator.

Let us consider x ≤ y ≤ z together with T 1 = T (yx), S 1 = S(yx), T 2 = T (zy), S 2 = S(zy), T 3 = T (zx) and S 3 = S(zx). Composition formulas can be deduced for the transmission and reflexion operators. It comes from the composition-commutation relations

e -A -1 Rd 3 = e -A -1 Rd 2 e -A -1 Rd 1 = e -A -1 Rd 1 e -A -1 Rd 2 , d 3 = d 1 + d 2 .
(63)

Proposition 12. One has the composition formulas

           T 3 = T 1 (I d -S 2 S 1 ) -1 T 2 = T 2 (I d -S 1 S 2 ) -1 T 1 , S 3 = S 2 + T 2 S 1 (I d -S 2 S 1 ) -1 T 2 = S 2 + T 2 (I d -S 1 S 2 ) -1 S 1 T 2 = S 1 + T 1 S 2 (I d -S 1 S 2 ) -1 T 1 = S 1 + T 1 (I d -S 2 S 1 ) -1 S 2 T 1 . (64)
where the matrices I d -S 2 S 1 and I d -S 1 S 2 can be shown non singular as consequence of the bound in (57).

Proof. The composition relation (63) yields

T 3 -S 3 T -1 3 S 3 S 3 T -1 3 -T -1 3 S 3 T -1 3 = T 2 -S 2 T -1 2 S 2 S 2 T -1 2 -T -1 2 S 2 T -1 2 T 1 -S 1 T -1 1 S 1 S 1 T -1 1 -T -1 1 S 1 T -1 1 . (65)
• It yields

T -1 3 = -T -1 2 S 2 S 1 T -1 1 + T -1 2 T -1 1 = T -1 2 (I d -S 2 S 1 ) T -1 1 .
It shows the first identity of the claim is deduced after inversion. Since the matrix exponentials commute, the other identity T 3 = T 2 (I d -S 1 S 2 ) -1 T 1 is immediate after switching the indices.

• One also has from (65)

S 3 T -1 3 = T 2 -S 2 T -1 2 S 2 S 1 T -1 1 + S 2 T -1 2 T -1 1 = T 2 S 1 T -1 1 + S 2 T -1 2 (I d -S 2 S 1 ) T -1 1 .
Therefore

S 3 = T 2 S 1 T -1 1 + S 2 T -1 2 (I d -S 2 S 1 ) T -1 1 T 1 (I d -S 2 S 1 ) -1 T 2 = T 2 S 1 (I d -S 2 S 1 ) -1 T 2 + S 2
which is the third identity of the claim. The fourth one comes from the identity S 1 (

I d -S 2 S 1 ) -1 = (I d -S 1 S 2 ) -1 S 1 .
The two last ones are obtained by switching the indices.

The two-states solver for S n model

The two-states solver ( 22)-( 25),

U n+1 j -U n j ∆t + 1 ∆x j e RA -1 ∆x + j AU * * j+ 1 2 -e RA -1 ∆x - j AU * * j-1 2
= 0, can be rewritten under the general form

U n+1 j -U n j ∆t + 1 ∆x j A e A -1 R∆x + j U * * j+ 1 2 -e A -1 R∆x - j U * * j-1 2 = 0. ( 66
)
With some additional algebra, one can prove the following elegant and compact result, which shows this method is also equal to the one published in [START_REF] Gosse | Transient radiative transfer in the grey case: Well-balanced and asymptotic-preserving schemes built on Case's elementary solutions[END_REF] and [START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws, Exponential-Fit, Well-Balanced and Asymptotic-Preserving[END_REF][formula (9.18) page 173], but is different from the one from [START_REF] Jin | Numerical schemes for hyperbolic conservation laws with stiff relaxation terms[END_REF].

Proposition 13. The two-states solver for the S n model can be written under the form

                     f n+1 j -f n j ∆t + D f n j -T j-1 2 3 f n j-1 -S j-1 2 3 g n j ∆x j = 0, g n+1 j -g n j ∆t -D S j+ 1 2 3 f n j + T j+ 1 2 3 g n j+1 -g n j ∆x j = 0. ( 67
)
where

T j-1 2 3 = T 3 (x j -x j-1 ) and S j-1 2 3 = S 3 (x j -x j-1
) for all j.

Proof. In the proof we do not note the reference to the interface for the transmission and reflection operators. Let us write the flux as U * *

j+ 1 2 = f * * j+ 1 2 , g * * j+ 1 2 t
, defined for all j by (46). We observe also that the eigenvectors u p in (46) are extremely simple to determine. Indeed they are of type u p = ( f p , 0) for λ p > 0, and of type u p = (0, g p ) for λ p < 0. With natural notations the system (46) writes

T -1 1 -T -1 1 S 1 -T -1 2 S 2 T -1 2         f * * j+ 1 2 g * * j+ 1 2         = f j g j+1
.

By comparison with (62) one has e

A -1 R∆x + j = T -1 1 -T -1 1 S 1 S 1 T -1 1 T 1 -S 1 T -1 1 S 1
with the notation T 1 = T (∆x + j ) and S 1 = S(∆x + j ). Therefore we define (with natural notations)

        f - j+ 1 2 g - j+ 1 2         = e A -1 R∆x + j U * * j+ 1 2 = T -1 1 -T -1 1 S 1 S 1 T -1 1 T 1 -S 1 T -1 1 S 1 =A T -1 1 -T -1 1 S 1 -T -1 2 S 2 T -1 2 -1 =B f j g j+1 .
One has the algebra

AB = T 1 -S 1 T -1 1 S 1 S 1 T -1 1 -T -1 1 S 1 T -1 1 -1 B = T -1 1 -T -1 1 S 1 -T -1 2 S 2 T -1 2 T 1 -S 1 T -1 1 S 1 S 1 T -1 1 -T -1 1 S 1 T -1 1 -1 = I 0 -S T -1 2 (I -S 2 S 1 ) T -1 1 -1 = I 0 -S T -1 3 -1 with S = T -1 2 S 2 T 1 + T -1 2 (I -S 2 S 1 ))T -1 1 S 1 . So AB = I 0 T 3 S T 3
. One can check that

T 3 S = T 1 (I d -S 2 S 1 ) -1 T 2 T -1 2 S 2 T 1 + T -1 2 (I -S 2 S 1 ))T -1 1 S 1 = T 1 (I d -S 2 S 1 ) -1 S 2 T 1 + S 1 = S 3
using the last identity of (64). Therefore one can write

        f - j+ 1 2 g - j+ 1 2         = f j S 3 f j + T 3 g j+1
. By symmetry one gets

        f + j-1 2 g + j-1 2         = S 3 f j-1 + S 3 g j g j
. The proof is ended.

A similar algebra can be performed for the one-state solver, but the results are less interesting for the moment. One starts from

U n+1 j -U n j ∆t + 1 ∆x j A e A -1 R∆x + j U * j+ 1 2 -e A -1 R∆x - j U * j-1 2 = 0. ( 68
)
It remains to identify the flux U *

j+ 1 2 = f * j+ 1 2 , g * j+ 1 2 t
which is defined for all j by (29). The operator that shows up in the first line of ( 29) is simply e -A -1 R∆x + j for which we can use the direct formula (61). Same kind of remarks hold for the second line of (29) but with the second line of (62) used with -x = -∆x - j+1 > 0. So one finds out that

         f * j+ 1 2 = T 1 -S 1 T -1 1 S 1 f j + S 1 T -1 1 g j , g * j+ 1 2 = S 2 T -1 2 f j+1 + T 2 -S 2 T -1 2 S 2 g j+1 ,
where T 2 = T (-∆x - j+1 ) and S 1 = S(-∆x - j+1 ). Therefore we define with natural notations

        f - j+ 1 2 g - j+ 1 2         = e A -1 R∆x + j U * j-1 2 = T -1 1 -T -1 1 S 1 S 1 T -1 1 T 1 -S 1 T -1 1 S 1        T 1 -S 1 T -1 1 S 1 f j + S 1 T -1 1 g j S 2 T -1 2 f j+1 + T 2 -S 2 T -1 2 S 2 g j+1        .
It yields the expression of the fluxes in terms of the reflexion and transmission operators. However this formulation is still quite involved and will not be pursued in this work.

Theorem 2. The scheme (67) is positive, L ∞ and L 1 stable under CFL Proof. The positivity of the first line of (67) is a consequence of the explicit Euler reformulation

f n+1 j = I - ∆t ∆x j D f n j + ∆t ∆x j DT j-1 2 3 f n j-1 + ∆t ∆x j DS j-1 2 3 g n j . (69)
Assuming the CFL condition ∆t ∆x j × spectral radius(D) ≤ 1, all matrices are non negative. Assume the vectors f n j , f n j-1 and g n j have non negative coefficients. Therefore the vector f n+1 j has also non negative coefficients. A similar algebra yields for the second line of (67)

g n+1 j = I - ∆t ∆x j D g n j + ∆t ∆x j DS j+ 1 2 3 f n j + ∆t ∆x j DT j+ 1 2 3 g n j . (70)
Let z ∈ R n be the vector with positive coefficients defined in proposition 9. Using the property of the proposition, one can rewrite (69-70) as

           cz -f n+1 j = I -∆t ∆x j D cz -f n j + ∆t ∆x j DT j-1 2 3 cz -f n j-1 + ∆t ∆x j DS j-1 2 3 cz -g n j , cz -g n+1 j = I -∆t ∆x j D cz -g n j + ∆t ∆x j DS j+ 1 2 3 cz -f n j + ∆t ∆x j DT j+ 1 2 3
czg n j .

(71)

Assuming that data are bounded in L ∞ at step t n , it is possible to find c > 0 large enough so that all quantities are non negative at t n , and remain non negative at time t n+1 . It shows an upper bound which implies the L ∞ stability (under CFL).

To show the stability in L 1 it is sufficient to consider that non negative solutions of (69-70) satisfy one conservation law. In view of proposition 9, one readily proves that (69-70) imply

j ∆x j D -1 z, f n+1 j + g n+1 j = j ∆x j D -1 D -1 z, f n j + g n j . (72)
Therefore non negative initial data which are bounded in L 1 remain bounded in L 1 . The case of a general initial data bounded in L 1 is as usual treated with a splitting of between positive and negative parts.

Remark 7. To prove the stability in L 2 , it is possible to redo the analysis using identities like T t DT + S t DS ≤ D which are consequences of the energy identity (55). This strategy is more involved and so is not developed in this work. Another possibility is to use the Riesz-Thorin theorem in interpolation theory [START_REF] Zygmund | Trigonometric series[END_REF] which automatically yields that the stability in L 1 and L ∞ yields the stability in any L p , 1 ≤ p ≤ ∞.

The conservation law (72) in the proof corresponds to the preservation of the total energy d dt R

1

-1 Idxdµ = 0 that comes from (49). Additional ones stemming from the general structure and rewriting in terms of the unknown α are given in proposition 14.

As stressed previously, all these schemes have a fully conservative interpretation using the variable α. A verification for the two-states S n solver is as followed. Proposition 14. Consider the two-states solver (67). Assume for simplicity that the initialization has compact support, that is f 0 j , g 0 j = 0 for | j| ≥ M ∈ N . Then the solution satisfies 2n conservation which can be expressed as Proof. The equality rewrites as j α n j = j α 0 j which the conservation laws for the variable α. It is possible to express the powers e RA -1 j∆x = e RA -1 ∆x j in function of the transmission and reflexion operators. To prove the equality we consider the difference of two iterates. One has

∆x j e RA -1 j∆x f n+1 j -f n j g n+1 j -g n j = - j e RA -1 j∆x A f n j -T 3 f n j-1 -S 3 g n j S 3 f n j + T 3 g n j+1 -g n j = -A j e A -1 R j∆x f n j -S 3 g n j S 3 f n j -g n j + -T 3 f n j-1 0 + 0 T 3 g n j+1 = -A j e A -1 R j∆x                     f n j -S 3 g n j S 3 f n j -g n j + e A -1 R∆x -T 3 f n j 0 + e -A -1 R∆x 0 T 3 g n j =W n j                     .
The representation formulas (61-62) show that

W n j = f n j -S 3 g n j S 3 f n j -g n j + e A -1 R∆x -T 3 f n j 0 + e -A -1 R∆x 0 T 3 g n j = f n j -S 3 g n j S 3 f n j -g n j + T -1 3 -T -1 3 S 3 S 3 T -1 3 T 3 -S 3 T -1 3 S 3 -T 3 f n j 0 + T 3 -S 3 T -1 3 S 3 S 3 T -1 3 -T -1 3 S 3 T -1 3 0 T 3 g n j = f n j -S 3 g n j S 3 f n j -g n j + -f n j -S 3 f n j + S 3 g n j g n j = 0.
Therefore ∆x j e RA -1 j∆x f n+1 j f n j g n+1 j g n j = 0. It shows the result by recurrence. The proof is ended.

One of these 2n conservation laws is actually the physical one j∈Z f n j + g n j , D -1 z = j∈Z f 0 j + g 0 j , D -1 z . It is easily seen using the definition of the scheme (67) and the identity (59).

The singular case

For many applications, the matrix A may be singular. We will detail some consequences on a example, which is the P 1 model coupled a linear temperature equation. It writes

         ∂ t p +∂ x u = τ(T -p), ∂ t u +∂ x p = -σu, ∂ t T = τ(p -T ), (73)
for which

A =           0 1 0 1 0 0 0 0 0           and R =           τ 0 -τ 0 σ 0 -τ 0 τ           .
We notice that A = A t is singular and R = R t ≥ 0. Assuming that σ, τ > 0, the solutions of the adjoint stationary equation satisfy

∂ x u = 0, ∂ x p = σ u, T = p.
Notice the additional inequality σ ≥ τ for physically motivated problems. One can take V 1 = (1, 0, 1) and V 2 = (σx, 1, σx) which yields two and only two linearly independent functions. So dim(V) = 2 < 3 which indicates a degeneracy. Therefore the general method developed previously cannot be used directly and must be adapted. It can be done at least in two directions.

Shifting of the spectrum

The idea is that stationary solutions (of the adjoint equation) are just particular special solutions. We can use more general particular solutions (A and R being symmetric)

∂ t V + A∂ x V = RV, ∂ t V + ξ∂ x V = 0 (74)
where ξ ∈ R is an arbitrary real number. It yields A ξ ∂ x V = RV where A ξ = A -ξI is a shift of the matrix A with det(A ξ ) 0. The spectrum is shifted as well. One has to take care that the test functions V ξ 1 (x, t), . . . , V ξ n (x, t) particular solutions of (74) depend also on the time variable.

However one can expect some kind of degeneracy as ξ → 0. Since it is difficult to determine in advance the behavior of the method in this regime, we do not pursue in this direction.

Ad-hoc integration in time

We present the idea on the example (73). Since dim(V) = 2 < 3, one can write two conservative equations for

α 1 = (U, V 1 ) = p + T and α 2 = (U, V 2 ) = σx(p + T ) + u ∂ t α 1 + ∂ x u = 0, ∂ t α 2 + ∂ x (uxσ + p) = 0.
Defining arbitrarily α 3 = T , one obtains the system

           ∂ t α 1 + ∂ x (-σxα 1 + α 2 ) = 0, ∂ t α 2 + ∂ x (1 -σ 2 x 2 )α 1 + σxα 2 -α 3 = 0, ∂ t α 3 = τ(α 1 -2α 3 ). (75)
The last equation is more an ODE and can be integrated in time since (e 2τt α 3 ) = τe 2τt α 1 . The solution is α 3 (x, t) = e -2τt t 0 τe 2τs α 1 (x, s)ds + e -2τt α 3 (x, 0). One can now perform a partial discretization in time during a time step of length ∆t. With an explicit Euler formulation for example, it yields α 3 (x, ∆t) ≈ e -2τ∆t ∆t 0 τe 2τs ds α 1 (x, 0) + e -2τ∆t α 3 (x, 0) = 1 2 1e -2τ∆t α 1 (x, 0) + e -2τ∆t α 3 (x, 0).

Plugging this approximation in (75) yields the system of two conservation laws

∂ t α 1 + ∂ x (-σxα 1 + α 2 ) = 0, ∂ t α 2 + ∂ x (1 -σ 2 x 2 )α 1 + σxα 2 -1 2 1 -e -2τ∆t α 1 -e -2τ∆t α 3 (0) = 0. ( 76
)
This approximation is valid during the time step ∆t. Within this time step it is possible to use a standard Riemann solver for systems of conservation laws to get a well-balanced scheme with explicit contribution of the time step, as in [START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws, Exponential-Fit, Well-Balanced and Asymptotic-Preserving[END_REF][page 185].

MultiD

This section is devoted to show that the definition of a multidimensional well-balanced scheme displays a much richer structure than in the one dimensional case. Indeed the set of adjoint stationary states may have an infinite dimension, while the dimension is only finite in one dimension. A solution will be detailed on the example of the 2D hyperbolic heat equation with the change-of-basis matrix P(x) which seems to be the correct object to manipulate in higher dimensions.

2D hyperbolic heat equation

We consider the two dimensional hyperbolic heat equation. The primal formulation writes

         ∂ t p +∂ x u +∂ y v = 0, ∂ t u +∂ x p = -σu, ∂ t v +∂ y p = -σv, (77)
where we take that σ > 0 is constant. The adjoint stationary states ( p, u, v) that correspond to (4-5) are solutions of

         ∂ x u +∂ y v = 0, ∂ x p = σ u, ∂ y p = σ v. So ( u, v) = 1 σ ∇ p and ∆ p = 0. Therefore V = ( p, u, v); p is harmonic and ( u, v) = 1 σ ∇ p .
In one dimension, the set of harmonic functions reduces to affine functions and so coincides with [START_REF] Coquel | Asymptotic preserving scheme for Euler system with large friction[END_REF], as observed in [START_REF] Birkhoff | Harmonic solutions of transport equations[END_REF] in the context of stationary transport equations. In two dimension, we define 6.2. Discretization of the 2D hyperbolic heat equation Using the structure (82), we show how to recover a well-balanced (and asymptotic diffusion preserving) scheme with corner fluxes recently designed in [START_REF] Buet | Design of asymptotic preserving schemes for the hyperbolic heat equation on unstructured meshes Numerish Mathematik[END_REF]. We start from a nodal based finite volume scheme for the homogeneous equation. The corresponding semi-discrete scheme writes with original notations on a 2D grid s j p j (t) + r l jr n jr • u r = 0, s j u j (t) + r l jr n jr p jr = 0.

(83)

Denoting n jr = cos θ jr , sin θ jr , it admits the reformulation more adapted to our purposes s j U j (t) + r l jr A jr U jr = 0, (

where U j = (p j , u j ), U jr = (p jr , u r ) and A jr = cos θ jr A + sin θ jr B = A t jr is symmetrix matrix. We now desire to modify this scheme in order to discretized the 2D conservative formulation (78), by means of an extension of the so-called two states solver.

We first write the analogue of the multiplicative scheme [START_REF] Jin | Asymptotic preserving (AP) schemes for multiscale kinetic and hyperbolic equations: a review[END_REF][START_REF] Jin | Numerical schemes for hyperbolic conservation laws with stiff relaxation terms[END_REF] as

The product of matrices can be rewritten using the composition formula (81). It writes under a more local formula P(x j ) -1 P(x r ) = P(x rx j ). These two formulations are the same since we assume the generalization of (24) under the form P(x r )A jr U jr = β jr .

(87)

It remains to define the fluxes to close the system.

To this end we decide first to preserve the structure of the fluxes defined in [START_REF] Buet | Design of asymptotic preserving schemes for the hyperbolic heat equation on unstructured meshes Numerish Mathematik[END_REF]. That is we keep U jr = (p jr , u r ) which means that the first component p jr ∈ R is delocalized around the node x r while the two other components u r ∈ R 2 are "attached" to the node. For convenience we generalize the two states solver starting from (47).

Proposition 15. The corner based 2D scheme based on the two-states solver is equal to the one published in [START_REF] Buet | Design of asymptotic preserving schemes for the hyperbolic heat equation on unstructured meshes Numerish Mathematik[END_REF]. The corner based linear solver writes p j r + σ x rx j , u rp j + n jr , u ru j = 0, j l jr n jr p jr = 0.

(88) Remark 8. What we call the linear vertex-based solver (around x r ) is made of the linear equations ( 89)-(90). The unknowns are u r and the p jr s for all cells j around the node x r . Unfortunately one cannot rely on theorem 1 (in 1D) to show this linear system is well posed. However it is been shown in [START_REF] Buet | Design of asymptotic preserving schemes for the hyperbolic heat equation on unstructured meshes Numerish Mathematik[END_REF] that this system is in practice non singular, for a wide range of meshes. Moreover the final scheme has been proved to be diffusion asymptotic preserving.

Proof. Indeed the matrix P(x j ) -t P(x j+ 1 2 ) t in one dimension becomes P(x j ) -t P(x r ) t = P(x r -

x j ) t =          
1 σ(x rx j ) σ(y ry j ) 0 1 0 0 0 1

         
. Notice that l jr n jr is equal to the orthogonal vector to the half of the vector that starts at x r-1 and finish at x r+1 . The center of the cell is an arbitrary point inside the cell.

The eigenvector u p for λ p > 0 is now in this context the eigenvector u jr of the matrix A jr associated a positive eigenvalue. There is no ambiguity since there is only one positive eigenvalue (equal to 1). One obtains the eigenvector u jr = (1, n jr ) t since A more general family of Friedrichs type is A 0 (x)∂ t U + ∂ x (A(x)U) + ∂ y (B(x)U) = S (x) -R(x)U for which the same ideas of duality and change of basis should apply identically. It can be used to treat the example evoked in remark 1. But the method can be developed in many other directions which are briefly detailed below. First directions already evoked in this work are: asymptotic preserving methods which usually start from wellbalanced techniques plus small parameters; and 2D formulations which display an interesting structure where the matrix P plays a major role. An open problem is to determine wether the commutation-composition relations (81) are necessary or if there are true only for our example. The effective calculation of the matrix exponentials for the S n model can be performed starting from the references [START_REF] Chandraseckhar | Radiative transfer[END_REF][START_REF] Case | Linear transport theory[END_REF][START_REF] Gosse | Computing Qualitatively Correct Approximations of Balance Laws, Exponential-Fit, Well-Balanced and Asymptotic-Preserving[END_REF][START_REF] Waterman | Matrix exponential description of radiative transfer[END_REF]. This technical difficulty does not show up for the hyperbolic heat equation since the matrix exponential of a nilpotent matrix has only a finite number of terms.

Other directions concern the numerical development of high order discretization Finite Volume Method or of completely different discretization technics such as Finite Element Method for example. The numerical analysis of the singular case is needed to understand the influence of low wave velocities for the S n model and similar models. Since any hyperbolic system can be locally linearized, it is reasonable to extend the techniques developed in this work to more general non linear problems. In the same vein, problems with variable coefficients can be addressed systematically with this formulation.

A fully open problem is the use of the conservative structure for the convergence analysis although stability is evoked briefly in remark 7, § 4.

Figure 1 :

 1 Figure 1: Schematic of a Riemann solver. The intermediate state is computed in function of a left state L and a right state R

Figure 2 :

 2 Figure 2: Schematic of the modified Riemann solver (31). The New-left and New-right states are modification of the initial Left and Right states.

Figure 4 :

 4 Figure 4: Schematic of the modified Riemann solver (46)-(48). The New-left and New-right internal states are modification of the intermediate state by the backward propagator. This illustration is to be compared with figure 2.

  ), and -R = 1 ⊗ w -I d . The notations are 1 = (1, . . . 1) ∈ R 2n , w = (w 1 , . . . , w n , w 1 , . . . , w n ) ∈ R 2n and I d = diag(1) ∈ R n×n the identity matrix.Let us define E = diag(w i ) ∈ R n×n . The system (52) is easily rewritten in a symmetric form with the definition of U

  s j U j (t) + P x j -1 r l jr β jr = 0. (85)The associated additive formulation (25) writes s j U j (t) + r l jr A jr U jr + r P(x j ) -1 P(x r ) -I A jr U jr = 0.

Figure 5 :

 5 Figure5: Notation for nodal formulation: the corner length l jr and the corner normal n jr . Notice that l jr n jr is equal to the orthogonal vector to the half of the vector that starts at x r-1 and finish at x r+1 . The center of the cell is an arbitrary point inside the cell.

Acknowledgments: The authors warmly thank Laurent Gosse for valuable comments and discussions and for pointing out some references in the literature about well-balanced schemes for various equations.

and the commutation and composition property

The identity (78) is also

Starting from p ∈ V 2 is equivalent to the decomposition

This set yields 5 conservation laws. It means that two additional algebraic conservational laws are satisfied by the solutions of the system (78). We do not know for the moment how to use this information for the development of numerical methods, even if it seems clear that discontinuous Galerkin methods [START_REF] Hindenlang | Explicit discontinuous Galerkin methods for unsteady problems[END_REF][START_REF] Tavelli | A high order semi-implicit discontinuous Galerkin method for the two dimensional shallow water equations on staggered unstructured meshes[END_REF] may take advantage of it.


Yogesh Karpate
email: yogesh.karpate@irisa.fr

Olivier Commowick

Christian Barillot

Robust Detection of Multiple Sclerosis Lesions from Intensity-Normalized Multi-Channel MRI

Keywords: Multiple Sclerosis, Intensity Normalization, Statistsics, MRI 1. DESCRIPTION OF PURPOSE

Multiple sclerosis (MS) is a disease with heterogeneous evolution among the patients. Quantitative analysis of longitudinal Magnetic Resonance Images (MRI) provides a spatial analysis of the brain tissues which may lead to the discovery of biomarkers of disease evolution. Better understanding of the disease will lead to a better discovery of pathogenic mechanisms, allowing for patient-adapted therapeutic strategies. To characterize MS lesions, we propose a novel paradigm to detect white matter lesions based on a statistical framework. It aims at studying the benefits of using multi-channel MRI to detect statistically significant differences between each individual MS patient and a database of control subjects. This framework consists in two components. First, intensity standardization is conducted to minimize the inter-subject intensity difference arising from variability of the acquisition process and different scanners. The intensity normalization maps parameters obtained using a robust Gaussian Mixture Model (GMM) estimation not affected by the presence of MS lesions. The second part studies the comparison of multi-channel MRI of MS patients with respect to an atlas built from the control subjects, thereby allowing us to look for differences in normal appearing white matter, in and around the lesions of each patient. Experimental results demonstrate that our technique accurately detects significant differences in lesions consequently improving the results of MS lesion detection.

Multiple sclerosis (MS) is a chronic inflammatory demyelinating disease affecting the central nervous system (CNS). The disease usually affects young adults, who suffer from a variety of neurological symptoms, with periods of relapse and remission (relapsing-remitting MS), or progressively leading to irreversible disabilities (primary or secondary progressive MS). The pathophysiology of the disease is complex, and not perfectly understood at this time. Quantitative analysis of Magnetic Resonance Images (MRI) of different patients provides a spatial analysis of the brain tissues, which may lead to the discovery of putative biomarkers of disease evolution. The origin and evolution of this disease are still not well understood, and numerous studies have been conducted to evaluate its evolution and its influence on neighboring brain structures. Although manual lesion detection by experts is the Gold Standard, the objective evaluation of lesion becomes difficult for the radiologist when the number of MR sequences grows dramatically. Consequently, several studies investigated the automatic/semi-automatic segmentation of MS lesions using MRI. [START_REF] Lladó | Segmentation of multiple sclerosis lesions in brain mri: A review of automated approaches[END_REF] One breed of MS lesion segmentation includes Gaussian Mixture Modeling (GMM) on multi-spectral MR images, where each multivariate Gaussian probability density function represents a tissue: e.g. cerebrospinal fluid (CSF), gray matter (GM) and white matter (WM). The GMM enables characterization of the image intensities with a reduced number of parameters which are estimated by a maximum likelihood estimator (MLE) with the Expectation-Maximization (EM) algorithm and its variants. [START_REF] Garca-Lorenzo | Trimmed-likelihood estimation for focal lesions and tissue segmentation in multisequence mri for multiple sclerosis[END_REF] The other approach contain detection based the group based comparison between controls and MS patients. In their pioneering work of usage of Diffusion Tensor Imaging (DTI) in MS, Filipi et al. [START_REF] Filippi | Diffusion tensor magnetic resonance imaging in multiple sclerosis[END_REF] evaluated the differences in fractional anisotropy (FA) between MS patients and controls on normal appearing white matter regions. Commowick et al. [START_REF] Commowick | Detection of DTI white matter abnormalities in multiple sclerosis patients[END_REF] proposed the use of the whole diffusion tensor information (DTI) to detect statistically significant differences between MS patients and controls. We propose a framework that relies on the multi-channel MRI information to detect statistically significant differences between MS patients and controls. The proposed technique is built on ideas similar to Commowick et al. [START_REF] Commowick | Detection of DTI white matter abnormalities in multiple sclerosis patients[END_REF] but taking into account conventional multi-channel MRI. Such an approach provides consistent and reliable lesion detections. The proposed framework consists of two major parts. First, we describe a new intensity standardization technique for minimizing the inter subject intensity differences which may occur due to protocol variations across the scanners. Since in MS, White Matter (WM) lesions are also present in addition to healthy brain tissues, the intensities from multichannel MRI are modeled with a parametric transform using a robust GMM estimation based on γ divergence. Such a formulation allows for keeping the lesions unaffected. It provides a technique that (1) uses tissue-specific intensity information by modeling them using a robust GMM; (2) provides a consistent intensity normalization between inter subject images. The second part concerns the computation and use of an unbiased atlas of multi-channel MRI from a database of controls. Based on this atlas, we show how to compute statistical differences between a patient and the atlas using a Mahalanobis distance derived from the multi-channel MRI. Subsequently, we demonstrate its crucial role for further lesion analysis.

METHODOLOGY 2.1 Intensity Standardization

Given two MR images of a single MS patient at time instant t 1 and t 2 , we seek to estimate a correction factor such that corresponding anatomical tissues adopt the same intensity profile. We model the image intensities of a healthy brain with a 3-class GMM, where each Gaussian represents one of the brain tissues White Matter (WM), Gray Matter(GM) and Cerebrospinal fluid (CSF). We consider the m MR sequences as a multidimensional image with n voxels. Each voxel i is represented as

x i = [x i1 ...x im ]
. The probability of intensity x i is calculated as follows:

f (x i |θ) = 3 k=1 π k N (µ k , Σ k ) (1)
where the mean µ k and covariance Σ k define the parameters of each Gaussian of the model N (µ k , Σ k ) along with their mixing proportions π k merged into parameter θ. If the proportions were known, θ could be estimated through the Maximum Likelihood Estimator (MLE):

θ = argmax θ L(θ) = argmax θ n i=1 f (x i |θ) (2)
Where x i are considered as i.i.d. samples. However, as π k are unknown, an Expectation Maximization (EM) algorithm 5 is used to estimate the parameters.

γ-loss Function for the Normal Distribution

The parameter estimation with classic MLE for GMM can deviate from its true estimation in presence of outliers.

In MS patients, such outliers may be of crucial importance as they may denote appearing or disappearing lesions.

Recently, to address these shortcomings, a modification of the MLE was proposed in order to make it more robust to outliers. The basic idea is to maximize Equation 2 in the form of γ divergence. Notsu et al. [START_REF] Notsu | Spontaneous clustering via minimum gamma-divergence[END_REF] proposed the γ-loss function for the Normal distribution with mean vector µ and covariance matrix Σ.

L γ (µ, Σ) = Σ -γ 2(1+γ) n i=1 exp - γ 2 (x i -µ) T Σ -1 (x i -µ) (3)
Where |.| indicates the determinant. The bounded influence function of an estimator is an indicator of robustness to outliers. The influence function for GMM with γ loss function is bounded whereas the one for regular GMM is unbounded. As γ grows larger, bounds become tighter. For a sufficiently large γ, (γ ≥ 0.1), the estimating equation has little impact from outliers contaminated in the data set. Equation 3 can be casted to yield an EM style algorithm as follows.

Expectation

Step. In the case of a GMM, the latent variables are the point-to-cluster assignments k i , i = 1, ..., n, one for each of n data points. The auxiliary distribution q(k i |x i ) = q ik is a matrix with n×K entries. Each row of q i can be thought of as a vector of soft assignments of the data points x i to each of the Gaussian modes.

q ik = π k exp -γ 2 (x i -µ k ) T Σ -1 k (x i -µ k ) K l=1 π l exp -γ 2 (x i -µ l ) T Σ -1 l (x i -µ l ) (4)

Maximization

Step. The maximization step estimates the parameters of the Gaussian mixture components and the mixing proportions π k , given the auxiliary distribution on the point-to-cluster assignments computed in the expectation step. The mean µ k of a Gaussian mode is obtained as the mean of the data points assigned to it (accounting for the strength of the soft assignments). The other quantities are obtained in a similar manner, yielding to:

µ k = n i=1 q ik x i n i=1 q ik (5)
Σ k = (1 + γ) n i=1 q ik (x i -µ k )(x i -µ k ) T n i=1 q ik (6)
π k = n i=1 q ik n i=1 K l=1 q il (7)

Selection of Parameter γ

The estimation of power index γ plays a critical role in our approach, since γ affects the estimated parameters in presence of outliers. The selection of γ was studied as a model selection problem based on Akaike information criterion (AIC). [START_REF] Notsu | Spontaneous clustering via minimum gamma-divergence[END_REF] Let K be the number of clusters, p be the total numbers of parameters of a model and (µ k , Σ k ), k = 1, .., K be the means and the covariance matrices of the clusters respectively. From Equation 2, the AIC is defined as follows:

AIC γ = -2 n i=1 log f γ (x i |θ) + 2 K p(p + 3) 2 + K -1 (8)
The value of γ which minimizes AIC is used as the optimal γ. For various values of γ, Equation 8 is evaluated in cross validation manner and the γ which results in minimum value is chosen for the experiment.

Intensity Correction

We obtain the means and covariances of tissues for the source and target images using the procedure mentioned above. We chose a linear correction function such that g(x) = Σ i β i x i . The coefficients β i are estimated to minimize the following cost function: Σ l=n l=1 (g(µ source,k ) -µ target,k ) 2 . l denotes the number of voxels in class k. This function can be solved by linear regression. Using the results of the linear regression, the intensity profiles of the two images are normalized by mapping the intensity of the source image to the target image. The resulting correction function is smooth and interpolates the intensity correction.

Multiple Sclerosis Lesions (MSL) Detection

The second part of our framework relies on the voxel-wise comparison of multi-channel MRI of MS patient against a set of controls. Here, we assume that all images (patients and controls) are aligned into a common space where each voxel of each image describes the exact same spatial position. Commowick et al. [START_REF] Commowick | Detection of DTI white matter abnormalities in multiple sclerosis patients[END_REF] introduced a methodology for the comparison of diffusion tensors between a patient and a group of controls.

We aim at comparing, for each voxel, the patient vector of intensities x q and intensity vectors from control subjects x q , q = 1, ..., M . The control population X = x q , q = 1, ..., M is assumed to follow a multivariate Normal distribution N ( X, Σ X ), where X, Σ X denote respectively the average and covariance matrix of population X. Following this assumption, we compute the difference statistic between x q and N ( X, Σ X ) as a Mahalanobis distance as

d 2 (x q ) = (x q -X) T Σ -1 X (x q -X) (9)
where d 2 will vary between 0 and +∞, getting smaller as the patient vector gets likely to belong the multivariate Normal distribution of the controls. A p-value can then be computed from this distance as the statistic T = M (M -h) h(M 2 -1) d 2 follows a Fisher distribution with parameters h and M -h: T ∼ F (h, M -h)where h is the vector dimension (1 for scalar images). The difference test p-value is therefore written as:

p(x q ) = 1 -F h,M -h M (M -h) h(M 2 -1) d 2 (x q ) ( 10
)
where F h,M -h is the cumulative distribution function of a Fisher distribution with parameters h and M -h.

Finally, as we are computing voxel-wise comparisons, a final step is to correct for multiple comparisons utilizing FDR correction as introduced by Benjamini et al. [START_REF] Benjamini | Controlling the false discovery rate: A practical and powerful approach to multiple testing[END_REF]

EXPERIMENTS & RESULTS

Whole-brain MR images were acquired on 16 MS patients and 20 controls. T1-w MPRAGE, T2-w and FLAIR modalities were chosen for the experiment. Expert annotations of lesions were carried out by an expert radiologist on all MS patients. The volume size for T1-w MPRAGE and FLAIR is 256×256×160 and voxel size is 1×1×1 mm 3 . For T2-w, the volume size is 256×256×44 and voxel size is 1×1×3 mm 3 . All imaging experiments for this study were performed on a 3T Siemens Verio (VB17) scanner with a 32-channel head coil. MR images from each patient are de-noised, 8 bias field corrected [START_REF] Tustison | N4ITK: improved N3 bias correction[END_REF] and registered with respect to T1-MPRAGE volume. [START_REF] Ourselin | Block matching: A general framework to improve robustness of rigid registration of medical images[END_REF][START_REF] Commowick | Block-matching strategies for rigid registration of multimodal medical images[END_REF] All the images are processed to extract intra-cranial region. [START_REF] Smith | Fast robust automated brain extraction[END_REF] We built a geometrically unbiased atlas for each sequence from dataset of controls. [START_REF] Guimond | Average brain models: A convergence study[END_REF] For each MS patient, nonlinear registration of the T1-w MPRAGE image onto the atlas was done. [START_REF] Commowick | Automated diffeomorphic registration of anatomical structures with rigid parts: Application to dynamic cervical MRI[END_REF] Patient to atlas registration consists in two steps. First, a global affine transformation was computed between the T1 images. Then, nonlinear registration was performed to locally align the anatomies of the atlas and the patient. The final transformations were then applied to T2-w and FLAIR images.

Lesion Detection Evaluation

We built a composite vector image for each patient and control from T1-w MPRAGE, T2-w and FLAIR. The atlas is considered as the reference image to which all controls and patients images are aligned using intensity normalization as described in Section 2.1.3. Finally, we applied our framework (Section 2.2) to compare the differences in images and report the results in the next section.

We performed experiments of lesion detection, to enhance the interest of our approach with respect to (1) intensity normalization and (2) to combine several modalities for lesion detection. To this end, we compared detections using each individual sequence i.e. T1-w MPRAGE, T2-w or FLAIR individually with and without intensity normalization; and then on composite vector images formed from all sequences with and without intensity normalization. The detections obtained are compared with the ground truth using Dice scores.

Quantitative Results

We report in Table 1, the Dice score average across 16 patients for each test variant. The lesion detection on different sequences are compared against each other. We report the best mean Dice score which is of T2-w as 0.76(±0.021). From Table 1, it is clear that intensity normalization plays crucial role in our framework. The lesion detection of T2-w achieves good performance because some lesions may not be found on T1-w MPRAGE and FLAIR. In Table 1, we report the mean Dice score of 0.79(±0.020) on vector images with intensity normalization and 0.68(±0.024) without intensity normalization, which shows a clear improvement on the best detection score from individual modalities.

Qualitative Results

Figure 1 shows the lesion detection results for T2-w and composite vector image. Each row represents the image and each column depicts the anatomical slice of patient MRI, corresponding ground truth, lesion detection results without and with intensity normalization respectively. The white labels show the true lesions detected while red ones indicate the brain tissues which are detected as lesion but actually are not part of the lesions. This figure demonstrates visually the ability of our approach to detect lesions. As seen from the last two columns, there is considerable improvement of lesion detection, thanks to intensity normalization.

NEW OR BREAKTHROUGH WORK TO BE PRESENTED

A novel paradigm for detection of MS lesions has been proposed. We devised a two layers strategy, first intensity standardization of images with respect to an atlas; and then computation of difference scores of a specific patient to a population of controls. Such comparisons allowed the accurate detection of lesions or other diffuse pathology-related changes.

CONCLUSION

We proposed a new framework for MS lesion detection based on differences between multi-channel MRI of patients and controls. This framework fused with intensity standardization was applied to the study of MS and highlighted the great interest of quantitative MRI measurements made on MRI for a better understanding and characterization of the disease. The efficacy of our method was evaluated though detection with and without intensity correction. Compared to un-normalized images, detection with intensity correction is a better choice for MS lesion analysis thanks to its ability to preserve the intensity variations caused by pathological changes while normalizing healthy tissue intensities. The resulting system is both efficient and accurate. This performance suggests that it can provide valuable assistance in detecting MS lesions in clinical routine with high reliability.

Our models are already capable of detecting highly variable lesion patterns, but we would like to move towards richer models. The framework described here allows for exploration of additional MR sequences with or without contrast agents. For example, one can consider infusing T1-w Gadolinium.

Table 1 :Figure 1 :

 11 Figure 1: Top row : A slice of MRI from T2-w Sequence (a), MSL detection obtained by without intensity normalization (b), MSL detection with intensity normalization (c) and corresponding ground truth (d). Bottom row : A slice of MRI from composite vector image (e), MSL detection obtained without intensity normalization (f), MSL detection with intensity normalization (g) and corresponding ground truth (h). White label shows the true positives where as red label indicates the false positives.

ACKNOWLEDGMENTS

This work was supported by the Region Council of Brittany of France and France Life Imaging (FLI) through the ANR-11-INBS-00006 grant.