T Pichard
email: pichard@mathcces.rwth-aachen.de

G W Alldredge

Stéphane Brull

B Dubroca

M Frank

M Frank Mathcces

An approximation of the M 2 closure: application to radiotherapy dose simulation

Keywords: Moment models, Entropy-based closure, Radiotherapy dose computation Mathematics Subject Classification (2000) 82C70, 82D99, 94A17

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

The aim of radiation treatments is to destroy tumor cells by prescribing a certain quantity of energy, called the dose, to the tumor cells. This dose is produced by radiation which can be modelled by the transport of particles (photons, electrons, protons, hadrons, depending on the type of radiation).

A large range of numerical approaches has been proposed in the literature to compute the dose. Dose distributions are typically numerically computed using Monte-Carlo algorithms (see e.g. [START_REF] Sempau | PENELOPE-2011: A Code System for Monte Carlo Simulation of Electron and Photon Transport[END_REF]) or discrete-ordinate methods (see e.g. [START_REF] Lewis | Computational methods of neutron transport[END_REF]). However such direct solution methods often require more computing resources than are typically available in medical centers. Resource effective alternatives to those approaches (see e.g. [START_REF] Mayles | Handbook of radiotherapy physics: Theory and practice[END_REF] and references therein), used in medical centers, include semi-empirical methods (e.g. Fermi-Eyges methods), probabilistic methods (fast Monte Carlo simulations, see e.g. [START_REF] Spezi | An overview of Monte Carlo treatment planning for radiotherapy[END_REF][START_REF] Chetty | Report of the AAPM task group no. 105: Issues associated with clinical implementation of Monte Carlo-based photon and electron external beam treatment planning[END_REF][START_REF] Zankowski | Fast electron monte carlo for eclipse TM[END_REF] and references therein) and PDE-based methods ( [START_REF] Wareing | Acuros XB advanced dose calculation for the Eclipse TM treatement planning system[END_REF]). However those alternatives may be unprecise for certain applications. The present method is an in-between alternative with an accuracy comparable to Monte-Carlo and much lower computational costs.

This paper is a follow-up to [START_REF] Duclous | A deterministic partial differential equation model for dose calculation in electron radiotherapy[END_REF][START_REF] Pichard | Relaxation schemes for the M 1 model with space-dependent flux: application to radiotherapy dose calculation[END_REF][START_REF] Caron | Deterministic model for the transport of energetic particles. application in the electron radiotherapy[END_REF]. The aim is to propose a PDE-based numerical approach which is fast and accurate enough for practical applications. We study in the next section a moment approach, i.e. a PDE-based approach, which has a much lower numerical cost than the Monte-Carlo methods.

The transport of particles for radiotherapy problem can be described by kinetic models ( [START_REF] Hensel | Deterministic model for dose calculation in photon radiotherapy[END_REF][START_REF] Olbrant | Generalized Fokker-Planck theory for electron and photon transport in biological tissues: application to radiotherapy[END_REF]) for the fluence ψ of the particles. Due to the high dimensionality of the fluence (it depends on position x ∈ R 3 , energy ε ∈ R + and direction of flight Ω ∈ S 2 where S 2 is the unit sphere), solving the kinetic transport equations is numerically expensive. The kinetic model can be reduced by extracting angular moments. The resulting models retain the major properties of the kinetic models.

Those models were applied in a large range of physics, such as fluid dynamics ([22,34, 41]), plasma physics ( [START_REF] Guisset | Asymptotic-preserving scheme for the Fokker-Planck-Landau-Maxwell system in the quasi-neutral regime[END_REF][START_REF] Guisset | Limits of the M 1 and M 2 angular moments models for kinetic plasma physics studies[END_REF][START_REF] Mallet | An entropic scheme for an angular moment model for the classical Fokker-Planck-Landau equation of electrons[END_REF][START_REF] Mallet | General moment system for plasma physics based on minimum entropy principle[END_REF]), semi-conductors ( [START_REF] Anile | Thermodynamic derivation of the hydrodynamical model for charge transport in semiconductors[END_REF][START_REF] Rosa | Exact maximum entropy closure of the hydrodynamical model for Si semiconductors: The 8-moment case[END_REF]) or radiative transfer ( [START_REF] Chandrasekhar | Radiative transfer[END_REF][START_REF] Dubroca | Hiérarchie des modèles aux moments pour le transfert radiatif[END_REF][START_REF] Levermore | Relating Eddington factors to flux limiters[END_REF][START_REF] Berthon | An HLLC scheme to solve the M 1 model of radiative transfer in two space dimensions[END_REF][START_REF] Brunner | One-dimensional riemann solvers and the maximum entropy closure[END_REF]).

The main difficulty arising when deriving moment models is computing a closure. Indeed moments equations have more unknowns than equations, therefore a closure needs to be computed. This closure is generally chosen to retrieve the basic features of the underlying kinetic models. Generally, one constructs an ansatz ψ R for the exact fluence so that the ansatz satisfies certain integral constraints. The closure is then computed by replacing the exact fluence with the ansatz.

In the present paper, an entropy-based closure is chosen which leads to the hierarchy of moment models colloquially known as M N , where N indicates the highest order of the moments in the model. This closure is based on the physics of the collisions and leads to a hyperbolic system of moment equations with an entropy dissipation property ( [START_REF] Levermore | Moment closure hierarchies for kinetic theories[END_REF]). However, computing such a closure directly requires solving (numerically) a minimization problem at every point in space and energy. Furthermore, the solution of this minimization problem requires repeated expensive quadrature computations ([27,2,1]). The main goal of this paper is to provide an approximation of the M 1 and the M 2 closures which avoids the optimization problem and the quadrature computations therein. We obtain the approximation by a careful study of the domain of realizability and by invariance properties of the entropy minimizer. There exist many approximate moment closures, e.g. the main alternative being the P N closures. However, they present drawbacks conflicting with the applications we have in mind. Especially, several of them are only derived in one space dimension (see e.g. [START_REF] Schneider | Kershaw closures for linear transport equations in slab geometry I: Model derivation[END_REF][START_REF] Vikas | Radiation transport modeling using extended quadrature method of moments[END_REF][START_REF] Alldredge | Approximating the M 2 method by the extended quadrature method of moments for radiative transfer in slab geometry[END_REF]), and it is not entirely clear how to generalize the ideas to multiple dimensions. In this paper we deal with the case of multiple dimensions.

In the next section, a simplified kinetic model of the transport of electrons is described, the procedure of moment extraction is presented and illustrated through the M 1 model and the advantages of the M 2 model are shown. The main result is presented in Section 3, it consists of an approximation of the M 2 closure for three dimensional problems. In order to validate our approach, the kinetic and moment models are compared on numerical test cases from medical applications in Section 4. The last section is devoted to conclusions and perspectives.

Models

The transport of electrons in the field of radiotherapy can be modeled using kinetic theory. We first recall a kinetic model used in the field of medical physics ( [START_REF] Hensel | Deterministic model for dose calculation in photon radiotherapy[END_REF][START_REF] Duclous | A deterministic partial differential equation model for dose calculation in electron radiotherapy[END_REF]), then the procedure of moment extraction, which is afterward illustrated through the first model in the hierarchy of entropy based moment models, the M 1 model.

Kinetic model

In this study, only the transport of electrons is considered. The transport of photons with the approach described in this paper has been studied in [START_REF] Pichard | The M 2 model for dose simulation in radiation therapy[END_REF]. We consider here only Møller's and Mott's collisions for electrons (see [START_REF] Hensel | Deterministic model for dose calculation in photon radiotherapy[END_REF] and references therein). Mott's collisions are elastic. Møller's collisions are inelastic and are ionizing interactions, i.e. two electrons emerge from this collision, a primary (the more energetic one) and a secondary. The electron transport can be modelled by the following kinetic equation [START_REF] Hensel | Deterministic model for dose calculation in photon radiotherapy[END_REF][START_REF] Duclous | A deterministic partial differential equation model for dose calculation in electron radiotherapy[END_REF] Ω .∇ x ψ(x, Ω , ε) = ρ(x)Q(ψ)(x, Ω , ε).

(

The unknown ψ is the fluence of electrons depending on position x ∈ Z, energy ε ∈ [0, ε max ], direction of flight Ω ∈ S 2 , and Q is the collision operator. For Mott's and Møller's collisions, the collision can be represented by linear Boltzmann (LB) gain and loss terms

Q(ψ)(x, Ω , ε) = S 2 ∞ ε (σ M,1 + σ M,2 )(ε ′ , ε, Ω ′ .Ω )ψ(x, Ω ′ , ε ′ )dε ′ dΩ ′ -σ T,M (ε)ψ(x, Ω , ε) + S 2 σ Mott (ε, Ω ′ .Ω )ψ(x, Ω ′ , ε)dΩ ′ -σ T,Mott (ε)ψ(x, Ω , ε). (2)
The gain terms are characterized by the differential cross sections for Mott's and Møller's primary and secondary electrons σ Mott , σ M,1 , σ M,2 , and the loss terms by the total cross sections σ M,T and σ T,Mott . The superscript ′ refers to the state of the particle before collision; the absence of this superscript refers to the post-collisional state. ρ(x) is the density of atomic cores in the medium at position x. Møller's cross section for primary electrons σ M,1 is very peaked in energy, meaning that most of the particles lose small energy during those collisions. Due to that peak in the cross section, the continuous-slowing down approximation (CSD, see [START_REF] Pomraning | The Fokker-Planck operator as an asymptotic limit[END_REF][START_REF] Olbrant | Generalized Fokker-Planck theory for electron and photon transport in biological tissues: application to radiotherapy[END_REF]) can be applied here, and the deflection due to Møller's effect is negligible compared to the one due to Mott's effect. This leads to approximating ( [START_REF] Larsen | Electron dose calculations using the method of moments[END_REF][START_REF] Duclous | A deterministic partial differential equation model for dose calculation in electron radiotherapy[END_REF][START_REF] Olbrant | Generalized Fokker-Planck theory for electron and photon transport in biological tissues: application to radiotherapy[END_REF])

Q(ψ)(x, Ω , ε) ≈ Q CSD (ψ)(x, Ω , ε) = ∂ ε (Sψ)(x, Ω , ε) (3)
+ S 2 ∞ ε σ M,2 (ε ′ , ε, Ω ′ .Ω )ψ(x, Ω ′ , ε ′ )dε ′ dΩ ′ + S 2 σ Mott (ε, Ω ′ .Ω )ψ(x, Ω ′ , ε)dΩ ′ -σ T,Mott (ε)ψ(x, Ω , ε),
where the stopping power S characterizes this energy loss.

Similarily, the remaining elastic cross section σ Mott is forward-peaked, meaning that most of the particles are slightly deflected during those collisions. Due to that peak in the cross section, the collision operator can be approximated by a Fokker-Planck (FP, see [START_REF] Pomraning | The Fokker-Planck operator as an asymptotic limit[END_REF] although the validity of this approximation was discussed in [START_REF] Olbrant | Generalized Fokker-Planck theory for electron and photon transport in biological tissues: application to radiotherapy[END_REF]) operator

Q CSD (ψ)(x, Ω , ε) ≈ Q FP (ψ)(x, Ω , ε) = ∂ ε (Sψ)(x, Ω , ε) +T (ε) ∂ µ (1 -µ 2 )∂ µ ψ (x, Ω , ε) + 1 1 -µ 2 ∂ 2 φ ψ(x, Ω , ε) (4) + S 2 ∞ ε σ M,2 (ε ′ , ε, Ω ′ .Ω )ψ(x, Ω ′ , ε ′ )dε ′ dΩ ′ ,
where the transport coefficient T (ε) characterizes this deflection, µ and φ are such that Ω = (µ, 1 -µ 2 cos φ , 1 -µ 2 sin φ ).

Moment models

A moment model is a reduction of a kinetic model that requires much lower computational times (see e.g. comparisons in Section 4). One can reduce the number of variables by approximating the angular distribution by an ansatz ψ ≈ ψ R (Ω ) satisfying integral constraints.

Let us define ψ i as the moment of order i of ψ, namely

ψ i := Ω ⊗ • • • ⊗ Ω i times ψ = S 2 Ω ⊗ • • • ⊗ Ω ψdΩ .
Here ⊗ denotes the tensor product. Instead of working with ψ which depends on ε, x and Ω , the moments ψ i of order 0 to N are studied.

One obtains equations for the moments ψ i by extracting moments of (1)

∇ x .ψ i+1 (x, ε) = ρ(x)Q i (ψ i )(x, ε), (5)
where the moments of order 0 to 2 of the collision operator read

Q 0 (ψ 0 )(x, ε) = ∞ ε σ 0 M (ε ′ , ε)ψ 0 (x, ε ′ )dε ′ -(σ T,M + σ T,Mott -σ 0 Mott )(ε), (6a)
Q 1 (ψ 1 )(x, ε) = ∞ ε σ 1 M (ε ′ , ε)ψ 1 (x, ε ′ )dε ′ -(σ T,M + σ T,Mott -σ 1 Mott )(ε)ψ 1 (x, ε ′ ) (6b) Q 2 (ψ 2 )(x, ε) = ∞ ε σ 0 M -σ 2 M 2 (ε ′ , ε)tr(ψ 2 )(x, ε ′ )Id + 3σ 2 M -σ 0 M 2 (ε ′ , ε)ψ 2 (x, ε ′ )dε ′ (6c) + σ 0 Mott -σ 2 Mott 2 (ε)tr(ψ 2 )(x, ε)Id + 3σ 2 Mott -σ 0 Mott 2 (ε)ψ 2 (x, ε) -(σ T,M + σ T,Mott )(ε)ψ 1 (x, ε ′ )
where σ M = σ M,1 + σ M,2 and σ i are scalars given by

σ i (ε ′ , ε) = 2π +1 -1 µ i σ (ε ′ , ε, µ)dµ.
Similarily the moments of the CSD operator (3) and the FP operator (4) yield

Q 0 CSD (ψ 0 )(x, ε) =∂ ε (Sψ 0 )(x, ε) (7a) + ∞ ε σ 0 M,2 (ε ′ , ε)ψ 0 (x, ε ′ )dε ′ -(σ T,Mott -σ 0 Mott )(ε)ψ 0 (x, ε), Q 1 CSD (ψ 1 )(x, ε) =∂ ε (Sψ 1 )(x, ε) (7b) + ∞ ε σ 1 M,2 (ε ′ , ε)ψ 1 (x, ε ′ )dε ′ -(σ T,Mott -σ 1 Mott )(ε)ψ 1 (x, ε), Q 2 CSD (ψ 2 )(x, ε) =∂ ε (Sψ 2 )(x, ε) (7c) + ∞ ε σ 0 M,2 -σ 2 M,2 2 (ε ′ , ε)tr(ψ 2 )(x, ε ′ )Id + 3σ 2 M,2 -σ 0 M,2 2 (ε ′ , ε)ψ 2 (x, ε ′ )dε ′ + σ 0 Mott -σ 2 Mott 2 (ε)tr(ψ 2 )(x, ε)Id + 3σ 2 Mott -σ 0 Mott 2 (ε)ψ 2 (x, ε) -σ T,Mott (ε)ψ 2 (x, ε), Q 0 FP (ψ 0 )(x, ε) =∂ ε (Sψ 0 )(x, ε) + ∞ ε σ 0 M,2 (ε ′ , ε)ψ 0 (x, ε ′ )dε ′ , (8a)
Q 1 FP (ψ 1 )(x, ε) =∂ ε (Sψ 1 )(x, ε) -2T (ε)ψ 1 (x, ε) + ∞ ε σ 1 M,2 (ε ′ , ε)ψ 1 (x, ε ′ )dε ′ , (8b)
Q 2 FP (ψ 2 )(x, ε) =∂ ε (Sψ 2 )(x, ε) -2T (ε) 3ψ 2 (x, ε) -tr(ψ 2 )(x, ε)Id (8c) + ∞ ε σ 0 M,2 -σ 2 M,2 2 (ε ′ , ε)tr(ψ 2 )(x, ε ′ )Id + 3σ 2 M,2 -σ 0 M,2 2 (ε ′ , ε)ψ 2 (x, ε ′ )dε ′ .
Remark that the only difference between the moments of the CSD and FP operators are the scalars before ψ 2 (x, ε) and tr(ψ 2 )(x, ε)Id.

The system (5) requires a closure, as it has more unknown than equations. In practice, this is done by approximating ψ by an ansatz ψ R , and then computing the higher-order term using this ansatz, i.e.

ψ N+1 ≈ Ω ⊗ • • • ⊗ Ω N+1 times ψ R ,
where ψ R is an ansatz having the moments (ψ 0 , ..., ψ N ).

For N = 2, choosing an ansatz ψ R provides an approximation of the flux ψ 3 depending on (ψ 0 , ψ 1 , ψ 2 ). For each possible set of moments (ψ 0 , ψ 1 , ψ 2 ), one needs to find an ansatz ψ R having the right moments. This problem can be written as

find ψ R (Ω ), s.t. Ω ⊗ • • • ⊗ Ω i times ψ R = ψ i , i = 0, ..., 2.
One possibility is simply to choose the ansatz ψ R as a polynomial of degree

N ψ R (Ω ) = λ . m(Ω ).
Here m(Ω ) is a basis of polynomial of degree N over S 2 , in particular, we chose

for N = 1, m(Ω ) = (1, Ω 1 , Ω 2 , Ω 3 ) , for N = 2, m(Ω ) = Ω 1 , Ω 2 , Ω 3 , Ω 2 1 , Ω 2 2 , Ω 2 3 , Ω 1 Ω 2 , Ω 1 Ω 3 , Ω 2 Ω 3
in the rest of the paper and λ is the unique vector of scalars such that the moments of ψ R are (ψ 0 , ..., ψ N ). This leads to the so-called P N models. Computing the ansatz (and then the closure) is simple and only requires the solution of a linear systems. However, P N models present several drawbacks. First the obtained function ψ R can be negative (see e.g. [START_REF] Hauck | Positive P N closures[END_REF][START_REF] Mcclarren | Robust and accurate filtered spherical harmonics expansion for radiative transfer[END_REF]), which is not physical. Second, in radiotherapy beams of particles are used. Perfect beams can be modeled by Dirac distribution in Ω . Such distributions are poorly approximated by polynomials.

Example 1 Consider a Dirac peak ψ = δ (Ω .e 1 -1). Its first moments read

ψ 0 = ψ = 1, ψ 1 = Ω ψ = e 1 .
The polynomial ψ R ∈ R 1 [X] of degree 1 having ψ 0 and ψ 1 for moments reads

ψ R (Ω ) = 1 + 3Ω 1 2 .
One first remark that the P 1 ansatz differs from the beam distribution. Then computing the second order moment of this ansatz, i.e. the P 1 closure, read

Ω ⊗ Ω ψ R (Ω ) = Diag 1 3 , 1 3 , 1 3 = e 1 ⊗ e 1 = Ω ⊗ Ω ψ ,
which differs from the second order moment of the Dirac peak, and therefore the P 1 model is not able to represent a beam. Similar computations show that for any N the P N model is also not able to capture the exact closure. One solution to obtain a good accuracy with the P N model when considering beams of particles consists in chosing the number of moments N large, which deteriorates the time efficiency of the method.

Among the possible candidates for ψ R (having the moments (ψ 0 , ..., ψ N )), we chose the one that minimizes Boltzmann entropy function

H ( f ) = S 2 ( f log( f ) -f ) (Ω )dΩ , (9)
which leads to the hierarchy of the so-called M N models.

Theorem 1 ([42,8,9,10,30,51,28]) Consider a vector of polynomials m(Ω ) and a vector ψ such that there exists at least one positive function ψ > 0 satisfying

ψ = mψ .
Then there exists a unique minimizer ψ R to [START_REF] Borwein | Partially finite convex programming, Part I: Quasi relative interiors and duality theory[END_REF] which has the form ψ R (Ω ) = exp( λ . m(Ω )).

Furthermore, the function that sends the moments ψ onto the ansatz (10) is a smooth bijection.

Beam-like distributions can be correctly approximated by ansätze of the form [START_REF] Borwein | Partially finite convex programming: Part II[END_REF]. Indeed, a perfect beam in the direction Ω 0 can be modeled by the Dirac distribution, which can be interpreted as the limit of a sequence of M 1 ansätze

δ (Ω .Ω 0 -1) = lim n→+∞ 1 2π exp ((log(n) -n) + nΩ .Ω 0 ) .
In this sense beam-like distributions are in the closure of the set of distributions of the form [START_REF] Borwein | Partially finite convex programming: Part II[END_REF].

In the next section, we focus on the following two properties of the M N model (and especially of the M 2 model) when constructing the approximation of the M 2 closure.

Hyperbolicity The M N equations are known to be symmetric hyperbolic ( [START_REF] Levermore | Moment closure hierarchies for kinetic theories[END_REF][START_REF] Friedriechs | Systems of conservation equations with a convex extension[END_REF]). This means that the Jacobian of the flux J = ∇ (ψ 0 ,...,ψ N ) (ψ 1 , ..., ψ N+1 ) is diagonalizable with real eigenvalues.

Realizability

The solution ψ to the kinetic equation 1, with the collision operator (2), (3) or (4), is positive. The realizability property corresponds to requiring that the solution (ψ 0 , ..., ψ N ) to the moments system [START_REF] Berntsen | Algorithm 698: Dcuhre: An adaptive multidemensional integration routine for a vector of integrals[END_REF] with the associated collision operator ( 6), [START_REF] Berthon | Numerical methods for balance laws with space dependent flux: application to radiotherapy dose calculation[END_REF] or [START_REF] Borwein | Duality relationships for entropy-like minimization problems[END_REF] are the moments of one positive distribution ψ, i.e.

∃ψ(Ω ) > 0, s.t. (ψ 0 , ..., ψ N ) = S 2 (1, ..., Ω ⊗ • • • ⊗ Ω N times )ψ(Ω )dΩ .

Definition 1

The set of all realizable moments of order up to N is called the realizability domain R N of order N

R N =    S 2 (1, ..., Ω ⊗ • • • ⊗ Ω N times )ψ(Ω )dΩ , ∀ψ > 0    . ( 11
)
Remark 1 The notion of realizability can be extended to general measures. In this paper, we avoid the technicalities associated to measures and always write non-negative densities, which formally includes the case of a sum of Diracs to represent a discrete measure.

Remark 2

The realizability domain is a convex cone. This means that positive combinations of two realizable moments is realizable. Indeed suppose α 1 > 0 and α 2 > 0 are two positive scalars and ψ 1 and ψ 2 are two positive functions of Ω and

ψ1 = (1, ..., Ω ⊗ • • • ⊗ Ω N times )ψ 1 , ψ2 = (1, ..., Ω ⊗ • • • ⊗ Ω N times )ψ 2 .
Then the sum

α 1 ψ1 + α 2 ψ2 = (1, ..., Ω ⊗ • • • ⊗ Ω N times )(α 1 ψ 1 + α 2 ψ 2 )
is obviously realizable since α 1 ψ 1 + α 2 ψ 2 is positive. This property will be used below in Subsection 3.3 to enforce realizability of the approximated closure.

The coefficients of λ are generally determined by solving the dual minimization problem ([27,2,1])

λ = argmin ā exp( ā. m(Ω )) -ā. ψ, ( 12
)
where ψ = m(Ω )ψ are the moments ψ i associated to m(Ω ) arranged as a column vector. Solving this minimization problem is however computationally expensive (see e.g. comparison in Section 4). Therefore we propose an alternative for the first two models in the hierarchy, the M 1 and M 2 models. It consists in smooth approximations preserving several important properties of the exact closure.

The M 1 model

We illustrate here the procedure described above for the first model in the hierarchy and propose an approximation of the closure. Considering only equations [START_REF] Berntsen | Algorithm 698: Dcuhre: An adaptive multidemensional integration routine for a vector of integrals[END_REF] for i = 1, 2, we need to express ψ 2 as a function of ψ 0 and ψ 1 . In that case, using (10)

, ψ R yields ψ R (Ω ) = exp(λ 0 + λ 1 Ω 1 + λ 2 Ω 2 + λ 3 Ω 3 ). ( 13
)
Then the coefficients λ 0 ∈ R and (λ 1 , λ 2 , λ 3 ) ∈ R 3 are such that ψ R has the right moments

ψ R (Ω ) = ψ 0 , Ω ψ R (Ω ) = ψ 1 .
Let us write n = ψ 1 /|ψ 1 |. Using rotational invariance ( [START_REF] Levermore | Relating Eddington factors to flux limiters[END_REF]), one can prove that there exists

α ∈ R such that (λ 1 , λ 2 , λ 3 ) = αn.
Using this form of the ansatz ψ R leads to

||ψ 1 || 2 ψ 0 = |α|coth(|α|) -1 |α| = f (|α|). (14)
This function f is bijection between R + and [0, 1[. Computing the closure ψ 2 leads to

ψ 2 = ψ 0 1 -χ 2 2 Id + 3χ 2 -1 2 n ⊗ n , (15)
χ 2 ||ψ 1 || 2 ψ 0 = +1 -1 µ 2 exp(α µ)dµ +1 -1 exp(α µ)dµ = 1 + 2 |α| (1 + coth (|α|)) ,
where |α| is a function of ||ψ 1 || 2 ψ 0 given by ( 14), i.e.

|α| = f -1 ||ψ 1 || 2 ψ 0 .
The Eddington factor χ 2 does not have an analytical formula but can be approximated (see below in Section 3.3).

The M 1 model is often used because it is simple to implement and covers a large range of physical phenomena. However it presents also several drawbacks. These drawbacks are presented in the next section in order to highlight the advantages of the M 2 model compared to the M 1 model. Then an approximation of the M 2 closure is proposed. Section 4 illustrates the efficiency of the M 2 model through numerical test cases.

The advantages of the M 2 model

Now we motivate the advance to the M 2 model, that is, increasing the moment order to two. The M 2 model is able to model a larger range of physical phenomena than the M 1 model.

First, the physics are modeled more accurately. Indeed decomposing the differential cross sections into polynomials reads

σ (ε ′ , ε, µ) = ∞ ∑ i=0 σ i (ε ′ , ε)µ i .
Then extracting the N first moments of the collisional operator is equivalent to truncating this expansion at degree N. So clearly the collisions are better modelled as N increases. This phenomenon is illustrated through the test cases in Sections 4.2 and 4.4 below.

Furthermore, M 1 is not able to distinguish certain multiple-beam cases which are of major importance in the field of radiotherapy. Indeed in external radiotherapy, the source of particles creating the dose are beams of particles applied on the boundary of a medium (see e.g. the numerical test cases in Section 4).

Example 2 Let us consider two perfect beams of opposite direction ±e 1 crossing each other. This is modelled by a distribution composed of two Dirac peaks

ψ = δ (Ω .e 1 -1) + δ (Ω .(-e 1 ) -1).
Extracting the moments of this distribution yields

ψ 0 = ψ = 2, ψ 1 = Ω ψ = 0 R 3 , ψ 2 = Ω ⊗ Ω ψ = 2e 1 ⊗ e 1 .
Working with the M 1 model means working with (ψ 0 , ψ 1 ). In that case, the first two mo- ments are the same of those of an isotropic distribution ψ = 1 2π . This means that the M 1 model is unable to distinguish two beams from a isotropic distribution. This produces an overestimation of the diffusion at the point where the beams cross each other. However with ψ 2 available, the M 2 model is able to recognize that the underlying distribution is not isotropic.

This problem also appears more generally for two-beam distributions.

ψ 0 = ψ = 2, ψ 1 = Ω ψ = e 1 + e 2 , ψ 2 = Ω ⊗ Ω ψ = e 1 ⊗ e 1 + e 2 ⊗ e 2 .
Here the M 1 model sees a single beam in the direction e 1 + e 2 . But because ψ 2 does not have the moments of a single beam at e 1 + e 2 (these would be (e 1 + e 2 ) ⊗ (e 1 + e 2 )), the M 2 ansatz can distinguish these cases.

Remark 3 This problem could be circumvented by exploiting the linearity of the underlying kinetic equation [START_REF] Alldredge | Adaptive change of basis in entropy-based moment closures for linear kinetic equations[END_REF][START_REF] Alldredge | High-order entropy-based closures for linear transport in slab geometry II: A comutational study of the optimization problem[END_REF][START_REF] Alldredge | Approximating the M 2 method by the extended quadrature method of moments for radiative transfer in slab geometry[END_REF][START_REF] Anile | Thermodynamic derivation of the hydrodynamical model for charge transport in semiconductors[END_REF]. Indeed, suppose that for each i, ψ i is the solution of (1) with the collision operator (2), ( 3) or ( 4) with a boundary condition ψ| ∂ Z = ψ b i , where ψ b i contains (for example) only a single i-th beam entering the domain, and an initial condition of zero. Then the solution of (1-4) with the boundary condition ψ| ∂ Z = ∑ i ψ b i and the same initial condition would be ∑ i ψ i . Each solution ψ i can be approximated by solving the M N system (5) with the associated collision operator. This way one obtains an approximation of the solution ∑ i ψ i with multiple beams by solving the M N system for each beam separately.

Approximation of the M 2 closure

The computation of the M 2 closure in three dimensions of space cannot be simplified to one dimension by using symmetry arguments as in the M 1 case (Section 2.3), therefore typically it would be computed by numerically solving the minimization problem [START_REF] Caron | Deterministic model for the transport of energetic particles. application in the electron radiotherapy[END_REF]. However, this method would be too expensive for the applications we have in mind, so we propose an approximation which takes advantage of the underlying hierarchical structure of the closure and is exact for certain special boundary cases.

A naive idea to approximate the M 2 closure consists in solving [START_REF] Caron | Deterministic model for the transport of energetic particles. application in the electron radiotherapy[END_REF] for a large amount of values of (ψ 0 , ψ 1 , ψ 2 ) and directly constructing a polynomial fitting those values. Such an approximation is not very satisfactory as it does not satisfy basic properties. For instance, when choosing a 1D function ψ(Ω 1 ), the closure ψ 3 should be the third order moment of a 1D distribution ψ R (Ω 1 ). This would not be enforced by such a naive polynomial approxi- mation.

Instead we propose the following: First we use the same method as in Section 2.3, to approximate the closure ψ 3 in the M 1 case. Then we extend this progressively to more general cases until obtaining an approximation of the general M 2 closure.

In the first subsection, we introduce the realizability domain for the M 2 model, which is the domain of physically relevant moments. In the second subsection, an approximation of the closure is proposed in some subsets of the realizability domain. In the third subsection, a polynomial fit is proposed which is correct in these subsets. Finally, we consider the properties of the obtained approximation.

The realizability domain for the M 2 model

The solution of the kinetic equation (1-4) is positive. This implies that its moments (11)). The solution of the moment system (5-8) with the M 2 closure needs to be inside this set. Indeed the M 2 closure exists only if there exists an ansatz of the form [START_REF] Borwein | Partially finite convex programming: Part II[END_REF] whose moments are (ψ 0 , ψ 1 , ψ 2 ). For moments on the unit sphere S 2 , this is equivalent (according to Theorem 1) to requiring

(1, Ω , Ω ⊗ Ω )ψ evolve in the realizability domain R 2 ⊂ R × R 3 × R 3×3 (defined in
(ψ 0 , ψ 1 , ψ 2 ) ∈ R 2 .
The realizability domain R 2 for moments of order up to 2 is characterized as follows.

Proposition 1

The realizability domain for second order moments can be written

R 2 = (ψ 0 , ψ 1 , ψ 2 ) ∈ R × R 3 × R 3×3 , s.t. tr(ψ 2 ) = ψ 0 > 0 ( 16
)
and ψ 0 ψ 2 -ψ 1 ⊗ ψ 1 is symmetric positive definite (s.p.d.) .
Proof This results follows directly from [START_REF] Kershaw | Flux limiting nature's own way[END_REF] which provided the following result

R2 = (ψ 0 , ψ 1 , ψ 2 ) ∈ R × R 3 × R 3×3 , s.t. tr(ψ 2 ) = ψ 0 ≥ 0 and ψ 0 ψ 2 -ψ 1 ⊗ ψ 1 is symmetric non-negative .
Then, using Theorem 1, one obtains the existence of a positive (instead of non-negative for the previous characterization) representing distribution for all moments (ψ 0 , ψ 1 , ψ 2 ) ∈ int(R 2 ) in the interior of the realizability domain. ( 16) is a characterization of this interior.

Finally, by replacing one of the inequality in ( 16) by an equality, one can prove that there is no strictly positive representing distribution for such moments (see e.g. [START_REF] Curto | Recusiveness, positivity, and truncated moment problems[END_REF][START_REF] Curto | A duality prood to Tchakaloff's theorem[END_REF]), which prove the equality [START_REF] Curto | A duality prood to Tchakaloff's theorem[END_REF].

In order to simplify the computation of ψ 3 , we use a transformation of the realizability domain. It consists of a normalization and a rotation.

Notation 1 -N i is the i-th moment normalized by the zeroth-order moment, i.e.

N i := ψ i ψ 0 .
-Let R be the rotation matrix that diagonalizes

ψ 0 ψ 2 -ψ 1 ⊗ ψ 1 . This also diagonalizes N 2 -N 1 ⊗ N 1 . -R T

is the set of realizable moments after these transformations (rotation and normalization)

R T := (N 1 , N 2 ), s.t. tr(N 2 ) = 1 and N 2 -N 1 ⊗ N 1 diagonal positive .
One passes easily from one set to the other

R 2 ∋ ψ 0 , ψ 1 , ψ 2 = ψ 0 , ψ 0 R.N 1 , ψ 0 R.N 2 .R T , with (N 1 , N 2 ) ∈ R T . (17)
Similarly, applying this transformation to ψ 3 yields

ψ 3 = ψ 0 Rot 3 (R, N 3 ), Rot 3 (R, N 3 ) = 3 ∑ i ′ =1 3 ∑ j ′ =1 3 ∑ k ′ =1 R i,i ′ R j, j ′ R k,k ′ N 3 i ′ , j ′ ,k ′ , (18)
where Rot 3 (R, N 3 ) is the tensor N 3 rotated using the matrix R. After transformation, N 3 only depends on

(N 1 , N 2 ) ∈ R T . Notation 2 For (N 1 , N 2 ) ∈ R T , as tr(N 2 ) = 1 we have tr(N 2 -N 1 ⊗ N 1 ) = 1 -||N 1 || 2 2 .
One can rewrite

N 2 = N 1 ⊗ N 1 + (1 -||N 1 || 2 2 )Diag(γ 1 , γ 2 , 1 -γ 1 -γ 2 ).
Then we can parametrize R T by

R P = (N 1 , γ 1 , γ 2 ) ∈ B(0 R 3 , 1)×]0, 1[×]0, 1 -γ 1 [ .
In the next subsection, we exhibit values of N 3 in subsets of R P .

Special values of the closure

Now we define a hierarchy of subdomains of R T and compute N 3 in each of them. Then we recall results about the boundary of R T .

N 3 in subdomains of R T

First we consider the relation between moments (N 1 , N 2 ) ∈ R T and functions of the form [START_REF] Borwein | Partially finite convex programming: Part II[END_REF] and vice versa. From these relations, we can give the form of the closure N 3 in certain subsets of R T .

Notation 3

We use the following monomial basis for polynomials up to degree two over the unit sphere

m(Ω ) = Ω 1 , Ω 2 , Ω 3 , Ω 2 1 , Ω 2 2 , Ω 2 3 , Ω 1 Ω 2 , Ω 1 Ω 3 , Ω 2 Ω 3 ,
and we write

ψ R (Ω ) = exp λ . m(Ω ) , ( 19
)
with λ ∈ R 9 .
Using Theorem 1, there exists a bijection between the Lagrange multipliers λ ∈ R 9 and the moments (ψ 0 , ψ 1 , ψ 2 ) ∈ int(R 2 ) in the interior of the realizability domain (remark that this interior is int(R 2 ) = R 2 , see ( 16)). We consider the following hierarchy of subdomains of R 9

L 1 := (λ 1 , 0, 0, λ 4 , λ 5 , λ 6 , 0, 0, 0), s.t.

(λ 1 , λ 4 , λ 5 , λ 6 ) ∈ R 4 ⊂ R 9 , (20a)
L 2 := (λ 1 , 0, 0, λ 4 , λ 5 , λ 5 , 0, 0, 0), s.t. (λ 1 , λ 4 , λ 5 ) ∈ R 3 ⊂ L 1 , (20b)
L 3 := (λ 1 , 0, 0, λ 4 , λ 4 , λ 4 , 0, 0, 0), s.t. (λ 1 , λ 4 ) ∈ R 2 ⊂ L 2 . (20c)
Choosing λ in one of those sets in [START_REF] Duclous | A deterministic partial differential equation model for dose calculation in electron radiotherapy[END_REF] gives ansätze ψ R of the form

ψ 1 (Ω ) = exp(λ 1 Ω 1 + λ 4 Ω 2 1 + λ 5 Ω 2 2 + λ 6 Ω 2 3 ), (21a) ψ 2 (Ω ) = exp(λ 5 + λ 1 Ω 1 + (λ 4 -λ 5 )Ω 2 1 ), (21b)
ψ 3 (Ω ) = exp(λ 4 + λ 1 Ω 1 ) (21c)
respectively.

Definition 2

We denote by H i the set of normalized and rotated moments of functions of the form [START_REF] Friedriechs | Systems of conservation equations with a convex extension[END_REF]:

H i := (Ω , Ω ⊗ Ω ) exp( λ . m(Ω )) , λ ∈ L i ∩ R T . ( 22
)
Remark 4 Choosing the sets L 3 ⊂ L 2 ⊂ L 1 ⊂ R 9 was motivated by two reason. First the hierarchical character of these sets is necessary to the construction of the expansion of N 3 described in the next section. Second each scale of this hierarchy corresponds to a particular type of problems:

-H 1 corresponds to the case where ψ 1 is an eigenvector of ψ 2 (see computations in Appendix A.1). -H 2 corresponds to a 1D problem. Indeed, one can see that the distribution (21b) depends only on one scalar Ω 1 . -H 3 corresponds to a first-order model (i.e. M 1 model). Indeed, one can see that the distribution (21c) corresponds to a distribution function obtained using a first-order model (i.e. when working with (ψ 0 , ψ 1 )).

In order to illustrate this hierarchy, possible distribution functions associated to each of those sets are depicted in Fig. 1, 2, 3, 4, 5 and 6. On those plots, the color on the unit sphere corresponds to the value of the distribution function (where blue corresponds to the lowest value and red the highest) given in [START_REF] Friedriechs | Systems of conservation equations with a convex extension[END_REF]. As H 3 and H 2 are related to 1D distribution, those distributions can be represented along the preferred axis, i.e. ψ R in (21a) and (21b) as a function of Ω 1 . Computing the moments of order one and two of the functions (21) reads Fig. 4 Distribution function associated to a vector of

H 2 \H 3 along its preferred axis in H 1 , N 1 = N 1 1 e 1 , N 2 = |N 1 1 | 2 e 1 ⊗ e 1 + (1 -|N 1 1 | 2 )Diag(γ 1 , γ 2 , 1 -γ 1 -γ 2 ), ( 23a
) in H 2 , N 1 = N 1 1 e 1 , N 2 = |N 1 1 | 2 e 1 ⊗ e 1 + (1 -|N 1 1 | 2 )Diag(γ 1 , 1 -γ 1 2 , 1 -γ 1 2 ), ( 23b
) in H 3 , N 1 = N 1 1 e 1 , N 2 = 3χ 2 (|N 1 1 |) -1 2 e 1 ⊗ e 1 + 1 -χ 2 (|N 1 1 |) 2
Id, (23c) where χ is the Eddington factor (see e.g. [START_REF] Levermore | Relating Eddington factors to flux limiters[END_REF]). This leads to the following parametrization of H i

(N 1 , γ 1 , γ 2 ) ∈ H 1 := (N 1 , γ 1 , γ 2 ) ∈ R P s.t. N 1 = N 1 1 e 1 , (N 1 , γ 1 , γ 2 ) ∈ H 2 := (N 1 , γ 1 , γ 2 ) ∈ H 1 s.t. γ 2 = 1 -γ 1 2 , (N 1 , γ 1 , γ 2 ) ∈ H 3 := (N 1 , γ 1 , γ 2 ) ∈ H 2 s.t. γ 1 = χ 2 (|N 1 1 |) -|N 1 1 | 2 1 -|N 1 1 | 2 .
Similarly, computing the third-order moments of the distributions [START_REF] Friedriechs | Systems of conservation equations with a convex extension[END_REF] reads

in H 1 , N 3 = κ 2 1 1,1,1 + κ 3 T 1,2,2 + (N 1 1 -κ 2 -κ 3 )T 1,3,3 , (24a)
in

H 2 , N 3 = κ 1 1 1,1,1 + N 1 1 -κ 1 2 (T 1,2,2 + T 1,3,3 ), (24b)
in

H 3 , N 3 = χ 3 1 1,1,1 + N 1 1 -χ 3 2 (T 1,2,2 + T 1,3,3 ), (24c)
T i, j, j = 1 i, j, j + 1 j,i, j + 1 j, j,i , 1 i, j,k = e i ⊗ e j ⊗ e k , where χ 3 , κ 1 , κ 2 and κ 3 are scalar coefficients depending on (N 1 1 , γ 1 , γ 2 ) respectively in H 3 for χ 3 , H 2 for κ 1 , and H 1 for κ 2 and κ 3 .

Limits of N 3 on the boundary of R T

The boundary of R T is characterized by

N 2 -N 1 ⊗ N 1 has a zero eigenvalue. ( 25
)
At the boundary, the problem [START_REF] Caron | Deterministic model for the transport of energetic particles. application in the electron radiotherapy[END_REF] has no solution (see e.g. [START_REF] Hauck | Convex duality and entropy-based moment closures: Characterizing degenerate densities[END_REF][START_REF] Borwein | Duality relationships for entropy-like minimization problems[END_REF]). However, we can deduce the following two results, first on the boundary of H 1 and then in the particular case when ||N 1 || 2 = 1, which we will use in the next subsection. 1 and N 2 is given by (23a) and

Proposition 2 Suppose (N 1 , N 2 ) ∈ ∂ H 1 , i.e. N
∃V ∈ S 2 , s.t. V. N 2 -N 1 ⊗ N 1 .V = 0, ( 26
)
or equivalently γ 1 = 0, γ 2 = 0 or 1 -γ 1 -γ 2 = 0 in (23a). N 1 1 N 2 1,1 N 2 2,2 1
1 Fig. 7 Representation of H 3 (red line), H 2 (green plane) and H 1 (blue volume) in the space

(N 1 1 , N 2 1,1 , N 2 2,2 ) ∈ R 3 .
Then the closure yields

N 3 i, j, j = N 1 i N 2 j, j , N 3 1,2,3 = 0. ( 27
)
Proof The moments in

H 1 satisfy N 2 = Diag(N 2 1,1 , N 2 2,2 , N 2 3,3 ) and N 1 = N 1 1 e 1 . Therefore the eigenvectors of N 2 -N 1 ⊗ N 1 are the cartesian axis V = e i .
This result follows from [START_REF] Kershaw | Flux limiting nature's own way[END_REF], where the author showed that the possible representing distributions ψ R (Ω ) for moments satisfying [START_REF] Hauck | Positive P N closures[END_REF] are zero everywhere except on the line Ω ∈ S 2 , s.t. (Ω -N 1 ).V = 0 . Computing the following moments provides part of the equalities ( 27)

S 2 Ω ⊗ Ω [(Ω -N 1 ).V ]ψ R (Ω )dΩ = N 3 .V -(N 1 .V )N 2 = 0 R 3×3 .
The other equalities are simply obtained by remarking that some of the values of M 2 closure N 3 are zero when (N 1 , N 2 ) ∈ H 1 , e.g. using the eveness of (21a) according to Ω 2 , one finds that N 3 2,2,2 = N 3 2,3,3 = N 3 1,2,3 = 0. As the M 2 closure is a continuous function of N 1 and N 2 , one retreives those zeros of N 3 on the boundary ∂ H 1 . Fig. 7 depicts the hierarchy H 3 ⊂ H 2 ⊂ H 1 and its boundary in the space (N 1 1 , N 2 1,1 , N 2 2,2 ) ∈ R 3 (these three components are indeed sufficient to parametrize H 1 ). H 3 is the red line included in H 2 the green plane, itself included in H 1 the blue volume.

Proposition 3 ([31]) Suppose

(N 1 , N 2 ) ∈ ∂ R T , such that ||N 1 || 2 = 1, then N 2 = N 1 ⊗ N 1 ,
and the closure yields

N 3 i, j,k = N 1 i N 1 j N 1 k . ( 28
)

Approximation of the M 2 closure: the expansion

We first start by approximating N 3 in the set H 3 . Then we extend it progressively into H 2 , H 1 , and finally R 2 .

The approximation is based on polynomial interpolations. For convenience, we introduce the following notation Notation 4 The polynomial of degree two interpolating the values A, B and C at the points a, b and c is denoted E, and Z denotes the polynomial of degree three which is zero in a, b and c:

E ((A, a), (B, b), (C, c)) (x) := A x -b a -b x -c a -c + B x -a b -a x -c b -c +C x -a c -a x -b c -b , Z (a, b, c) (x) := (x -a)(x -b)(x -c).
The hyperbolicity and the realizability properties (see Subsection 2.2) are considered at each step of the construction of approximated closure. The realizability, the hyperbolicity and the precision of the approximated closure are finally studied in the next subsection.

Initialization: the closure in H 3

First, one needs to approximate the Eddington factor χ 2 to compute the M 1 closure [START_REF] Curto | Recusiveness, positivity, and truncated moment problems[END_REF]. As in the characterization of ( 16), the moments (ψ 0 , ψ 1 ) are the moments of a positive function if and only if ( [START_REF] Kershaw | Flux limiting nature's own way[END_REF])

||ψ 1 || 2 ψ 0 < 1.
Similarly using [START_REF] Curto | A duality prood to Tchakaloff's theorem[END_REF], simple computations show that (ψ 0 , ψ 1 , ψ 2 ), where ψ 2 has the form (15), are the moments of a positive function if and only if

x 2 < χ 2 (x) ≤ 1, for x ∈ [0, 1[.
In order to construct a realizable closure (see Remark 2), we propose to approximate χ 2 by the convex combination

χ 2 (x) ≈ x 2 θ 1 (x) + (1 -θ 1 (x))1.
The coefficient θ 1 is chosen to be an even function of x so that θ 1 (x) ∈ [0, 1] and the following exact value of the Eddington factor χ 2 and its derivative are satisfied on the boundary of the realizability domain (i.e. at x = 1) and for the isotropic case (i.e. in x = 0). The choice of fixing the exact values of the derivative χ ′ 2 is motivated by the hyperbolic character of the M 2 closure (see the "hyperbolicity" paragraphs of Subsections 2.2 and 3.4). Those values are obtained by rewriting

χ 2 (|N 1 |) = Ω 2 1 exp(α(|N 1 |)Ω 1 ) exp(α(|N 1 |)Ω 1 ) , χ ′ 2 (|N 1 |) = dα d|N 1 | d dα Ω 2 1 exp(α(|N 1 |)Ω 1 ) exp(α(|N 1 |)Ω 1 )
,

where dα d|N 1 | = d|N 1 | dα -1
. Then one can remark that

|N 1 |(α = 0) = 0, lim α→+∞ |N 1 | = 1,
and since |N 1 | is a bijection of α (Theorem 1), this implies that

α(|N 1 | = 0) = 0, lim |N 1 |→1 α = +∞,
then computing χ 2 and its derivative at those values using a symbolic computation software such as Maple TM ( [START_REF] Maple | Technical report[END_REF]) leads to

χ 2 (1) = 1, χ 2 (0) = 1 3 , χ ′ 2 (1) = 2, χ ′ 2 (0) = 0. ( 29
)
In the end, we chose

θ 1 (x) = x 2 + 2 3 (1 -x 2 ) + x 2 (1 -x 2 ) c 0 + c 1 x 2 + c 2 x 4 ,
where the coefficients c i are fitted to approximate the exact χ 2 for 10 3 values of x equally dis- tributed in [0, 1] using Maple TM ( [START_REF] Maple | Technical report[END_REF]): c 0 = -0.0954823981432433, c 1 = 0.229069986304953 and c 2 = -0.0344846229504588. One can check that θ 1 (x) ∈ [0, 1] for all x ∈ [-

-(N 1 1 , N 2 1,1 ) < N 3 1,1,1 < b + (N 1 1 , N 2 1,1 ), b -(N 1 1 , N 2 1,1 ) := -N 2 1,1 + (N 1 1 +N 2 1,1 ) 2 (1+N 1 1 ) , b + (N 1 1 , N 2 1,1 ) := N 2 1,1 - (N 1 1 -N 2 1,1 ) 2 (1-N 1 1 ) .
Then similarly to above, we construct a realizable closure in H 3 as a convex combination (see Remark 2)

χ 3 (x) = b -(x, χ 2 (x))θ 2 (x) + b + (x, χ 2 (x))(1 -θ 2 (x)).
The coefficient θ 2 is chosen so that θ 2 (x) ∈ [0, 1], such that χ 3 is odd and the following exact value of χ 3 and its derivative are satisfied, on the boundary of the realizability domain (i.e. at x = 1) and for the isotropic case (i.e. at x = 0)

χ 3 (1) = 1, χ 3 (0) = 0, χ ′ 3 (1) = 3, χ ′ 3 (0) = 1 2 .
Those values are obtained from the same method as for χ 2 in (29). In the end, we chose

θ 2 (x) = 1 2 + x - 1 2 + (1 -x 2 )(d 0 + d 1 x 2 + d 2 x 4 ) ,
where the coefficients d i are fitted to approximate the exact χ 3 for 10 3 values of x equally dis- tributed in [0, 1] using Maple TM ( [START_REF] Maple | Technical report[END_REF]): d 0 = 0.386143553495150, d 1 = 0.488034553677475 and d 2 = -0.681343955348390. One can check that θ 2 (x) ∈ [0, 1] for all x ∈ [-1, 1].

Extension to H 2 : the 1D approximated closure

Using the above results and the form of N 3 in (24b), we define an approximation of N 3 in H 2 by approximating κ 1 as follows

κ 1 (x, y) = b -(x, y)θ 3 (x, y) + b + (x, y)(1 -θ 3 (x, y)).
The coefficient θ 3 is chosen such that θ 3 (x, y) ∈ [0, 1] so that the resulting closure is realiz- able (see Remark 2). Furthermore it is chosen to have the exact value on the boundary of the realizability domain (i.e. when (N 1 , N 2 ) ∈ ∂ H 2 as in Proposition 2) and the approximated value in H 3 given above, that is

κ 1 (N 1 1 , 0) = (N 1 1 ) 3 , κ 1 N 1 1 , χ 2 (N 1 1 ) -|N 1 1 | 2 1 -|N 1 1 | 2 = χ 3 (N 1 1 ), κ 1 (N 1 1 , 1) = N 1 1 .
In the end, we chose

T 1 := b + -κ 1 b + -b -(x, 0), T 2 := b + -κ 1 b + -b -x, χ 2 (x)-x 2 1-x 2 , T 3 := b + -κ 1 b + -b -(x, 1), θ 3 (x, y) = E (T 1 , 0) , T 2 , χ 2 (x)-x 2 1-x 2 , (T 3 , 1) (y) + Z 0, χ 2 (x)-x 2 1-x 2 , 1 (y)Q 1 (x, y)
, where Q 1 is a polynomial of x and y of degree sixteen. Its coefficients are chosen such that the discrete L 2 distance between the approximated and the exact κ 1 (computed by solving [START_REF] Caron | Deterministic model for the transport of energetic particles. application in the electron radiotherapy[END_REF] for 10 4 values of (N 1 , N 2 ) ∈ H 2 , given by 100 values of x equally distributed in [0, 1] and 100 of y equally distributed in [0, 1]) is minimized. The degree of this polynomial approximation is chosen very high because a high precision is required for the next extension (i.e. to H 1 and then to R T ). The discrete L ∞ error compared to the solution of the minimization problem [START_REF] Caron | Deterministic model for the transport of energetic particles. application in the electron radiotherapy[END_REF] for those 10 4 values of (N 1 , N 2 ) ∈ H 2 is 8.43 × 10 -3 . Our polynomial satisfies θ 3 (x, y) ∈ [0, 1] for all x ∈ [-1, 1] and y ∈ [0, 1].

Extension to H 1

This procedure can be repeated to approximate N 3 in H 1 .

Using Proposition 2 and the previous approximation, we aim to write an approximation of κ 2 and κ 3 (in (24a)) satisfying

κ 2 (N 1 1 , γ 1 , 0) = (N 1 1 ) 3 + N 1 1 (1 -|N 1 1 | 2 )γ 1 , κ 2 (N 1 1 , γ 1 , 1 -γ 1 2 ) = κ 1 (N 1 1 , γ 1 ), κ 2 (N 1 1 , γ 1 , 1 -γ 1 ) = (N 1 1 ) 3 + N 1 1 (1 -|N 1 1 | 2 )γ 1 , κ 3 (N 1 1 , γ 1 , 0) = 0, κ 3 (N 1 1 , γ 1 , 1 -γ 1 2 ) = 1 2 N 1 1 -κ 1 (N 1 1 , γ 1 ) , κ 3 (N 1 1 , γ 1 , 1 -γ 1 ) = (N 1 1 ) 3 + N 1 1 (1 -|N 1 1 | 2 )(1 -γ 1
). Similarly to above, we propose an approximation of the form

κ 2 (x, y, z) = E (κ 2 (x, y, 0), 0) , κ 2 (x, y, 1 -y 2 ), 1 -y 2 , (κ 2 (x, y, 1 -y), 1 -y) , (z) +Z 0, 1 -y 2 , 1 -y (z)Q 2 (x, y, z), κ 3 (x, y, z) = E (κ 3 (x, y, 0), 0) , κ 3 (x, y, 1 -y 2 ), 1 -y 2 , (κ 3 (x, y, 1 -y), 1 -y) , (z) +Z 0, 1 -y 2 , 1 -y (z)Q 3 (x, y, z),
where Q 2 and Q 3 are polynomials in x, y and z of degree eight. The coefficients of those polynomials are chosen such that the discrete L 2 distance between the approximated and the exact κ 2 and κ 3 (computed by solving (12) for 8 × 10 3 values of (N 1 , N 2 ) ∈ H 1 , given by 20 values of x equally distributed in [0, 1], and 20 values of y in [0, 1] and 20 of z in [0, 1y]) is minimized. The discrete L ∞ error compared to the solution of the minimization problem [START_REF] Caron | Deterministic model for the transport of energetic particles. application in the electron radiotherapy[END_REF] for those 8 × 10 3 values of (N 1 , N 2 ) ∈ H 1 is of 2.09.10 -2 .

Extension to the whole realizability domain R T

We now aim to extend this approximation to R T which will provide us with an approximation of ψ 3 for any (ψ 0 , ψ 1 , ψ 2 ) ∈ R 2 through (17). This last extension again consists of an interpolation.

The previous approximation provides a closure when N 1 is along one of the Cartesian axes. One can also compute the closure when |N 1 | = 1 through Proposition 3. Now suppose we want to compute N 3 at the point

P 0 = (N 1 , N 2 ) ∈ R T , parametrized by (N 1 , γ 1 , γ 2 ) = (xe 1 + ye 2 + ze 3 , β 1 , β 2 ).
Let us define the following points (see Fig. 8 and9)

P 1 ≡ (N 1 , γ 1 , γ 2 ) = (xe 1 , β 1 , β 2 ), P 2 ≡ (N 1 , γ 1 , γ 2 ) = (ye 2 , β 1 , β 2 ), P 3 ≡ (N 1 , γ 1 , γ 2 ) = (ze 3 , β 1 , β 2 ),
corresponding to projections of P 0 onto each Cartesian axis. At those points, N 1 is an eigenvalue of N 2 . Then we can use the approximation techniques in H 1 of the previous paragraph (see remark 4). Now let us define the lines and points (see Fig. 8 and9) L 1 = (P 1 , P 0 ), L 2 = (P 2 , P 0 ), L 3 = (P 3 , P 0 )

P 4 = L 1 ∩ {||N 1 || 2 = 1} ≡ (xe 1 + a 1 ye 2 + a 1 ze 3 , β 1 , β 2 ), a 1 = 1-x 2 y 2 +z 2 , P 5 = L 2 ∩ {||N 1 || 2 = 1} ≡ (a 2 xe 1 + ye 2 + a 2 ze 3 , β 1 , β 2 ), a 2 = 1-y 2 x 2 +z 2 , P 6 = L 3 ∩ {||N 1 || 2 = 1} ≡ (a 3 xe 1 + a 3 ye 2 + ze 3 , β 1 , β 2 ), a 3 = 1-z 2
x 2 +y 2 .

The closure N 3 is either known (at P 4 , P 5 and P 6 ) or approximated (at P 1 , P 2 and P 3 ) at each of those points. Thus for (N 1 , N 2 ) ∈ R T the different components of the closure N 3 are approximated by convex combinations of the values of N 3 computed or approximated at the points P i for i ∈ {1, 2, . . . 6}. In particular, for i ∈ {1, 2, 3}, we approximate N 3 i, j, j with a convex combination of its approximate value for P i and its exact value for P 3+i . Similarly, the value of N 3 1,2,3 is known at the points

(N 1 , γ 1 , γ 2 ) = (0 R 3 , γ 1 , γ 2 ) and (N 1 , γ 1 , γ 2 ) = (N 1 /||N 1 || 2 , γ 1 , γ 2
), and we simply approximate it by a convex combination of the value at those points. In the

Properties of the approximation

We study here the accuracy and the numerical cost to compute the approximated closure, the hyperbolicity and the realizability property with the approximated closure.

Precision and numerical cost

Our approach consists of a polynomial approximation. Evaluating the polynomials is very fast compared to solving the minimization problem [START_REF] Caron | Deterministic model for the transport of energetic particles. application in the electron radiotherapy[END_REF], as the numerical methods generally require iterative algorithms. The approximate closure is compared to the one obtained by solving this minimization problem using the routines HUMSL of MINPACK [START_REF] Moré | User Guide for MINPACK-1[END_REF] and DCUHRE from [START_REF] Berntsen | Algorithm 698: Dcuhre: An adaptive multidemensional integration routine for a vector of integrals[END_REF]. The routine HUMSL is based on an iterative algorithm using the gradient and the Hessian of ψ 3 . The routine DCUHRE is an adaptive quadrature algorithm for functions of several variables. The precision for this minimization algorithm (both for the minimization and the quadrature) can be fixed. In order to have a fair comparison of computational times and precision, we first ran the code with high precision (L ∞ error of 10 -9 ) in order to have an accurate reference (it ran for 3h 33min), and then we ran it again with a precision equivalent to that of our approximation. We fixed the maximum number of iterations for both the quadrature and the minimization routine at 10 6 iterations, and we checked that this was not reached during the computations (meaning the desired precision was always obtained). Computing the approximation is much faster than solving the minimization problem with this method as shown in Table 1.

Minimization solver Approximation Computation times

1654 sec = 27min 34sec 0.434 sec

L ∞ error < 3 × 10 -2 3.12 × 10 -2
Table 1 Comparing computation times for the closure with 3.2 × 10 6 points in R T between the minimization solver and our approximation.

Hyperbolicity

It is well-known [START_REF] Levermore | Moment closure hierarchies for kinetic theories[END_REF] that the M 2 system is a hyperbolic system of conservation laws, as long as (ψ 0 , ψ 1 , ψ 2 ) ∈ R 2 . This means that the Jacobian ∇ (ψ 0 ,ψ 1 ,ψ 2 ) (ψ 1 , ψ 2 , ψ 3 ) is diago- nalizable with real eigenvalues. In this paragraph, we study the hyperbolicity of our approximation of the M 2 closure at each step of the expansion of the previous subsection, i.e. when (N 1 , N 2 ) is in H 3 , H 2 , H 1 , and R T .

For the M 1 approximation, i.e. for moments in H 3 , we can work in a reference frame in which ψ 1 = ||ψ 1 ||e 1 . Using this simple rotation, one can show that the Jacobian reduces into

∇ (ψ 0 ,||ψ 1 || 2 ) ||ψ 1 || 2 ψ 0 χ 2 ( ||ψ 1 || 2 ψ 0 ) = 0 1 χ 2 ( ||ψ 1 || 2 ψ 0 ) -||ψ 1 || ψ 0 χ ′ 2 ( ||ψ 1 || 2 ψ 0 ) χ ′ 2 ( ||ψ 1 || ψ 0 )
. This leads to the following requirement for hyperbolicity

∀x ∈ [0, 1], χ ′ 2 (x) -4(χ 2 (x) -xχ ′ 2 (x)) ≥ 0,
which is satisfied by our approximation. In H 2 (i.e. the 1D case), using the same method, one obtains the hyperbolicity if the following matrix is diagonalizable with real eigenvalues

  0 1 0 0 0 1 N 3 1,1,1 -N 1 1 ∂ N 1 1 N 3 1,1,1 -N 2 1,1 ∂ N 2 1,1 N 3 1,1,1 ∂ N 1 1 N 3 1,1,1 ∂ N 2 1,1 N 3 1,1,1   .
Studying the roots of the characteristic polynomial of this matrix (i.e. roots of a cubic function), one can verify that this matrix is indeed diagonalizable with real eigenvalues. This means that the M 2 system of equations for 1D problems with the approximated closure is hyperbolic.

In H 1 and in R T , the problem is 3D. The Jacobian can not be easily reduced into a smaller matrix, so the full Jacobian is studied. As studying moments in H 1 does not provide any simplification of the Jacobian compared to moments in R T , we directly study the hyperbolicity in the general case, for moments in R T . Checking that the eigenvalues of the full Jacobian ∇ (ψ 0 ,ψ 1 ,ψ 2 ) (ψ 1 , ψ 2 , ψ 3 ) are real for all (N 1 , N 2 ) ∈ R T is complicated, and we have not been to verify analytically that the eigenvalues are always real. Instead, we verified this property was satisfied for a finite number of points in R T , i.e. for 10 5 values of (N 1 , N 2 ) ∈ R T . Those values are obtained from 10 values of

N 1 1 equally distributed in [0, 1], 10 of N 1 2 in [0, 1 -|N 1 1 | 2 ], 10 of N 1 2 in [0, 1 -|N 1 1 | 2 -|N 1 2 | 2 ], 10 of γ 1 in [||N 1 || 2 2 , 1] and 10 of γ 2 in [0, 1].
Remark 5 Since solving the minimization problem [START_REF] Caron | Deterministic model for the transport of energetic particles. application in the electron radiotherapy[END_REF] numerically also introduces errors, this may also result in a loss of hyperbolicity.

Realizability

For any set of realizable moments (ψ 0 , ψ 1 , ψ 2 ) ∈ R 2 , we have constructed a closure ψ 3 . One may expect that the set of moments (ψ 0 , ψ 1 , ψ 2 , ψ 3 ) ∈ R 3 is also realizable, i.e. are the moments of one positive function. The realizability condition on ψ 3 in 1D (see e.g. [START_REF] Kershaw | Flux limiting nature's own way[END_REF][START_REF] Curto | Recusiveness, positivity, and truncated moment problems[END_REF]), i.e. in the particular case when (ψ 0 , ψ 1 , ψ 2 ) ∈ H 2 , can easily be verified. In our framework, this condition can be rewritten (see Remark 2)

θ 3 ∈ [0, 1],
which is satisfied by our approximation. This means that the proposed closure is realizable when working with 1D problems. In multi-D, there is, to the authors' knowledge, no similar characterization of the realizability to check. As our approach approximates a realizable closure ψ 3 and has the right value in limit case when (ψ 0 , ψ 1 , ψ 2 ) ∈ ∂ R T , we may expect the approximation to be realizable.

Remark 6

Numerically solving the minimization problem [START_REF] Caron | Deterministic model for the transport of energetic particles. application in the electron radiotherapy[END_REF] introduces errors which may also result in a loss of realizability.

Numerical results

Now we compare the solutions computed with the M 1 and M 2 models (both with approximated closure and by using the minimization algorithm previously mentioned for [START_REF] Caron | Deterministic model for the transport of energetic particles. application in the electron radiotherapy[END_REF]) to those results obtained by solving the kinetic model directly on several test cases. First, the numerical approaches used to solve the kinetic and moment equations are described. Then the results for each test cases are presented. Those tests consist in computing the dose for given boundary conditions which correspond to: a single beam in 1D, two opposite beams in 1D, a single beam in 2D whose direction is not aligned with the mesh, two beams in 2D and finally a single beam in an inhomogeneous medium corresponding to a human chest.

Numerical approaches

In each test case, the moments results are compared to the results of a reference code.

In 1D, this reference code is a kinetic solver for (1). However, discretizing (1) directly with the LB operator (2) or CSD operator (3) is not efficient. Indeed the differential cross sections in ( 2) and ( 3) are very peaked, and such a discretization would require a meshgrid in angle and energy fine enough to represent those peaks. We could not afford to use such a fine mesh. For the same reason, the proposed moment solver is only used for the equations ( 5) with the CSD [START_REF] Berthon | Numerical methods for balance laws with space dependent flux: application to radiotherapy dose calculation[END_REF] or FP [START_REF] Borwein | Duality relationships for entropy-like minimization problems[END_REF] operator. In 1D, the kinetic and moment solvers are based on the equation ( 1) and ( 5) with a FP operator ( 4) and [START_REF] Borwein | Duality relationships for entropy-like minimization problems[END_REF].

In 1D, we also compared the results with the approximated closures to the one using the closure obtained by solving the minimization problem [START_REF] Caron | Deterministic model for the transport of energetic particles. application in the electron radiotherapy[END_REF]. This closure was obtained by using the routines "HUMSL" of MINPACK [START_REF] Moré | User Guide for MINPACK-1[END_REF] and "DQAGS" of QUADPACK [START_REF] Piessens | QUADPACK: a subroutine package for automatic integration[END_REF], using a maximum tolerance of 3.10 -2 (equivalent to the precision as our approximation) and a maximum number of iterations of 100. For the purposes of dose calculation these routines provide an accurate M N closure which is easy to incorporate into a numerical solver. However, for cases where the moments are closer to the boundary of realizability (and the optimization problem becomes numerically ill-conditioned) special techniques have been developed [START_REF] Hauck | High-order entropy-based closures for linear transport in slab geometry[END_REF][START_REF] Alldredge | High-order entropy-based closures for linear transport in slab geometry II: A comutational study of the optimization problem[END_REF][START_REF] Alldredge | Adaptive change of basis in entropy-based moment closures for linear kinetic equations[END_REF] which make the entropy closure problem tractable.

We used the Monte-Carlo solver PENELOPE ( [START_REF] Sempau | PENELOPE-2011: A Code System for Monte Carlo Simulation of Electron and Photon Transport[END_REF]) as reference. PENELOPE is one of the state-of-the-art codes for electron transport, and has been validated against experiments. Note that PENELOPE, when compared to our model, also takes into account more physical effects (like e.g. pair production and bremsstrahlung).

The numerical schemes are constructed as follows: Due to the integral in (4), in order to compute ψ at energy ε, one requires the knowledge of the value of ψ for all energies between ε and ε max . In practice, we solve it from a maximum energy ε max to zero.

The 1D kinetic equation is solved using an upwind scheme for the spatial and energy derivatives. The angular diffusion operator is discretized using a central scheme, and the integral operator with a midpoint quadrature rule.

The moment solver was obtained by adapting the relaxation approach proposed in [START_REF] Pichard | Relaxation schemes for the M 1 model with space-dependent flux: application to radiotherapy dose calculation[END_REF] to the equations ( 5) with the collisions operators ( 8) and [START_REF] Berthon | Numerical methods for balance laws with space dependent flux: application to radiotherapy dose calculation[END_REF]. For the sake of clarity, the numerical scheme for moment equations is written here for a 1D problem with a Fokker-Planck operator [START_REF] Borwein | Duality relationships for entropy-like minimization problems[END_REF]. It can be easily extended to 2D or 3D problems and to the CSD colllision operator [START_REF] Berthon | Numerical methods for balance laws with space dependent flux: application to radiotherapy dose calculation[END_REF]. First let us rewrite the moment system 5 with the FP operator (8) under the form (30) where ψ is a vector of moments and F( ψ) its associated flux; σ is the matrix composed of the moments σ i of the cross section; and Ā is a matrix such that Ā ψ corresponds to the moments of the angular diffusion term, e.g. for

∂ x F( ψ)(x, ε) = ρ(x) ∂ ε (S ψ)(x, ε) -T (ε) Ā ψ(x, ε ′ ) + ε max ε σ (ε ′ , ε) ψ(x, ε ′ )dε ′ ,
M 1 ψ = (ψ 0 , ψ 1 1 ) T , F( ψ) = (ψ 1 1 , ψ 2 1,1 ) T , σ = Diag(σ 0 , σ 1 ), Ā = Diag(0, 2).
In the following, the subscript l corresponds to the spatial variable x l and the superscript p to energy variable ε p . In order not to be constrained by a Courant-Friedrichs-Lewy (CFL) condition which could be very restrictive for our problem (see e.g. [START_REF] Berthon | Numerical methods for balance laws with space dependent flux: application to radiotherapy dose calculation[END_REF][START_REF] Pichard | Relaxation schemes for the M 1 model with space-dependent flux: application to radiotherapy dose calculation[END_REF]), we aim to construct an implicit numerical scheme. We propose the following scheme, which was obtained by discretizing in [START_REF] Junk | Maximum entropy for reduced moment problems[END_REF], the integral with a quadrature formula, the energy derivative with an Euler discretization and the spatial flux with an implicit form of HLL ( [START_REF] Harten | On upstream differencing and Gudonov-type schemes for hyperbolic conservation laws[END_REF]) scheme

F( ψ) p l+1 -F( ψ) p l-1 2∆ x - ψ p l-1 -2 ψ p l + ψ p l+1 2∆ x +ρ l - (S ψ) p-1 l -(S ψ) p l ∆ ε p + T p Ā -σ p,p ∆ ε p Id ψ p l - p-1 ∑ p ′ =1 σ p ′ ,p ψ p ′ l ∆ ε p ′ = 0.
This is an implicit equation for ψ p . We use an iterative solver inspired by [START_REF] Dubroca | An iterative method for transport equations in radiotherapy[END_REF], which consists of sweeping in the l variable. Let us write ψ p,k l , the intermediate solution after sweeping k times. The solution is initialized by choosing ψ p,0 l = ψ p-1 l , then the following systems are solved iteratively: First

F( ψ) p,k l+1 -F( ψ) p,k+1 l-1 2∆ x - ψ p,k+1 l-1 -2 ψ p,k+1 l + ψ p,k l+1 2∆ x ( 31
)
+ρ l - (S ψ) p-1 l -(S ψ) p,k+1 l ∆ ε p + T p Ā -σ p,p ∆ ε p Id ψ p,k+1 l - p-1 ∑ p ′ =1 σ p ′ ,p ψ p ′ l ∆ ε p ′ = 0.
This corresponds to letting l run from 0 to l max , i.e. from left to right, and implicitly using what has already been computed (i.e. ψ p,k+1 l-1 ) and explicitly using what is not (i.e. ψ p,k l+1 ). Then we have l run from l max to zero, i.e. from right to left

F( ψ) p,k+1 l+1 -F( ψ) p,k l-1 2∆ x - ψ p,k l-1 -2 ψ p,k+1 l + ψ p,k+1 l+1 2∆ x ( 32
)
+ρ l - (S ψ) p-1 l -(S ψ) p,k+1 l ∆ ε p + T p Ā -σ p,p ∆ ε p Id ψ p,k+1 l - p-1 ∑ p ′ =1 σ p ′ ,p ψ p ′ l ∆ ε p ′ = 0.
Equations [START_REF] Kershaw | Flux limiting nature's own way[END_REF] and [START_REF] Larsen | Electron dose calculations using the method of moments[END_REF] are solved iteratively until reaching a converged state, i.e. until k satisfies

ψ p,k+1 l -ψ p,k l ψ p,k l < r
for all l, and for a chosen maximum residual r. In practice, we chose r = 10 -1 , and a minimum number of iterations k of 3. The function of interest for medical physicists is the dose given by

D(x) = ε max 0 S(ε)ψ 0 (x, ε)dε.
For each test case we compute the dose produced by beams of electrons prescribed on the boundary of the medium by for Ω .n < 0, ψ(x, ε, Ω ) = 10 10 exp -c e (ε 0 -ε)

2 exp -c o (1 -Ω p .Ω ) 2 , ( 33
)
where ε 0 , Ω p are the energy and direction, respectively, of the prescribed beam, and n is a vector normal to the boundary in the outgoing direction.

If no beam is applied on one boundary we apply a zero distribution and assume that all particles injected at the other borders are stopped inside the medium.

The boundary conditions for the moment systems are defined by extracting moments from these boundary conditions.

The following test cases provide qualitative analysis of the approximation of the M 2 closure. Further analysis of the numerical schemes for moments and quantitative study (convergence results) will be provided in another paper. In all test cases, we used mesh cells ∆ x = 0.01 cm and 100 cells in energy, with a tolerance r = 10 -1 in the scheme for moments models. This precision was sufficient to observe the considered phenomena.

Simple beam in 1D

In this test case a 6 cm long uniform water phantom is irradiated with a beam of electrons. This is modelled with fixed ρ = 1 and boundary conditions [START_REF] Levermore | Relating Eddington factors to flux limiters[END_REF] at x = 0 cm : ε 0 = 10 MeV, Ω p = e 1 , c e = 200 and c o = 10 3 .

The domain was meshed with 600 cells in position and 100 cells in energy. For this test case, we normalize the dose by the quantity of injected particles, i.e.

ψ i = Ω 1 >0 ε max 0 ψ(0, ε, Ω )dεdΩ .
The doses obtained with the kinetic solver, the M 1 and M 2 solver with the closures obtained from the minimization procedure and approximations are represented on Fig. the M N models have the same qualitative behaviour as the kinetic reference. The dose obtained with the M 1 model is however imprecise. Its derivative is too high at the entry and raised slowly and gradually until reaching its maximum, while the dose obtained with the M 2 and the kinetic models have a lower derivative at the entry, a higher maximum, and faster decrease after the maximum. These differences are due to the error produced when approximating the kinetic solution of (1) by a M N ansatz, as described in Section 2.4, but the results of the M 2 model already follow the kinetic reference with good precision. The approximations of the M 1 and the M 2 closures give good agreement with the closure obtained from the minimization procedure. The accuracy of the approximated closure is characterized by the errors on the doses gathered in Table 3. As described in Table 2, the approximations of the M N closures significantly accelerate (by around 200 times) the computation compared to using a minimization solver, which are themselves faster (between 5 to 10 times faster) than the kinetic computations. Table 3 Discrete L ∞ and L 2 error between the doses obtained with the approximated and the minimization closure respectively for the M 1 and M 2 models.

Double beam in 1D

The multi-beam instability in 1D is studied through this test case (see Section 2.4): Two beams of same energy and intensity but with opposite direction are prescribed on both ends of a 8 cm long homogeneous water phantom, characterized by [START_REF] Levermore | Relating Eddington factors to flux limiters[END_REF] The domain was meshed with 800 cells in position and 100 cells in energy. The dose normalized by ψ i obtained with the Monte-Carlo solver, and the M 1 , M 2 , double-M 1 (i.e. computing the dose of each beam seperately with M 1 ) and double-M 2 models is given on Fig. 11, and the computational times are presented in Table 4. The dose with the double M N model was simply obtained from the dose of the previous test case by summing D(x) + D(8 cmx) (see Remark 3). As described in [START_REF] Hauck | High-order entropy-based closures for linear transport in slab geometry[END_REF], when using a M N model for this problem, a shock appears in the middle of the medium. For M 1 this gives a high dose peak in the center of the medium, and for M 2 by a significantly smaller drop. It has been observed in [START_REF] Hauck | High-order entropy-based closures for linear transport in slab geometry[END_REF] that the higher is the order of model, the smaller this shock becomes. This artificial shock does not appear when considering the two beams separately. Through this case, we see that the non-physical effects of the angular approximation are smaller as closer to the ones of the Monte-Carlo reference, one can remark that it is also somewhat overly diffusive.

The absolute errors between the doses obtained with the PENELOPE Monte-Carlo solver and the approximated M 1 and M 2 models are represented on Fig. 13.

The maximum error in the dose when using the approximated M 1 or M 2 model compared to the reference Monte-Carlo results are located at the entry of the medium and on both sides of the beam, and, when using the approximated M 1 model, in the middle of the beam between 1 and 2 cm depth.

Double beam in 2D

A 2D version of the test case 4.3 is now studied. It consists of two beams crossing each other with at an angle of π/2. As described in Section 2.4, M 1 model is known to fail to represent this phenomenon (also in 2D). The same mesh as the previous case was used (i.e. 100 cells in energy and 600×600 cells in position), so the computation times for this case are identical to the previous one (see Table 6). The dose obtained with the different models are represented on Fig. 14. As expected, when using the M 1 model, the two incoming beams turn into one of direction (e 1 + e 2 ). This effect is artificial and is due to M 1 moments extraction. It does not appear when using the M 2 model. The dose obtained with the M 2 model is slightly more diffusive than the dose obtained with the Monte-Carlo solver.

The absolute errors between the doses obtained with the PENELOPE Monte-Carlo solver and the approximated M 1 , M 2 , double-M 1 and double-M 2 models are represented on Fig. 15.

With the M 2 model, no visible shock appears, in contrast to the 1D double beam case of Section 4.3. However, the M 2 model gives a dose slightly more diffused than that of the double-M 2 model. The dose as computed by the double-M 2 model is indeed very close to the Monte-Carlo reference.

As in 1D, the double-beam instability does not appear when using the double-M N models.

Chest geometry

This last test case is a 2D simple beam case with a density ρ corresponding to the density of a human chest. Here we show that our approach is valid with more realistic geometries. The isodose curves at 5%, 10%, 25%, 50%, 70% and 80% of the maximum dose obtained with the different models are plotted on Fig. 16, and the computational times are given in Table 7. As in the previous test cases, the dose computed with the M 1 solver is more diffusive 7 Characteristics of the computations with the different models than the others. In particular, the dose in the lungs (the dark regions) is overestimated. With both moment models, the dose is underestimated at the entry in the region corresponding to the backbone.

To characterize the errors when considering heterogeneous media, the 3% -3 mm distanceto-agreement is often used. A voxel is within this agreement if the absolute error in the dose compared to the reference at this point is lower than 3% of the maximum dose or if the dose obtained with the reference code at this point is also obtained with the moment code in a radius of less than 3mm around this voxel. The percentage of voxels not satisfying the 3% -3 mm distance-to-agreement with the approximated M 1 and M 2 model compared to the reference Monte-Carlo solver for this test case are gathered in Table 8. For both models these voxels are located at the entry on both sides of the beam. For the M 1 model such voxels are also found in the middle of the medium between 4 and 5 cm depth (in the backbone).

Model approximated M 1 approximated M 2 Percentage of voxels 0.12% 0.054% Table 8 Percentage of voxels not satisfying the 3% -3 mm distance-to-agreement compared to the reference Monte-Carlo dose for the 2D simple beam test case in a chest.

See also [START_REF] Pichard | The M 2 model for dose simulation in radiation therapy[END_REF] for more applications of those models in photon transport for radiotherapy.

Conclusion

We have proposed an approximation of the M 2 closure and have shown this model to be significantly better than the M 1 model. This approximation is based on the construction of entropy-based closures and the hierarchical structure of such models. In particular, the approximate closure is consistent with the exact closure for the moments associated with isotropic, M 1 , 1D, or 2D distributions. Numerical tests show that the M 2 model is much more accurate than the M 1 model and is valid for a larger range of physical phenomena. The dose computed from the M 2 model is close to the one provided by a reference Monte-Carlo code and required a much lower computational time (between fifty and two hundred times faster). Future work includes a comparison against optimization methods [27,2,1] which were specifically designed for entropy closures. Furthermore, we plan to generalize our physical model to the coupled system of photons, electrons and positrons, and include bremsstrahlung as well as pair production. We expect that these refinements make the dose calculation accurate enough for clinical purposes.

A Computation of moments in H i

A.1 Computation of moments in H 1

Consider λ ∈ L 1 , and the associated ansatz ψ 1 (defined in (21a)) and its moments (ψ 1 , ψ 2 ) given by

ψ 1 (Ω ) = exp λ 1 Ω 1 + λ 4 Ω 2 1 + λ 5 Ω 2 2 + λ 6 Ω 2 3 , (ψ 1 , ψ 2 ) = (Ω , Ω ⊗ Ω )ψ 1 ∈ H 1 .
One can remark that ψ 1 is an even function of Ω 2 or Ω 3 , and therefore the moment of ψ 1 according to any odd polynomial is zero, in particular

ψ 1 2 = ψ 1 3 = ψ 2 1,2 = ψ 2 1,3 = ψ 2 2,3 = 0.
With those computations, the moments ψ 1 and ψ 2 actually reads

ψ 1 = ψ 1 1 e 1 , ψ 2 = Diag ψ 2 1,1 , ψ 2 2,2 , ψ 2 3,3 ,
and ψ 1 is therefore an eigenvector of ψ 2 . Using Notations 1 leads to write N 1 and N 2 under the form (23a), and one may obsereve that the eigenvectors of N 2 -N 1 ⊗ N 1 are along the cartesian axis e i . Using again evenness of ψ 1 , one obtains Ω i tr(Ω ⊗ Ω )ψ 1 (Ω )dΩ = ψ 1 i .

This leads to write N 3 under the form (24a).

Proposition 4 Consider realizable moments (ψ 0 , ψ 1 , ψ 2 ) ∈ R 2 such that ψ 1 is an eigenvector of ψ 2 . Then the rotated normalized moments (N 1 , N 2 ) given by [START_REF] Dubroca | Hiérarchie des modèles aux moments pour le transfert radiatif[END_REF] are in H 1 .

Proof Under those hypothesis, the decomposition (17) can be simplified. Indeed, since ψ 1 is an eigenvector of ψ 2 , then N 1 is an eigenvector of N 2 . So a rotation R diagonalizing N 2 will send N 1 onto one of the cartesian axis (chose R such that N 1 is along e 1 ). Then this rotation also diagonalizes N 2 -N 1 ⊗ N 1 since it diagonalizes both N 2 and N 1 ⊗ N 1 and one can write N 1 and N 2 under the form (21a).

Finally one can prove that the unique exponential representation for moments (N 1 , N 2 ) satisfying (23a) is (21a) by using Theorem 1 with m(Ω ) = (Ω 1 , Ω 2 1 , Ω 2 2 , Ω 2 3 ). Indeed this theorem provide the existence of a unique function ψ of the form (21a) satisfying

(N 1 1 , N 2 1,1 , N 2 2,2 , N 2 3,3 ) = mψ .
Computing the other moments of such a function (21a) read

ψ 1 2 = ψ 1 3 = ψ 2 1,2 = ψ 2 1,3 = ψ 2 2,3 = 0,
i.e. it satisfies the other moment constraints. Then the unique function [START_REF] Borwein | Partially finite convex programming: Part II[END_REF] satisfying all the moment constraints has the form (21a).

A.2 Computation of moments in H 2

Consider λ ∈ L 2 , and the associated ansatz ψ 2 (defined in (21b)) and its moments (ψ 1 , ψ 2 ) given by ψ 2 (Ω ) = exp λ 5 + λ 1 Ω 1 + (λ 4 -λ 5 )Ω 2 1 , (ψ 1 , ψ 2 ) = (Ω , Ω ⊗ Ω )ψ 2 ∈ H 2 .

The computations of the previous subsection hold. Since ψ 2 does not depend on Ω 2 nor Ω 3 , one deduces that This leads to write N 3 under the form (24b).

Proposition 5 Consider realizable moments (ψ 0 , ψ 1 , ψ 2 ) ∈ R 2 such that ψ 1 is an eigenvector of ψ 2 and ψ 2 2,2 = ψ 2 3,3 . Then the rotated normalized moments (N 1 , N 2 ) given by [START_REF] Dubroca | Hiérarchie des modèles aux moments pour le transfert radiatif[END_REF] are in H 2 .

The proof is identical to the one of Proposition 4 with m(Ω ) = (1, Ω 1 , Ω 2 1 ).

A.3 Computation of moments in H 3

Consider λ ∈ L 3 , and the associated ansatz ψ 3 (defined in (21c)) and its moments (ψ 1 , ψ 2 ) given by ψ 3 (Ω ) = exp (λ 4 + λ 1 Ω 1 ) , (ψ 1 , ψ 2 ) = (Ω , Ω ⊗ Ω )ψ 3 ∈ H 2 .

The computations of the previous subsections hold. In this case, the ansatz ψ 3 is the M 1 ansatz defined in (13), and therefore ψ 2 1,1 is the Eddington factor χ 2 defined in Subsection 2.3. Using Notations 1 leads to write N 1 and N 2 under the form (23c).

The form of N 3 is not simplified compared to the previous case.

Proposition 6 Consider realizable moments (ψ 0 , ψ 1 , ψ 2 ) ∈ R 2 such that ψ 1 is an eigenvector of ψ 2 , ψ 2 2,2 = ψ 2 3,3 and ψ 2 1,1 = ψ 0 χ 2 |ψ 1 |/ψ 0 . Then the rotated normalized moments (N 1 , N 2 ) given by [START_REF] Dubroca | Hiérarchie des modèles aux moments pour le transfert radiatif[END_REF] are in H 2 .

The proof is identical to the one of Proposition 4 with m(Ω ) = (1, Ω 1 ).

Example 3 A

 3 distribution of two beams in directions e 1 and e 2 ψ = δ (Ω .e 1 -1) + δ (Ω .e 2 -1) has moments

Fig. 1

 1 Fig. 1 Unit sphere colored by a distribution function associated to a vector of H 3

Fig. 2 Fig. 3

 23 Fig. 2 Distribution function associated to a vector of H 3 along its preferred axis

Fig. 5 2 Fig. 6

 526 Fig. 5 Unit sphere colored by a distribution function associated to a vector of H 1 \H 2

Fig. 10

 10 Fig.[START_REF] Borwein | Partially finite convex programming: Part II[END_REF] Normalized dose produced by a 10 MeV electron beam in water using a kinetic, M 1 and M 2 solvers with the closures obtained from the minimization procedure and approximations.

  with at x = 0 cm : ε 0 = 10 MeV, Ω p = e 1 , c e = 200 and c o = 10 3 ; at x = 8 cm : ε 0 = 10 MeV, Ω p = -e 1 , c e = 200 and c o = 10 3 .

Fig. 11

 11 Fig.11Normalized dose produced by two 10 MeV electron beams in water using kinetic, M 1 and M 2 solvers with the closures obtained from the minimization procedure and approximations.

Fig. 12

 12 Fig. 12 Dose produced by a 10 MeV electron beam in water using the PENELOPE Monte-Carlo solver (top) and the M 1 (below left) and M 2 (below right) solvers using the approximated closures.

Fig. 13

 13 Fig.[START_REF] Chandrasekhar | Radiative transfer[END_REF] Absolute error between the normalized doses obtained with Monte-Carlo solver and the M 1 (left) and M 2 (right) solvers using the approximated closures.

Fig. 14

 14 Fig. 14 Dose produced by two 10 MeV electron beams in water using the PENELOPE Monte-Carlo solver (top), a M 1 solver (middle left) and a M 2 solver (middle right), a double-M 1 solver (below left) and a double-M 2 solver (below right), with the approximated closures.

Fig. 15

 15 Fig. 15 Absolute error between the doses obtained with Monte-Carlo solver and the M 1 (top left), M 2 (top right), double-M 1 (below left) and double-M 2 solvers (below right) using the approximated closures.

Dose Monte-Carlo Dose M 1 Dose M 2 Fig. 16

 216 Fig. 16 Dose produced by a 10 MeV electron beam in a cut of a chest using PENELOPE Monte-Carlo solver (left), a M 1 solver (middle) and a M 2 solver (right), with the approximated closures.

  that tr(Ω ⊗ Ω ) = 1, one obtains that 3

λ 5 + λ 1 Ω 1 + (λ 4 -λ 5 λ 5 + λ 1 Ω 1 + (λ 4 -λ 5 3 . 2 Ω 1 2 Ω 1

 511455114532121 Using Notations 1 leads to write N 1 and N 2 under the form (23b). Similarily, one has Ω 2 2 exp (λ4 + λ 1 Ω 1 ) dΩ = Ω 2 3 exp (λ 4 + λ 1 Ω 1 ) dΩ =

Table 2

 2 10, and the computation times for this test case are gathered in Table2. This simple case shows that Computation times with the different models for the simple 10 MeV beam test case in 1D

	Model	Kinetic	minimization M 1	minimization M 2	approximated M 1	approximated M 2
	Times	43.76 sec	4.58 sec	8.84 sec	0.016 sec	0.047 sec

Acknowledgements

The authors are grateful to Gabriele Birindelli (Université de Bordeaux) and Kerstin Küpper (RWTH Aachen) for performing the PENELOPE simulations, and to Nuria Escobar Corral (Uniklinikum RWTH Aachen) for providing the density map in the chest for the last test case. T. Pichard's PhD was funded by IdEx Bordeaux and Aquitaine Region. This work was also partially funded by Aquitaine Region and FEDER fund through IOPRA interface.

end, those approximations are

For N 3 1,2,3 , several other linear combinations were possible, but we found that this simple expression gave the best approximation of the M 2 closure. The discrete L ∞ error compared to the solutions of the minimization problem (12) for 3.2 × 10 6 values of (N 1 , N 2 ) ∈ R T is of 3.12.10 -2 . Those values are obtained from 20 values of the order of the model raises. Similarly to the previous case, the M 2 results are close to the kinetic reference, except in the middle of medium, while the M 1 results are less precise. Finally the double-M N models gives satisfactory results when high precision is required. The accuracy of the approximations of the closures of the M 1 , M 2 , double M 1 and double M 2 models is characterized by the errors on the doses gathered in Table 5.

Simple inclined beam in 2D

In 2D, the computations need more computational power. We could not afford to use a minimization procedure in 2D to obtain the M N closure, so only results using the approximated closures are presented. This test consists of a beam of 10 MeV electrons in water, whose direction forms an angle of π/3 with the normal to the boundary. It is modeled by [START_REF] Levermore | Relating Eddington factors to flux limiters[END_REF] with the parameters at x = 0 cm, for y ∈ [4.75 cm, 5.25 cm] :

, c e = 200 and c o = 10 3 .

In order to compare our results with those of PENELOPE Monte-Carlo code, the dose is normalized by the maximum dose, i.e. we computed the percentage depth dose (PDD)

.

The domain was meshed with 600×600 cells in position and 100 cells in energy.

The doses obtained with the different models are given on Fig. 12, and the computation times for that test case are given in Table 6 is overly diffusive. However, the effect in multi-D is lower than in 1D. The M 2 results are


Chloé Pasin

Laura Villain

François Dufour

Daniel Commenges

Mélanie Prague

Rodolphe Thiébaut

Use of mathematical modeling for optimizing and adapting immunotherapy protocols in HIV-infected patients

come

Use of mathematical modeling for optimizing and adapting immunotherapy protocols in HIV-infected patients

Context

Infection by the Human Immunodeficiency Virus (HIV) compromises the immune system: depletion of CD4 `T lymphocytes Combination antiretroviral therapies (cART) ùñ improvment of patient's survival Some HIV-infected patients under cART are unable to recover normal CD4 `T cell levels Immune therapy using interleukine 7 (IL7) complements cART: § Phase I trials and observational studies ùñ safety and beneficial effect of IL7 injections § Phase I/II human clinical trials (INSPIRE 1, 2 and 3 studies) ùñ repeated cycles of 3 weekly IL7 injections help maintaining CD4 Main issue: criterion in the original protocol (decision to begin a new cycle when CD4 levels ă 550 cells/µ) inadequate for some patients Next step: optimizing and adapting protocols of injections Aim: maintaining the CD4 levels ą 500 cells/µL without using too many injections

Pipeline

Inclusion of a new patient

Estimation of its individual parameters using the first observations.

Parameters estimation

Estimation of the ODE parameters θ using INSPIRE data ùñ determination of a posterior distribution.

Optimal control approach

Use of optimal control theory to determine an optimal strategy of injections from a cost function.

Stochasticity attributed to the biological process.

Bayesian approach

Treatment decision is adapted at each new observation, by updating individual parameters estimation with MCMC algorithm.

Stochasticity attributed to the uncertainty on the parameters estimation.

Assessment of the methods

Pseudo-patients parameters generation

Sample N parameter sets from the posterior distribution of θ.

Parameters estimation

Estimation of the ODE parameters θ using INSPIRE data ùñ determination of a posterior distribution.

Optimal control approach

Use of optimal control theory to determine an optimal strategy of injections from a cost function.

Stochasticity attributed to the biological process.

N=50, horizon = 1 year

Bayesian approach

Treatment decision is adapted at each new observation, by updating individual parameters estimation with MCMC algorithm.

Stochasticity attributed to the uncertainty on the parameters estimation.

N=150

  determine actions (injection of not and choice of dose) at given time points over an horizon of time ùñ impulse control problem in the optimal control theory Two-steps method: § Developing an adapted mathematical model for the process: Piecewise Deterministic Markov model[Davis (1984)] § Using the theory of optimal control from[Costa et al. (2016)] to solve the impulse control problem by iteration of an integro-differential operator computation of the optimal cost determination of an optimal strategy of injections comparison of the obtained optimal strategy to other clinical protocols C. Pasin, F. Dufour, L. Villain, H. Zhang, R. Thiébaut. Controlling IL-7 injections in HIV-infected patients. To appear in Bulletin of Mathematical Biology Pasin, Villain, Dufour, Commenges, Prague, Thiébaut Optimizing and adapting immunotherapy protocols 7 / 15 Optimal control: method Jump in the process : modification of parameter π § Deterministic: CD4 levels ă 500 cells/µL ùñ IL7 injection § Random: stochastic time τ of effect of an IL7 injection π " π0 `βpnq π d 0.25 1 ttPrt inj ,t inj `τ suCost function: combination of two criteria § Time spent with CD4 levels ă 500 cells/µL § Number of injections realized Application of the method to determine the optimal strategy of 50 pseudo-patients on a reduced model: § 2 possible doses (10,20) : 3-injections cycles. P2: 3-injections cycle then 2-injections cycles. P3: 2-injections cycles. P4: 2-injections cycle then 1-injection cycles. P5: 1-injection cycles. Pasin, Villain, Dufour, Commenges, Prague, Thiébaut Optimizing and adapting immunotherapy protocols Optimal control: conclusion Development of a dynamic programming method to apply the theory of optimal control in a biological framework Intuitive optimal strategy: first cycles of 2 injections to increase the number of CD4 `then cycles of 1 injection to maintain the CD4 levels Limitations: § One year horizon only due to computational time § Main hypothesis: known individual parameters Pasin, Villain, Dufour, Commenges, Prague, Thiébaut Optimizing and adapting immunotherapy protocols 10 / 15

  Regular measurements of CD4 `levels and proliferation marker Ki67 Visits every three months. CD4 `levels ă 550 cells/µL ùñ new cycle Total duration of the studies: between 12 and 24 months

	Mechanistic model
	Dynamics of CD4 `T lymphocytes, resting (R) or proliferating (P) [Thiébaut et al.
	(2014)]
		µ R		µ P
					$ ' &	dR dt	" λ ´µR R ´πR `2ρP
	λ	R	π ρ	P	' % initial condition = equilibrium dP " ´µP P ´ρP `πR dt
	Parameters estimation with a population approach [Prague et al. (2013)]
				πi " π0	`βpnq π d 0.25 i	1 ttPrt i inj ,t i inj	`τ i su

Total of 128 HIV-infected patients under cART with CD4 levels between 100-400 cells/µL and undetectable viral load for at least 6 months Cycles of 3 weekly IL7 injections of dose 20µg/kg (21 received dose 10 or 30) Identification of random effects on parameters λ and ρ IL-7 injection: effect on the proliferation rate π (differs for each injection within a cycle) [Jarne et al. (2017)] d : injected dose, n P t1, 2, 3u : injection of the cycle, τ : time of effect of the injection Pasin, Villain, Dufour, Commenges, Prague, Thiébaut Optimizing and adapting immunotherapy protocols 3 / 15

Pasin, Villain, Dufour, Commenges, Prague, ThiébautOptimizing and adapting immunotherapy protocols

Acknowledgments

INRIA SISTM / CQFD teams and other collaborators Vaccine Research Institute (Y Levy, JD Lelievre) INSPIRE Study group (JP Routy, I Sereti, M Fischl, P Ive, RF Speck, G D'offizi, S Casari, S Foulkes, V Natarajan, T Croughs, Jean-F Delfraissy, G Tambussi, MM Lederman)