Existence analysis of a degenerate diffusion system for heat-conducting fluids

The existence of global weak solutions to a parabolic energy-transport system in a bounded domain with no-flux boundary conditions is proved. The model can be derived in the diffusion limit from a kinetic equation with a linear collision operator involving a non-isothermal Maxwellian. The evolution of the local temperature is governed by a heat equation with a source term that depends on the energy of the distribution function. The limiting model consists of cross-diffusion equations with an entropy structure. The main difficulty is the nonstandard degeneracy, i.e., ellipticity is lost when the fluid density or temperature vanishes. The existence proof is based on a priori estimates coming from the entropy inequality and the $H^{-1}$ method and on techniques from mathematical fluid dynamics (renormalized formulation, div-curl lemma).

The first three authors acknowledge partial support from the FWF, the Austrian Science Fund (FWF), grants F65 and W1245. The second author has been additionally supported by the grants P30000 and P33010 of the FWF. The fourth author acknowledges support from the Alexander von Humboldt Foundation. a heat equation for the background temperature governed by a Fourier law. We refer to Section 2 for modeling details.
A major difficulty of system (1) is the derivation of suitable a priori estimates. This issue will be tackled by exploiting the entropy structure of the system. This means that equations (1) can be written in the cross-diffusion form (4) ∂ t u = div(M∇ q), .
Still, this approach is not sufficient. Indeed, because of the degeneracy at θ = 0, we cannot expect to achieve any control on the gradient of ρ, and moreover, the bounds from the entropy estimate are not sufficient to conclude. Our idea, detailed below, is to apply well-known tools from mathematical fluid dynamics like H −1 estimates and compensated compactness. The originality of this work consists in the combination of these tools and entropy methods, which allows us to treat non-standard degeneracies.
1.1. State of the art. Equations (1) belong to the class of energy-transport models which have been investigated particularly in semiconductor theory [13]. The first energy-transport model for semiconductors was presented by Stratton [16]. First existence results were concerned with models with very particular diffusion coefficients (being not of the form (1)) [1,2] or with uniformly positive definite diffusion matrices [8]. Existence results for physically more realistic diffusion coefficients were shown in [6], but only for situations close to equilibrium. A degenerate energy-transport system with a simplified temperature equation was analyzed in [14]. Energy-transport models do not only appear in semiconductor theory. For instance, they have been used to model self-gravitating particle clouds [3] and the dynamics in optical lattices [5].
In [18], the global existence of weak solutions to the model in a bounded domain Ω with no-flux boundary conditions was proved. At first glance, equations (1) look simpler than (6) because of the additional diffusion in the energy equation. However, the ideas in [18] cannot be easily applied to (1). Indeed, the key idea in [18] was to introduce the variables u = ρθ and v = ρθ 2 and to apply the Stampacchia trunction method to a time-discretized version of The functionals´Ω ρ 2 θ b dx turn out to be Lyapunov functionals along solutions to (7) for suitable values of b ∈ R, leading to uniform gradient estimates. However, the additional term in the energy equation of (1) complicates the derivation of a priori estimates. Thus, the proof in [18] seems to be rather specific to system (6) and is not generalizable. Our idea is to treat (1) by combining entropy methods and tools from mathematical fluid dynamics, which may be also applied to other cross-diffusion systems.
Clearly, the entropy estimates are not sufficient to pass to the de-regularization limit in the approximate scheme. Further bounds are derived from the H −1 (Ω) method, i.e., we use basically (−∆) −1 ρ and (−∆) −1 E, respectively, as test functions in the weak formulation of (1) (second key idea). This method gives estimates for Combining these bounds with those coming from the entropy inequality and the conservation laws leads to estimates for ∇(ρθ) = √ ρθ∇ √ ρθ, ∇θ = θ∇ log θ and consequently for E in W 1,1 (Ω). Moreover, ∂ t E is bounded in some dual Sobolev space. This allows us to apply the Aubin-Lions lemma to E. Unfortunately, we do not obtain gradient estimates for ρ.
To overcome this issue, we use tools from mathematical fluid dynamics (third key idea). Let (ρ δ , θ δ ) be approximate solutions to (1) (in a sense made precise in Section 3). First, we write the mass balance equation in the renormalized form in the sense of distributions for smooth functions f with bounded derivatives. Let g another smooth function with bounded derivatives and introduce the vectors We deduce from the properties of f and g and the a priori estimates that div (t,x) U δ and curl (t,x) V δ are uniformly bounded in L 1 (Ω × (0, T )) and hence relatively compact in W −1,r (Ω) for some r > 1. The div-curl lemma implies that U δ · V δ = U δ · V δ a.e., where the bar denotes the weak limit of the corresponding sequence. Thus, f (ρ δ )g(θ δ ) = f (ρ δ ) g(θ δ ) a.e. A truncation procedure yields that ρ δ θ δ = ρθ, where ρ and θ are the weak limits of (ρ δ ) and (θ δ ), respectively. As (E δ ) converges strongly, by the Aubin-Lions lemma, we are able to prove that θ δ → θ and eventually ρ δ → ρ a.e. These limits allow us to identify the weak limits and to pass to the limit δ → 0 in the approximate equations. The approximate scheme contains additional terms which need to be treated carefully such that our arguments are more technical than presented here. In fact, we need three approximation levels; see Section 3 for details.

Formal derivation from a kinetic model
We consider a gas which is rarefied enough such that collisions between gas particles can be neglected, but there are thermalizing collisions at a fixed rate with a nonmoving background. This is modeled by sampling post-collisional velocities from a Maxwellian distribution with zero mean velocity and with the background temperature, which is determined from the assumptions of energy conservation as well as heat transport in the background governed by the Fourier law. These assumptions lead to the equations which are written in dimensionless form with a diffusive macroscopic scaling with the scaled Knudsen number 0 < ε ≪ 1. The gas is described by the distribution function f ε (x, v, t) with the velocity v ∈ R 3 , and the temperature of the background is θ ε (x, t). The gradient and Laplace operators are meant with respect to the position variable x, and the Maxwellian is given by Finally, the position density of the gas is defined by The right-hand side of the heat equation (12) has been chosen such that the sum of the kinetic energy of the gas and the thermal energy of the background is conserved. In [11], the energy-transport system (1) has been derived formally from (11)-(12) in the macroscopic limit ε → 0. We repeat the argument here for completeness.
In the computations, the moments of the Maxwellian up to order 4 will be needed: where v i , v j denote the components of v (i, j = 1, 2, 3). From (11)- (12), the local conservation laws for mass and energy, can be derived by integration of (11) with respect to v and, respectively, by integration of (11) against |v| 2 /2 and adding to (12). In a formal convergence analysis, we assume f ε → f , ρ ε → ρ, and θ ε → θ as ε → 0 and deduce from (11) that f = ρM(θ). With (14), we obtain for the kinetic energy density The limit of the mass flux is obtained by multiplication of (11) by v/ε, integration with respect to v, and passing to the limit, using again (14): Analogously, we compute the flux of the kinetic energy, for i = 1, 2, 3. Using these results in the limits of the conservation laws leads to (1).

Proof of Theorem 1
We approximate equations (1) in the following way. The time derivative is replaced by the implicit Euler discretization with parameter τ > 0. This is needed to avoid issues related to the time regularity. A higher-order H 4 regularization for φ = ∂h/∂ρ in the mass balance equation with parameter ε > 0 gives H 2 (Ω) regularity and compactness in W 1,4 (Ω). Furthermore, H 2 (Ω) and W 1,4 (Ω) regularizations for log θ with the same parameter are added to the energy balance equation. The W 1,4 (Ω) regularization is needed to derive estimates when using both log θ and −1/θ as test functions in (1). Furthermore, we add an additional H 2 (Ω) regularization for φ in the mass balance equation with parameter δ > 0, which removes the degeneracy of the diffusion matrix M in (4). Finally, we add the artificial heat flux ∆θ 3 in the energy density equation with the same parameter δ to obtain gradient estimates for the temperature, and we add the term θ −N log θ for some N > 0 to achieve an estimate for θ −(N +1) .
After having proved the existence of solutions to the approximate problem and some a priori estimates coming from the entropy inequality, we perform the limits ε → 0, τ → 0, and δ → 0 (in this order).
3.1. Solution of the approximate problem. We wish to solve a system which approximates (1) and is formulated in the variables φ and w = log θ, similarly as in (4). We interpret ρ and E = θ(1 + 3 2 ρ) as functions of (φ, w), i.e.
In this notation, the diffusion coefficients become Let T > 0 and let the approximation parameters τ > 0 (such that T /τ ∈ N), ε > 0, and δ > 0 be given. Furthermore, let 0 < N < 5 be a number needed for the approximation θ −N log θ in the energy balance equation.
Let σ ∈ (0, 1] and let (φ, w) be a fixed point of S(·, ·, σ). It is a solution to (16) (16) and (17), respectively, and add both equations. (We use 1 − e −w instead of −e −w as a test function in order to be able to treat the term ε´Ω(1 + e w )wψ 2 dx and to obtain the entropy and energy balance in one single equation.) Then To estimate the first integral I 1 , we use the entropy density (8), formulated in terms of the variables (ρ, E), The function h in the variables (ρ, E) is convex, since the determinant of its Hessian, for any (ρ 1 , E 2 ), (ρ 2 , E 2 ) > 0, and consequently, The second integral I 2 is nonnegative since The integrals I 3 , I 7 , and I 8 are estimated according to Therefore, we obtain from (19) where C > 0 is here and in the following a generic constant independent of τ , ε, and δ. This gives a uniform H 2 (Ω) estimate for φ and w, independent of σ (but depending on ε and δ), and hence the desired uniform estimate for (φ, w) in W 1,4 (Ω; R 2 ). By the Leray-Schauder fixed-point theorem, there exists a solution (φ k , w k ) : We reformulate equations (16)-(17) by inserting definition (15) of the diffusion coefficients and computing (we drop the superindex k) for test functions ψ 1 , ψ 2 ∈ H 2 (Ω).
Next, we consider the following term appearing in J 3 : where the last inequality follows from the fact that z → z log z is bounded from below. Furthermore, we deduce from Lemma 2, bound (33) for φ τ , and the Poincaré-Wirtinger inequality that This shows that The first integral on the right-hand side can be controlled by the last integral on the lefthand side of (38). The last integral on the right-hand side is controlled after applying Gronwall's inequality. The integrals J 4 , J 5 , and J 6 can be controlled by the expressions on the left-hand side of (38). We conclude that We deduce from this estimate and Young's inequality that This proves the lemma.
Step 3: Strong convergence of (ρ τ ) and (θ τ ). First, we prove a gradient bound for the particle density.
Lemma 5 (Bounds for the discrete time derivative). There exists a constant C > 0 which does not depend on τ such that Proof. We infer from (42) and (33) that Furthermore, Taking into account Lemma 3, the first three terms on the right-hand side are uniformly bounded. Since N < 5, the last term can be estimated from above by δ θ which is bounded because of (33). This finishes the proof.
Concerning (E τ ), Lemmas 4 and 5 allow us to apply the Aubin-Lions lemma in the version of [9] (or Theorem 3 in [7] with m = 1) to obtain a subsequence of (E τ ) (not relabeled) such that, as τ → 0, E τ → E strongly in L 1 (Ω T ).
The previous bounds and the strong convergences of (ρ τ ) and (θ τ ) allow us to pass to the limit τ → 0 in (28)-(29). For this, we observe that, by (42), Furthermore, by Lemma 3, The strong convergence of (θ τ ) to θ, the uniform bounds on (θ τ ), and the a.e. positivity of θ imply that Then (28)-(29) become in the limit τ → 0, for any test functions ψ 1 , ψ 2 ∈ C 2 0 (Ω T ). 3.5. Limit δ → 0. In this subsection, we need some tools from mathematical fluid dynamics, in particular the concept of renormalized solutions and the div-curl lemma. In the following, we denote by u δ the weak or distributional limit of a sequence (u δ ) whenever it exists. Let (ρ δ , E δ ) be a weak solution to (43)-(44) and set φ δ = log(ρ δ /θ Step 1: Renormalized mass balance equation. We compute the renormalized form of (43).
The div-curl lemma [12,Theorem 10.21] implies that U δ · V δ = U δ · V δ a.e. in Ω T , which means that Step 3: Proof of ρ δ θ δ = ρθ. We wish to relax the assumptions on the functions f and g. To this end, we introduce the truncation function T 1 ∈ C 2 ([0, ∞)) by T 1 (s) = s for 0 ≤ s < 1, T 1 (s) = 2 for s > 3, and T 1 is nondecreasing and concave in [0, ∞). Then we define T k (s) = kT 1 (s/k) for s > 0 and k ∈ N. It is possible to choose f = T k in (47). Together with Fatou's lemma and the boundedness of g, we infer that Furthermore, we deduce from (30) that such that we obtain for any k ≥ 2, Then the limit k → ∞ implies that We claim that both terms on the right-hand side converge to zero as k → ∞. Indeed, it follows from Fatou's lemma and the L 1 (Ω T ) bound for (ρ δ θ 2 δ ) from Lemma 3 that ρ δ (θ δ − T k (θ δ )) 1 + ρ while we deduce from Fatou's lemma and the L 2 (Ω T ) bound for (θ δ ), again from Lemma 3, that We infer from (49) that for any k ≥ 1, which implies, in the limit k → ∞, that (50) ρ δ θ δ = ρθ a.e. in Ω T .
The positivity of ρ δ implies that ρ ≥ 0 a.e. in Ω T . Note, however, that we cannot conclude that ρ > 0 a.e., since the control on φ δ is now lost.