\documentclass[reqno]{amsart} \usepackage{hyperref} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2016 (2016), No. 245, pp. 1--31.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu} \thanks{\copyright 2016 Texas State University.} \vspace{8mm}} \begin{document} \title[\hfilneg EJDE-2016/245\hfil Low Mach and P\'eclet number limit] {Low Mach and P\'eclet number limit for a model of stellar tachocline and upper radiative zones} \author[D. Donatelli, B. Ducomet, M. Kobera, S. Ne\v{c}asov\'a \hfil EJDE-2016/245\hfilneg] {Donatella Donatelli, Bernard Ducomet, \\ Marek Kobera, \v{S}\'arka Ne\v{c}asov\'a} \address{Donatella Donatelli \newline Dipartimento di Ingegneria e Scienze dell'Informazione e Matematica, Universita degli Studi dell'Aquila, 67100 L'Aquila, Italy} \email{donatella.donatelli@univaq.it} \address{Bernard Ducomet \newline CEA, DAM, DIF, F-91297 Arpajon, France} \email{bernard.ducomet@cea.fr} \address{Marek Kobera \newline Mathematical Institute of the Charles University at Prague, Sokolovsk\'a 83, 186 75 Praha 8, Czech Republic} \email{kobera@centrum.cz} \address{\v{S}\'arka Ne\v{c}asov\'a \newline Institute of Mathematics of the Academy of Sciences of the Czech Republic, \v{Z}itn\'a 25, 115 67 Praha 1, Czech Republic} \email{matus@math.cas.cz} \thanks{Submitted June 21, 2016. Published September 10, 2016.} \subjclass[2010]{35Q35, 76N15, 67W05} \keywords{Navier-Stokes-Fourier-Poisson system; radiation transfer; \hfill\break\indent compressible magnetohydrodynamics; rotation; stellar radiative zone; weak solution; \hfill\break\indent elliptic-parabolic initial boundary value problem; vanishing P\'eclet number; \hfill\break\indent vanishing Mach number; vanishing Alfv\'en number; classical physics; plasma} \begin{abstract} We study a hydrodynamical model describing the motion of internal stellar layers based on compressible Navier-Stokes-Fourier-Poisson system. We suppose that the medium is electrically charged, we include energy exchanges through radiative transfer and we assume that the system is rotating. We analyze the singular limit of this system when the Mach number, the Alfv\'en number, the P\'eclet number and the Froude number approache zero in a certain way and prove convergence to a 3D incompressible MHD system with a stationary linear transport equation for transport of radiation intensity. Finally, we show that the energy equation reduces to a steady equation for the temperature corrector. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \newtheorem{proposition}[theorem]{Proposition} \newtheorem{definition}[theorem]{Definition} \allowdisplaybreaks \section{Introduction}\label{i} Our motivation in this work is the rigorous analysis of the equations describing parts of stars called radiative zones which are one of the most basic structures constituting stars among cores, convection zones, photospheres and atmospheres. Our model can be also applied to tachoclines which are transition layers between convection and radiative zones of stars. In this context it is conjectured that magnetic field of stars arises when poloidal orientation of magnetic fields changes to toroidal and that a dynamo effect is present in tachoclines \cite{Sh}. Tachoclines are not homogeneous and stable structures and they move steadily. In their upper parts the P\'eclet numbers are high (of the order 600), but in the vicinity of the radiative part they drop below 1. Their distinctive feature concerns rotation, na\"ively speaking the convective zone behaves in this respect as a fluid and rotates differentially, whereas the radiative zone more like a solid and rotates as a rigid body. The origin of these rotational changes has preocuppied astrophysicists and astronomers, particularly in connection with helioseismological observations \cite{Bro}. Gravitational forces in these regions are high, however the fluid is no longer strongly stratified as show non-dimensional numbers associated to the \emph{solar tacho-cline}. Namely \emph{the Froude number Fr} measuring the strength of gravitational interactions (see Section \ref{fsa} below for precise definitions) is $Fr = 3.11 \times 10^{-3} U$, where $U$ is the referential speed of flow in SI units. \emph{The Mach number Ma\/} measuring the compressibility is $Ma = 1.49 \times 10^{-7} U$, i.e. the fluid is almost incompressible for sufficiently slow motions and one has $Fr^2 \sim Ma$ (Mach number is due to high temperatures when radiation dominates). Finally \emph{P\'eclet number Pe\/} need not to be sufficiently small in the solar tachocline (we assume $Ma^2 = Pe$), but thermal diffusivity in \textbf{giant stars} can be seven orders of magnitude larger than that of the Sun (see \cite[page 22]{Gar}). Notice in conclusion that our low stratification model can be applied to other compact stellar objects, as the fraction of $Fr$ and $Ma$ depends on the ratio of temperature, density and is inversely proportional to the square of characteristic length. Therefore white dwarfs are too cold to be described by low stratification models, but neutron stars, especially newly born are not. Validity of classical MHD may be restricted to their (outer) crusts though; in their superfluid cores a quantum description is inevitable. Let us complete this physical introduction by drawing the reader's attention to the fact that models in stellar physics are computationally time consuming. Rieutord \cite{Rie} has estimated for example that modelling a single supergranule on the Sun would require having more than power of the Sun at our disposal! That is why Ligni\`eres \cite{Lig} has initiated studies of models at small P\'eclet number as through the Boussinesq-Oberbeck approximation density variations with temperature enter through the buoyancy force only and moreover temperature can be expressed by the velocity field. In our previous work \cite{ddn} we analyzed a thick disk model for the Mach number of order $\varepsilon$, $\varepsilon\to 0$ whereas the Peclet number was of order 1. Instead as in \cite{NRT}, in the present one we consider a model where the Peclet number is of order $\epsilon^2$ and the domain is general. The mathematical model we consider is the compressible heat conducting MHD system \cite{DF} describing the motion of a viscous plasma confined in $ \Omega $, a 3D domain, moreover as we suppose a global rotation of the system, some new terms appear due to the change of frame and we also suppose that the fluid exchanges energy with radiation through radiative cooling/heating (see \cite{DF}, \cite{dkn}), but neglect radiative accelerations. More precisely, the non-dimensional system of equations giving the evolution of the mass density $\varrho = \varrho (t,x)$, the velocity field $\vec{u} = \vec{u}(t,x)$, the (divergence-free) magnetic field $\vec B=\vec B(x,t)$, and the radiative intensity $I=I(x,t,\vec \omega,\nu)$ as functions of the time $t\in (0,T)$, the spatial coordinate $x=(x_1,x_2,x_3) \in \Omega \subset \mathbb{R}^3$, and (for $I$) the angular and frequency variables $(\vec \omega,\nu)\in \mathcal{S}^2\times \mathbb{R}_+$, reads as follows \begin{gather} \label{i1} \partial_t \varrho + \operatorname{div}_x (\varrho \vec{u}) = 0\quad \text{in } (0,T) \times \Omega, \\ \label{i2} \begin{aligned} &\partial_t (\varrho \vec{u}) + \operatorname{div}_x (\varrho \vec{u} \otimes \vec{u}) + \nabla_x p +2\varrho \vec \chi \times \vec{u} \\ &= \operatorname{div}_x \mathbb{S}+\varrho \nabla \Psi +\frac{1}{2}\varrho \nabla_x |\vec \chi \times \vec x|^2+ \vec{j}\times \vec B\quad \text{in } (0,T) \times \Omega, \end{aligned}\\ \label{i3} \begin{aligned} &\partial_t ( \varrho e ) + \operatorname{div}_x ( \varrho e \vec{u} ) + \operatorname{div}_x \vec{q} \\ &= \mathbb{S}:\nabla_x \vec{u} -p\operatorname{div}_x \vec{u}+\vec j\cdot \vec E - S_E\quad \text{in } (0,T) \times \Omega, \end{aligned}\\ \label{i4} \frac{1}{c} \partial_t I + \vec \omega \cdot \nabla_x I = S\quad \text{in } (0,T) \times \Omega \times (0, \infty) \times \mathcal{S}^2, \\ \label{i5} \partial_t \vec{B} + \operatorname{curl}_x( \vec{B} \times \vec{u} ) + \operatorname{curl}_x( \lambda\ \operatorname{curl}_x \vec{B} )= 0\quad \text{in } (0,T) \times \Omega,\\ \label{i6bis} -\Delta \Psi = 4\pi G (\tilde{\eta}\varrho+g) \quad \text{in } (0,T) \times \Omega. \end{gather} In the electromagnetic source terms, electric current $\vec j$ and electric field $\vec E$ are interrelated by \emph{Ohm's law} \[ \vec j=\sigma(\vec E+\vec{u}\times \vec B), \] and \emph{Amp\`ere's law} \[ \zeta\vec j=\operatorname{curl}_x \vec B, \] where $\zeta>0$ is the (constant) magnetic permeability { and $\sigma$ is the coefficient of electric conductivity. This is a simplified version of Ohm's law for plasmas as both the Hall effect and the ambipolar diffusion from density gradients and the electron inertia are neglected}. Moreover in \eqref{i5} $\lambda=\lambda(\vartheta)>0$ is the magnetic diffusivity of the fluid. In \eqref{i6bis} $\Psi$ is the gravitational potential and the corresponding source term in \eqref{i2} is the Newton force $\varrho \nabla \Psi$. $G$ is the Newton constant and $g$ is a given function, modelling an external gravitational effect. Supposing that $\varrho$ is extended by 0 outside $\Omega$ and solving \eqref{i6bis}, we have \[ \Psi(t,x)= G\int_{\Omega}K(x-y)(\tilde{\eta}\varrho(t,y)+g(y))\,dy, \] where $K(x)=\frac{1}{|x|}$, and the parameter $\tilde{\eta}$ may take the values 0 or 1: for $\tilde{\eta}=1$ selfgravitation is present and for $\tilde{\eta}=0$ gravitation only acts as an external field. We also assume that the system is globally rotating at uniform velocity $\chi$ around the vertical direction $\vec e_3$ and we note $\vec \chi=\chi\vec e_3$. Then {Coriolis acceleration term 2$\varrho \vec \chi \times \vec{u}$ appears in the system, together with} the centrifugal force term $\varrho \nabla_x |\vec \chi \times \vec x|^2$ (see \cite{C}). We consider here the simplified model studied in \cite{DuNe} where radiation does not appear in the momentum equation (see also \cite{TSGKS}): only the source term $S_E$ is present, in the energy equation $$ S_E(t,x) = \int_{\mathcal{S}^2} \int_0^\infty S(t,x, \vec \omega, \nu) \,d \vec \omega \,d \nu. $$ The symbol $p=p(\varrho,\vartheta)$ denotes the thermodynamic pressure and $e=e(\varrho,\vartheta)$ is the specific internal energy, interrelated through {Maxwell's relation} \begin{equation} \label{i5bis} \frac{\partial e}{\partial \varrho} = \frac{1}{\varrho^2} \Big( p(\varrho,\vartheta) - \vartheta \frac{\partial p}{\partial \vartheta} \Big). \end{equation} Furthermore, $\mathbb{S}$ is the Newtonian viscous stress tensor determined by \begin{equation} \label{i6} \mathbb{S} = \mu \Big( \nabla_x \vec{u} + \nabla_x^t \vec{u} - \frac{2}{3} \operatorname{div}_x \vec{u}\, \mathbb{I}\Big) + \eta \operatorname{div}_x \vec{u} \mathbb{I}, \end{equation} where the shear viscosity coefficient $\mu = \mu(\vartheta) > 0$ and the bulk viscosity coefficient $\eta = \eta(\vartheta) \geq 0$ are effective functions of the temperature. Similarly, $\vec q$ is the heat flux given by Fourier's law \begin{equation} \label{i7} \vec q = - \kappa \nabla_x \vartheta, \end{equation} with the heat conductivity coefficient $\kappa = \kappa(\vartheta) > 0$. Finally, \begin{equation} \label{i8} S = S_{a,e} + S_s , \end{equation} where \begin{equation} \label{i9} S_{a,e} = \sigma_a \big( B(\nu, \vartheta) - I \big), \quad S_s = \sigma_s ( \tilde I - I ). \end{equation} In this formula $\tilde I:= \frac{1}{4 \pi} \int_{\mathcal{S}^2} I(\cdot , \vec \omega) \,d \vec \omega$ and $B(\nu, \vartheta)=2h\nu^3c^{-2}(e^{\frac{h\nu}{k\vartheta}}-1)^{-1}$ is the radiative equilibrium function where $h$ and $k$ are the Planck and Boltzmann constants, $\sigma_a = \sigma_a(\nu, \vartheta) \geq 0$ is the absorption coefficient and $\sigma_s = \sigma_s(\nu, \vartheta) \geq 0$ is the scattering coefficient. More restrictions on these structural properties of constitutive quantities will be imposed in Section \ref{m} below. System \eqref{i1}--\eqref{i6bis} is supplemented with the ``no-slip, thermal insulation, perfect conductor, no reflection, no radiative entropy flux'' boundary conditions: \begin{gather} \label{i11} \vec u|_{\partial \Omega} = 0, \quad \vec{q} \cdot \vec{n} |_{\partial \Omega} = 0,\quad \vec B\cdot \vec n|_{\partial \Omega} = 0,\quad \vec E\times \vec n|_{\partial \Omega} = \vec 0, \\ \label{i12} I (t,x,\nu, \vec \omega) = 0 \quad \text{on } \Gamma_{\_},\quad \vec q^R\cdot\vec n (x) = 0 \quad \text{for } x \in \partial \Omega, \end{gather} where $\vec n$ denotes the outer normal vector to $\partial \Omega$, $\Gamma_{-} := \{ (x, \vec\omega) \in \partial \Omega\times\mathcal{S}^2 : \vec \omega \cdot \vec{n}_x \leq 0\}$ and the radiative entropy flux $\vec q^R$ will be defined in the next Section. Similarly we define $\Gamma_{+} : = \partial \Omega\times\mathcal{S}^2 \setminus \Gamma_{-}$. Let us mention that previous works have been achieved in the previous framework but, to our knowledge, not in the case of rotating fluid with radiation (with the exception of \cite{ddn}). Among them: Kuku\v{c}ka \cite{K} studied the case when Mach and Alfv\'en number go to zero in the case of a bounded domain and Novotn\'y and collaborators \cite{NRT} investigated the problem in the case of strong stratification. Let us also mention the works of Trivisa et al. \cite{TK} and Wang et al.\cite {HW}, and related articles of Jiang et al. \cite{JJL,JJL1,JJLX}. Our work differs from theirs in that we take a larger Froude number and add radiation and non-inertial effects. This paper is organized as follows. In Section \ref{m}, we list the principal hypotheses imposed on constitutive relations, introduce the concept of weak solution to problem \eqref{i1}--\eqref{i12}, and state the existence result for our model. In Section \ref{fsa} we compute the formal asymptotics of the problem. Uniform bounds imposed on weak solutions by the data are derived in Section \ref{ge}. The convergence theorem is proved in Section \ref{pconv}. Existence of a solution for the target system is briefly given in the Appendix. \section{Hypotheses and stability result}\label{m} As in \cite {DFN} we consider a pressure law in the form \begin{equation} \label{m1} p(\varrho, \vartheta) = \vartheta^{5/2} P ( \frac{\varrho}{\vartheta^{3/2}} ) + \frac{a}{3} \vartheta^4, \quad a > 0, \end{equation} where $P: [0, \infty) \to [0, \infty)$ is a given function with the following properties: \begin{gather} \label{m2} P \in C^2([0, \infty)), \quad P(0) = 0, \quad P'(Z) > 0 \quad \text{for all } Z \geq 0, \\ \label{m3} 0 < \frac{ \frac{5}{3} P(Z) - P'(Z) Z }{Z} < c \quad \text{for all } Z \geq 0, \\ \label{m4} \lim_{Z \to \infty} \frac{P(Z)}{Z^{5/3}} = p_\infty > 0. \end{gather} According to {Maxwell's equation \eqref{i5bis}}, the specific internal energy $e$ is \begin{equation} \label{m5} e(\varrho, \vartheta) = \frac{3}{2} \vartheta \frac{\vartheta^{3/2}}{\varrho} P ( \frac{\varrho}{\vartheta^{3/2}} ) +{a} \frac{\vartheta^4}{\varrho}, \end{equation} % ( ) and the associated specific entropy reads \begin{equation} \label{m6} s(\varrho, \vartheta) = M ( \frac{\varrho}{\vartheta^{3/2}} ) + \frac{4a}{3} \frac{\vartheta^3}{\varrho}, \end{equation} with \[ M'(Z) = - \frac{3}{2} \frac{ \frac{5}{3} P(Z) - P'(Z) Z } {Z^2} < 0. \] To ensure positivity of the total entropy production rate, as in \cite{ddn}, in this paper we explicitly introduce the entropy for the photon gas in the sequel. The transport coefficients $\mu$, $\eta$, $\kappa$ and $\lambda$ are continuously differentiable and Lipschitz functions of the absolute temperature with the properties, \begin{gather} \label{m7} c_1 (1 + \vartheta) \leq \mu (\vartheta), \quad \mu'(\vartheta) < c_2, \quad 0 \leq \eta(\vartheta) \leq c_3(1 + \vartheta), \\ \label{m8} c_1 (1 + \vartheta^r) \leq \kappa(\vartheta)\leq c_2 (1 + \vartheta^r) \\ \label{m8z} c_5(1 + \vartheta)<\lambda(\vartheta)\le c_4(1 + \vartheta^p) \end{gather} for any $\vartheta \geq 0$, for a $1\le p<\frac{17}{6}$ and $r = 3$. Moreover we assume that $\sigma_a$, $\sigma_s$, $B$ are continuous functions of $\nu$, $\vartheta$ such that \begin{gather} \label{m9} 0 < \sigma_a (\nu, \vartheta)\le c_1\,\quad 0 \leq \sigma_s (\nu, \vartheta),\quad |\partial_{\vartheta} \sigma_a (\nu, \vartheta)|, \quad |\partial_{\vartheta} \sigma_s (\nu, \vartheta)| \leq c_1, \\ \label{m9bis} 0 \leq \sigma_a (\nu, \vartheta) B (\nu, \vartheta), \quad |\partial_{\vartheta} \{ \sigma_a (\nu, \vartheta) B (\nu, \vartheta)\}|\leq c_2, \\ \label{m10} \sigma_a (\nu, \vartheta), \sigma_s (\nu, \vartheta), \quad \sigma_a (\nu, \vartheta) B(\nu, \vartheta) \leq h(\nu) ,\quad h \in L^1(0,\infty). \end{gather} for all $\nu \geq 0$, $\vartheta \geq 0$, where $c_{1,2,3,4,5}$ are positive constants. Let us recall some definitions introduced in \cite{DFN}. $\bullet$ In the weak formulation of the Navier-Stokes-Fourier system the \emph{equation of continuity} \eqref{i1} is replaced by its (weak) renormalized version \cite{DL} represented by the family of integral identities \begin{equation} \label{m12} \begin{aligned} &\int_0^T\!\! \int_\Omega{ \Big[ \Big( \varrho + b(\varrho) \Big) \partial_t \varphi + \Big(\varrho + b(\varrho) \Big) \vec{u} \cdot \nabla_x \varphi + \Big( b(\varrho) - b'(\varrho) \varrho \Big) \operatorname{div}_x \vec{u}\ \varphi \Big]}\,dx\,dt \\ & = -\int_\Omega{\! \Big( \varrho_0 + b(\varrho_0) \Big) \varphi (0, \cdot)}\,dx \end{aligned} \end{equation} satisfied for any $\varphi \in C^\infty_c([0, T) \times \overline{\Omega})$, and any $b \in C^\infty([0,\infty))$, $b' \in C^\infty_c([0, \infty))$, where \eqref{m12} implicitly includes the initial condition $\varrho(0, \cdot) = \varrho_0$. $\bullet$ Similarly, the \emph{momentum equation} \eqref{i2} is replaced by \begin{equation} \label{m13} \begin{aligned} &\int_0^T \int_{\Omega} \Big( (\varrho \vec{u}) \cdot \partial_t \vec{\varphi} + (\varrho \vec{u} \otimes \vec{u}): \nabla_x \vec{\varphi} + p \operatorname{div}_x \vec{\varphi} + 2\varrho \vec \chi \times \vec{u} \cdot \vec{\varphi}\Big) \,dx \, dt\\ &=\int_0^T \int_{\Omega} \Big( \mathbb{S} : \nabla_x \vec{\varphi} - \varrho \nabla_x \Psi \cdot \vec{\varphi} -\vec{j} \times \vec{B} \cdot \vec{\varphi} -\frac12\varrho \nabla_x |\vec \chi \times \vec x|^2 \cdot \vec{\varphi} \Big) \,dx \,dt \\ &\quad - \int_{\Omega} (\varrho \vec{u})_0 \cdot \vec{\varphi} (0,\cdot )\,dx \end{aligned} \end{equation} for any $\vec{\varphi} \in C^\infty_c([0, T) \times \Omega; \mathbb{R}^3 )$. As usual, for \eqref{m13} to make sense, the field $\vec{u}$ must belong to a certain Sobolev space with respect to the spatial variable and we require that \begin{equation} \label{m14} \vec{u} \in L^2(0,T; W^{1,2}_0 (\Omega; \mathbb{R}^3)), \end{equation} where \eqref{m14} already includes the no-slip boundary condition \eqref{i11}. $\bullet$ The \emph{magnetic equation} \eqref{i5} is replaced by \begin{equation} \label{v9} \int_0^T \int_{\Omega} \Big( \vec{B} \cdot \partial_t \varphi - ( \vec{B} \times \vec{u} + \lambda \operatorname{curl}_x \vec{B}) \cdot \operatorname{curl}_x \varphi \Big) \,dx \,dt + \int_{\Omega} \vec{B}_0 \cdot \varphi (0, \cdot) \,dx=0, \end{equation} to be satisfied for any vector field $\varphi \in \mathcal{D}([0,T) \times \Omega;\mathbb{R}^3)$. Here, according to the boundary conditions, one has to take \begin{equation} \label{v10} \vec{B}_0 \in L^2(\Omega),\quad \operatorname{div}_x \vec{B}_0 = 0 \quad \text{in } \mathcal{D}'(\Omega),\quad \vec{B}_0 \cdot \vec{n}|_{\partial \Omega} = 0. \end{equation} Following \cite[Theorem 1.4]{TEM}, $\vec{B}_0$ belongs to the closure of all solenoidal functions from $\mathcal{D}(\Omega)$ with respect to the $L^2$-norm. Anticipating (see \eqref{m18} below) we see that \[ \vec{B} \in L^{\infty}(0,T; L^2(\Omega; \mathbb{R}^3)),\quad \operatorname{curl}_x \vec{B} \in L^2 (0,T; L^2(\Omega; \mathbb{R}^3)) \] and we deduce from \eqref{v9} that \[ \operatorname{div}_x \vec{B}(t) = 0 \text{ in } \mathcal{D}'(\Omega),\quad \vec{B}(t) \cdot \vec{n}|_{\partial \Omega} = 0 \text{ for a.a. } t \in (0,T). \] In particular, using \cite[Theorem 6.1]{DULI}, we conclude that \begin{equation} \label{v11} \vec{B} \in L^2(0,T; W^{1,2}(\Omega; \mathbb{R}^3)),\quad \operatorname{div}_x \vec{B} (t) = 0,\quad \vec{B} \cdot \vec{n}|_{\partial \Omega} = 0 \quad \text{for a.a. } t \in (0,T). \end{equation} $\bullet$ From \eqref{i2} and \eqref{i3} we find the \emph{energy conservation law} \begin{equation} \label{b15} \begin{aligned} &\partial_t \Big( { \frac12} \varrho |\vec{u}|^2 + \varrho e + {\frac1{2 \zeta}} |\vec{B}|^2 \Big) +\operatorname{div}_x \Big( ({ \frac12} \varrho |\vec{u}|^2 + \varrho e + p) \vec{u} + \vec{E} \times \vec{B} - \mathbb{S} \vec{u}+\vec{q} \Big) \\ & = \varrho \nabla_x \Psi \cdot \vec{u}+\frac12\varrho \nabla_x | \vec \chi\times\vec x|^2\cdot \vec{u}-S_E. \end{aligned} \end{equation} As the gravitational potential $\Psi$ is determined by equation \eqref{i6bis} considered on the whole space $\mathbb{R}^3$, the density $\varrho$ being extended to be zero outside $\Omega$ we observe that $$ \int_{\Omega} \varrho \nabla_x \Psi \cdot \vec{u} \,dx = \frac{{\rm d}}{{\rm d}t} {\frac12} \int_{\Omega} \varrho \Psi \,dx, $$ in the same stroke $$ \frac12\int_{\Omega} \varrho \nabla_x |\vec \chi\times\vec x|^2\cdot \vec{u} \,dx = \frac{{\rm d}}{{\rm d}t}\frac12 \int_{\Omega} \varrho |\vec \chi\times\vec x|^2 \,dx. $$ Denoting now by $E^R$ the radiative energy given by \begin{equation} \label{i0} E^R(t,x) = \frac{1}{c} \int_{\mathcal{S}^2} \int_0^\infty I(t,x, \vec \omega, \nu) \,d \vec \omega \,d \nu, \end{equation} and integrating the radiative transfer equation \eqref{i4}, we obtain \[ \partial_t \int_{\Omega} E^R \,dx + \int\int_{\partial \Omega \times \mathcal{S}^2, \ \vec \omega \cdot \vec{n} \geq 0 } \int_0^\infty I (t,x,\vec \omega, \nu)\, \vec \omega \cdot \vec{n} \, {\rm d}\nu \, {\rm d} \vec \omega \,dS_x \ =\int_{\Omega} S_E\,dx. \] so, by using boundary conditions, we can integrate \eqref{b15}, as follows, \begin{align*} &\frac{{\rm d}}{{\rm d}t} \int_{\Omega}\Big( { \frac12 } \varrho |\vec{u}|^2 + \varrho e + {\frac1{2 \zeta}} |\vec{B}|^2\Big) + \int_{\partial\Omega}( ({ \frac12} \varrho |\vec{u}|^2 + \varrho e + p) \vec{u} + \vec{E} \times \vec{B} - \mathbb{S} \vec{u}+\vec{q} )\cdot \vec{n}\,dS \\ &=\int_{\Omega}(\varrho \nabla_x \Psi \cdot \vec{u}+\frac12\varrho \nabla_x |\vec \chi\times\vec x|^2\cdot \vec{u}-S_E)\,\,dx\\ \end{align*} \begin{equation} \label{b17} \begin{aligned} &\frac{{\rm d}}{{\rm d}t} \int_{\Omega}\Big( { \frac12 } \varrho |\vec{u}|^2 + \varrho e + {\frac1{2 \zeta}} |\vec{B}|^2 - {\frac12 } \varrho \Psi - {\frac12}\varrho |\vec \chi\times\vec x|^2+E^R\Big) \,dx \\ &+\int\int_{\Gamma_{+}} \int_0^\infty I (t,x,\vec \omega, \nu) \, \vec \omega \cdot \vec{n} \,d\nu \,d \vec \omega \,dS_x = 0 \end{aligned} \end{equation} by \eqref{i11} and \eqref{i12}. $\bullet$ Finally, dividing \eqref{i3} by $\vartheta$ and using Maxwell's relation \eqref{i5bis}, we obtain the \emph{entropy equation} \begin{equation} \label{m16} \partial_t ( \varrho s ) + \operatorname{div}_x ( \varrho s \vec{u} )+\operatorname{div}_x (\frac{\vec q}{\vartheta})=\varsigma, \end{equation} where \begin{equation} \label{iprod} \varsigma=\frac{1}{\vartheta} \Big(\mathbb{S}:\nabla_x \vec{u} -\frac{\vec q\cdot \nabla_x\vartheta}{\vartheta}+\frac{\lambda}{\zeta}|\operatorname{curl}_x \vec B|^2\Big) -\frac{S_E}{\vartheta}, \end{equation} where the first term $\varsigma_m:=\frac{1}{\vartheta} (\mathbb{S}:\nabla_x \vec{u} -\frac{\vec q\cdot \nabla_x\vartheta}{\vartheta}+\frac{\lambda}{\zeta}|\operatorname{curl}_x \vec B|^2)$ is the (positive) electromagnetic matter entropy production. To identify the second term in \eqref{iprod}, let us recall from \cite{BAL} the formula for the entropy of a photon gas \begin{equation} \label{iprodrad} s^R=-\frac{2k}{c^3}\int_0^{\infty}\int_{\mathcal{S}^2} \nu^2 [ n\log n-(n+1)\log(n+1)]d\vec \omega d\nu, \end{equation} where $n=n(I)=\frac{c^2I}{2h\nu^3}$ is the occupation number. Defining the radiative entropy flux \begin{equation} \label{ifluxrad} \vec q^R=-\frac{2k}{c^2}\int_0^{\infty}\int_{\mathcal{S}^2} \nu^2 [ n\log n-(n+1)\log(n+1)]\vec \omega\,d\vec \omega d\nu, \end{equation} and using the radiative transfer equation, we obtain the equation \begin{equation} \label{ierad} \partial_t s^R+\operatorname{div}_x \vec q^R =-\frac{k}{h}\int_0^{\infty}\int_{\mathcal{S}^2} \frac{1}{\nu} \log \frac{n}{n+1} S\,d\vec \omega d\nu=:\varsigma^R. \end{equation} With the identity $\log \frac{n(B)}{n(B)+1}=-\frac{h\nu}{k \vartheta}$ with $B=B(\vartheta,\nu)$ denoting Planck's function, and using the definition of $S$, the right-hand side of \eqref{ierad} rewrites \begin{align*} \varsigma^R&=\frac{S_E}{\vartheta} -\frac{k}{h}\int_0^{\infty}\int_{\mathcal{S}^2} \frac{1}{\nu} \big[\log \frac{n(I)}{n(I)+1}-\log \frac{n(B)}{n(B)+1}\big] \sigma_a(B-I)\,d\vec \omega d\nu \\ &\quad -\frac{k}{h}\int_0^{\infty}\int_{\mathcal{S}^2} \frac{1}{\nu} \big[\log \frac{n(I)}{n(I)+1}-\log \frac{n(\tilde I)}{n(\tilde I)+1} \big]\sigma_s(\tilde I-I)\,d\vec \omega d\nu, \end{align*} where we used the hypothesis that the transport coefficients $\sigma_{a,s}$ do not depend on $\vec \omega $. So we obtain finally \begin{equation} \label{ient} \partial_t ( \varrho s +s^R) + \operatorname{div}_x ( \varrho s \vec{u} +\vec q^R) +\operatorname{div}_x (\frac{\vec q}{\vartheta}) =\varsigma+\varsigma^R. \end{equation} and equation \eqref{m16} is replaced in the weak formulation by the inequality \begin{equation} \label{m17} \begin{aligned} &\int_0^T \int_{\Omega} \Big( (\varrho s+s^R) \partial_t \varphi + \varrho s \vec{u} \cdot \nabla_x \varphi + \big(\frac{ \vec q }{\vartheta}+\vec q^R\big) \cdot \nabla_x \varphi \Big) \,dx \,dt \\ &\leq - \int_{\Omega} (\varrho s+s^R)_0 \varphi (0, \cdot) \,dx\\ &\quad - \int_0^T \int_{\Omega} \frac{1}{\vartheta} \Big( \mathbb{S} : \nabla_x \vec{u} - \frac{\vec{q} \cdot \nabla_x \vartheta }{\vartheta} {+\frac{\lambda}{\zeta}|\operatorname{curl}_x \vec B|^2}\Big) \varphi \,dx \,dt \\ &\quad - \frac{k}{h}\int_0^T \int_{\Omega} \Big[\int_0^{\infty}\int_{\mathcal{S}^2} \frac{1}{\nu} \big[\log \frac{n(I)}{n(I)+1}-\log \frac{n(B)}{n(B)+1}\big] \sigma_a(B-I)\,d\vec \omega d\nu \\ &\quad +\int_0^{\infty}\int_{\mathcal{S}^2} \frac{1}{\nu} \big[\log \frac{n(I)}{n(I)+1}-\log \frac{n(\tilde I)}{n(\tilde I)+1}\big] \sigma_s(\tilde I-I)\,d\vec \omega d\nu \Big] \varphi \,dx \,dt \end{aligned} \end{equation} for any $\varphi \in C^\infty_c([0,T) \times \overline{\Omega})$, $\varphi \geq 0$, where the sign of all the terms in the right hand side may be controlled. $\bullet$ Since replacing equation \eqref{i3} by inequality \eqref{m17} would result in a formally under-determined problem, system \eqref{m12}, \eqref{m13}, \eqref{m17} must be supplemented with the \emph{total energy balance} \begin{equation} \label{m18} \begin{aligned} &\int_{\Omega} ( \frac{1}{2} \varrho |\vec{u}|^2 + \varrho e(\varrho, \vartheta)+ {\frac1{2 \zeta}} |\vec{B}|^2 - {\frac12 } \varrho \Psi - {\frac12 }\varrho |\vec \chi\times\vec x|^2 + E^R ) (\tau, \cdot) \,dx \\ & + \int_0^\tau \int\int_{\Gamma_{+}} \int_0^\infty I (t,x,\vec \omega, \nu)\,\vec \omega \cdot \vec{n} \,d\nu \,d \vec \omega \,dS_x \,dt \\ &= \int_{\Omega} \Big( \frac{1}{2 \varrho_0}|(\varrho \vec{u})_0 |^2 + (\varrho e)_0 + {\frac1{2 \zeta}} |\vec{B}_0|^2 - {\frac12 } \varrho_0 \Psi_0 - {\frac12 }\varrho_0 |\vec \chi\times\vec x|^2+ E^R_0 \Big) \,dx, \end{aligned} \end{equation} where \[ E^R_0(x) = \frac{1}{c} \int_{\mathcal{S}^2} \int_0^\infty I(0,x, \vec \omega, \nu) \,d \vec \omega \,d \nu. \] The transport equation \eqref{i4} can be extended to the whole physical space $\mathbb{R}^3$ provided we set $\sigma_a (x, \nu, \vartheta) = \mathbb{I}_\Omega \sigma_a (\nu, \vartheta)$ and $\sigma_s (x, \nu, \vartheta) = \mathbb{I}_\Omega \sigma_s (\nu, \vartheta)$, where $\mathbb{I}_{A}$ is the characteristic function of a set $A$ and take the initial distribution $I_0 (x, \vec \omega, \nu)$ to be zero for $x \in \mathbb{R}^3 \setminus \Omega$. Accordingly, for any fixed $\vec \omega \in \mathcal{S}^2$, equation \eqref{i4} can be viewed as a linear transport equation defined in $(0,T) \times \mathbb{R}^3$, with a right-hand side $S$. With the above mentioned convention, extending $\vec{u}$ to be zero outside $\Omega$, we may therefore assume that both $\varrho$ and $I$ are defined on the whole physical space $\mathbb{R}^3$. \begin{definition} \label{def2.1} \rm We say that $\varrho, \vec{u}, \vartheta, \vec B, I$ is a weak solution of problem \eqref{i1}--\eqref{i6bis} if \begin{gather*} \varrho \geq 0, \; \vartheta > 0 \text{ for a.a. } (t,x) \times \Omega,\quad I \geq 0 \text{ a.a. in } (0,T) \times \Omega \times \mathcal{S}^2 \times (0, \infty), \\ \varrho \in L^\infty(0,T; L^{5/3}(\Omega)), \quad \vartheta \in L^\infty(0,T; L^4(\Omega)), \\ \vec{u} \in L^2(0,T; W^{1,2}_0(\Omega; \mathbb{R}^3)), \quad \vartheta \in L^2(0,T; W^{1,2}(\Omega)), \\ \vec B \in L^2(0,T; W^{1,2}(\Omega; \mathbb{R}^3)),\quad \vec B \cdot \vec n\big|_{\partial\Omega} = 0, \\ I \in L^\infty ((0,T) \times \Omega \times \mathcal{S}^2 \times (0,\infty)),\quad I \in L^\infty(0,T; L^1(\Omega \times \mathcal{S}^2 \times (0,\infty)), \end{gather*} and if $\varrho$, $\vec{u}$, $\vartheta$, $\vec B$, $I$ satisfy the integral identities \eqref{m12}, \eqref{m13}, \eqref{m17}, \eqref{v9}, \eqref{m18}, together with the transport equation \eqref{i4}. \end{definition} The stability result of \cite{dkn} reads as follows \begin{theorem} \label{thm1} Let $\Omega \subset \mathbb{R}^3$ be a bounded $C^{2,\alpha}$ with $\alpha>0$ domain. Assume that the thermodynamic functions $p$, $e$, $s$ satisfy hypotheses \eqref{m1}--\eqref{m6}, and that the transport coefficients $\mu$, $\lambda$, $\kappa$, $\sigma_a$, and $\sigma_s$ comply with \eqref{m7}--\eqref{m10}. Let $\{ \varrho_\varepsilon, \vec{u}_\varepsilon, \vartheta_\varepsilon, \vec B_\varepsilon, I_\varepsilon \}_{\varepsilon > 0}$ be a family of weak solutions to problem \eqref{i1}--\eqref{i12} in the sense of Definition \ref{def2.1} such that \begin{gather} \label{m19} \varrho_\varepsilon(0, \cdot) \equiv \varrho_{\varepsilon,0} \to \varrho_0 \quad \text{in } L^{5/3}(\Omega), \\ \label{m20} \begin{aligned} &\int_{\Omega} \Big( \frac{1}{2} \varrho_\varepsilon |\vec{u}_\varepsilon|^2 + \varrho_\varepsilon e(\varrho_\varepsilon, \vartheta_\varepsilon) + {\frac1{2 \zeta}} |\vec{B}_\varepsilon|^2 - {\frac12 } \varrho_\varepsilon \Psi_\varepsilon \\ &- {\frac12 }\varrho_\varepsilon |\vec \chi\times\vec x|^2 + E_{\varepsilon}^R \Big) (0, \cdot) \,dx \\ &\equiv \int_{\Omega} ( \frac{1}{2 \varrho_{0,\varepsilon}} |(\varrho \vec{u})_{0,\varepsilon} |^2 + (\varrho e)_{0,\varepsilon} + E_{0,\varepsilon}^R + {\frac1{2 \zeta}} |\vec{B}_{0,\varepsilon}|^2 - {\frac12 }\varrho_{\varepsilon,0} |\vec \chi\times\vec x|^2 \\ &\quad - {\frac12 } \varrho_{\varepsilon,0} \Psi_{\varepsilon,0}) \,dx \leq E_0, \end{aligned} \\ \nonumber \int_{\Omega} [\varrho_\varepsilon s(\varrho_\varepsilon, \vartheta_\varepsilon)+s^R(I_\varepsilon)](0, \cdot) \,dx \equiv \int_{\Omega} (\varrho s+s^R)_{0, \varepsilon} \,dx \geq S_0, \\ \nonumber 0 \leq I_\varepsilon(0, \cdot) \equiv I_{0, \varepsilon} (\cdot) \leq I_0,\quad |I_{0, \varepsilon} (\cdot, \nu) | \leq h(\nu) \quad \text{for a certain } h \in L^1(0,\infty). \end{gather} Then \begin{gather*} \varrho_\varepsilon \to \varrho \quad \text{in } C_{\rm weak}([0,T]; L^{5/3}(\Omega)), \\ \vec{u}_\varepsilon \to \vec{u} \quad \text{weakly in } L^2(0,T; W^{1,2}_0 (\Omega; \mathbb{R}^3)), \\ \vartheta_\varepsilon \to \vartheta \quad \text{weakly in } L^2(0,T; W^{1,2}(\Omega)), \\ \vec B_\varepsilon \to \vec B \quad \text{weakly in } L^2(0,T; W^{1,2} (\Omega; \mathbb{R}^3)), \vec B \cdot \vec n\big|_{\partial\Omega} = 0, \\ I_\varepsilon \to I \ \text{weakly-(*) in} \ L^\infty((0,T) \times \Omega \times \mathcal{S}^2 \times (0, \infty) ), \end{gather*} at least for suitable subsequences, where $\{ \varrho, \vec{u}, \vartheta, \vec B, I \}$ is a weak solution of problem \eqref{i1}--\eqref{i12}. \end{theorem} \section{Formal scaling analysis} \label{fsa} To identify the appropriate limit regime we perform a general scaling, denoting by $L_{\rm ref}$, $T_{\rm ref},U_{\rm ref}$, $\rho_{\rm ref}$, $\vartheta_{\rm ref}$, $p_{\rm ref}$, $e_{\rm ref}$, $\mu_{\rm ref}$, $ \lambda _{\rm ref}$, $\ \kappa_{\rm ref}$, the reference hydrodynamical quantities (length, time, velocity, density, temperature, pressure, energy, viscosity, conductivity), by $I_{\rm ref}$, $\nu_{\rm ref}$, $\sigma_{a,{\rm ref}}$, $\sigma_{s,{\rm ref}}$, the reference radiative quantities (radiative intensity, frequency, absorption and scattering coefficients), by $\chi_{\rm ref}$ the reference rotation velocity, and by $\zeta_{\rm ref},\ B_{\rm ref}$ the reference electrodynamic quantities (permeability and magnetic induction). We also assume the compatibility conditions $ p_{\rm ref}\equiv \rho_{\rm ref}e_{\rm ref}$, $\nu_{\rm ref}=\frac{k\vartheta_{\rm ref}}{h}$, $I_{\rm ref}=\frac{2h\nu_{\rm ref}^3}{c^2}$, $\tilde\lambda = \frac{\lambda_{\rm ref}}{L_{\rm ref}U_{\rm ref}}$ and we denote by $Sr:=\frac{L_{\rm ref}}{T_{\rm ref}U_{\rm ref}}$, $Ma:=\frac{U_{\rm ref}}{\sqrt{p_{\rm ref}/\rho_{\rm ref}}}$, $Re:=\frac{U_{\rm ref}\rho_{\rm ref}L_{\rm ref}}{\mu_{\rm ref}}$, $Pe:=\frac{U_{\rm ref}p_{\rm ref}L_{\rm ref}}{\vartheta_{\rm ref}\kappa_{\rm ref}}$, $Fr:=\frac{U_{\rm ref}}{\sqrt{G\rho_{\rm ref}L^2_{\rm ref}}}$, ${\mathcal C}:=\frac{c}{U_{\rm ref}}$, the Strouhal, Mach, Reynolds, P\'eclet, Froude and ``infrarelativistic" dimensionless numbers corresponding to hydrodynamics, by $Ro:=\frac{U_{\rm ref}}{\chi_{\rm ref}L_{\rm ref}}$ the Rossby number, by $Al:=\frac{U_{\rm ref}\rho_{\rm ref}^{1/2}\zeta_{\rm ref}^{1/2}}{B_{\rm ref}}$ the Alfv\'en number and by ${\mathcal L}:=L_{\rm ref}\sigma_{a,{\rm ref}}$, ${\mathcal L}_s:=\frac{\sigma_{s,{\rm ref}}}{\sigma_{a,{\rm ref}}}$, ${\mathcal{P}}:=\frac{2k^4\vartheta_{\rm ref}^4}{h^3c^3\ \rho_{\rm ref}e_{\rm ref}}$, various dimensionless numbers corresponding to radiation. Using these scalings and using carets to symbolize renormalized variables we obtain $$ S=\frac{I_{\rm ref}}{L_{\rm ref}}\ \hat S, $$ where $$ \hat S={\mathcal L}\hat \sigma_a(B(\hat \nu,\hat \vartheta)-\hat I) +{\mathcal L}{\mathcal L}_s\hat \sigma_s\Big(\frac{1}{4\pi} \int_{\mathcal{S}^2} \hat I(\cdot , \vec \omega) \,d \vec \omega -\hat I\Big). $$ Omitting the carets in the following, we obtain first the scaled equation for $I$, in the region $(0,T) \times \Omega \times (0, \infty)\times \mathcal{S}^2$, \begin{equation} \label{i4bis} \frac{Sr}{{\mathcal C}}\ \partial_tI+\vec \omega \cdot \nabla_x I =s={\mathcal L} \sigma_a(B- I) +{\mathcal L}{\mathcal L}_s \sigma_s \Big(\frac{1}{4\pi} \int_{\mathcal{S}^2} I \,d \vec \omega - I\Big), \end{equation} where we used the same notation $B$ for the dimensionless Planck function \[ B(\nu,\vartheta)=\frac{\nu^3}{e^{\nu/\vartheta}-1}. \] Denoting also by $ E^R= \int_{\mathcal{S}^2} \int_0^\infty I\,d \nu\,d\vec \omega$, the (renormalized) radiative energy, by \[ \vec F^R = \int_{\mathcal{S}^2} \int_0^\infty I \vec \omega{\rm d}\nu \,d \vec \omega, \] the renormalized radiative momentum, by $s_E = \int_{\mathcal{S}^2} \int_0^\infty s \,d \nu \,d \vec \omega$, the renormalized radiative energy source, by $\vec s^R=-\int_0^{\infty}\int_{\mathcal{S}^2} \nu^2 [ n\log n-(n+1)\log(n+1)]\,d\vec \omega d\nu$, the renormalized radiative entropy with $n=n(I)=\frac{I}{\nu^3}$, by \[ \vec q^R=-\int_0^{\infty}\int_{\mathcal{S}^2} \nu^2 [ n\log n-(n+1)\log(n+1)]\vec \omega\,d\vec \omega d\nu, \] the renormalized radiative entropy flux, and taking the first moment of \eqref{i4bis} with respect to $\vec \omega$, we obtain first an equation for $E^R$ \begin{equation} \label{erfr} \frac{Sr}{{\mathcal C}}\ \partial_t E^R+\nabla_x \cdot\vec F^R=s_E. \end{equation} The continuity equation is now \begin{equation} \label{i1bis} Sr\ \partial_t \varrho + \operatorname{div}_x (\varrho \vec{u}) = 0, \end{equation} and the momentum equation reads \begin{equation} \label{i2bis} \begin{aligned} & Sr\ \partial_t (\varrho \vec{u}) + \operatorname{div}_x (\varrho \vec{u} \otimes \vec{u}) + \frac{1}{Ma^2}\ \nabla_x p(\varrho,\vartheta) + \frac{2}{Ro}\varrho \vec \chi \times \vec{u} \\ & = \frac{1}{Re}\ \operatorname{div}_x \mathbb{S} +\frac{1}{Fr^2}\varrho \nabla \Psi +\frac{1}{2 Ro^2}\varrho \nabla_x |\vec \chi \times \vec x|^2 + \frac{1}{Al^2}\vec{j}\times \vec B. \end{aligned} \end{equation} The balance of internal energy rewrites \begin{align*} &Sr \partial_t \big( \varrho e +\frac1{\mathcal C} E^R\big) + \operatorname{div}_x \big( \varrho e \vec{u} + \vec F^R\big)+\frac{1}{Pe}\ \operatorname{div}_x \vec q\\ &= \frac{Ma^2}{Re}\ \mathbb{S}:\nabla_x \vec{u} -p\operatorname{div}_x \vec{u} + Sr\frac{Ma^2}{Al^2}\vec j\cdot \vec E, \end{align*} and we obtain the balance of matter (fluid) entropy \begin{equation} \label{i3bis} Sr \partial_t (\varrho s ) + \operatorname{div}_x (\varrho s \vec{u}) + \frac{1}{Pe}\ \operatorname{div}_x ( \frac{\vec{q}}{\vartheta}) = \varsigma, \end{equation} with \[ \varsigma= \frac{1}{\vartheta} \Big( \frac{Ma^2}{Re} \mathbb{S}:\nabla_x \vec{u} -\frac{1}{Pe}\frac{\vec q\cdot \nabla_x\vartheta}{\vartheta} +\frac{Ma^2}{Al^2}\frac{\lambda}{\zeta}|\operatorname{curl}_x \vec B|^2\Big)-\frac{S_E}{\vartheta}, \] and the balance of radiative entropy \begin{equation} \label{i3ter} \frac{Sr}{{\mathcal C}}\ \partial_t s^R+ \operatorname{div}_x \vec q^R= \varsigma^R, \end{equation} with \begin{align*} \varsigma^R &=\mathcal{P}\mathcal{L}\int_0^{\infty}\int_{\mathcal{S}^2} \frac{1}{\nu} \big[\log \frac{n(I)}{n(I)+1} -\log \frac{n(B)}{n(B)+1}\big] \sigma_a(I-B)\,d\vec \omega d\nu \\ &\quad +\mathcal{P}\mathcal{L}\mathcal{L}_s\int_0^{\infty} \int_{\mathcal{S}^2} \frac{1}{\nu} \big[\log \frac{n(I)}{n(I)+1} -\log \frac{n(\tilde I)}{n(\tilde I)+1}\big]\sigma_s(I-\tilde I)\, d\vec \omega d\nu+\frac{S_E}{\vartheta}. \end{align*} The scaled equation for the electromagnetic field is \begin{equation} \label{i5ter} Sr\partial_t \vec{B} + \operatorname{curl}_x( \vec{B} \times \vec{u} ) + \operatorname{curl}_x( \lambda \operatorname{curl}_x \vec{B} )= 0. \end{equation} The scaled equation for total energy gives finally the total energy balance \begin{equation} \label{toten} \begin{aligned} &Sr \frac{d}{dt} \int_{\Omega} \Big(\frac{Ma^2}{2}\ \varrho |\vec{u}|^2+ \varrho e +\frac{1}{{\mathcal C}}\ E^R +\frac{Ma^2}{2 Al^2}\frac{1}{\zeta}|\vec B|^2 - {\frac12 } \frac{Ma^2}{Fr^2}\varrho \Psi \\ &- {\frac12}\frac{Ma^2}{Ro^2}\varrho |\vec \chi\times\vec x|^2 \Big)\, dx + \int_0^{\infty}\int_{\Gamma_+} I \vec \omega \cdot \vec n \, d\Gamma_+ d\nu=0. \end{aligned} \end{equation} In the sequel we analyze the asymptotic regime defined by \[ Ma=\varepsilon, \quad Al=\varepsilon,\quad Fr=\varepsilon^{1/2},\quad {\mathcal C}=\varepsilon^{-1}, \quad Pe = \varepsilon^2 \] where $\varepsilon>0$ is small and we put $Sr=1$, $Re=1$, $Ro=1$, ${\mathcal{P}}=1$, ${\mathcal L}={\mathcal L}_s=1$ in the previous system. Plugging this scaling into the previous system gives \begin{gather} \label{radeq} \varepsilon\partial_t I+\vec \omega \cdot \nabla_x I = \sigma_a(B- I) + \sigma_s(\frac{1}{4\pi} \int_{\mathcal{S}^2} I \,d \vec \omega - I),\\ \label{masseq} \partial_t \varrho + \operatorname{div}_x (\varrho \vec{u}) = 0, \\ \label{momeq} \begin{aligned} &\partial_t (\varrho \vec{u}) + \operatorname{div}_x (\varrho \vec{u} \otimes \vec{u}) + \frac{1}{\varepsilon^2}\ \nabla_x p(\varrho,\vartheta) + 2\varrho \vec \chi \times \vec{u}\\ &= \operatorname{div}_x \mathbb{S}+\frac{1}{\varepsilon}\varrho \nabla \Psi + \frac12 \varrho \nabla_x |\vec \chi \times \vec x|^2 + \frac{1}{\varepsilon^2}\vec{j}\times \vec B, \end{aligned} \\ \label{eninteq} \begin{aligned} &\partial_t ( \varrho e+\varepsilon E^R) + \operatorname{div}_x ( \varrho e \vec{u} + \vec F^R)+ \frac{1}{\varepsilon^2}\operatorname{div}_x \vec q \\ &=\varepsilon^2 \mathbb{S}:\nabla_x \vec{u} -p\operatorname{div}_x \vec{u}+ \vec j\cdot \vec E, \end{aligned} \\ \label{enteq} \partial_t (\varrho s+\varepsilon s^R) + \operatorname{div}_x (\varrho s \vec{u}+ \vec q^R) + \ \frac{1}{\varepsilon^2}\operatorname{div}_x ( \frac{\vec{q}}{\vartheta})\geq \varsigma_{\varepsilon}, \end{gather} with \begin{gather} \nonumber \begin{aligned} \varsigma_{\varepsilon} &=\frac{1}{\vartheta} ( \varepsilon^2 \mathbb{S}:\nabla_x \vec{u} -\frac{\vec q\cdot \nabla_x\vartheta}{\varepsilon^2\vartheta}+\frac{\lambda}{\zeta}|\operatorname{curl}_x \vec B|^2) \\ &\quad +\int_0^{\infty}\int_{\mathcal{S}^2} \frac{1}{\nu} \big[\log \frac{n(I)}{n(I)+1}-\log \frac{n(B)}{n(B)+1}\big] \sigma_a(I-B)\,d\vec \omega d\nu \\ &\quad +\int_0^{\infty}\int_{\mathcal{S}^2} \frac{1}{\nu} \big[\log \frac{n(I)}{n(I)+1}-\log \frac{n(\tilde I)}{n(\tilde I)+1}\big] \sigma_s(I-\tilde I)\,d\vec \omega d\nu, \end{aligned} \\ \label{magneq} \partial_t \vec{B} + \operatorname{curl}_x( \vec{B} \times \vec{u} ) + \operatorname{curl}_x( \lambda\ \operatorname{curl}_x \vec{B} )= 0, \end{gather} and finally \begin{equation} \label{etoteq} \begin{aligned} &\frac{d}{dt} \int_{\Omega} \Big(\frac{1}{2} \varepsilon^2\varrho|\vec{u}|^2 + \varrho e +\varepsilon E^R +\frac{1}{2\zeta}|\vec B|^2 - {\frac12 } \varepsilon\varrho \Psi - {\frac12 }\varrho\varepsilon^2 |\vec \chi\times\vec x|^2\Big)\,dx \\ &+ \int_0^{\infty}\int_{\Gamma_+} \vec \omega \cdot \vec n I \, d\Gamma_+ d\nu=0. \end{aligned} \end{equation} To compute the limit system, we consider now the formal expansions \begin{equation} { ( I,\varrho,\vec{u},\vartheta,p,\vec B)=( \breve{I_0},\breve{\varrho_0},\breve{u_0},\breve{\vartheta_0}, \breve{p_0},\breve{B_0})+\varepsilon ( I_1,\varrho_1,\vec{u}_1,\vartheta_1,p_1,\vec B_1) +O(\varepsilon^2)}. \label{scal} \end{equation} $\bullet$ We first observe from \eqref{momeq} that $\breve{\varrho_0}=const := \overline\varrho$ and $\breve{\vartheta_0}=const := \overline\vartheta$, moreover \begin{equation} \nabla_x p_1=\overline{\varrho} \nabla_x\Psi(\overline{\varrho}). \label{rt02} \end{equation} Let us fix the constants in the Neumann problem for perturbations of the temperature \begin{equation} \int_\Omega\vartheta_i\,dx = 0\quad\text{for any } i\ge1. \label{zeromeantemp} \end{equation} From \eqref{masseq} we derive the incompressibility condition \begin{equation} \operatorname{div}_x \breve{u_0}=0, \label{wed1bis} \end{equation} and \begin{equation} \partial_t \varrho_1+\operatorname{div}_x ( \overline{\varrho} \vec{u}_1+\varrho_1 \breve{u_0})=0. \label{v1} \end{equation} $\bullet$ From \eqref{radeq} we obtain now two stationary linear transport equations for the two moments $\breve{I_0}$ and $I_1$ \begin{gather} \label{radeq0} \vec \omega \cdot \nabla_x \breve{I_0} = \sigma_{a,0}(B_0- \breve{I_0}) + \sigma_{s,0}(\tilde{\breve{I_0}} - \breve{I_0}), \\ \label{radeq1} \begin{aligned} \vec \omega \cdot \nabla_x I_1 &= \sigma_{a,0}(\partial_{\vartheta}B_0\vartheta_1- I_1) +\partial_{\vartheta}\sigma_{a,0}(B_0- \breve{I_0})\vartheta_1 \\ &\quad +\partial_{\vartheta} \sigma_{s,0}(\tilde{\breve{I_0}} - \breve{I_0})\vartheta_1 +\sigma_{s,0}(\tilde I_1 - I_1), \end{aligned} \end{gather} where $\tilde I:=\frac{1}{4\pi} \int_{\mathcal{S}^2} I \,d \vec \omega$, $\sigma_{a,0}=\sigma_a(\nu,\breve{\vartheta_0})$, $\sigma_{s,0}=\sigma_s(\nu,\breve{\vartheta_0})$ and $B_{0}=B(\nu,\breve{\vartheta_0})$. $\bullet$ The limit momentum equation reads \begin{equation} \label{momeqbis} \overline{\varrho}(\partial_t \breve{u_0} + \operatorname{div}_x ( \breve{u_0} \otimes \breve{u_0}) + \nabla_x \Pi + 2\overline{\varrho}\vec \chi \times \breve{u_0} = \operatorname{div}_x \mathbb{S}(\breve{u_0}) + \frac{1}{\zeta}\operatorname{curl}_x \vec B_1\times \vec B_1 +\vec F, \end{equation} where $\mu_0=\mu(\breve{\vartheta_0})$ is used in $\mathbb{S}(\breve{u}_0)$, $\vec F=\varrho_1 \nabla_x \Psi(\overline{\varrho})$ and $\Pi$ is an effective pressure for which it holds $\nabla_x\Pi = \frac12\overline{\varrho} \nabla_x |\vec \chi \times \vec x|^2+ \overline{\varrho}\nabla_x\Pi(\varrho_1)+p_{\varrho,\varrho}(\overline{\varrho},\overline{\vartheta})\varrho_1\nabla_x\varrho_1$. Here we set $\vartheta_1=0$ which is consistent with the $O(\varepsilon^{-1})-$order of the internal energy equation \eqref{eninteq} and the additional zero mean of $\vartheta - \breve{\vartheta_0}$ requirement. $\bullet$ The limit magnetic field $\vec B_1$ solves \begin{equation} \label{mhdeq} \partial_t \vec{B}_1 + \operatorname{curl}_x( \vec{B}_1 \times \breve{u}_0 ) + \operatorname{curl}_x( \overline{\lambda}\ \operatorname{curl}_x \vec{B}_1 )= 0, \end{equation} for $\overline{\lambda}=\lambda(\breve{\vartheta}_0)$. $\bullet$ At the lowest order ($O(\varepsilon^0)$) the energy equation \eqref{eninteq} gives \begin{equation} \label{enerzero} \overline{\kappa}\Delta\vartheta_2 = {s_E}_0 \end{equation} where $-{s_E}_0 = \int_0^\infty\int_{\mathcal{S}^2}\sigma_{a,0}(\breve{I}_0-B_0) \,d\vec\omega\,d\nu$ and $\overline{\kappa} = \kappa(\overline{\vartheta})$. $\bullet$ At the order ($O(\varepsilon)$) we simplify the energy equation \eqref{eninteq}. Observing that from \eqref{rt02} we have \begin{equation} \label{constraint} \partial_{\varrho}p(\overline{\varrho},\overline{\vartheta}) D \varrho_1 + \overline{\varrho}\breve{u}_0\cdot\nabla_x\Psi(\overline{\varrho})=0, \end{equation} where $D:=\partial_t +\breve{u}_0\cdot\nabla_x$, and from \eqref{v1} \[ \overline{\varrho}\operatorname{div}_x \vec u_1=-D\varrho_1, \] and after \eqref{radeq1} \[ {S_E}_1=-\int_0^{\infty}\int_{\mathcal{S}^2}\sigma_{a,0}I_1\,d\vec \omega\,d\nu, \] and simplifying by \eqref{i5bis} we end up with \[ \partial_t \varrho_1+\operatorname{div}_x( \varrho_1 \breve{u}_0) = -\overline{\alpha} (\overline{\kappa}\Delta\vartheta_3 +\int_0^{\infty}\int_{\mathcal{S}^2}\sigma_{a,0}I_1\,d\vec \omega\,d\nu) \] where $\overline{\alpha} :=\frac{\overline{\varrho}}{\overline{\vartheta}}\partial_{\theta}p(\overline{\varrho},\overline{\vartheta})$. Putting $\vec U=\breve{u}_0$, $\Theta=\vartheta_3$, $\vec B=\vec B_1$, $\overline{\varrho}=\breve{\varrho}_0$, $\overline{\vartheta}=\breve{\vartheta}_0$, $\overline{\mu}=\mu(\breve{\vartheta}_0)$, $\sigma_a=\sigma_{a,0}$, $\sigma_s=\sigma_{s,0}$, $ B=B_0$, $\mathbb{D}(\vec U)= \frac{1}{2}(\nabla \breve{u}_0+\nabla^T \breve{u}_0)$, and \[ G= \frac{\int_0^\infty \int_{\mathcal{S}^2} \sigma_{a,0}(\breve{I}_0-B_0)\,d \vec\omega\,d\nu}{\overline\kappa} \] we observe that the solution of the equation \eqref{radeq0} is up to the boundary condition \eqref{i11}$_2$ $\breve{I}_0 = B_0$ which in turn entails that the equation for $\vartheta_2$ turns in $\Omega$ into the Laplace homogeneous equation ($G = 0$) and therefore $\vartheta_2 = 0$ and we obtain the limit system in $ (0,T) \times \Omega$ \begin{gather} \label{bous1} \operatorname{div}_x \vec U = 0, \\ \label{bous2} \overline{\varrho}(\partial_t \vec U + \operatorname{div}_x ( \vec U \otimes \vec U)) + \nabla_x \Pi = \operatorname{div}_x ( 2\overline{\mu}\ \mathbb{D}(\vec U)) + \frac{1}{\zeta}\operatorname{curl}_x \vec B\times \vec B+\vec F, \\ \label{bous1bis} \partial_t \vec{B} + \operatorname{curl}_x( \vec{B} \times \vec U ) + \operatorname{curl}_x( \overline{\lambda}\ \operatorname{curl}_x \vec{B} )= 0, \\ \label{bous1ter} \operatorname{div}_x \vec B = 0, \\ \label{bous3} -\Delta\Theta =\frac{1}{\overline{\alpha}\overline{\kappa}}\vec U\cdot\nabla_x \tilde r - \frac{1}{\overline{\kappa}}\int_0^{\infty}\sigma_{a}\int_{\mathcal{S}^2}I_1\,d\vec \omega\,d\nu + \tilde{h}(t), \\ \label{bous5} \vec \omega \cdot \nabla_x I_1 = -\sigma_{a} I_1 +\sigma_{s}(\tilde I_1 - I_1), \end{gather} together with the \emph{Boussinesq relation} \eqref{rt02} \begin{equation} \label{bous5bis} \nabla_x \tilde{r}=\frac{\overline{\varrho} \nabla_x\Psi(\overline{\varrho})}{\partial_{\varrho}p(\overline{\varrho},\overline{\vartheta})}, \end{equation} where $\tilde r := \varrho_1-\overline\varrho$ and $\tilde{h}$ is an undetermined function which allows satisfaction of \eqref{bous10}$_2$. We finally consider the boundary conditions \begin{equation} \label{bous6} \vec U|_{\partial \Omega} = 0, \quad \nabla \Theta \cdot \vec{n} |_{\partial \Omega} = 0, \quad \vec B \cdot \vec{n}|_{\partial \Omega} = 0, \quad \operatorname{curl}_x\vec B \times \vec{n}|_{\partial \Omega} = 0 \end{equation} for \eqref{bous1}-\eqref{bous3} and \begin{equation} \label{bous8} I_1 (x,\nu, \vec \omega) = 0 \quad \text{for } x \in \partial \Omega, \quad \vec \omega \cdot \vec{n} \leq 0 \end{equation} for \eqref{bous5}, and the initial conditions \begin{equation} \label{bous9} \vec U|_{t=0} = \vec U_0,\quad \vec B|_{t=0} = \vec B_0. \end{equation} Moreover, we endow system \eqref{bous1}--\eqref{bous5bis} with the additional conditions \begin{equation} \label{bous10} \operatorname{div}_x \vec B_0 = 0, \quad \int_{\Omega} \Theta \,dx= 0. \end{equation} For this system we have the following existence result (see the Appendix for a short proof) \begin{theorem} \label{ROB} Let $\Omega \subset \mathbb{R}^3$ be a bounded Lipschitz domain. For any $T>0$ the initial-boundary value problem \eqref{bous1}-\eqref{bous10} has at least a weak solution $(\vec U,\Theta,\vec B,I_1)$ such that \begin{enumerate} \item \[ \vec U\in L^{\infty}(0,T;{\mathcal H}(\Omega))\cap L^2(0,T;{\mathcal U}(\Omega)), \quad \vec B\in L^{\infty}(0,T;{\mathcal V}(\Omega))\cap L^2(0,T;{\mathcal W}(\Omega)), \] with ${\mathcal H}(\Omega)=\{ \vec U\in L^2(\Omega;\mathbb{R}^3), \operatorname{div}_x \vec U=0\ \text{in } \Omega,\ \vec U\big|_{\partial \Omega}=0\}$, ${\mathcal U}(\Omega)={\mathcal H}(\Omega)\cap W^{1,2}_0(\Omega;\mathbb{R}^3))$, ${\mathcal V}(\Omega)=\big\{\vec b\in L^2(\Omega;\mathbb{R}^3),\ \operatorname{div}_x \vec b=0, \vec b\cdot\vec n\big|_{\partial \Omega}=0\big\}$ and $ {\mathcal W}(\Omega)={\mathcal V}(\Omega)\cap W^{1,2}_0(\Omega;\mathbb{R}^3)$, \item \[ \Theta\in L^{\infty}((0,T; W^{2,2}(\Omega)) \cap L^{2}((0,T; W^{q,2}(\Omega))\quad\text{for any $q<\frac{5}{2}$}, \] \item \[ I_1 \in L^{\infty}((0,T)\times \Omega\times\mathcal{S}^2\times\mathbb{R}_+), \] with \[ \vec \omega\cdot\nabla_x I_1 \in L^p((0,T)\times \Omega\times\mathcal{S}^2\times\mathbb{R}_+), \] for any $p>1$ and any $\vec \omega\in\mathcal{S}^2$. \end{enumerate} \end{theorem} The remaining part of the paper is devoted to the proof of the convergence of the primitive system \eqref{i1}--\eqref{i12} to the target system \eqref{bous1}--\eqref{bous10}. \section{Global existence for the primitive system and uniform estimates} \label{ge} For system \eqref{i1}--\eqref{i12} we prepare the initial data as follows \begin{equation} \label{data1} \begin{gathered} \varrho(0,\cdot)=\varrho_{0,\varepsilon}=\overline{\varrho}+\varepsilon\varrho_{0,\varepsilon}^{(1)},\\ \vec u(0,\cdot)=\vec u_{0,\varepsilon},\\ \vartheta(0,\cdot)=\vartheta_{0,\varepsilon}=\overline{\vartheta}+\varepsilon^3\vartheta_{0, \varepsilon}^{(3)},\\ I(0,\cdot,\cdot,\cdot)=I_{0,\varepsilon} =\overline{I}+\varepsilon I_{0,\varepsilon}^{(1)},\\ \vec B(0,\cdot)={B_{0,\varepsilon}}=\varepsilon \vec B_{0,\varepsilon}^{(1)}, \end{gathered} \end{equation} where $\overline{\varrho}>0$, $\overline{\vartheta}>0$, $\overline{I}>0$ are spacetime constants and $\int_{\Omega}\varrho_{0,\varepsilon}^{(1)}\,dx=0 =\int_{\Omega}\vartheta_{0,\varepsilon}^{(3)}\,dx$ for any $\varepsilon>0$. As in \cite{FEINOV}, for any locally compact Hausdorff metric space $X$ we denote by ${\mathcal M}(X)$ the set of signed Borel measures on $X$ and by ${\mathcal M}^+(X)$ the cone of non-negative elements of ${\mathcal M}(X)$. From Theorem \ref{thm1} we obtain immediately (by combining the approximating schemes introduced in \cite{DFN} and \cite{DF}) the existence of a weak solution $(\varrho_{\varepsilon},\vec u_{\varepsilon},\vartheta_{\varepsilon}, I_{\varepsilon}, \vec B_{\varepsilon})$ to the radiative MHD system \eqref{i1}--\eqref{i12}. \begin{theorem} \label{thm2} Let $\Omega \subset \mathbb{R}^3$ be a bounded $C^{2,\alpha}$ with $\alpha>0$ domain. Assume that the thermodynamic functions $p$, $e$, $s$ satisfy hypotheses \eqref{m1}--\eqref{m6}, and that the transport coefficients $\mu$, $\lambda$, $\kappa$, $\sigma_a$, $\sigma_s$ and the equilibrium function $B$ comply with \eqref{m7}--\eqref{m10}. Let the initial data $(\varrho_{0,\varepsilon},\vec u_{0,\varepsilon},\vartheta_{0,\varepsilon}, I_{0,\varepsilon}, \vec B_{0,\varepsilon})$ be given by \eqref{data1}, where $(\varrho_{0,\varepsilon}^{(1)}, \vartheta_{0,\varepsilon}^{(3)}, I_{0,\varepsilon}^{(1)}, \vec B_{0,\varepsilon}^{(1)})$ are uniformly bounded measurable functions. Then for any $\varepsilon>0$ small enough (to maintain positivity of $\varrho_{0,\varepsilon}$ and $\vartheta_{0,\varepsilon}$), there exists a weak solution $(\varrho_{\varepsilon},\vec u_{\varepsilon},\vartheta_{\varepsilon}, I_{\varepsilon},\vec B_{\varepsilon})$ to the radiative Navier-Stokes system \eqref{i1}--\eqref{i9} for $(t,x,\vec \omega,\nu)\in (0,T)\times \Omega \times \mathcal{S}^2 \times \mathbb{R}_+$, supplemented with the boundary conditions \eqref{i11}--\eqref{i12} and the initial conditions \eqref{data1}. More precisely we have \begin{equation} \label{mass} \begin{aligned} &\int_0^T\int_{\Omega} \varrho_{\varepsilon} b(\varrho_{\varepsilon}) (\partial_t\phi+\vec u_{\varepsilon}\cdot \nabla_x \phi)\,dx\,dt \\ &=\int_0^T\int_{\Omega} \beta(\varrho_{\varepsilon}) \operatorname{div}_x u_{\varepsilon}\ \phi\,dx\,dt -\int_{\Omega} \varrho_{0,\varepsilon} b(\varrho_{0,\varepsilon}) \phi(0,\cdot)\,dx, \end{aligned} \end{equation} for any $\beta$ such that $\beta\in (L^{\infty}\cap C)([0,\infty))$, $b(\varrho)=b(1)+\int_1^{\varrho}\frac{\beta(z)}{z^2}\,dz$ and any $\phi\in C^{\infty}_c([0,T)\times \overline{\Omega})$, \begin{equation} \label{mom} \begin{aligned} &\int_0^T\int_{\Omega} \Big(\varrho_{\varepsilon} \vec u_{\varepsilon}\cdot \partial_t\vec{\varphi} +\varrho_{\varepsilon} \vec u_{\varepsilon}\otimes \vec u_{\varepsilon}:\nabla_x \vec{\varphi} +\frac{p_{\varepsilon}}{\varepsilon^2}\ \operatorname{div}_x \vec{\varphi} -2\varrho_{\varepsilon} \vec \chi \times \vec{u}_{\varepsilon} \cdot \vec{\varphi}\Big)\,dx\,dt\\ &=\int_0^T\int_{\Omega} \Big(\mathbb{S}_{\varepsilon}:\nabla_x \vec{\varphi} -\frac{1}{\varepsilon} \varrho_{\varepsilon} \nabla_x \Psi_{\varepsilon} \cdot \vec{\varphi} -\frac{1}{\varepsilon^2}(\vec{j}_{\varepsilon} \times \vec{B}_{\varepsilon}) \cdot \vec{\varphi} \\ &\quad -\frac12\varepsilon^2\varrho_{\varepsilon} \nabla_x |\vec \chi \times \vec x|^2 \cdot \vec{\varphi}\Big)\,dx\,dt -\int_{\Omega} \varrho_{0,\varepsilon} \vec u_{0,\varepsilon} \cdot \vec{\varphi}(0,\cdot)\,dx, \end{aligned} \end{equation} for any $\vec{\varphi}\in C^{\infty}_c([0,T)\times \overline{\Omega};\mathbb{R}^3)$ with $p_{\varepsilon}=p(\varrho_{\varepsilon},\vartheta_{\varepsilon})$, $\mathbb{S}_{\varepsilon}=\mathbb{S}(\vec u_{\varepsilon},\vartheta_{\varepsilon})$, and $\vec{j}_{\varepsilon}=\frac{1}{\zeta}\operatorname{curl}_x \vec B_{\varepsilon}$, \begin{equation} \label{enertot} \begin{aligned} &\int_{\Omega} (\frac{\varepsilon^2}{2}\ \varrho_{\varepsilon} |\vec u_{\varepsilon}|^2 +\varrho_{\varepsilon}e_{\varepsilon}+\varepsilon E^R_{\varepsilon} +\frac{1}{2\zeta}|\vec B_{\varepsilon}|^2 - {\frac12 } \varepsilon\varrho_{\varepsilon} \Psi_{\varepsilon} - {\frac12 }\varrho_{\varepsilon}\varepsilon^2 |\vec \chi\times\vec x|^2 )\,dx\,dt\\ & +\int_0^T \int_0^{\infty}\int_{\Gamma_+} \vec \omega\cdot \vec n_x I_{\varepsilon}(t,x,\vec \omega,\nu)\,d\Gamma_{+}\,d\nu\,dt \\ &=\int_{\Omega} \Big(\frac{\varepsilon^2}{2}\ \varrho_{0,\varepsilon} |\vec u_{0,\varepsilon}|^2 +\varrho_{0,\varepsilon}e_{0,\varepsilon}+\varepsilon E^R_{0,\varepsilon} +\frac{1}{2\zeta}|\vec B_{0,\varepsilon}|^2- {\frac12 } \varepsilon\varrho_{0,\varepsilon} \Psi_{0,\varepsilon} \\ &\quad - {\frac12 }\varepsilon^2\varrho_{0,\varepsilon} |\vec \chi\times\vec x|^2\Big)\,dx, \end{aligned} \end{equation} for a.a. $t\in(0,T)$ with $e_{\varepsilon}=e(\varrho_{\varepsilon},\vartheta_{\varepsilon})$, $\Psi_{\varepsilon}=\Psi(\varrho_{\varepsilon})$, $\Psi_{0,\varepsilon}=\Psi(\varrho_{0,\varepsilon})$ and \newline $E^R_{\varepsilon}(t,x)=\int_0^{\infty}\int_{\mathcal{S}^2} I_{\varepsilon}(t,x,\vec \omega,\nu)\,d\vec \omega\,d\nu$ \begin{equation} \label{magfield} \begin{aligned} &\int_0^T \int_{\Omega} \Big( \vec{B}_{\varepsilon} \cdot \partial_t \vec{\varphi} - (\vec{B}_{\varepsilon} \times \vec{u}_{\varepsilon} + \lambda_{\varepsilon} \operatorname{curl}_x \vec{B}_{\varepsilon}) \cdot \operatorname{curl}_x \vec{\varphi} \Big) \,dx \,dt \\ &+ \int_{\Omega} \vec{B}_{0,\varepsilon} \cdot \vec{\varphi} (0, \cdot) \,dx=0, \end{aligned} \end{equation} for any vector field $\vec{\varphi} \in \mathcal{D}([0,T) \times \mathbb{R}^3, \mathbb{R}^3)$, with $\lambda_{\varepsilon}=\lambda(\vartheta_{\varepsilon})$. \begin{equation} \label{entro} \begin{aligned} &\int_0^T \int_{\Omega} \Big(( \varrho_{\varepsilon} s_{\varepsilon}+\varepsilon s^R_{\varepsilon}) \partial_t \varphi +(\varrho_{\varepsilon} s_{\varepsilon}\vec{u}_{\varepsilon} +{{\vec q}^R}_{\varepsilon}) \cdot \nabla_x \varphi \Big) \,dx \,dt \\ &+\int_0^T \int_{\Omega}\frac{ {\vec q}_{\varepsilon} }{\varepsilon^2\vartheta_{\varepsilon}} \cdot \nabla_x \varphi \,dx \,dt +\langle\varsigma^m_{\varepsilon}+\varsigma^R_{\varepsilon} ; \varphi\rangle_{[{\mathcal M};C]([0,T)\times \overline{\Omega})} \\ &= - \int_{\Omega} ( ((\varrho s)_{0,\varepsilon} +\varepsilon s^R_{0,\varepsilon}) \varphi (0, \cdot))\,dx, \end{aligned} \end{equation} where \[ \varsigma^m_{\varepsilon} \geq \frac{1}{\vartheta_{\varepsilon}} ( \varepsilon^2\mathbb{S}_{\varepsilon} : \nabla_x \vec{u}_{\varepsilon} - \frac{\vec{q}_{\varepsilon} \cdot \nabla_x \vartheta_{\varepsilon} }{\varepsilon^2\vartheta_{\varepsilon}} +\frac{\lambda_{\varepsilon}}{\zeta}|\operatorname{curl}_x \vec B_{\varepsilon}|^2), \] and \begin{align*} \varsigma^R_{\varepsilon} &\geq \int_0^{\infty}\int_{\mathcal{S}^2} \frac{1}{\nu} \big[\log \frac{n(I_{\varepsilon})}{n(I_{\varepsilon})+1} -\log \frac{n(B_{\varepsilon})}{n(B_{\varepsilon})+1}\big] {\sigma_a}_{\varepsilon}(B_{\varepsilon}-I_{\varepsilon})\,d\vec \omega d\nu \\ &\quad +\int_0^{\infty}\int_{\mathcal{S}^2} \frac{1}{\nu} \big[\log \frac{n(I_{\varepsilon})}{n(I_{\varepsilon})+1} -\log \frac{n(\tilde I_{\varepsilon})}{n(\tilde I_{\varepsilon})+1}\big] {\sigma_s}_{\varepsilon}(\tilde I_{\varepsilon}-I_{\varepsilon})\,d\vec \omega d\nu, \end{align*} for any $\varphi\in C^{\infty}_c([0,T)\times \overline{\Omega})$ with $\varsigma^m_{\varepsilon}\in {\mathcal M}^+([0,T)\times \overline{\Omega})$ and $\varsigma^R_{\varepsilon}\in {\mathcal M}^+([0,T)\times \overline{\Omega})$, and with ${\sigma_a}_{\varepsilon}=\sigma_a(\nu,\vartheta_{\varepsilon})$, ${\sigma_s}_{\varepsilon}=\sigma_s(\nu,\vartheta_{\varepsilon})$, $B_{\varepsilon}=B(\nu,\vartheta_{\varepsilon})$, ${\vec q}_{\varepsilon}=\kappa(\vartheta_{\varepsilon})\nabla_x \vartheta_{\varepsilon}$, $s_{\varepsilon}=s(\varrho_{\varepsilon},\vartheta_{\varepsilon})$, $s^R_{\varepsilon}=s^R(I_{\varepsilon})$, ${{\vec q}^R}_{\varepsilon}={\vec q}^R(I_{\varepsilon})$ and $\tilde I_{\varepsilon}:=\frac{1}{4\pi} \int_{\mathcal{S}^2} I_{\varepsilon}(t,x,\vec \omega, \nu) \,d \vec \omega$, \begin{align} &\int_0^T\int_{\Omega}\int_0^{\infty}\int_{\mathcal{S}^2} (\varepsilon\partial_t\psi +\vec \omega \cdot \nabla_x \psi)I_{\varepsilon}\,d\vec \omega\,d\nu\,dx\,dt \nonumber \\ &+\int_0^T \int_{\Omega}\int_0^{\infty}\int_{\mathcal{S}^2} \big[{\sigma_a}_{\varepsilon}(B_{\varepsilon}-I_{\varepsilon}) +{\sigma_s}_{\varepsilon}(\tilde I_{\varepsilon}-I_{\varepsilon})\big] \psi\,d\vec \omega\,d\nu\,dx\,dt \label{transfer} \\ &=-\int_{\Omega}\int_0^{\infty}\int_{\mathcal{S}^2} \varepsilon I_{0,\varepsilon}\psi(0,x,\vec \omega,\nu)\,d\vec \omega\,d\nu\,dx +\int_0^T \int_{\Gamma_+} \int_0^{\infty} I_{\varepsilon} \vec \omega\cdot \vec n_x\psi\,d\Gamma_{+}\,d\nu\,dt, \nonumber \end{align} for any $\psi\in C^{\infty}_c([0,T)\times \overline{\Omega}\times \mathcal{S}^2\times \mathbb{R}_+)$. \end{theorem} \subsection{Uniform estimates} We recall from \cite{FEINOV} the necessary definitions in the formalism of essential and residual sets (see \cite{DuNe}). Given three numbers $\overline{\varrho}\in\mathbb{R}_+$, $\overline{\vartheta}\in\mathbb{R}_+$ and $\overline{E}\in\mathbb{R}_+$ we define ${\mathcal O}_{\rm ess}^H$ the set of hydrodynamical essential values \begin{equation} \label{essH} {\mathcal O}_{\rm ess}^H:=\big\{ (\varrho,\vartheta)\in\mathbb{R}^2\ :\ \frac{\overline{\varrho}}{2}<\varrho<2\overline{\varrho}, \frac{\overline{\vartheta}}{2}<\vartheta<2\overline{\vartheta} \big\}, \end{equation} and ${\mathcal O}_{\rm ess}^R$ the set of radiative essential values \begin{equation} \label{essR} {\mathcal O}_{\rm ess}^R:=\big\{ E^R\in\mathbb{R} : \frac{\overline{E}}{2}0}$ be a family of solutions of the scaled radiative Navier-Stokes system given in Theorem \ref{thm2}. We call ${\mathcal M}_{\rm ess}^{\varepsilon}\subset (0,T)\times \Omega$ the set \[ {\mathcal M}_{\rm ess}^{\varepsilon}=\left\{(t,x)\in (0,T)\times \Omega : (\varrho_{\varepsilon}(t,x),\vartheta_{\varepsilon}(t,x),E^R_{\varepsilon}(t,x)) \in{\mathcal O}_{\rm ess}\right\}, \] and ${\mathcal M}_{\rm res}^{\varepsilon}=(0,T)\times \Omega \backslash {\mathcal M}_{\rm ess}^{\varepsilon}$ the corresponding residual set. To any measurable function $h$ we associate its decomposition into essential and residual parts \[ h=[h]_{\rm ess}+[h]_{\rm res}, \] where $[h]_{\rm ess}=h\cdot \mathbb{I}_{{\mathcal M}_{\rm ess}^{\varepsilon}}$ and $[h]_{\rm res}=h\cdot \mathbb{I}_{{\mathcal M}_{\rm res}^{\varepsilon}}$. Denoting by $H_{\overline{\vartheta}}$ the Helmholtz function for matter \[ H_{\overline{\vartheta}}(\varrho,\vartheta)=\varrho e-\overline{\vartheta}\varrho s, \] and for radiation \[ {H^R}_{\overline{\vartheta}}(I)= E^R-\overline{\vartheta}s^R, \] and using \eqref{entro} we rewrite \eqref{enertot} as \begin{align*} & \int_{\Omega} \Big(\frac{\varepsilon^2}{2}\ \varrho_{\varepsilon} |\vec u_{\varepsilon}|^2 + H_{\overline{\vartheta}}(\varrho_{\varepsilon},\vartheta_{\varepsilon}) +\varepsilon {H^R}_{\overline{\vartheta}}(I_{\varepsilon}) +\frac{1}{2\zeta}|\vec B_{\varepsilon}|^2- {\frac12 } \varepsilon\varrho_{\varepsilon} \Psi_{\varepsilon} - {\frac12 }\varepsilon^2\varrho_{\varepsilon} |\vec \chi\times\vec x|^2 \Big)\,dx \\ &+\int_0^T\int_0^\infty \int_{\Gamma_+} \vec \omega\cdot \vec n_x I_{\varepsilon}(t,x,\vec \omega,\nu)\,d\Gamma\,d\nu\,dt +\overline{\vartheta}(\varsigma^m_{\varepsilon} +\varsigma^R_{\varepsilon})\big[[0,T]\times \overline{\Omega}\big] \\ &=\int_{\Omega} \Big(\frac{\varepsilon^2}{2}\ \varrho_{0,\varepsilon} |\vec u_{0,\varepsilon}|^2 +\varrho_{0,\varepsilon}e_{0,\varepsilon}+\varepsilon E^R_{0,\varepsilon} +\frac{1}{2\zeta}|\vec B_{0,\varepsilon}|^2 - {\frac12 } \varepsilon\varrho_{0,\varepsilon} \Psi_{0,\varepsilon} \\ &\quad - {\frac12 }\varepsilon^2\varrho_{0,\varepsilon} |\vec \chi\times\vec x|^2 \Big)\,dx. \end{align*} Observing that the total mass is a constant of motion $M=\int_{\Omega}\varrho_{\varepsilon}\,dx=\overline{\varrho}|\Omega|$ and using Hardy-Littlewood-Sobolev inequality, we obtain \[ \frac{\varepsilon}{2}\int_{\Omega}\varrho_{\varepsilon} \Psi_{\varepsilon}\,dx \leq \frac{G\varepsilon}{2}CM^{2/3} \|\varrho_{\varepsilon}\|_{L^{4/3}(\Omega)}^{4/3}. \] By \eqref{m1} and \eqref{m5} we have also $\varrho_{\varepsilon}e(\varrho_{\varepsilon},\vartheta_{\varepsilon}) \geq a\vartheta_{\varepsilon}^4+\frac{3p_{\infty}}{2}\varrho_{\varepsilon}^{5/3}$, so we have the lower bound \[ \int_{\Omega}[H_{\overline{\vartheta}}(\varrho_{\varepsilon},\vartheta_{\varepsilon}) -\frac{1}{2}\varepsilon\varrho_{\varepsilon} \Psi_{\varepsilon}]\,dx \geq c\int_{\Omega}H_{\overline{\vartheta}}(\varrho_{\varepsilon},\vartheta_{\varepsilon})\,dx, \] for $\varepsilon$ small and a $c(\varepsilon)<1$ and we deduce finally the dissipation energy-entropy inequality \begin{equation} \label{energy} \begin{aligned} &\int_{\Omega} \Big( \frac{\varepsilon^2}{2}\ \varrho_{\varepsilon} |\vec u_{\varepsilon}|^2 +H_{\overline{\vartheta}}(\varrho_{\varepsilon},\vartheta_{\varepsilon}) -(\varrho_{\varepsilon}-\overline{\varrho})\partial_{\varrho}H_{\overline{\vartheta}}(\overline{\varrho},\overline{\vartheta}) -H_{\overline{\vartheta}}(\overline{\varrho},\overline{\vartheta}) +\frac{1}{2\zeta}|\vec B_{\varepsilon}|^2 \\ &- \frac{\varepsilon^2}{2}\varrho_{\varepsilon} |\vec \chi\times\vec x|^2 +\varepsilon {H^R}_{\overline{\vartheta}}(I_{\varepsilon}) \Big) dx \\ &\quad +\int_0^T \int_0^\infty\int_{\Gamma_+} I_{\varepsilon} (t,x,\vec \omega,\nu)\ \vec \omega\cdot \vec n_x\, d\Gamma\, d\nu\, dt + \overline{\vartheta}(\varsigma^m_{\varepsilon} + \varsigma^R_{\varepsilon})\big[[0,T]\times \overline{\Omega}\big] \\ &\leq C \int_{\Omega} \Big(\frac{\varepsilon^2}{2}\,\varrho_{0,\varepsilon} |\vec u_{0,\varepsilon}|^2 + H_{\overline{\vartheta}}(\varrho_{0,\varepsilon},\vartheta_{0,\varepsilon}) -(\varrho_{0,\varepsilon}-\overline{\varrho})\partial_{\varrho}H_{\overline{\vartheta}} (\overline{\varrho},\overline{\vartheta}) \\ &\quad -H_{\overline{\vartheta}}(\overline{\varrho},\overline{\vartheta}) +\frac{1}{2\zeta}|\vec B_{0,\varepsilon}|^2 +\varepsilon {H^R}_{\overline{\vartheta}}(I_{0,\varepsilon}) \Big)\,dx. \end{aligned} \end{equation} Now, according to \cite[Lemma 4.1]{DuNe} (see \cite{FEINOV}) we have the following properties for material and radiative Helmholtz functions. \begin{lemma} \label{lem1} Let $\overline{\varrho}>0$ and $\overline{\vartheta}>0$ two given constants and let \[ H_{\overline{\vartheta}}(\varrho,\vartheta)=\varrho e-\overline{\vartheta}\varrho s, \quad {H^R}_{\overline{\vartheta}}(I)= E^R-\overline{\vartheta}s^R. \] Let ${\mathcal O}_{\rm ess}$ and ${\mathcal O}_{\rm res}$ be the sets of essential and residual values introduced in \eqref{essH}--\eqref{res}. There exist positive constants $C_j=C_j(\overline{\varrho},\overline{\vartheta})$ for $j=1,\dots ,4$ and positive constants $C_j=C_j(\overline{E},\overline{\vartheta})$ for $j=5,\dots ,8$ such that \begin{equation} \label{ess1} \begin{aligned} C_1(|\varrho-\overline{\varrho}|^2+|\vartheta-\overline{\vartheta}|^2) &\leq H_{\overline{\vartheta}}(\varrho,\vartheta) -(\varrho-\overline{\varrho})\partial_{\varrho}H_{\overline{\vartheta}}(\overline{\varrho},\overline{\vartheta}) -H_{\overline{\vartheta}}(\overline{\varrho},\overline{\vartheta}) \\ &\leq C_2(|\varrho-\overline{\varrho}|^2+|\vartheta-\overline{\vartheta}|^2), \end{aligned} \end{equation} for all $(\varrho,\vartheta)\in{\mathcal O}_{\rm ess}^H$; \begin{equation} \label{res1} \begin{aligned} &H_{\overline{\vartheta}}(\varrho,\vartheta) -(\varrho-\overline{\varrho})\partial_{\varrho}H_{\overline{\vartheta}}(\overline{\varrho},\overline{\vartheta}) -H_{\overline{\vartheta}}(\overline{\varrho},\overline{\vartheta}) \\ &\geq \inf_{\tilde\varrho,\tilde\vartheta \in{\mathcal O}_{\rm res}^H} \left\{H_{\overline{\vartheta}}(\tilde\varrho,\tilde\vartheta) -(\tilde\varrho-\overline{\varrho})\partial_{\varrho}H_{\overline{\vartheta}}(\overline{\varrho},\overline{\vartheta}) -H_{\overline{\vartheta}}(\overline{\varrho},\overline{\vartheta})\right\} = C_3, \end{aligned} \end{equation} for all $(\varrho,\vartheta)\in{\mathcal O}_{\rm res}^H$; \begin{equation} \label{res2} H_{\overline{\vartheta}}(\varrho,\vartheta) -(\varrho-\overline{\varrho})\partial_{\varrho}H_{\overline{\vartheta}}(\overline{\varrho},\overline{\vartheta}) -H_{\overline{\vartheta}}(\overline{\varrho},\overline{\vartheta}) \geq C_4( \varrho e(\varrho,\vartheta)+\varrho |s(\varrho,\vartheta)|), \end{equation} for all $(\varrho,\vartheta)\in{\mathcal O}_{\rm res}^H$; \begin{equation} \label{ess2} C_5 |E^R-\overline{E}|^2 \leq {H^R}_{\overline{\vartheta}}(I) \leq C_6|E^R-\overline{E}|^2, \end{equation} for all $E\in{\mathcal O}_{\rm ess}^R$; \begin{equation} \label{res3} {H^R}_{\overline{\vartheta}}(I) \geq \inf_{\tilde I \in{\mathcal O}_{\rm res}^R} {H^R}_{\overline{\vartheta}}(\tilde I) = C_7, \end{equation} for all $E\in{\mathcal O}_{\rm res}^R$; \begin{equation} \label{res4} {H^R}_{\overline{\vartheta}}(I) \geq C_8( E^R(I)+ |s^R(I)|) \end{equation} for all $E\in{\mathcal O}_{\rm res}^R$. \end{lemma} Using \eqref{energy} and Lemma \ref{lem1}, we obtain the following energy estimates \begin{lemma} \label{lem2} Suppose that the initial data satisfy \begin{gather*} \|[ \varrho_{0,\varepsilon}-\overline{\varrho}]_{\rm ess}\|^2_{L^2(\Omega)} \leq C\varepsilon^2, \quad \|[ \vartheta_{0,\varepsilon}-\overline{\vartheta}]_{\rm ess}\|^2_{L^2(\Omega)} \leq C\varepsilon^2, \quad \| E^R_{0,\varepsilon}-\overline{E}\|^2_{L^2(\Omega)}\leq C{\varepsilon}, \\ \| \vec B_{0,\varepsilon}\|_{L^2(\Omega;{\mathbb{R}^3)}}^2\leq C\varepsilon^2, \quad \| \sqrt{\varrho_{0,\varepsilon}} \vec u_{0,\varepsilon}\|_{L^2(\Omega;\mathbb{R}^3)}\leq C. \end{gather*} Then the following estimates hold \begin{gather} \label{Mres} \operatorname{ess\,sup}_{t \in (0,T)} | {\mathcal M}^{\varepsilon}_{\rm res}(t)| \leq C\varepsilon^2, \\ \label{ress} \operatorname{ess\,sup}_{t \in (0,T)} \|[ \varrho_{\varepsilon} -\overline{\varrho}]_{\rm ess}(t)\|^2_{L^2(\Omega)}\leq C\varepsilon^2, \\ \label{tess} \operatorname{ess\,sup}_{t \in (0,T)} \|[ \vartheta_{\varepsilon} -\overline{\vartheta}]_{\rm ess}(t)\|^2_{L^2(\Omega)}\leq C\varepsilon^2, \\ \label{ERess} \operatorname{ess\,sup}_{t \in (0,T)} \|[ E^R_{\varepsilon} -\overline{E}]_{\rm ess}(t)\|^2_{L^2(\Omega)}\leq C\varepsilon, \\ \label{reres} \operatorname{ess\,sup}_{t \in (0,T)} \|[ \varrho_{\varepsilon}e(\varrho_{\varepsilon}, \vartheta_{\varepsilon})]_{\rm res}(t)\|_{L^1(\Omega)}\leq C\varepsilon^2, \\ \label{rsres} \operatorname{ess\,sup}_{t \in (0,T)} \|[ \varrho_{\varepsilon}s(\varrho_{\varepsilon}, \vartheta_{\varepsilon})]_{\rm res}(t)\|_{L^1(\Omega)}\leq C\varepsilon^2, \\ \label{ERres} \operatorname{ess\,sup}_{t \in (0,T)} \|[ E^R(I_{\varepsilon})]_{\rm res}(t) \|_{L^1(\Omega)}\leq C\varepsilon, \\ \label{sRres} \operatorname{ess\,sup}_{t \in (0,T)} \|[ s^R(I_{\varepsilon})]_{\rm res}(t) \|_{L^1(\Omega)}\leq C\varepsilon, \\ \label{sigres} (\varsigma^m_{\varepsilon}+\varsigma^R_{\varepsilon}) [[0,T]\times \overline{\Omega}] \leq C\varepsilon^2, \\ \label{rB} \operatorname{ess\,sup}_{t \in (0,T)} \| \frac{\vec B_{\varepsilon}(t)}{\varepsilon}\|_{L^2(\Omega;\mathbb{R}^3)}\leq C, \\ \label{ru} \operatorname{ess\,sup}_{t \in (0,T)} \| \sqrt{\varrho_{\varepsilon}} \vec u_{\varepsilon}(t)\|_{L^2(\Omega;\mathbb{R}^3)}\leq C, \\ \label{rtpower} \operatorname{ess\,sup}_{t \in (0,T)} \int_{\Omega}( [\varrho_{\varepsilon}]_{\rm res}^{\frac{5}{3}} +[ \vartheta_{\varepsilon}]_{\rm res}^4)(t)\,dx\leq C\varepsilon^2, \\ \label{uW} \|\vec u _{\varepsilon}\|_{L^2(0,T;W^{1,2}(\Omega;\mathbb{R}^3))}\leq C, \\ \label{tW} \| \frac{\vartheta_{\varepsilon}-\overline{\vartheta}}{\varepsilon^2} \|_{L^2(0,T;W^{1,2}(\Omega))}\leq C, \\ \label{logtW} \| \frac{\log(\vartheta_{\varepsilon})-\log(\overline{\vartheta})}{\varepsilon^2} \|_{L^2(0,T;W^{1,2}(\Omega))}\leq C, \\ \label{BW} \|\frac{\vec B _{\varepsilon}}{\varepsilon}\|_{L^2(0,T;W^{1,2}(\Omega;\mathbb{R}^3))} \leq C. \end{gather} \end{lemma} \begin{proof} Estimate \eqref{Mres} follows from \eqref{res1}. Bounds \eqref{ress}, \eqref{tess} and \eqref{ERess} follow from \eqref{ess1} and \eqref{ess2}. Estimates \eqref{reres} and \eqref{rsres} follow from \eqref{res2}. Bounds \eqref{ERres} and \eqref{sRres} follow from \eqref{res4}. Estimates \eqref{sigres}, \eqref{rB} and \eqref{ru} follow from the dissipation energy-entropy inequality \eqref{energy}. Bound \eqref{rtpower} follows from \eqref{reres} and \eqref{m5} (cf. a lower bound for $\varrho e$ before \eqref{energy}). From \eqref{sigres} we see that \begin{equation} \label{Du} \|\nabla_x \vec{u}_{\varepsilon}+\nabla_x^T \vec{u}_{\varepsilon}-\frac{2}{3} \operatorname{div}_x \vec{u}_{\varepsilon}\mathbb{I} \|_{L^2(0,T; L^2(\Omega;\mathbb{R}^{3\times 3}))}\leq C. \end{equation} From \eqref{m7}, \eqref{ru} and \eqref{Du} we obtain \eqref{uW}. Details can be found in \cite{dkn,FEINOV}. From \eqref{sigres} we obtain \[ \| \nabla_x(\frac{\vartheta_{\varepsilon}}{\varepsilon^2})\|_{L^2(0,T; L^2(\Omega;\mathbb{R}^{3}))} +\| \nabla_x(\frac{\log\vartheta_{\varepsilon}}{\varepsilon^2}) \|_{L^2(0,T; L^2(\Omega;\mathbb{R}^{3}))}\leq C, \] which, using Poincar\'e inequality, gives \eqref{tW} and \eqref{logtW}. Finally by \eqref{m8z}, \eqref{iprod} and \eqref{sigres} one gets \[ \|\frac{\operatorname{curl}_x \vec B_{\varepsilon}}{\varepsilon}\|_{L^2(0,T;L^2(\Omega;\mathbb{R}^3))} \leq C, \] and \eqref{BW} follows by using \cite[Theorem 6.1]{DULI}. \end{proof} Our goal in the next Section will be to prove that the incompressible system \eqref{bous1}-\eqref{bous9} is the limit of the primitive system \eqref{mass}-\eqref{transfer} in the following sense \begin{theorem} \label{conv} Let $\Omega\subset \mathbb{R}^3$ be a bounded domain of class $C^{2,\nu}$. Assume that the thermodynamic functions $p$, $e$, $s$ satisfy hypotheses \eqref{m1}--\eqref{m6} with $P\in C^1([0,\infty))\cap C^2(0,\infty)$, and that the transport coefficients $\mu$, $\eta$, $\kappa$, $\lambda$, $\sigma_a$, $\sigma_s$ and the equilibrium function $B$ comply with \eqref{m7}--\eqref{m10}. Let $(\varrho_{\varepsilon},\vec u_{\varepsilon},\vartheta_{\varepsilon}, \vec B_{\varepsilon}, I_{\varepsilon})$ be a weak solution of the scaled system \eqref{i1}--\eqref{i9} for $(t,x,\vec \omega,\nu)\in [0,T]\times \Omega \times \mathcal{S}^2 \times \mathbb{R}_+$, supplemented with the boundary conditions \eqref{i11}--\eqref{i12} and initial conditions $(\varrho_{0,\varepsilon},\vec u_{0,\varepsilon},\vartheta_{0,\varepsilon}, \vec B_{0,\varepsilon},I_{0,\varepsilon})$ \begin{gather*} \varrho_{\varepsilon}(0,\cdot)=\overline{\varrho}+\varepsilon \varrho_{0,\varepsilon}^{(1)},\quad \vec u_{\varepsilon}(0,\cdot)= \vec u_{0,\varepsilon}, \quad \vartheta_{\varepsilon}(0,\cdot)=\overline{\vartheta}+\varepsilon^3 \vartheta_{0,\varepsilon}^{(3)},\\ I_{\varepsilon}(0,\cdot)=\overline{I}+\varepsilon I_{0,\varepsilon}^{(1)}, \quad \vec B_{\varepsilon}(0,\cdot)= \varepsilon \vec B_{0,\varepsilon}^{(1)}, \end{gather*} where $\overline{\varrho}>0,\ \overline{\vartheta}>0,\ \overline{I}>0$ are constants in $(0,T)\times\Omega$ and \[ \int_{\Omega}\varrho_{0,\varepsilon}^{(1)}\,dx=0,\quad \int_{\Omega}\vartheta_{0,\varepsilon}^{(3)}\,dx=0, \quad \int_{\Omega}I_{0,\varepsilon}^{(1)}\,dx=0,\quad \int_{\Omega}\vec B_{0,\varepsilon}^{(1)}\,dx=0\quad \text{for all } \varepsilon>0. \] Assume that \begin{gather*} \varrho_{0,\varepsilon}^{(1)}\to \varrho_0^{(1)}\quad\text{weakly-(*) in } L^{\infty}(\Omega),\\ \vec u_{0,\varepsilon}\to \vec U_0\quad\text{weakly-(*) in } L^{\infty}(\Omega;\mathbb{R}^3),\\ \vartheta_{0,\varepsilon}^{(3)}\to \vartheta_0^{(3)}\quad\text{weakly-(*) in } L^{\infty}(\Omega),\\ I_{0,\varepsilon}^{(1)}\to I_0^{(1)}\quad\text{weakly-(*) in } L^{\infty}(\Omega\times\mathcal{S}^2\times \mathbb{R}_+),\\ \vec B_{0,\varepsilon}^{(1)}\to \vec B_0^{(1)}\quad\text{weakly-(*) in } L^{\infty}(\Omega;\mathbb{R}^3), \end{gather*} Then \begin{equation} \label{conv1} \operatorname{ess\,sup}_{t \in (0,T)} \| \varrho_{\varepsilon}(t)-\overline{\varrho}\|_{L^{\frac{5}{3}}(\Omega)} \leq C\varepsilon, \end{equation} and up to subsequences \begin{gather} \label{conv2} \vec u_{\varepsilon}\to \vec U\quad\text{weakly in } L^2(0,T;W^{1,2}(\Omega; \mathbb{R}^3)),\\ \label{conv3} \frac{\vartheta_{\varepsilon}-\overline{\vartheta}}{\varepsilon^3}=:\vartheta^{(3)}_\varepsilon\to \Theta \quad\text{weakly in } L^\frac{4}{3}(0,T;W^{1,\frac{4}{3}}(\Omega))\\ \label{conv4} I_{\varepsilon}\to \overline{I} = B_0\quad\text{weakly in } L^2(0,T;L^2(\Omega\times\mathcal{S}^2\times \mathbb{R}_+)), \\ \label{conv5} \frac{\vec B_{\varepsilon}}{\varepsilon}=\vec B^{(1)}_\varepsilon\to \vec B \quad\text{weakly in } L^2(0,T;W^{1,2}(\Omega; \mathbb{R}^3)), \\ \label{conv6} \frac{I_{\varepsilon}-\overline{I}}{\varepsilon}=I^{(1)}_\varepsilon\to I_1 \quad\text{weakly in } L^2(0,T;L^2(\Omega\times\mathcal{S}^2\times \mathbb{R}_+)), \end{gather} where $(\vec U,\Theta,\vec B, I_1)$ solves the system \eqref{bous1}-\eqref{bous5}. \end{theorem} \section{Proof of Theorem \ref{conv}} \label{pconv} Let us first quote the following result of \cite{DuNe} (see \cite{FEINOV}). \begin{proposition} \label{prop2} Let $\{\varrho_{\varepsilon}\}_{\varepsilon>0}, \{\vartheta_{\varepsilon}\}_{\varepsilon>0}, \{I_{\varepsilon}\}_{\varepsilon>0}$ be three sequences of non-negative measurable functions such that \begin{gather*} [\varrho_{\varepsilon}^{(1)}]_{\rm ess}\to \varrho^{(1)}\quad\text{weakly-(*) in } L^{\infty}(0,T;L^2(\Omega)), \\ [\vartheta_{\varepsilon}^{(1)}]_{\rm ess}\to \vartheta^{(1)}\quad\text{weakly-(*) in } L^{\infty}(0,T;L^2(\Omega)), \\ [I_{\varepsilon}^{(1)}]_{\rm ess}\to I^{(1)}\quad\text{weakly-(*) in } L^{\infty}(0,T;L^2(\Omega)),\text{ a.e. in } \mathcal{S}^2\times \mathbb{R}_+, \end{gather*} where \[ \varrho_{\varepsilon}^{(1)}=\frac{\varrho_{\varepsilon}-\overline{\varrho}}{\varepsilon},\quad \vartheta_{\varepsilon}^{(1)}=\frac{\vartheta_{\varepsilon}-\overline{\vartheta}}{\varepsilon},\quad I_{\varepsilon}^{(1)}=\frac{I_{\varepsilon}-\overline{I}}{\varepsilon}. \] Suppose that \begin{equation} \label{Mresbis} \operatorname{ess\,sup}_{t \in (0,T)} | {\mathcal M}^{\varepsilon}_{\rm res}(t)| \leq C\varepsilon^2. \end{equation} Let $G,G^R \in C^1(\overline{{\mathcal O}_{\rm ess}})$ be given functions. Then \[ \frac{[G(\varrho_{\varepsilon},\vartheta_{\varepsilon})]_{\rm ess} -G(\overline{\varrho},\overline{\vartheta})}{\varepsilon} \to \frac{\partial G(\overline{\varrho},\overline{\vartheta})}{\partial \varrho}\ \varrho^{(1)} +\frac{\partial G(\overline{\varrho},\overline{\vartheta})}{\partial \vartheta}\ \vartheta^{(1)}, \] weakly-(*) in $L^{\infty}(0,T;L^2(\Omega))$, and if we denote \[ [G^R(I_{\varepsilon})]_{\rm ess}:=[G^R(I_{\varepsilon}(\cdot,\cdot,\vec \omega, \nu))]_{\rm ess}=G^R(I_{\varepsilon})\cdot \mathbb{I}_{{\mathcal M}^{\varepsilon}_{\rm ess}}, \quad \text{for a.a. } (\vec \omega,\nu)\in \mathcal{S}^2\times \mathbb{R}_+, \] we have \[ \frac{[G^R(I_{\varepsilon})]_{\rm ess}-G^R(\overline{I})}{\varepsilon} \to \frac{\partial G(\overline{I})}{\partial I}\ I^{(1)}, \] weakly-(*) in $L^{\infty}(0,T;L^2(\Omega))$, a.e. in $ \mathcal{S}^2\times \mathbb{R}_+$. Moreover if $G,G^R \in C^2(\overline{{\mathcal O}_{\rm ess}})$ then \[ \big\| \frac{[G(\varrho_{\varepsilon},\vartheta_{\varepsilon})]_{\rm ess} -G(\overline{\varrho},\overline{\vartheta})}{\varepsilon} -\frac{\partial G(\overline{\varrho},\overline{\vartheta})}{\partial \varrho}\ [\varrho^{(1)}]_{\rm ess} -\frac{\partial G(\overline{\varrho},\overline{\vartheta})}{\partial \vartheta}\ [\vartheta^{(1)}]_{\rm ess} \big\|_{L^{\infty}(0,T;L^1(\Omega))}\leq C\varepsilon, \] and \[ \big\|\frac{[G^R(I_{\varepsilon})]_{\rm ess}-G^R(\overline{I})}{\varepsilon} -\frac{\partial G(\overline{I})}{\partial I}\ [I^{(1)}]_{\rm ess} \big\|_{L^{\infty}(0,T;L^1(\Omega))}\leq C\varepsilon, \] for a.a. $(\vec \omega,\nu)\in \mathcal{S}^2\times \mathbb{R}_+$. \end{proposition} Clearly, this result provides us with the convergence properties \eqref{conv1}--\eqref{conv2}, \eqref{conv5}--\eqref{conv6}. The convergence of radiative intensity \eqref{conv4} follows from \eqref{ERres}, \eqref{ERess}, and the linearity of \eqref{radeq}, cf.~the section \ref{subsec:rte} Radiative transfer equation. The equilibrium Planck function $B_0$ does not satisfy the boundary condition \eqref{i12}${}_1$ however since it is isotropic; therefore has to be modified at the boundary $\partial\Omega$. The last convergence \eqref{conv3} is postponed to Section \ref{subsec:eb}. To conclude the proof of Theorem \ref{conv}, let us prove that the limit quantities $(\vec U,\Theta,\vec B,I_1)$ solve the target system \eqref{bous1}-\eqref{bous5}. As number of terms in the equations of our model are similar to those of the radiative Navier-Stokes-Fourier analyzed in \cite{DuNe} we focus on the new contributions only. \subsection{Continuity and momentum equations} For the continuity equation, one expects that in the low Mach number limit, it reduces to the incompressibility constraint. In fact, from Lemma \ref{lem2} we know that $\int_0^T \|\vec u _{\varepsilon}(t)\|^2_{W^{1,2}(\Omega;\mathbb{R}^3)}\,dt\leq C$ so passing to the limit after possible extraction of a subsequence, we deduce that \begin{equation}\label{convu} \vec u _{\varepsilon} \to \vec U , \quad\text{weakly in } L^2(0,T;W^{1,2}(\Omega;\mathbb{R}^3)). \end{equation} In the same stroke $\varrho_{\varepsilon} \to \overline{\varrho}$, weakly in $L^{\infty}(0,T;L^{5/3}(\Omega;\mathbb{R}^3))$. So we can pass to the limit in the weak continuity equation \eqref{mass} which gives $\int_0^T\int_{\Omega} \vec U\cdot \nabla_x \phi\,dx\,dt=0$ for all $\phi\in {\mathcal D}((0,T)\times \overline{\Omega})$, which rewrites \[ \operatorname{div}_x \vec U=0,\quad \text{a.e. in } (0,T)\times \Omega,\quad \vec U\big|_{\partial \Omega}=0, \] provided $\partial\Omega$ is regular. For the momentum equation one knows that due to possible strong time oscillations of the gradient component of velocity, one has only $\varrho_{\varepsilon} \vec u _{\varepsilon}\otimes \vec u _{\varepsilon}\to\overline{\varrho \vec U\otimes \vec U}$ weakly in $L^2(0,T;L^{\frac{30}{29}}(\Omega;\mathbb{R}^3))$. However one can show by the analysis in \cite{FEINOV} that one can pass to the limit in the convective term and obtain \[ \int_0^T \int_{\Omega} \overline{ \varrho\ \vec U\otimes \vec U}:\nabla_x \vec{\phi}\,dx\,dt \to \int_0^T \int_{\Omega} \overline{\varrho}\ \vec U\otimes \vec U:\nabla_x \vec{\phi}\,dx\,dt. \] According to the hypotheses on the pressure law, the temperature $\vartheta_{\varepsilon}$ is bounded in $L^{\infty}((0,T);L^4(\Omega))\cap L^2(0,T;L^6(\Omega))$, which together with the strong convergence of $\vartheta_\varepsilon$ by \eqref{tW} imply that $\mathbb{S}_{\varepsilon}\to\mu(\overline{\vartheta})(\nabla_x \vec U+\nabla_x^T \vec U)$ weakly in $L^{\frac{34}{23}}(0,T;L^{\frac{34}{23}}(\Omega;\mathbb{R}^3))$. So taking a divergence free test vector field $\vec{\phi}$ in \eqref{mom}, we have \begin{equation} \label{mombis} \begin{aligned} &\int_0^T\int_{\Omega}\big( \varrho_{\varepsilon} \vec u_{\varepsilon}\cdot \partial_t\vec{\phi}+\varrho_{\varepsilon} \vec u_{\varepsilon}\otimes \vec u_{\varepsilon}:\nabla_x \vec{\phi} -2\varrho_{\varepsilon} \vec \chi \times \vec{u}_{\varepsilon} \cdot \vec{\phi} \big)\,dx\,dt \\ &=\int_0^T\int_{\Omega} \Big( \mathbb{S}_{\varepsilon}:\nabla_x \vec{\phi} -\frac{\varrho_{\varepsilon}-\overline{\varrho}}{\varepsilon} \nabla_x \Psi_{\varepsilon} \cdot {\vec{\phi}} -\frac{1}{\zeta}\frac{{\operatorname{curl}_x \vec{B}_{\varepsilon}}}{\varepsilon} \times \frac{\vec{B}_{\varepsilon}}{\varepsilon} \cdot {\vec{\phi}} \\ &\quad -\frac12\varepsilon^2\varrho_{\varepsilon} \nabla_x |\vec \chi \times \vec x|^2 \cdot {\vec{\phi}}\Big)\,dx\,dt -\int_{\Omega} \varrho_{0,\varepsilon} \vec u_{0,\varepsilon} \cdot \vec{\phi}(0,\cdot)\,dx. \end{aligned} \end{equation} Moreover, using \eqref{v9} together with estimates \eqref{rB}, \eqref{BW} and Aubin-Lions lemma we obtain \begin{equation} \label{convB} \begin{gathered} \frac{\vec{B}_{\varepsilon}}{\varepsilon} \to \vec{B} \quad \text{weakly in $L^2(0,T;W^{1,2}(\Omega;\mathbb{R} ^3))$ and strongly in } L^2(0,T; L^q(\Omega, \mathbb{R}^3)), \\ \frac{1}{\zeta}\frac{{\operatorname{curl}_x \vec{B}_{\varepsilon}}}{\varepsilon} \times \frac{\vec{B}_{\varepsilon}}{\varepsilon} \to \frac{1}{\zeta}{{\operatorname{curl}_x \vec{B}}} \times \vec{B}\quad \text{ weakly in } L^{5/4}((0,T)\times \Omega;\mathbb{R} ^3), \end{gathered} \end{equation} for any $1\le q<6$. Then passing to the limit and using \eqref{conv2}-\eqref{conv6}, we obtain \begin{align*} &\int_0^T \int_{\Omega} \big( \overline{\varrho}\vec U\cdot \partial_t\vec{\phi}+\overline{\varrho} \vec U\otimes \vec U:\nabla_x \vec{\phi} -2\overline{\varrho} \vec \chi\times \vec U\cdot \vec{\phi} \big)\,dx\,dt \\ &=\int_0^T \int_{\Omega} \Big( \mu(\overline{\vartheta})(\nabla_x \vec U+\nabla_x^T \vec U):\nabla_x \vec{\phi}-\varrho_1\nabla_x \Psi(\overline{\varrho})\cdot \vec{\phi} -\frac{1}{\zeta}\operatorname{curl}_x \vec B\times\vec B\cdot \vec{\phi} \Big)\,dx\,dt \\ &\quad -\int_{\Omega} \overline{\varrho}\vec U_0\cdot \vec{\phi}\,dx, \end{align*} provided that $ \vec u _{0,\varepsilon}\to \vec U_0$ weakly-$*$ in $L^{\infty}(\Omega;\mathbb{R}^3)$. Moreover as in \cite{FEINOV}, the formal relation between $\varrho^{(1)}$ and $\overline{\varrho}$ is recovered by multiplying the momentum equation by $\varepsilon$. One gets, using Proposition \ref{prop2} and passing to the limit $\varepsilon\to0$, \begin{equation} \label{Gp} \int_0^T \int_{\Omega} \big(\nabla_x p^{(1)} -\overline{\varrho}\nabla_x \Psi(\overline{\varrho}) \big)\cdot \phi\,dx\,dt=0, \end{equation} which is the weak formulation of \begin{equation} \label{press1} \partial_{\varrho}p(\overline{\varrho},\overline{\vartheta})\nabla_x \varrho^{(1)} +\partial_{\vartheta}p(\overline{\varrho},\overline{\vartheta})\nabla_x \vartheta^{(1)} -\overline{\varrho}\nabla_x \Psi(\overline{\varrho})=0. \end{equation} We rewrite this as \begin{equation} \label{press2} \partial_{\varrho}p(\overline{\varrho},\overline{\vartheta})\nabla_x \varrho_1 -\overline{\varrho}\nabla_x \Psi(\overline{\varrho})=0, \end{equation} once we establish that $\vartheta^{(1)}=\vartheta_1=\vartheta_2=0$ in the section \ref{subsec:eb}. That means we have got an explicit formula for $\varrho_1$ \begin{equation} \label{rhoone} \varrho_1 = \frac{{\overline{\varrho}}\Psi(\overline{\varrho})}{\partial_{\varrho}p(\overline{\varrho}, \overline{\vartheta})} + h(t), \end{equation} where $h$ is an undetermined function. \subsection{Radiative transfer equation} \label{subsec:rte} Using the $L^{\infty}$ bound shown in (\cite{dkn}) for $I_\varepsilon$, based on the initial data bound \eqref{data1}, it is clear that $I_{\varepsilon} \to I_0$ weakly in $L^2((0,T)\times\Omega\times\mathcal{S}^2\times\mathbb{R}_+)$, and we have also by \eqref{tW} $\vartheta_{\varepsilon} \to \overline{\vartheta}$ weakly in $L^2(0,T;W^{1,2}(\Omega))$. By using the cut-off hypotheses \eqref{m9}, \eqref{m10} and the same notation for any time-independent test function $\psi\in C^{\infty}_c(\overline{\Omega}\times \mathcal{S}^2\times \mathbb{R}_+)$, we can pass to the limit in \eqref{transfer} and obtain \begin{align*} &\int_{\Omega}\int_0^{\infty}\int_{\mathcal{S}^2} \vec \omega \cdot \nabla_x \psi\ I_0\,d\vec \omega\,d\nu\,dx \\ &+ \int_{\Omega}\int_0^{\infty}\int_{\mathcal{S}^2} [\sigma_a(\nu,\overline{\vartheta})(B(\nu,\overline{\vartheta})-I_0) +\sigma_s(\nu,\overline{\vartheta})(\tilde I_0-I_0)] \psi\,d\vec \omega\,d\nu\,dx \\ &= \int_{\Gamma_+} \int_0^{\infty}I_0\,\vec \omega\cdot \vec n_x \ \psi\,d\Gamma\,d\nu, \end{align*} which is the weak formulation of the stationary problem \begin{equation} \label{transfer0} \vec \omega \cdot \nabla_x I_0=S_0, \end{equation} with the boundary condition \begin{equation} \label{transfer0BC} I_0=0\quad \text{on } \Gamma_{-}, \end{equation} where $S_0=\sigma_a(\nu,\overline{\vartheta})(B(\nu,\overline{\vartheta})-I_0) +\sigma_s(\nu,\overline{\vartheta})(\tilde I_0-I_0)$. The solution of \eqref{transfer0}--\eqref{transfer0BC} is the function equal to $B(\nu,\overline{\vartheta}) = B_0$ in $\Omega$ and $0$ on $\Gamma_{-}$. This solution is unique at least for a particular type of domains thanks to the linearity of \eqref{transfer0}. Subtracting from \eqref{transfer} and dividing by $\varepsilon$ gives \begin{align*} &\int_0^T\int_{\Omega}\int_0^{\infty}\int_{\mathcal{S}^2} (\varepsilon\partial_t\psi +\vec \omega \cdot \nabla_x \psi) \frac{I_{\varepsilon}-I_0}{\varepsilon} \,d\vec \omega\,d\nu\,dx\,dt \\ &+\int_0^T \int_{\Omega}\int_0^{\infty}\int_{\mathcal{S}^2} [ \frac{S_{\varepsilon}-S_0}{\varepsilon}] \psi\,d\vec \omega\,d\nu\,dx\,dt, \\ &=-\int_{\Omega}\int_0^{\infty}\int_{\mathcal{S}^2} \varepsilon \frac{I_{0,\varepsilon}-I_0}{\varepsilon} \psi(0,x,\vec \omega,\nu)\,d\vec \omega\,d\nu\,dx \\ &\quad +\int_0^T \int_{\Gamma_+} \int_0^{\infty}\vec \omega\cdot \vec n_x \frac{I_{\varepsilon}-I_0}{\varepsilon} \psi\,d\Gamma\,d\nu\,dt, \end{align*} for any $\psi\in C^{\infty}_c([0,T)\times \overline{\Omega}\times \mathcal{S}^2\times \mathbb{R}_+)$, with $S_{\varepsilon}-S_0 := S(I_{\varepsilon})-S(I_0)$. From Proposition \ref{prop2}, we obtain \begin{align*} \frac{S_{\varepsilon}-S_0}{\varepsilon} \to S_1:=&\partial_{\vartheta}(\sigma_a B)(\nu,\overline{\vartheta})\vartheta^{(1)} -\partial_{\vartheta}\sigma_a(\nu,\overline{\vartheta})\vartheta^{(1)}I_0 -\sigma_a(\nu,\overline{\vartheta})I_1 \\ & +\partial_{\vartheta}\sigma_s(\nu,\overline{\vartheta})\vartheta^{(1)}\tilde I_0 +\sigma_s(\nu,\overline{\vartheta})\tilde I_1 -\partial_{\vartheta}\sigma_s(\nu,\overline{\vartheta})\vartheta^{(1)} I_0 -\sigma_s(\nu,\overline{\vartheta}) I_1 \\ &= -\sigma_a(\nu,\overline{\vartheta})I_1 + \sigma_s(\nu,\overline{\vartheta}) (\tilde I_1 - I_1) \end{align*} weakly-* in $L^{\infty}(0,T;L^2(\Omega\times\mathcal{S}^2\times\mathbb{R}_+))$ with $I_1:=I^{(1)}$. Passing to the limit we find the limit equation based on the assumption $I_{0,\varepsilon}^{(1)}\to I_{0}^{(1)}$ weakly-* in $L^\infty({\Omega}\times \mathcal{S}^2\times \mathbb{R}_+)$ \begin{equation} \label{trans2} \begin{aligned} &\int_{\Omega}\int_0^{\infty}\int_{\mathcal{S}^2} \vec \omega \cdot \nabla_x \psi\ I_1\,d\vec \omega\,d\nu\,dx + \int_{\Omega}\int_0^{\infty}\int_{\mathcal{S}^2} S_1 \psi\,d\vec \omega\,d\nu\,dx \\ &= \int_{\Gamma_+} \int_0^{\infty}I_1\,\vec \omega\cdot \vec n_x \psi\,d\Gamma\,d\nu, \end{aligned} \end{equation} using the same notation for any time-independent test function $\psi\in C^{\infty}_c(\overline{\Omega}\times \mathcal{S}^2\times \mathbb{R}_+)$ which is the weak formulation of the stationary problem \begin{equation} \label{transfer1} \vec \omega \cdot \nabla_x I_1=S_1, \end{equation} with the boundary condition \begin{equation} \label{transfer1BC} I_1=0\quad \text{on } \Gamma_{-}. \end{equation} \subsection{Entropy balance} \label{subsec:eb} First of all we analyze the weak limit of \eqref{entro}, then we subtract it (with weak limits denoted by bars) from \eqref{entro} and divide by $\varepsilon$ as in the last section. We follow the ideas of \cite{FN} and \cite{K2}. The most obvious convergence in \eqref{entro} is in the entropy production rate measures. By virtue of \eqref{sigres} it holds \begin{gather} \label{mechradentrprodconv0} \langle \varsigma^m_{\varepsilon}+\varsigma^R_{\varepsilon} ; \varphi\rangle _{[{\mathcal M};C]([0,T)\times \overline{\Omega})}\to0\quad \text{as } \varepsilon\to 0^+, \\ \label{mechradentrprodconv1} \frac{1}{\varepsilon}\langle \varsigma^m_{\varepsilon}+\varsigma^R_{\varepsilon} ; \varphi\rangle _{[{\mathcal M};C]([0,T)\times \overline{\Omega})}\to 0\quad \text{as } \varepsilon\to 0^+. \end{gather} We split the heat flux term into residual and essential parts as follows: \begin{equation} \label{heatfluxconv0} \begin{aligned} -\int_0^T \int_{\Omega}\frac{ {\vec q}_{\varepsilon} }{\varepsilon^2\vartheta_{\varepsilon}} \cdot \nabla_x \varphi\,dx\,dt & =\int_0^T \int_{\Omega}\frac{ { \kappa([\vartheta_\varepsilon]_{\rm res})} }{[\vartheta_{\varepsilon}]_{\rm res}} \frac{\nabla_x\vartheta_\varepsilon}{\varepsilon^2}\cdot \nabla_x \varphi\,dx\,dt\\ &\quad + \int_0^T \int_{\Omega}\frac{ { \kappa([\vartheta_\varepsilon]_{\rm ess})} } {[\vartheta_{\varepsilon}]_{\rm ess}}\frac{\nabla_x\vartheta_\varepsilon}{\varepsilon^2}\cdot \nabla_x \varphi\,dx\,dt. \end{aligned} \end{equation} The first term on the right-hand side vanishes. The argument is as follows: Firstly, from \eqref{sigres}, \eqref{iprod} and \eqref{m8} we obtain an exact estimate \begin{equation} \label{tempsobolev} \int_0^T \int_{\Omega}\vartheta_\varepsilon|\nabla_x\frac{\vartheta_\varepsilon-\overline{\vartheta}}{\varepsilon^2}|^2\,dx\,dt\le c. \end{equation} From \eqref{tW} we know that $\|\frac{\vartheta_\varepsilon-\overline{\vartheta}}{\varepsilon}\|_{L^2(0,T;W^{1,2}(\Omega))}\le c$, thus \begin{equation} \label{tempstrong} \vartheta_\varepsilon\to\overline{\vartheta}\quad\text{in } L^2(0,T;W^{1,2}(\Omega)) \end{equation} strongly. On the residual set we now apply the Sobolev embedding and interpolate \eqref{tempstrong} with the information in \eqref{rtpower}. (Similarly we obtain $ \vartheta_\varepsilon^{(1)}\to 0$ in $L^2(0,T;W^{1,2}(\Omega))$ as well.) This leads to the convergence \begin{equation} \label{tempresstronghom} [\vartheta_\varepsilon]_{\rm res}\to[\overline{\vartheta}]_{\rm res} = 0\quad\text{in } L^{\frac{14}{3}}(0,T;L^{\frac{14}{3}}(\Omega)), \end{equation} meaning that the first integral in \eqref{heatfluxconv0} converges. With the intention that its limit is zero we apply \eqref{m8} and split the integral into two parts. The second part, namely, \begin{equation} \label{resheatconv2} \int_0^T \int_{\Omega}([\vartheta_\varepsilon]_{\rm res}^{\frac{3}{2}} -\overline{\vartheta}^{\frac{3}{2}}+\overline{\vartheta}^{\frac{3}{2}})\sqrt{[\vartheta_\varepsilon]_{\rm res}} \nabla_x\frac{ \vartheta_\varepsilon }{\varepsilon^2}\cdot \nabla_x \varphi\,dx\,dt \end{equation} converges to zero as $\varepsilon\to0$ with the rate $\varepsilon^2$ by the Poincar\'e inequality \begin{equation} \label{temppoincare} \|\vartheta_\varepsilon^{\frac{3}{2}}-\overline{\vartheta}^{\frac{3}{2}}\|_{L^2(0,T;L^{2}(\Omega))}\le c \|\sqrt{\vartheta_\varepsilon}\nabla_x\vartheta_\varepsilon\|_{L^2(0,T;L^{2}(\Omega))} \le C\varepsilon^2 \end{equation} by \eqref{tempsobolev}, \eqref{Mres} and \eqref{rtpower}. The first part, namely, \begin{equation} \label{resheatconv1} \int_0^T \int_{\Omega} [1]_{\rm res}\vartheta_\varepsilon^{-1}\nabla_x\frac{ \vartheta_\varepsilon }{\varepsilon^2}\cdot \nabla_x \varphi\,dx\,dt \end{equation} converges to zero as $\varepsilon\to0$ with the rate $\varepsilon$ by Cauchy-Schwarz inequality, \eqref{logtW} and \eqref{Mres}. The second term on the right-hand side of \eqref{heatfluxconv0} converges by virtue of \eqref{tempsobolev} and \eqref{tempstrong}, at least for a subsequence, to \begin{equation} \label{heatres} \int_0^T \int_{\Omega} \kappa(\overline{\vartheta})\overline{\vartheta}^{-1}\nabla_x\vartheta_2\cdot \nabla_x \varphi\,dx\,dt. \end{equation} For the convergence of the initial entropies in \eqref{entro} we use Proposition \ref{prop2} and we obtain \begin{equation} \label{initent1} \begin{aligned} &- \int_{\Omega} \big\{\big( \frac{[(\varrho s)_{0,\varepsilon}]_{\rm ess}-\overline{\varrho} s(\overline{\varrho},\overline{\vartheta})}{\varepsilon} +\varepsilon\frac{[s^R_{0,\varepsilon}]_{\rm ess}-{s^R(\overline{I})}}{\varepsilon} \big) \varphi (0, \cdot) \big\}\,dx \\ & \to - \int_{\Omega} \overline{\varrho} (\partial_{\varrho}s(\overline{\varrho},\overline{\vartheta})\varrho_0^{(1)}) \phi(0,\cdot)\,dx. \end{aligned} \end{equation} In particular \begin{gather} \label{initentmech0} [(\varrho s)_{0,\varepsilon}]_{\rm ess}\to\overline{\varrho} s(\overline{\varrho},\overline{\vartheta})\quad\varepsilon\to 0^+, \\ \label{initentrad0} [s^R_{0,\varepsilon}]_{\rm ess}\to {s^R(\overline{I})}\quad\varepsilon\to 0^+, \end{gather} weakly-(*) in $ L^{\infty}(0,T;L^2(\Omega))$. Residual parts of the initial conditions disappear thanks to the $L^{\infty}$ weak star convergence of the initial data in Theorem \ref{conv} for $\varepsilon$ sufficiently small. For the convergence in advective part of the entropy balance \eqref{entro} we use \eqref{tempstrong} and the fact that \begin{equation} \label{densstrong} \varrho_\varepsilon\to\overline{\varrho} \quad\varepsilon\to 0^+ \end{equation} in $L^{\infty}(0,T;L^{5/3}(\Omega))$. This allows to make the limit of entropy to a constant for a subsequence \begin{equation} \label{sconv0} s(\varrho_\varepsilon,\vartheta_\varepsilon)\to s(\overline{\varrho},\overline{\vartheta}) \quad\varepsilon\to 0^+ \text{ a. e. in }(0,T)\times\Omega. \end{equation} The convergence of entropy of a photon gas follows from Proposition \ref{prop2} as \begin{gather} \label{sRconv1} \varepsilon \frac{[ s^R_{\varepsilon}]_{\rm ess}-\overline{s^R}}{\varepsilon} \to 0,\\ \label{sRconv0} s^R_{\varepsilon}\to\overline{s^R} \end{gather} weakly-(*) in $L^{\infty}(0,T;L^2(\Omega))$ as $\varepsilon\to0+$ according to \eqref{sRres} and \eqref{Mres} again. The convergence of the next term containing $\varrho_\varepsilon s_\varepsilon \vec u_\varepsilon$ is again split into two terms, first one on the residual, second one on the essential set. For the second one we use again Proposition \ref{prop2}, the first one \begin{equation} \label{entfluxconv0res} \int_0^T \int_{\Omega}[\varrho_\varepsilon s_\varepsilon(\varrho_\varepsilon,\vartheta_\varepsilon)]_{\rm res}\vec{u}_\varepsilon\cdot \nabla_x \varphi\,dx\,dt \end{equation} converges to 0 just in $L^1((0,T)\times\Omega)$ as $\varepsilon\to0+$ because of the estimates \eqref{Mres}, \eqref{rsres}, \eqref{rtpower}. While the convergence of the equilibrial radiative entropy flux can be readily improved, e.~g.~to the space $L^\frac{12}{11}((0,T)\times\Omega)$ because of the Gibbs' relation between specific entropy and energy, cf. \eqref{m5} and \eqref{m6}, {the} integral with the material entropy flux part does not seem to have a right regularity to be meaningful. However, we can use usual cut-off functions $T_K(z) := \min(K,z)$, choose $K$ large enough, e.~g.~ $K =\varepsilon^{-\frac16}$ and split the integral to two parts \begin{align*} &\int_0^T \int_{\Omega}|[\varrho_\varepsilon s_\varepsilon(\varrho_\varepsilon,\vartheta_\varepsilon)]_{\rm res}| |\vec{u}_\varepsilon| |\nabla_x \varphi|\,dx\,dt \\ &= \int_0^T \int_{\Omega}|[\varrho_\varepsilon s_\varepsilon(\varrho_\varepsilon,\vartheta_\varepsilon)]_{\rm res} |T_{\varepsilon^{-\frac16}}(|\vec{u}_\varepsilon|) |\nabla_x \varphi|\,dx\,dt \\ &\quad + \int_0^T \int_{\Omega}|[\varrho_\varepsilon s_\varepsilon(\varrho_\varepsilon,\vartheta_\varepsilon)]_{\rm res} |[|\vec{u}_\varepsilon|-T_{\varepsilon^{-\frac16}}(|\vec{u}_\varepsilon|)]| \nabla_x \varphi|\,dx\,dt. \end{align*} The first part converges to 0 by \eqref{rsres}, the second one is of order $O(\varepsilon)$ by Sobolev embedding, estimate \eqref{Mres} and Markov-Chebyshev inequality. The limiting part of this estimate is the first part, where the need to improve the regularity the material part of the entropy flux faces the problem that we have not got generally a better estimate than \eqref{rsres}. Previous works \cite{NRT,K2,FN} rely on the closedness of the equation of state to the ideal gas law so that $\varrho_\varepsilon s_\varepsilon$ is estimated essentially by $\varrho_\varepsilon|\log\varrho_\varepsilon|$, $\vartheta_\varepsilon^3$ and $\varrho_\varepsilon|\log\vartheta_\varepsilon|$, the last one being the most restrictive, leading to the convergence in \eqref{entfluxconv0res} in $L^2(0,T;L^{\frac{30}{29}}(\Omega))$. Without such an assumption we would estimate the entropy by $\varrho_\varepsilon^2\vec{u}_\varepsilon$ which is not tractable in view of \eqref{rtpower}. Nevertheless, in our case of low stratification we do not need to identify the limit of the entropy flux on the essential set since it vanishes after an integration by parts. From \eqref{tempstrong} and \eqref{densstrong}, \begin{equation} \label{entfluxconv0ess} \int_0^T \int_{\Omega}[\varrho_\varepsilon s_\varepsilon(\varrho_\varepsilon,\vartheta_\varepsilon)]_{\rm ess}\vec{u}_\varepsilon\cdot \nabla_x \varphi\,dx\,dt\to \int_0^T \int_{\Omega}\overline{\varrho} s(\overline{\varrho},\overline{\vartheta})\vec U\cdot \nabla_x \varphi\,dx\,dt = 0. \end{equation} The last term contains the non-equlibirial radiative entropy flux $\vec{q}^R_\varepsilon$. Let us recall \[ {{\vec q}^R}_{\varepsilon}={\vec q}^R(I_{\varepsilon}) =-\int_0^{\infty}\int_{\mathcal{S}^2} \nu^2 \{n_{\varepsilon}\log n_{\varepsilon}-(n_{\varepsilon}+1)\log(n_{\varepsilon}+1) \} \vec \omega\,d\vec \omega\,d\nu, \] with $n_{\varepsilon}=n(I_{\varepsilon}) =\frac{I_{\varepsilon}}{\nu^3}$. We claim \begin{equation} \label{nonentfluxconv0} J_\varepsilon := \int_0^T\int_\Omega\vec{q}^R_\varepsilon\cdot\nabla_x\varphi\,dx\,dt \to \int_0^T\int_\Omega\overline{\vec{q}}^R\cdot\nabla_x\varphi\,dx\,dt =: J_0 \end{equation} because of the convergence on the essential set $\mathcal M_{\rm ess}^\varepsilon$ that follows from \eqref{nonentflux1conv} and on the residual set $\mathcal M_{\rm res}^\varepsilon$ we use \eqref{sRres}. Collecting now all the aforementioned convergences in this section we readily get the weak formulation of \eqref{enerzero}. With \eqref{zeromeantemp} we see that $\vartheta_2 \equiv 0$ and search for $\nabla_x\Theta = \nabla_x\vartheta_3 := w-\lim_{L^{4/3}((0,T)\times\Omega), \varepsilon\to0+}\nabla_x\frac{\vartheta_\varepsilon-\overline{\vartheta}}{\varepsilon^3}$. Let us recall that $r=3$ that is needed for the proof of existence and try to relax it for the current proof of convergence $\varepsilon\to0+$ and realize that we can extract from \eqref{sigres} the bound \begin{equation} \label{tempgradwithr} \int_0^T\int_\Omega \vartheta_\varepsilon^{r-2}\frac{|\nabla_x(\vartheta_\varepsilon-\overline{\vartheta})|^2}{\varepsilon^4}\,dx\,dt < c \end{equation} with a constant $c$ independent of $\varepsilon$. Therefore for $r\ge2$, $|\frac{\nabla_x(\vartheta_\varepsilon-\overline{\vartheta})}{\varepsilon^3}|$ is bounded in $L^{4/3}((0,T)\times\Omega)$ and $\Theta$ exists. We subtract equation \eqref{entro} from its limit and divide by $\varepsilon$, \begin{equation} \label{entconv1} \begin{aligned} &\int_0^T \int_{\Omega} \big\{ \varrho_{\varepsilon}\ \frac{ s_{\varepsilon} -\overline{s}}{\varepsilon}(\partial_t {\varphi}+ \vec{u}_{\varepsilon} \cdot \nabla_x \varphi ) + \frac{s^R_{\varepsilon}-\overline{s^R}}{\varepsilon} \varepsilon\partial_t \varphi +\frac{ {{\vec q}^R}_{\varepsilon} - \overline{{\vec q}^R}}{\varepsilon}\cdot \nabla_x \varphi\big\} \,dx \,dt\\ &+\int_0^T \int_{\Omega}\frac{\kappa(\vartheta_{\varepsilon})} {\vartheta_{\varepsilon}}\nabla_x \frac{ \vartheta_{\varepsilon} }{\varepsilon^3}\cdot \nabla_x \varphi \,dx \,dt +\frac{1}{\varepsilon}\langle \varsigma^m_{\varepsilon}+\varsigma^R_{\varepsilon} ; {\varphi}\rangle _{[{\mathcal M};C]([0,T]\times \overline{\Omega})} \\ &= - \int_{\Omega} \big\{ \big( \varrho_{0,\varepsilon}\frac{ s_{0,\varepsilon}-\overline{s}}{\varepsilon} +\varepsilon\frac{s^R_{0,\varepsilon}-\overline{s_0^R}}{\varepsilon} \big) \varphi (0, \cdot)\big\}\,dx. \end{aligned} \end{equation} We claim that all the terms in \eqref{entconv1} are uniformly bounded, especially \begin{equation} \label{heatdiffconv} \frac{\kappa(\vartheta_{\varepsilon})}{\vartheta_{\varepsilon}} \nabla_x\frac{\vartheta_{\varepsilon}-\overline{\vartheta}}{\varepsilon^3} \to \frac{\kappa(\overline{\vartheta})}{\overline{\vartheta}}\ \nabla_x \vartheta_3, \end{equation} weakly in $L^{\frac{6r+16}{6r+15}}((0,T)\times\Omega;\mathbb{R}^3)$ which gives for $r=3$ the summability with the exponent of $\frac{34}{33}$. To show this we restrict ourselves to the residual set $\mathcal M_{\rm res}^\varepsilon$, since on the essential set $\mathcal M_{\rm ess}^\varepsilon$ the boundedness is easy. For the set $A_\varepsilon := \{(t,x) : |\nabla_x\vartheta_\varepsilon(t,x)|<1\}$ we use the estimates \eqref{Mres}, \eqref{rtpower} with H\"older's inequality and $r\in[3,5]$, \begin{equation} \label{estmresanda} \begin{aligned} K_0:=&\int_{\mathcal M_{\rm res}^\varepsilon\cap A_\varepsilon} \int|[\vartheta_{\varepsilon}]_{\rm res}^{r-1} \nabla_x \frac{ \vartheta_{\varepsilon} }{\varepsilon^3}|\,dx\,dt \\ \le& \varepsilon^{-3}\int\int_{\mathcal M_{\rm res}^\varepsilon} \vartheta_\varepsilon^{r-1}\,dx\,dt\\ \le& T\varepsilon^{-3}\|[\vartheta_\varepsilon]_{\rm res}\|_{L^4(\Omega)} ^{r-1} \operatorname{ess\,sup}\|\mathbb{I}_{\mathcal M_{\rm res}^\varepsilon(t)}\|_{L^{\frac{4}{5-r}} (\Omega)} \\ \le& C \varepsilon^{-3}\varepsilon^{\frac{r-1}2}\varepsilon^{2} = C \varepsilon^{\frac{r-3}2} \le c \end{aligned} \end{equation} with $c$ independent of $\varepsilon$. In the opposite case (the complement of this set in $\mathcal M_{\rm res}^\varepsilon$) we estimate as follows \begin{equation} \label{heatdiffbound} \begin{aligned} K_1:=&\int_{\mathcal M_{\rm res}^\varepsilon\setminus A_\varepsilon} \int|[\vartheta_{\varepsilon}]_{\rm res}^{r-1} \nabla_x \frac{ \vartheta_{\varepsilon} }{\varepsilon^3}|\,dx\,dt\\ =& \int_{\mathcal M_{\rm res}^\varepsilon\setminus A_\varepsilon} \int\Big|[\vartheta_{\varepsilon}]_{\rm res}^{\frac{r}4+\frac12} |\nabla_x\vartheta_\varepsilon|^{-\frac12}\underbrace{\vartheta_\varepsilon^{\frac34 r-\frac32} \frac{|\nabla_x \vartheta_{\varepsilon}|^{\frac32} }{\varepsilon^3}}_{\in L^{4/3} ({\mathcal M_{\rm res}^\varepsilon\setminus A_\varepsilon})}\Big|\, dx\,dt\\ \le & %\label{heatdiffbound2} \int_0^T\int_{\Omega}|[\vartheta_{\varepsilon}]_{\rm res}^{\frac{r}4+\frac12} \vartheta_\varepsilon^{\frac34 r-\frac32} \frac{|\nabla_x \vartheta_{\varepsilon}|^{\frac32} }{\varepsilon^3}|\, dx\,dt < c \end{aligned} \end{equation} provided $[\vartheta_\varepsilon]_{\rm res}^{\frac{r}4+\frac12}$ is uniformly bounded in $L^{4}((0,T)\times\Omega)$; that is, $[\vartheta_\varepsilon]_{\rm res}^{{r}+2}$ is uniformly bounded in the space $L^{1}((0,T)\times\Omega)$. However we know that $[\vartheta_\varepsilon]_{\rm res}$ is bounded in $L^{\infty}(0,T;L^{4}(\Omega))\cap L^{r}(0,T;L^{3r}(\Omega))$ as in \cite{FEINOV}. By interpolation we obtain $[\vartheta_\varepsilon]_{\rm res}$ is uniformly bounded in $L^{r+\frac83}((0,T)\times\Omega)$ and that is why $K_1$ converges; moreover when we reiterate the same argument with a $s$-power of its integrand, we obtain the bound $s\le\frac{6r+16}{6r+15}$ for H\"older's inequality. Similarly to \cite{FEINOV}, using Proposition \ref{prop2} and energy estimates, we see that \[ \varrho_{\varepsilon}\ \frac{ s_{\varepsilon}-\overline{s}}{\varepsilon} \to \overline{\varrho} (\partial_{\varrho}s(\overline{\varrho},\overline{\vartheta}) \varrho^{(1)}+\partial_{\vartheta}s(\overline{\varrho},\overline{\vartheta})\vartheta^{(1)}) = \overline{\varrho}\varrho_1\partial_{\varrho}s(\overline{\varrho},\overline{\vartheta}) \] weakly-* in $L^{\infty}(0,T;L^2(\Omega;\mathbb{R}^3))$, and remind \eqref{mechradentrprodconv1}, \eqref{sRconv1}, \eqref{sRres}, \eqref{Mres} and \eqref{initent1}. Moreover the advective part weakly converges according to Proposition \ref{prop2} again \[ \varrho_{\varepsilon}\ \frac{ s_{\varepsilon}-\overline{s}}{\varepsilon} \vec u_{\varepsilon} \to \overline{\varrho} (\partial_{\varrho}s(\overline{\varrho},\overline{\vartheta}) \varrho^{(1)}+\partial_{\vartheta}s(\overline{\varrho},\overline{\vartheta})\vartheta^{(1)}) \vec U = \overline{\varrho} \partial_{\varrho}s(\overline{\varrho},\overline{\vartheta}) \varrho_1 \vec U, \] weakly in $L^2(0,T;L^{3/2}(\Omega;\mathbb{R}^3))$. This allows to pass to the limit in all terms of \eqref{entconv1} except the nonequilibrial radiative entropy flux term \begin{equation} \label{nonentflux1} \int_0^T \int_{\Omega}\frac{ {{\vec q}^R}_{\varepsilon}- \overline{{\vec q}^R}}{\varepsilon}\cdot \nabla_x \varphi \,dx \,dt. \end{equation} Let us compute the limit of $\frac{ {{\vec q}^R}_{\varepsilon}- \overline{{\vec q}^R}}{\varepsilon}$. Applying once more Proposition \ref{prop2} with $G^R(I)=n(I)\log n(I)-(n(I)+1)\log(n(I)+1)$ and integrating on $\mathcal{S}^2\times\mathbb{R}_+$, we find \[ \frac{ {{\vec q}^R}_{\varepsilon}- \overline{{\vec q}^R}}{\varepsilon} \to \int_0^{\infty}\int_{\mathcal{S}^2} \frac{1}{\nu} \log\frac{n(\overline{I})+1}{n(\overline{I})} \vec \omega\ I^{(1)}\,d\vec \omega\,d\nu, \] weakly-* in $L^{\infty}(0,T;L^2(\Omega;\mathbb{R}^3))$ on the essential set $\mathcal M_{\rm ess}^\varepsilon$ and as $\log[\frac{n(\overline{I})+1}{n(\overline{I})}]=\frac{\nu}{\overline{\vartheta}}$, we have \[ \frac{ {{\vec q}^R}_{\varepsilon}-\overline{ {\vec q}^R } }{\varepsilon} \to \frac{1}{\overline{\vartheta}}\ {\vec F}^R(I^{(1)}), \] with the radiative momentum ${\vec F}^R(I^{(1)})=\int_0^{\infty}\int_{\mathcal{S}^2} \vec \omega\ I^{(1)}\,d\vec \omega\,d\nu$. So \begin{equation} \label{nonentflux1conv} \int_0^T \int_{\Omega} \big(\frac{ {{\vec q}^R}_{\varepsilon}- \overline{{\vec q}^R}}{\varepsilon} \big)\cdot \nabla_x \varphi\,\,dx\,dt \to \int_0^T \int_{\Omega} \frac{\operatorname{div}_x {\vec F}^R(I^{(1)})}{\overline{\vartheta}} \varphi\,\,dx\,dt \end{equation} by the Proposition \ref{prop2}, \eqref{sRres} and \eqref{Mres} on $\mathcal M_{\rm res}^\varepsilon$. As we have from \eqref{transfer1}, \[ \operatorname{div}_x {\vec F}^R =\int_0^{\infty}\int_{\mathcal{S}^2} [\partial_{\vartheta}\sigma_a(\nu,\overline{\vartheta})(B(\nu,\overline{\vartheta})-I_0)\vartheta^{(1)} +\sigma_a(\nu,\overline{\vartheta})(\partial_{\vartheta}B(\nu,\overline{\vartheta})\vartheta^{(1)}-I_1) ]\,d\vec \omega\,d\nu, \] the limit contribution in \eqref{entconv1} becomes \[ \int_0^T \int_{\Omega}\int_0^{\infty}\int_{\mathcal{S}^2} \frac{-\sigma_a(\nu,\overline{\vartheta}) I_1(t,x,\vec \omega,\nu)}{\overline{\vartheta}} \varphi \,d\vec \omega\,d\nu\ \,dx \,dt. \] Gathering all these terms, we find the limit equation for entropy \begin{align*} &\overline{\varrho} \partial_{\varrho}s(\overline{\varrho},\overline{\vartheta}) \int_0^T \int_{\Omega} \varrho_1 \big(\partial_t \phi+\vec U\cdot \nabla_x \phi \big)\,dx\,dt +\frac{\overline{\kappa}}{\overline{\vartheta}}\int_0^T \int_{\Omega} \nabla_x \Theta\cdot\nabla_x \varphi\,dx\,dt \\ &\quad -\frac1{\overline{\vartheta}}\int_0^T \int_{\Omega}\int_0^{\infty} \sigma_a(\nu,\overline{\vartheta})\int_{\mathcal{S}^2} I_1(t,x,\vec \omega,\nu) \varphi\,d\vec \omega\,d\nu \,dx \,dt \\ &=- \overline{\varrho} \partial_{\varrho}s(\overline{\varrho},\overline{\vartheta})\int_{\Omega} \varrho_0^{(1)} \varphi(0,\cdot)\,dx. \end{align*} Using \eqref{press2} we easily verify that we finally obtained the thermal equation \eqref{bous3} once we take the Maxwell relation $\partial_\vartheta p = \partial_\varrho s$ into account. \subsection{Maxwell equation} From \eqref{convu} and \eqref{convB} we obtain \begin{gather*} \frac{\vec{B}_{\varepsilon}}{\varepsilon} \times \vec u \to \vec{B}\times \vec U \quad\text{weakly in $L^q(0,T; L^q(\Omega, \mathbb{R} ^3))$ for } q\in [1,\frac53), \\ \lambda(\vartheta_{\varepsilon}) \operatorname{curl}_x \frac{\vec{B}_{\varepsilon}}{\varepsilon} \to \overline{\lambda} \operatorname{curl}_x\vec{B} \quad\text{weakly in } L^{\frac{34}{6p+17}}(0,T,L^{\frac{34}{6p+17}}(\Omega, \mathbb{R} ^3)). \end{gather*} Then it is easy to pass to the limit in \eqref{magfield}, realizing that ${\frac{34}{6p+17}}>1$ for {$1\le p<\frac{17}{6}$}. This completes the proof of Theorem \ref{conv}. \section{Appendix: Proof of Theorem \ref{ROB}} Consider now the linearly coupled problem for the remaining equations \begin{gather} \label{b1} \operatorname{div}_x \vec U = 0, \\ \label{b2} \partial_t \vec U + (\vec U\cdot\nabla_x)\vec U +\ \nabla_x \Pi - \overline{\mu}\Delta\vec U+ \frac{1}{\zeta} \nabla_x(\frac{\vec B^2}{2})- \frac{1}{\zeta}( \vec B\cdot\nabla_x) \vec B=\vec F, \\ \label{b3} \partial_t \vec{B} + (\vec U\cdot\nabla_x) \vec B+(\vec B\cdot\nabla_x) \vec U - \overline{\lambda}\Delta \vec{B}= 0, \\ \label{b4} \operatorname{div}_x \vec B = 0, \\ \label{b5} -\Delta\Theta =\vec U\cdot\vec\beta - \frac{1}{\overline{\kappa}}\int_0^{\infty}\sigma_{a} \int_{\mathcal{S}^2}I_1\,d\vec \omega\,d\nu + \tilde{h} \\ \label{b6} \vec \omega \cdot \nabla_x I_1+\sigma_{a} I_1-\sigma_{s}(\tilde I_1 - I_1)=0, \end{gather} where $\vec\beta\in (L^{\infty}(\Omega))^3$, together with the boundary conditions \begin{equation} \label{b7} \vec U|_{\partial \Omega} = 0, \quad \nabla \Theta \cdot \vec{n} |_{\partial \Omega} = 0, \quad \vec B \cdot \vec{n}|_{\partial \Omega} = 0, \quad \operatorname{curl}_x\vec B \times \vec{n}|_{\partial \Omega} = 0 \end{equation} for \eqref{b1}--\eqref{b5} and \begin{equation} \label{b8} I_1 (x,\nu, \vec \omega) = 0 \quad \text{for } x \in \partial \Omega, \; \vec \omega \cdot \vec{n} \leq 0 \end{equation} for \eqref{b6}, and the initial conditions \begin{equation} \label{b9} \vec U|_{t=0} = \vec U_0,\quad \vec B|_{t=0} = \vec B_0. \end{equation} We first consider the solution $(\vec U,\vec B, I_1)$ of the ``radiative-MHD problem'' \begin{gather} \label{b11} \operatorname{div}_x \vec U = 0, \\ \label{b21} \partial_t \vec U + (\vec U\cdot\nabla_x)\vec U +\ \nabla_x \Pi - \overline{\mu}\Delta\vec U=\frac{1}{\zeta} \operatorname{curl}_x\vec B\times \vec B + \vec F, \\ \label{b31} \partial_t \vec{B} + (\vec U\cdot\nabla_x) \vec B + \vec B\cdot\nabla_x\vec U - \overline{\lambda}\Delta \vec{B}= 0, \\ \label{b41} \operatorname{div}_x \vec B = 0, \\ \label{b61} \vec \omega \cdot \nabla_x I_1+\sigma_{a} I_1-\sigma_{s}(\tilde I_1 - I_1)=0, \end{gather} with \begin{gather*} \vec U|_{\partial \Omega} = 0, \quad \vec B \cdot \vec{n}|_{\partial \Omega} = 0, \quad \operatorname{curl}_x\vec B \times \vec{n}|_{\partial \Omega} = 0, \\ \vec U|_{t=0} = \vec U_0,\quad \vec B|_{t=0} = \vec B_0,\quad I_1 (x,\nu, \vec \omega) = 0 \quad \text{for } x \in \partial \Omega, \; \vec \omega \cdot \vec{n} \leq 0. \end{gather*} The MHD part has a weak solution $\vec U\in L^2(0,T;{\mathcal U}(\Omega))$, $\vec B\in L^2(0,T;{\mathcal W}(\Omega))$ thanks to an extension of the Leray-Hopf Theorem (see \cite{ST}). Moreover the stationary radiative equation \eqref{b61} also has a weak solution $I_1\in L^2((0,T)\times \Omega\times\mathcal{S}^2\times\mathbb{R}_+) $ according to \cite[Theorem 1 and Proposition 2]{BGPS}. Then we consider the solution $\Theta$ of the stationary diffusion equation \begin{equation} \label{b51} -\Delta\Theta =\vec U\cdot\vec\beta - \frac{1}{\overline{\kappa}}\int_0^{\infty}\sigma_{a}\int_{\mathcal{S}^2}I_1\,d \vec \omega\,d\nu + \tilde{h} \end{equation} with \[ \nabla \Theta \cdot \vec{n} \big|_{\partial \Omega} = 0 \] subject to $\int_\Omega \Theta\,dx = 0$ for all times. This problem admits a weak solution $\Theta\in L^{\infty}((0,T; W^{2,2}(\Omega)) \cap L^{2}((0,T; W^{q,2}(\Omega))$ for all $q<5/2$, thanks to classical elliptic regularity theory and due to regularity of the right-hand side due to \cite{GLPS}. Since the ``radiative-MHD problem'' does not depend on the temperature perturbation $\Theta$ the proof is complete. \subsection*{Acknowledgments} \v{S}. N. was supported by the Czech Agency of the Czech Republic No. 16-03230S and by RVO 67985840. Bernard Ducomet is partially supported by the ANR project INFAMIE (ANR-15-CE40-0011). \begin{thebibliography}{00} \bibitem{BAL} R. Balian: \newblock \emph{From microphysics to macrophysics. Methods and applications of statistical physics Vol. II}. \newblock Springer Verlag, Berlin, Heidelberg, New York, 1992. \bibitem{BGPS} C. Bardos, F. Golse, B. Perthame, R. Sentis: \newblock \emph{The nonaccretive radiative transfer equations: Existence of solutions and Rosseland approximation}. \newblock J. Funct. Anal. 77 (2) (1988) 434--460. \bibitem{Bro} A.-M.~Broomhall, P. Chatterjee, R. Howe, A. A. Norton, M. J. Thompson: \newblock \emph{The Sun's interior structure and dynamics, and the solar cycle}. \newblock \emph{ArXiv:1411.5941v2 [astro-ph.SR] 25 Nov 2014}. \bibitem{C} A.R. Choudhuri: \newblock \emph{The physics of fluids and plasmas, an introduction for astrophysicists}. \newblock Cambridge University Press, 1998. \bibitem{ddn} D. Donatelli, B. Ducomet, \v{S}. Ne\v{c}asov\'a: \newblock \emph{Low Mach number limit for a model of accretion disk}, \newblock Submitted \bibitem{DL} R.J. DiPerna, P.-L. Lions: \newblock \emph{Ordinary differential equations, transport theory and {S}obolev spaces}. \newblock Invent. Math. 98 (1989), 511--547. \bibitem{DF} B. Ducomet, E. Feireisl: \newblock \emph{The equations of magnetohydrodynamics: on the interaction between matter and radiation in the evolution of gaseous stars}. \newblock Commun. Math. Phys., 266 (2006), 595--629. \bibitem{DFN} B. Ducomet, E. Feireisl, \v{S}. Ne\v{c}asov\'a: \newblock \emph{On a model of radiation hydrodynamics}. \newblock Ann. I. H. Poincar\'e AN-28 (2011), 797--812. \bibitem{DFPS} B. Ducomet, E. Feireisl, H. Petzeltov\'a, I. Stra\v skraba: \newblock \emph{Global in time weak solutions for compressible barotropic self-gravitating fluids}. \newblock Discrete and Continuous Dynamical Systems, 11 (2004), 113--130. \bibitem{dkn} B. Ducomet, M. Kobera, \v{S}. Ne\v{c}asov\'a: \newblock \emph{Global existence of a weak solution for a model in radiation magnetohydrodynamics}. \newblock Submitted. \bibitem{DuNe} B. Ducomet, \v{S}. Ne\v{c}asov\'a: \newblock \emph{Low Mach number limit for a model of radiative flow}. \newblock J. of Evol. Equ., 14 (2014), 357--385. \bibitem{dnn} B. Ducomet, M. Caggio, \v{S}. Ne\v{c}asov\'a, M. Pokorn\'y: \newblock \emph{The rotating Navier-Stokes-Fourier-Poisson system on thin domains}. \newblock Submitted. \bibitem{dnnp} B. Ducomet, \v{S}. Ne\v{c}asov\'a, M. Pokorn\'y: \newblock \emph{Thin domain limit for a polytropic model of radiation hydrodynamics}. \newblock In preparation. \bibitem{DULI} G.~Duvaut, J.-L. Lions: \newblock \emph{Inequalities in mechanics and physics}. \newblock Springer-Verlag, Heidelberg, 1976. \bibitem{FEINOV} E.~Feireisl, A.~Novotn\'y: \newblock \emph{Singular limits in thermodynamics of viscous fluids}. \newblock Birkhauser, Basel, 2009. \bibitem{FN} E.~Feireisl, A.~Novotn\'y: \newblock \emph{Small P\'eclet approximation as a singular limit of the full Navier-Stokes-Fourier equation with radiation}. \newblock Preprint. \bibitem{Gar} P.~Garaud, L.~Kulenthirarajah: \newblock \emph{Turbulent transport in a strongly stratified forced shear layer with thermal diffusion}. \newblock ArXiv:1512.08774v1 [astro-ph.SR] 29 Dec 2015. \bibitem{GLPS} F.~Golse, P. L. Lions, B. Perthame, R. Sentis: \emph{Regularity of the moments of the solution of a transport equation}, J. Funct. Anal. 16 (1988) 110--125. \bibitem {HW} X Hu, D. Wang: \newblock \emph{Low Mach number limit of viscous compressible magnetohydrodynamic flows}. \newblock SIAM J. Math. Anal. 41 (2009), no. 3, 1272--1294. \bibitem{JJLX} S. Jiang, Q. Ju, F. Li, Z. Xin: \newblock \emph{Low Mach number limit for the full compressible magnetohydrodynamic equations with general initial data}. \newblock Adv. Math. 259 (2014), 384--420. \bibitem{JJL} S. Jiang, Q. Ju, F. Li: \newblock \emph{Low Mach number limit for the multi-dimensional full magnetohydrodynamic equations}. \newblock Nonlinearity 25 (2012), no. 5, 1351--1365. \bibitem{JJL1} S. Jiang, Q. Ju, F. Li: \newblock \emph{Incompressible limit of the compressible magnetohydrodynamic equations with vanishing viscosity coefficients}. \newblock SIAM J. Math. Anal. 42 (2010), no. 6, 2539--2553. \bibitem{K} P. Kuku\v{c}ka: \newblock \emph{Singular limits of the equations of magnetohydrodynamics}. \newblock J. Math. Fluid Mech. 13 (2011) 173--189. \bibitem{K2} P. Kuku\v{c}ka: \newblock \emph{Incompressible limits for the Navier-Stokes-Fourier Systems on Unbounded Domains under Strong Stratification}. \newblock Preprint. \bibitem{TK} Y. Kwon, K. Trivisa: \newblock \emph{On the incompressible limits for the full magnetohydrodynamics flows}. \newblock J. Differential Equations 251 (2011), no. 7, 1990--2023. \bibitem{Lig} F.~Ligni\`eres: \newblock \emph{The small-P\'eclet-number approximation in stellar radiative zones}. \newblock ArXiv: 9908.182v1 [astro-ph] 17 Aug 1999. \bibitem{LSU} O. A. Lady\v{z}enskaja, V. A. Solonnikov, N. N. Ural'ceva: \newblock \emph{Linear and quasilinear equations of parabolic type}. \newblock AMS 1968. \bibitem{LL} E. H. Lieb, M. Loss: \newblock \emph{Analysis}, \newblock Graduate Studies in Mathematics, Vol. 14, A.M.S., 1997. \bibitem{NR} J. Ne\v{c}as, T. Roub\'i\v{c}ek: \newblock \emph{Buoyancy-driven viscous flow with $L^1$-data}. \newblock Nonlinear Analysis 46 (2001) 737--755. \bibitem{NS} A. Novotn\'y, I. Stra\v skraba: \newblock \emph{Introduction to the mathematical theory of compressible flows}, \newblock Oxford Lectures Series in Mathematics and Its Applications, Vol. 27, Oxford University Press, 2004. \bibitem{NRT} A. Novotn\'y, M. R\accent23 u\v{z}i\v{c}ka, G. Th\"ater: \newblock \emph{Singular limit of the equations of magnetohydrodynamics in the presence of strong stratification}. \newblock Math. Models Methods Appl. Sci. 21 (2011), no. 1, 115--147. \bibitem{P} T. Padmanabhan: \newblock \emph{Theoretical astrophysics, Vol II: Stars and stellar systems}, \newblock Cambridge University Press, 2001. \bibitem{Rie} M.~Rieutord: \newblock \emph{Magnetohydrodynamics and solar physics}. \newblock ArXiv:1410.3725v2 [astro-ph.SR] 15 Oct 2014. \bibitem{Rou} T. Roub\'i\v{c}ek: \newblock \emph{Nonlinear partial differential equations with applications}, \newblock Bikh\"auser, 2005. \bibitem{S} S. Schochet: \newblock \emph{Fast singular limits of hyperbolic PDE's}, \newblock J. Diff. Eq., 114 (1994) 476 -- 512. \bibitem{ST} M. Sermange, R. Temam: \newblock \emph{Some mathematical questions related to magnetohydrodynamics}, \newblock Comm. Pure Appl. Math., 36 (1983) 635--664. \bibitem{Sh} S.N. Shore: \newblock \emph{An introduction to astrophysical hydrodynamics}, \newblock Academic Press, 1992. \bibitem{TSGKS} I. Teleaga, M. Sea\"id, I. Gasser, A. Klar, J. Struckmeier, \emph{Radiation models for thermal flows at low Mach number}, J. of Comput. Phys. 215 (2006) 506--525. \bibitem{TEM} R. Temam: \newblock \emph{{N}avier-{S}tokes equations}. \newblock North-Holland, Amsterdam, 1977. \end{thebibliography} \end{document}