Pauline Godillon-Lafitte

Thierry Goudon
email: thierry.goudon@math.univ-lille1.fr

UMR Labo Paul Painlevé

A Coupled Model for Radiative Transfer: Doppler Effects, Equilibrium and Non-Equilibrium Diffusion Asymptotics

Keywords: Hydrodynamic limits, Diffusion approximation, Radiative transfer, Doppler correction, AMS Subject classification, 35Q99 35B25

published or not. The documents may come

A coupled model for radiative transfer: Doppler effects, equilibrium, and nonequilibrium diffusion asymptotics

Introduction

This paper is devoted to the asymptotic analysis of a system of PDEs arising in the modeling of radiative transfer. The unknowns are the (nonnegative) specific intensity of radiation f (t, x, v), which depends on variables of time t ≥ 0, space x ∈ R 3 and direction v ∈ S 2 , and the (nonnegative) material temperature θ(t, x). The former satisfies a kinetic equation, the latter a drift-diffusion equation, and the coupling arises from source terms. Let us introduce the problem by dealing with a simple dimensionless example :

   ε∂ t f ε + v • ∇ x f ε = Q ε (f ε , θ ε ) in R + t × R 3 x × S 2 v , ∂ t θ ε + div x (uθ ε ) -D∆ x θ ε = - S 2 Λ ε γ ε Q ε (f ε , θ ε ) dv in R + t × R 3 x , (1)
where D is the positive diffusion constant, and u : (t, x) ∈ R × R 3 → R 3 is a given velocity field. We wish to investigate the behavior of the solution (f ε , θ ε ) as the positive small parameter ε, which derives from physical quantities, tends to 0. The weights Λ ε and γ ε which arise in these formulae are related to Doppler corrections. They take into account the fact that certain quantities are measured in a mobile frame, instead of being evaluated directly in the reference frame. They depend on the field u as follows

Λ ε (t, x, v) = 1 -εu(t, x) • v 1 -ε 2 |u(t, x)| 2 , γ ε (t, x) = 1 1 -ε 2 |u(t, x)| 2 . ( 2
)
The source term in (1) splits into two parts

Q ε (f ε , θ ε ) = L s Q s,ε (f ε ) + L a Q a,ε (f ε , θ ε ),
1 involving parameters L s and L a which depend on ε as it will be specified later on. Indeed, photons are submitted to two kinds of interaction processes, embodied in the operators Q s,ε and Q a,ε respectively: scattering and absorption/emission. They are defined as follows

       Q s,ε (f )(v) = σ s 1 Λ ε (v) 3 S 2 Λ ε (v ′ ) 2 f (v ′ ) dv ′ -Λ ε (v)f (v) , Q a,ε (f, θ) = σ a Aθ 4 Λ 3 ε -Λ ε f , (3)
with given constant coefficients A > 0, σ s > 0 and σ a > 0. (Non constant kernels will be considered in the following sections.) Moreover, initial data are prescribed

f ε,|t=0 = f ε,0 ≥ 0, θ ε,|t=0 = θ ε,0 ≥ 0.
Details on the physical background and discussion on the scaling will be given in Section 2. This work relies on the analysis of diffusion asymptotics of kinetic equations. There exists a wide literature on the topics, with various viewpoints and fields of applications. Let us mention the papers of Bardos-Golse-Perthame [START_REF] Bardos | The Rosseland approximation for the radiative transfer equations[END_REF] and Bardos-Golse-Perthame-Sentis [START_REF] Bardos | The nonaccretive radiative transfer equations: existence of solutions ans Rosseland approximations[END_REF] for applications to radiative transfer, Poupaud [START_REF] Poupaud | Diffusion approximation of the linear semiconductor Boltzmann equation: analysis of boundary layers[END_REF] and Golse-Poupaud [START_REF] Golse | Limite fluide des équations de Boltzmann des semi-conducteurs pour une statistique de Fermi-Dirac[END_REF] for applications to semi-conductors theory, Lions-Toscani [START_REF] Lions | Diffuse limit for finite velocity Boltzmann kinetic models[END_REF] for special arguments when dealing with discrete velocity models, Degond-Goudon-Poupaud [START_REF] Degond | Diffusion limit for non homogeneous and non reversible processes[END_REF] when considering collision operators without microreversibilty assumptions, Chalub-Markowich-Perthame-Schmeiser [START_REF] Chalub | Kinetic models for chemotaxis and their drift-diffusion limits[END_REF] for applications to chemotaxis... We also refer to the survey of Golse [START_REF] Golse | From kinetic to macroscopic models in Kinetic equations and asymptotic theory[END_REF]. In our work, the interesting features are the following:

-at first, the collision operator splits into two parts. Then, depending on the dominant effect -scattering or absorption, depending on how L s and L a depend on ε -it can relax to different equilibrium states. Accordingly, we can be led to a limit process where the radiations have a temperature which differs from the material temperature. This is the so-called non-equilibrium regime. In turn, we obtain a coupled system of diffusion equations, see [START_REF] Battaner | Astrophysical fluid dynamics[END_REF] below. Note that this aspect has also been investigated in a different framework by Dogbe [START_REF] Dogbe | The radiative transfer equations: diffusion approximation under accretiveness and compactness assumptions[END_REF]; -secondly, the Doppler corrections can induce in the limit ε → 0 some additional convection terms. This is due to the fact that, for ε > 0, the equilibrium states of the scattering operator have a non-vanishing flux of order O(ε);

-thirdly, our proof strategy differs slightly from those in [START_REF] Bardos | The nonaccretive radiative transfer equations: existence of solutions ans Rosseland approximations[END_REF], [START_REF] Dogbe | The radiative transfer equations: diffusion approximation under accretiveness and compactness assumptions[END_REF]. We use in a more systematic way the dissipative properties of the system to obtain useful estimates and compactness properties. Furthermore, strong convergence properties are established via a compensated compactness argument, namely the Div-Curl lemma, see [START_REF] Tartar | Remarks on homogenization in Homogenization and effective moduli of material and media[END_REF], instead of applying average lemma techniques (see, e.g. [START_REF] Bardos | The nonaccretive radiative transfer equations: existence of solutions ans Rosseland approximations[END_REF], [START_REF] Golse | Regularity of the moments of the solution of a transport equation[END_REF], [START_REF] Perna | L p regularity of velocity averages[END_REF]...).

To be more precise, we shall prove the following results:

-Equilibrium Regime: if L a = 1/ε and L s = ε, then f ε and θ ε converge (in some sense...) to Aθ 4 and θ, respectively, where θ is a solution of the following non linear drift-diffusion equation

∂ t (θ + Aθ 4 ) + div x u(θ + Aθ 4 ) -D∇ x θ - A 3σ a ∇ x θ 4 + A 3 θ 4 div x (u) = 0. (4)
-Non-Equilibrium Regime: if L s = 1/ε and L a = ε, then f ε (respectively θ ε ) converges (in some sense...) to ρ (respectively θ), which depends only on the time and space variables, with (ρ, θ) solution of the following coupled system    ∂ t ρ -1 3σ s ∆ x ρ + div x (ρu) + 1 3 ρdiv x u = σ a (Aθ 4 -ρ), ∂ t θ + div x (uθ) -D∆ x θ = σ a (ρ -Aθ 4 ).

(

) 5
This kind of behavior (additional convection terms, non-equilibrium asymptotics...) is known by physicists (see e.g. Lowrie-Morel-Hittinger [START_REF] Lowrie | The coupling of radiation and hydrodynamics[END_REF]). It has been pointed out more recently on a mathematical viewpoint by Buet-Després [START_REF] Ch | Asymptotic analysis of fluid models for the coupling of radiation and hydrodynamics[END_REF], who describe much more complicated models than [START_REF] Aregba-Driollet | Convergence of relaxation schemes for conservation laws[END_REF], involving a coupling to the Euler equations. In particular, they discuss the approximations by moment closure and their connections to shock relations. Here, our goal is to obtain a complete and rigorous convergence result illustrating these interesting features, at least for a somewhat simplified situation.

The paper is organized as follows. In Section 2, we introduce in detail the coupling of radiative transfer with hydrodynamics. Then, we make some drastic simplifications in order to obtain a model we are able to deal with, from a mathematical viewpoint. This work is a new step towards more realistic models in radiative transfer. In Section 3, we discuss the dimension analysis, identifying relevant scaling parameters by means of physical quantities. Section 4 is devoted to a preliminary discussion of the properties of the scattering operator, which, in turn, leads to dissipative properties of the coupled system of PDEs. There, we give also a complete statement of the main results of the paper. The proof for the non-equilibrium regime is postponed to Section 5 and in Section 6, we prove the results for the equilibrium regime. Section 7 is devoted to the numerical simulations that we obtain for the non-equilibrium regime using splitting techniques, that fit this situation, since the (stiff) source term can be split into a vanishing mean value part and a nonzero mean value part : we obtain very satisfactory results when the bulk velocity is a rarefaction wave and stretch our computations to cases that do not satisfy the theoretical assumptions. We also study the influence of the variations of the coefficients σ a and σ s . The paper ends with an appendix containing some technical details.

Physical Background and Motivation

First of all, let us introduce the equations of radiative transfer, taking into account relativistic effects, Doppler corrections and the coupling of radiation with hydrodynamics. We refer for further details on the physics to the classical treatises of Battaner [START_REF] Battaner | Astrophysical fluid dynamics[END_REF], Mihalas-Mihalas [START_REF] Mihalas | Foundations of radiation hydrodynamics[END_REF] (see also the notes of Rutten [START_REF] Rutten | Radiative transfer in stellar atmospheres[END_REF]) and also to the recent papers by Lowrie-Morel-Hittinger [START_REF] Lowrie | The coupling of radiation and hydrodynamics[END_REF] and Buet-Després [START_REF] Ch | Asymptotic analysis of fluid models for the coupling of radiation and hydrodynamics[END_REF]. We adopt an Eulerian viewpoint. Accordingly, we take into account the motion of the fluid with respect to the reference frame which induces corrective terms when dealing with quantities measured in a comobile frame. Secondly, we will present some simplifications of the model.

Coupling of Radiation and Hydrodynamics

Let f (t, x, v, ν) stand for the specific intensity of radiation. It depends on the time t, the position x ∈ R 3 , the solid angle v ∈ S 2 and the frequency ν ∈ [0, +∞) and has the dimension of an energy per surface, time and frequency units. Precisely, let Ω ⊂ R 3 , D ⊂ S 2 and N ⊂ R + . Denote by c the speed of light. Then, the integral 1 c Ω D N f dν dv dx gives, at time t, the radiant energy corresponding to the set of photons with frequencies in N , directions within a solid angle in D, and located at positions in Ω. It can be related to the distribution function in the phase space F (t, x, p), which gives the number of photons per volume unit of the phase space. For the momentum variable, we have p = hν c v (and dp = ( h c ) 3 ν 2 dν dv), with h the Planck constant. It follows that

hν F (t, x, p) dp dx = hν F t, x, hν c v h 3 c 3 ν 2 dν dv dx = 1 c f (t, x, v, ν) dν dv dx
holds. Throughout the paper dv stands for the normalized Lebesgue measure on S 2 . In other words, we set

f (t, x, v, ν) = h 4 ν 3 c 2 F t, x, hν c v .
Note that all photons are travelling with velocities having the same modulus c. The quantity f satisfies the following evolution equation

1 c ∂ t f + v • ∇ x f = Q(f )
which links the transport of photons with velocity cv along the straight lines x + ctv to the interaction processes submitted by the photons, that is scattering, absorption and emission phenomena. When describing such phenomena, we want to take into account Doppler corrections. To this end, we introduce the velocity of the surrounding fluid u(t, x) and the Lorentz factor defined by

γ(t, x) = 1 1 -|u(t, x)| 2 /c 2 .
In what follows, we denote by the 0 superscript the quantities which are evaluated in the comobile frame.

In particular, the formulae

ν 0 = νγ 1 - u • v c , v 0 = ν ν 0 v - γ c u 1 - u • v c γ γ + 1 (6)
link the frequency and the direction, respectively, of the photons measured in the comobile frame to the quantities ν and v measured in the reference frame, see [START_REF] Mihalas | Foundations of radiation hydrodynamics[END_REF]. A key property is the invariance relation

ν 0 dν 0 dv 0 = ν dν dv. ( 7
)
In the sequel, we will also use the notation

Λ(t, x, v) = ν 0 ν = 1 -v • u(t, x)/c 1 -|u(t, x)| 2 /c 2 ,
and we will often omit the time and space variables when no confusion can arise. Of course, we come back to (2) by considering a scaling where c is large compared to the characteristic velocity of the fluid.

Scattering Operator

In the comobile frame, scattering interactions only produce a change in the direction of the trajectories of the photons. Let ℓ s be the scattering mean free path. It describes the average distance between successive scattering events. Let f 0 stand for the intensity of radiation measured in the comobile frame. Then, considering given frequency ν 0 and direction v 0 , scattering effects are described by a gain term

1 ℓ s Q 0,+ s (f 0 ) = 1 ℓ s S 2 σ 0 s (ν 0 , v 0 , v 0 ′ )f 0 (ν 0 , v 0 ′ ) dv 0 ′ ,
and a loss term

1 ℓ s Q 0,- s (f 0 ) = 1 ℓ s Σ 0 s (ν 0 , v 0 )f (ν 0 , v 0 ).
These definitions involve nonnegative (dimensionless) coefficients σ 0 s (ν 0 , v 0 , v 0 ′ ), which characterize the change of direction from v 0 ′ to v 0 , and Σ 0 s . These coefficients may also depend on the position x. The gain term accounts for the photons that change from the states (v 0 ′ , ν 0 ′ ) to (v 0 , ν 0 ) due to the scattering phenomena, while the loss term accounts for the change of photons in the state (v 0 , ν 0 ). The crucial assumption is that the scattering operator

Q 0 s = Q 0,+ s -Q 0,- s is conservative that is ∞ 0 S 2 Q 0 s dv 0 dν 0 = 0.
In other words, when seen in the comobile frame, scattering has no contribution to the total energy balance. Obviously, this is satisfied by imposing the relation

Σ 0 s (ν 0 , v 0 ) = S 2 σ 0 s (ν 0 , v 0 ′ , v 0 ) dv 0 ′ ,
which will be assumed from now on.

Then, we go back to the reference frame by using the formulae (6) and the conversion relations for which we refer to [START_REF] Mihalas | Foundations of radiation hydrodynamics[END_REF])

1 ν 3 f (ν, v) = 1 (ν 0 ) 3 f 0 (ν 0 , v 0 ), Q s (ν, v) = ν ν 0 2 Q 0 s (ν, v). (8)
Let us give at once the expression of the collision operator Q s , as a linear operator acting on f , postponing the detailed computations to Appendix B. Consider a triple (ν, v, v ′ ) ∈ R + × S 2 × S 2 . First, we associate to (ν, v) the pair (ν 0 , v 0 ) ∈ R + × S 2 by using [START_REF] Bouchut | Diffusive BGK approximations for nonlinear multidimensional parabolic equations[END_REF]. Then, we set

ν ′ = ν 0 γ(1 -v ′ • u/c) = ν 1 -v • u/c 1 -v ′ • u/c , (9)
and we simply denote σ s (ν, v, v ′ ) = σ 0 s (ν 0 , v 0 , v 0 ′ ) with v 0 ′ defined by ( 6) from (ν ′ , v ′ ). Finally, we obtain

Q s (f )(ν, v) = 1 Λ(v) 2 S 2 σ s (ν, v, v ′ )Λ(v ′ )f (ν ′ , v ′ ) dv ′ - S 2 σ s (ν, v ′ , v) Λ(v ′ ) 2 dv ′ Λ(v)f (ν, v). ( 10
)
Remark 1 The special case where the scattering kernel σ s is constant leads to some simplifications. Indeed, we note that S 2 dv Λ 2 = 1. Therefore, for such a kernel, we get

Q s (f )(ν, v) = σ s 1 Λ(v) 2 S 2 Λ(v ′ )f (ν ′ , v ′ ) dv ′ -Λ(v)f (ν, v) . ( 11
)
Remark 2 It is worth pointing out that the conservation property of the scattering phenomena only holds in the comobile frame: we have ∞ 0 S 2 Q s dv dν = 0. We will come back to this aspect later on (see Proposition 1).

Absorption-Emission Operator

Absorption consists of a loss of photons. It involves a coefficient σ a (ν, v) (possibly depending also on the space variable) and a mean free path ℓ a ; namely it gives the following loss term

- 1 ℓ a σ a (ν, v) ν 0 ν f (ν, v).
Moreover, photons are produced when the energy stored into atoms or molecules is transferred to the radiation field. It gives the following gain term

1 ℓ a σ a (ν, v) ν ν 0 2 B(ν 0 , θ)
where the function B characterizes the emission law, depending on the temperature of the material θ, which depends on the time and space variables. For instance, when photons are emitted in a thermodynamic equilibrium with a black body system, B is given by the following Planck function

B(ν, θ) = 2hν 3 c 2 1 e hν/(kθ) -1 ( 12
)
where k is the Boltzmann constant.

Coupled Equations

To sum up, f satisfies the following equation

1 c ∂ t f + v • ∇ x f = Q(f, θ) = 1 ℓ s Q s (f ) + 1 ℓ a Q a (f, θ) = 1 ℓ s 1 Λ(v) 2 S 2 σ s (ν, v, v ′ )Λ(v ′ ) f (ν ′ , v ′ ) dv ′ - S 2 σ s (ν, v ′ , v) Λ(v ′ ) 2 dv ′ Λ(v) f (ν, v) + σ a (ν, v) ℓ a 1 Λ(v) 2 2hν 3 c 2 B(ν 0 , θ) -Λ(v) f (ν, v) . (13)
Next, this equation has to be considered as coupled to hydrodynamic equations for the fluid. The unknowns are the density n, the bulk velocity u and the energy E. These quantities satisfy the relativistic Euler equations

         ∂ t n + div x (nu) = 0, ∂ t (γHnu) + Div x (γHnu ⊗ u + pI) = - 1 c ∞ 0 S 2 vQ(f ) dv dν, ∂ t (γHn -p/c 2 ) + div x (γHnu) = - ∞ 0 S 2 Q(f ) dv dν, ( 14
)
where p is the pressure and the enthalpy H is defined by

H = 1 + e + γp/n c 2 ,
with e the specific internal energy of the fluid. Assuming the gas follows a simple perfect gas pressure law, we get p = Γen/γ and H = 1 + (1 + Γ)e/c 2 , with Γ the adiabatic constant. We are usually interested in physical situations where |u|/c ≪ 1; hence terms of order 2 with respect to |u|/c are neglected. In such a case, the left-hand side in ( 14) is replaced by this of the classical Euler equations, see [START_REF] Lowrie | The coupling of radiation and hydrodynamics[END_REF], [START_REF] Ch | Asymptotic analysis of fluid models for the coupling of radiation and hydrodynamics[END_REF], and we get

         ∂ t n + div x (nu) = 0, ∂ t (nu) + Div x (nu ⊗ u + pI) = - 1 c ∞ 0 S 2 vQ(f ) dv dν, ∂ t (nE) + div x (nEu + pu) = - ∞ 0 S 2 Q(f ) dv dν, (15)
with E = e + u 2 /2 the total energy. The coupling between hydrodynamics and radiation arises from the right-hand sides in [START_REF] Golse | Regularity of the moments of the solution of a transport equation[END_REF] (or [START_REF] Golse | Limite fluide des équations de Boltzmann des semi-conducteurs pour une statistique de Fermi-Dirac[END_REF]). They describe momentum and energy exchanges between the material and the photons, since ∞ 0 S 2 Q(f ) dv dν is the energy balance of the radiation while

∞ 0 S 2 v c Q(f )
dv dν is the momentum balance of the radiation. In particular, it is worth pointing out that the total energy of the system is conserved d dt

1 c R 3 ∞ 0 S 2 f dv dν dx + R 3 nE dx = 0. ( 16
)
2.2 Simplified Models

Grey Assumption

The first assumption that simplifies the model consists in assuming that the coefficients do not depend on the frequency variable ν. This is the so-called "grey assumption". This can be reasonable for the scattering coefficient -considering for instance Thomson scattering, see e.g. [START_REF] Battaner | Astrophysical fluid dynamics[END_REF] -but it is more questionable for the absorption coefficient. Therefore, we can integrate equation [START_REF] Golse | From kinetic to macroscopic models in Kinetic equations and asymptotic theory[END_REF] with respect to ν. Let us introduce the following notations

       f (t, x, v) = ∞ 0 f (t, x, ν, v) dν, 1 Λ B(θ) = ∞ 0 B(ν 0 , θ) dν = ∞ 0 B(ν 0 , θ) dν 0 Λ .
For instance, considering the black body emission law [START_REF] Dogbe | The radiative transfer equations: diffusion approximation under accretiveness and compactness assumptions[END_REF], we get

B(θ) = Λ ∞ 0 2h(ν 0 ) 3 c 2 dν e hν 0 /(kθ) -1 = Λ 2h c 2 ∞ 0 kθ h 3 µ 3 e µ -1 kθ hΛ dµ = σ π θ 4
where σ = 2π 5 k 4 15h 3 c 2 is the Stefan-Boltzmann constant. Furthermore, denoting Λ(t, x, v ′ ) by Λ ′ for the sake of simplicity and recalling [START_REF] Degond | Diffusion limit for non homogeneous and non reversible processes[END_REF], integrating the scattering gain term leads to

1 Λ 2 ∞ 0 S 2 σ s (v, v ′ )Λ ′ f (ν ′ , v ′ ) dv ′ dν = 1 Λ 2 S 2 σ s (v, v ′ )Λ ′ ∞ 0 f (ν ′ , v ′ ) Λ ′ Λ dν ′ dv ′ = 1 Λ 3 S 2 σ s (v, v ′ )Λ ′2 f (v ′ ) dv ′ .
Hence, from now on we drop the bars and we denote by • the integration with respect to the angular variable v. We are thus led to consider the following kinetic equation

1 c ∂ t f + v • ∇ x f = Q(f ) = 1 ℓ s σ s Λ 2 f Λ 3 - σ s Λ 2 Λf + 1 ℓ a σ a B(θ) Λ 3 -Λf := 1 ℓs Q s (f ) + 1 ℓa Q a (f, θ). (17)

A Simpler Coupling

Next, it is very difficult to deal with a coupled model involving the Euler equations. Instead, we assume that n and u are given (with a constant density), and we only consider a diffusion equation for the material temperature θ. Hence, the kinetic equation ( 17) is simply coupled to the following drift-diffusion equation, whose unknown is the energy 3 2 kθ,

3 2 k ∂ t θ + ∇ x • (uθ) -d∆ x θ = - Λ γ Q(f ) . (18)
In particular, let us look at the evolution of the following energy functional

E(t) = 1 c R 3 S 2 f dv dx + 3 2 k R 3 θ dx.
We have

d dt E = R 3 S 2 1 - Λ γ 1 ℓ s Q s + 1 ℓ a Q a dv dx = R 3 u c • S 2 v 1 ℓ s Q s + 1 ℓ a Q a dv dx. ( 19
)
The right-hand side keeps track of the evolution of the kinetic energy associated to the Euler equation [START_REF] Golse | Limite fluide des équations de Boltzmann des semi-conducteurs pour une statistique de Fermi-Dirac[END_REF], and of the conservation of the total energy in that system, see [START_REF] Goudon | Approximation by homogeneization and diffusion of kinetic equations[END_REF]. This remark explains the role of the corrective term Λ/γ in the right-hand side of (18): the energy production in [START_REF] Lions | Diffuse limit for finite velocity Boltzmann kinetic models[END_REF] is exactly what would come from the kinetic energy balance in [START_REF] Golse | Limite fluide des équations de Boltzmann des semi-conducteurs pour une statistique de Fermi-Dirac[END_REF] if we had consider an evolution equation for the bulk velocity u.

Remark 3 Considering the isotropic case [START_REF] Perna | L p regularity of velocity averages[END_REF] with the grey assumption, we have

Q s (f ) = σ s Λ 2 f Λ 3 -Λf .
With the black body emission law [START_REF] Dogbe | The radiative transfer equations: diffusion approximation under accretiveness and compactness assumptions[END_REF], it leads to the interaction operator that was stated in the introduction, see (1).

Dimension Analysis

We aim at writing the equation in a dimensionless form. To this end, we introduce an observation length unit L, characteristic of the flow behavior. Next, we use the material speed of sound u ∞ as the velocity unit. Thus, we naturally define a time scale by T = L/u ∞ . The material temperature will be scaled according to n ∞ u 2 ∞ /(3k), with n ∞ the density of the fluid, while radiation terms will be evaluated in comparison to a reference temperature τ ∞ . Precisely, we will use the quantity στ 4 ∞ as a dimension unit for f . Therefore, we define starred dimensionless quantities by the following relations

             t = T t ⋆ = L u ∞ t ⋆ , x = L x ⋆ , u(t, x) = u ∞ u ⋆ (t ⋆ , x ⋆ ), θ(t, x) = n ∞ u 2 ∞ 3k θ ⋆ (t ⋆ , x ⋆ ), f (t, x, v) = στ 4 ∞ f ⋆ (t ⋆ , x ⋆ , v), B(θ) = στ 4 ∞ B ⋆ (θ ⋆ ).
For the diffusion coefficient, we assume it scales to d = u ∞ L D, with D > 0.

The rescaled equations will depend on the following dimensionless parameters

C = c u ∞ , P = 2στ 4 ∞ n ∞ u 3 ∞ , L s = L ℓ s , L a = L ℓ a , D = d u ∞ L .
Indeed, performing the change of variables in ( 17) and ( 18) leads to the following dimensionless system, where we dropped the stars,

               1 C ∂ t f + v • ∇ x f = L s σ s Λ 2 f Λ 3 - σ s Λ 2 Λf + L a σ a B(θ) Λ 3 -Λf ∂ t θ + ∇ x • (uθ) -D∆θ = -PL s Λ γ σ s Λ 2 f Λ 3 - σ s Λ 2 Λf -PL a Λ γ σ a B(θ) Λ 3 -Λf (20)
with

Λ(t, x, v) = γ(t, x) 1 - v • u(t, x) C , γ(t, x) = 1 1 -|u(t, x)| 2 /C 2 .
Remark 4 If the momentum evolution is taken into account, the dimensionless form of (15) reads

∂ t (nu) + Div x (nu ⊗ u + pI) = -L s σ s Λ 2 f v Λ 3 - σ s Λ 2 vΛf + L a σ a B(θ) v Λ 3 -vΛf . ( 21
)
In what follows, we deal with asymptotic regimes where

C = 1 ε ≫ 1, P C = 1, D = 1. ( 22
)
It means that we are interested in nonrelativistic flows with a moderate amount of radiation in the flow. Next, it remains to discuss the scaling for the mean free paths. Actually, we can discuss two different asymptotic regimes. The former consists of assuming

L s = ε ≪ 1, L a = 1 ε ≫ 1.
Here, radiations adapt to the temperature of the material since we can naturally guess that f ∼ B(θ). This is the equilibrium regime which yields a diffusion equation for the common temperature. This kind of regime is in the spirit of results presented in [START_REF] Bardos | The Rosseland approximation for the radiative transfer equations[END_REF], [START_REF] Bardos | The nonaccretive radiative transfer equations: existence of solutions ans Rosseland approximations[END_REF]. For the latter we assume

L s = 1 ε ≫ 1, L a = ε ≪ 1.
This is the non-equilibrium regime since f relaxes to a state characterized by a temperature which differs from the material temperature. Accordingly, we are led to a coupled diffusion system (see [START_REF] Dogbe | The radiative transfer equations: diffusion approximation under accretiveness and compactness assumptions[END_REF]). Furthermore, Doppler corrections induce additional convective terms in the limit equation.

Remark 5 Note that the definition of the parameters by means of ε = 1/C can be seen as a constraint on the observation scales with respect to the physical constants. For instance the relation L a,s = 1/ε imposes L 2 /T = cℓ a,s . Note also that our scaling assumptions imply that

n∞u 2 ∞ 2 = στ 4 ∞ c
, which means that the kinetic energy of the fluid balances the radiative energy.

Dissipation Properties; Main Results

Properties of the Scattering Operator; Dissipation

Our analysis relies crucially on the dissipative properties of System [START_REF] Lowrie | The coupling of radiation and hydrodynamics[END_REF]. First, let us consider the scattering operator

Q s (f ) = σ s Λ 2 f Λ 3 - σ s Λ 2 Λf. Proposition 1 Let σ s (v, v ′ ) be a nonnegative function such that σ s (v, v ′ ) = σ s (v ′ , v).
Then, the following properties hold: i) let G : R + → R + be a nondecreasing function. We have

S 2 Q s (f ) G(Λ 4 f ) Λ dv = - S 2 S 2 σ s (v, v ′ ) 2Λ(v) 2 Λ(v ′ ) 2 Λ(v) 4 f (v) -Λ(v ′ ) 4 f (v ′ ) × G(Λ 4 f )(v) -G(Λ 4 f )(v ′ ) dv ′ dv ≤ 0; ii) Ker(Q s ) = Span(Λ -4 ); iii) Λ Q s = 0. Proof. Of course, iii) is a direct consequence of i) when choosing G(z) = 1.
In particular it is worth pointing out that Q s is not a conservative operator, in the sense that Q s = 0. The conservation property of the scattering holds only in the comobile frame. Besides, it is clear that Q s (Λ -4 ) = 0 and the converse is also a direct consequence of i) (choose for instance G(s) = s). Therefore, only i) needs to be proved. We have

S 2 Q s (f ) G(Λ 4 f ) Λ dv = S 2 S 2 σ s (v, v ′ ) Λ(v ′ ) 2 f (v ′ ) Λ(v) 3 - Λ(v)f (v) Λ(v ′ ) 2 dv ′ G(Λ 4 f )(v)Λ(v) dv = S 2 S 2 σ s (v, v ′ ) Λ(v) 2 Λ(v ′ ) 2 Λ(v ′ ) 4 f (v ′ ) -Λ(v) 4 f (v) G(Λ 4 f )(v) dv ′ dv.
We conclude by using the change of variables

(v, v ′ ) → (v ′ , v) which leaves the kernel σs(v,v ′ ) Λ(v) 2 Λ(v ′ ) 2 invariant.
Remark 6 Let us point out that the equilibrium state Λ -4 has a nonvanishing flux

S 2
vΛ -4 dv = 0.

This remark will lead to additional terms in the limit equations.

Remark 7 Since we consider the sum Q(f ) of the interaction processes with weight Λ/γ in the right-hand side of the diffusion equation ( 18), then, by iii), only the absorption term actually contributes :

Λ γ Q(f ) = Λ γ Q a (f, θ) .
Remark 8 If one deals only with the functions G(z) = 1 or G(z) = z it is not necessary to assume the symmetry of the kernel σ s (we refer to [START_REF] Degond | Diffusion limit for non homogeneous and non reversible processes[END_REF] on this aspect).

We aim at deducing some a priori estimates, independent of the parameters C, P, L a , L s . Let us collect here the set of requirements that will be imposed on the data of the problem.

(H0) The velocity field

u belongs to W 1,∞ (R + × R 3 ). (H1) There exist constants σ * , σ * > 0 such that 0 < σ * ≤ σ a (x, v, v ′ ) ≤ σ * , 0 < σ * ≤ σ s (x, v, v ′ ) = σ s (x, v ′ , v) ≤ σ * for almost all (x, v, v ′ ) ∈ R 3 × S 2 × S 2 .
(H2) The function B is nonnegative and increasing. Let us introduce the following auxiliary functions, associated to B

Γ(z) = z 0 B(θ) dθ, G(z) = z 0 B ′ (θ) dθ, J(z) = z 0 B (-1) (s) ds,
where B (-1) stands for the inverse function of B. We suppose that there exist constants p ∈ (1, 2), M > 0 such that

B(θ) p ≤ M Γ(θ), θ 2 B ′ (θ) ≤ M Γ(θ), zJ ′ (z) ≤ M J(z), J ′ (z) p ′ ≤ M J(z), 1/p + 1/p ′ = 1.
(H3) The initial data f 0 ≥ 0 and θ 0 ≥ 0 are such that the following quantities

R 3 θ 0 dx, R 3 |θ 0 | 2 dx, R 3 Γ(θ 0 ) dx R 3 S 2 f 0 dv dx, R 3 S 2 |f 0 | 2 dv dx, R 3 S 2 J(f 0 ) dv dx
are bounded uniformly with respect to C, P, L a , L s .

Remark 9 Restricting to the special case [START_REF] Dogbe | The radiative transfer equations: diffusion approximation under accretiveness and compactness assumptions[END_REF], the functions of the temperature involved in these relations are defined as follows

B(θ) = 1 π θ 4 , Γ(θ) = 1 5π θ 5 , G(θ) = 2 5π θ 5/2 .
Hence, assumption (H2) holds in this case (with p = 5/4).

The crucial properties of the complete system (20) can be summarized in the following claim (recall that we assume D = 1).

Proposition 2 Let (f, θ) be a solution of [START_REF] Xin | The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF][START_REF] Klar | Numerical Passage from Radiative Heat Transfer to Nonlinear Diffusion Models[END_REF]. We assume that (H0) and (H1) hold. We also assume that C ∈ (C * , +∞), with a large enough C * > 0 depending on (H0). Let us consider two nonnegative and C 2 convex functions Φ, Ψ : R + → R + . We assume that there exists a constant M such that for any s ≥ 0, we have |Φ ′ (s)|s ≤ M Φ(s). We set

H(f, θ)(t) = R 3 S 2 Φ(Λ 4 f ) 1 Λ 3 dv dx + C P R 3 Ψ(θ) dx ≥ 0. ( 23
)
Then, we have

d dt H(f, θ) + C P R 3 Ψ"(θ)|∇ x θ| 2 + D s + D a = C 1 + C 2 + R 1 + R 2 , (24)
with

D s = CL s R 3 S 2 S 2 σ s (v, v ′ ) 2Λ(v) 2 Λ(v ′ ) 2 Λ(v) 4 f (v) -Λ(v ′ ) 4 f (v ′ ) × Φ ′ (Λ 4 f )(v) -Φ ′ (Λ 4 f )(v ′ ) dv ′ dv dx ≥ 0, D a = CL a R 3 S 2 σ a (v) γΛ(v) 2 Λ(v) 4 f (v) -B(θ) × Φ ′ (Λ 4 f )(v) -Ψ ′ (θ) dv dx.
The right-hand side terms in ( 24) can be estimated as follows

|C 1 | ≤ K R 3 S 2 Φ(Λ 4 f )(v) 1 Λ 3 dv dx, (25)
|C 2 | ≤ C 2P R 3 Ψ"(θ)|∇ x θ| 2 dx + C u 2 L ∞ 2P R 3 θ 2 Ψ"(θ) dx, (26)
|R 1 | ≤ K C R 3 S 2 Φ(Λ 4 f )(v) 1 Λ 3 dv dx, (27)
|R 2 | ≤ K L a C R 3 S 2 |Λ 4 f -B(θ)| |Φ ′ (Λ 4 f )| Λ 2 dv dx, ( 28
)
where the constant K > 0 does not depend on the parameters C, P, L s and L a .

Remark 10

The term D a becomes nonnegative provided that Φ and Ψ can be linked by the relation

Φ ′ (B(θ)) = Ψ ′ (θ) or equivalently Φ ′ (z) = Ψ ′ (B (-1) (z)).
Note that this implies, by convexity, that Φ"(B(θ))B ′ (θ) = Ψ"(θ) ≥ 0. In such a case, (23) and (24) can be seen as generalized entropy and entropy dissipation relations. For the black body law [START_REF] Dogbe | The radiative transfer equations: diffusion approximation under accretiveness and compactness assumptions[END_REF], the physical choice Ψ(θ) = ln(θ) leads to Φ(s) = 4 3 π -1/4 s 3/4 (coming back to a model with frequencies we would obtain, up to some constants, the usual entropy function Φ(s) = (1 + s) ln(1 + s) -s ln(s), see [START_REF] Ch | Asymptotic analysis of fluid models for the coupling of radiation and hydrodynamics[END_REF]).

We postpone the proof to Appendix C. In consequence, we are able to derive the following useful a priori estimates, whose proofs will also be found in Appendix D and E.

Corollary 1 Assume that (H0-3) are fulfilled and consider C ∈ (C * , +∞) as in Proposition 2. We assume that there exists a constant M > 0 such that

P C , L a C , 1 CL a ∈ (0, M ]. ( 29
)
Let T ∈ (0, ∞). Then, the following quantities

R 3 θ 2 dx, R 3 Γ(θ) dx (30) R 3 S 2 f 2 dv dx, R 3 S 2 J(f ) dv dx (31) t 0 R 3 |∇ x θ| 2 dx ds, t 0 R 3 |∇ x G(θ)| 2 dx ds (32) CL s t 0 R 3 S 2 |Λ 4 f (v ′ ) -Λ 4 f (v)| 2 dv ′ dv dx ds, ( 33
)
CL a t 0 R 3 S 2 |Λ 4 f -B(θ)| 2 dv dx ds, ( 34
)
are bounded uniformly with respect to C, P, L a , L s and t ∈ [0, T ].

Remark 11 Of course, the restriction ( 29) is fulfilled for the asymptotic regimes we are interested in. Precisely, we always assume P/C = 1 and either

CL a = 1/ε 2 ≫ 1, L a /C = 1 (equilibrium regime) or CL a = 1, L a /C ≪ 1 (nonequilibrium regime).
Corollary 2 Assume that (H0-3) are fulfilled and consider C ∈ (C * , +∞) as in Proposition 2. Then, the quantity

R 3 S 2 f dv dx + C P R 3 θ dx
as well as B(θ) dx are bounded uniformly with respect to C, P, L a ,L s and t ∈ [0, T ].

Statements of the Results

We can now write the statements of our main results.

Theorem 1 We consider the non-equilibrium regime C = 1/ε, P = 1/ε, L s = 1/ε, L a = ε. Assume that (H0-3) are fulfilled. Then, up to a subsequence, ρ ε = f ε and θ ε converge to ρ and θ respectively, strongly in L 2 loc ((0, T ) × R 3 ) and in C 0 ([0, T ]; L 2 (R 3 ) -weak), f ε converges to ρ weakly in L 2 ((0, T ) × R 3 × S 2 ) and the limit satisfies the following system of drift-diffusion equations1

∂ t ρ -div x (D∇ x ρ -3 v ⊗ v uρ) = -(D x u) T v • v ρ + σ a (B(θ) -ρ), (35)
∂ t θ + div x (uθ) -∆θ = σ a (ρ -B(θ)). ( 36
)
These equations involve the positive diffusion coefficient

D(x) = v ⊗ χ > 0 with χ ∈ L ∞ (R 3 ; L 2 (S 2 )) 3 that solves S 2 σ s (x, v, v ′ ) χ(x, v ′ ) -χ(x, v) dv ′ = -v
(we refer to Proposition 3 below for a precise statement on the definition of χ). The system is completed by the initial data ρ |t=0 = lim ε→0 f ε,0 and θ |t=0 = lim ε→0 θ ε,0 , the limits being understood as weakly in L 2 (R 3 ).

Theorem 2 We consider the equilibrium regime C = 1/ε, P = 1/ε, L s = ε, L a = 1/ε. Assume that (H0-3) are fulfilled. Then, up to a subsequence, ρ ε = f ε and θ ε converge to B(θ) and θ respectively, strongly in

L 2 loc ((0, T ) × R 3 ), (ρ ε + θ ε ) converges to B(θ) + θ in C 0 ([0, T ]; L 2 (R 3 ) -weak)
, and f ε converges to B(θ) weakly in L 2 ((0, T ) × R 3 × S 2 ). The limit satisfies the following non linear drift-diffusion equation

∂ t θ + B(θ) + div x u θ + 3 v ⊗ v B(θ) -∇ x θ - v ⊗ v σ a ∇ x B(θ) = -(D x u) T v • v B(θ).
The system is completed by the initial data (B(θ) + θ) |t=0 = lim ε→0 ( f ε,0 + θ ε,0 ), the limit being understood weakly in L 2 (R 3 ).

Remark 12 Since we are able to prove strong compactness properties, it is worth pointing out that our proofs also apply to more non linear problems. Namely, we can consider coefficients σ s,a which depend continuously on ρ or θ, without requiring any monotonicity property.

Remark 13 Note that the limit equations, in both cases, contain drift terms depending on u, as well as zero-th order terms. In particular, the limit equations are not conservative.

5 Non-equilibrium Regime: Proof of Theorem 1

Theorem 1 deals with the scaling [START_REF] Marchuk | Splitting and alternating direction methods in Handbook of numerical analysis[END_REF] and

L s = 1 ε , L a = ε.
Accordingly, we are studying the behavior as ε tends to 0, of the sequence (f ε , θ ε ) of solutions of the following system

                         ε∂ t f ε + v • ∇ x f ε = 1 ε Q s,ε (f ε ) + εQ a,ε (f ε , θ ε ), Q s,ε (f ε ) = σ s Λ 2 ε f ε Λ 3 ε - σ s Λ 2 ε Λ ε f ε , Q a,ε (f ε , θ ε ) = σ a B(θ ε ) Λ 3 ε -Λ ε f ε , ∂ t θ ε + ∇ x • (uθ ε ) -∆θ ε = - 1 ε Λ ε γ ε 1 ε Q s,ε (f ε ) + εQ a,ε (f ε , θ ε ) = - Λ ε γ ε Q a,ε (f ε , θ ε ) , (37)
where we used the notation

Λ ε = γ ε (1 -εu • v), γ ε = 1/ √ 1 -ε 2 u 2
(note also that we used Proposition 1-iii)). It is worth starting with the following preliminary lemma (see (65) and (66)).

Lemma 1 Assume (H0) is satisfied. There exist constants λ * , λ * and K which do not depend on ε ∈ (0, ε * )

such that 0 ≤ λ * ≤ Λ ε (t, x, v) ≤ λ * , |∂ t Λ ε (t, x, v)| ≤ Kε hold for every (t, x, v) ∈ (0, T ) × R 3 × S 2 .
Moreover, the sequence Λ ε satisfies the following properties when

ε tends to 0 Λ ε → 1, 1 -Λ α ε ε → αu • v (α ∈ R), ∇ x Λ ε ε → -(D x u) T v, uniformly on (0, T ) × R 3 × S 2 .
Let us define the following macroscopic quantities

ρ ε (t, x) = f ε , J ε (t, x) = v ε f ε .
Since the scattering operator is penalized, we guess that f ε behaves for small ε as an element of the kernel of the operator. Recalling Proposition 1, we thus set

f ε (t, x, v) = Λ 4 ε f ε Λ -4 ε + εg ε . (38)
Now, the proof of Theorem 1 splits into four steps.

A priori Estimates

First of all, let us discuss the uniform estimates satisfied by the sequence of solutions (f ε , θ ε ). They can be deduced as consequences of Corollaries 1 and 2.

Lemma 2

The following properties are satisfied i)

f ε ε>0 is bounded in L ∞ (0, T ; L 1 ∩ L 2 (R 3 × S 2 )), ii) ρ ε ε>0 is bounded in L ∞ (0, T ; L 1 ∩ L 2 (R 3 )), iii) g ε ε>0 is bounded in L 2 ((0, T ) × R 3 × S 2 ), iv) J ε ε>0 is bounded in L 2 ((0, T ) × R 3 ) 3 , v) θ ε ε>0 is bounded in L ∞ (0, T ; L 1 ∩ L 2 (R 3 )) and in L 2 (0, T ; H 1 (R 3 )), vi) B(θ ε ) -Λ 4 ε f ε ε>0 is bounded in L 2 ((0, T ) × R 3 × S 2 ).
Remark 14 This claim justifies the expansion of f ε since εg ε is of order O(ε). Note also that, coming back to a physical quantity, it is natural to interpret (ρ/σ) 1/4 as the radiation temperature.

Proof. Properties i), iii), v) and vi) are direct consequences of ( 31), ( 33), ( 30), ( 32) and (34) respectively. Estimate ii) follows from [START_REF] Tartar | Remarks on homogenization in Homogenization and effective moduli of material and media[END_REF] and a simple application of the Cauchy-Schwarz inequality. A similar argument shows that Λ 4

ε f ε is bounded in L ∞ (0, T ; L 2 (R 3
)). Finally, we note that

J ε = v ε Λ 4 ε f ε Λ -4 ε + εg ε = Λ 4 ε f ε v ε (Λ -4 ε -1) + v ε + vg ε = Λ 4 ε f ε v ε (Λ -4 ε -1) + 0 + vg ε , (39)
since v = 0. We conclude by using Lemma 1 combined to the L 2 estimate on g ε .

Moment Equations

Now, let us turn to the moment equations satisfied by the moments of f ε . We split the result in the following two Lemmata.

Lemma 3

We have

∂ t ρ ε + div x J ε - (1 -Λ ε ) v ε f ε = ε(∂ t S ε + V ε ) + U ε where S ε and V ε are bounded in L ∞ (0, T ; L 2 (R 3 )) and U ε = v ε •∇ x Λ ε f ε + Λ ε Q a,ε is bounded in L 2 ((0, T )× R 3 ).
Proof. Let us make the following computation

∂ t ρ ε + div x J ε = ∂ t Λ ε f ε + div x v ε Λ ε f ε +∂ t (1 -Λ ε ) f ε + div x v ε (1 -Λ ε ) f ε .
However, we check that

∂ t Λ ε f ε = ∂ t Λ ε f ε -div x v ε Λ ε f ε + v ε • ∇ x Λ ε f ε + Λ ε Q a,ε ,
since Λ ε Q s,ε = 0. Therefore, we obtain the announced formula with

S ε = 1 ε (1 -Λ ε ) f ε , V ε = 1 ε ∂ t Λ ε f ε which are both bounded in L ∞ (0, T ; L 2 (R 3
)) by virtue of Lemma 1 and 2-i). On the other hand, by using Lemma 1 and 2-i) again, we realize that v

• ∇xΛε ε f ε is also bounded in L ∞ (0, T ; L 2 (R 3 )). Eventually, writing B(θ) Λ 3 -Λf = 1 Λ 3 (B(θ) -Λ 4 f ), we see that Λ ε Q a,ε is bounded in L 2 ((0, T ) × R 3 × S 2
) by using Lemma 2 vi).

Lemma 4 We have ε 2 ∂ t J ε + Div x P ε = Q ε where the right-hand side Q ε = v ε Q s,ε (f ε ) + εvQ a,ε (f ε ) is bounded in L 2 ((0, T ) × R 3 )
3 and

P ε = v ⊗ vf ε = v ⊗ v ρ ε + εK ε K ε being bounded in L 2 ((0, T ) × R 3 ) 3×3 .
Proof. Multiplying Equation (37) by v and integrating leads immediately to the formula. What remains to be discussed are the uniform bounds. At first, we simply expand

v ⊗ vf ε = v ⊗ v ρ ε + v ⊗ v(f ε -ρ ε ) .
Then, we note that f ε -

ρ ε = εg ε + (Λ 4 ε -1)f ε Λ -4 ε + ρ ε (Λ -4 ε -1
) which proves that the components of the matrix

K ε = 1 ε v ⊗ v(f ε -ρ ε ) are bounded in L 2 ((0, T ) × R 3 ). Next, note that 1 ε Q s,ε (f ε ) = Λ 4 ε f ε ε Q s,ε (Λ -4 ε ) + Q s,ε (g ε ) = Q s,ε (g ε ).
It is bounded in L 2 ((0, T ) × R 3 ) since we have

S 2 |Q s,ε (g ε )| 2 dv ≤ K S 2 | g ε | 2 dv + S 2 g 2 ε dv ≤ K S 2 g 2 ε dv
where the constant K depends only on (H1) and the bounds in Lemma 2-iii). This ends the proof since we already proved that Q a,ε is bounded in L 2 ((0, T ) × R 3 ).

Compactness Properties

These relations allow us to deduce strong compactness properties on the sequence ρ ε ε>0 . To this end we use a compensated compactness argument which relies on the structure of the moment equations. This has been remarked first by Marcati-Milani [START_REF] Marcati | The one-dimensional Darcy's law as the limit of a compressible Euler flow[END_REF], and then used in various contexts when dealing with diffusion approximations [START_REF] Lions | Diffuse limit for finite velocity Boltzmann kinetic models[END_REF], [START_REF] Degond | Diffusion limit for non homogeneous and non reversible processes[END_REF], [START_REF] Goudon | Approximation by homogeneization and diffusion of kinetic equations[END_REF]...

Lemma 5

The sequence ρ ε ε>0 lies in a (strong) compact set of

L 2 loc ((0, T ) × R 3 ). It is also compact in C 0 ([0, T ]; L 2 (R 3 ) -weak).
Proof. The second part of the statement means that we can extract a subsequence, still labelled ρ ε ε>0 , such that for any test function ϕ ∈ L 2 (R 3 ), ρ ε ϕ dx converges to ρϕ dx in C 0 ([0, T ]) for some limit function ρ ∈ L ∞ (0, T ; L 2 (R 3 )). This follows from the bound on ρ ε (Lemma 2-ii)) combined to Lemma 3 which tells us that

∂ t (ρ ε + εS ε ) is bounded in L 2 (0, T ; H -1 (R 3 )), with S ε bounded in L ∞ (0, T ; L 2 (R 3
)). This property allows us to recover the initial data when passing to the limit ε → 0.

To prove the strong L 2 -compactness, we use the celebrated Div-Curl Lemma of Murat-Tartar [START_REF] Tartar | Remarks on homogenization in Homogenization and effective moduli of material and media[END_REF]. Indeed, Lemma 3 tells us that

div t,x ρ ε , J ε - (1 -Λ ε ) v ε f ε , ε > 0
lies in a compact set of H -1 loc ((0, T )×R 3 ) by using the Rellich embedding theorem. Of course, the components (ρ ε ) and (J ε -(1-Λε) v ε f ε ) of this vector field are bounded in L 2 ((0, T ) × R 3 ). Next, we rewrite Lemma 4 as

∇ x ρ ε = S 2 v ⊗ v dv -1 -ε 2 ∂ t J ε -εDiv x K ε + Q ε , (40)
which implies that the curl of the R 4 vector field (ρ ε , 0, 0, 0) belongs to a compact set of H -1 loc ((0, T )×R 3 )

4×4

since the matrix v ⊗ v dv is invertible (note the loss of one derivative with respect to time). Assume we are dealing with weakly convergent sequences. Then, the Div-Curl lemma implies that the product (ρ ε , J ε -(1-Λε) v ε f ε ) • (ρ ε , 0) passes to the limit in the D ′ ((0, T ) × R 3 ) sense, which in turn proves that ρ ε ) ε>0 converges strongly in L 2 loc ((0, T ) × R 3 ).

Remark 15

The key argument relies on the invertibility of the matrix v ⊗ v dv. It is worth pointing out that this argument allows us to consider discrete velocity models, which is not the case if we use arguments based on average lemma techniques (see [START_REF] Golse | Regularity of the moments of the solution of a transport equation[END_REF], [START_REF] Perna | L p regularity of velocity averages[END_REF]). Indeed, our results apply considering the variable v in some measured space (V, dv),

V ⊂ R 3 , such that:    V is bounded, 1 1 = 1, v = 0, meas{v ∈ V, v • ξ = 0} > 0,
for any ξ ∈ S 2 .

We refer to [START_REF] Goudon | Approximation by homogeneization and diffusion of kinetic equations[END_REF] and [START_REF] Degond | Diffusion limit for non homogeneous and non reversible processes[END_REF] for further comments and applications on this aspect.

Now, let us go back to the material temperature θ ε .

Lemma 6

The sequence θ ε ε>0 lies in a (strong) compact set of

L 2 loc ((0, T ) × R 3 ). It is also compact in C 0 ([0, T ]; L 2 (R 3 ) -weak).
Proof. The result follows from the relation

∂ t θ ε = -div x (uθ ε ) + ∆ x θ ε - Λ ε γ ε Q a,ε ,
where the right-hand side is bounded in L 2 (0, T ; H -1 (R 3 )). Hence, the strong compactness follows from an application of the Aubin lemma, see e.g. [START_REF] Simon | Compact sets in the space L p (0, T ; B)[END_REF].

Passage to the limit

We are now ready to pass to the limit in Equation (37). By Lemma 2, we can assume, possibly at the cost of extracting a subsequence, that

               f ε ⇀ f weakly in L 2 ((0, T ) × R 3 × S 2 ), ρ ε → ρ = S 2 f dv strongly in L 2 ((0, T ) × R 3 ) and in C 0 ([0, T ]; L 2 (R 3 ) -weak), J ε ⇀ J weakly in L 2 ((0, T ) × R 3 ) 3 , θ ε → θ strongly in L 2 ((0, T ) × R 3 ) and in C 0 ([0, T ]; L 2 (R 3 ) -weak).
Actually, the limit f is the macroscopic quantity ρ since we have

f ε = Λ 4 ε f ε Λ -4 ε + εg ε ⇀ f = f + 0 = ρ weakly in L 2 ((0, T ) × R 3 × S 2 ).
Furthermore, we can also assume that (ρ ε ) and (θ ε ) converge almost everywhere and are dominated. Then, by using Lemma 2-i) and vi), we realize that

B(θ ε ) = Λ 4 ε f ε + B(θ ε ) -Λ 4 ε f ε is bounded in L 2 ((0, T ) × R 3 )
. By using classical tricks of integration theory we can thus show that B(θ ε ) → B(θ) strongly in L 1 loc ((0, T ) × R 3 ) and a.e. (

(the convergence holds at least in L p loc ((0, T ) × R 3 ), 1 ≤ p < 2 and can certainly be improved up to some assumptions on the behavior of the function B). Therefore, passing to the limit in the zero-th moment equation yields

∂ t ρ + div x (J -v ⊗ v uρ) = -v (D x u) T v ρ + σ a (B(θ) -ρ) (42)
by using Lemma 1.

Besides, let us remark that the well-known regularizing effects of diffusion asymptotics also apply in this context: the macroscopic limit ρ has better regularity properties than ρ ε itself.

Lemma 7 The limit ρ of ρ ε satisfies ∇ x ρ ∈ L 2 ((0, T ) × R 3 ). Proof. For any test function ϕ ∈ C ∞ c ((0, T ) × R 3 , (40) leads to ∞ 0 R 3 ρ∇ x ϕ dx dt = lim ε→0 ∞ 0 R 3 ρ ε ∇ x ϕ dx dt = lim ε→0 ∞ 0 R 3 ε 2 v ⊗ v -1 J ε ∂ t ϕ + ε v ⊗ v -1 K ε ∇ x ϕ dx dt + ∞ 0 R 3 v ⊗ v -1 Q ε ϕ dx dt ≤ sup ε>0 Q ε L 2 ϕ L 2 ,
since all derivated terms vanish as

ε → 0. Since (Q ε ) is bounded in L 2 ((0, T ) × R 3 ), we conclude that ∞ 0 R 3 ρ∇ x ϕ dx dt ≤ K ϕ L 2 , which thus implies that ∇ x ρ ∈ L 2 ((0, T ) × R 3 ).
Next, we immediately pass to the limit in the equation for the material temperature and we get

∂ t θ + div x (uθ) -∆θ = σ a (ρ -B(θ)).
We are thus left with the task of identifying the limit J of J ε . We remark that J is related to a weighted average of g, the (weak) limit of g ε . Indeed, letting ε → 0 in (39) yields J = 4 v u • v ρ + vg , by using Lemma 1. Therefore, our aim becomes identifying g. To this end, we need the following claim.

Proposition 3 Let us set Q s,0 (f ) = S 2 σ s (v, v ′ ) f (v ′ ) -f (v) dv ′ .
It defines a self-adjoint bounded operator on L 2 (S 2 ). The kernel of this operator is Ker(Q s,0 ) = Span{1 1} and for any ψ ∈ L 2 (S 2 ) such that ψ dv = 0, there exists a unique f ∈ L 2 (S 2 ) satisfying

Q s,0 (f ) = ψ, S 2 f dv = 0.
Taking into account the dependence with respect to x of the collision kernel, if ψ belongs to L ∞ (R 3 ; L 2 (S 2 )) then f does too.

The proof relies on the application of the Fredholm alternative for the operator Q s,0 . It will be detailed in Appendix F. In particular, this result allows us to define χ ∈ L ∞ (R 3 ; L 2 (S 2 ))

3 as the solution (with vanishing integral) of Q s,0 (χ) = -v, as in Theorem 1.

Let ϕ ∈ C ∞ c (R 3 × S 2 ). We have

ε d dt R 3 S 2 f ε ϕ dv dx - R 3 S 2 vf ε • ∇ x ϕ dv dx = 1 ε R 3 S 2 Q s,ε (f ε ) ϕ dv dx + ε R 3 S 2 Q a,ε (f ε , θ ε ) ϕ dv dx.
We rewrite

1 ε Q s,ε (f ε ) = Q s,ε (g ε ) so that we obtain, as ε → 0, - R 3 S 2 vρ • ∇ x ϕ dv dx = R 3 S 2 Q s,0 (g) ϕ dv dx.
This is also a simple consequence of Lemma 1 combined with the weak convergence

g ε ⇀ g in L 2 ((0, T ) × R 3 × S 2 ). It means that g is solution of Q s,0 (g) = v • ∇ x ρ, which belongs to L 2 (R 3 × S 2 ), as noted in Lemma 7. It follows that g = -χ • ∇ x ρ
(up to an element in Ker(Q s,0 ) which is irrelevant in the following results). We conclude that

J(t, x) = 4 v u • v ρ -v ⊗ χ ∇ x ρ. (43)
Inserting this relation into (42) ends the proof of Theorem 1.

Remark 16

The diffusion coefficient is positive. This is a consequence of Proposition 1 (applied with u = 0, i.e. Λ = 1). Indeed, for any ξ ∈ R N \{0}, we have

v ⊗ χ ξ • ξ = -Q s,0 (χ • ξ) χ • ξ ≥ σ * χ • ξ 2 > 0.
The right-hand side cannot vanish, since otherwise Q s,0 (χ • ξ) = -v • ξ would vanish a.e.

Remark 17 Considering the isotropic case [START_REF] Perna | L p regularity of velocity averages[END_REF], we simply have χ = v σs and J = 4 3 ρu -1 3σs ∇ x ρ. We are thus led to [START_REF] Battaner | Astrophysical fluid dynamics[END_REF].

Equilibrium Regime: Proof of Theorem 2

Theorem 2 deals with the scaling [START_REF] Marchuk | Splitting and alternating direction methods in Handbook of numerical analysis[END_REF] and

L s = ε, L a = 1 ε .
Accordingly, we are studying the behavior as ε tends to 0, of the sequence (f ε , θ ε ) of solutions of the following system

         ε∂ t f ε + v • ∇ x f ε = εQ s,ε (f ε ) + 1 ε Q a,ε (f ε , θ ε ), ∂ t θ ε + ∇ x • (uθ ε ) -∆θ ε = - 1 ε Λ ε γ ε εQ s,ε (f ε ) + 1 ε Q a,ε (f ε , θ ε ) = - 1 ε 2 Λ ε γ ε Q a,ε (f ε , θ ε ) , (44)
where all notations are left unchanged with respect to those of the previous section. The proof follows the same scheme as the proof of Theorem 1; hence we skip some details and only emphasize the main changes in the arguments.

Step

A priori Estimates

We start with the following direct consequence of Corollaries 1 and 2:

-

f ε ε>0 is bounded in L ∞ (0, T ; L 1 ∩ L 2 (R 3 × S 2 )), -ρ ε ε>0 is bounded in L ∞ (0, T ; L 1 ∩ L 2 (R 3 )), -θ ε ε>0 is bounded in L ∞ (0, T ; L 1 ∩ L 2 (R 3 )) and in L 2 (0, T ; H 1 (R 3 )); -g ε = fε-Λ -4 ε B(θε) ε ε>0 is bounded in L 2 ((0, T ) × R 3 × S 2
). Furthermore, we remark that

ρ ε -B(θ ε ) L 2 ((0,T )×R 3 ) = (1 -Λ 4 ε )f ε + ε Λ 4 ε g ε L 2 ((0,T )×R 3 ) ≤ Kε (45)
by using Lemma 1. Finally, let us set J ε = v ε f ε . We rewrite it as

J ε = vg ε + v Λ -4 ε -1 ε B(θ ε ) (46) which shows that (J ε ) is bounded in L 2 ((0, T ) × R 3 )
3 , by using Lemma 1 again.

Step 2. Moment Equations Combining the zero-th moment equation and the temperature equation leads to

∂ t (ρ ε + θ ε ) + div x J ε -v 1 -Λ ε ε f ε + uθ ε -∇ x θ ε = ε(∂ t S ε + V ε ) + U ε ( 47
)
where S ε and V ε are defined as in Lemma 3 and are bounded in L ∞ (0, T ; L 2 (R 3 )) while

U ε = v ε • ∇ x Λ ε f ε + Λ ε 1 -1/γ ε ε 2 Q a,ε (f ε , θ ε ) . Obviously, v ε ∇ x Λ ε f ε is bounded in L 2 ((0, T ) × R 3
). Then, we note that

Q a,ε (f ε , θ ε ) = -εσ a Λ ε g ε ---→ ε→0 0 in L 2 ((0, T ) × R 3 × S 2 ) while we check readily that 1-1/γε ε 2
is bounded in L ∞ ((0, T ) × R 3 ). Hence, the last term in the definition of U ε vanishes as ε → 0.

Besides, we get

ε 2 ∂ t J ε + Div x P ε = 1 ε vQ a,ε + ε vQ s,ε (f ε ) = -vσ a Λ ε g ε + ε vQ s,ε . (48)
The right-hand side is bounded in L 2 ((0, T ) × R 3 × S 2 ) 3 . As a consequence of (45) we obtain

P ε (t, x) = v ⊗ v ρ ε + ε K ε , (49)
with

K ε = -v⊗v ε (ρ ε -B(θ ε )) -εΛ 4 ε g ε + (Λ 4 ε -1)f ε bounded in L 2 ((0, T ) × R 3 ) 3 .
Step 3. Compactness Properties Extracting subsequences if necessary, we can assume that -

f ε ⇀ f weakly in L 2 ((0, T ) × R 3 × S 2 ), -g ε ⇀ g weakly in L 2 ((0, T ) × R 3 × S 2 ), -ρ ε ⇀ ρ weakly in L 2 ((0, T ) × R 3 ), -J ε ⇀ J weakly in L 2 ((0, T ) × R 3 ) 3 , -θ ε ⇀ θ weakly in L 2 ((0, T ) × R 3 ). Since f ε = B(θ ε )Λ -4
ε + εg ε , we deduce from (45) that f ε ⇀ f = ρ = f which does not depend on the variable v and coincides with the weak limit in L 2 of the sequence B(θ ε ). Thus, we might expect that ρ = B(θ), but the justification of this relation requires some strong compactness. To this end, we use a compensated compactness argument again. Indeed, from (47) we deduce that

div t,x ρ ε + θ ε , J ε -v 1 -Λ ε ε f ε + uθ ε -∇ x θ ε , ε > 0 lies in a compact set of H -1 loc ((0, T ) × R 3
). Note that this relation also proves that (ρ ε + θ ε ) lies in a compact set of C 0 ([0, T ]; L 2 (R 3 ) -weak).

On the other hand, we rewrite ( 48) and (49) as

∇ x ρ ε = v ⊗ v -1 -ε 2 ∂ t J ε -εDiv x K ε -v(σ a Λ ε g ε -εQ s,ε (f ε ))
which implies that:

-∇ x ρ ε ε>0 belongs to a compact set of H -1 loc ((0, T ) × R 3 ), -the gradient with respect to x of the limit ρ belongs to L 2 ((0, T ) × R 3 ). Furthermore, the a priori estimates guarantee that (∇ x θ ε ) is bounded in L 2 ((0, T ) × R 3 ), thus it also lies in a compact set of H -1 loc ((0, T ) × R 3 ). We interpret these facts by saying that curl t,x (ρ ε + θ ε , 0, 0, 0) belongs to a compact set of H -1 loc ((0, T ) × R 3 ) 4×4 . Therefore, applying the Div-Curl lemma yields

(ρ ε + θ ε ) 2 ⇀ (ρ + θ) 2 , in D ′ ((0, T ) × R 3 ). We deduce that ρ ε + θ ε converges strongly to ρ + θ in L 2 loc ((0, T ) × R 3 ). Let us set F (z) = z + B(z). Using (45), we realize that F (θ ε ) = ρ ε + θ ε + (B(θ ε ) -ρ ε ) converges strongly to ρ + θ in L 2 loc ((0, T ) × R 3
). Extracting more subsequences if necessary, we can assume that this convergence holds a.e. so that θ ε = F (-1) (F (θ ε )) tends to F (-1) (ρ + θ) a.e. Consequently, θ ε converges to θ = F (-1) (ρ + θ) a.e. and strongly in L p loc ((0, T ) × R 3 ), 1 ≤ p < 2. Coming back to (45), it follows that ρ ε converges to F (θ) -θ = B(θ) a.e. and strongly in L p loc ((0, T ) × R 3 ), 1 ≤ p < 2.

Step 4. Passage to the Limit Letting ε tend to 0 in (47) yields

∂ t (ρ + θ) + div x J -(v ⊗ v)u ρ + uθ -∇ x θ = -v D x u T v ρ.
We have seen that ρ = B(θ) and it only remains to identify J. Let ϕ ∈ C ∞ c (R 3 × S 2 ). We have

ε d dt R 3 S 2 f ε ϕ dv dx - R 3 S 2 vf ε • ∇ x ϕ dv dx = 1 ε R 3 S 2 Q a,ε (f ε , θ ε ) ϕ dv dx + ε R 3 S 2 Q s,ε (f ε ) ϕ dv dx. We rewrite 1 ε Q a,ε (f ε , θ ε ) = -σ a Λ ε g ε . Letting ε → 0 leads to - R 3 S 2 vρ • ∇ x ϕ dv dx = - R 3 S 2
σ a gϕ dv dx.

Since this relation holds for any test function and

∇ x ρ ∈ L 2 ((0, T ) × R 3 ), we deduce that v • ∇ x ρ = -σ a g
holds. Hence, coming back to (46), we obtain the formula

J(t, x) = - v ⊗ v σ a ∇ x ρ + 4 v ⊗ v uρ,
which ends the proof of Theorem 2.

Numerical Results in One Dimension

In order to be able to compute solutions of (37) numerically, we need to simplify the expressions of the source-terms so that we can recognize in the kinetic equation the non-stiff and the stiff parts that allow the use of classical splitting schemes [START_REF] Strang | Approximating semigroups and the consistency of difference schemes[END_REF], [START_REF] Marchuk | Splitting and alternating direction methods in Handbook of numerical analysis[END_REF], [START_REF] Descombes | Operator splitting for nonlinear reaction-diffusion systems with an entropic structure: singular perturbation and order reduction[END_REF]. These methods have been well known since decades and are particularly efficient for evolution equations; they have been used very successfully to treat relaxation of hyperbolic problems, see e.g. [START_REF] Xin | The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF][START_REF] Aregba-Driollet | Convergence of relaxation schemes for conservation laws[END_REF]..., the survey [START_REF] Natalini | Recent results on hyperbolic relaxation problems in Analysis of systems of conservation laws[END_REF] and the references therein. Let us consider the non-equilibrium regime described in equations (37) in the one-dimensional case, that is, the space variable x belongs to R and the direction variable v belongs to (-1, 1). For numerical simulations of the equilibrium regime, we refer to [START_REF] Klar | Numerical Passage from Radiative Heat Transfer to Nonlinear Diffusion Models[END_REF]. In what follows we restrict to the situations with constant coefficients σ a,s and to the simple coupling with the convection-diffusion equation for the material temperature. Further numerical investigations with the full hydrodynamic system are postponed to a forthcoming work. Let us also mention the recent work [START_REF] Buet | Asymptotic preserving scheme and numerical methods for radiative hydrodynamic models[END_REF] which uses a different numerical approach to treat a macroscopic version of the problem. The aims of the numerical investigation are: -to check on numerical grounds the convergence as ε → 0 to the solutions of the limit equations, -to discuss the role of the assumptions, in particular the regularity of the velocity field u and the heat diffusion D.

Expansion of the equations

At first, let us expand equations (37) with respect to ε. Expanding the kinetic unknown as

f ε (t, x, v) = ρ ε (t, x) + εg ε (t, x, v), we rewrite Q s,ε as 1 σ s Q s,ε (f ε ) = -εg ε + T 1 (v)ρ ε + T 2 ρ ε + εS 1 (v, gε ), (50)
where

T 1 (v) = Λ 2 Λ -3 -Λ 2 Λ -3 -(Λ -Λ ), T 2 = Λ 2 Λ -3 -Λ , S 1 (v, gε ) = Λ 2 gε Λ -3 -(Λ -1)g ε ,
keeping in mind that we will get rid of the higher oder terms in (37) (see Remark 18). We get

• T 1 (v) = 0 and T 1 (v) = 4εuv + 2ε 2 u 2 (3v 2 -1) + O(ε 3 ), • T 2 = 4 3 ε 2 u 2 + O(ε 4 ), • S 1 (v, gε ) = εu(vg ε -2 vg ε ) + O(ε 2 ).
In the same way, we have 1

σ a Q a,ε (f ε , θ ε ) = B(θ ε ) -ρ ε + O(ε). (51)
Let us now rewrite (37) according to the previous Taylor expansions, adding -v∂ x ρ on the left-and righthand sides of the kinetic equation to introduce g, so that

∂ t f ε + v∂ x gε = -σ s gε -v∂ x ρ ε + σ s 4vuρ ε ε + σ s u 2 3v 2 - 1 3 ρ ε u + vg ε -2 vg ε +σ a (B(θ ε ) -ρ ε ) + O(ε), ∂ t θ ε + ∂ x (uθ ε ) -D∂ 2 xx θ ε = -σ a (B(θ ε ) -ρ ε ) + O(ε).
(52)

Numerical scheme

We need to describe precisely the discrete setting we are going to use : let the direction variables (v k ) k∈{1,...,2p} belong to (-1, 1) \ {0} for p ∈ N * and satisfy

• 1 2p 2p k=1 v k = 0, • 1 2p 2p k=1 v 2 k = 1 3 ,
in accordance with the choice of a Lebesgue measure in the continuous setting. Let us denote the discretized quantities by f n ε,k,j (resp. θ n ε,j , u n j ) where j ∈ Z and n ∈ N, which represents an approximate value of f ε (n∆t, j∆x, v k ) (resp. θ ε (n∆t, j∆x), u(n∆t, j∆x)) where ∆t and ∆x are the time and space steps : f is now the double sequence of vectors f ε = (f n ε,k,j ) 1≤k≤2p,n∈N,j∈Z . For the sake of simplicity, we will often omit the subscripts k and j . In order to compute numerically the solutions of (37), we will use a splitting-type scheme for the intensity f . Indeed, the kinetic equations for f contain source terms that can be split into a non-stiff part and a stiff part, the mean value of which vanishes, so that the second part of the splitting scheme can be solved explicitly up to a spatial discretization. This is close to the approach in [START_REF] Bouchut | Diffusive BGK approximations for nonlinear multidimensional parabolic equations[END_REF]. We set up the numerical resolution of (52) as follows:

• assume the initial data f 0 (respectively

θ 0 ) is discretized in (f 0 j ) j∈Z (resp. (θ 0 j ) j∈Z ) such that f 0 k,j = (∆x) -1 [(j-1/2)∆x,(j+1/2)∆x) f (0, x, v k ) dx, j ∈ Z, k ∈ {1, . . . , 2p} (resp.
θ 0 j = (∆x) -1 [(j-1/2)∆x,(j+1/2)∆x) θ(0, x) dx, j ∈ Z ;
• iterate the following three-step scheme :

1. compute θ n+1 ε using a finite difference scheme that takes into account correctly the convection term. For a scheme of order one, one can take an explicit three-point scheme for the Laplacian and a upwind-type discretization for the convection term. Indeed, this is safer than a centered discretization, which can be expected to be strongly oscillatory, since the Peclet number Pe = ||u|| ∞ ∆x/D is of order 1 (the critical value to ensure stability being Pe crit = 2, see [START_REF] Natalini | Metodi Numerici per Problemi di Evoluzione[END_REF] p. 44, [START_REF] Richtmyer | Difference methods for initial-value problems[END_REF] p. 197): choosing a upwind-scheme can be interpreted as smoothing since it increases artificially the diffusion. The Courant number µ = ∆t/(∆x) 2 is then chosen so that µ = 1/4. 2. solve, for each k ∈ {1, . . . , 2p},

∂ t f ε,k + v k ∂ x gε,k = σ s u 2 3v 2 k - 1 3 ρ ε u + v k gε,k -2 vg ε,k + σ a (B(θ ε ) -ρ ε ), (53)
thanks to a convection-type spatial discretization : for example, one can take an upwind-type D k for the transport term -v k ∂ x gε and one gets, for each k ∈ {1, . . . , 2p},

f n+1/2 ε,k = f n ε,k + ∆t D k (g n ε,k ) +∆t σ s u n 2 3v 2 k - 1 3 ρ n ε u n + v k gn ε,k -2 vg n ε,k + σ a (B(θ n ε ) -ρ n ε ) , ρ n+1/2 ε = 2p k=1 f n+1/2 ε,k , g n+1/2 ε,k = f n+1/2 ε,k -ρ n+1/2 ε ε , k ∈ {1, . . . , 2p}; 3. solve, for each k ∈ {1, . . . , 2p},      ∂ t g ε,k = -σ s g ε,k ε 2 - v k ∂ x ρ ε ε 2 + σ s 4v k uρ ε ε 2 , ∂ t f ε,k = -σ s f ε,k -ρ ε ε 2 - v k ∂ x ρ ε ε + σ s 4v k uρ ε ε . (54)
Note that, since

ρ ε = f ε , ∂ t ρ ε vanishes in Step 3 by (54), that is ρ n+1 ε = ρ n+1/2 ε
! Therefore, one has a semi-discrete explicit expression for f ε and gε :

gn+1 ε,k = e -σs∆t/ε 2 gn+1/2 ε,k + (1 -e -σs∆t/ε 2 ) 1 σ s Dk (ρ n+1/2 ε ) + 4u n v k ρ n+1/2 ε , (55)
f n+1 ε,k = e -σs∆t/ε 2 f n+1/2 ε,k + (1 -e -σs∆t/ε 2 ) ρ n+1/2 ε + ε 1 σ s Dk (ρ n+1/2 ε ) + 4u n v k ρ n+1/2 ε , ( 56
)
where Dk is a spatial discretization of -v k ∂ x . The numerical approximation of the convective term -v k ∂ x ρ ε that appears in the expression of gn+1 ε is indeed really intricate. However tempting it may seem to take Dk = D k , this naive choice can lead to misleading approximations. This can be readily remarked at once in the basic case u = 0 and σ a = 0 which leads to a system of non-coupled pure heat equations

   ∂ t ρ - 1 3σ s ∂ 2 xx ρ = 0 ∂ t θ -D∂ 2 xx θ = 0.
Choosing an upwind-type discretization for D k leads to choosing the opposite direction for Dk in order to have a three point centered scheme for the final heat equations. Indeed, the limit-scheme at ε = 0 for ρ is

ρ n+1 = ρ n + ∆t σ s k Dk D k (ρ n ),
so that, if Dk = D k for all k, one finds the scheme

ρ n+1 j = ρ n j + ∆t 3σ s ∆x 2 (ρ n j -2ρ n j-1 + ρ n j+1 ),
which is unconditionnally L 2 -unstable. Moreover, taking a close look at the expression of f n+1 ε , we note that the maximum principle is not true in general for (f n ε ) n : Formula (56) can produce negative values. Note that the mean value of the O(ε) term in ε in (56) vanishes, so that it has no influence on our final computation of ρ ε and θ ε at order O(ε). So, neglecting the terms of order ε, f n+1 ε clearly appears as a convex combination of non-negative quantities :

f n+1 ε = e -σs∆t/ε 2 f n+1/2 ε + (1 -e -σs∆t/ε 2 )ρ n+1/2 ε . (57)
Finally, we use Equation (57) to compute f ε and Equation (55) for gε .

Remark 18 Let us point out the following facts:

• The important ratio here is of course σ s ∆t/ε 2 : if σ s and ∆t are fixed, then, as ε tends to zero, one finds that gε,k tends to -v k ∂ρ + 4uv k ρ and this limit is compatible with the limit (43) of g ε which was defined by (38) when taking into account Lemma 1. Numerically, we will of course consider that there exists a constant C such that ε ≤ C∆x. Since we are interested in the limit ε → 0, this assumption is satisfactory.

• Having got rid of the O(ε) terms in the computations of f ε does not contradict the fact we also work with gε : indeed, in the end, we are only interested in viewing the results for the macroscopic quantities ρ ε and θ ε .

Let us now take a look at the numerical results : at first, we take data that satisfy the hypotheses of Theorem 1 at t > 0; then we stretch the computations to a bulk velocity that violates (H0), namely a shock wave, and to a vanishing thermal diffusion coefficient D (see Equation ( 20)).

Results in the case of a rarefaction wave

Let us now focus on the results we obtained using this scheme (movies are available at the address http://math.univ-lille1.fr/lafitte/rt).

Consider at first that the bulk velocity u is a rarefaction wave, solution of the Burgers equation

∂ t u + ∂ x u 2 /2 = 0, (58)
such that u(0, x) = -10 if x < 0 and u(0, x) = 10 if x > 0.

We represent hereafter the results for ε = 0, computed through the following steps :

Step 1 is the same as before,

Step 2' is

f n+1/2 0,k = f n 0,k + ∆t D k (g n 0,k ) +∆t σ s u 2 3v 2 k - 1 3 ρ n 0 u n + (v k gn 0,k -2 vg n 0,k >) + σ a (B(θ n 0 ) -ρ n 0 ) , (59)
ρ n+1/2 0 = 2p k=1 f n+1/2 0,k , (60) (61)
Step 3' is merely gn+1

0,k = 1 σ s Dk (ρ n+1/2 0 ) + 4u n v k ρ n+1/2 0 , (62)
f n+1 0,k = ρ n+1/2 0 , (63)
We choose the number of positive directions v to be 5.

• Thermal equilibrium

Here, the initial data ρ 0 , f 0 and θ 0 are chosen such that the initial system is at the thermal equilibrium, that is Q a,ε (f 0 , θ 0 ) = 0 :

θ 0 = ρ 0 = 1 [-2,2]
. Figure 1 represents the speed u, the temperature θ and the intensity ρ at four different times. The interesting aspects here are the facts that the thermal equilibrium is obviously unstable, that is the absorption phenomenon is clearly not the main one, and that the diffusion coefficients are very different for ρ and θ : the slopes of ρ are steeper than those of θ and a gap appears in ρ.

• Disjoint supports and different amplitudes Assume now the temperature of the fluid is high on a compact interval and that a high intensity is observed at some other place. Figure 2 represents the speed u, the temperature θ and the intensity ρ at four different times. One notes at once three interesting phenomena that are quite in accordance with our model (37) :

the convection in the opposite directions following the bulk velocity and the fact that the energy stays localized,

the rising of temperature (resp. intensity) where initially only the intensity (resp. the temperature) was high,

the diffusive aspects : the temperature θ is, quite in accordance with (37), more spread out than the intensity ρ, since the diffusion coefficient D = 1 is larger than the limit coefficient 1/3 = v 2 dv.

Remark 19 Numerically, the energy ρ + θ is conserved when ε = 0, but not otherwise : of course, the difference is of order O(ε), in accordance with our assumptions (see [START_REF] Lions | Diffuse limit for finite velocity Boltzmann kinetic models[END_REF] and Corollary 2). This is due to the fact that the bulk velocity u is given, so that there is indeed a hidden pressure term in [START_REF] Golse | Regularity of the moments of the solution of a transport equation[END_REF].

If we had built our model such that the total energy was conserved, then we would have had a non zero RHS in (58).

• Comparing the results for different values of ε Figure 3 represents the L 1 spatial distance between the solutions θ ε (resp, ρ ε ) and the relaxed solution θ 0 (resp. ρ 0 ) at the same time tmax, using five different discretizations. The different values that we chose for ε are proportional to √ σ s ∆t, since the important ratio is σ s ∆t/ε 2 (see Remarks 18). In our case, since we chose ∆t = ∆x 2 /4, ε is simply proportional to ∆x : the computation were made using 2ε √ σ s ∆x ∈ {1, 1.5, 2.3, 3, 4, 5, 6, 7, 10, 20, 30, 50, 60}.

At first, one sees that the convergence as ε tends to 0 seems to be quite fast : the error is constant after a short while. In fact, we find that the L 1 -error for ρ is of order O((∆x) 2 ). This is not surprising and is due to the splitting method. Moreover, at ε = 0, computing ρ 0 through (59), (60) and then (62) and ( 63) is equivalent to treating directly the drift-diffusion limit equation for ρ with an upwinddiscretization for the convection term and the classical centered three-point scheme for the diffusion term. What is more interesting is that the L 1 -error for θ is also of order O((∆x) 2 ) if σ a > 0, that is as soon as the coupling induces an intervention of ε in the equation for θ. Finally, the order of θ is better also because we chose a splitting method.

Results in the case of a shock wave

One can at once guess that the term ρ∂ x u in (35) will create numerical instabilities in the case of a shock wave. Hereafter are two cases that illustrate this phenomenon.

• Thermal equilibrium Figure 4 shows that the fact that the initial data is nonzero around x = 0 implies that there is an immediate blow-up for ρ around x = 0, but this blow-up is delayed for θ because of the effects of the diffusion.

• Disjoint supports

We witness in Figure 5 the fact that it is indeed precisely at x = 0 that the blow-up occurs : the energy travels towards x = 0 since the bulk velocity is positive on the left and negative on the right. When the two parts of the solution get at x = 0, the energy concentrates and there is an blow-up.

7.5 Changing the parameters σ a and σ s

• Changing σ a Figure 6 allows us to compare different evolutions depending on different values for σ a in the case of a rarefaction wave: one notes at once that, although the initial data are symmetric, the solution is not when σ a = 1 but gets more and more symmetric as σ a increases. Indeed, the coupling gets bigger and consequently the phenomena tend to balance each other.

In Figure 7, in the case of a shock wave, we note that increasing σ a delays the blow-up for ρ. Moreover, the energy of ρ transfers to θ a lot more quickly when σ a is bigger.

• Changing σ s Increasing σ s means that the coefficient diffusion in the limit equation (35) decreases, and this is the reason why we note an obvious sharpening of the slopes of ρ in Figure 8 and oscillations. What we also note is that our increasing σ s speeds up the blow-up for θ.

Results when there is no thermal diffusion

The major effect that one notes is the sharpening of the slopes of θ, as shown in Figure 9. This is of course not surprising since the equation for θ is now a purely convective one. Recall however that the assumption D > 0 was used in the proof of convergence (see Subsection 5.3) as a means to obtain some strong compactness propreties. From a physical viewpoint, the fact that the scheme is able to reproduce the nondiffusive case is very satisfactory.

The expression of the loss term follows immediately. By using (8), we obtain

Q - s (ν, v) = ν ν 0 2 S 2 σ 0 s (ν 0 , v 0 ′ , v 0 ) dv 0 ′ f 0 (ν 0 , v 0 ) = ν ν 0 2 S 2 σ 0 s (ν 0 , v 0 ′ , v 0 ) 1 Λ(v ′ ) 2 dv ′ ν 0 ν 3 f (ν, v) = S 2 σ s (ν, v ′ , v) 1 Λ(v ′ ) 2 dv ′ Λf (ν, v)
where we have set σ s (ν, v, v ′ ) = σ 0 s (ν 0 , v 0 , v 0 ′ ), ν 0 and v 0 being associated to ν, v by ( 6) as v 0 ′ is defined from v ′ and ν 0 .

Then, we turn to the gain term by using the same tricks. We get

ν 0 ν 2 Q + s (ν, v) = S 2 σ 0 s (ν 0 , v 0 , v 0 ′ ) f 0 (ν 0 , v 0 ′ ) dv 0 ′ = lim k→∞ ∞ 0 S 2 σ 0 s (µ 0 , v 0 , v 0 ′ ) f 0 (µ 0 , v 0 ′ )ζ k (µ 0 -ν 0 ) dv 0 ′ dµ 0 = lim k→∞ ∞ 0 S 2 σ 0 s (µ 0 , v 0 , v 0 ′ ) µ 0 µ ′ 3 f (µ ′ , v ′ )ζ k (µ 0 -ν 0 ) dv 0 ′ dµ 0
(use of ( 8), (µ ′ , v ′ ) functions of (µ 0 , v 0 ′ ) by ( 6))

= lim k→∞ ∞ 0 S 2 σ 0 s (µ 0 , v 0 , v 0 ′ ) µ 0 µ ′ 3 f (µ ′ , v ′ )ζ k (µ 0 -ν 0 ) µ ′ µ 0 dv ′ dµ ′ (change of variables (µ 0 , v 0 ′ ) → (µ ′ , v ′ ) with (7)) = S 2 lim k→∞ ∞ 0 σ 0 s (µ 0 , v 0 , v 0 ′ ) µ 0 µ ′ 2 f (µ ′ , v ′ )ζ k (µ 0 -ν 0 ) µ ′ µ 0 dµ 0 dv ′ (change of variables µ ′ → µ 0 with v ′ fixed) = S 2 σ 0 s (ν 0 , v 0 , v 0 ′ )Λ(v ′ )f (ν ′ , v ′ ) dv ′ = S 2 σ s (ν, v, v ′ )Λ(v ′ )f (ν ′ , v ′ ) dv ′ .

C Proof of Proposition 2

As a preliminary remark, note that we can work with large enough values of the parameter C so that for some constants λ * , λ * ,

0 < λ * ≤ Λ = 1 -v • u/C 1 -u 2 /C 2 ≤ λ * (65)
holds for any C ≥ C * and almost all (t, x) ∈ R + × R 3 , the starred quantities depending on (H0). Furthermore, we check readily that in this context we have also

     |∇ x Λ| ≤ K C , |∂ t Λ| ≤ K C 0 < k ≤ γ ≤ 1, 0 < k C 2 ≤ 1 - 1 γ ≤ K C 2 (66)
for some constants K, k > 0 which depend only on (H0), but do not depend on the scaling parameters C, P, L s and L a . We adopt from now on the convention that K stands for a constant which depends only on (H0), (H1), but does not depend on the scaling parameters disregarding the fact that the value of the constant may change from a line to another.

As usual, we discuss the result by using the equations and performing formal integration by parts. A rigorous and complete proof can be obtained through standard regularization arguments and in the construction of the solution. Indeed, we write

d dt H(f, θ) = R 3 S 2 Φ ′ (Λ 4 f )∂ t f Λ dv dx + C P R 3 Ψ ′ (θ)∂ t θ dx + R 3 S 2 4Φ ′ (Λ 4 f )f -3Φ(Λ 4 f ) 1 Λ 4 ∂ t Λ dv dx := -R 1 = C 1 + CL s R 3 S 2 Φ ′ (Λ 4 f )Q s (f ) Λ dv dx +CL a R 3 S 2 Φ ′ (Λ 4 f )Q a (f, θ) Λ dv dx -CL a R 3 Ψ ′ (θ) Λ γ Q a (f, θ) dx +C 2 - C P R 3 Ψ"(θ)|∇ x θ| 2 dx -R 1 (67)
with the notation

C 1 = -C R 3 S 2 Φ ′ (Λ 4 f )v • ∇ x f Λ dv dx C 2 = C P R 3 Ψ"(θ)∇ x θ • uθ dx.
(Note that we have used Remark 7.) By using (65), we estimate easily

|R 1 | ≤ K R 3 S 2 Φ(Λ 4 f ) + Φ ′ (Λ 4 f )Λ 4 f |∂ t Λ| 1 Λ 3 dv dx.
It yields ( 27) by using (66) and the assumption on the function Φ.

The second term in the right-hand side of (67) has been treated in Proposition 1: it is exactly -D s . The third and forth terms combine. We obtain

CL a R 3 S 2 Φ ′ (Λ 4 f ) -Ψ ′ (θ) 1 γ Q a (f, θ) Λ dv dx = CL a R 3 S 2 σ a γΛ 2 Φ ′ (Λ 4 f ) -Ψ ′ (θ) B(θ) -Λ 4 f dv dx := -D a +CL a R 3 S 2 σ a 1 - 1 γ B(θ) -Λ 4 f Φ ′ (Λ 4 f ) Λ 2 dv dx := R 2 .
The remaining term R 2 is readily estimated by using (66).

It remains to deal with the convective terms C 1 and C 2 . At first, we simply use the Cauchy-Schwarz inequality to dominate

|C 2 | ≤ C P R 3 Ψ"(θ)|∇ x θ| Ψ"(θ)|uθ| dx ≤ C 2P R 3 Ψ"(θ)|∇ x θ| 2 dx + u 2 L ∞ R 3 Ψ"(θ) θ 2 dx .
Eventually, integrations by part yield

C 1 = C R 3 S 2 Φ ′ (Λ 4 f ) v • - ∇ x (Λ 4 f ) Λ 4 + f Λ 4 ∇ x Λ 4 Λ dv dx = C R 3 S 2 Φ(Λ 4 f ) v • ∇ x (Λ -3 ) dv dx +C R 3 S 2 Φ ′ (Λ 4 f )Λ 4 f v • ∇ x (Λ 4 ) Λ 4 1 Λ 3 dv dx.
Then, we conclude that (25) holds by using (66). This ends the proof.

D Proof of Corollary 1

Our task consists of dominating the right-hand side in [START_REF] Natalini | Metodi Numerici per Problemi di Evoluzione[END_REF] so that we will be able to use the Gronwall lemma. It follows from Proposition 2 with suitable choices of the functions Φ and Ψ. First, notice that C 1 and R 1 are clearly dominated by H(f, θ). Next, the first term in the right-hand side of ( 26) is compensated by the second integral in the left-hand side of [START_REF] Natalini | Metodi Numerici per Problemi di Evoluzione[END_REF]. Besides, assuming that θ 2 Ψ"(θ) ≤ M Ψ(θ) allows us to dominate the second term by H(f, θ). The treatment of the terms coming from R 2 is slightly more delicate. Using the Young inequality yields

|R 2 | ≤ KνCL a R 3 S 2 |Λ 4 f -B(θ)| 2 dv dx + KL a 4νC 3 R 3 S 2 |Φ ′ (Λ 4 f )| 2 dv dx Λ 4 , (68)
for any ν ∈ (0, 1). Note that, up to now, we do not use any relation between the functions Φ and Ψ. Then, we first use these relations with Φ(z) = z 2 /2, and, according to Remark 10, Ψ(z) = z 0 B(θ) dθ, which has been denoted by Γ. With this choice, we simply have

H(f, θ) = 1 2 R 3 S 2 f 2 Λ 5 dv dx + C P R 3 Γ(θ) dx.
Next, using (H1) and (65), we get

D a ≥ σ * (λ * ) 2 CL a R 3 S 2 |Λ 4 f -B(θ)| 2 dv dx, while D s ≥ σ * 2(λ * ) 4 CL s R 3 S 2 |Λ(v) 4 f (v) -Λ(v ′ ) 4 f (v ′ )| 2 dv ′ dv dx.
On the other hand, we have

|Φ ′ (Λ 4 f )| 2 Λ 4 ≤ Kf 2 Λ 5 ,
in view of (65). Coming back to (68), we thus obtain

|R 2 | ≤ KνCL a R 3 S 2
|Λ 4 f -B(θ)| 2 dv dx + KL a 4νC 3 H(f, θ).

Combining these informations and choosing, say, ν = σ * 2(λ * ) 2 , we are finally led to

d dt H(f, θ) + C 2P R 3 |∇ x G(θ)| 2 dx + σ * 2(λ * ) 2 CL a R 3 S 2 |Λ 4 f -B(θ)| 2 dv dx + σ * 2(λ * ) 4 L s C R 3 S 2 |Λ(v) 4 f (v) -Λ(v ′ ) 4 f (v ′ )| 2 dv ′ dv dx ≤ K + 2(λ * ) 2 L a 4σ * C 2 H(f, θ). (69)
A direct application of the Gronwall lemma, with ( 29) and (H3), proves the second estimates in ( 30) and (32), the L 2 -estimate in [START_REF] Tartar | Remarks on homogenization in Homogenization and effective moduli of material and media[END_REF] and it also justifies (33) and (34).

We go one step further by using now another choice for the functions Φ and Ψ. Namely, let us use Ψ(θ) = θ 2 /2 and Φ(z) = z 0 B (-1) (s) ds (which has been denoted by J(z)). As explained above, the main difficulty relies on the term denoted by R 2 . We simply write

|R 2 | ≤ K L a C R 3 S 2 Λ 4 f |J ′ (Λ 4 f )| dv dx Λ 2 + R 3 S 2 B(θ) |J ′ (Λ 4 f )| dv dx Λ 2 . ( 70
)
Assumption (H2) allows to dominate the first term by H(f, θ). Finally, the Young inequality yields

R 3 S 2 B(θ) J ′ (Λ 4 f ) dv dx Λ 2 ≤ 1 p R 3 B(θ) p dx + 1 p ′ R 3 S 2 J ′ (Λ 4 f ) Λ 2 p ′ dv dx ≤ K R 3 Γ(θ) dx + R 3 S 2 J(f Λ 4 ) 1 Λ 3 dv dx ,
where the last line follows by using (H2). Therefore, we can now complete the estimates in ( 31) and ( 30) by means of the Gronwall lemma, using ( 29), (H3) and the bound on Γ(θ) which has been justified in the previous step.

E Proof of Corollary 2

It is more convenient to deal with a modified energy involving the weight Λ. Precisely, we have d dt

R 3 S 2 f Λ dv dx + C P R 3 θ dx = R 3 S 2 f ∂ t Λ dv dx -C R 3 S 2 v • ∇ x f Λ dv dx +CL a R 3 S 2 1 - 1 γ ΛQ a dv dx = R 3 S 2 f ∂ t Λ dv dx + C R 3 S 2 vf • ∇ x Λ dv dx +CL a R 3 S 2 1 - 1 γ B(θ) Λ 2 dv dx -CL a R 3 S 2 1 - 1 γ Λ 2 f dv dx.
We use the assumption (H2) which implies Γ(θ) Furthermore, coming back to (66)and ( 29), we note that all the quantities C∇ x Λ, ∂ t Λ and CL a (1 -1/γ) ≤ KL a /C are bounded uniformly with respect to the scaling parameters. We thus obtain d dt

R 3 S 2 f Λ dv dx + R 3 θ dx ≤ K R 3 S 2 f Λ dv dx + R 3 θ dx + R 3 Γ(θ) dx
for some constant K which does not depend on the parameters C, P, L s , L a . An application of the Gronwall lemma ends the proof.

F Proof of Proposition 3

Proposition 1 applies for u = 0, Λ = 1 1. Thus, we have Ker(Q s,0 ) = Span{1 1}, Q s,0 (f ) = 0, which shows that Ran(Q s,0 ) ⊂ Span{1 1} ⊥ . In particular the vanishing condition h = 0 is a necessary condition to solve the equation Q s,0 (f ) = h. Moreover, the following coerciviy estimate

- S 2 Q s,0 (f )f dv ≥ σ * 2 S 2 S 2 |f (v ′ ) -f (v)| 2 dv ′ dv = σ * S 2 |f (v) -f | 2 dv
holds. Therefore, a(f, g) = -S 2 Q s,0 (f )g dv defines a bilinear continuous form on L 2 (S 2 ), which is coercive on the closed set {f ∈ L 2 (S 2 ), f = 0}. Hence, the equation Q s,0 (f ) = h can be solved on this set by a direct application of the Lax-Milgram lemma. Estimates on the solution are immediate consequences of the coercivity relation. (Much more general collision operators are dealt with in [START_REF] Degond | Diffusion limit for non homogeneous and non reversible processes[END_REF].)

3 B

 3 ) 1/p ds ≤ M θ 0 ds = M θ. Then, we use the Hölder inequality with exponent q = 1/(2 -p), q ′ = 1/(p -1) and we get R

Figure 1 :Figure 2 :Figure 3 :Figure 4 :Figure 5 : 10 Figure 6 :

 12345106 Figure 1: Speed, temperature and intensity -tmax = 0.3, σ s = σ a = 1, ∆x = 0.02

Figure 7 : 20 Figure 8 :Figure 9 :

 72089 Figure 7: Temperature and intensity -tmax = 0.35, σ s = 1, ∆x = 0.02, u ∞ = 10, shock wave

We denote by Dxu the jacobian matrix of the vector field u whose components are ∂x j u i .

Acknowledgments. We are gratefully indebted to Ch. Buet and B. Després who gave us the preliminary version of ref. [7] and for very fruitful discussions. We also warmly thank D. Ramon, from society Hygeos, for many advice and interesting comments on the physics of the problem.

Appendix

A Table of Physical Constants

For the sake of information, let us recall here the value of some physical constants used in the paper.

B Expression of the Scattering Operator

The computations start with the following preliminary result which shows how some weights arise from the change of variables [START_REF] Bouchut | Diffusive BGK approximations for nonlinear multidimensional parabolic equations[END_REF].

Lemma 8 For any (smooth enough) function ϕ, we have

where in the right-hand side v 0 ′ is understood as the function of ν 0 and v ′ defined by the formulae [START_REF] Bouchut | Diffusive BGK approximations for nonlinear multidimensional parabolic equations[END_REF] and

Proof. We follow the arguments given in [START_REF] Ch | Asymptotic analysis of fluid models for the coupling of radiation and hydrodynamics[END_REF]. Namely, we use the invariance property ( 7) by considering a sequence of mollifiers ζ k : R → R + . Let us denote the integral in the left-hand side of (64) by I 0 . Then, we compute as follows

ϕ(ν 0 , v 0 ′ ) 1 Λ(v ′ ) 2 dv ′ by using the definition Λ(v ′ ) = µ 0 /µ ′ .