Kazuo Aoki
email: aoki.kazuo.7a@kyoto-u.ac.jp

Shingo Kosuge
email: kosuge.shingo.6r@kyoto-u.ac.jp

Prof Tai-Ping Liu

Academia Sinica (New Series) THE STEADY BOLTZMANN AND NAVIER-STOKES EQUATIONS

Keywords: AMS Subject Classification: 35Q30, 35Q20 (76P05, 76D05, 82C40) Steady Boltzmann equation, Steady Navier-Stokes equation, Heat diffusion, Viscous heating, Periodic solutions

The paper discusses the similarities and the differences in the mathematical theories of the steady Boltzmann and incompressible Navier-Stokes equations posed in a bounded domain. First we discuss two different scaling limits in which solutions of the steady Boltzmann equation have an asymptotic behavior described by the steady Navier-Stokes Fourier system. Whether this system includes the viscous heating term depends on the ratio of the Froude number to the Mach number of the gas flow. While the steady Navier-Stokes equations with smooth divergence-free external force always have at least one smooth solutions, the Boltzmann equation with the same external force set in the torus, or in a bounded domain with specular reflection of gas molecules at the boundary may fail to have any solution, unless the force field is identically zero. Viscous heating seems to be of key importance in this situation. The nonexistence of any steady solution of the Boltzmann equation in this context seems related to the increase of temperature for the evolution problem, a phenomenon that we have established with the help of numerical simulations on the Boltzmann equation and the BGK model.

Introduction

The Boltzmann equation and its hydrodynamic limits have been widely studied in the time-dependent regime. The Cauchy problem for the Boltz-mann equation is discussed in [START_REF] Ukai | On the existence of global solutions of a mixed problem for the nonlinear Boltzman equation[END_REF][START_REF] Kaniel | The Boltzmann Equation: Uniqueness and Local Existence[END_REF][START_REF] Illner | The Boltzmann equation: Global existence for a rare gas in an infinite vacuum[END_REF][START_REF] Diperna | On the Cauchy problem for the Boltzmann equation: Global existence and weak stability results[END_REF]. Hydrodynamic limits of the Boltzmann equation are analyzed by various methods in [START_REF] Nishida | Fluid dynamical limit of the nonlinear Boltzmann equation to the level of the compressible Euler equation[END_REF][START_REF] Caflisch | The fluid dynamic limit of the nonlinear Boltzmann equation[END_REF][START_REF] Bardos | Sur les limites asymptotiques de la théorie cinétique conduisant à la dynamique des fluides incompressibles[END_REF][START_REF] Demasi | Incompressible Navier-Stokes and Euler Limits of the Boltzmann Equation[END_REF][START_REF] Bardos | Fluid dynamic limits of the Boltzmann equation I[END_REF][START_REF] Bardos | Fluid dynamic limits of kinetic equations II: Convergence proofs for the Boltzmann equation[END_REF][START_REF] Bardos | The classical incompressible NavierStokes limit of the Boltzmann equation[END_REF][START_REF] Bardos | The acoustic limit for the Boltzmann equation[END_REF][START_REF] Lions | From Boltzmann equation to the NavierStokes and Euler equations I[END_REF][START_REF] Lions | From Boltzmann equation to the NavierStokes and Euler equations II[END_REF][START_REF] Golse | The StokesFourier and acoustic limits for the Boltzmann equation[END_REF][START_REF] Saint-Raymond | Convergence of solutions to the Boltzmann equation in the incompressible Euler limit[END_REF][START_REF] Golse | The NavierStokes limit of the Boltzmann equation for bounded collision kernels[END_REF][START_REF] Golse | The incompressible Navier-Stokes limit of the Boltzmann equation for hard cutoff potentials[END_REF][START_REF] Levermore | From the Boltzmann equation to an incompressible Navier-Stokes-Fourier system[END_REF] -see also [START_REF] Villani | Limites hydrodynamiques de léquation de Boltzmann[END_REF] for a nice introduction to the mathematical analysis of hydrodynamic limits of the Boltzmann equation. Some of the results known in the case of the Cauchy problem set on a spatial domain which either the Euclidean space R 3 or the periodic box T 3 have been extended to the case of the initial-boundary value problem in a domain Ω of R 3 : see [START_REF] Mischler | Kinetic equations with Maxwell boundary conditions[END_REF][START_REF] Masmoudi | From the Boltzmann equation to the Stokes-Fourier system in a bounded domain[END_REF][START_REF] Bardos | The incompressible Euler limit of the Boltzmann equation with accommodation boundary condition[END_REF][START_REF] Golse | From the Boltzmann equation to the Euler equations in the presence of boundaries[END_REF] -see also [START_REF] Saint-Raymond | Hydrodynamic Limits of the Boltzmann Equation[END_REF] for a more synthetic presentation of this material as well as further results.

By comparison, the mathematical literature on the analogous problems in the case of steady solutions is much more scarce. A great collection of asymptotic and numerical results on the Boltzmann equation in the steady regime can be found in the books [START_REF] Sone | Kinetic Theory and Fluid Dynamics[END_REF][START_REF] Sone | Molecular Gas Dynamics: Theory, Techniques and Applications[END_REF]. As for the mathematical analysis of the boundary value problem for the Boltzmann equation, the main references are [START_REF] Guiraud | Problème aux limites intérieur pour léquation de Boltzmann linéaire[END_REF][START_REF] Guiraud | Problème aux limites intérieur pour léquation de Boltzmann en régime stationnaire, faiblement non linéaire[END_REF] (see also [START_REF] Guiraud | Problème aux limites intérieures pour léquation de Boltzmann[END_REF]), together with the more modern references [START_REF] Arkeryd | The stationary Boltzmann equation in R n with given indata[END_REF][START_REF] Esposito | Non-isothermal boundary in the Boltzmann theory and Fourier law[END_REF]. Steady solutions of the Boltzmann equation for a gas flow past an obstacle have been investigated in detail in [START_REF] Ukai | Steady solutions of the Boltzmann equation for a gas flow past an obstacle. I. Existence[END_REF][START_REF] Ukai | Steady solutions of the Boltzmann equation for a gas flow past an obstacle. II. Stability[END_REF].

There are striking analogies between the Boltzmann and the incompressible Navier-Stokes equations in three space dimensions -in the words of P.-L. Lions [START_REF] Lions | Compactness in Boltzmanns equation via Fourier integral operators and applications II[END_REF] "[...] the global existence result of [renormalized] solutions [...] can be seen as the analogue for Boltzmann's equation to the pioneering work on the Navier-Stokes equations by J. Leray". This analogy is at the origin of the program outlined in [START_REF] Bardos | Sur les limites asymptotiques de la théorie cinétique conduisant à la dynamique des fluides incompressibles[END_REF][START_REF] Bardos | Fluid dynamic limits of the Boltzmann equation I[END_REF][START_REF] Bardos | Fluid dynamic limits of kinetic equations II: Convergence proofs for the Boltzmann equation[END_REF] and carried out in [START_REF] Golse | The NavierStokes limit of the Boltzmann equation for bounded collision kernels[END_REF][START_REF] Golse | The incompressible Navier-Stokes limit of the Boltzmann equation for hard cutoff potentials[END_REF] -see also [START_REF] Levermore | From the Boltzmann equation to an incompressible Navier-Stokes-Fourier system[END_REF] for the extension to a more general class of collision kernels.

It has been known for a long time that the regularity theory of solutions of the incompressible Navier-Stokes equations is much easier in the steady than in the time-dependent regime. For instance, in space dimension 3, steady solutions of the Dirichlet problem for the incompressible Navier-Stokes equations have the same regularity as their boundary data and the external force field driving them: see Proposition 1.1 and Remark 1.6 in chapter II, §1 of [START_REF] Temam | Navier-Stokes equations[END_REF]. At variance, it is still unknown at the time of this writing whether Leray solutions of the Cauchy problem for the Navier-Stokes equations in space dimension 3 propagate the regularity of their initial data -see Problems A-B in [START_REF] Fefferman | Existence and smoothness in the Navier-Stokes equations[END_REF].

One striking difference between the steady and the time-dependent problems for the incompressible Navier-Stokes equations is the uniqueness theory, and its relation to the regularity of solutions. Smooth solutions of the timedependent incompressible Navier-Stokes equations in space dimension 3 are known to be uniquely determined by their initial data (and driving force field) within the class of Leray weak solutions, a remarkable result proved by Leray himself (see §32 in [START_REF] Leray | Essai sur le mouvement dun liquide visqueux emplissant lespace[END_REF]). On the contrary, it can be proved that bifurcations do occur on the steady problem for the incompressible Navier-Stokes equations, leading to nonuniqueness results. Such a nonuniqueness result for the Taylor-Couette problem has been proved in [START_REF] Velte | Stabilität und Verzweigung stationärer Lösungen der Navier-Stokesschen Gleichungen beim Taylor Problem[END_REF] -see also Chapter II, §4 in [START_REF] Temam | Navier-Stokes equations[END_REF]. Analogous bifurcations for the Boltzmann equation have been observed numerically in [START_REF] Sone | Analytical study of bifurcation of a flow of a gas between coaxial circular cylinders with evaporation and condensation[END_REF], and mathematically in [START_REF] Arkeryd | On a Taylor-Couette type bifurcation for the stationary nonlinear Boltzmann equation[END_REF].

These considerations suggest studying the mathematical theory of existence and regularity for steady solutions of the Boltzmann equation driven by an external force field. In particular, does the assumption of a steady regime simplify regularity issues, as in the case of the Navier-Stokes equations?

As we shall see, there are strong similarities between the Boltzmann and the Navier-Stokes equations in the steady regime. After reviewing the basic structure of the Boltzmann equation in section 1, we propose in section 2 a formal derivation of two variants of the incompressible Navier-Stokes-Fourier system from the Boltzmann equation under appropriate scaling assumptions. Section 3 discusses the steady Navier-Stokes and Boltzmann equation in the periodic setting, which leads to a striking difference between both models. Section 4 provides a physical explanation for this difference, based on numericla simulations of the evolution problem. After a brief section 5 summarizing our conclusions, some computations involving Gaussian averages of certain vector and tensor fields have been put together in an appendix.

Professor Tai-Ping Liu is at the origin of some of the most striking results on the mathematical analysis of the equations fluid dynamics -his work [START_REF] Liu | Solutions in the large for the equations of nonisentropic gas dynamics[END_REF] on the compressible Euler system, and his analysis of the stability of the Boltzmann shock profile in collaboration with S.-H. Yu [START_REF] Liu | Boltzmann equation: Micro-macro decompositions and positivity of shock profiles[END_REF] have had a lasting impact on our field. We are pleased to offer him this modest contribution on the occasion of his 70th birthday.

The Steady Boltzmann Equation with External Force Field

The Boltzmann equation with external force field f ≡ f (x) ∈ R 3 is posed in a bounded, convex spatial domain Ω ⊂ R 3 with smooth boundary ∂Ω. The outward unit normal vector field on ∂Ω is denoted by n x . According to §1.9 in [START_REF] Sone | Molecular Gas Dynamics: Theory, Techniques and Applications[END_REF], its dimensionless form is

v • ∇ x F + Ma 2 Fr 2 f • ∇ v F = 1 Kn C(F ) , x ∈ Ω , v ∈ R 3 , (1.1)
where

C(F ) := R 3 ×S 2 (F F * -F F * )|(v -v * ) • ω|dv * dω ,
with the usual notation

F = F (x, v) , F * = F (x, v * ) , F = F (x, v ) , and F * = F (x, v * ) , assuming that v = v -(v -v * ) • ωω , and v * = v * + (v -v * ) • ωω .
The dimensionless numbers Ma, Fr and Kn are respectively the Mach, Froude and Knudsen numbers. We recall the definitions of these dimensionless numbers:

Ma = U 0 √ RT 0 , Fr = U 0 √ F 0 L 0 , Kn = 0 L 0 ,
where R is the specific gas constant, henceforth set to one for simplicity. In these formulas, U 0 , T 0 , F 0 , L 0 are respectively the reference speed, temperature and external force in the gas, while L 0 is the reference length scale and 0 is the mean free path of the gas molecules at the reference state.

We recall that the local conservation laws of mass, momentum and energy for the collision integral are

R 3 C(F )dv = R 3 v i C(F )dv = R 3
|v| 2 C(F )dv = 0 a.e. on Ω for i = 1, 2, 3, provided that F ≥ 0 is measurable and decays rapidly enough as |v| → ∞ -so that, for instance

Ω×R 3 (1 + |v| 3 )(1 + |f |)F dvdx < ∞ .
(See for instance chapter I, §4, Corollary 1 in [START_REF] Glassey | The Cauchy Problem in Kinetic Theory[END_REF].) This implies the local mass, momentum and energy balance identities, which are the differential identities

div x R 3 vF dv = 0 , (mass) div x R 3 v ⊗ vF dv = Ma 2 Fr 2 f R 3 F dv , (momentum)
div x R 3 v 1 2 |v| 2 F dv = Ma 2 Fr 2 f • R 3 vF dv , (energy)
are satisfied by all the solutions of the Boltzmann equation having the decay property mentioned above. Notice that this property implies that

R 2 n Ω |v|=Rn F |f |ds(v)dx → 0
for some sequence R n → ∞. (The notation ds designates the surface element on the sphere of radius R n centered at the origin.)

Integrating further in x and applying Green's formula leads to the identities:

∂Ω R 3 F v • n x dvdS(x) = 0 , ( mass
) ∂Ω R 3 vF (v • n x )dvdS(x) = Ma 2 Fr 2 Ω×R 3 f F dxdv , ( momentum
) ∂Ω R 3 1 2 |v| 2 F (v • n x )dvdS(x) = Ma 2 Fr 2 Ω×R 3 v • f F dxdv , ( energy
)
where dS is the surface element on ∂Ω, which are the global balance laws of mass, momentum and energy.

Multiplying both sides of the Boltzmann equation by ln F + 1 leads to the identity

v • ∇ x (F ln F ) + Ma 2 Fr 2 f • ∇ v (F ln F ) = 1 Kn C(F )(ln F + 1) .
Integrating in v leads to the local form of Boltzmann's H theorem (see for instance chapter I, §4, Corollary 1 in [START_REF] Glassey | The Cauchy Problem in Kinetic Theory[END_REF])

div x R 3 vF ln F dv = 1 Kn R 3 C(F ) ln F dv ≤ 0 ,
assuming again that F decays rapidly enough as |v| → ∞ -for instance

Ω×R 3 (1 + |v|)(1 + |f |)F ln F dvdx < ∞ , so that R n Ω |v|=Rn F | ln F ||f |ds(v)dx → 0 for some sequence R n → ∞.
Integrating further in x, one obtains the global form of Boltzmann's H theorem

∂Ω R 3 F ln F v • n x dvdS(x) = 1 Kn Ω×R 3 C(F ) ln F dvdx ≤ 0 .
Boltzmann's H theorem asserts that the inequality above is an equality if and only if F is a local Maxwellian, i.e. is of the form

F (x, v) = M (ρ(x),u(x),θ(x)) (v) ,
with the notation

M (ρ,u,θ) = ρ (2πθ) 3/2 e -|v-u| 2 2θ
.

(

(See for instance chapter I, § §5 and 7, Corollary 2 in [START_REF] Glassey | The Cauchy Problem in Kinetic Theory[END_REF].) .

The Navier-Stokes Limit for the Boltzmann Equation

In this section, we explain how two variants of the Navier-Stokes-Fourier system can be derived from the Boltzmann equation. The exposition is formal and follows the style adopted in [START_REF] Bardos | Sur les limites asymptotiques de la théorie cinétique conduisant à la dynamique des fluides incompressibles[END_REF][START_REF] Bardos | Fluid dynamic limits of the Boltzmann equation I[END_REF]. The difference between the limiting systems comes from the different scaling assumptions on the Froude number. Specifically, we are concerned with those variants of the incompressible Navier-Stokes-Fourier system which do, or do not include the vis-cous heating term. This particular feature of the fluid dynamic limit of the Boltzmann equation is of key importance for the discussion in the present work.

From Boltzmann to Navier-Stokes-Fourier

In this section, we assume the following scaling, where > 0 is a small parameter: Ma = Kn = , and Fr = √ .

(2.1)

Hence, the scaled Boltzmann equation takes the form

v • ∇ x F + f • ∇ v F = 1 C(F ) , x ∈ Ω , v ∈ R 3 .
In addition, write the Helmholtz decomposition of the external force field as

f (x) = -∇Φ(x) + f s (x) , with div f s = 0 .
(In other words, we assume that the divergence free component of the external force is small compared to its curl free component.) With this additional assumption, the scaled Boltzmann equation becomes

v • ∇ x F -∇Φ(x) • ∇ v F + 2 f s (x) • ∇ v F = 1 C(F ) , x ∈ Ω , v ∈ R 3 . (2.2)
Seek F in the form

F (x, v) = M (v) Z e Φ(x) + g (x, v) , (2.3)
assuming that

Ω×R 3 g M dxdv = 0 , (2.4)
with the notation

M := M (1,0,1)
,

and 1 Z := 1 |Ω| Ω e Φ(x) dx . (2.5)
In other words, the distribution function is sought in the form of a perturbation of the order of the Mach number Ma = about the uniform Maxwellian M . Except for the scaling assumption on the external force, this is exactly the same scaling assumption as in [START_REF] Bardos | Sur les limites asymptotiques de la théorie cinétique conduisant à la dynamique des fluides incompressibles[END_REF][START_REF] Bardos | Fluid dynamic limits of the Boltzmann equation I[END_REF].

In terms of g , the scaled Boltzmann equation (2.2) takes the form

v • ∇ x g -2 M -1 ∇ x Φ • ∇ v (M g ) + 2 M -1 f s (x) • ∇ v M (Z e Φ(x) + g ) = -Z e Φ(x) Lg + Q(g , g ) (2.6)
where

Lg := -M -1 DC(M ) • (M g) , Q(g, g) := M -1 C(M g) , (2.7
Q(f, g) := 1 2 (Q(f + g, f + g) -Q(f, f ) -Q(g, g)) .
In other words,

Lg = R 3 ×S 2 (g + g * -g -g * )|(v -v * ) • ω|M * dv * dω , while Q(f, g) = 1 2 R 3 ×S 2 (f g * + f * g -f g * -f * g)|(v -v * ) • ω|M * dv * dω .
Henceforth, the integration with respect to the Gaussian weight M is denoted as follows:

φ := R 3 φ(v)M (v)dv .
Theorem 1. Let F be a family of solutions of the scaled Boltzmann equation (2.2), whose relative fluctuation g defined in (2.3) satisfies

g → g in L 2 (Ω × R 3 ; (1 + |v| 2 )M dvdx) weak, and
Q(g , g ) → Q(g, g) in L 1 (Ω; L 2 (R 3 ; (1 + |v| 2 )M dv)) weak, while g ∇φ(v) → 0 in L 2 (Ω) weak for each φ ∈ L 2 (R 3 ; (1 + |v| 2 )M dv). Then g(x, v) := θ + u(x) • v + θ(x) 1 2 (|v| 2 -5)
where

θ := 1 |Ω| Ω θ(x)dx
and (u, θ) is a solution of the incompressible Navier-Stokes-Fourier system

                 div u = 0 , div(u ⊗2 ) + ∇p = ν∆u + θ∇Φ + f s , 5 2 div(uθ) = κ∆θ -u • ∇Φ . (2.8)
The values of the viscosity ν and heat diffusivity κ are determined implicitly in terms of the collision integral, by formulas (2.15) and (2.16).

We recall the following fundamental result.

Lemma 1. The operator L is an unbounded self-adjoint operator on the

Hilbert space L 2 (R 3 ; M dv) with domain L 2 (R 3 ; (1 + |v| 2 )M dv). Moreover L ≥ 0 and Ker L = span{1, v 1 , v 2 , v 3 , |v| 2 } .
Finally, L satisfies the Fredholm alternative: one has

Ran L = (Ker L) ⊥ .
This has been proved by Hilbert in 1912 (see for instance [START_REF] Glassey | The Cauchy Problem in Kinetic Theory[END_REF], chapter III, § §4-5 and [START_REF] Cercignani | Theory and applications of the Boltzmann equation[END_REF], chapter IV, §6).

Proof. The proof of Theorem 1 is rather involved; we follow the discussion in [START_REF] Bardos | Sur les limites asymptotiques de la théorie cinétique conduisant à la dynamique des fluides incompressibles[END_REF][START_REF] Bardos | Fluid dynamic limits of the Boltzmann equation I[END_REF].

Step 1: asymptotic form of the number density fluctuations

Assuming that g → g weakly in L 2 (Ω × R 3 ; M dxdv), we pass to the limit in the sense of distributions in both sides of the equality (2.6), and obtain

Lg = 0 .
Thus g is of the form

g(x, v) = ρ(x) + u(x) • v + θ(x) 1 2 (|v| 2 -3) .
(2.9)

Step 2: divergence-free and hydrostatic conditions

Multiplying both sides of the scaled Boltzmann equation by 1 M and integrating in v leads to

div x vg = R 3 div v (M (g ∇Φ(x) -(Z e Φ(x) + g )f s (x)))dv = 0 .
Passing to the limit in both sides of this equality in the sense of distributions as → 0, we get div x u = div x vg = 0 .

(2.10)

Mutiplying both sides of the scaled Boltzmann equation by 1 M v and integrating in v leads to div x v ⊗2 g = -∇Φ(x) g + f s (x) (Z e Φ(x) + g ) .

Passing to the limit as → 0 shows that

∇ x (ρ + θ) = div x v ⊗2 g = 0 .
(2.11)

Passing to the limit in (2.4) shows that

Ω ρ(x)dx = Ω g (x)dx = 0 , so that ρ(x) + θ(x) = θ , with θ := 1 |Ω| Ω θ(x)dx .
(2.12)

Step 3: motion equation

Mutiplying both sides of the scaled Boltzmann equation by

1 2 M v and integrating in v shows that div x 1 A(v)g + ∇ x 1 1 3 |v| 2 g = -∇Φ g + f s (Z e Φ + g ) → -ρ∇Φ(x) + f s
in the sense of distributions as → 0, with the notation

A(v) := v ⊗2 -1 3 |v| 2 .
Using (2.12), we recast the limit above as

div x 1 A(v)g + ∇ x 1 1 3 |v| 2 g → (θ -θ)∇Φ + f s . (2.13)
The elementary properties of the tensor field A used in the paper are recalled in the Appendix. In particular, elementary computations and Hilbert's Lemma 1 show that A ∈ (Ker L) ⊥ = Ran(L) , so that there exists a unique tensor field denoted  such that L  = A and Â⊥ Ker L .

(See Lemma 3.) Thus

1 A(v)g = Â(v) 1 Lg .
Returning to the scaled Boltzmann equation, we observe that

Z e Φ(x) 1 Lg =Q(g , g ) -v • ∇ x g + M -1 ∇ x Φ • ∇ v (M g ) -M -1 f s (x) • ∇ v (M (Z e Φ(x) + g )) →Q(g, g) -v • ∇ x g as → 0.
At this point, we recall the following useful result.

Lemma 2. For each φ, ψ ∈ Ker L, one has

Q(φ, ψ) = 1 2 L(φψ) .
See formula (60) in [START_REF] Bardos | Fluid dynamic limits of the Boltzmann equation I[END_REF] (one should notice the slightly different definitions of L and Q in formula (20) of [START_REF] Bardos | Fluid dynamic limits of the Boltzmann equation I[END_REF], which account for the different sign and normalizing factor 1 2 ). Since g(x, •) ∈ Ker L, Lemma 2 implies that Q(g, g) = 1 2 L(g 2 ), so that

Z e Φ(x) 1 Lg → 1 2 L(g 2 ) -v • ∇ x g .
On the other hand, (2.9) and (2.12) imply that

g = θ + u • v + θ 1 2 (|v| 2 -5) (2.14) so that g 2 =A(u) : A(v) + 2θu • B(v) + 1 3 |v| 2 |u| 2 + 1 4 θ 2 (|v| 2 -5) 2 + θ 2 + 2θu • v + θθ(|v| 2 -5) , while v • ∇ x g = A(v) : ∇ x u + B(v) • ∇ x θ ,
since u is divergence-free. Here, the notation B(v) designates the vector field

B(v) := 1 2 (|v| 2 -5)v ,
whose properties are recalled in the Appendix.

Therefore

1 A(v)g = Z -1 e -Φ(x) Â(v)Z e Φ(x) 1 Lg → ( 1 2 A(v)g 2 -Â(v)v • ∇ x g) =A(u) : 1 2 A(v) ⊗2 -∇u : Â(v) ⊗ A(v) =A(u) -ν(∇u + (∇u) T ) , on account of (2.

10) and of the identities

A ij A kl = δ ik δ jl + δ il δ jk -2 3 δ ij δ kl , Âij A kl = ν(δ ik δ jl + δ il δ jk -2 3 δ ij δ kl ) , with ν := 1 10 Â : A . (2.15) (See statement (2) in Lemma 4.) Therefore div(A(u) -ν(∇u + (∇u) T )) + ∇q = (θ -θ)∇Φ + f s , which can be recast as div(u ⊗ u) -ν∆u + ∇p = θ∇Φ(x) + f s (x)
where p = q -1 3 |u| 2 + θΦ. This is precisely the motion equation in the Navier-Stokes-Fourier system. Notice that the term

1 1 3 |v| 2 g
appearing on the left hand side of the equality (2.13) does not converge in the sense of distributions in general, but its gradient does. Define

T := lim →0 ∇ 1 1 3 |v| 2 g .
For each compactly supported, divergence free test vector field ξ ≡ ξ(x), one has

T, ξ = lim →0 R 3 ξ • ∇ 1 1 3 |v| 2 g dx = -lim →0 R 3 div ξ 1 1 3 |v| 2 g dx = 0 .
By Theorem 17' in [START_REF] De Rham | Differentiable manifolds. Forms, Currents, Harmonic Forms[END_REF], T viewed as a 1-current is homologous to 0, which means precisely that T = ∇q for some distribution q.

Step 4: the heat conduction equation

div x 1 B(v)g = -∇Φ(x) • vg + f s (x) • v(Z e Φ(x) + g ) → -u • ∇Φ(x) .
On the other hand

1 B(v)g = B(v) 1 Lg = Z -1 e -Φ(x) B(v)Z e Φ(x) 1 Lg → B(v) 1 2 L(g 2 ) -v • ∇ x g = 1 2 B(v)g 2 -B(v)v • ∇ x g = B(v) ⊗2 • uθ -B(v) ⊗ B(v) • ∇θ = 5 2 uθ -κ∇θ .
Indeed,

B i B j = 5 2 δ ij , Bi B j = κδ ij , with κ = B • B .
(2.16)

(See statement (1) in Lemma 4.) Therefore 5 2 div(uθ) -κ∆θ = -u • ∇Φ ,
which is precisely the heat conduction equation in the Navier-Stokes-Fourier system.

From Kinetic Theory to Viscous Heating

In the asymptotic Navier-Stokes-Fourier regime discussed above, the fluctuations of velocity field and of temperature are small and of the same order O( ). In this case, the fluctuation of kinetic energy is negligible when compared to the fluctuation of internal energy. In order to keep both fluctuations small and of the same order, it is natural to scale the fluctuation of velocity field as O( ), while the fluctuation of temperature should be of order O( 2 ).

At the level of the Boltzmann equation, this scaling assumption is obtained by choosing the distribution function of the form

F (x, v) = M (1 + g (x, v) + 2 h (x, v)) (2.17)
where

g (x, v) = -g (x, -v) , while h (x, v) = h (x, -v) (2.18)
for a.e. (x, v). Here, we assume that there is no conservative force, i.e. we take the potential Φ identically 0. The total external force acting on the gas is therefore f s ≡ f s (x) such that div f s = 0.

The dimensionless Boltzmann equation (1.1) is scaled as follows:

Ma = Kn = , and Fr = 1 .

(2. [START_REF] Glassey | The Cauchy Problem in Kinetic Theory[END_REF] In other words, the scaled Boltzmann equation takes the form

v • ∇ x F + 2 f s (x) • ∇ v F = 1 C(F ) . (2.20)
The exposition in this section follows closely [START_REF] Bardos | Kinetic equations : fluid dynamical limits and viscous heating[END_REF], where the idea of the even-odd decomposition of the distribution function seems to have been used for the first time.

First, we express the local balance laws of mass, momentum and energy in terms of the fluctuations g and h . The odd contributions of either g or h vanish after integration in v, so that

R 3 F dv = 1 + 2 h , R 3 vF dv = vg , while R 3 v ⊗ vF dv = I + 2 v ⊗ vh , R 3 v|v| 2 F dv = v|v| 2 g .
Hence, the local balance laws of mass, momentum and energy implied by the Boltzmann equation take the form

div x vg = 0 , (mass)
while

div x v ⊗ vh = f s + 2 f s h , (momentum)
and

div x v 1 2 |v| 2 g = 2 f s • vg . (energy)
Likewise, both sides of the Boltzmann equation are decomposed into even and odd components, observing that, for each rapidly decaying F

C(F • R) = C(F ) • R , for all R ∈ O 3 (R) .
Since the Maxwellian M is a radial function, one has M • R = M , and therefore

L(φ • R) = (Lφ) • R , Q(φ • R, ψ • R) = Q(φ, ψ) • R
for all R ∈ O 3 (R) and all rapidly decaying φ, ψ. These identities are satisfied in particular for R = -I. Hence the even and odd components of C(F ) are respectively

C(F ) even = -2 M Lh + 2 M Q(g , g ) + 4 M Q(h , h ) , C(F ) odd = -M Lg + 2 3 M Q(g , h ) .
Therefore, the scaled Boltzmann equation is equivalent to the system

       v • ∇ x g + 2 M -1 f s • ∇ v (M g ) = -Lh + Q(g , g ) + 2 Q(h , h ) , 2 v • ∇ x h -2 v • f s + 4 M -1 f s • ∇ v (M h ) = -Lg + 2 2 Q(g , h ) . (2.21)
Theorem 2. Let F be a family of solutions of the scaled Boltzmann equation (2.20), whose relative fluctuations g , h defined in (2.17)-(2.18) satisfy

g → g , h → h in L 2 (Ω × R 3 ; (1 + |v| 2 )M dvdx) weak, and
Q(g , g ) → Q(g, g) , Q(g , h ) → Q(g, h) in L 1 (Ω; L 2 (R 3 ; (1 + |v| 2 )M dv)) weak, while 2 h ∇φ(v) → 0 in L 2 (Ω) weak for each φ ∈ L 2 (R 3 ; (1 + |v| 2 )M dv). Then g(x, v) = u(x) • v , (2.22) while h(x, v) = 1 2 A(u(x)) : A(v) -Â(v) : ∇ x u(x) + ρ(x) + (θ(x) + 1 3 |u(x)| 2 ) 1 2 (|v| 2 -3) (2.23)
where (u, θ) is a solution of the incompressible Navier-Stokes-Fourier system with viscous heating

                 div u = 0 , div(u ⊗2 ) + ∇p = ν∆u + f s , 5 2 div(uθ) -u • ∇ x p = κ∆θ + 1 2 ν ∇ x u + (∇ x u) T 2 .
(2.24)

The values of the viscosity ν and heat diffusivity κ are determined implicitly in terms of the collision integral, by formulas (2.15) and (2.16), as in Theorem 1.

Proof. The proof follows more or less the same lines as that of Theorem 1; see also [START_REF] Bardos | Kinetic equations : fluid dynamical limits and viscous heating[END_REF].

Step 1: asymptotic form of g and divergence-free condition

We first deduce from the second equation in (2.21) and the assumption that g → g in the sense of distributions that Lg = 0 .

Hence g(x, •) ∈ Ker L and v → g(x, •) is odd for a.e. x, so that g is of the form (2.22). The local conservation of mass implies that

0 = div x vg → div x vg = div x u
in the sense of distributions, so that div x u = 0 .

(2.25)

Step 2: asymptotic form of h and divergence-free condition

Next we deduce from the first equation in (2.21) that

h → h with v • ∇ x g = -Lh + Q(g, g) .
Since g(x, •) ∈ Ker L, applying Lemma 2 shows that

L(h -1 2 g 2 ) = -v • ∇ x g = -v ⊗ v : ∇ x u = -A(v) : ∇ x u ,
since u is divergence free by (2.25). Hence

h(x, v) = 1 2 g 2 (x, v) -Â(v) : ∇ x u + h 0 (x, v) , with h 0 (x, •) ∈ Ker L for a.e. x. Since h 0 is even in v, it is of the form h 0 (x, v) = (x) + θ(x) 1 2 (|v| 2 -3) . (2.26)
On the other hand

g 2 (x, v) = (u(x) • v) 2 = u(x) ⊗ u(x) : v ⊗ v = A(u(x)) : A(v) + 1 3 |u(x)| 2 |v| 2 , so that h is of the form h(x, v) = 1 2 A(u(x)) : A(v) -Â(v) : ∇ x u(x) + ρ(x) + (θ(x) + 1 3 |u(x)| 2 ) 1 2 (|v| 2 -3) , (2.27
)

with ρ = + 1 6 |u| 2 .
Step 3: motion equation Passing to the limit in the sense of distributions in the local balance law of momentum shows that

div x v ⊗ vh = f s .
We insert the expression for h found above in this identity. Observe that

v ⊗ vh = A(v)( 1 2 A(u(x)) : A(v) -Â(v) : ∇ x u(x)) + 1 3 |v| 2 ρ(x) + (θ(x) + 1 3 |u(x)| 2 ) 1 2 (|v| 2 -3) I because A and  ⊥ Ker L. By statements (1)-(2) in Lemma 4 A(v)( 1 2 A(u(x)) : A(v) -Â(v) : ∇ x u(x)) = A(u(x)) -ν(∇ x u + (∇ x u) T ) , while 1 3 |v| 2 ρ(x) + (θ(x) + 1 3 |u(x)| 2 ) 1 2 (|v| 2 -3) = ρ(x) + θ(x) + 1 3 |u(x)| 2 .
Hence

div x A(u) -ν(∇ x u + (∇ x u) T ) + ∇ x (ρ + θ + 1 3 |u| 2 ) = f s , or equivalently div x (u ⊗ u) -ν∆ x u + ∇ x p = f s (2.28)
with p = ρ + θ.

Step 4: heat equation Finally, we combine the local balance of mass and energy so that

div x 1 2 Bg = f s • vg → f s • u in the sense of distributions.
Next, we transform the left hand side of the equality above by exactly the same method as in the previous section:

1 2 Bg = B 1 2 Lg = B(2Q(g , h ) -v • ∇ x h -2 M -1 f s • ∇ v (M h )) → B(2Q(g, h) -v • ∇ x h) .
Notice that

R 3 Bf s • ∇ v M dv = -Bf s • v = 0 , because B ⊥ Ker L. Hence div x B(2Q(g, h) -v • ∇ x h) = f s • u ,
and we insert in this last expression the explicit formulas for g and h.

First

2Q(g, h) =2Q(u • v, 1 2 A(u) : A(v) + ρ + (θ + 1 3 |u| 2 ) 1 2 (|v| 2 -3)) -2Q(u • v, Â(v) : ∇ x u) .
(2.29)

Applying formula (2.18c) in [START_REF] Bardos | Kinetic equations : fluid dynamical limits and viscous heating[END_REF] shows that

2Q(u • v, 1 2 A(u) : A(v) + ρ + (θ + 1 3 |u| 2 ) 1 2 (|v| 2 -3)) = L(θu • B(v) + 1 3 C(v) : u ⊗ u ⊗ u) ,
where

C(v) := 1 2 (v ⊗ v ⊗ v -3v ⊗ I) . Hence 2 B(v)Q(u • v, 1 2 A(u) : A(v) + ρ + (θ + 1 3 |u| 2 ) 1 2 (|v| 2 -3)) = B(v)L(θu • B(v) + 1 3 C(v) : u ⊗ u ⊗ u) = B(v)(θu • B(v) + 1 3 C(v) : u ⊗ u ⊗ u) = 1 2 (|u| 2 + 5θ)u (2.30)
by Lemma 5 and statement (1) in Lemma 4.

By Proposition 2.6 in [START_REF] Bardos | Kinetic equations : fluid dynamical limits and viscous heating[END_REF] and statements (1) and (3) of Lemma 4,

2 B(v)Q(u • v, Â(v) : ∇ x u) = (A(v) • u) Â(v) : ∇ x u -B(v)(u • v)A(v) : ∇ x u =(ν -2 5 κ)(∇ x u + (∇ x u) T ) • u .
(2.31)

Finally, we compute

B(v) ⊗ vh = 1 2 B(v) ⊗ v ⊗ A(v) : A(u) -B(v) ⊗ v ⊗ Â(v) : ∇ x u + (θ + 1 3 |u| 2 ) B(v) ⊗ B(v) = 2 5 κA(u) -c(∇ x u + (∇ x u) T ) + κ(θ + 1 3 |u| 2 )I (2.32)
according to statements (1) and (3) of Lemma 4. Notice that

ρ B ⊗ v = 0 ,
and that

B ⊗ v 1 2 (|v| 2 -3) = B ⊗ B , because B ⊥ Ker L.
Therefore, putting together (2.30)-(2.31)-(2.32), we arrive at the identity

div x B(2Q(g, h) -v • ∇ x h) = div x ( 1 2 (|u| 2 + 5θ)u -(ν -2 5 κ)(∇ x u + (∇ x u) T ) • u) -∇ x ⊗ ∇ x : ( 2 5 κA(u) -c(∇ x u + (∇ x u) T ) + κ(θ + 1 3 |u| 2 )I) .
(2.33)

This identity can be substantially simplified, as follows. First,

∇ x ⊗ ∇ x : ∇ x u + (∇ x u) T = 0 .
because of the divergence-free condition (2.25). On the other hand,

div x (∇ x u + (∇ x u) T ) • u) = ∇ x ⊗ ∇ x : u ⊗ u + ∆ x 1 2 |u| 2 , while ∇ x ⊗ ∇ x : (A(u) + 5 6 |u| 2 I) = ∇ x ⊗ ∇ x : u ⊗ u -1 3 ∆ x |u| 2 + 5 6 ∆ x |u| 2 = ∇ x ⊗ ∇ x : u ⊗ u + ∆ x 1 2 |u| 2 . Hence div x B(2Q(g, h) -v • ∇ x h) = div x 1 2 (|u| 2 + 5θ)u -ν ∇ x u + (∇ x u) T • u -κ∆ x θ ,
so that the limiting form of the local energy balance is

div x 1 2 (|u| 2 + 5θ)u = ν div x ((∇ x u + (∇ x u) T ) • u) + κ∆ x θ + f s • u .
On the other hand, multiplying both sides of the Navier-Stokes motion equation by u, we arrive at the identity

div x u 1 2 |u| 2 + u • ∇ x p = f s • u + νu • ∆ x u = f s • u + νu • div x ∇ x u + (∇ x u) T .
We further simplify the right hand side of the equality above as follows:

u • div x ∇ x u + (∇ x u) T = div x ∇ x u + (∇ x u) T u -∇ x u + (∇ x u) T : ∇ x u = div x ∇ x u + (∇ x u) T u -1 2 ∇ x u + (∇ x u) T 2 .
Combining these identities with (2.33) above leads to

5 2 div x (uθ) -u • ∇ x p = κ∆ x θ + 1 2 ν ∇ x u + (∇ x u) T 2 .
(2.34)

Observe that the Navier-Stokes motion equation is exactly the same when Φ = 0 for both scaling assumptions (2.1) and (2.19). The temperature equation, however, is very different according to whether the Froude number is O( √ ) as in (2.1), or O(1) as in (2.19).

The viscous heating term on the right hand side of the equation governing the temperature field appears in Sone's asymptotic analysis of the hydrodynamic limits of the Boltzmann equation in the weakly nonlinear regime -see [START_REF] Sone | Molecular Gas Dynamics: Theory, Techniques and Applications[END_REF]. Sone's original work [START_REF] Sone | Asymptotic theory of flow of rarefied gas over a smooth boundary II[END_REF] on the weakly nonlinear hydrodynamic limits of kinetic theory was written in the case of the BGK model; see [START_REF] Sone | Steady gas flows past bodies at small Knudsen numbers -Boltzmann and hydrodynamic systems[END_REF] for the extension to the Boltzmann equation. Sone's argument is based on the Hilbert expansion. One should pay attention to a particular feature of Sone's theory: the pressure field p in (2.8) and the velocity and temperature fields u and θ in (2.8) do not appear at the same order in Sone's expansion; in fact p appears at order O( 2 ) in the expansion of the distribution function in powers of , while u and θ appear at order O( ) in that same expansion. See formulas (3.77), (3.79b-c), (3.80d) and (3.88a-c) in section 3.2.2 of [START_REF] Sone | Molecular Gas Dynamics: Theory, Techniques and Applications[END_REF]. The limit leading to (2.24) corresponds to a situation where u appears at order O( ) in Sone's expansion, while the leading order temperature fluctuation (denoted τ S1 in Sone's analysis, see formula (3.79c) in [START_REF] Sone | Molecular Gas Dynamics: Theory, Techniques and Applications[END_REF]) is identically zero. The temperature fluctuation appears at order O( 2 ) in Sone's expansion, together with the pressure field p, and the temperature equation in (2.24) coincides with formula (3.89c) in section 3.2.2 of [START_REF] Sone | Molecular Gas Dynamics: Theory, Techniques and Applications[END_REF]. (Sone's analysis in [START_REF] Sone | Molecular Gas Dynamics: Theory, Techniques and Applications[END_REF] does not involve an external force, but the work of the external force disappears from the temperature equation when combining the motion equation and the energy equation as explained above.)

Boundary Conditions

The discussion of the incompressible Navier-Stokes-Fourier limit of the Boltzmann equation presented above would remain incomplete without discussing the boundary condition. In this section, we briefly describe the simplest imaginable situation.

Assume that the scaled Boltzmann equation (2.2) is supplemented with a diffuse reflection condition at the boundary of the spatial domain Ω on which the Boltzmann equation (1.1) is posed. In other words, for each x ∈ ∂Ω, one has

F (x, v) = √ 2πM R 3 F (x, v)(v • n x ) + dv , v • n x < 0 , (2.35)
where x → n x is the unit normal field defined on the boundary ∂Ω of the spatial domain. Here, we have assumed for simplicity that there is no temperature gradient on ∂Ω. The constant temperature at the boundary (i.e. the temprature 1 in the Maxwellian state M ) defines the scale of the speed of sound in the interior of the domain.

By construction

R 3 F (x, v)v • n x dv = 0 , x ∈ ∂Ω ,
which means that the net mass flux at each point x ∈ ∂Ω is identically 0. This suggests that the boundary condition (2.35) should be supplemented with the additional condition

Ω×R 3 F (x, v)dxdv = |Ω| , (2.36)
that is consistent to leading order with the normalization of Z in (2.3)-(2.17), and is equivalent to the condition (2.4) already introduced above in the case Ma = Kn = Fr 2 = , and

Ω h dx = 0
in the case Ma = Kn = Fr = . (Notice that, in the latter case, g = 0 a.e. on Ω since g is odd in v.)

Besides, we assume that the force f s satisfies both

div x f s = 0 on Ω , f s • n x = 0 on ∂Ω . (2.37) Define Λ x (φ) := √ 2π φ(v • n x ) + . Theorem 3. Let Ma = Kn = Fr 2 = ,
∂Ω × R 3 satisfies g ∂Ω×R 3 → g ∂Ω×R 3 weakly in L 2 (∂Ω × R 3 ; |v • n x |M dvdS(x))
where g is such that g → g weakly in L 2 (Ω × R 3 , M dvdx). Then

g(x, v) = θ + u(x) • v + θ(x) 1 2 (|v| 2 -5)
where

θ = 1 |Ω| Ω θ(x)dx ,
and (u, θ) is a solution of the Navier-Stokes-Fourier system (2.8), with Dirichlet boundary condition

u ∂Ω = 0 , θ ∂Ω = 0 .
Proof. Indeed, the diffuse reflection condition implies that

g (x, v) = Λ x (g (x, •)) , v • n x < 0 . (2.38)
Thus one can pass to the limit as → 0 in (2.38). One arrives at

g(x, v) = Λ x (g(x, •)) , x ∈ ∂Ω , v ∈ R 3 .
Since we already know from Theorem 1 that g of the form

g(x, v) = θ + u(x) • v + θ(x)
u ∂Ω = 0 , θ ∂Ω = 0 .
Proof. Specializing the equality above to the case where v is tangential to the boundary, we find that

g (x, v) = Λ x ((g + h )(x, •)) -h (x, v) , x ∈ ∂Ω , v • n x = 0 ,
and observe that the left hand side of this equality is odd in v = v-(v•n x )n x , while the right hand side is even. Therefore both sides vanish, so that

g (x, v) = 0 , x ∈ ∂Ω , v • n x = 0 . while h (x, v) = Λ x 1 g + h (x, •) , x ∈ ∂Ω , v • n x = 0 .
Passing to the limit on both sides of the first equality as → 0 shows that

g(x, v) = Λ x (g(x, •)) , x ∈ ∂Ω , v • n x = 0 ,
and we conclude from (2.22) that u ∂Ω = 0 .

Next we consider the differential operator

P (x, D v ) := (I -n x ⊗ n x )∇ v
-which is the orthogonal projection of ∇ v on the tangential direction of ∂Ω at x, and observe that

P (x, D v )h (x, v) = 0 , x ∈ ∂Ω , v • n x = 0 .
Passing to the limit in both sides of this identity, we conclude that

P (x, D v )h(x, v) = 0 , x ∈ ∂Ω , v • n x = 0 .
Substituting the expression (2.27) in this equality, we find that

P (x, D v )h(x, v) = -P (x, D v ) Â(v) : ∇u(x) + v τ θ(x) , x ∈ ∂Ω , v • n x = 0 with v τ = v -(v • n x )n x .
Observe that

P (x, D v ) Â(v) : ∇u(x) = v τ • ∇u(x) + (I -n x ⊗ n x )∇(u(x) • v τ ) = 0 ,
since we already know that u ∂Ω = 0 and all the derivatives of u appearing in the expression above are taken in directions tangential to ∂Ω. Hence

v τ θ(x) = 0 , for each x ∈ ∂Ω , v • n x = 0 ,
which implies that θ ∂Ω = 0.

Spatially Periodic Steady Solutions

The Navier-Stokes Equations

In this section, we assume that the spatial domain is Ω = T 3 . Consider the system (2.8) posed on T 3 , and seek solutions satisfying

T 3 u(x)dx = T 3 θ(x)dx = 0 . (3.1)
For simplicity, we assume further that Φ ≡ 0.

Multiplying both sides of the last equation in (2.8) by θ, we see that

u • ∇( 1 2 θ 2 ) = div(u 1 2 θ 2 ) = κθ∆θ , so that κ T 3 |∇θ(x)| 2 dx = 0 . Therefore θ ≡ 0.
Conversely, if u is a solution of the motion equation and Φ ≡ 0, then (u, 0) is a solution of the Navier-Stokes-Fourier system without viscous heating.

There exist indeed nontrivial solutions of the motion equation with nonzero solenoidal external force f s . The simplest example is the case of a shear flow

u(x) = (0, 0, U (x 1 , x 2 )) , f s (x) = (0, 0, a(x 1 , x 2 ))
with -ν∆U

T 2 U (x 1 , x 2 )dx 1 dx 2 = T 2 a(x 1 , x 2 )dx 1 dx 2 = 0 . Obviously div u = ∂ x 3 U = 0 , div f s = ∂ x 3 a = 0 , and div(u ⊗2 ) = (u • ∇)u = U ∂ x 3 u = 0 ,
(x 1 , x 2 ) = a(x 1 , x 2 ) , (x 1 , x 2 ) ∈ T 2 .
For each zero-mean a ∈ L 2 (T 3 ), there exists a unique zero-mean U ∈ H 2 (T 2 ) of the Poisson equation above.

More generally, the following result is classical.

Theorem 5. For each f s ∈ L 2 (T 3 ) satisfying

T 3 f s (x)dx = 0 ,
there exists at least one solution u ∈ H 2 (T 3 ) of the Navier-Stokes equations with external force f s such that

T 3 u(x)dx = 0 .
Besides, there exists c > 0 such that the solution u ∈ H 2 (T 3 ) and is unique

if f s L 2 ≤ cν 2 .
The proof of this classical result is given below -see chapter II, §1 in [START_REF] Temam | Navier-Stokes equations[END_REF] for a similar result in a slightly different (nonperiodic) setting. First, we recall some elements of notation. We denote by H the subspace of L 2 (T 3 ; R 3 ) of vector fields v such that

T 3 v(x)dx = 0 ,
and we set H 1 = H 1 (T 3 ) ∩ H. It will be convenient to use the norm

v H 1 := ∇v L 2 .
We denote by Π the L 2 -orthogonal projection on divergence-free vector fields. In other words, if v ∈ H, its Fourier decomposition is

v(x) = k∈Z 3 \{0} v k e i2πk•x , and
Πv(x) = k∈Z 3 \{0} (v k -|k| -2 k(k • v k ))e i2πk•x .
Likewise, for each zero-mean φ ∈ H we define

(-∆) -s φ(x) := k∈Z 3 \{0} |2πk| -2s φ k e i2πk•x where φ k = T 3 φ(x)e -i2πk•x dx .
Proof. Consider the map

T : v → (-∆) -1 Π div(v ⊗2 ) .
First, observe that T maps H 1 into itself:

T v H 1 = (-∆) -1/2 Π div(v ⊗2 ) L 2 = Π(-∆) -1/2 div(v ⊗2 ) L 2 ≤ (-∆) -1/2 div(v ⊗2 ) L 2 ≤ v ⊗2 L 2 ≤ v 2 L 4 ≤ C 2 v 2 H 1 by Sobolev embedding (H 1 ⊂ L p (T 3 ) for 1 ≤ p ≤ 6). Similarly T u -T v H 1 = (-∆) -1/2 Π div(u ⊗2 -v ⊗2 ) L 2 ≤ u ⊗2 -v ⊗2 L 2 ≤ u ⊗ (u -v) L 2 + (u -v) ⊗ v L 2 ≤ ( u L 4 + v L 4 ) u -v L 4 ≤ C 2 ( u H 1 + v H 1 ) u -v H 1 .
Hence T is continuous from H 1 into itself, and Lipschitz continuous on balls of H 1 .

The Navier-Stokes equations can be put in the form

u = (-ν∆) -1 f s -1 ν T u ,
and is embedded into the family of equations

F α (u) = (-ν∆) -1 f s , where F α = I + α ν T parametrized by α ∈ [0, 1]. If u, v ∈ H 1 with u H 1 ≤ R and v H 1 ≤ R, then the map G α defined by G α (w) = (-ν∆) -1 f s + w -F α (w) so that ∇u L 2 ≤ 1 ν (-∆) -1/2 f s L 2 . Setting R = 1 ν (-∆) -1/2 f s L 2 + 1, we see that F α maps B(0, R ) ⊂ H 1 into H 1 and that u ∈ ∂B(0, R ) ⊂ H 1 implies that F α (u) = (-ν∆) -1 f s for all α ∈ [0, 1]. Therefore degree(F 1 , B(0, R ), (-ν∆) -1 f s ) = degree(F 0 , B(0, R ), (-ν∆) -1 f s ) = 1
because F 0 is the identity. Hence the equation

F 1 (u) = (-ν∆) -1 f s ,
which is equivalent to the steady Navier-Stokes equations in T 3 , has at least one solution in H 1 .

However, the only solution of the Navier-Stokes-Fourier system with viscous heating (2.24) satisfying (3.1) is (0, 0). Indeed, integrating both sides of the last equation in (2.24) shows that

T 3 |∇u + (∇u) T | 2 (x)dx = 0 . Hence ∇u + (∇u) T = 0 .
In particular, div u = 0 and div(∇u) = -∆u = div((∇u) T ) = ∇(div u) = 0 , so that u is a harmonic vector field on T 3 satisfying (3.1). Hence u = 0. Returning to the heat equation in (2.24), we see that θ is a harmonic function on T 3 , and (3.1) implies that θ = 0.

The Boltzmann Equation

Theorem 6. Let F ≡ F (x, v) be a C 1 solution of the steady Boltzmann equation

v • ∇ x F + f (x) • ∇ v F = C(F ) , x ∈ T 3 , v ∈ R 3 ,
and assume that F is rapidly decaying in v while ln F has polynomial growth in v as |v| → ∞.

Then f is a gradient field, while F is a Maxwellian distribution with constant temperature. More precisely, there exists Φ ∈ C 1 (T 3 ), a vector u ∈ R 3 , and two constants θ > 0 and C ≥ 0 such that f = -∇Φ , and

F (x, v) = C exp - 1 θ 1 2 |v -u| 2 + Φ(x) with u • ∇Φ = 0 .
In particular, if div f = 0, then f = 0 and F is a uniform Maxwellian.

Proof. The global form of Boltzmann's H theorem shows that

T 3 ×R 3 C(F ) ln F dxdv = 0 . Hence F (x, v) is a local Maxwellian satisfying v • ∇ x F + f (x) • ∇ v F = C(F ) = 0 , x ∈ T 3 , v ∈ R 3 . Setting F (x, v) = M (ρ(x),u(x),θ(x)) (v), the Boltzmann equation reduces to v • ∇ρ(x) ρ(x) + (v -u) ⊗ v : ∇u(x) θ(x) + v • ∇θ(x) 2θ(x) 2 (|v -u(x)| 2 -3θ(x)) = f (x) • v -u(x) θ(x) .
Setting

V := v -u(x) θ(x) ,
the equality above is recast as

u(x) • ∇ ln ρ(x) + θ(x) -1/2 V • (θ(x)∇ ln ρ(x) + (u • ∇)u) + 1 2 (|V | 2 -3)u(x) • ∇ ln θ(x) + V ⊗2 : ∇u(x) +(|V | 2 -3)V • ∇ θ(x) = θ(x) -1/2 f (x) • V Hence ∇ √ θ = 0 so that θ = Const.
and the equality above reduces to

u(x) • ∇ ln ρ(x) + θ -1/2 V • (θ∇ ln ρ(x) + (u • ∇)u) + V ⊗2 : ∇u(x) = θ(x) -1/2 f (x) • V . Therefore u • ∇ρ = 0 , θ∇ ln ρ + (u • ∇)u = f , ∇u + (∇u) T = 0 .
Using the third equation, the second equation is recast as

θ∇ ln ρ + (u • ∇)u = θ∇ ln ρ -(∇u) T u = θ∇ ln ρ -∇ 1 2 |u| 2 = f so that f = ∇(θ ln ρ -1 2 |u| 2 )
must be a gradient field.

Next ∇u + (∇u) T = 0 which implies as above that u is harmonic on T 3 , and therefore is a constant.

We conclude that

f (x) = -∇Φ(x) , F (x, v) = Ce -1 θ ( 1 2 |v-u| 2 +Φ(x))
with u • ∇Φ(x) = 0 .

Hence u = 0 unless Φ takes its values in an affine space orthogonal to u.

Physical discussion and study of a numerical example

As we observed in Sec. 3.2, there is no spatially periodic steady solution of the Boltzmann equation for an external force that is not derived from a potential. In contrast, the Navier-Stokes equations have a spatially periodic and steady solution for such an external force. This discrepancy seems to be contradictory, since the Navier-Stokes equations are derived from the Boltzmann equation, as shown in Sec. 2.1. In this section, we will further examine this seemingly contradictory results.

A possible physical explanation of the paradox

Let us consider a gas in a periodic box or in a specularly reflecting box. As we have seen in Sec. 3.2, the steady solution of the Boltzmann equation with an external force derived from a potential has the following properties: the temperature is uniform, the macroscopic flow of the gas vanishes except for a special case, and the density is distributed according to the potential (i.e., a stratified gas at rest). Now, let us consider a time-dependent problem starting from a given initial state. If the external force has a potential, the time-dependent solution should approach the steady solution mentioned above, i.e., the solution without a gas flow and with a density stratification, in the long-time limit. Then, what will happen when the force does not have a potential? At the initial stage, a gas flow is caused by the external force. But, since the force does not have a potential, the density stratification that blocks the gas flow cannot be formed. Therefore, the flow remains forever, or at least, for much longer time. The induced gas flow, in general, has a shear. If we consider the case with relatively small Knudsen number, this shear gives rise to the viscous heating. However, because the boundary is periodic or adiabatic (in the case of specularly reflecting box), the heat generated by the shear in the gas cannot escape through the boundary. This means that the temperature in the gas will increase indefinitely. This is the reason why there is no steady solution for the Boltzmann equation.

On the other hand, the Navier-Stokes equations (2.8) are derived in the limit where the effect of the viscous heating is negligibly small. Therefore, the mass and momentum equations (the first two equations in (2.8)) are decoupled with the energy equation (the last equation) when Φ = 0. The former equations give a velocity field with a shear, and the latter equation gives a constant temperature field. However, if we include the effect of viscous heating, as was done in Sec. 2.2, the energy equation is changed to the last equation in (2.24), the mass and momentum equations being unchanged. Therefore, the flow velocity field is the same, but the new energy equation does not have a solution for this velocity field. If we consider a time-dependent version of this energy equation, we can easily see that the temperature increases indefinitely. In conclusion, the difference between the Boltzmann equation and the Navier-Stokes equations is due to the fact that the effect of viscous heating is neglected in the Navier-Stokes limit with Ma = Kn = Fr 2 = .

Numerical example

For the purpose of understanding the phenomenon of heating predicted above, we consider a simple numerical example. Let us consider a gas in a two-dimensional square box -1/2 < x 1 < 1/2, -1/2 < x 2 < 1/2 with periodic condition on each side. We assume that the external force is of the form f = (0, f 0 sin 2πx 1 , 0), which is divergence free and does not have a potential. Initially, the gas is in a uniform equilibrium state at rest with density 1 and temperature 1. We pursue the time evolution of the solution and observe whether the temperature increases indefinitely or not.

We analyze the problem mainly using the Bhatnagar-Gross-Krook (also known as BGK) model, but some results based on the original Boltzmann equation will also be presented. Because of the form of the external force, we can assume that the flow field is spatially one dimensional depending only on x 1 and is periodic in x 1 with period 1. Moreover, the external force is symmetric with respect to x 1 = -1/4 and x 1 = 1/4, so that we can also assume the same symmetry for the flow field. Therefore, placing specularly reflecting boundaries at x 1 = -1/4 and x 1 = 1/4, we can analyze the problem in the finite interval -1/4 ≤ x 1 ≤ 1/4.

Here, we formulate the problem using the BGK model. In the present problem, the BGK model in an appropriate dimensionless form can be written as

∂ t F + v 1 ∂ x 1 F + f 0 sin 2πx 1 ∂ v 2 F = 1 Kn C BGK (F ),
where F (x 1 , t, v) is the velocity distribution function including the time variable t, and

C BGK (F ) = 8 π 1/2 ρ (M -F ) , M = M (ρ(x 1 ,t),v(x 1 ,t),θ(x 1 ,t)) = ρ(x 1 , t) [2πθ(x 1 , t)] 3/2 exp - |v -u(x 1 , t)| 2 2θ(x 1 , t) , ρ(x 1 , t) = R 3 F dv, u i (x 1 , t) = 1 ρ R 3 v i F dv, (i = 1, 2, 3; u 3 = 0) θ(x 1 , t) = 1 3ρ R 3 |v -u| 2 F dv.
Here, Kn is the Knudsen number, i.e., the mean free path in the initial equilibrium state at rest divided by the length of the period (note that Kn is denoted by in Sec. 2 because the limit as → 0 is discussed there).

The factor (8/π) 1/2 appears because of the manner of nondimensionalization used here in consistency with the form of the Boltzmann equation in earlier sections. The specular reflection condition at x 1 = -1/4 and 1/4 is as follows:

F (±1/4, t, v) = F (±1/4, t, Rv), for ∓ v 1 > 0,
where R is the reflection operator: Rv = (-v 1 , v 2 , v 3 ). The initial condition is given by

F (x 1 , 0, v) = M (1,0,1) .
We solve this initial-boundary-value problem by the finite-difference method.

Now we show some of the numerical results for Kn = 0.1 and for f 0 = 2 (Fig. Finally, we present some results based on the Boltzmann equation. Let us replace the right-hand side, (1/Kn)C BGK (F ), of the basic equation with (1/Kn)C(F ), where C is the dimensionless Boltzmann collision operator for hard-sphere molecules, as in Sec. 1. The mean free path used here to define Kn is ( √ 2πd 2 n 0 ) -1 with d the diameter of a molecule and n 0 the reference molecular number density, which is related to 0 in Sec. 1 [cf. (1.1)] as (2 √ 2π) -1 0 . We employ the direct simulation Monte Carlo (DSMC) method as the solution method. The method is a particle and stochastic one, so that it contains the inherent statistical fluctuations. For steady problems, we can take the time average over a long interval of time to reduce the fluctuations. But the time-averaging does not work for time-dependent problems as the present one. The only possible way to reduce them is to perform many independent runs and take an ensemble average over the runs. The method is also not appropriate for describing small quantities because they are hidden in the fluctuations. In fact, it is impossible to obtain θθ av , ρ, and u 1 . Nevertheless, we show in Figs. 2) because the longer computation is practically impossible. However, the convergence of α seems to be faster for hard-sphere molecules. It is likely that α tends to approach about 0.665 in Fig. 4.3(d) and about 0.66 in Fig. 4.4(c). Therefore, irrespective of the model of the collision term, we have the asymptotic behavior like θ av ≈ C θ t 0.66 for large t.

Interpretation of the numerical results

Next, in order to understand the slow increase of the temperature and the slow decrease of the flow speed, we try a rough discussion on the basis of the compressible Navier-Stokes equations. The numerical results show that the flow is almost unidirectional, i.e., u ≈ (0, u 2 , 0). Therefore, let us assume that the flow is unidirectional, i.e., u = (0, u 2 , 0) and the problem in spatially one-dimensional (∂/∂x 2 = ∂/∂x 3 = 0). Then, the compressible Navier-Stokes equations reduce to the following equations: ∂p ∂x 1 = 0,

∂u 2 ∂t = 1 ρ ∂ ∂x 1 µ(θ) ∂u 2 ∂x 1 + g 0 sin 2πx 1 , ∂θ ∂t = 2 3 1 ρ ∂ ∂x 1 κ(θ) ∂θ ∂x 1 + 2 3 µ(θ) ρ ∂u 2 ∂x 1 2
, where a suitable nondimensionalization has been made, and g 0 is a constant.

In addition, µ(θ) and κ(θ) are, respectively, dimensionless forms of the vis- Suppose that µ(θ av ) = C µ θ δ av . Then, for the initial condition θ = 1 at t = 0, we obtain θ 1+δ av ≈ C 1 t + 1 with C 1 = (1 + δ)g 2 0 ρ/12π 2 C µ , or for large t, θ av ≈ (C 1 t + 1)

1 1+δ ≈ C 1 1+δ 1 t 1 1+δ .
Since ∂θ/∂t = O(∂θ av /∂t) = O(t -δ/(1+δ) ) and µ(∂u 2 /∂x 1 ) 2 = O(1/µ(θ av )) = O(t -δ/(1+δ) ), they are small in the energy equation. If we neglect these terms in the energy equation, we obtain ∂θ/∂x 1 ≈ 0, i.e., θ ≈ θ av . Then, from the expression ∂u 2 /∂x 1 ≈ (g 0 ρ/2πµ) cos 2πx 1 and the boundary condition, u 2 is obtained as u 2 ≈ (g 0 ρ/4π 2 µ) sin 2πx 1 . To summarize, we obtain θ ≈ C θ t 1 1+δ , u 2 ≈ C u t -δ 1+δ sin 2πx 1 , (for t 1), with new constants C θ and C u . Because δ = 1 for the BGK model, we have θ ≈ C θ t 1/2 and u 2 ≈ C u t -1/2 sin 2πx 1 for large t. On the other hand, since δ = 0.5 for hard-sphere molecules, θ and u 2 behave as θ ≈ C θ t 2/3 and u 2 ≈ C u t -1/3 sin 2πx 1 . Although this conclusion for hard-sphere molecules is consistent with the numerical result based on the Boltzmann equation, it does not coincide precisely with the numerical result based on the BGK model. It is natural because the argument is too sketchy. However, it provides qualitative explanation for the slow increase of the temperature and its uniformity and for the slow decrease of the flow speed and the sinusoidal flow-velocity profile. More specifically, as the result of the temperature rise caused by the viscous heating, the viscosity and the thermal conductivity increase. The high thermal conductivity leads to a uniform θ. On the other hand, the high viscosity tends to prevent the external force from causing the gas flow, so that the flow speed decreases as the temperature increases. However, the decrease of the flow speed results in the decrease of the viscous heating. The rough estimate based on the compressible Navier-Stokes equations show that, as time proceeds, the amount of viscous heating decreases, but the total amount of heat produced from the initial time increases indefinitely. Thus, the temperature continues to increase indefinitely.

Conclusion

We have discussed the numerous similarities between the steady problem for the Boltzmann equation and for the Navier-Stokes-Fourier system. In particular, the presence of the viscous heating term in the temperature equation depends on the scaling of the external divergence-free force field. However, we have observed a significant difference between the steady Navier-Stokes equation and the steady Boltzmann equation: while the steady Navier-Stokes equation with prescribed divergence-free external force always have a solution, either in the periodic setting or in the case of a bounded spatial domain with Dirichlet boundary condition, the steady Boltzmann equation with nonzero divergence-free external force field cannot have a nonzero solution. We have proposed a physical explanation for this difference, based on the long time behavior of the evolution Boltzmann equation with nonzero divergence-free external force field in the periodic setting. Our numerical computations based on the BGK model suggest that the temperature field in the gas increases indefinitely, so that the gas cannot reach a steady state.

While we have presented our results on the hydrodynamic limit of the steady Boltzmann equation in the case of a hard sphere gas, the same results should remain true in the case of hard cutoff potentials (in the sense of Grad).

Finally, our numerical simulations suggest the following problem, which we believe is open at the time of this writing.

Problem: Consider the evolution Boltzmann equation set in the periodic box, with a prescribed nonzero divergence-free external force field (for instance f (x 1 , x 2 , x 3 ) = (0, f 0 sin(2πx 1 ), 0)). What is the asymptotic behavior of the temperature field averaged in the space variable in the long time limit? In particular, does there exist α > 0 such that the average temperature field is asymptotically equivalent to Ct α for some C > 0 as the time variable t tends to infinity?

  ) are respectively the linearized collision integral at M and the collision integral intertwined with the multiplication by M . (The notation DC(M ) • (M g) designates the differential of the collision integral C evaluated at M , and applied to the variation M g of distribution function.) The quadratic operator Q defines a unique bilinear symmetric operator, also denoted Q, by the polarization formula

  Next we explain how to derive the heat conduction equation. Mutiplying both sides of the scaled Boltzmann equation by 1 M 1 2 (|v| 2 -5) and integrating in v shows that

  so that the Navier-Stokes equation reduces to the Poisson equation in T 2 :

Figure 4 . 1 : 2 0

 412 Figure 4.1: BGK model with Kn = 0.1 and f 0 = 2

Figure 4

 4 Figure 4.1: Continued

Figure 4 . 2 . 2 1/ 4 - 1 / 2 . 2 1/ 4 - 1 / 4

 4224122414 Figure 4.1(a) shows the profile of the x 2 component of the flow velocity u 2 from t = 0 to 1000/ √ 2 in the half interval 0 ≤ x 1 ≤ 1/4. Since u 1 and u 2 are odd functions of x 1 , and ρ and θ are even functions of x 1 , we show the profiles of these quantities only in the half interval here and in what follows. The sinusoidal external force parallel to the x 2 direction induces u 2 at the very early stage, but u 2 starts decreasing after t = 1.5/ √ 2. Figure 4.1(b) shows the time evolution of the average temperature θ av (t) = 2 1/4 -1/4 θ(x 1 , t)dx 1 until t = 1000/ √ 2, and Fig. 4.1(c) the corresponding evolution of the profile of the deviation θθ av . The θ becomes nonuniform at the early stage but tends to become uniform after t = 1.5/ √ 2. On the other hand, θ av increases and reaches 50 at t = 1000/ √ 2. As shown in Fig. 4.1(d), the density ρ is nonuniform only at the very early stage and is almost uniform (i.e., almost ρ = 1) at t = 30/ √ 2. Corresponding to the nonuniformity of the density, the flow-velocity component u 1 perpendicular to the external force arises at the very early stage [Fig. 4.1(e)], but its magnitude is very small and practically vanishes at t = 10/ √ 2. Figures 4.1(f) and 4.1(g) show the long-time behavior of θ av (t) and the average speed, |u 2 | av (t) = 2 1/4 -1/4 |u 2 (x 1 , t)|dx 1 , up to t = 24000/ √ 2. The left figures show the double-logarithmic plot of θ av versus t and that of |u 2 | av versus t. In the right figures, the gradients of the curves in the left figures, i.e., α = d ln θ av /d ln t and β = d ln |u 2 | av /d ln t are plotted. If α and β approach constant values, say α 0 and β 0 , respectively, then we have the long-time behavior as θ av ≈ C θ t α 0 and |u 2 | av ≈ C u t β 0 with positive constants C θ and C u . From Figs. 4.1(f) and 4.1(g), it is still not clear whether α and β converge to finite values or not. But, if it is the case, it is likely that α 0 ≈ 0.66 and β 0 ≈ -0.34 ≈ -(1α 0 ).Figures 4.2(a)-4.2(g) show the behavior, corresponding to Figs. 4.1(a)-4.1(g), for a weaker external force (f 0 = 0.2). The tendency of the time evolution of the solution is similar to Fig. 4.1. However, since the magnitude of the force is 1/10, the resulting flow and the temperature rise are smaller. As Fig. 4.2(a) shows, the flow speed |u 2 |, which is smaller by one order of magnitude, takes the maximum at around t = 4/ √ 2 and decreases more slowly than in Fig. 4.1(a).Figure 4.2(b) shows that the increase of |θ av | is much slower compared to Fig. 4.1(b). One sees from Figs. 4.2(c)-4.2(e) that the nonuniformity of θ and ρ and the magnitude of u 1 are smaller by two orders of magnitude. In Figs. 4.2(f) and 4.2(g), we show the long-time behavior of θ av and |u 2 | up to an extremely large time, t = 336000/ √ 2. As

Figure 4 .

 4 Figures 4.2(a)-4.2(g) show the behavior, corresponding to Figs. 4.1(a)-4.1(g), for a weaker external force (f 0 = 0.2). The tendency of the time evolution of the solution is similar to Fig. 4.1. However, since the magnitude of the force is 1/10, the resulting flow and the temperature rise are smaller. As Fig. 4.2(a) shows, the flow speed |u 2 |, which is smaller by one order of magnitude, takes the maximum at around t = 4/ √ 2 and decreases more slowly than in Fig. 4.1(a).Figure 4.2(b) shows that the increase of |θ av | is much slower compared to Fig. 4.1(b). One sees from Figs. 4.2(c)-4.2(e) that the nonuniformity of θ and ρ and the magnitude of u 1 are smaller by two orders of magnitude. In Figs. 4.2(f) and 4.2(g), we show the long-time behavior of θ av and |u 2 | up to an extremely large time, t = 336000/ √ 2. As

  4.3 and 4.4 the results for the same cases as Figs. 4.1 and 4.2, respectively. That is, Fig. 4.3 is for Kn = 0.1, f 0 = 2, and Fig. 4.4 for Kn = 0.1, f 0 = 0.2. Figures 4.3(a) and 4.4(a) correspond respectively to Figs. 4.1(a) and 4.2(a); Figs. 4.3(b) and 4.4(b) to Figs. 4.1(b) and 4.2(b); Fig. 4.3(c) to Fig. 4.1(c) (it is impossible to obtain the corresponding figure for f 0 = 0.2); Figs. 4.3(d) and 4.4(c) to Figs. 4.1(f) and 4.2(f); and Figs. 4.3(e) and 4.4(d) to Figs. 4.1(g) and 4.2(g).

Figure 4 .

 4 [START_REF] Bardos | Sur les limites asymptotiques de la théorie cinétique conduisant à la dynamique des fluides incompressibles[END_REF] shows the result of the ensemble average over 96 independent runs. In addition, in Fig.4.3(e), β, which is the average of β over the time interval [t -500/ √ 2, t], is shown instead of β itself, since it is impossible to obtain the plot of β in a reasonable form. In Fig.4.

Figure 4 . 2 :

 42 Figure 4.2: BGK model with Kn = 0.1 and f 0 = 0.2

Figure 4 . 2 :

 42 Figure 4.2: Continued

Figure 4 . 3 : 2 ,Figure 4 . 3 :

 43243 Figure 4.3: Boltzmann equation with Kn = 0.1 and f 0 = 2. The ensemble average over 96 independent runs is shown. The β(t) is the time average of β over the interval [t -500/ √ 2, t].

Figure 4 . 4 :

 44 Figure 4.4: Boltzmann equation with Kn = 0.1 and f 0 = The result shown here is based on the ensemble average over N independent runs and time average over the interval [t -C/ √ 2, t + C/ √ 2], where N and C are chosen appropriately depending on t. The details are given in the main text. The β is the average of β over the time interval [t -10000/ √ 2, t].

4 with 4 .

 44 cosity and thermal conductivity and are functions of θ. Corresponding to the initial-boundary-value problem of the BGK model solved numerically, the above equations should be considered in the interval -1/4 < x 1 < 1/We should keep in mind that the compressible Navier-Stokes equations listed above hold only approximately because u 1 is not exactly zero in the numerical solution based on the BGK model. Now, we assume that ρ ≈ const and ∂u 2 /∂t is negligibly small in the third equation. Integrating this equation with respect to x 1 and taking into account the boundary condition, we have ∂u 2 /∂x 1 ≈ (g 0 ρ/2πµ) cos 2πx 1 . We insert the expression of ∂u 2 /∂x 1 in the forth equation and integrate it with respect to x 1 from x 1 = -1/4 to 1/4 assuming that ρ ≈ const and µ(θ) ≈ µ(θ av ). Then, we obtain µ(θ av )

  and consider a family of solutions of the scaled Boltzmann equation (2.2) supplemented with the diffuse reflection condition (2.35) and with the total mass condition (2.36). Assume that

F = M (Z e Φ + g ) as in (2.3) and satisfies the same assumptions as in Theorem 1. Assume moreover that the family of traces of g on the boundary

  Let Ma = Kn = while Fr = 1, and consider a family of solutions of the scaled Boltzmann equation (2.20) supplemented with the diffuse reflection condition (2.35) and with the total mass condition (2.36). Assume that F satisfies the same assumptions as in Theorem 2, and that the odd and even part of the relative fluctuation of distribution function, resp.

	1 2 (|v| 2 -5) ,
	this implies that
	u ∂Ω = 0 , and θ ∂Ω = 0 .
	Theorem 4.

g and h defined in (2.17) are continuous in v and satisfy the condition

g ∂Ω×R 3 → g ∂Ω×R 3 and h ∂Ω×R 3 → h ∂Ω×R 3

locally uniformly in x, v, where we recall that g and h are the weak limits of g and h in L 2 (Ω × R 3 ; (1 + |v| 2 )M dvdx) as → 0. Then g and h are given by the expressions (2.22) and

(2.27)

, where (u, θ) is a solution of the system (2.24) with the Dirichlet boundary condition

The left figures show the double-logarithmic plot of θ av versus t and that of |u 2 | av versus t. In the right figures, the gradients of the curves in the left figures, i.e., α = d ln θ av /d ln t and β = d ln |u 2 | av /d ln t are plotted. If α and β approach constant values, say α 0 and β 0 , respectively, then we have the long-time behavior as θ av ≈ C θ t α 0 and |u 2 | av ≈ C u t β 0 with positive constants C θ and C u . From Figs. 4.1(f) and 4.1(g), it is still not clear whether α and β converge to finite values or not. But, if it is the case, it is likely that α 0 ≈ 0.66 and β 0 ≈ -0.34 ≈ -(1α 0 ).

Acknowledgement. This work was started during a visit of the first author at Ecole polytechnique, and completed while the second author was visiting Kyoto University. We are grateful to both institutions for their hospitality and support.

satisfies the bound

Thus, if

then G α maps the closed ball B(0, R) ⊂ H 1 into itself and is a strict contraction on B(0, R). Hence G α has a unique fixed point in B(0, R)

In particular, for α = 1, this unique fixed point of G α is the unique solution of the Navier-Stokes equation in B(0, R) ⊂ H 1 .

The estimate

and Rellich's theorem imply that the map T is compact in

) by the Rellich compactness theorem, and the inequality above with v = u n shows that T u n → T u strongly in H 1 . On the other hand, for each u ∈ H 1 , the equation

Multiplying both sides of the motion equation by u and integrating over T 3 , one finds that

Appendix: Properties of A and B Lemma 3. The tensor field A and the vector field B defined by the formulas

satisfy the following properties.

(1) The orthogonality relations

(2) There exists a unique tensor field  and a unique vector field B in

(3) The tensor field  and the vector field B are of the form

Statement ( 1) is Lemma 5.3 in [START_REF] Golse | Fluid Dynamic Limits of the Kinetic Theory of Gases[END_REF], while statement (3) is Lemma 5.4 in [START_REF] Golse | Fluid Dynamic Limits of the Kinetic Theory of Gases[END_REF]. Statement (3) has been used systematically in the literature on the Boltzmann equation, referring to §7.31 in [START_REF] Chapman | The Mathematical Theory of Non-Uniform Gases[END_REF]. However, the discussion in [START_REF] Chapman | The Mathematical Theory of Non-Uniform Gases[END_REF] is incomplete. The key argument leading to the structure of  and B in statement ( 3) is the invariance of the linearized collision operator under the orthogonal group O 3 (R). It seems that the first complete proof of statement (3) following this idea is in [START_REF] Desvillettes | A remark concerning the Chapman-Enskog asymptotics[END_REF]. Statement (2) follows from Hilbert's lemma (the fact that the linearized collision operator L satisfies the Fredholm alternative).

Lemma 4. The tensor fields A and Â, and the vector fields B and B satisfy the following properties.

(1) For each i, j = 1, 2, 3, one has

(2) For each i, j, k, l = 1, 2, 3, one has

and

(3) For each i, j, k, l = 1, 2, 3, one has

and

Statement ( 2) is Lemma 3.4 in [START_REF] Golse | Fluid Dynamic Limits of the Kinetic Theory of Gases[END_REF] (proved in Appendix 2 of [START_REF] Golse | Fluid Dynamic Limits of the Kinetic Theory of Gases[END_REF]). See formulas (4.10) and (4.13b) in [START_REF] Bardos | Fluid dynamic limits of kinetic equations II: Convergence proofs for the Boltzmann equation[END_REF] for statement [START_REF] Arkeryd | The stationary Boltzmann equation in R n with given indata[END_REF], which is proved by a similar (though simpler) argument as statement [START_REF] Arkeryd | On a Taylor-Couette type bifurcation for the stationary nonlinear Boltzmann equation[END_REF]. (Notice the slight difference in normalization in the definitions of κ in [START_REF] Bardos | Fluid dynamic limits of kinetic equations II: Convergence proofs for the Boltzmann equation[END_REF]and here).

Proof of statement (3). By statement (3) in the previous lemma Bi

where the second equality follows form the fact that B ⊥ Ker L. Then, according to Lemma 4.3 in [START_REF] Golse | Fluid Dynamic Limits of the Kinetic Theory of Gases[END_REF], one has

for some λ ∈ R, while, by statement (1) of the previous lemma,

and since

we conclude that λ = 2 5 κ, which immediatly implies the first equality in statement [START_REF] Bardos | Sur les limites asymptotiques de la théorie cinétique conduisant à la dynamique des fluides incompressibles[END_REF]. The second equality is obtained in exactly the same manner. Contracting i, k and j, l, one finds the announced formula for c. Lemma 5. Let C be the tensor field defined by the formula

Then, for each i, j, k, l = 1, 2, 3, one has

for some µ ∈ R, by Lemma 4.3 in [START_REF] Golse | Fluid Dynamic Limits of the Kinetic Theory of Gases[END_REF]. By contraction on the indices k, l

where the first equality follows from the orthogonality relation B ⊥ Ker L by statement (1) in Lemma 3. Hence µ = 1 2 , and this immediatly implies the desired identity.


Taghreed Mohammadridha
email: taghreed.mohammadridha@irccyn.ec-nantes.fr

Claude H Moog
email: claude.moog@irccyn.ec-nantes.fr

Emmanuel Delaleau
email: emmanuel.delaleau@isen-bretagne.fr

Michel Fliess
email: michel.fliess@polytechnique.edu

Cédric Join
email: cedric.join@univ-lorraine.fr

A Variable Reference Trajectory for Model-Free Glycemia Regulation

Keywords: 1 diabetes, artificial pancreas, control algorithms, model-free control, intelligent proportional controller, variable reference trajectory, algebraic estimation techniques

The control design of an artificial pancreas, a hot research topic in diabetology, is tackled via the newly introduced model-free control and its corresponding "intelligent" proportional controller, which were already quite successful in many concrete and diverse situations. It results in an insulin injection for type 1 diabetes which displays via constant references a good nocturnal/fasting response, but unfortunately a poor postprandial behavior due to long hyperglycemia. When a variable reference is introduced, that switches between a constant one, when glycemia is more or less normal or moderate, and an exponential decay reference path, when a high glycemia rate indicates a meal intake, the results in silico, which employ real clinical data, become excellent. We obtain a bolus-shaped insulin injection rate during postprandial phases. The hyperglycemic peaks are therefore lowered a lot.

INTRODUCTION

Type-1 diabetes mellitus (T1DM), or insulin-dependent diabetes, is a chronic disease that results from the autoimmune destruction of the insulin-producing beta cells in the pancreas (see, e.g., [START_REF] Atkinson | Type 1 diabetes[END_REF], and the references therein). A functional insulin therapy (FIT) is a most effective pedagogical treatment (see, e.g., [START_REF] Howorka | Funktionelle Insulintherapie[END_REF], and the references therein) where exogenous insulin is injected, or infused subcutaneously via a pump, according to the patient's everyday life.

Rapid technological advancements justify the research of an artificial pancreas, i.e., of an automated device which consists of:

• a blood glucose sensor, or a continuous glucose monitor (CGM), {michel.fliess, cedric.join}@alien-sas.com

Projet Non-A, INRIA Lille -Nord-Europe, France.

• an insulin pump,

• a control algorithm for injecting the appropriate insulin dose.

Most diverse control techniques have been proposed and often tested in silico, i.e., via computer simulators. See, e.g., [START_REF] Bequette | Challenges and recent progress in the development of a closed-loop artificial pancreas[END_REF][START_REF] Chee | Closed-Loop Control of Blood Glucose[END_REF][START_REF] Cobelli | Artificial pancreas: Past, present, future[END_REF][START_REF] Doyle | Closed-loop artificial pancreas systems: Engineering the algorithms[END_REF][START_REF] Lunze | Blood glucose control algorithms for type 1 diabetic patients: A methodological review[END_REF][START_REF] Youssef | A review of closedloop algorithms for glycemic control in the treatment of type 1 diabetes[END_REF] for reviews. Obtaining a suitable mathematical modeling is therefore an essential issue (see, e.g., the previous reviews, and [START_REF] Join | Développement d'un nouveau modèle dédié à la commande du métabolisme glucidique appliqué aux patients diabétiques de type 1[END_REF][START_REF] Cobelli | Advancing our understanding of the glucose system via modeling: A perspective[END_REF][START_REF] Magdelaine | A long-term model of the glucose-insulin dynamics of type 1 diabetes[END_REF]). Like in many concrete applications, PID control (see, e.g., [START_REF] Cobelli | Artificial pancreas: Past, present, future[END_REF][START_REF] Marchetti | An improved PID switching control strategy for type 1 diabetes[END_REF][START_REF] Mohammad Ridha | Back stepping-based-PID controller designed for an artificial pancreas model[END_REF][START_REF] Weinzimer | Fully automated closed-loop insulin delivery versus semiautomated hybrid control in pediatric patients with type 1 diabetes using an artificial pancreas[END_REF]) and model-based predictive control (see, e.g., [START_REF] Del Favero | First use of model predictive control in outpatient wearable artificial pancreas[END_REF][START_REF] Hovorka | Nonlinear model predictive control of glucose concentration in subjects with type 1 diabetes[END_REF][START_REF] Lee | A closed-loop artificial pancreas based on MPC: human-friendly identification and automatic meal disturbance rejection[END_REF][START_REF] Magni | Model-predictive control of type 1 Diabetes: An in silico trial[END_REF][START_REF] Plank | Multicentric, randomized, controlled trial to evaluate blood glucose control by the model predictive control algorithm versus routine glucose management protocols in intensive care unit patients[END_REF]) are quite often encountered (see, e.g., [START_REF] Hovorka | Continuous glucose monitoring and closed-loop systems[END_REF]). 1 They are however far from being completely satisfactory:

1. Although PIDs do not require a precise mathematical model, their tuning (see, e.g., [START_REF] Åstrom | Advanced PID Control. Instrument Soc. Amer[END_REF][START_REF] O'dwyer | Handbook of PI and PID Controller Tuning Rules[END_REF]) may be extremely difficult and quite sensitive to "disturbances" like meals or stress.

2.

The reference trajectory which is connected to modelbased predictive control (see, e.g., [START_REF] Camacho | Model Predictive control[END_REF][START_REF] Grüne | Nonlinear Model Predictive Control -Theory and Algorithms[END_REF][START_REF] Richalet | Predictive Functional Control: Principles and Industrial Applications[END_REF]) is a most attractive and useful feature. It is however paid by the need of a mathematical modeling and by cumbersome optimization procedures.

Let us also mention [START_REF] Russell | Blood glucose control in type 1 diabetes with a bihormonal bionic endocrine pancreas[END_REF] which presents a wearable, automated, bihormonal, bionic pancreas that improved mean glycemic levels, with less frequent hypoglycemic episodes.

Our purpose is to show that the new model-free control (MFC) [START_REF] Fliess | Model-free control[END_REF], which has already been successfully applied in many concrete and diverse situations (see [START_REF] Fliess | Model-free control[END_REF][START_REF] Lafont | A modelfree control strategy for an experimental greenhouse with an application to fault accommodation[END_REF][START_REF] Menhour | A new model-free design for vehicle control and its validation through an advanced simulation platform[END_REF], and the references therein), might get closer to a real artificial pancreas. As a matter of fact the benefits of this control algorithm are numerous:

• As for PIDs there is no need of any mathematical modeling, which is anyway extremely difficult to write down in the case of diabetes, and, more generally, for biological systems.

• Contrarily to PIDs the tuning of the corresponding intelligent controller is straightforward.

• This controller, which may be easily implemented on a cheap programmable device [START_REF] Join | Intelligent" controllers on cheap and small programmable devices[END_REF], is moreover quite robust with respect to disturbances and system modifications.

• Our feedforward standpoint, i.e., the choice of a reference trajectory, is similar to what is done in flatnessbased control (see [START_REF] Fliess | Flatness and defect of non-linear systems: introductory theory and examples[END_REF][START_REF] Fliess | Continuous-time linear predictive control and flatness: a module-theoretic setting with examples[END_REF], and [START_REF] Lévine | Analysis and Control of Nonlinear Systems -A Flatness-Based Approach[END_REF][START_REF] Sira-Ramírez | Differentially Flat Systems[END_REF]). Contrarily to predictive control, optimal control becomes pointless.

An excellent description of the tremendous advantages of a model-free setting is provided by the following citation, which is extracted from Section 1.4.3 in [START_REF] Bastani | Model-free intelligent diabetes management using machine learning[END_REF], where a machine learning viewpoint is developed: Model-free methods, in either open-loop or closed-loop systems, have some advantages over their model-based counterparts. . . . One problem with the model-based methods is that they need an accurate model, and building such a model usually requires special and costly physiological studies of the processes involved in glucose regulation. However, model-free control algorithms rely solely on the "simple" type of data that patients normally collect. Another problem with the model-based methods is their dependence on a patient specific model. Thus any modification to a component of the mathematical model will require the parameters of the model-based method to be re-calculated. . . . We believe model-free methods are an important avenue of research that has not been exploited to its full advantage by the diabetes community.

Our contributions, which seem to be quite new in the huge, and, sometimes, old literature on automated pancreas, might be summarized as follows:

1. The design of a fully automated intelligent controller for blood glucose (BG) regulation including postprandial phase.

2. The design of an ad hoc time-varying glycemic reference trajectory, based on the results of [START_REF] Marchetti | An improved PID switching control strategy for type 1 diabetes[END_REF]. It is modified to reproduce an artificial insulin bolus, i.e., regular to rapid-acting insulin dose, when meals are taken.

3. This variable reference yields a postprandial bolusshaped injection rate which is more effective in handling hyper peaks of glycemia.

4. Clinical data are used to test the design.

Our paper is organized as follows. Model-free control and the corresponding controllers are briefly reviewed in Section 2. Section 3 and 4 present respectively the virtual patient simulator, and the system constraints and limitations. Implementation issues are discussed in Section 5. Promising simulation results are displayed in Section 6. Concluding remarks are given in Section 7. • the derivation order of y is 1 like in most concrete situations,

• α ∈ R is chosen by the practitioner such that αu and ẏ are of the same magnitude.

The following explanations on F might be useful:

• F is estimated via the measure of u and y,

• F subsumes not only the unknown system structure but also any perturbation.

Intelligent controllers

The loop is closed by intelligent proportional controller, or iP,

(2.2) u = - F -ẏ * + K P e α where
• y is the reference trajectory,

• e = yy is the tracking error,

• K P is the usual tuning gain.

Combining Equations (2.1) and (2.2) yields:

ė + K P e = 0
where F does not appear anymore. The tuning of K P becomes therefore quite straightforward. This is a major benefit when compared to the tuning of "classic" PIDs (see, e.g., [START_REF] Åstrom | Advanced PID Control. Instrument Soc. Amer[END_REF][START_REF] O'dwyer | Handbook of PI and PID Controller Tuning Rules[END_REF], and the references therein), which

• necessitate a "fine" tuning in order to deal with the poorly known parts of the pant,

• exhibit a poor robustness with respect to "strong" perturbations and/or system alterations.

See Figure 1 for the corresponding block diagram.

REMARK 2.1. See [START_REF] Fliess | Model-free control[END_REF] for academic comparisons with classic PIDs, and [START_REF] Gédouin | Experimental comparison of classical pid and model-free control: position control of a shape memory alloy active spring[END_REF] for a concrete one.

REMARK 2.2. See [START_REF] Fliess | Model-free control[END_REF] for the connection with the classic integral-proportional (PI) controllers.

REMARK 2.3. See [START_REF] Gao | On the centrality of disturbance rejection in automatic control[END_REF][START_REF] De Larminat | Automatique appliquée[END_REF] for a slightly different presentation of model-free control. See [START_REF] Fliess | Model-free control[END_REF] for other references.

Estimation of F

The calculations below stem from new estimation techniques (see [START_REF] Fliess | An algebraic framework for linear identification[END_REF][START_REF] Fliess | Closed-loop parametric identification for continuous-time linear systems via new algebraic techniques[END_REF], and [START_REF] Sira-Ramírez | Algebraic Identification and Estimation Methods in Feedback Control Systems[END_REF]).

First approach

The term F in Equation (2.1) may be assumed to be "well" approximated by a piecewise constant function F est . Rewrite then Equation (2.1) in the operational domain (see, e.g., [START_REF] Yosida | Operational Calculus[END_REF]): sY = Φ s +αU +y(0), where Φ is a constant. We get rid of the initial condition y(0) by multiplying both sides on the left by d ds : Y + s dY ds = -Φ s 2 + α dU ds . Noise attenuation is achieved by multiplying both sides on the left by s -2 . It yields in the time domain the realtime estimate, thanks to the equivalence between d ds and the multiplication by -t,

(2.3) F est (t) = - 6 τ 3 t t-τ [(τ -2σ )y(σ ) + ασ (τ -σ )u(σ )] dσ

Second approach

Close the loop with the iP (2.2):

(2.4)

F est (t) = 1 τ t t-τ ( ẏ -αu -K P e) dσ
2.4. Note the following facts:

• integrals (2.3) and (2.4) are low pass filters,

• τ > 0 might be quite small, i.e., online estimation is achieved,

• the integrals may of course be replaced in practice by classic digital filters.

Virtual patient

Our control algorithm is tested and evaluated in silico via a long-term linear time-invariant model, which was recently developed [START_REF] Magdelaine | A long-term model of the glucose-insulin dynamics of type 1 diabetes[END_REF]:

(3.5)      Ġ İ Ï Ḋ D     =       0 -k si 0 1 0 0 0 1 0 0 0-1 T 2 u -2 Tu 0 0 0 0 0 0 1 0 0 0 -1 T 2 r -2 Tr            G I İ D Ḋ     +       0 0 Ku V i T 2 u 0 0       u+       0 0 0 0 Kr V B T 2 r       r+      k l -128 M 0 0 0 0     
• G(t) is the glycemia (mg/dL),

• I(t) is the insulinemia (U/dL),

• D(t) is the increase of glycemia due to digestion of carbohydrate (mg/dL).min -1 ,

• u(t) is the insulin infusion rate (U/min),

• r(t) is the carbohydrate in the meal (mg).

The parameters are given in Table 1. Recall some advantage of this modeling with respect to existing ones. It was shown, for instance in [START_REF] Magdelaine | A long-term model of the glucose-insulin dynamics of type 1 diabetes[END_REF], that the parameters in [START_REF] Cobelli | Diabetes: Models, signals, and control[END_REF] cannot be estimated from standard clinical data (CGM, injection and meal intake). Most existing models introduce, however, some apparent equilibria which are not consistent with real life. Equation (3.5), though it is linear, or, more exactly, affine, appears to be a significant scientific alternative as it displays a long-term fit with clinical data unlike other short-term predictions models such as [START_REF] Cescon | Continuous-time interval model identification of blood glucose dynamics for type 1 diabetes[END_REF][START_REF] Kirchsteiger | Continuous-time interval model identification of blood glucose dynamics for type 1 diabetes[END_REF]. Fundamental concrete quantities can be deduced from Equation (3.5), that are not found in [START_REF] Cescon | Continuous-time interval model identification of blood glucose dynamics for type 1 diabetes[END_REF][START_REF] Cobelli | Diabetes: Models, signals, and control[END_REF][START_REF] Kirchsteiger | Continuous-time interval model identification of blood glucose dynamics for type 1 diabetes[END_REF]. Let us mention here:

• The basal rate: it is the constant insulin infusion rate that maintains glycemia at a constant value during fasting.

• The insulin sensitivity factor (ISF): it corresponds to the glycemic drop per unit of extra insulin when the basal infusion rate is correctly set.

• Raise: it is the glycemic increase due to the digestion of 15g of carbohydrate (CHO) while a basal insulin rate is infused. Note that insulinemia is then at equilibrium.

• U/P: additional insulin units, above the infused basal rate, per portion of CHO required to steer glycemia back to its initial level prior to digestion.

The parameters of the simulator are identified on the model thanks to the data provided by the Centres hospitaliers universitaires (CHU) de Nantes et Rennes3 (see Table 2). After the identification process there is still a modeling error (including measurement fluctuations) between the model output and the corresponding CGM data. This error signal is computed for each subject and added to the model. In this manner, the clinical CGM data are reconstructed for each case to simulate the corresponding virtual patient. Furthermore, it allows a better comparison with the openloop control as the CGM record of FIT controlled glycemia can now be reconstructed by the simulator.

For a non-diabetic subject during fasting/nocturnal period, where r = 0, G is maintained steady with respect to liver production and brain consumption of glucose (k l -128 M ) in normal range by a basal (equilibrium) level of u = u eq and hence basal insulinemia I b . In the T1DM case the main challenge is that with a very low constant rate (u < u eq ) during fasting/nocturnal conditions G diverges. It is due to the continuous liver glucose production resulting in fasting hyperglycemia [START_REF] Sheehan | Fasting hyperglycemia: etiology, diagnosis, and treatment[END_REF]. It yields Ġ > 0.

Thus, for an insulin-dependent diabetic, the control objective is to maintain G within the interval [70, 120] mg/dL relatively fast. Two main constraints are emphasized:

• hypoglycemic limit of 70 mg/dL,

• the input is nonnegative (u ≥ 0).
Insulin infusion rate cannot be reversed once it is spread out in the blood stream. If BG drops into hypoglycemia, the only possible control is to shut off the insulin pump temporarily and to have some carbohydrates to raise BG back to euglycemia, i.e., to a normal level of sugar in the blood. Therefore, the system is externally positive [START_REF] Farina | Positive linear systems: theory and applications[END_REF].

Control implementation

The first step in the design is to define the reference input to the controller. The time-varying reference trajectory which was introduced in [START_REF] Marchetti | An improved PID switching control strategy for type 1 diabetes[END_REF] for the purpose of PID switching control, is modified in order to be utilized here:

G ref (t) =      (G(t s ) -G * ) exp -t-t s τ ref + G * if G(t) > G * G r if G(t) ≤ G *
• t s is the switching time,

• τ ref is a design parameter,

• G r is the constant set point, G * is the switching threshold.

The following properties hold:

1. The time-varying trajectory starts decaying directly at G(t s ).

2. The switching threshold G * is set to 140 mg/dL while G r = 120 mg/dL to avoid hypoglycemia.

3. t s is the switching/reset time: to re-start the trajectory if (tt s ) > 45 min and Ġ > 0 (BG is still increasing).

When G is above the postprandial hyperglycemia level (G ≥ 140mg/dL), the controller switches from a constant reference G r = 120 mg/dL to an exponential-decay trajectory initiated at G ref (t) = G(t s ). It settles eventually at G r . The resulting discontinuity produces an impulsive control due to ė that appears in the iP control law. This bolus-shaped insulin rate permits to have a fast response in order to regulate G towards the normal level and to avoid extended postprandial hyperglycemia.

The control will first be designed and tested on IF2 subject, the parameters are specified as follows:

• τ = nT s , where T s is the sampling time, and n ≥ 1 is an integer,

• T s = 1 min, n = 30, τ = 30 min for smooth estimation,

• τ ref = 5 min,

• since G responds inversely to the control variable, α is negative and: α -1 = -0.07.

• K P is tuned individually aiming to have a minimum number of hypoglycemic events.

Once they are designed, the parameters τ, τ ref , α are kept constant. Only K P is adjusted per patient for an optimal performance. The closed-loop design, with constant or variable reference, will be compared to the open-loop FIT therapy with the same meal protocol.

Simulation results

The proportional gain is set respectively to 0.023 and to 0.018 for the fixed set point and for the variable reference trajectory (see Section 5). Figure 2 displays the difference between the two strategies in the postprandial phase. Postprandial peaks are reduced by (25 ± 23) mg/dL. IF2 has a one-day fasting phase, i.e., a non-glucidic regime, during which our control algorithm has a much better performance than the open-loop FIT control with a fasting hyperglycemia time lapse greater than 11 hours (see Figure 2b). This design is tested on the two other virtual patients IF3 and BE as shown in Figures 3 and4. The simulator is considered as a black box where the model parameters and meals are all unknown. The proportional gain is scaled to K P = 0.0065 and 0.01 for IF3 and BE respectively for variable-reference control. These values represent the upper bound with minimum hypoglycemic episodes during the simulated period. Tables 3, 4 and The frequent initialization of the exponential decay reference path appeared as a series of consecutive impulsive insulin rates. These bolus-shaped rapid pulses are infused only when a postprandial phase is detected, i.e., BG ≥ 140 mg/dL, positive rate. It yields an Ersatz, i.e., a replacement, for the missing pre-meal boluses. The postprandial peaks under closed-loop insulin administration occur more often than under FIT control. It might be due to the FIT injections anticipatory character, i.e., the boluses are infused before meals as shown in Figures 2, 3 and4 (parts d ande). For the fully automatic control, meals are unknown disturbances. As a result, G grows until the retarded plasma insulin reaches its peak (see Figures 2c,3c,4c). The insulin absorption rate is slow. According to the American Diabetes Association (ADA), the minimum insulin onset 4 is 15 minutes. It takes about 60 minutes to reach its peak with a minimum overall duration of 2-4 hours.

Our control design is tested on three different virtual patients: IF2, IF3, and BE. They are stable diabetics according to the Mean Amplitude Glycemic Excursion (MAGE) index. Unstable or brittle diabetics 5 are not included in this study. Severe BG swings from hyper to hypoglycemia, as a matter of fact, are hard to control and usually need extensive care.

The following two achievements summarize the benefits of our approach:

1. When postprandial phase is detected, artificial rapid insulin bolus are produced.

2. The overall BG mean value is improved and postprandial peaks are lowered. 4 It is the time lapse before insulin reaches the bloodstream and begins lowering blood glucose. 5 T1DM is an intrinsically unstable condition. The terminology brittle diabetes corresponds to an instability, causing life disruption, or prolonged hospitalization [START_REF] Vantyghem | Management strategies for brittle diabetes[END_REF].

Conclusion

Let us emphasize once again the following features of our intelligent proportional controller:

• a most promising in silico behavior,

• the few parameters which need to be calibrated are easy to tune,

• the implementation on cheap programmable devices is straightforward.

Some fault accommodations may also be easily handled via this setting (see [START_REF] Fliess | Model-free control[END_REF], and [START_REF] Lafont | A modelfree control strategy for an experimental greenhouse with an application to fault accommodation[END_REF] for a concrete application). 6Some other future investigations should also be added:

• test the controller on UVa/Padova simulator [START_REF] Cobelli | Diabetes: Models, signals, and control[END_REF],

• apply it on a large cohort of diabetics,

• compare it with other control schemes,

• relate K P with a known characteristic of the patient like the body mass.

Model-free control yields a quite elementary control synthesis, whereas questions stemming from computer simulations remain most complex. This fact, which was already mentioned in [START_REF] Fliess | Model-free control[END_REF], where several other references are provided, might be an "epistemological break" (rupture épistémologique in French) in the sense of Bachelard [START_REF] Bachelard | La Formation de l'esprit scientifique[END_REF].

2 2 2. 1

 21 Model-free control and intelligent controllers The ultra-local model Replace the unknown global description by the ultra-local model: (2.1) ẏ = F + αu where • the control and output variables are respectively u and y,

Figure 1 :

 1 Figure 1: MFC block diagram .

  5 present simulation statistics for open and closed-loop control, including a conventional iP controller with a constant reference trajectory. Hyperglycemia is then reduced in time. Hypoglycemia remains small. The sudden one-minute drop in Figure 3b at t = 12.6 h of ∆BG = 29 mg/dL is too fast to be the response to a bolus or an overdose. It might be a measurement issue. An unavoidable hypoglycemia episode follows. It happens for all three controllers.

Figure 2 :Figure 3 :Figure 4 :

 234 Figure 2: IF2 response for FIT, where K p = 0.023 or 0.018. (a): Meal intake and the corresponding increase of BG. (b): BG behavior. (c): Insulinemia. (d): The integral of the control rate per hour. (e): Controlled insulin injection minute rate.

Other mathematical techniques have also been investigated, like delay differential equations and partial differential equations (see, e.g.,[START_REF] Makroglou | Mathematical models and software tools for the glucose-insulin regulatory system and diabetes: an overview[END_REF], and the references therein).

See[START_REF] Fliess | Model-free control[END_REF] for more details.

Nantes and Rennes University-Hospital Centers.

This important topic seems to have been neglected to some extent. See however[START_REF] Vega-Hernandez | Actuator fault tolerant control for an artificial pancreas, 6th Int. Conf. Elec. Engin. Comput. Sci. Automat. Control[END_REF].