\documentclass[reqno]{amsart} \usepackage{hyperref} \AtBeginDocument{{\noindent\small \emph{Electronic Journal of Differential Equations}, Vol. 2013 (2013), No. 232, pp. 1--29.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu} \thanks{\copyright 2013 Texas State University - San Marcos.} \vspace{9mm}} \begin{document} \title[\hfilneg EJDE-2013/232\hfil Motion of compressible magnetic fluids] {Motion of compressible magnetic fluids in $\mathbb{T}^3$} \author[W. Yan \hfil EJDE-2013/232\hfilneg] {Weiping Yan} % in alphabetical order \address{Weiping Yan \newline College of of Mathematics, Jilin University, Changchun 130012, China} \email{yan8441@126.com} \thanks{Submitted August 8, 2013. Published October 18, 2013.} \subjclass[2000]{76W05, 36D30, 35B10} \keywords{Magnetohydrodynamics; weak solution; periodic solution} \begin{abstract} This article shows the existence of weak time-periodic motion of a three-dimensional system of compressible magnetic fluid driven by time-dependent external forces in a torus $\mathbb{T}^3$. The model consists of the mass conservation equation, the linear momentum equation, the angular momentum equation, the Bloch-Torrey type equation and the magnetostatic equation. This analysis is based on the Faedo-Galerkin method and weak compactness techniques. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \newtheorem{proposition}[theorem]{Proposition} \allowdisplaybreaks \section{Introduction}\label{intro} Magnetic fluids (also called ferrofluids) have found a wide variety of applications: magnetic liquid seals, cooling and resonance damping for loudspeaker coils, printing with magnetic inks, magnetic resonance imaging contrast enhancement, magnetic separation, rotating shaft seals in vacuum chambers used in the semiconductor industry, for more applications see \cite{Grief,Pankhurst,Zahn}. Magnetic fluids are colloidal suspensions of fine magnetic mono domain nanoparticles in nonconducting liquids. A set of such system is a combination of the compressible Navier-Stokes equations, and the angular momentum equation, the Bloch-Torrey equation and the magnetostatic equations, see \cite{Ami,Rosensweig}. The Bloch-Torrey equation is proposed by Torrey \cite{Torrey} as a generalization of the Bloch equations, which is the magnetization equation with a diffusion term describing the situations when the diffusion of the spin magnetic moment is not negligible. More precisely, we consider the following model: the continuity equation, \begin{equation}\label{E1-1} \partial_t\rho+\operatorname{div}(\rho\mathbf{u})=0, \end{equation} the linear momentum equation \begin{equation} \label{E1-2} \partial_t(\rho\mathbf{u})+\operatorname{div}(\rho\mathbf{u}\otimes \mathbf{u})-\mu\Delta\mathbf{u}-(\lambda+\mu) \nabla(\operatorname{div}\mathbf{u})+\nabla(P(\rho,\mathbf{M})) =R+\rho\mathbf{f}(t,x),\ \end{equation} the angular momentum equation \begin{equation}\label{E1-3} \partial_t(\rho\Omega)+\operatorname{div}(\rho\mathbf{u}\otimes \Omega)-\mu'\Delta\Omega-(\lambda'+\mu')\nabla(\operatorname{div}\Omega) =S+\rho\mathbf{g}(t,x), \end{equation} the Bloch-Torrey equation \begin{equation}\label{E1-R1} \partial_t\mathbf{M}+\operatorname{div}(\mathbf{u}\otimes \mathbf{M})-\sigma\Delta\mathbf{M}+\frac{1}{\tau} (\mathbf{M}-\chi_0\mathbf{H})=\Omega\times\mathbf{M}+\rho\mathbf{l}(t,x), \end{equation} the magnetostatic equation \begin{equation}\label{E1-R5} \mathbf{H}=\nabla\psi,\quad \operatorname{div}(\mathbf{H}+\mathbf{M})=\mathbb{F}, \end{equation} where $(t,x)\in\mathbb{R}^1\times\mathbb{T}^3$, $\rho\in\mathbb{R}^1_+$ denotes the density, $\mathbf{u}\in\mathbb{R}^3$ denotes the fluid velocity, $\mathbf{M}\in\mathbb{R}^3$ denotes the magnetization, $\mathbf{H}\in\mathbb{R}^3$ denotes the magnetic field, $\Omega\in\mathbb{R}^3$ denotes the angular velocity, $\mathbb{F}$ is a given function defined on $\mathbb{R}^1\times\mathbb{T}^3$ with $\int_{\mathbb{T}^3}\mathbb{F}dx=0$, for all $t\in\mathbb{R}^1$, \begin{gather}\label{E1-R2} R=\mu_0\mathbf{M}\cdot\nabla\mathbf{H}-\zeta\nabla\times (\nabla\times\mathbf{u}-2\Omega), \\ \label{E1-R3} S=\mu_0\mathbf{M}\times\mathbf{H}+2\zeta(\nabla\times\mathbf{u}-2\Omega), \end{gather} the pressure $P(\rho,\mathbf{M})$ obeys the state law \cite{Rosensweig}: \begin{equation}\label{E1-R4} P(\rho,\mathbf{M})=P_e(\rho)+P_m(\mathbf{M}) \end{equation} with the isentropic pressure $P_e(\rho)=a\rho^{\gamma}$, $a>0$ and the adiabatic exponent $\gamma>\frac{3}{2}$ are constants, the magnetic pressure $P_m(\mathbf{M})=\frac{\mu_0}{2}|\mathbf{M}|^2$. The parameters $\lambda$, $\mu$, $\lambda'$, $\mu'$, $\chi_0$, $\mu_0$, $\zeta$ and $\tau$ are positive and their physical meaning can be found in \cite{Rosensweig}. $\mathbf{f}(t,x)$, $\mathbf{g}(t,x)$ and $\mathbf{l}(t,x)$ denote the time periodic external forces with a period $\omega>0$. Moreover, the functions $\mathbf{f}(t,x)$, $\mathbf{g}(t,x)$ and $\mathbf{l}(t,x)$ have bounded and measurable components with no restriction on their amplitudes. The symbol $\otimes$ denotes the Kronecker tensor product. In recent years, compressible magnetohydrodynamics (MHD) has attracted much attention. Such a system (MHD) is a combination of the compressible Navier-Stokes equations of fluid dynamics and Maxwell's equations of electromagnetism. Duvaut and Lions \cite{Duv}, Sermange and Temam \cite{Ser} obtained some existence and long time behavior results; Ducomet and Feireisl \cite{Fei4} proved existence of global in time weak solutions to a multi-dimensional nonisentropic MHD system for gaseous stars coupled with the Poisson equation with all the viscosity coefficients and the pressure depend on temperature and density asymptotically, respectively. Hu and Wang \cite{Hu1} studied the global variational weak solution to the three-dimensional full magnetohydrodynamic equations with large data by an approximation scheme and a weak convergence method. In \cite{Hu2}, by using the Faedo-Galerkin method and the vanishing viscosity method, they also studied the existence and large-time behavior of global weak solutions for the three-dimensional equations of compressible magnetohydrodynamic isentropic flows \eqref{E1-1}-\eqref{E1-3}. They \cite{Hu3} showed that the convergence of weak solutions of the compressible MHD system to a weak solution of the viscous incompressible MHD system. Jiang, et all. \cite{Jiang,Jiang1} obtained that the convergence towards the strong solution of the ideal incompressible MHD system in the whole space and periodic domain, respectively. But the study of magnetic fluids differs from MHD that concerns itself with nonmagnetizable but electrically conducting fluids. In fact, it is more complicated and can be treated as homogeneous monophase fluid, see \cite{Rosensweig}. As the authors known, for compressible case without external forces, only Amirat and Hamdache \cite{Ami1,Ami} obtained the existence of global in time weak solution with finite energy to magnetic fluids equations \eqref{E1-1}-\eqref{E1-R4} under suitable the boundary conditions. A natural question of the existence of time-periodic solution arises when the external forces are time-periodic. Feireisl \cite{Fei1} first studied three dimensional compressible Navier-Stokes equations driven by a time-periodic external force. They imposed so-called no-stick boundary condition. For the three dimensional flat boundary case, this condition means that the vorticity is perpendicular (see \cite{Fei1}). Using the Faedo-Galerkin method and the vanishing viscosity method, they obtained the existence of time-periodic weak solution for three dimensional compressible Navier-Stokes equations. But the uniqueness is highly nontrivial and far from being solved. For MHD driven by the time periodic external forces, Yan and Li \cite{Yan} showed that such system has the time periodic weak solution. In what follows, we use the ideas and techniques in the books of Lions \cite{Lions3} and Feireisl \cite{Fei3} to construct a time-periodic weak solution of the problem \eqref{E1-1}-\eqref{E1-R5}. The main difficulty is how to deal with the strong coupling of the hydrodynamic motion and the magnetic field. Here, the density $\rho(t)$ is a renormalized solution of the continuity equation \eqref{E1-1} in domain $\mathcal{D}'(\mathbb{R}^1)$, which is introduced by DiPerna and Lions \cite{Duv}, it says that if for any function $b\in\mathbb{C}^1(\mathbb{R}^+)$, there holds \[ \frac{\partial b(\rho)}{\partial t}+\operatorname{div}(b(\rho) \mathbf{u})+(b'(\rho)\rho-b(\rho))\operatorname{div}\mathbf{u}=0. \] Assume that \[ \mathbf{f}(t+\omega,x)=\mathbf{f}(t,x),\quad \mathbf{g}(t+\omega,x)=\mathbf{g}(t,x),\quad \mathbf{l}(t+\omega,x)=\mathbf{l}(t,x), \] for a.e. $x\in\mathbb{T}^3,~t\in\mathbb{R}^1$, with some constant $\omega>0$ (periodic). The corresponding periodic solution of equations \eqref{E1-1}-\eqref{E1-R5} should satisfy \begin{gather*} \rho(t+\omega,x)=\rho(t,x),\quad \mathbf{u}(t+\omega,x)=\mathbf{u}(t,x),\quad \Omega(t+\omega,x)=\Omega(t,x),\\ \mathbf{M}(t+\omega,x)=\mathbf{M}(t,x),\quad \mathbf{H}(t+\omega,x)=\mathbf{H}(t,x), \end{gather*} for all $t\in\mathbb{R}^1$, $x\in\mathbb{T}^3$. We impose so-called no-stick boundary conditions to the fluid velocity, the angular velocity, the magnetization field and the magnetic field: \begin{gather*} \mathbf{u}(t,x)\cdot\nu(x)=0,\quad \Omega(t,x)\cdot\nu(x)=0,\quad \forall t\in\mathbb{R}^1,\; x\in\partial\mathbb{D},\\ \mathbf{M}(t,x)\cdot\nu(x)=0,\quad \mathbf{H}(t,x)\cdot\nu(x)=0\quad \forall t\in\mathbb{R}^1,\;x\in\partial\mathbb{D},\\ (\nabla\times\mathbf{M})\times\nu(x)=0,\quad \forall t\in\mathbb{R}^1,\; x\in\partial\mathbb{D}, \end{gather*} where $\nu(x)$ denotes the outer normal vector. Here, the domain we considered is a three dimensional torus $\mathbb{T}^3$, this no-stick boundary conditions on the fluid velocity, the angular velocity, the magnetization field and the magnetic field is equal to \begin{equation}\label{E1-A} \begin{gathered} u_i=0\quad\text{on the opposite faces}\\ \{x_i=0,x_j\in[0,\pi],j\neq i\}\cup\{x_i=\pi,x_j\in[0,\pi],j\neq i\},\\ \Omega_i=0\quad\text{on the opposite faces}\\ \{x_i=0,x_j\in[0,\pi],j\neq i\}\cup\{x_i=\pi,x_j\in[0,\pi],j\neq i\},\\ M_i=0\quad\text{on the opposite faces}\\ \{x_i=0,x_j\in[0,\pi],j\neq i\}\cup\{x_i=\pi,x_j\in[0,\pi],j\neq i\},\\ \frac{\partial M_j}{\partial x_i}=0\quad\text{for $i\neq j$ on}\\ \{x_i=0,x_j\in[0,\pi],j\neq i\}\cup\{x_i=\pi,x_j\in[0,\pi],j\neq i\},\\ H_i=0\quad\text{on the opposite faces}\\ \{x_i=0,x_j\in[0,\pi],j\neq i\}\cup\{x_i=\pi,x_j\in[0,\pi],j\neq i\}. \end{gathered} \end{equation} The density satisfies the physical requirements: \[ \rho\geq0,\quad \int_{\mathbb{T}^3}\rho(t)dx=m>0,\quad \forall t\in\mathbb{R}^1, \] for a given positive mass $m$. The rest of this article is organized as follows: In section 2, the concept of the finite energy weak time-periodic solution is introduced and the main result of this paper is also given. In section 3, using the Faedo-Galerkin method, the time periodic solution of corresponding approximation equations is obtained. By employing compactness arguments developed by Lions, et all. \cite{Lions1} and Feireisl et all. \cite{Fei1,Fei3} of compressible barotropic flows, we pass to the limit, in section 4 as $\epsilon\to 0$ and then in section 5 as $\delta\to 0$. We denote $C$ as a generic constant, depending only on some bounds of the physical data, which can take different values in different occurrences. \section{Preliminary results} In this section, we introduce some notations and the main result. Because this paper is devoted to finding the existence of weak time-periodic solution, it is convenient to consider the time $t$ belonging to the one dimensional sphere \[ t\in\mathbb{S}^1=[0,\omega]|_{[0,\omega]}. \] Moreover, we set \[ \mathbb{Q}=\mathbb{S}^1\times\mathbb{T}^3. \] The time-periodic external forces $\mathbf{f}(t,x)$, $\mathbf{g}(t,x)$ and $\mathbf{l}(t,x)$ are \begin{gather*} \mathbf{f}(t,x)=[f^1(t,x),f^2(t,x),f^3(t,x)],\\ \mathbf{g}(t,x)=[g^1(t,x),g^2(t,x),g^3(t,x)],\\ \mathbf{l}(t,x)=[l^1(t,x),l^2(t,x),l^3(t,x)]. \end{gather*} Let $\mathbb{L}^p(\mathbb{T}^3)$ and $\mathbb{W}^{p,q}(\mathbb{T}^3)$ $(1\leq p\leq\infty,1\leq q\leq\infty)$ be the usual Lebesgue and Sobolev spaces. $\mathbb{C}([0,\omega];\mathbb{X}'_{\rm weak})$ is the space of function $u:[0,\omega]\to \mathbb{X}'$ which are continuous with respect to the weak topology, where $\mathbb{X}$ is a Banach space and $\mathbb{X}'$ is the the dual space of $\mathbb{X}$. We denote \[ \mathcal{M}=\{\mathbf{u}\in\mathbb{L}^2(\mathbb{T}^3): \operatorname{div}\mathbf{u}\in\mathbb{L}^2(\mathbb{T}^3),\; \nabla\times\mathbf{u}\in(\mathbb{L}^2(\mathbb{T}^3))^3,\; \mathbf{u}\cdot\nu(x)=0\text{ on }\partial\mathbb{D}\}, \] the Hilbert space equipped with the scalar product: \[ (\mathbf{u},\mathbf{v})=\int_{\mathbb{T}^3}\mathbf{u}\cdot\mathbf{v}dx +\int_{\mathbb{T}^3}(\operatorname{div}\mathbf{u})(\operatorname{div}\mathbf{v})dx +\int_{\mathbb{T}^3}(\nabla\times\mathbf{u})(\nabla\times\mathbf{v})dx \] and the associated norm. From \cite{Dau}, we know that if we also denote \[ \mathcal{M}=\{\mathbf{u}\in\mathbb{W}^{1,2}(\mathbb{T}^3):\mathbf{u}\cdot\nu(x)=0 \text{ on }\partial\mathbb{D}\}, \] then $\|\cdot\|_{\mathcal{M}}$ and $\|\cdot\|_{\mathbb{W}^{1,2}(\mathbb{T}^3)}$ are two equivalent norms on the space $\mathcal{M}$. Now, we introduce the concept of finite energy weak solution $(\rho,\mathbf{u},\Omega,\mathbf{M},\mathbf{H})$ to the problem \eqref{E1-1}-\eqref{E1-R4}, which is taken from the strategy in \cite{Ami,Fei1,Fei3,Lions3} as (1) The non negative density function $\rho$ and the momentum $\rho\mathbf{u}$ satisfy \begin{gather*} \rho\in\mathbb{L}^{\infty}\big(\mathbb{S}^1;\mathbb{L}^{\gamma}(\mathbb{T}^3)\big) \cap\mathbb{L}^{\gamma}(\mathbb{Q}), \quad \gamma>2,\\ \rho\mathbf{u}\in\mathbb{C}\big(\mathbb{S}^1; \mathbb{L}_{\rm weak}^{\frac{2\gamma}{\gamma+1}}(\mathbb{Q})\big). \end{gather*} (2) The fluid velocity $\mathbf{u}$, the angular velocity $\Omega$ and the magnetization $\mathbf{M}$ satisfy \begin{gather*} \mathbf{u}\in\left(\mathbb{L}^2(\mathbb{S}^1;\mathbb{W}^{1,2}(\mathbb{T}^3))\right)^3, \quad \Omega\in\left(\mathbb{L}^2(\mathbb{S}^1;\mathbb{W}^{1,2} (\mathbb{T}^3))\right)^3,\\ \rho\Omega\in\mathbb{C}\big(\mathbb{S}^1;\mathbb{L}_{\rm weak} ^{\frac{2\gamma}{\gamma+1}}(\mathbb{Q})\big), \quad \mathbf{M}\in\left(\mathbb{L}^2(\mathbb{S}^1;\mathcal{M})\right)^3\cap \left(\mathbb{C}(\mathbb{S}^1;\mathbb{L}^2_{\rm weak}(\mathbb{T}^3))\right)^3. \end{gather*} Moreover, $\rho\mathbf{u}\otimes\mathbf{u}$ is integrable on $\mathbb{Q}$. (3) The magnetic field $\mathbf{H}$ is so that $\mathbf{H}=\nabla\psi$ where \[ \psi\in\mathbb{L}^{\infty}(\mathbb{S}^1;\mathbb{W}^{1,2}(\mathbb{T}^3)) \cap\mathbb{L}^2(\mathbb{S}^1;\mathbb{W}^{2,2}(\mathbb{T}^3)), \] solving the problem \begin{equation}\label{E1-m1} -\Delta\psi=\operatorname{div}\mathbf{M}-\mathbb{F}, \end{equation} with $\int_{\mathbb{T}^3}\psi dx=0$ and \begin{gather*} \frac{\partial\psi_i}{\partial x_j}=0\quad\text{on the opposite faces}\\ \{x_i=0,x_j\in[0,\pi],\,j\neq i\}\cup\{x_i=\pi,\,x_j\in[0,\pi],\,j\neq i\}. \end{gather*} (4) For any $(t,x)\in\mathbb{Q}$ and $i=1,2,3$, there holds \begin{gather}\label{E2-4} \int_{\mathbb{T}^3}\rho(t)dx=m>0,\quad \rho(t,Y_i(x))=\rho(t,x),\\ \label{E2-4R1} u^i(t,Y_i(x))=-u^i(t,x),\quad u^i(t,Y^j(x))=u^j(t,x),\quad j\neq i, \\ \label{E2-4R2} \Omega^i(t,Y_i(x))=-\Omega^i(t,x),\quad \Omega^i(t,Y^j(x))=\Omega^j(t,x),\quad j\neq i, \end{gather} where $Y_i$ satisfies \[ Y_i[x_1,\dots,x_i,\dots,x_3]=[x_1,\dots,-x_i,\dots,x_3]. \] (5) The continuity system \eqref{E1-1}-\eqref{E1-3} is satisfied in the sense of renormalized solutions; i.e., the equation \begin{equation}\label{E1-R} \frac{\partial b(\rho)}{\partial t}+\operatorname{div}(b(\rho)\mathbf{u})+(b'(\rho)\rho-b(\rho))\operatorname{div}\mathbf{u}=0 \end{equation} holds in $\mathcal{D}'(\mathbb{R}^1)$ and for any function $b\in\mathbb{C}^1(\mathbb{R}^+)$. (6) The energy inequality \begin{align*} E(t)+C\int_0^tE_f(s)ds &\leq E(0)+C\int_0^t(1+\|\mathbb{F}(s)\|^2+\|\partial_t\mathbb{F}(s)\|^2)ds\\ &\quad +C\int_0^t(1+\|\mathbf{f}(s)\|^2_{\mathbb{L}^{\infty}(\mathbb{Q})} +\|\mathbf{g}(s)\|^2_{\mathbb{L}^{\infty}(\mathbb{Q})} +\|\mathbf{l}(s)\|^2_{\mathbb{L}^{\infty}(\mathbb{Q})})ds, \end{align*} holds for almost everywhere $t\in\mathbb{S}^1$, where \begin{gather*} E(t)=\int_{\mathbb{T}^3}\Big(\frac{1}{2}\rho(|\mathbf{u}|^2+|\Omega|^2) +\frac{a}{\gamma-1}\rho^{\gamma}+\frac{\mu_0}{2}(|\mathbf{M}|^2 +|\mathbf{H}|^2)\Big)dx, \\ \begin{aligned} E_f(t) &=\int_{\mathbb{T}^3}(\mu|\nabla\mathbf{u}|^2+\mu'|\nabla\Omega|^2 +(\mu+\lambda)|\operatorname{div}\mathbf{u}|^2 +(\mu'+\lambda')|\operatorname{div}\Omega|^2)dx\\ &\quad +\int_{\mathbb{T}^3}(\frac{\mu_0}{\tau}|\mathbf{M}|^2 +\frac{\mu_0(2\chi_0+1)}{\tau}|\mathbf{H}|^2 +\mu_0\sigma|\nabla\times\mathbf{M}|^2)dx\\ &\quad +\int_{\mathbb{T}^3}2\mu_0\sigma|\operatorname{div}\mathbf{M}|^2 +\zeta|\nabla\times\mathbf{u}-2\Omega|^2dx\,, \end{aligned}\\ E(0)=\int_{\mathbb{T}^3}\Big(\frac{1}{2}\rho_0(\mathbf{u}_0^2+\Omega^2_0) +\frac{a}{\gamma-1}\rho_0^{\gamma}+\frac{\mu_0}{2} (|\mathbf{M}_0|^2+|\mathbf{H}_0|^2)\Big)dx. \end{gather*} The main result of this paper is stated as follows. \begin{theorem} \label{thm1} Assume that $f^i,g^i,l^i\in\mathbb{L}^{\infty}(\mathbb{Q})$ and satisfy \begin{gather*} f^i(t,Y_i(x))=-f^i(t,x),\quad f^i(t,Y^j(x))=f(t,x),\quad j\neq i,\quad i,j=1,2,3,\\ g^i(t,Y_i(x))=-g^i(t,x),\quad g^i(t,Y^j(x))=g(t,x),\quad j\neq i,\quad i,j=1,2,3, \end{gather*} and $\mathbb{F}\in\mathbb{H}^1(\mathbb{S}^1;\mathbb{L}^2(\mathbb{T}^3))$, $\int_{\mathbb{T}^3}\mathbb{F}(t)dx=0$, the pressure function $P(\rho,\mathbf{M})=P_e(\rho)+P_m(\mathbf{M})$ with $P_e(\rho)=a\rho^{\gamma}$ and $P_m(\mathbf{M})=\frac{\mu_0}{2}|\mathbf{M}|^2$, where $a,\mu_0>0$ and $\gamma>\frac{3}{2}$. Then, for a given $m\geq0$, equation \eqref{E1-1}-\eqref{E1-R5} has a time-periodic weak solution \[ (\rho(t,x),\mathbf{u}(t,x),\Omega(t,x),\mathbf{M}(t,x),\mathbf{H}(t,x)) \quad\text{in the sense of (1)-(6) above}. \] \end{theorem} The proof of this theorem will be carried on by means of a three-level approximation scheme based on a modified system under boundary conditions \eqref{E1-A}: the continuity equation \begin{equation}\label{E2-1} \partial_t\rho+\operatorname{div}(\rho\mathbf{u}) =\epsilon\Delta\rho-2\epsilon\rho+M(\int_{\mathbb{T}^3}\rho dx), \end{equation} the linear momentum equation \begin{equation} \label{E2-2} \begin{aligned} &\partial_t(\rho u^i)+\operatorname{div}(\rho u^i\otimes \mathbf{u})-\mu\Delta u^i-(\mu+\lambda)\partial_{x_i} (\operatorname{div}\mathbf{u})+\partial_{x_i}P(\rho,\mathbf{M}) \\ &=R-\delta\partial_{x_i}\rho^{\beta}-\epsilon\nabla u^i\nabla\rho -2\epsilon\rho u^i+\rho f^i, \end{aligned} \end{equation} the angular momentum equation \begin{equation} \label{E2-R1} \begin{aligned} &\partial_t(\rho\Omega^i)+\operatorname{div}(\rho\Omega^i\otimes\mathbf{u} )-\mu'\Delta\Omega^i-(\lambda'+\mu')\nabla(\operatorname{div}\Omega) \\ &=S-\epsilon\nabla\Omega^i\nabla\rho-2\epsilon\rho\Omega^i+\rho g^i(t,x), \end{aligned} \end{equation} the Bloch-Torrey equation, \begin{equation}\label{E2-3} \partial_t\mathbf{M}+\operatorname{div}(\mathbf{u}\otimes \mathbf{M})-\sigma\Delta\mathbf{M}+\frac{1}{\tau}(\mathbf{M}-\chi_0\mathbf{H}) =\Omega\times\mathbf{M}+\rho\mathbf{l}(t,x), \end{equation} the magnetostatic equation \begin{equation}\label{E2-R2} \mathbf{H}=\nabla\psi,\quad \operatorname{div}(\mathbf{H}+\mathbf{M})=\mathbb{F}, \end{equation} where \begin{gather}\label{E2-R3} R=\mu_0\mathbf{M}\cdot\nabla\mathbf{H}-\zeta\nabla\times (\nabla\times u^i-2\Omega), \\ \label{E2-R4} S=\mu_0\mathbf{M}\times\mathbf{H}+2\zeta(\nabla\times\mathbf{u}-2\Omega^i), \end{gather} $\epsilon, \delta>0$ are small, $\beta>0$ sufficiently large and $M(t)\in\mathbb{C}^{\infty}(\mathbb{R}^1)$, \[ M(t)=\begin{cases} 1, & t\in(-\infty,0],\\ M'(t)<0, & t\in(0,m),\\ 0, & t\in[m,\infty). \end{cases} \] In following sections, we will prove our main result by three steps. Note that the continuity equation \eqref{E2-1} can be solved as done in \cite{Fei1}. Hence, this result is taken from \cite[Lemma 1]{Fei1}. At the first step, by using a modified Faedo-Galerkin approximation and a standard fixed point theorem, a finite-dimensional system is solved and the existence of periodic weak solution (in the sense of (1)-(6)) $(\rho_{n,\epsilon,\delta},\mathbf{u}_{n,\epsilon,\delta}, \Omega_{n,\epsilon,\delta},\mathbf{M}_{n,\epsilon,\delta}, \mathbf{H}_{n,\epsilon,\delta})$ is obtained, i.e., a fixed point of the corresponding Poincar\'e map on a bounded invariant set is found. The rest steps are: let the artificial viscosity terms and artificial pressure term represented by the $\epsilon$-quantities and the $\delta$-quantity go to zero, respectively. To eliminate the $\epsilon$-quantities, we use the weak continuity of the effective viscous flux $P(\rho)-(\lambda+2\mu)\operatorname{div}\mathbf{u}$ which used by Lions \cite{Lions3} and Feireisl \cite{Fei1,Fei3} to remove $\epsilon$-quantities. More specially, it is shown that \[ (P(\rho_n)-(\lambda+2\mu)\operatorname{div}\mathbf{u}_n)b(\rho_n)\to (\overline{P(\rho)}-(\lambda+2\mu)\operatorname{div}\mathbf{u}) \overline{b(\rho)} \] weakly as $n\to 0$, where the over bar and $(\rho_n,\mathbf{u}_n)$ denote the corresponding weak limit and a suitable sequence of approximate solutions, respectively. The process of elimination of the $\delta$-quantity is similar with removing $\epsilon$-quantities, with only difference that the loss of integrability of the density must be compensated for by replacing the function $\rho$ by $\mathcal{T}_k[\rho]$ (defined in \eqref{E5-13}) in \[ \int_{\mathbb{Q}}(P(\rho_n)-(\lambda+2\mu)\operatorname{div} \mathbf{u}_n)b(\rho_n)\,dx\,dt\to \int_{\mathbb{Q}}(\overline{P(\rho)}-(\lambda+2\mu) (\operatorname{div}\mathbf{u})\overline{b(\rho)}\,dx\,dt, \] weakly as $\epsilon\to 0$. \section{Faedo-Galerkin approximation scheme} Let \begin{gather*} \begin{aligned} \mathbb{X}_n &=\Big\{\mathbf{w}=[w^1,w^2,w^3]: w^j=\sum_{|\mathbf{k}|\leq n} a_{\mathbf{k}}[w^j]e^{i\mathbf{k}\cdot x}, \text{ where } a_{Y_i[\mathbf{k}]}[w^i]=-a_{\mathbf{k}}[w^i],\\ &\quad a_{Y_j[\mathbf{k}]}[w^i]=a_{\mathbf{k}}[w^i],\, j\neq i\text{ for all } i=1,2,3\Big\}, \end{aligned}\\ \mathbb{Y}_n=\{\mathbf{M}_n|\mathbf{M}_n(x,t)=\sum_{|k|\leq n} \mathbf{m}_k(t)\varphi_k(x),\, n\in\mathbb{N}\}, \\ \mathbb{Z}_n= \{\mathbf{H}_n|\mathbf{H}_n(x,t) =\sum_{|k|\leq n}\mathbf{h}_k(t)\nabla\varphi_k(x),\, n\in\mathbb{N}\} \end{gather*} be the finite dimensional space with the $\mathbb{L}^2$, $\mathbb{H}_0^1$ and $\mathbb{H}^2$ Hilbert space structure, respectively, where $\{\psi_k(x)\}_{k\in\mathbb{N}}$ is an orthonormal basis of $\mathbb{H}_0^1(\mathbb{T}^3)$ and $\{\nabla\psi_k(x)\}_{k\in\mathbb{N}}$ be an orthonormal basis of $\mathbb{H}^2(\mathbb{T}^3)$, the symbols $a_{\mathbf{k}},h_{\mathbf{k}},\mathbf{k}\in\mathbb{Z}^3$ denote the Fourier coefficients. We define an approximate solution $(\rho_n,\mathbf{u}_n,\Omega_n,\mathbf{M}_n,\mathbf{H}_n)$ of the problem \eqref{E2-1}-\eqref{E2-R2} by the following scheme: The equation of $\rho_n$: \begin{equation}\label{E3-1} \partial_t\rho_n+\operatorname{div}(\rho_n\mathbf{u}_n)=\epsilon\Delta\rho_n-2\epsilon\rho_n+M(\int_{\mathbb{T}^3}\rho_n dx), \end{equation} with the initial data $\rho_n(0)$ satisfying \begin{equation}\label{E3-4} 0<\bar{\rho}(0)\leq\rho_n(0)\leq\tilde{\rho}(0),\quad \rho_n(0)\in\mathbb{C}(\mathbb{T}^3)\cap\mathbb{W}^{1,2}(\mathbb{T}^3),\quad \rho_0(Y_i(x))=\rho_0(x), \end{equation} where $\bar{\rho}(0)$ and $\tilde{\rho}(0)$ are given constants. The equation of $\mathbf{u}_n$: \begin{align} &\frac{d}{dt}\int_{\mathbb{T}^3}\rho_n\mathbf{u}_n\cdot\mathbf{w}dx \nonumber\\ &= \int_{\mathbb{T}^3}\left(\rho_n(\mathbf{u}_n\otimes \mathbf{u}_n)\cdot\nabla\mathbf{w}-\mu\nabla \mathbf{u}_n\cdot\nabla\mathbf{w}-(\mu+\lambda)\operatorname{div} \mathbf{u}_n\cdot\operatorname{div}\mathbf{w}\right)dx \nonumber\\ &\quad +\int_{\mathbb{T}^3}\left(a(\rho_n)^{\gamma}\operatorname{div}\mathbf{w}+\delta(\rho_n)^{\beta}\operatorname{div}\mathbf{w} +\frac{\mu_0}{2}|\mathbf{M}_n|^2\operatorname{div}\mathbf{w} -\epsilon\nabla\mathbf{u}_n\cdot\nabla\rho_n\cdot\mathbf{w}\right)dx \nonumber\\ &\quad +\int_{\mathbb{T}^3}\left(-2\epsilon\rho_n\mathbf{u}_n \cdot\mathbf{w}+\mu_0\mathbf{M}_n\cdot\nabla\mathbf{H}_n\cdot\mathbf{w} +2\zeta(\nabla\times\Omega_n)\cdot\mathbf{w}\right)dx \nonumber\\ &\quad +\int_{\mathbb{T}^3}\left(-\zeta(\nabla\times\mathbf{u}_n) \cdot(\nabla\times\mathbf{w}) +\rho_n\mathbf{f}_n\cdot \mathbf{w}\right)dx, \label{E3-2} \end{align} with the initial conditions $\mathbf{u}_n(0)\in\mathbb{X}_n$. The equation of $\Omega_n$: \begin{equation}\label{E3-2R} \begin{aligned} &\frac{d}{dt}\int_{\mathbb{T}^3}\rho_n\Omega_n\cdot\mathbf{w}dx \\ &= \int_{\mathbb{T}^3}\left(\rho_n(\mathbf{u}_n\otimes \Omega_n)\cdot\nabla\mathbf{w}-\mu'\nabla \Omega_n\cdot\nabla\mathbf{w}-(\mu'+\lambda')\operatorname{div} \Omega_n\cdot\operatorname{div}\mathbf{w}\right)dx \\ &\quad +\int_{\mathbb{T}^3}\left(-2\epsilon\rho_n\Omega_n\cdot \mathbf{w}-\epsilon\nabla\Omega_n\cdot\nabla\rho_n\cdot\mathbf{w} +\mu_0(\mathbf{M}_n\times\mathbf{H}_n)\cdot\mathbf{w}\right)dx \\ &\quad +\int_{\mathbb{T}^3}\left(2\zeta(\nabla\times\mathbf{u}_n -2\Omega_n)\cdot\mathbf{w}+\rho_n\mathbf{g}_n\cdot \mathbf{w}\right)dx, \end{aligned} \end{equation} with the initial conditions $\Omega_n(0)\in\mathbb{X}_n$. The equation of $\mathbf{M}_n$: \begin{equation}\label{E3-3} \partial_t\mathbf{M}_n+\operatorname{div}(\mathbf{u}_n\otimes \mathbf{M}_n)-\sigma\Delta\mathbf{M}_n+\frac{1}{\tau}(\mathbf{M}_n -\chi_0\mathbf{H}_n)=\Omega_n\times\mathbf{M}_n+\rho_n\mathbf{l}_n(t,x), \end{equation} with the initial conditions $\mathbf{M}_n(0)\in\mathbb{Y}_n$. The equation of $\mathbf{H}_n$: \begin{equation}\label{E3-3R} \mathbf{H}_n=\nabla\psi_n,\quad \operatorname{div}(\mathbf{H}_n+\mathbf{M}_n) =\mathbb{F}_n, \end{equation} with the initial conditions $\mathbf{H}_n(0)\in\mathbb{Z}_n$. Moreover, for any fixed constant $n$, $\mathbf{f}_n,\mathbf{g}_n,\mathbf{l}_n\in\mathbb{C}^{\infty}(\mathbb{Q})$, $\mathbf{f}_n,\mathbf{g}_n$ satisfies \eqref{E2-4R1} and \begin{equation}\label{ER} \begin{gathered} \|f_n^i\|_{\mathbb{L}^{\infty}(\mathbb{Q})} \leq\|f^i\|_{\mathbb{L}^{\infty}(\mathbb{Q})}, \quad \|g_n^i\|_{\mathbb{L}^{\infty}(\mathbb{Q})} \leq\|g^i\|_{\mathbb{L}^{\infty}(\mathbb{Q})},\quad \|\mathbf{l}_n\|_{\mathbb{L}^{\infty}(\mathbb{Q})} \leq\|\mathbf{l}\|_{\mathbb{L}^{\infty}(\mathbb{Q})},\\ f_n^i\to f^i,\quad g_n^i\to g^i,\quad \mathbf{l}_n\to \mathbf{l}\quad \text{strongly in }\mathbb{L}^2(\mathbb{Q}),\quad i=1,2,3. \end{gathered} \end{equation} Now, we have the following result: \begin{proposition} \label{prop1} Assume that $\beta>4$. Then there exits a time periodic solution $(\rho,\mathbf{u},\Omega,\mathbf{M},\mathbf{H})$ of the problem \begin{gather}\label{E3-R1} \partial_t\rho+\operatorname{div}(\rho\mathbf{u}) =\epsilon\Delta\rho-2\epsilon\rho+M(\int_{\mathbb{T}^3}\rho dx),\quad \text{a.e. on }\mathbb{Q},\\ \label{E3-R2} \begin{aligned} &\partial_t(\rho u^i)+\operatorname{div}(\rho u^i\otimes \mathbf{u})-\mu\Delta u^i-(\mu+\lambda)\partial_{x_i} (\operatorname{div}\mathbf{u})+\partial_{x_i}P(\rho,\mathbf{M}) \\ &=R-\delta\partial_{x_i}\rho^{\beta}-\epsilon\nabla u^i\nabla\rho -2\epsilon\rho u^i+\rho f^i,\quad i=1,2,3,\text{ in } \mathcal{D}'(\mathbb{Q}), \end{aligned} \\ \label{E6-4} \begin{aligned} &\partial_t(\rho\Omega^i)+\operatorname{div}(\rho\Omega^i\otimes\mathbf{u} )-\mu'\Delta\Omega^i-(\lambda'+\mu')\nabla(\operatorname{div}\Omega) \\ &=S-\epsilon\nabla\Omega^i\nabla\rho-2\epsilon\rho\Omega^i+\rho g^i(t,x), \quad\text{in }\mathcal{D}'(\mathbb{Q}), \end{aligned} \\ \label{E3-R3} \partial_t\mathbf{M}+\operatorname{div}(\mathbf{u}\otimes \mathbf{M})-\sigma\Delta\mathbf{M}+\frac{1}{\tau}(\mathbf{M}-\chi_0\mathbf{H}) =\Omega\times\mathbf{M}+\rho\mathbf{l}(t,x),~in~\mathcal{D}'(\mathbb{Q}),\\ \label{E6-5} \mathbf{H}=\nabla\psi,\quad \operatorname{div}(\mathbf{H}+\mathbf{M}) =\mathbb{F},\quad\text{in }\mathcal{D}'(\mathbb{Q}), \end{gather} where \begin{gather*} R=\mu_0\mathbf{M}\cdot\nabla\mathbf{H}-\zeta\nabla \times(\nabla\times u^i-2\Omega), \\ S=\mu_0\mathbf{M}\times\mathbf{H}+2\zeta(\nabla\times\mathbf{u}-2\Omega^i), \end{gather*} function $\rho\geq0$ belongs to the class corresponding to \eqref{E3-28}-\eqref{E3-29} and satisfies \eqref{E2-4} and \begin{equation}\label{E5-14} \int_{\mathbb{T}^3}\rho(t)dx=m_{\epsilon}, \quad \forall t\in\mathbb{S}^1, \quad 2\epsilon m_{\epsilon}=|\mathbb{T}^3|M(m_{\epsilon}). \end{equation} The fluid velocity and angular velocity $\mathbf{u},\Omega\in[\mathbb{L}^2(\mathbb{S}^1; \mathbb{W}^{1,2}(\mathbb{T}^3))]^3$ satisfies \eqref{E2-4R1} and \eqref{E2-4R2} a.e. on $\mathbb{Q}$. The Magnetization $\mathbf{M}\in\mathbb{C}(\mathbb{S}^1;\mathbb{L}^2_{\rm weak} (\mathbb{T}^3))\cap\mathbb{L}^2(\mathbb{S}^1;\mathbb{H}^1(\mathbb{T}^3))$ and $\mathbf{H}=\nabla\psi$, $\psi\in\mathbb{L}^2(\mathbb{S}^1; \mathbb{H}^2(\mathbb{T}^3))\cap\mathbb{L}^{\infty}(\mathbb{S}^1; \mathbb{H}^1(\mathbb{T}^3))$. The energy $E_{\delta}(t)\in\mathbb{L}^{\infty}(\mathbb{S}^1)$ such that \begin{equation}\label{E3-30} \begin{aligned} &E_n(t)+C\int_0^tE_{nf}(s)ds\\ &\leq E(0)+C\int_0^t(1+\|\mathbb{F}(s)\|^2+\|\partial_t\mathbb{F}(s)\|^2)ds \\ &\quad +C\int_0^t(1+\|\mathbf{f}(s)\|^2_{\mathbb{L}^{\infty}(\mathbb{Q})} +\|\mathbf{g}(s)\|^2_{\mathbb{L}^{\infty}(\mathbb{Q})} +\|\mathbf{l}(s)\|^2_{\mathbb{L}^{\infty}(\mathbb{Q})})ds \end{aligned} \end{equation} holds on $\mathbb{S}^1$, where \begin{align*} E_n(t)=\int_{\mathbb{T}^3}\Big(\frac{1}{2}\rho(|\mathbf{u}|_n^2 +|\Omega|_n^2)+\frac{a}{\gamma-1}\rho_n^{\gamma} +\frac{\mu_0}{2}(|\mathbf{M}_n|^2+|\mathbf{H}_n|^2)\Big)dx, \\ \begin{aligned} E_{nf}(t) &=\int_{\mathbb{T}^3}(\mu|\nabla\mathbf{u}_n|^2 +\mu'|\nabla\Omega_n|^2+(\mu+\lambda)|\operatorname{div}\mathbf{u}_n|^2 +(\mu'+\lambda')|\operatorname{div}\Omega_n|^2)dx\\ &\quad +\int_{\mathbb{T}^3}(\frac{\mu_0}{\tau}|\mathbf{M}_n|^2 +\frac{\mu_0(2\chi_0+1)}{\tau}|\mathbf{H}_n|^2 +\mu_0\sigma|\nabla\times\mathbf{M}_n|^2)dx\\ &\quad +\int_{\mathbb{T}^3}2\mu_0\sigma|\operatorname{div}\mathbf{M}_n|^2 +\zeta|\nabla\times\mathbf{u}_n-2\Omega_n|^2dx\,, \end{aligned} \\ E(0)=\int_{\mathbb{T}^3}\left(\frac{1}{2}\rho_0(\mathbf{u}_0^2+\Omega^2_0) +\frac{a}{\gamma-1}\rho_0^{\gamma}+\frac{\mu_0}{2}(|\mathbf{M}_0|^2 +|\mathbf{H}_0|^2)\right)dx, \end{align*} where $C$ is independent of $\epsilon$ and $\delta$. \end{proposition} The following lemma is taken from \cite{Fei1} which shows that the local solvability of \eqref{E3-1}-\eqref{E3-4}. \begin{lemma} \label{lem1} Assume that the initial data $\rho_n(0)$ satisfies \eqref{E3-4} and fix $\mathbf{u}_n$ in the space $\mathbb{C}([0,T],\mathbb{X}_n)$. Then \eqref{E3-1}-\eqref{E3-4} has a unique solution $\rho_n(t)$ on the interval $(0,\omega)$ satisfying \begin{itemize} \item[(1)] $\rho_n\in\mathbb{C}\left([0,T]; \mathbb{C}(\mathbb{T}^3)\cap\mathbb{W}^{1,2}(\mathbb{T}^3)\right)$ is a classical solution of \eqref{E3-1}-\eqref{E3-4}. \item[(2)] For any $t\in[0,\omega]$, \[ 0<\bar{\rho}(t)\leq\rho_n(t)\leq\tilde{\rho}(t), \] where \begin{gather*} \bar{\rho}(t)=\bar{\rho}(0)\exp(-\epsilon t-\int_{0}^t\|\operatorname{div} \mathbf{u}_n\|_{\mathbb{L}^{\infty}(\mathbb{T}^3)}ds),\\ \tilde{\rho}(t)=\tilde{\rho}(0)\exp(-\epsilon t+\int_{0}^t\|\operatorname{div} \mathbf{u}_n\|_{\mathbb{L}^{\infty}(\mathbb{T}^3)}ds)+t. \end{gather*} \end{itemize} \end{lemma} This lemma tells us that there exists a locally Lipschitz continuous function $\mathbf{u}_n\to \rho_n[\mathbf{u}_n]$ acting between the spaces $\mathbb{C}([0,T];\mathbb{X}_n)$ and $\mathbb{C}([0,T];\mathbb{W}^{1,2}(\mathbb{T}^3)])$. Before showing the solvability of \eqref{E3-2}-\eqref{E3-3R}, we first derive the energy equation and some a priori estimates. Let $\mathbf{w}=\mathbf{u}_n(t)$ in \eqref{E3-2}. We obtain \begin{equation}\label{E3-13} \begin{aligned} &\frac{d}{dt}\int_{\mathbb{T}^3}\Big(\frac{1}{2}\rho_n|\mathbf{u}_n|^2 +\frac{a}{\gamma-1}\rho_n^{\gamma}+\frac{\delta}{\beta-1}\rho_n^{\beta}\Big)dx \\ &+\int_{\mathbb{T}^3}\Big(\epsilon\rho_n|\mathbf{u}_n|^2+\mu|\nabla\mathbf{u}_n|^2 +(\mu+\lambda)|\operatorname{div}\mathbf{u}_n|^2 +\zeta|\nabla\times\mathbf{u}_n|^2+\frac{2a\epsilon\gamma}{\gamma-1} \rho_n^{\gamma}\\ &\quad +\frac{2\delta\epsilon\beta}{\beta-1}\rho_n^{\beta}\Big)dx \\ &\leq \int_{\mathbb{T}^3}|\mathbf{f}_n||\mathbf{u}_n|\rho_ndx +\int_{\mathbb{T}^3}M(\int_{\mathbb{T}^3}\rho_ndx) \Big(\frac{a\gamma}{\gamma-1}\rho_n^{\gamma-1}+\frac{\delta\beta}{\beta-1} \rho_n^{\beta-1}\Big)dx \\ &\quad +\int_{\mathbb{T}^3}\Big(\frac{\mu_0}{2}|\mathbf{M}_n|^2\operatorname{div}\mathbf{u}_n+\mu_0\mathbf{M}_n\cdot\nabla\mathbf{H}_n\cdot\mathbf{u}_n +2\zeta(\nabla\times\Omega_n)\cdot\mathbf{u}_n\Big)dx. \end{aligned} \end{equation} Taking $\mathbf{w}=\Omega_n(t)$ in \eqref{E3-2R} yields \begin{equation}\label{E3-13R} \begin{aligned} &\frac{d}{dt}\int_{\mathbb{T}^3}\frac{1}{2}\rho_n|\Omega_n|^2dx + \int_{\mathbb{T}^3}\epsilon\rho_n|\Omega_n|^2+\mu'|\nabla\Omega_n|^2 +(\mu'+\lambda')|\operatorname{div}\Omega_n|^2)dx\\ &+\int_{\mathbb{T}^3}\zeta|\nabla\times\Omega_n|^2dx+4\zeta|\Omega_n|^2dx \\ &\leq \int_{\mathbb{T}^3}|\mathbf{g}_n||\Omega_n|\rho_n +\mu_0(\mathbf{M}_n\times\mathbf{H}_n)\cdot\Omega_n +2\zeta(\nabla\times\Omega_n)\cdot\Omega_ndx. \end{aligned} \end{equation} Multiplying \eqref{E3-3} by $\mathbf{M}_n$, integrating on $\mathbb{T}^3$ and using the equation \eqref{E3-3R}, we have \begin{equation}\label{E3-13RR} \begin{aligned} &\frac{1}{2}\frac{d}{dt}\int_{\mathbb{T}^3}|\mathbf{M}_n|^2dx +\int_{\mathbb{T}^3}\frac{1}{\tau}(|\mathbf{M}_n|^2+\chi_0|\mathbf{H}_n|^2) +\sigma(|\nabla\mathbf{M}_n|^2+|\operatorname{div}\mathbf{M}_n|^2)dx \\ &=-\int_{\mathbb{T}^3}\frac{\mu_0}{2}|\mathbf{M}_n|^2\operatorname{div} \mathbf{u}_n-\frac{\chi_0}{\tau}\mathbb{F}_n\psi_n +\rho_n\mathbf{l}_n\mathbf{M}_ndx. \end{aligned} \end{equation} Multiplying \eqref{E3-3} by $\mathbf{H}_n$, and integrating on $\mathbb{T}^3$, we have \begin{equation}\label{E3-14} \begin{aligned} &\int_{\mathbb{T}^3}\partial_t\mathbf{M}_n\cdot\mathbf{H}_ndx -\sigma\int_{\mathbb{T}^3}\Delta\mathbf{M}_n\cdot\mathbf{H}_ndx +\frac{1}{\tau}\int_{\mathbb{T}^3}\mathbf{M}_n\cdot\mathbf{H}_ndx \\ &= \int_{\mathbb{T}^3}\left(\operatorname{div}(\mathbf{u}_n\otimes \mathbf{M}_n)\cdot\mathbf{H}_n+(\Omega_n\times\mathbf{M}_n) \cdot\mathbf{H}_n+\frac{\chi_0}{\tau}|\mathbf{H}_n|^2 +|\mathbf{l}_n\mathbf{H}_n|\rho_n\right)dx. \end{aligned} \end{equation} First using the relation $\nabla\times\mathbf{H}_n=0$, we have \begin{equation}\label{E3-14'} \int_{\mathbb{T}^3}\operatorname{div}(\mathbf{u}_n\otimes\mathbf{M}_n) \cdot\mathbf{H}_ndx =-\int_{\mathbb{T}^3}\mathbf{M}_n\cdot\nabla\mathbf{H}_n\cdot\mathbf{u}_n dx =\int_{\mathbb{T}^3}\mathbf{u}_n\cdot\nabla\mathbf{H}_n\cdot\mathbf{M}_ndx. \end{equation} Multiplying \eqref{E3-3R} by $\psi_n$, then differentiating, we obtain \begin{equation}\label{E3-14''} \int_{\mathbb{T}^3}\partial_t\mathbf{M}_n\cdot\mathbf{H}_ndx =-\frac{1}{2}\frac{d}{dt}\int_{\mathbb{T}^3}\|\mathbf{H}_n\|^2dx -\int_{\mathbb{T}^3}\partial_t\mathbb{F}\psi_n dx. \end{equation} Multiplying the identity \[ \nabla\times(\nabla\times\mathbf{M}_n)=\nabla \operatorname{div}\mathbf{M}_n-\Delta\mathbf{M}_n, \] by $\mathbf{H}_n$, using the boundary condition $(\nabla\times\mathbf{M})\times\nu(x)=0$ and \eqref{E3-3R}, we have \begin{equation}\label{E3-14R} -\int_{\mathbb{T}^3}\Delta\mathbf{M}_n\cdot\mathbf{H}_ndx =-\int_{\mathbb{T}^3}|\operatorname{div}\mathbf{M}_n|^2dx +\int_{\mathbb{T}^3}\mathbb{F}_n\operatorname{div}\mathbf{M}_ndx \end{equation} Substituting \eqref{E3-14'}-\eqref{E3-14R} into \eqref{E3-14}, we obtain \begin{equation}\label{E3-14RR} \begin{aligned} &\frac{1}{2}\frac{d}{dt}\int_{\mathbb{T}^3}|\mathbf{H}_n|^2dx +\sigma\int_{\mathbb{T}^3}|\operatorname{div}\mathbf{M}_n|^2dx\\ &=-\int_{\mathbb{T}^3}\partial_t\mathbb{F}\psi_n +\sigma\mathbb{F}_n\operatorname{div}\mathbf{M}_ndx -\int_{\mathbb{T}^3}\mathbf{u}_n\cdot\nabla\mathbf{H}_n\cdot\mathbf{M}_n -(\Omega_n\times\mathbf{M}_n)\cdot\mathbf{H}_ndx \\ &\quad -\int_{\mathbb{T}^3}\frac{\chi_0}{\tau}|\mathbf{H}_n|^2 -|\mathbf{l}_n\mathbf{H}_n|\rho_ndx. \end{aligned} \end{equation} Hence, summing up \eqref{E3-13}, \eqref{E3-13R}, \eqref{E3-13RR} and \eqref{E3-14RR}, we obtain the energy inequality \begin{equation}\label{E3-15} \begin{aligned} \frac{d}{dt}E_n(t)+E_f(t) &\leq\int_{\mathbb{T}^3}(|\mathbf{f}_n||\mathbf{u}_n|+|\mathbf{g}_n||\Omega_n|+|\mathbf{l}_n|(|\mathbf{M}_n|+|\mathbf{H}_n|))\rho_ndx \\ &\quad +\int_{\mathbb{T}^3}M(\int_{\mathbb{T}^3}\rho_ndx) \Big(\frac{a\gamma}{\gamma-1}\rho_n^{\gamma-1} +\frac{\delta\beta}{\beta-1}\rho_n^{\beta-1}\Big)dx \\ &\quad +C\int_{\mathbb{T}^3}(1+|\mathbb{F}|^2+|\partial_t\mathbb{F}|^2)dx, \end{aligned} \end{equation} where \begin{gather}\label{E3-16} E_n(t)=\int_{\mathbb{T}^3}\Big(\frac{1}{2}\rho_n(|\mathbf{u}_n|^2 +|\Omega_n|^2)+\frac{a}{\gamma-1}\rho_n^{\gamma}+\frac{\delta}{\beta-1} \rho_n^{\beta}+\frac{\mu_0}{2}(|\mathbf{M}_n|^2+|\mathbf{H}_n|^2)\Big)dx,\\ \label{E3-16'} \begin{aligned} E_f(t) &=\int_{\mathbb{T}^3}(\mu|\nabla\mathbf{u}_n|^2+\mu'|\nabla\Omega_n|^2 +(\mu+\lambda)|\operatorname{div}\mathbf{u}_n|^2+(\mu'+\lambda') |\operatorname{div}\Omega_n|^2)dx \\ &\quad +\int_{\mathbb{T}^3}(\frac{\mu_0}{\tau}|\mathbf{M}_n|^2 +\frac{\mu_0(2\chi_0+1)}{\tau}|\mathbf{H}_n|^2 +\mu_0\sigma|\nabla\times\mathbf{M}_n|^2)dx \\ &\quad +\int_{\mathbb{T}^3}\frac{2\epsilon a\gamma}{\gamma-1}\rho_n^{\gamma} +\frac{2\epsilon \delta\beta}{\beta-1}\rho_n^{\beta} +2\mu_0\sigma|\operatorname{div}\mathbf{M}_n|^2 +\zeta|\nabla\times\mathbf{u}_n-2\Omega_n|^2dx. \end{aligned} \end{gather} In the following result, we show that if $(\mathbf{u}_n,\Omega_n)$ is given, then there exists $\omega$-periodic solution $(\mathbf{M}_n,\mathbf{H}_n)$ for \eqref{E3-3} and \eqref{E3-3R}. \begin{lemma} \label{lem2} Assume that the initial data $(\mathbf{M}_n(0),\mathbf{H}_n(0))\in\mathbb{Y}_n\times\mathbb{Z}_n$ and fix a velocity field $\mathbf{u}_n,\Omega_n\in\mathbb{C}([0,T],\mathbb{X}_n)$. Then \eqref{E3-3} and \eqref{E3-3R} have solutions \[ (\mathbf{M}_n(t,x),\mathbf{H}_n(t,x))\in\mathbb{C}^1([0,\omega]; \mathbb{H}^1_0(\mathbb{T}^3))\times\mathbb{C}^1([0,\omega]; \mathbb{H}^1_0(\mathbb{T}^3))\cap\mathbb{C}^1([0,\omega]; \mathbb{H}^2(\mathbb{T}^3)) \] such that $(\mathbf{M}_n(0,x),\mathbf{H}_n(0,x)) =(\mathbf{M}_n(\omega,x),\mathbf{H}_n(\omega,x))$. Moreover, the operators $(\mathbf{u}_n,\Omega_n)\to \mathbf{M}_n(\mathbf{u}_n,\Omega_n)$ and $(\mathbf{u}_n,\Omega_n)\to \mathbf{H}_n(\mathbf{u}_n,\Omega_n)$ map bounded sets from $\mathbb{C}([0,T],\mathbb{X}_n\times\mathbb{X}_n)$ into bounded subsets of $\mathbb{Y}_n$ and $\mathbb{Z}_n$, and the solution operators are continuous operators. \end{lemma} \begin{proof} Multiplying \eqref{E3-3} by $\nabla\mathbf{h}$, integrating on $\mathbb{T}^3$, we have \begin{equation}\label{E4-A} \begin{aligned} &(\partial_t\mathbf{M}_n,\nabla\mathbf{h}) -\sigma(\Delta\mathbf{M}_n,\nabla\mathbf{h}) +\frac{1}{\tau}(\mathbf{M}_n,\nabla\mathbf{h})\\ &=(\operatorname{div}(\mathbf{u}_n\otimes \mathbf{M}_n),\nabla\mathbf{h}) +((\Omega_n\times\mathbf{M}_n),\nabla\mathbf{h}) +\frac{\chi_0}{\tau}(\mathbf{H}_n,\nabla\mathbf{h}) +(\rho_n\mathbf{g}_n,\nabla\mathbf{h}). \end{aligned} \end{equation} Now using the relation $\nabla\times(\nabla\mathbf{h})=0$, we obtain \begin{equation}\label{E4-A1} (\operatorname{div}(\mathbf{u}_n\otimes\mathbf{M}_n),\nabla\mathbf{h}) =(\mathbf{u}_n\cdot\mathbf{M}_n,\operatorname{div}(\nabla\mathbf{h})). \end{equation} Multiplying \eqref{E3-3R} by $\mathbf{h}_n$, then differentiating, we obtain \begin{equation}\label{E4-A2} (\partial_t\mathbf{M}_n,\nabla\mathbf{h}) =-\frac{d}{dt}(\mathbf{H}_n,\nabla\mathbf{h}) -(\partial_t\mathbb{F},\nabla\mathbf{h}). \end{equation} Substituting \eqref{E4-A1}-\eqref{E4-A2} into \eqref{E4-A}, we obtain \begin{equation}\label{E4-A4} \begin{aligned} &-\frac{d}{dt}(\mathbf{H}_n,\nabla\mathbf{h}) -\sigma(\Delta\mathbf{M}_n,\nabla\mathbf{h}) +\frac{1}{\tau}(\mathbf{M}_n,\nabla\mathbf{h})\\ &=(\partial_t\mathbb{F},\nabla\mathbf{h}) +(\mathbf{u}_n\cdot\mathbf{M}_n,\operatorname{div}(\nabla\mathbf{h})) +((\Omega_n\times\mathbf{M}_n),\nabla\mathbf{h})\\ &\quad +\frac{\chi_0}{\tau}(\mathbf{H}_n,\nabla\mathbf{h}) +(\rho_n\mathbf{g}_n,\nabla\mathbf{h}). \end{aligned} \end{equation} Multiplying \eqref{E3-3} by $\mathbf{h}_n$, integrating on $\mathbb{T}^3$ and using the equation \eqref{E3-3R}, we have \begin{equation}\label{E4-A5} \begin{aligned} &\frac{d}{dt}(\mathbf{M}_n,\mathbf{h}) -(\mathbf{u}_n\otimes\mathbf{M}_n,\nabla\mathbf{h}) +\sigma(\nabla\mathbf{M}_n,\nabla\mathbf{h})\\ &=(\Omega_n\times\mathbf{M}_n,\mathbf{h}) -\frac{1}{\tau}(\mathbf{M}_n-\chi_0\mathbf{H}_n,\mathbf{h}) +(\rho_n\mathbf{l}_n,\mathbf{h}). \end{aligned} \end{equation} Take $\mathbf{h}(x)=\varphi_k(x)$ from an orthonormal basis of $\mathbb{Y}_n$ and $\mathbb{Z}_n$ in \eqref{E4-A4}-\eqref{E4-A5}, respectively. Note that $\mathbf{M}_n(t,x)\in\mathbb{Y}_n$ and $\mathbf{H}_n(t,x)\in\mathbb{Z}_n$. We can write \[ \mathbf{M}_n(x,t)=\sum_{|k|\leq n}\mathbf{m}_k(t)\varphi_k(x),\quad \mathbf{H}_n(x,t)=\sum_{|k|\leq n}\mathbf{h}_k(t)\nabla\varphi_k(x), \] where the coefficients $\mathbf{m}_k(t)$ and $\mathbf{h}_k(t)$ are required to solve the coupled ordinary differential equations \begin{gather}\label{E3-6R} \frac{d\mathbf{m}_{k}}{dt}+\sum_{|k|\leq n}A^1_{j,k}(t)\mathbf{m}_{k} +\sum_{|k|\leq n}A^2_{j,k}(t)\mathbf{h}_k=B^1_k(t),\quad |k|\leq n,\\ \label{E3-6} \frac{d\mathbf{h}_{k}}{dt}+\sum_{|k|\leq n}A^3_{j,k}(t)\mathbf{h}_{k} +\sum_{|k|\leq n}A^4_{j,k}(t)\mathbf{m}_{k}=B^2_k(t),\quad |k|\leq n, \end{gather} where \begin{gather*} A^1_{j,k}(t)=-(\mathbf{u}_n\otimes\psi_j,\nabla\varphi_k) +\sigma(\nabla\varphi_j,\nabla\varphi_k)-(\Omega_n\times\varphi_j,\varphi_k) +\frac{1}{\tau}(\varphi_j,\varphi_k), \\ A^2_{j,k}(t)=\frac{\chi_0}{\tau}(\nabla\varphi_k,\varphi_j), \quad A^2_{j,k}(t)=\frac{\chi_0}{\tau}(\nabla\varphi_j,\nabla\varphi_k), \\ A^4_{j,k}(t)=\sigma(\Delta\varphi_j,\nabla\varphi_k) -\frac{1}{\tau}(\varphi_j,\nabla\varphi_k) +(\mathbf{u}_n\cdot\varphi_j,\operatorname{div}\nabla\varphi_k) +(\Omega_n\times\varphi_j,\nabla\varphi_k), \\ B^1_k(t)=(l_n\rho_n,\varphi_{k}),\quad B^2_k(t)=(l_n\rho_n,\nabla\varphi_{k}). \end{gather*} Following \cite{Prouse}, for given initial data $\mathbf{M}_n(0)\in\mathbb{Y}_n$ and $\mathbf{H}_n(0)\in\mathbb{Z}_n$, the system \eqref{E3-6R}-\eqref{E3-6} has a unique solution $(\mathbf{m}_k,\mathbf{h}_k)\in\mathbb{C}^1((0,T); \mathbb{Y}_n)\times\mathbb{C}^1((0,T);\mathbb{Z}_n)$ for some $T\leq\omega$. Multiplying both sides \eqref{E3-6R}-\eqref{E3-6} by $\mathbf{m}_k$ and $\mathbf{h}_k$, summing over $k$, integrating by parts, we have \begin{equation}\label{E3-8} \begin{aligned} &\frac{1}{2}\frac{d}{dt}(\|\mathbf{M}_n\|^2+\|\mathbf{H}_n\|^2) +\frac{1}{\tau}\|\mathbf{M}_n\|^2+(\chi_0 +\frac{\chi_0}{\tau})\|\mathbf{H}_n\|^2 +\sigma\|\nabla\mathbf{M}_n\|^2\\ &+2\sigma\|\operatorname{div}\mathbf{M}_n\|^2 \\ &=-\frac{\mu_0}{2}(\mathbf{M}_n\cdot\operatorname{div}\mathbf{u}_n,\mathbf{M}_n) -\frac{\chi_0}{\tau}(\mathbb{F}_n,\psi_n)-(\partial_t\mathbb{F}_n,\psi_n) +\sigma(\mathbb{F}_n,\operatorname{div}\mathbf{M}_n)\\ &\quad -(\mathbf{u}_n\cdot\mathbf{M}_n,\nabla\mathbf{H}_n) +(\Omega_n\times\mathbf{M}_n,\mathbf{H}_n) +(\rho_n\mathbf{l}_n,\mathbf{M}_n+\mathbf{H}_n). \end{aligned} \end{equation} Note that by the Young inequality and $\operatorname{div}\mathbf{H}_n=\mathbb{F}_n-\operatorname{div}\mathbf{M}_n$, \begin{gather}\label{E3-9R} -\frac{\mu_0}{2}(\mathbf{M}_n\cdot\operatorname{div}\mathbf{u}_n,\mathbf{M}_n) \leq \frac{\mu_0}{2}\|\operatorname{div}\mathbf{u}_n\|\|\mathbf{M}_n\|^2,\\ \label{E3-9R1} -\frac{\chi_0}{\tau}(\mathbf{F}_n,\psi_n)-(\partial_t\mathbf{F}_n,\psi_n) \leq \frac{\chi_0}{2\tau}\|\mathbb{F}_n\|^2+\frac{2\tau}{\chi_0} \|\partial_t\mathbb{F}_n\|^2+\frac{\chi_0}{\tau}\|\psi_n\|^2,\\ \label{E3-9R2} \sigma(\mathbf{F}_n,\operatorname{div}\mathbf{M}_n) \leq \frac{\sigma}{2}\|\operatorname{div}\mathbf{M}_n\|^2 +\frac{2}{\sigma}\|\mathbb{F}_n\|^2,\\ \label{E3-9R3} -(\mathbf{u}_n\cdot\mathbf{M}_n,\nabla\mathbf{H}_n) \leq \frac{\sigma}{2}\|\operatorname{div}\mathbf{M}_n\|^2 +\frac{\sigma}{2}\|\mathbb{F}_n\|^2+\frac{2}{\sigma}\|\mathbf{u}_n\|^2 \|\mathbf{M}_n\|^2,\\ \label{E3-9} (\Omega_n\times\mathbf{M}_n,\mathbf{H}_n) \leq \frac{\chi_0}{4}\|\mathbf{H}_n\|^2+\frac{4}{\chi_0}\|\Omega_n\|^2 \|\mathbf{M}_n\|^2,\\ \label{E3-10} (\rho_n\mathbf{l}_n,\mathbf{M}_n+\mathbf{H}_n) \leq \frac{\chi_0}{2}\|\mathbf{H}_n\|^2+\frac{1}{2\tau}\|\mathbf{M}_n\|^2 +(\frac{2}{\chi_0}+2\tau)\|\rho_n\|^2\|\mathbf{l}_n\|^2. \end{gather} Using the Poincar\'{e} inequality and $\mathbf{H}_n=\nabla\psi_n$, \begin{gather}\label{E3-11} \|\nabla\mathbf{M}_n\|_{\mathbb{L}^2(\mathbb{T}^2)} \geq C\|\mathbf{M}_n\|_{\mathbb{L}^2(\mathbb{T}^2)},\\ \label{E3-11'} \|\mathbf{M}_n\|= \|\nabla\psi_n\|\geq C\|\psi_n\|. \end{gather} From \eqref{E3-8}-\eqref{E3-11'}, we derive \begin{align*} &\frac{d}{dt}(\|\mathbf{M}_n(t)\|^2+\|\mathbf{H}_n(t)\|^2) + C(\sigma,\tau,\chi_0)(\|\mathbf{M}_n\|^2+\|\mathbf{H}_n\|^2) \\ &\leq C(\tau,\chi_0)(\|\mathbb{F}_n\|^2+\|\partial_t\mathbb{F}_n\|^2) +(\frac{2}{\chi_0}+2\tau)\|\rho_n\|^2\|\mathbf{l}_n\|^2. \end{align*} Thus, for $t\in[0,T]$, we obtain \begin{equation}\label{E3-12} \begin{aligned} &\|\mathbf{M}_n(t)\|^2+\|\mathbf{H}_n(t)\|^2 \\ &\leq (\|\mathbf{M}_n(0)\|^2+\|\mathbf{H}_n(0)\|^2)e^{-C(\sigma,\tau,\chi_0)t} \\ &\quad+\int_0^te^{-C\nu(t-s)}(C(\tau,\chi_0)(\|\mathbb{F}_n\|^2 +\|\partial_t\mathbb{F}_n\|^2) +(\frac{2}{\chi_0}+2\tau)\|\rho_n\|^2\|\mathbf{l}_n\|^2)ds. \end{aligned} \end{equation} This implies that for $t\in[0,T]$, \begin{align*} &\|\mathbf{M}_n(t)\|^2+\|\mathbf{H}_n(t)\|^2\\ &\leq (\|\mathbf{M}_n(0)\|^2+\|\mathbf{H}_n(0)\|^2)e^{-C(\sigma,\tau,\chi_0)t} \\ &\quad +\int_0^{\omega}e^{-C\nu(t-s)}(C(\tau,\chi_0)\|\mathbb{F}_n\|^2 +(\frac{2}{\chi_0}+2\tau)\|\rho_n\|^2\|\mathbf{l}_n\|^2)ds. \end{align*} Since $\{\varphi_k(x)\}_{k\in\mathbb{N}}$ and $\{\nabla\varphi_k(x)\}_{k\in\mathbb{N}}$ are orthonormal bases of $\mathbb{H}_0^1(\mathbb{T}^3)$ and $\mathbb{H}^2(\mathbb{T}^3)$, respectively, we have $|\mathbf{h}_k(t)|^2=\|\mathbf{H}_n(t)\|^2$ and $|\mathbf{m}_k(t)|^2=\|\mathbf{M}_n(t)\|^2$, from which we conclude that $T=\omega$. Define the ball $B_{R}$ of radius $R$ and the map $\Pi:B_{R}\to B_{R}$ such that \[ \Pi((\mathbf{M}_n(0),\mathbf{H}_n(0)))=(\mathbf{M}_n(\omega), \mathbf{H}_n(\omega)), \] where the radius $R$ satisfies \[ R\geq\Big(\frac{\int_0^{\omega}e^{-C\nu(t-s)}(C(\tau, \chi_0)(\|\mathbb{F}_n\|^2+\|\partial_t\mathbb{F}_n\|^2) +(\frac{2}{\chi_0}+2\tau)\|\rho_n\|^2\|\mathbf{l}_n\|^2)ds} {1-e^{-C(\sigma,\tau,\chi_0)\omega}}\Big)^{1/2}. \] By the same process as in \cite{Prouse}, we can prove that the map $\Pi$ is continuous. Hence, it has a fixed point. Moreover, from \eqref{E3-12}, we know that the solution operators $(\mathbf{u}_n,\Omega_n)\to \mathbf{M}_n(\mathbf{u}_n,\Omega_n)$ and $(\mathbf{u}_n,\Omega_n)\to \mathbf{H}_n(\mathbf{u}_n,\Omega_n)$ maps bounded sets in $\mathbb{C}([0,T],\mathbb{X}_n)$ into bounded subsets of the set $\mathbb{Y}_n$ and $\mathbb{Z}_n$. Then, as done in \cite{Hu2}, the solution operators $(\mathbf{u}_n,\Omega_n)\to \mathbf{M}_n(\mathbf{u}_n,\Omega_n)$ and $(\mathbf{u}_n,\Omega_n)\to \mathbf{H}_n(\mathbf{u}_n,\Omega_n)$ are continuous. This completes the proof. \end{proof} The following result shows the local solvability of problem \eqref{E3-2} and \eqref{E3-2R}. The proof is inspired by the work of \cite{Ami,Fei1,Fei3}. \begin{lemma} \label{lem3} For initial $(\mathbf{u}_n(0),\Omega_n(0))\in\mathbb{X}_n\times\mathbb{X}_n$, system \eqref{E3-2}-\eqref{E3-2R} has a unique solution $(\mathbf{u}_n(t,x),\Omega_n(t,x))$ in $\mathbb{C}([0,\omega]; \mathbb{X}_n\times\mathbb{X}_n)$. \end{lemma} \begin{proof} Define a family of linear self-adjoint operators \begin{gather*} \mathcal{F}_1[\rho_n]:\mathbb{X}_n\to \mathbb{X}^*_n,\quad (\mathcal{F}_1[\rho_n]\mathbf{u}_n,\mathbf{w}) =\int_{\mathbb{T}^3}\rho_n\mathbf{u}_n\cdot\mathbf{w}dx, \\ \mathcal{F}_2[\rho_n]:\mathbb{X}_n\to \mathbb{X}^*_n,\quad (\mathcal{F}_2[\rho_n]\Omega_n,\mathbf{w}) =\int_{\mathbb{T}^3}\rho_n\Omega_n\cdot\mathbf{w}dx. \end{gather*} By \eqref{E3-4}, we can get that those operators are invertible and the inverse operators $\mathcal{F}_1^{-1}[\rho_n], \mathcal{F}_2^{-1}[\rho_n]\in\mathcal{L}(\mathbb{X}_n^*,\mathbb{X}_n)$ are well defined and locally Lipschitz continuous on the set $\{\rho_n\geq\frac{\tilde{\rho}(0)}{2}\}$. Furthermore, problem \eqref{E3-2}-\eqref{E3-2R} can rewritten as the integral equations \begin{gather}\label{E3-5} \mathbf{u}_n(t)=\mathcal{F}_1^{-1}[\rho_n(t)] \Big(\mathcal{F}_1[\rho_n(0)]\mathbf{u}_n(0)+\int_{0}^tG_1[\rho_n, \mathbf{u}_n,\Omega_n,\mathbf{M}_n,\mathbf{H}_n](s)ds\Big), \\ \label{E3-5R} \Omega_n(t)=\mathcal{F}_2^{-1}[\rho_n(t)] \Big(\mathcal{F}_2[\rho_n(0)]\Omega_n(0)+\int_{0}^tG_2[\rho_n, \mathbf{u}_n,\Omega_n,\mathbf{M}_n,\mathbf{H}_n](s)ds\Big), \end{gather} where \begin{align*} &(G_1[\rho_n,\mathbf{u}_n,\Omega_n,\mathbf{M}_n,\mathbf{H}_n],\mathbf{w}) \\ &=\int_{\mathbb{T}^3}\left(\rho_n(\mathbf{u}_n\otimes \mathbf{u}_n)\cdot\nabla\mathbf{w}-\mu\nabla \mathbf{u}_n\cdot\nabla\mathbf{w}-(\mu+\lambda)\operatorname{div} \mathbf{u}_n\cdot\operatorname{div}\mathbf{w}\right)dx \\ &\quad+\int_{\mathbb{T}^3}\left(a(\rho_n)^{\gamma}\operatorname{div} \mathbf{w}+\delta(\rho_n)^{\beta}\operatorname{div}\mathbf{w} +\frac{\mu_0}{2}|\mathbf{M}_n|^2\operatorname{div}\mathbf{w} -\epsilon\nabla\mathbf{u}_n\cdot\nabla\rho_n\cdot\mathbf{w}\right)dx \\ &\quad +\int_{\mathbb{T}^3}\left(-2\epsilon\rho_n\mathbf{u}_n\cdot\mathbf{w} +\mu_0\mathbf{M}_n\cdot\nabla\mathbf{H}_n\cdot\mathbf{w} +2\zeta(\nabla\times\Omega_n)\cdot\mathbf{w}\right)dx \\ &\quad +\int_{\mathbb{T}^3}\left(-\zeta(\nabla\times\mathbf{u}_n) \cdot(\nabla\times\mathbf{w}) +\rho_n\mathbf{f}_n\cdot \mathbf{w}\right)dx, \end{align*} and \begin{align*} &(G_2[\rho_n,\mathbf{u}_n,\Omega_n,\mathbf{M}_n,\mathbf{H}_n],\mathbf{w}) \\ &=\int_{\mathbb{T}^3}\left(\rho_n(\mathbf{u}_n\otimes \Omega_n)\cdot\nabla\mathbf{w}-\mu'\nabla \Omega_n\cdot\nabla\mathbf{w}-(\mu'+\lambda')\operatorname{div} \Omega_n\cdot\operatorname{div}\mathbf{w}\right)dx \\ &\quad +\int_{\mathbb{T}^3}\left(-2\epsilon\rho_n\Omega_n\cdot\mathbf{w} -\epsilon\nabla\Omega_n\cdot\nabla\rho_n\cdot\mathbf{w} +\mu_0(\mathbf{M}_n\times\mathbf{H}_n)\cdot\mathbf{w}\right)dx \\ &\quad +\int_{\mathbb{T}^3}\left(2\zeta(\nabla\times\Omega_n-2\Omega_n) \cdot\mathbf{w}+\rho_n\mathbf{g}_n\cdot \mathbf{w}\right)dx. \end{align*} By Lemmas \ref{lem1}-\ref{lem2}, we know that the nonlinear term \[ G_j:\mathbb{W}^{1,2}(\mathbb{T}^3)\times\mathbb{X}_n\times\mathbb{X}_n \times\mathbb{Y}_n\times\mathbb{Z}_n\to \mathbb{X}_n^*,\quad j=1,2 \] are locally Lipschitz. This together with the properties of $\rho_n[\mathbf{u}_n]$, $\mathbf{M}_n(\mathbf{u}_n,\Omega_n)$, $\mathbf{H}_n(\mathbf{u}_n,\Omega_n)$, $\mathcal{F}_1^{-1}[\rho_n]$ and $\mathcal{F}_2^{-1}[\rho_n]$ implies that integral equations \eqref{E3-5}-\eqref{E3-5R} possess unique solutions $(\mathbf{u}_n,\Omega_n)\in\mathbb{C}([0,T);\mathbb{X}_n\times\mathbb{X}_n)$ by the Banach fixed point theorem, for sufficiently small $T$. Then, by Lemma \ref{lem3} and \eqref{E3-15}, we know that $\rho_n$, $\mathbf{u}_n$, $\Omega_n$, $\mathbf{M}_n$ and $\mathbf{H}_n$ are bounded in $[0,T)$. Hence, $(\mathbf{u}_n,\Omega_n)$ can be prolonged up to $T=\omega$, which implies that the solution $(\rho_n, \mathbf{u}_n, \Omega_n, \mathbf{M}_n, \mathbf{H}_n)$ is global. \end{proof} Define a bounded convex set \begin{align*} \Gamma&=\Big\{[\rho_n,\mathbf{u}_n, \Omega_n, \mathbf{M}_n, \mathbf{H}_n]|\|\rho_n\|_{\mathbb{W}^{1,2}(\mathbb{T}^3) \cap\mathbb{C}(\mathbb{T}^3)}\leq R_1,\, \rho_n\geq R_2,\, \int_{\mathbb{T}^3}\rho_n dx\leq m_1, \\ &\quad \frac{1}{2}\int_{\mathbb{T}^3}|\mathbf{u}_n|^2dx\leq\frac{E_1}{R_2},\, \frac{1}{2}\int_{\mathbb{T}^3}|\Omega_n|^2dx\leq\frac{E_1}{R_2},\, \|\mathbf{M}_n\|_{\mathbb{W}^{1,2}(\mathbb{T}^3)\cap\mathbb{C}(\mathbb{T}^3)} \leq R_3, \\ &\quad \|\mathbf{H}_n\|_{\mathbb{W}^{1,2}(\mathbb{T}^3) \cap\mathbb{C}(\mathbb{T}^3)}\leq R_3,\, \int_{\mathbb{T}^3}\frac{a}{\gamma-1}\rho^{\gamma} +\frac{\delta}{\beta-1}\rho^{\beta}dx\leq E_1\Big\}, \end{align*} a map $\mathcal{P}_2$ from $[\mathbb{C}(\mathbb{T}^3)\cap\mathbb{W}^{1,2}(\mathbb{T}^3)] \times\mathbf{X}_n\times\mathbf{X}_n\times\mathbf{Y}_n\times\mathbf{Z}_n$ to itself, and \[ \mathcal{P}_1([\rho,\mathbf{u}_n, \Omega_n, \mathbf{M}_n, \mathbf{H}_n]) =[\rho_n,r\mathbf{u}_n, r\Omega_n, r\mathbf{M}_n, r\mathbf{H}_n], \] with \begin{align*} r&=r(\rho_n,\mathbf{u}_n, \Omega_n, \mathbf{M}_n, \mathbf{H}_n)\\ &=\min\big\{1,\big(\frac{E_1+\varepsilon-(\int_{\mathbb{T}^3} \frac{a}{\gamma-1}\rho^{\gamma}+\frac{\delta}{\beta-1}\rho^{\beta}dx)} {\frac{1}{2}\int_{\mathbb{T}^3}(\rho_n|\mathbf{u}_n| +\rho_n|\Omega_n|+\mu_0(|\mathbf{M}_n|^2+|\mathbf{H}_n|^2))dx}\big)^{1/2}\big\}. \end{align*} Then, by the properties of $\rho_n$, $\mathbf{u}_n$, $\Omega_n$, $\mathbf{M}_n$ and $\mathbf{H}_n$ obtained in Lemmas \ref{lem1}-\ref{lem3}, we get that the continuous map $\mathcal{P}_1:\Gamma\to \mathbb{D}_{E_1}\subset\Gamma$ and $\mathcal{P}_1|_{\mathbb{D}_{E_1-\varepsilon}}=Id$, where \begin{align*} \mathbb{D}_{E_1} &=\Big\{[\rho_n,\mathbf{u}_n, \Omega_n, \mathbf{M}_n, \mathbf{H}_n]| \|\rho_n\|_{\mathbb{W}^{1,2}(\mathbb{T}^3)\cap\mathbb{C}(\mathbb{T}^3)} \leq R_1,\, \rho_n\geq R_2,\quad \int_{\mathbb{T}^3}\rho_n dx\leq m_1, \\ &\quad \|\mathbf{M}_n\|_{\mathbb{W}^{1,2}(\mathbb{T}^3)\cap\mathbb{C} (\mathbb{T}^3)}\leq R_3,\, \|\mathbf{H}_n\|_{\mathbb{W}^{1,2} (\mathbb{T}^3)\cap\mathbb{C}(\mathbb{T}^3)}\leq R_3,\, E(t)\leq E_1\Big\} \end{align*} and \begin{align*} \mathbb{D}_{E_1-\epsilon} &=\Big\{[\rho_n,\mathbf{u}_n, \Omega_n, \mathbf{M}_n, \mathbf{H}_n]| \|\rho_n\|_{\mathbb{W}^{1,2}\cap\mathbb{C}(\mathbb{T}^3)}\leq R_1,\, \rho_n\geq R_2,\, \int_{\mathbb{T}^3}\rho_n dx\leq m_1, \\ &\quad \|\mathbf{M}_n\|_{\mathbb{W}^{1,2}\cap\mathbb{C}(\mathbb{T}^3)}\leq R_3,\, \|\mathbf{H}_n\|_{\mathbb{W}^{1,2}(\mathbb{T}^3)\cap\mathbb{C}(\mathbb{T}^3)} \leq R_3,\, E(t)\leq E_1-\varepsilon\Big\}. \end{align*} Moreover, the map $\mathcal{P}_2$ is continuous and compact, $\mathcal{P}_2\circ\mathcal{P}_1$ maps $\mathbb{D}_{E_1}$ into $\mathbb{D}_{E_1-\varepsilon}$. Then, using Banach fixed theorem and Lemmas \ref{lem1}-\ref{lem3}, we summarize the existence of periodic solution for system \eqref{E3-1}-\eqref{E3-3R} as follows. \begin{lemma} \label{lem4} Assume that $\epsilon$, $\delta$ and $\beta$ are given positive parameters. Then, for any fixed $n$, systems \eqref{E3-1}-\eqref{E3-3} have time-periodic solutions $\rho_n$, $\mathbf{u}_n$, $\Omega_n$, $\mathbf{M}_n$ and $\mathbf{H}_n$. Moreover, $\rho_n\in\mathbb{C}^1(\mathbb{S}^1;\mathbb{C}^2(\mathbb{T}^3))$ is a classical solution of \eqref{E3-1} on $\mathbb{S}^1$, and there exists $R_2$ depending on $n$ such that \[ \rho_n\geq R_2(n)>0,\quad \int_{\mathbb{T}^3}\rho_n(t)dx=m_{\epsilon},\quad\text{with } |\mathbb{T}^3|M(m_{\epsilon})=2\epsilon m_{\epsilon}. \] The energy inequality \begin{equation}\label{E6-1} \begin{aligned} &E_n(t)+C\int_0^tE_{nf}(s)ds\\ &\leq E(0)+C\int_0^t(1+\|\mathbb{F}(s)\|^2+\|\partial_t\mathbb{F}(s)\|^2)ds \\ &\quad +C\int_0^t(1+\|\mathbf{f}(s)\|^2_{\mathbb{L}^{\infty}(\mathbb{Q})} +\|\mathbf{g}(s)\|^2_{\mathbb{L}^{\infty}(\mathbb{Q})} +\|\mathbf{l}(s)\|^2_{\mathbb{L}^{\infty}(\mathbb{Q})})ds \end{aligned} \end{equation} holds on $\mathbb{S}^1$, where \begin{gather*} E_n(t)=\int_{\mathbb{T}^3}\Big(\frac{1}{2}\rho(\mathbf{u}_n^2+\Omega_n^2) +\frac{a}{\gamma-1}\rho_n^{\gamma}+\frac{\mu_0}{2}(|\mathbf{M}_n|^2 +|\mathbf{H}_n|^2)\Big)dx,\\ \begin{aligned} E_{nf}(t) &= \int_{\mathbb{T}^3}(\mu|\nabla\mathbf{u}_n|^2+\mu'|\nabla\Omega_n|^2 +(\mu+\lambda)|\operatorname{div}\mathbf{u}_n|^2 +(\mu'+\lambda')|\operatorname{div}\Omega_n|^2)dx\\ &\quad +\int_{\mathbb{T}^3}(\frac{\mu_0}{\tau}|\mathbf{M}_n|^2 +\frac{\mu_0(2\chi_0+1)}{\tau}|\mathbf{H}_n|^2 +\mu_0\sigma|\nabla\times\mathbf{M}_n|^2)dx\\ &\quad +\int_{\mathbb{T}^3}2\mu_0\sigma|\operatorname{div}\mathbf{M}_n|^2 +\zeta|\nabla\times\mathbf{u}_n-2\Omega_n|^2dx\,, \end{aligned}\\ E(0)=\int_{\mathbb{T}^3}\Big(\frac{1}{2}\rho_0(\mathbf{u}_0^2+\Omega^2_0) +\frac{a}{\gamma-1}\rho_0^{\gamma}+\frac{\mu_0}{2}(|\mathbf{M}_0|^2 +|\mathbf{H}_0|^2)\Big)dx. \end{gather*} Moreover, there exits a constant $E_1$ depending on $\delta$, $\epsilon$ such that \[ E_n[\rho_n,\mathbf{u}_n, \Omega_n, \mathbf{M}_n, \mathbf{H}_n](0)\leq E_1. \] \end{lemma} In what follows, passing to the limit as $n\to \infty$ in the sequence of approximation solution $\{\rho_n,\mathbf{u}_n,\Omega_n, \mathbf{M}_n, \mathbf{H}_{n}\}$ is the main target, which is treated similarly to \cite{Ami,Fei1,Fei2,Fei3,Lions3}. The following result is directly from the energy inequality \eqref{E6-1} and some properties about the $L^p$-theory of parabolic equations. \begin{lemma} \label{lem5} Assume that $\beta\geq4$. Then $(\rho_n,\mathbf{u}_n, \Omega_n, \mathbf{M}_n, \mathbf{H}_n)$, the approximation solutions constructed above, satisfy the following estimates: \begin{gather*} \|\rho_n\|_{\mathbb{L}^{\infty}(\mathbb{S}^1;\mathbb{L}^{\beta}(\mathbb{T}^3))} \leq C,\quad \|\mathbf{u}_n\|_{\mathbb{L}^2(\mathbb{S}^1;\mathbb{H}_0^1 (\mathbb{T}^3))} \leq C, \quad \|\Omega_n\|_{\mathbb{L}^2 (\mathbb{S}^1;\mathbb{H}_0^1(\mathbb{T}^3))}\leq C \\ \|\rho_n^{1/2}|\mathbf{u}_n\|_{\mathbb{L}^{\infty} (\mathbb{S}^1;\mathbb{L}^{\beta}(\mathbb{T}^3))}\leq C,\quad \|\rho_n^{1/2}|\Omega_n\|_{\mathbb{L}^{\infty}(\mathbb{S}^1;\mathbb{L}^{\beta} (\mathbb{T}^3))}\leq C, \\ \|\mathbf{M}_n\|_{\mathbb{L}^{\infty}(\mathbb{S}^1;\mathbb{L}^2(\mathbb{T}^3))} \leq C,\quad \|\mathbf{M}_n\|_{\mathbb{L}^2(\mathbb{S}^1;\mathcal{M})}\leq C, \\ \|\mathbf{H}_n\|_{\mathbb{L}^{\infty}(\mathbb{S}^1;\mathbb{L}^2(\mathbb{T}^3))} \leq C,\quad \|\mathbf{H}_n\|_{\mathbb{L}^2(\mathbb{S}^1;\mathcal{M})}\leq C. \end{gather*} Furthermore, for $12$, $r_2>\beta$ and $q_1=\frac{6\beta}{4\beta+3}$, it follows that \begin{gather}\label{E3-28} \epsilon\int_{\mathbb{Q}}|\rho_n|^2+|\nabla\rho_n|^2\,dx\,dt\leq C,\quad \int_{\mathbb{Q}}|\nabla\rho_n|^{r_1}+|\rho_n|^{r_2}\,dx\,dt\leq C\\ \label{E3-29} \int_{\mathbb{S}^1}\|\partial_t\rho_n\|^p_{\mathbb{L}^{q_1}(\mathbb{T}^3)} +\|\Delta\rho_n\|_{\mathbb{L}^{q_1}(\mathbb{T}^3)}^pdt\leq C, \end{gather} where $C$ is a constant independent of $n$. \end{lemma} To obtain the weak solution of the approximation equations \eqref{E3-1}-\eqref{E3-3R}, the first step is necessary to pass to the limit, as $n\to \infty$, in the sequence of approximate solutions $(\rho_n,\mathbf{u}_n,\Omega_n,\mathbf{M}_n,\mathbf{H}_n)$ constructed Lemma \ref{lem4}. We use the Lions-Aubin lemma (see \cite{Lions4}), as done in \cite{Ami,Fei2}, and obtain some convergence results about $\rho_n$, $\mathbf{u}_n$ and $\Omega_n$: \begin{gather}\label{E6-2} \rho^{\gamma}_n\to \rho^{\gamma},\quad \rho^{\beta}_n\to \rho^{\beta},\quad \text{strongly in $\mathbb{L}^{r_3}(\mathbb{Q})$ for a certain }r_3>1, \\ \nabla\rho_n\nabla\mathbf{u}_n^i\to \nabla\rho\nabla\mathbf{u}^i,\quad \text{weakly in $\mathbb{L}^{r_3}(\mathbb{Q})$ for a certain }r_3>1, \nonumber\\ \nabla\rho_n\nabla\Omega_n\to \nabla\rho\nabla\Omega,\quad \text{weakly in $\mathbb{L}^{r_3}(\mathbb{Q})$ for a certain }r_3>1, \nonumber\\ \rho_n\mathbf{u}_n^i\to \rho\mathbf{u}^i,\quad \text{strongly in }\mathbb{L}^2(\mathbb{S}^1;\mathbb{W}^{-1,2}(\mathbb{T}^3)), \nonumber\\ \rho_n\Omega_n^i\to \rho\Omega^i,\quad \text{ strongly in } \mathbb{L}^2(\mathbb{S}^1;\mathbb{W}^{-1,2}(\mathbb{T}^3)), \nonumber\\ \rho_n\mathbf{u}_n^i\mathbf{u}_n^j\to \rho\mathbf{u}^i\mathbf{u}^j,\quad \rho_n\mathbf{u}_n^i\Omega_n^j\to \rho\mathbf{u}^i\Omega^j,\quad \text{in }\mathcal{D}'(\mathbb{Q}), \nonumber \end{gather} where $i,j=1,2,3$. Now we pass to the limit, as $n\to \infty$, in the sequence of the approximate solutions $(\mathbf{M}_n,\mathbf{H}_n)$ of approximation equations \eqref{E3-3}-\eqref{E3-3R}. By the energy inequality \eqref{E3-15}, it is easy to see that \begin{gather*} \mathbf{M}_n\to \mathbf{M}\quad \text{ weakly* in $\mathbb{L}^{\infty}(\mathbb{S}^1;\mathbb{L}^2(\mathbb{T}^3))$ and in $\mathbb{L}^2(\mathbb{S}^1;\mathcal{M})$ weak},\\ \psi_n\to \psi\quad \text{weakly* in $\mathbb{L}^2(\mathbb{S}^1;\mathbb{H}^1(\mathbb{T}^3))$ and in $\mathbb{L}^2(\mathbb{S}^1;\mathbb{H}^2(\mathbb{T}^3))$ weak}, \\ \mathbf{H}_n\to \mathbf{H}\quad \text{weakly* in $\mathbb{L}^{\infty}(\mathbb{S}^1;\mathbb{L}^2(\mathbb{T}^3))$ and in $\mathbb{L}^2(\mathbb{S}^1;\mathcal{M})$ weak}. \end{gather*} Due to the Sobolev imbedding $\mathbb{H}^1(\mathbb{T}^3)\hookrightarrow\mathbb{L}^6(\mathbb{T}^3)$, $\Omega_n\in\mathbb{L}^2(\mathbb{S}^1;\mathbb{H}^1(\mathbb{T}^3))$ and $\mathbf{M}_n\in\mathbb{L}^{\infty}(\mathbb{S}^1;\mathbb{L}^2(\mathbb{T}^3))$, we deduce that \begin{equation}\label{E6-3} \Omega_n\times\mathbf{M}_n\in\mathbb{L}^2(\mathbb{S}^1;\mathbb{L}^{3/2} (\mathbb{T}^3)). \end{equation} Consider the orthogonal projections: \[ \Lambda_n:[\mathbb{L}^2(\mathbb{T}^3)]^3\to \mathbb{Z}_n,\quad \Upsilon_n=Id-\Lambda_n. \] It follows from \eqref{E3-3} and \eqref{E6-3} that \begin{equation}\label{E3-25} \partial_t(\Lambda_n\mathbf{M}_n)\quad\text{are bounded in } [\mathbb{L}^2(\mathbb{S}^1;\mathbb{W}^{-1,2}(\mathbb{T}^3))]^3. \end{equation} Then, by the Lions-Aubin lemma, for $-13. \end{gathered} \end{equation} This is the main results in this section. \begin{proposition} \label{prop2} Assume that $\beta>5$. Then there exists a time-periodic solution $(\rho,\mathbf{u},\Omega,\mathbf{M},\mathbf{H})$ of the problem in $\mathcal{D}'(\mathbb{Q})$, \begin{gather}\label{E4-RR1} \partial_t\rho+\operatorname{div}(\rho\mathbf{u})=0, \\ \label{E4-RR2} \begin{aligned} &\partial_t(\rho u^i)+\operatorname{div}(\rho u^i\otimes \mathbf{u})-\mu\Delta u^i \\ &-(\lambda+\mu)\partial_{x_i}(\operatorname{div}\mathbf{u}) -\partial_{x_i}P(\rho) \\ &+\delta\partial_{x_i}\rho^{\beta} =R+\rho f^i(t,x), \end{aligned} \\ \label{E4-RRRR1} \partial_t(\rho\Omega^i)+\operatorname{div}(\rho\Omega^i\otimes\mathbf{u} )-\mu'\Delta\Omega^i-(\lambda'+\mu')\nabla(\operatorname{div}\Omega) =S+\rho g^i(t,x), \\ \label{E4-RR3} \partial_t\mathbf{M}+\operatorname{div}(\mathbf{u}\otimes \mathbf{M})-\sigma\Delta\mathbf{M}+\frac{1}{\tau}(\mathbf{M}-\chi_0\mathbf{H}) =\Omega\times\mathbf{M}+\rho\mathbf{l}(t,x), \\ \label{E4-RRRR2} \mathbf{H}=\nabla\psi,\quad \operatorname{div}(\mathbf{H}+\mathbf{M}) =\mathbb{F}, \end{gather} where $i=1,2,3$, \begin{gather*} R=\mu_0\mathbf{M}\cdot\nabla\mathbf{H}-\zeta\nabla \times(\nabla\times u^i-2\Omega), \\ S=\mu_0\mathbf{M}\times\mathbf{H}+2\zeta(\nabla\times\mathbf{u}-2\Omega^i), \\ F_{\epsilon}=(2\mu+\lambda)\operatorname{div}\mathbf{u}_{\epsilon} +\frac{\mu_0}{2}|\mathbf{M}_0|^2-a\rho_{\epsilon}^{\gamma} -\delta\rho^{\beta}_{\epsilon}, \\ F_{\epsilon}\to \overline{F}\quad \text{weakly in } \mathbb{L}^{\frac{\beta+1}{\beta}}(\mathbb{Q}). \end{gather*} Moreover, $\rho$, $\mathbf{u}$ and $\Omega$ satisfy \eqref{E2-4}-\eqref{E2-4R2}, and \begin{gather*} \rho\in\mathbb{L}^{\infty}(\mathbb{S}^1;\mathbb{L}^{\beta}(\mathbb{T}^3)) \cap\mathbb{L}^{\beta+1}(\mathbb{Q}), \quad \mathbf{u},\Omega\in\mathbb{L}^2(\mathbb{S}^1;\mathbb{W}^{1,2}(\mathbb{T}^3)), \\ \mathbf{M}\in\mathbb{L}^{\infty}(\mathbb{S}^1;\mathbb{L}^2(\mathbb{T}^3)) \cap\mathbb{L}^2(\mathbb{S}^1;\mathcal{M}), \\ \mathbf{H}=\nabla\psi,\quad \psi\in\mathbb{L}^{\infty} (\mathbb{S}^1;\mathbb{H}^1(\mathbb{T}^3))\cap\mathbb{L}^2 (\mathbb{S}^1;\mathbb{H}^2(\mathbb{T}^3)). \end{gather*} Equation \eqref{E4-RR1} holds in the sense of renormalized solutions and for any continuous sublinear function $h$ on $\mathbb{R}^1$, then it holds \begin{equation}\label{E4-18} \int_{\mathbb{Q}}h(\rho)\operatorname{div}\mathbf{u}\,dx\,dt=0. \end{equation} The energy $E_{\delta}[\rho,\mathbf{u},\Omega,\mathbf{M},\mathbf{H}]$ satisfies \begin{equation}\label{E4-19} \begin{aligned} \frac{d}{dt}E_{\delta}(t)+CE_{\delta f}(t) &\leq C(1+\|\mathbb{F}(t)\|^2 +\|\partial_t\mathbb{F}(t)\|^2)\\ &\quad +C\Big(1+\int_{\mathbb{T}^3}\rho(|\mathbf{f}||\mathbf{u}| +|\mathbf{g}||\Omega|+|\mathbf{l}||\mathbf{M}|)dx\Big), \end{aligned} \end{equation} in $\mathcal{D}'(\mathbb{S}^1)$, where \begin{align*} E_{\delta f}(t) &=\int_{\mathbb{T}^3}(\mu|\nabla\mathbf{u}|^2+\mu'|\nabla\Omega|^2 +(\mu+\lambda)|\operatorname{div}\mathbf{u}|^2+(\mu'+\lambda') |\operatorname{div}\Omega|^2)dx\\ &\quad +\int_{\mathbb{T}^3}(\frac{\mu_0}{\tau}|\mathbf{M}|^2 +\frac{\mu_0(2\chi_0+1)}{\tau}|\mathbf{H}|^2 +\mu_0\sigma|\nabla\times\mathbf{M}|^2)dx\\ &\quad +\int_{\mathbb{T}^3}2\mu_0\sigma|\operatorname{div}\mathbf{M}|^2 +\zeta|\nabla\times\mathbf{u}-2\Omega|^2dx. \end{align*} \end{proposition} \begin{lemma} \label{lem6} Let $\rho\geq0$, $\mathbf{u}$, $\Omega$, $\mathbf{M}$ and $\mathbf{H}$ satisfy \begin{gather*} \rho\in\mathbb{C}(\mathbb{S}^1;\mathbb{L}^{\gamma}(\mathbb{T}^3)),\quad \|\rho\|_{\mathbb{C}(\mathbb{S}^1;\mathbb{L}^{1}(\mathbb{T}^3))}\leq m,\quad \mathbf{u}\in[\mathbb{L}^2(\mathbb{S}^1;\mathbb{W}^{1,2}(\mathbb{T}^3))]^3, \\ \Omega\in[\mathbb{L}^2(\mathbb{S}^1;\mathbb{W}^{1,2}(\mathbb{T}^3))]^3,\quad \text{$\mathbf{u}$ and $\Omega$ satisfy \eqref{E2-4R1} and \eqref{E2-4R2}}, \\ \mathbf{M}\in\mathbb{L}^{\infty}(\mathbb{S}^1;\mathbb{L}^2(\mathbb{T}^3)) \cap\mathbb{L}^2(\mathbb{S}^1;\mathcal{M}), \\ \mathbf{H}=\nabla\psi,\quad \psi\in\mathbb{L}^{\infty}(\mathbb{S}^1;\mathbb{H}^1(\mathbb{T}^3)) \cap\mathbb{L}^2(\mathbb{S}^1;\mathbb{H}^2(\mathbb{T}^3)). \end{gather*} Then there holds \[ \|E_{\delta}(t)\|_{\mathbb{C}(\mathbb{S}^1)} \leq C(1+\int_{\mathbb{Q}}P(\rho(t))\,dx\,dt), \] where $C$ is a constant depending on $\mu$, $\nu$, $\|f\|_{\mathbb{L}^{\infty}(\mathbb{Q})}$, $\|g\|_{\mathbb{L}^{\infty}(\mathbb{Q})}$ and $\|l\|_{\mathbb{L}^{\infty}(\mathbb{Q})}$, $P$ denotes a convex function such that $P(\rho(t))\geq\frac{a}{\gamma-1}\rho^{\gamma}(t)$ for $\gamma>\frac{5}{3}$, \[ E_{\delta}(t)=\int_{\mathbb{T}^3} \Big(\frac{1}{2}\rho(\mathbf{u}^2+\Omega^2)+P(\rho(t)) +\frac{\mu_0}{2}(|\mathbf{M}|^2+|\mathbf{H}|^2)\Big)dx \in\mathbb{L}^1(\mathbb{S}^1) \] satisfying the energy inequality \begin{equation}\label{E4-7} \begin{aligned} \frac{d}{dt}E_{\delta}(t)+CE_{\delta f}(t) &\leq C\Big(1+\|\mathbb{F}(t)\|^2+\|\partial_t\mathbb{F}(t)\|^2\Big) \\ &\quad +C(1+\int_{\mathbb{T}^3}\rho(|\mathbf{f}||\mathbf{u}| +|\mathbf{g}||\Omega|+|\mathbf{l}||\mathbf{M}|)dx) \end{aligned} \end{equation} and \begin{align*} E_{\delta f}(t) &= \int_{\mathbb{T}^3}(\mu|\nabla\mathbf{u}|^2 +\mu'|\nabla\Omega|^2+(\mu+\lambda)|\operatorname{div}\mathbf{u}|^2 +(\mu'+\lambda')|\operatorname{div}\Omega|^2)dx\\ &\quad +\int_{\mathbb{T}^3}(\frac{\mu_0}{\tau}|\mathbf{M}|^2 +\frac{\mu_0(2\chi_0+1)}{\tau}|\mathbf{H}|^2 +\mu_0\sigma|\nabla\times\mathbf{M}|^2)dx\\ &\quad +\int_{\mathbb{T}^3}2\mu_0\sigma|\operatorname{div}\mathbf{M}|^2 +\zeta|\nabla\times\mathbf{u}-2\Omega|^2dx. \end{align*} \end{lemma} \begin{proof} Integrating \eqref{E4-7} over $\mathbb{S}^1$, we have \begin{align*} \int_{\mathbb{S}^1}\int_{\mathbb{T}^3}E_{\delta f}(t)dx\, dt &\leq C\int_{\mathbb{S}^1}(1+\|\mathbb{F}(s)\|^2 +\|\partial_t\mathbb{F}(s)\|^2)dt \\ &\quad+C\Big(1+\int_{\mathbb{S}^1}\int_{\mathbb{T}^3} (|\mathbf{f}||\mathbf{u}|+|\mathbf{g}||\Omega|+|\mathbf{l}||\mathbf{M}|) \rho dx\,dt\Big). \end{align*} By the Poincar\'{e} inequality, the standard Sobolev embedding theorem and the definition of $E_{\delta f}(t)$, \begin{align*} &\int_{\mathbb{S}^1}(\|\mathbf{u}\|^2_{\mathbb{L}^6(\mathbb{T}^3)} +\|\Omega\|^2_{\mathbb{L}^6(\mathbb{T}^3)} +\|\mathbf{M}\|^2_{\mathbb{L}^2(\mathbb{T}^3)} +\|\mathbf{H}\|^2_{\mathbb{L}^2(\mathbb{T}^3)})dt \\ &\leq C(1+\|\rho\|_{\mathbb{C}(\mathbb{S}^1;\mathbb{L}^{6/5} (\mathbb{T}^3))}\int_{\mathbb{S}^1}(\|\mathbf{u}\|_{\mathbb{L}^6 (\mathbb{T}^3)}+\|\Omega\|_{\mathbb{L}^6(\mathbb{T}^3)})dt\\ &\quad +\|\rho\|_{\mathbb{C}(\mathbb{S}^1;\mathbb{L}^2(\mathbb{T}^3))} \int_{\mathbb{S}^1}(\|\mathbf{M}\|_{\mathbb{L}^2(\mathbb{T}^3)} +\|\mathbf{H}\|_{\mathbb{L}^2(\mathbb{T}^3)})dt), \end{align*} which implies \begin{equation}\label{E4-8} \begin{aligned} &\int_{\mathbb{S}^1}(\|\mathbf{u}\|^2_{\mathbb{L}^6(\mathbb{T}^3)} +\|\Omega\|^2_{\mathbb{L}^6(\mathbb{T}^3)} +\|\mathbf{M}\|^2_{\mathbb{L}^2(\mathbb{T}^3)} +\|\mathbf{H}\|^2_{\mathbb{L}^2(\mathbb{T}^3)})dt\\ &\leq C\Big(1+\|\rho\|_{\mathbb{C}(\mathbb{S}^1;\mathbb{L}^{6/5} (\mathbb{T}^3))}\Big)^2. \end{aligned} \end{equation} It follows from the H\"{o}lder inequality and \eqref{E4-8} that \begin{equation}\label{E4-10} \begin{aligned} &\int_{\mathbb{Q}}\frac{1}{2}\rho(t)(|\mathbf{u}(t)|^2 +|\Omega(t)|^2)+\frac{\mu_0}{2}(|\mathbf{M}(t)|^2+|\mathbf{H}(t)|^2)\,dx\,dt \\ &\leq \|\rho\|_{\mathbb{C}(\mathbb{S}^1;\mathbb{L}^{3/2}(\mathbb{T}^3))} \int_{\mathbb{S}^1}(\|\mathbf{u}\|^2_{\mathbb{L}^6(\mathbb{T}^3)} +\|\Omega\|^2_{\mathbb{L}^6(\mathbb{T}^3)})dt\\ &\quad +\frac{\mu_0}{2}\int_{\mathbb{S}^1}(\|\mathbf{M}\|^2_{\mathbb{L}^2(\mathbb{T}^3)} +\|\mathbf{H}\|^2_{\mathbb{L}^2(\mathbb{T}^3)})dt \\ &\leq C\Big(1+\|\rho\|_{\mathbb{C}(\mathbb{S}^1;\mathbb{L}^{6/5} (\mathbb{T}^3))}\Big)^2 \Big(1+\|\rho\|_{\mathbb{C}(\mathbb{S}^1;\mathbb{L}^{3/2}(\mathbb{T}^3))}\Big). \end{aligned} \end{equation} By \eqref{E4-7}, we have \begin{equation}\label{E4-9} \|E_{\delta}(t)\|_{\mathbb{C}(\mathbb{S}^1)} \leq C\Big(1+\int_{\mathbb{Q}}\rho(|\mathbf{u}|^2+|\Omega|^2) +|\mathbf{M}|^2+|\mathbf{H}|^2+P(\rho)\,dx\,dt\Big). \end{equation} Then, by \eqref{E4-10}-\eqref{E4-9} and the interpolation \[ \|\rho\|_{\mathbb{L}^{3/2}(\mathbb{T}^3)},~\|\rho\|_{\mathbb{L}^{6/5}(\mathbb{T}^3)}\leq\|\rho\|_{\mathbb{L}^{1}(\mathbb{T}^3)}^{1-\varpi}\|\rho\|_{\mathbb{L}^{\gamma}(\mathbb{T}^3)}^{\varpi} \] for $\gamma>\frac{5}{3}$, we can complete the proof. \end{proof} \begin{lemma} \label{lem7} Assume that $\beta>4$. Then \begin{gather*} \epsilon\int_{\mathbb{S}^1}\|\rho_{\epsilon}\|^2_{\mathbb{W}^{1,2} (\mathbb{T}^3)}dt,\quad \int_{\mathbb{S}^1}\|\mathbf{u}_{\epsilon}\|^2_{\mathbb{W}^{1,2}(\mathbb{T}^3)}dt, \quad \int_{\mathbb{S}^1}\|\Omega_{\epsilon}\|^2_{\mathbb{W}^{1,2}(\mathbb{T}^3)}dt,\\ \int_{\mathbb{S}^1}\|\mathbf{M}_{\epsilon}\|^2_{\mathbb{W}^{1,2}(\mathbb{T}^3)}dt, \quad \int_{\mathbb{S}^1}\|\mathbf{H}_{\epsilon}\|^2_{\mathbb{W}^{1,2} (\mathbb{T}^3)}dt,\quad \|\rho\|^{\beta+1}_{\mathbb{L}^{\beta+1} (\mathbb{Q})},\\ \|E_{\delta}[\rho_{\epsilon},\mathbf{u}_{\epsilon}, \Omega_{\epsilon},\mathbf{M}_{\epsilon},\mathbf{H}_{\epsilon}]\|_{\mathbb{C}(\mathbb{S}^1)} \end{gather*} are bounded independently of $\epsilon>0$. \end{lemma} \begin{proof} Integrating the energy inequality \eqref{E4-7} over $\mathbb{S}^1$, we deduce \begin{align*} &\int_{\mathbb{S}^1}\|\mathbf{u}\|^2_{\mathbb{W}^{1,2}(\mathbb{T}^3)} +\|\Omega\|^2_{\mathbb{W}^{1,2}(\mathbb{T}^3)} +\|\mathbf{M}\|^2_{\mathbb{W}^{1,2}(\mathbb{T}^3)} +\|\mathbf{H}\|^2_{\mathbb{W}^{1,2}(\mathbb{T}^3)}dt \\ &\leq C\Big(1+\int_{\mathbb{Q}}(|\mathbf{u}|+|\Omega|+|\mathbf{M}| +|\mathbf{H}|)\rho \,dx\,dt\Big). \end{align*} It follows from the Sobolev embedding theorem that \begin{equation}\label{E4-1'} \begin{aligned} &\int_{\mathbb{S}^1}\|\mathbf{u}\|^2_{\mathbb{W}^{1,2}(\mathbb{T}^3)} +\|\Omega\|^2_{\mathbb{W}^{1,2}(\mathbb{T}^3)} +\|\mathbf{M}\|^2_{\mathbb{W}^{1,2}(\mathbb{T}^3)} +\|\mathbf{H}\|^2_{\mathbb{W}^{1,2}(\mathbb{T}^3)}dt \\ &\leq C\Big(1+\|\rho\|_{\mathbb{C}(\mathbb{S}^1;\mathbb{L}^{6/5}(\mathbb{T}^3))} \Big), \end{aligned} \end{equation} which implies \begin{gather*} \int_{\mathbb{S}^1}\|\mathbf{u}_{\epsilon}\|^2_{\mathbb{W}^{1,2}(\mathbb{T}^3)}dt,\quad \int_{\mathbb{S}^1}\|\Omega_{\epsilon}\|^2_{\mathbb{W}^{1,2}(\mathbb{T}^3)}dt,\\ \int_{\mathbb{S}^1}\|\mathbf{M}_{\epsilon}\|^2_{\mathbb{W}^{1,2}(\mathbb{T}^3)}dt,\quad \int_{\mathbb{S}^1}\|\mathbf{H}_{\epsilon}\|^2_{\mathbb{W}^{1,2}(\mathbb{T}^3)}dt \end{gather*} are bounded independently of $\epsilon>0$. Multiplying the continuity equation \eqref{E3-R1} by $\rho$ and integrating by parts, we obtain \[ \epsilon\int_{\mathbb{S}^1}\|\rho\|^2_{\mathbb{W}^{1,2}(\mathbb{T}^3)}dt \leq C\Big(1+\int_{\mathbb{S}^1}\|\mathbf{u}\|_{\mathbb{W}^{1,2} (\mathbb{T}^3)}\|\rho\|^2_{\mathbb{L}^{4}(\mathbb{T}^3)}dt\Big). \] The estimate of $\|\rho\|^{\beta+1}_{\mathbb{L}^{\beta+1}(\mathbb{Q})}$ is similar with \cite{Fei1}. Let $b$ be an convex function. For convenience, we take $b(x)=x\log x$. Multiplying the continuity equation \eqref{E3-R1} by $b'(\rho_{\epsilon})$ and integrating by parts, we have \[ \int_{\mathbb{Q}}(b'(\rho_{\epsilon})\rho_{\epsilon}-b(\rho_{\epsilon})) \operatorname{div}\mathbf{u}_{\epsilon}+2\epsilon\rho_{\epsilon} b'(\rho_{\epsilon})-\frac{2\epsilon m_{\epsilon}}{|\mathbb{T}^3|} b'(\rho_{\epsilon})\,dx\,dt\leq0, \] which leads to \begin{equation}\label{E4-6} \int_{\mathbb{Q}}\rho_{\epsilon}(\operatorname{div}\mathbf{u}_{\epsilon})\,dx\,dt \leq\epsilon C, \end{equation} where $C$ depends on $|\mathbb{Q}|$ and $m$. Denote $\mathcal{R}_{i,j}=\partial_{x_i}\Delta^{-1}\partial_{x_j}$. Note that \begin{equation}\label{E4-2} \begin{aligned} &\int_{\mathbb{Q}}a\rho_{\epsilon}^{\gamma+1}+\delta\rho_{\epsilon}^{\beta+1} +\frac{\mu_0}{2}\rho_{\epsilon}|\mathbf{M}_{\epsilon}|^2\,dx\,dt\\ &= \int_{\mathbb{Q}}\frac{m}{|\mathbb{T}^3|}(a\rho_{\epsilon}^{\gamma} +\delta\rho_{\epsilon}^{\beta}) +(\lambda+2\mu)\rho_{\epsilon}(\operatorname{div}\mathbf{u}_{\epsilon})\,dx\,dt +I_1+I_2+I_3+I_4, \end{aligned} \end{equation} where $\mathcal{A}_i[\rho_{\epsilon}]$ are a test functions of \eqref{E3-R2} and \begin{gather*} I_1=\int_{\mathbb{Q}}\rho_{\epsilon}\mathbf{u}_{\epsilon}^i\mathcal{A}_{i} [2\epsilon\rho_{\epsilon}-\epsilon\Delta\rho_{\epsilon}]\,dx\,dt, \\ I_2=\int_{\mathbb{Q}}\mathbf{u}_{\epsilon}^i(\rho_{\epsilon}\mathcal{R}_{i,j} [\rho_{\epsilon}\mathbf{u}_{\epsilon}^j]-\rho_{\epsilon}\mathbf{u}_{\epsilon}^j \mathcal{R}_{i,j}[\rho_{\epsilon}]) \,dx\,dt, \\ I_3=\int_{\mathbb{Q}}(\epsilon\nabla\rho_{\epsilon}\nabla\mathbf{u}_{\epsilon} +2\epsilon\rho_{\epsilon}\mathbf{u}_{\epsilon}-f^i\rho_{\epsilon})\mathcal{A}_{i} [\rho_{\epsilon}]\,dx\,dt, \\ I_4=\int_{\mathbb{Q}}\mu_0\mathbf{M}_{\epsilon}\cdot\nabla\mathbf{H}_{\epsilon} \cdot\mathcal{A}_{i}[\rho_{\epsilon}]\,dx\,dt, \\ I_5=\int_{\mathbb{Q}}-\zeta\nabla\times(\nabla\times \mathbf{u}_{\epsilon} -2\Omega_{\epsilon})\mathcal{A}_{i}[\rho_{\epsilon}]\,dx\,dt. \end{gather*} By the H\"{o}lder inequality, \eqref{E4-1} and \eqref{ER}, we obtain the following estimates: \begin{gather}\label{E4-3} |I_1|\leq\epsilon C\|\rho\|_{\mathbb{C}(\mathbb{S}^1;\mathbb{L}^3 (\mathbb{T}^3))}\Big(\int_{\mathbb{S}^1}\|\mathbf{u}\|^2_{\mathbb{L}^6 (\mathbb{T}^3)}dt\Big)^{1/2} \Big(\int_{\mathbb{S}^1}\|\rho\|^2_{\mathbb{W}^{1,2}(\mathbb{T}^3)}dt\Big)^{1/2},\\ \label{E4-4} |I_2|\leq C\|\rho\|^2_{\mathbb{C}(\mathbb{S}^1;\mathbb{L}^3(\mathbb{T}^3))} \Big(\int_{\mathbb{S}^1}\|\mathbf{u}\|^2_{\mathbb{L}^6(\mathbb{T}^3)}dt\Big),\\ \label{E4-5} \begin{aligned} |I_3|&\leq \epsilon C\|\rho\|_{\mathbb{C}(\mathbb{S}^1;\mathbb{L}^4 (\mathbb{T}^3))}\Big(\int_{\mathbb{S}^1}\|\mathbf{u}\|^2_{\mathbb{W}^{1,2} (\mathbb{T}^3)}dt\Big)^{1/2} \Big(\int_{\mathbb{S}^1}\|\rho\|^2_{\mathbb{W}^{1,2}(\mathbb{T}^3)}dt\Big)^{1/2} \\ &\quad +C\int_{\mathbb{S}^1}\|\rho\|^2_{\mathbb{L}^2(\mathbb{T}^3)}dt. \end{aligned}\\ \label{E4-5'} \begin{aligned} |I_4|&=\int_{\mathbb{Q}}\mu_0\mathbf{M}_{\epsilon}\cdot \nabla\mathbf{H}_{\epsilon}\cdot\mathcal{A}_{i}[\rho_{\epsilon}]dx \\ &\leq \mu_0\|\mathbf{M}_{\epsilon}\cdot\nabla\mathbf{H}_{\epsilon}\|_{\mathbb{L}^1 (\mathbb{S}^1;\mathbb{L}^{3/2}(\mathbb{T}^3))} \|\mathcal{A}_{i}[\rho_{\epsilon}]\|_{\mathbb{L}^{\infty} (\mathbb{S}^1;\mathbb{L}^3{(\mathbb{T}^3)})} \\ &\leq \mu_0\|\rho_{\epsilon}\|_{\mathbb{L}^{\infty} (\mathbb{S}^1;\mathbb{L}^5{(\mathbb{T}^3)})} \|\mathbf{M}_{\epsilon}\|_{\mathbb{L}^2 (\mathbb{S}^1;\mathcal{M})}\|\mathbf{H}_{\epsilon}\|_{\mathbb{L}^2 (\mathbb{S}^1;\mathcal{M})} \\ &\leq C\|\rho_{\epsilon}\|_{\mathbb{L}^{\infty} (\mathbb{S}^1;\mathbb{L}^5{(\mathbb{T}^3)})}, \end{aligned}\\ \label{E4-5R} \begin{aligned} |I_5|&= \zeta\int_{\mathbb{Q}}(\nabla\times \mathbf{u}_{\epsilon}) \cdot(\nabla\times\mathcal{A}_{i}[\rho_{\epsilon}])dx -2\zeta\int_{\mathbb{Q}}(\nabla\times \Omega_{\epsilon}) \cdot\mathcal{A}_{i}[\rho_{\epsilon}]dx \\ &\leq \zeta\|\nabla\times \mathbf{u}_{\epsilon}\|_{\mathbb{L}^2 (\mathbb{S}^1;\mathbb{L}^2(\mathbb{T}^3))}\|\nabla\times\mathcal{A}_{i} [\rho_{\epsilon}]\|_{\mathbb{L}^2(\mathbb{S}^1;\mathbb{L}^2(\mathbb{T}^3))} \\ &\quad +2\zeta\|\nabla\times\Omega_{\epsilon}\|_{\mathbb{L}^2 (\mathbb{S}^1;\mathbb{L}^2(\mathbb{T}^3))}\|\mathcal{A}_{i} [\rho_{\epsilon}]\|_{\mathbb{L}^2(\mathbb{S}^1;\mathbb{L}^2(\mathbb{T}^3))} \\ &\leq C\|\mathbf{u}_{\epsilon}\|_{\mathbb{L}^2(\mathbb{S}^1;\mathbb{H}^{1} (\mathbb{T}^3))}\|\mathcal{A}_{i}[\rho_{\epsilon}]\|_{\mathbb{H}^1 (\mathbb{S}^1;\mathbb{L}^2(\mathbb{T}^3))}\\ &\quad +C\|\Omega_{\epsilon}\|_{\mathbb{L}^2(\mathbb{S}^1;\mathbb{H}^1 (\mathbb{T}^3))}\|\mathcal{A}_{i}[\rho_{\epsilon}]\|_{\mathbb{L}^2 (\mathbb{S}^1;\mathbb{L}^2(\mathbb{T}^3))} \\ &\leq C\|\rho_{\epsilon}\|_{\mathbb{L}^{\infty}(\mathbb{S}^1;\mathbb{L}^4 {(\mathbb{T}^3)})}(\|\mathbf{u}_{\epsilon}\|_{\mathbb{L}^2 (\mathbb{S}^1;\mathbb{H}^{1}(\mathbb{T}^3))}+ \|\Omega_{\epsilon}\|_{\mathbb{L}^2(\mathbb{S}^1;\mathbb{H}^1(\mathbb{T}^3))}). \end{aligned} \end{gather} Thus, by \eqref{E4-1'}-\eqref{E4-5R} and the estimate about $\int_{\mathbb{S}^1}\|\mathbf{u}\|^2_{\mathbb{W}^{1,2}(\mathbb{T}^3)}$, $\int_{\mathbb{S}^1}\|\mathbf{M}\|^2_{\mathbb{W}^{1,2}(\mathbb{T}^3)}$ and $\int_{\mathbb{S}^1}\|\mathbf{H}\|^2_{\mathbb{W}^{1,2}(\mathbb{T}^3)}$, we have \[ \int_{\mathbb{Q}}a\rho^{\gamma+1}+\delta\rho^{\beta+1} +\frac{\mu_0}{2}\rho_{\epsilon}|\mathbf{M}_{\epsilon}|^2\,dx\,dt \leq C\Big(1+\|\rho\|_{\mathbb{L}^{\infty}(\mathbb{S}^1;\mathbb{L}^{\beta} (\mathbb{T}^3))}\Big)^4. \] This combining with Lemma \ref{lem6} implies \[ \int_{\mathbb{S}^1}\|\rho_{\epsilon}\|^{\beta+1}_{\mathbb{L}^{\beta+1} (\mathbb{Q})}dt\leq C,\quad \|E_{\delta}[\rho_{\epsilon},\mathbf{u}_{\epsilon},\Omega_{\epsilon}, \mathbf{M}_{\epsilon},\mathbf{H}_{\epsilon}]\|_{\mathbb{C}(\mathbb{S}^1)}\leq C, \] where $C$ is independent of $\epsilon$. \end{proof} In fact, from Lemma \ref{lem7}, using the standard Sobolev compact embedding and the Arzel\`{a}-Ascoli theorem, we can obtain the following convergence results about $(\rho_{\epsilon},\mathbf{u}_{\epsilon},\Omega_{\epsilon},\mathbf{M}_{\epsilon},\mathbf{H}_{\epsilon})$ as $\epsilon\to 0$. \begin{gather}\label{E4-13} \rho_{\epsilon}\to \rho\quad \text{in }\mathbb{C} (\mathbb{S}^1;\mathbb{L}^{\beta}_{\rm weak}(\mathbb{T}^3)), \\ \label{E4-13R1} \mathbf{u}^i_{\epsilon}\to \mathbf{u}^i,\quad \Omega^i_{\epsilon}\to \Omega^i,\quad \text{weakly in }\mathbb{L}^2(\mathbb{S}^1;\mathbb{W}^{1,2}(\mathbb{T}^3)),\; i=1,2,3, \\ \label{E4-13R2} \mathbf{M}_{\epsilon}\to \mathbf{M},\quad \text{weakly* in } \mathbb{L}^{\infty}(\mathbb{S}^1;\mathbb{L}^2(\mathbb{T}^3)) \cap\mathbb{L}^2(\mathbb{S}^1;\mathcal{M}),\\ \label{E4-13RR} \psi_{\epsilon}\to \psi\quad \text{weakly* in } \mathbb{L}^{\infty}(\mathbb{S}^1;\mathbb{H}^1(\mathbb{T}^3)) \cap\mathbb{L}^2(\mathbb{S}^1;\mathbb{H}^2(\mathbb{T}^3)),\\ \label{E4-13R21} \mathbf{H}_{\epsilon}\to \mathbf{H},\quad \text{weakly* in } \mathbb{L}^{\infty}(\mathbb{S}^1;\mathbb{L}^2(\mathbb{T}^3)) \cap\mathbb{L}^2(\mathbb{S}^1;\mathcal{M}),\\ \label{E4-13R3} \rho_{\epsilon}\mathbf{u}^i_{\epsilon}\to \rho\mathbf{u}^i,\quad \rho_{\epsilon}\Omega^i_{\epsilon}\to \rho\Omega^i,\quad \text{weakly in }\mathbb{L}^2(\mathbb{S}^1;\mathbb{L}^{\frac{6\beta}{\beta+6}} (\mathbb{T}^3)), \\ \label{E4-14} \rho_{\epsilon}\mathbf{u}^i_{\epsilon}\to \rho\mathbf{u}^i,\quad \rho_{\epsilon}\Omega^i_{\epsilon}\to \rho\Omega^i,\quad \text{weakly in } \mathbb{C}(\mathbb{S}^1;\mathbb{L}_{\rm weak}^{\frac{2\beta}{\beta+1}} (\mathbb{T}^3)),\quad i=1,2,3, \\ \label{E4-13R4} \rho_{\epsilon}\mathbf{u}^i_{\epsilon}\mathbf{u}^j_{\epsilon}\to \rho\mathbf{u}^i \mathbf{u}^j,\quad \text{weakly in } \mathbb{L}^2(\mathbb{S}^1;\mathbb{L}^{\frac{6\beta}{4\beta+3}}(\mathbb{T}^3))\quad i,j=1,2,3,\\ \label{E4-13RR1} \rho_{\epsilon}\mathbf{u}^i_{\epsilon}\Omega^j_{\epsilon}\to \rho\mathbf{u}^i\Omega^j,\quad \text{weakly in } \mathbb{L}^2(\mathbb{S}^1;\mathbb{L}^{\frac{6\beta}{4\beta+3}}(\mathbb{T}^3))\quad i,j=1,2,3. \end{gather} Letting $\epsilon\to 0$, we obtain the time periodic solution $(\rho_{\delta},\mathbf{u}_{\delta},\Omega_{\delta}, \mathbf{M}_{\delta}, \mathbf{H}_{\delta})$ of these equations in $\mathcal{D}'(\mathbb{Q})$, \begin{gather}\label{E4-R1} \partial_t\rho+\operatorname{div}(\rho\mathbf{u})=0, \\ \label{E4-R2} \partial_t(\rho u^i)+\operatorname{div}(\rho u^i\otimes \mathbf{u})-\mu\Delta u^i+\mu\partial_{x_i}(\operatorname{div}\mathbf{u}) -\partial_{x_i}\overline{\mathbb{F}}=R+f^i\rho, \\ \label{E4-RRR1} \partial_t(\rho\Omega^i)+\operatorname{div}(\rho\Omega^i\otimes\mathbf{u} )-\mu'\Delta\Omega^i-(\lambda'+\mu')\nabla(\operatorname{div}\Omega) =S+\rho g^i(t,x), \\ \label{E4-R3} \partial_t\mathbf{M}+\operatorname{div}(\mathbf{u}\otimes \mathbf{M})-\sigma\Delta\mathbf{M}+\frac{1}{\tau}(\mathbf{M} -\chi_0\mathbf{H})=\Omega\times\mathbf{M}+\rho\mathbf{l}(t,x), \\ \label{E4-RRR2} \mathbf{H}=\nabla\psi,\quad \operatorname{div}(\mathbf{H}+\mathbf{M}) =\mathbb{F},\text{in }\mathcal{D}'(\mathbb{Q}), \end{gather} where $i=1,2,3$, \begin{gather*} R=\mu_0\mathbf{M}\cdot\nabla\mathbf{H}-\zeta\nabla\times (\nabla\times u^i-2\Omega), \\ S=\mu_0\mathbf{M}\times\mathbf{H}+2\zeta(\nabla\times\mathbf{u}-2\Omega^i), \\ F_{\epsilon}=(2\mu+\lambda)\operatorname{div}\mathbf{u}_{\epsilon} +\frac{\mu_0}{2}|\mathbf{M}_0|^2-a\rho_{\epsilon}^{\gamma} -\delta\rho^{\beta}_{\epsilon},\\ F_{\epsilon}\to \overline{F}\quad \text{weakly in } \mathbb{L}^{\frac{\beta+1}{\beta}}(\mathbb{Q}). \end{gather*} Now, our attention turns to proving the strong convergence of $\rho$ by means of the classical Minty trick. This follows the idea of Lions\cite{Lions5} , Feireisl et al.\cite{Fei1} and Amirat et al.\cite{Ami}. Define the bilinear operator \begin{equation}\label{E4-11} \mathcal{Q}_{i,j}[u,v]=u\mathcal{R}_{i,j}[v]-v\mathcal{R}_{i,j}[u]. \end{equation} According to the work of Coifman et al \cite{Coi}, the bilinear operator $\mathcal{Q}_{i,j}$ is weakly continuous. Applying the discrete Fourier transform to \eqref{E4-11}, for $\mathbf{k}\neq0$ and $a_{\mathbf{0}}=0$, we have \[ a_{\mathbf{k}}[\mathcal{Q}_{i,j}[u,v]] =\sum_{\mathbf{n}\in\mathbb{Z}^3/\{\mathbf{0}\}} \Big(\frac{(k_i-n_i)(k_j-n_j)}{|\mathbf{k}-\mathbf{n}|^2} -\frac{n_in_j}{|\mathbf{n}|^2}\Big)a_{\mathbf{k}-\mathbf{n}}[u]a_{\mathbf{n}}[v]. \] Moreover, if we take $u,v\in\mathbb{L}^2(\mathbb{T}^3)$ such that $u_{\epsilon}\to u$ and $v_{\epsilon}\to v$ weakly in $\mathbb{L}^2(\mathbb{T}^3)$, the operator $\mathcal{Q}_{i,j}$ has the property \begin{equation}\label{E4-15} \mathcal{Q}_{i,j}[u_{\epsilon},v_{\epsilon}]\to \mathcal{Q}_{i,j}[u,v],\quad \text{in }\mathcal{D}'(\mathbb{T}^3). \end{equation} Following the idea of Lions\cite{Lions5}, letting $\epsilon\to 0$, using \eqref{E4-6}-\eqref{E4-5R}, we obtain \begin{equation}\label{E4-12} \begin{aligned} &\limsup_{\epsilon\to 0}\int_{\mathbb{Q}}a\rho^{\gamma+1} +\frac{\mu_0}{2}|\mathbf{M}|^2+\delta\rho^{\beta+1}\,dx\,dt \\ &\leq \int_{\mathbb{Q}}\frac{m}{|\mathbb{T}^3|} \overline{(a\rho^{\gamma}+\delta\rho^{\beta})} +\overline{(\mathbf{u}^i\mathcal{Q}_{i,j}[\rho,\rho\mathbf{u}^j])} +\rho f^{i}\mathcal{A}_i[\rho]+R\mathcal{A}_{i}[\rho]\,dx\,dt. \end{aligned} \end{equation} Moreover, by \eqref{E5-14}, we have \[ \lim_{\epsilon\to 0}m_{\epsilon} =\lim_{\epsilon\to 0}\int_{\mathbb{T}^3}\rho_{\epsilon}dx =\lim_{\epsilon\to 0}\int_{\mathbb{T}^3}\rho dx=m. \] \begin{lemma} \label{lem8} Assume that $\rho_{\epsilon}$ and $\rho_{\epsilon}\mathbf{u}^j_{\epsilon}$ satisfy \eqref{E4-13}-\eqref{E4-13R4}. Then \[ \mathcal{Q}_{i,j}[\rho_{\epsilon},\rho_{\epsilon}\mathbf{u}^j_{\epsilon}]\to \mathcal{Q}_{i,j}[\rho,\rho\mathbf{u}^j],\quad weakly~in~\mathbb{L}^{p}(\mathbb{T}^3),\quad \forall t\in\mathbb{S}^1, \] where $15$, then \[ \mathcal{Q}_{i,j}[\rho_{\epsilon},\rho_{\epsilon}\mathbf{u}^j_{\epsilon}]\to \mathcal{Q}_{i,j}[\rho,\rho\mathbf{u}^j],\quad strongly~in~\mathbb{W}^{-1,2}(\mathbb{T}^3),\quad \forall t\in\mathbb{S}^1. \] \end{lemma} \begin{proof} Denote the cut-off function \begin{equation}\label{E5-13} \mathcal{T}_{k}[x]=\min\{|x|,k\}sign(x). \end{equation} Then \begin{gather*} \|\mathcal{T}_{k}[\rho_{\epsilon}\mathbf{u}^j_{\epsilon}] -\rho_{\epsilon}\mathbf{u}^j_{\epsilon}\|_{\mathbb{L}^{p_1}} \leq k^{p_1-\frac{2\beta}{\beta+1}}\sup_{\epsilon}\|\rho_{\epsilon} \mathbf{u}^j_{\epsilon}\|^{\frac{2\beta}{\beta+1}}_{\mathbb{L} ^{\frac{2\beta}{\beta+1}}(\mathbb{T}^3)},\\ \|\overline{\mathcal{T}_{k}[\rho_{\epsilon}\mathbf{u}^j_{\epsilon}]} -\rho_{\epsilon}\mathbf{u}^j_{\epsilon}\|_{\mathbb{L}^{p_1}} \leq k^{p_1-\frac{2\beta}{\beta+1}}\sup_{\epsilon}\|\rho_{\epsilon} \mathbf{u}^j_{\epsilon}\|^{\frac{2\beta}{\beta+1}}_{\mathbb{L} ^{\frac{2\beta}{\beta+1}}(\mathbb{T}^3)}, \end{gather*} where $p_1=\beta/(\beta-1)$ and $\overline{\mathcal{T}_{k}[\rho_{\epsilon}\mathbf{u}^j_{\epsilon}]}$ denotes a weak limit of $\mathcal{T}_{k}[\rho_{\epsilon}\mathbf{u}^j_{\epsilon}]$. Note that \begin{align*} \mathcal{Q}_{i,j}[\rho_{\epsilon},\rho_{\epsilon}\mathbf{u}^j_{\epsilon}] -\mathcal{Q}_{i,j}[\rho,\rho_{\epsilon}\mathbf{u}^j] &=\mathcal{Q}_{i,j}[\rho_{\epsilon},\mathcal{T}_{k}[\rho_{\epsilon} \mathbf{u}^j_{\epsilon}]]-\mathcal{Q}_{i,j}[\rho,\mathcal{T}^*_{k} [\rho_{\epsilon}\mathbf{u}^j_{\epsilon}]]\\ &\quad +\mathcal{Q}_{i,j}[\rho_{\epsilon},\rho_{\epsilon}\mathbf{u}^j_{\epsilon} -\mathcal{T}_{k}[\rho_{\epsilon}\mathbf{u}^j_{\epsilon}]] +\mathcal{Q}_{i,j}[\rho_{\epsilon},\mathcal{T}_{k}^*[\rho\mathbf{u}^j] -\rho\mathbf{u}^j] \end{align*} Hence, for a sufficiently large $k$, by virtue of \eqref{E4-15}, we can obtain the result. This completes the proof. \end{proof} Multiplying both sides of \eqref{E4-R2} by $\mathcal{A}_i[\rho]$, we obtain \begin{equation}\label{E4-16} \begin{aligned} \int_{\mathbb{Q}}\overline{(a\rho^{\gamma}+\delta\rho^{\beta})}\rho \,dx\,dt &=\int_{\mathbb{Q}}\frac{m}{|\mathbb{T}^3|}\overline{(a\rho^{\gamma} +\delta\rho^{\beta})} +(\lambda+2\mu)\rho \operatorname{div}\mathbf{u}\\ &\quad +\mathbf{u}^i\mathcal{Q}_{i,j}[\rho,\rho\mathbf{u}^j] +\rho f^{i}\mathcal{A}_i[\rho]+R\mathcal{A}_{i}[\rho]\,dx\,dt, \end{aligned} \end{equation} where we use \[ \int_{\mathbb{Q}}\overline{(\mathbf{u}^i\mathcal{Q}_{i,j} [\rho,\rho\mathbf{u}^j])}\,dx\,dt =\int_{\mathbb{Q}}\mathbf{u}^i\mathcal{Q}_{i,j}[\rho,\rho\mathbf{u}^j]\,dx\,dt. \] By \eqref{E4-12} and \eqref{E4-16}, we have \[ \limsup_{\epsilon\to 0}\int_{\mathbb{Q}}a\rho_{\epsilon}^{\gamma+1} +\delta\rho_{\epsilon}^{\beta+1}\,dx\,dt \leq\int_{\mathbb{Q}}\rho\overline{(a\rho^{\gamma} +\delta\rho^{\beta})}-(\lambda+2\mu)\rho(\operatorname{div}\mathbf{u})\,dx\,dt, \] which implies the strong convergence of $\rho$ by means of the classical Minty trick. Next, we will follow \cite{Duv} and \cite{Fei1} to prove that $\rho$ and $\mathbf{u}$ are a renormalized time-periodic solution of \eqref{E4-R1}. Let $\theta_k=k^3\theta(kx)$ be the regularizing kernels with $\theta\in\mathbb{C}_{0}^{\infty}(\mathbb{T}^3)$ and $\int_{\mathbb{T}^3}\theta dx=1$, for $k>0$. Then $\rho_k=\theta_k\ast\rho$ is a regular solution of \begin{equation}\label{E4-17} \partial_t\rho_k+\operatorname{div}(\rho_k\mathbf{u}) =\operatorname{div}(\rho_k\mathbf{u})-\theta_k\ast(\operatorname{div} (\rho\mathbf{u}))\quad\text{on }\mathbb{Q}. \end{equation} Note that $\rho,\operatorname{div}\mathbf{u}\in\mathbb{L}^2(\mathbb{Q})$. Using \cite[Lemma 2.3]{Lions1}, we obtain that \[ \operatorname{div}(\rho_k\mathbf{u})-\theta_k\ast(\operatorname{div} (\rho\mathbf{u}))\to 0,\quad\text{as }k\to \infty. \] Multiplying \eqref{E4-17} by $b'(\rho_k)$ and passing to the limit, we obtain \[ \frac{\partial b(\rho)}{\partial t}+\operatorname{div} (b(\rho)\mathbf{u})+(b'(\rho)\rho-b(\rho))\operatorname{div}\mathbf{u}=0. \] Furthermore, taking $b'(\rho_k)=\rho_k\log\rho_k$ and integrating above equation over $\mathbb{Q}$, we obtain \[ \int_{\mathbb{Q}}\rho(\operatorname{div}\mathbf{u})\,dx\,dt=0. \] This shows that $\rho$ and $\mathbf{u}$ are a renormalized time-periodic solution of \eqref{E4-R1}. \section{Passing to the limit as $\delta\to 0$} In this section, our target is to pass to the limit as $\delta\to 0$ in \eqref{E4-RR1}-\eqref{E4-RRRR2}. This process of proof is similar to the proof of passing to the limit as $\epsilon\to 0$. Firstly, we need to show that \begin{gather*} \|E_{\delta}[\rho_{\delta},\mathbf{u}_{\delta},\Omega_{\delta}, \mathbf{M}_{\delta},\mathbf{H}_{\delta}]\|_{\mathbb{C}(\mathbb{S}^1)},\quad \|\mathbf{u}_{\delta}\|_{\mathbb{L}^2(\mathbb{S}^1; \mathbb{W}^{1,2}(\mathbb{T}^3))},\quad \|\Omega_{\delta}\|_{\mathbb{L}^2(\mathbb{S}^1;\mathbb{W}^{1,2}(\mathbb{T}^3))}\\ \|\mathbf{M}_{\delta}\|_{\mathbb{L}^2(\mathbb{S}^1;\mathcal{M}(\mathbb{T}^3))}, \quad \|\mathbf{H}_{\delta}\|_{\mathbb{L}^2(\mathbb{S}^1; \mathcal{M}(\mathbb{T}^3))},\quad \int_{\mathbb{Q}}a\rho^{\gamma+\vartheta}+\delta\rho^{\beta+\vartheta}\,dx\,dt \end{gather*} are bounded independently of $\delta$, for some $\vartheta>0$. Multiplying both side of \eqref{E4-RR1} by $\vartheta\rho^{\vartheta-1}$, for some $\vartheta>0$, we obtain \begin{equation}\label{E4-20} \partial_t\rho^{\vartheta}+\operatorname{div}(\rho^{\vartheta}\mathbf{u}) +(\vartheta-1)\rho^{\vartheta}\operatorname{div}\mathbf{u}=0,\quad \text{in }\mathcal{D}'(\mathbb{Q}). \end{equation} As in \cite{Fei1}, we take $\vartheta=1/5$ and a test function $\mathcal{A}_i[\rho^{\vartheta}]$ of \eqref{E4-RR2}. Then by \eqref{E4-1}, \eqref{E4-18} and \eqref{E4-20}, we obtain \begin{equation}\label{E5-6} \begin{aligned} &\int_{\mathbb{Q}}a\rho^{\gamma+\vartheta}+\delta\rho^{\beta+\vartheta} +\frac{\mu_0}{2}\rho|\mathbf{M}|^2\,dx\,dt \\ &= \int_{\mathbb{S}^1}(\int_{\mathbb{T}^3}a\rho^{\gamma} +\delta\rho^{\beta})\int_{\mathbb{T}^3}\rho^{\vartheta}\,dx\,dt +\int_{\mathbb{Q}}R\mathcal{A}_i[\rho^{\vartheta}]\,dx\,dt \\ &\quad +(1-\vartheta)\int_{\mathbb{Q}}\rho^{\vartheta} (\operatorname{div}\mathbf{u})\mathcal{A}_i[\rho\mathbf{u}^i] -f^i\mathcal{A}_i[\rho^{\vartheta}]\rho +\mathbf{u}^i\mathcal{Q}_{i,j}[\rho,\rho\mathbf{u}^j]\,dx\,dt. \end{aligned} \end{equation} The following estimates are taken from \cite{Fei1}. For $\gamma\geq 9/5$, \begin{gather}\label{E5-1} \|\mathcal{A}_i[\rho\mathbf{u}^i]\|_{\mathbb{L}^{18/7}(\mathbb{T}^3)} \leq C\|\rho\|_{\mathbb{L}^{\gamma}(\mathbb{T}^3)} \|\mathbf{u}\|_{\mathbb{W}^{1,2}(\mathbb{T}^3)}, \\ \label{E5-2} |(1-\vartheta)\int_{\mathbb{Q}}\rho^{\vartheta}(\operatorname{div}\mathbf{u})\mathcal{A}_i [\rho\mathbf{u}^i]\,dx\,dt| \leq C(1+\int_{\mathbb{S}^1}\|\mathbf{u}\|_{\mathbb{W}^{1,2}(\mathbb{T}^3)} \|\rho\|_{\mathbb{L}^{\gamma}(\mathbb{T}^3)}^{6/5}dt),\\ \label{E5-3} |\int_{\mathbb{Q}}f^i\mathcal{A}_i[\rho^{\vartheta}]\rho \,dx\,dt|\leq C, \\ \label{E5-4} |\int_{\mathbb{Q}}\mathbf{u}^i\mathcal{Q}_{i,j}[\rho,\rho\mathbf{u}^j]\,dx\,dt| \leq C(1+\int_{\mathbb{S}^1}\|\mathbf{u}\|_{\mathbb{W}^{1,2}(\mathbb{T}^3)} \|\rho\|_{\mathbb{L}^{\gamma}(\mathbb{T}^3)}^{6/5}dt). \end{gather} The rest is to estimate $\int_{\mathbb{Q}}R\mathcal{A}_i[\rho^{\vartheta}]\,dx\,dt$. Integrating the energy inequality \eqref{E4-19} over $\mathbb{S}^1$ and by \eqref{E2-4}, we obtain \begin{align*} &\int_{\mathbb{S}^1}(\|\mathbf{u}\|_{\mathbb{W}^{1,2}(\mathbb{T}^3)} + \|\Omega\|_{\mathbb{W}^{1,2}(\mathbb{T}^3)} +\|\mathbf{M}\|_{\mathbb{W}^{1,2}(\mathbb{T}^3)} +\|\mathbf{H}\|_{\mathbb{W}^{1,2}(\mathbb{T}^3)})\,dx\,dt\\ &\leq C(1+\int_{\mathbb{S}^1}\|\sqrt{\rho}\|_{\mathbb{L}^3(\mathbb{T}^3)}( \|\mathbf{u}\|_{\mathbb{L}^6(\mathbb{T}^3)}+\|\Omega\|_{\mathbb{L}^6 (\mathbb{T}^3)}+\|\mathbf{M}\|_{\mathbb{L}^6(\mathbb{T}^3)})\,dx\,dt), \end{align*} which implies \begin{equation}\label{E5-7} \begin{aligned} &\int_{\mathbb{S}^1}(\|\mathbf{u}\|_{\mathbb{W}^{1,2}(\mathbb{T}^3)} +\|\Omega\|_{\mathbb{W}^{1,2}(\mathbb{T}^3)} +\|\mathbf{M}\|_{\mathbb{W}^{1,2}(\mathbb{T}^3)} +\|\mathbf{H}\|_{\mathbb{W}^{1,2}(\mathbb{T}^3)})\,dx\,dt \\ &\leq C(1+\|\rho\|_{\mathbb{C}(\mathbb{S}^1;\mathbb{L}^{3/2}(\mathbb{T}^3))}). \end{aligned} \end{equation} By H\"{o}lder's inequality, the above inequality and \eqref{E4-1}, for $s\geq 3/4$, we have \begin{equation}\label{E5-5} |\int_{\mathbb{Q}}R\mathcal{A}_i[\rho^{\vartheta}]\,dx\,dt|\leq C(1+\|\rho\|_{\mathbb{C}(\mathbb{S}^1;\mathbb{L}^s(\mathbb{T}^3))}). \end{equation} By H\"{o}lder's inequality and the standard Sobolev theorem, for $\tau>2$, it follows from \eqref{E5-6}-\eqref{E5-5} that \[ \int_{\mathbb{Q}}a\rho^{\gamma+\vartheta}+\delta\rho^{\beta+\vartheta}\,dx\,dt \leq C(1+\|\rho\|_{\mathbb{C}(\mathbb{S}^1;\mathbb{L}^{\gamma} (\mathbb{T}^3))})^\tau. \] Furthermore, using the method of estimating \eqref{E5-6}, the boundness independent of $\delta$ of the values $\|E_{\delta}[\rho_{\delta},\mathbf{u}_{\delta},\Omega_{\delta}, \mathbf{M}_{\delta},\mathbf{H}_{\delta}]\|_{\mathbb{L}^{\infty}(\mathbb{S}^1)}$, $\|\mathbf{u}_{\delta}\|_{\mathbb{L}^2(\mathbb{S}^1;\mathbb{W}^{1,2} (\mathbb{T}^3))}$, $\|\Omega_{\delta}\|_{\mathbb{L}^2(\mathbb{S}^1; \mathbb{W}^{1,2}(\mathbb{T}^3))}$, $\|\mathbf{M}_{\delta}\|_{\mathbb{L}^2(\mathbb{S}^1;\mathcal{M}(\mathbb{T}^3))}$ and $\|\mathbf{H}_{\delta}\|_{\mathbb{L}^2(\mathbb{S}^1;\mathcal{M}(\mathbb{T}^3))}$ can be obtained. Letting $\delta\to 0$ in \eqref{E4-RR1}-\eqref{E4-RRRR2}, we obtain the time periodic solution $(\rho,\mathbf{u},\Omega,\mathbf{M},\mathbf{H})$ of the equations in $\mathcal{D}'(\mathbb{Q})$, \begin{gather}\label{E5-R1} \partial_t\rho+\operatorname{div}(\rho\mathbf{u})=0, \\ \label{E5-R2} \partial_t(\rho u^i)+\operatorname{div}(\rho u^i\otimes \mathbf{u})-\mu\Delta u^i+\mu\partial_{x_i}(\operatorname{div}\mathbf{u})-\partial_{x_i}\overline{\mathbb{F}} =R+\rho f^i(t,x), \\ \label{E5-R3} \partial_t(\rho\Omega^i)+\operatorname{div}(\rho\Omega^i\otimes\mathbf{u} )-\mu'\Delta\Omega^i-(\lambda'+\mu')\nabla(\operatorname{div}\Omega) =S+\rho g^i(t,x), \\ \label{E5-RR1} \partial_t\mathbf{M}+\operatorname{div}(\mathbf{u}\otimes \mathbf{M})-\sigma\Delta\mathbf{M}+\frac{1}{\tau}(\mathbf{M} -\chi_0\mathbf{H})=\Omega\times\mathbf{M}+\rho\mathbf{l}(t,x), \\ \label{E5-RR2} \mathbf{H}=\nabla\psi,\quad \operatorname{div}(\mathbf{H}+\mathbf{M})=\mathbb{F}, \end{gather} where $i=1,2,3$, \begin{gather*} R=\mu_0\mathbf{M}\cdot\nabla\mathbf{H}-\zeta\nabla\times (\nabla\times u^i-2\Omega), \\ S=\mu_0\mathbf{M}\times\mathbf{H}+2\zeta(\nabla\times\mathbf{u}-2\Omega^i), \\ F_{\epsilon}=(2\mu+\lambda)\operatorname{div}\mathbf{u}_{\epsilon} +\frac{\mu_0}{2}|\mathbf{M}_{\epsilon}|^2-a\rho_{\epsilon}^{\gamma} -\delta\rho^{\beta}_{\epsilon}, \\ \overline{F}=(2\mu+\lambda)\operatorname{div}\mathbf{u} +\frac{\mu_0}{2}|\mathbf{M}|^2-\overline{p(\rho)},\\ F_{\epsilon}\to \overline{F}\quad \text{weakly in } \mathbb{L}^{\frac{\gamma+\vartheta}{\gamma}}(\mathbb{Q}). \end{gather*} Furthermore, using the same process of deducing that $\rho$ is a renormalized solution of the continuous equation \eqref{E5-R1}, for sublinear continuous function $h$, it holds \[ \int_{\mathbb{Q}}h(\rho)\operatorname{div}\mathbf{u}\,dx\,dt=0. \] To prove our main result, Theorem \ref{thm1}, it remains to show the density $\rho_{\delta}\to \rho$ strongly on $\mathbb{Q}$. We introduce a uniformly bounded smooth function $h(x)\in\mathbb{C}^1(\mathbb{R})$ with $h(0)=0$ and uniformly bounded $h'(x)$ and a sequence of functions $\mathcal{H}_k=h(x)-\mathcal{I}_k(x-h(x))$, where $\mathcal{I}_k(x)\in\mathbb{C}^{\infty}(\mathbb{R}^+)$ is bounded nondecreasing function with $0\leq\mathcal{I}_k(x)\leq x$, $\mathcal{I}_k(x)=x$ on $[0,k]$ and $\mathcal{I}_{k+1}(x)\geq\mathcal{I}_k(x)$. According to the definition of $\mathcal{H}_k$, it is easy to check that \begin{equation}\label{E5-12} \int_{\mathbb{Q}}|\mathcal{H}_k(\rho_{\delta})-\rho_{\delta}|^2\,dx\,dt \leq k^{-\zeta}\int_{\mathbb{Q}}|\rho_{\delta}|^{\gamma+\vartheta}\,dx\,dt,\quad \zeta=\gamma+\vartheta-2>0. \end{equation} Multiplying both side \eqref{E4-RR2} by $h(\rho_{\delta})$ and using \eqref{E1-R} and \eqref{E4-18}, we obtain \begin{equation}\label{E5-9} \begin{aligned} &\int_{\mathbb{Q}}a\rho_{\delta}^{\gamma}h(\rho_{\delta}) +\delta\rho_{\delta}^{\beta}h(\rho_{\delta})\,dx\,dt \\ &= \int_{\mathbb{S}^1}(\int_{\mathbb{T}^3}a\rho_{\delta}^{\gamma} +\delta\rho_{\delta}^{\beta})\int_{\mathbb{T}^3}h(\rho_{\delta})\,dx\,dt +\int_{\mathbb{Q}}R\mathcal{A}_i[h(\rho_{\delta})]\,dx\,dt \\ &\quad +\int_{\mathbb{Q}}(h(\rho_{\delta})-\rho_{\delta}h'(\rho_{\delta})) (\operatorname{div}\mathbf{u}_{\delta})\mathcal{A}_i [\rho_{\delta}\mathbf{u}^i_{\delta}]-f^i\mathcal{A}_i[h(\rho_{\delta})] \rho_{\delta}\,dx\,dt \\ &\quad +\int_{\mathbb{Q}}\mathbf{u}^i\mathcal{Q}_{i,j} [h(\rho_{\delta}),\rho_{\delta}\mathbf{u}_{\delta}^j]\,dx\,dt. \end{aligned} \end{equation} Note that $h(\rho_{\delta})$ is bounded and satisfies \eqref{E1-R}. Letting $\delta\to 0$ implies that \begin{gather}\label{E5-8} h(\rho_{\delta})\to \overline{h(\rho)}\quad \text{in }\mathbb{C}(\mathbb{S}^1;\mathbb{L}^q(\mathbb{T}^3)),\quad \forall q\in(1,\infty), \\ \label{E5-10} \partial_t(\overline{h(\rho)})+\operatorname{div} (\overline{h(\rho)}\mathbf{u})+\overline{(h'(\rho)\rho-h(\rho))} \operatorname{div}\mathbf{u}=0. \end{gather} Combining \eqref{E4-13}-\eqref{E4-13R4}, Lemma \ref{lem8} and \eqref{E5-8}, we pass to the limit in \eqref{E5-9} to obtain \begin{align*} \int_{\mathbb{Q}}a\overline{(\rho^{\gamma}h(\rho))}\,dx\,dt &= \int_{\mathbb{S}^1}(\int_{\mathbb{T}^3}a\overline{\rho^{\gamma}}dx) \Big(\int_{\mathbb{T}^3}\overline{h(\rho)}dx\Big)dt +\int_{\mathbb{Q}}R\mathcal{A}_i[\overline{h(\rho)}]\,dx\,dt \\ &\quad +\int_{\mathbb{Q}}\overline{(h(\rho)-\rho h'(\rho))\operatorname{div} \mathbf{u}}\mathcal{A}_i[\rho\mathbf{u}^i]-f^i\mathcal{A}_i [\overline{h(\rho)}]\rho \\ &\quad +\mathbf{u}^i\mathcal{Q}_{i,j}[\overline{h(\rho)},\rho\mathbf{u}^j ]\,dx\,dt. \end{align*} Multiplying both sides of \eqref{E5-R2} by $\mathcal{A}_i(\overline{h(\rho)})$ and using \eqref{E5-10}, we have \begin{align*} \int_{\mathbb{Q}}a\overline{\rho^{\gamma}}\overline{h(\rho))}\,dx\,dt &=\int_{\mathbb{S}^1}(\int_{\mathbb{T}^3}a\overline{\rho^{\gamma}}dx) \Big(\int_{\mathbb{T}^3}\overline{h(\rho)}dx\Big)dt +\int_{\mathbb{Q}}R\mathcal{A}_i[\overline{h(\rho)}]\,dx\,dt \\ &\quad +\int_{\mathbb{Q}}\overline{h(\rho)-\rho h'(\rho))\operatorname{div} \mathbf{u}}\mathcal{A}_i[\rho\mathbf{u}^i] -f^i\mathcal{A}_i[\overline{h(\rho)}]\rho \,dx\,dt \\ &\quad +\int_{\mathbb{Q}}\mathbf{u}^i\mathcal{Q}_{i,j} [\overline{h(\rho)},\rho\mathbf{u}^j]\,dx\,dt. \end{align*} From the two equalities above, it follows that \begin{equation}\label{E5-11} \int_{\mathbb{Q}}\overline{\rho^{\gamma}h(\rho))}\,dx\,dt =\int_{\mathbb{Q}}\overline{\rho^{\gamma}}\overline{h(\rho)}\,dx\,dt -\frac{2\mu+\lambda}{a}\int_{\mathbb{Q}}(\overline{h(\rho)}-\rho) \operatorname{div}\mathbf{u}\,dx\,dt. \end{equation} Note that $h$ is nondecreasing and $\mathcal{H}_k$ is increasing about $k$. So, by \eqref{E5-11} and \eqref{E5-12}, it follows that \begin{align*} 0&\leq\int_{\mathbb{Q}}(\overline{\rho^{\gamma}h(\rho))} -\overline{\rho^{\gamma}}\overline{h(\rho)}\,dx\,dt\\ &\leq \int_{\mathbb{Q}}\overline{\rho^{\gamma}\mathcal{H}_k(\rho)} -\overline{\rho^{\gamma}}\overline{\mathcal{H}_k(\rho)}\,dx\,dt\\ &\leq -\frac{2\mu+\lambda}{a}\int_{\mathbb{Q}} \overline{\mathcal{H}_k(\rho)}-\rho)\operatorname{div}\mathbf{u}\,dx\,dt\\ &\leq Ck^{-\zeta}. \end{align*} Hence, for increasing $k$, we obtain \[ \int_{\mathbb{Q}}\overline{\rho^{\gamma}h(\rho)}\,dx\,dt =\int_{\mathbb{Q}}\overline{\rho^{\gamma}}\overline{h(\rho)}\,dx\,dt, \] which implies that the density $\rho_{\delta}\to \rho$ strongly on $\mathbb{Q}$. \subsection*{Acknowledgements} The author is supported by grants: 11201172 from the NSFC, 20120061120002 from SRFDP, and project 985 from Jilin University. \begin{thebibliography}{00} \bibitem{Ami1} Amirat, Y., Hamdache, K., Muiat, F.: Global weak solutions to the equations of motion for magnetic fuilds. \textit{J. Math. Fluid. Mech} \textbf{10}, 326-351 (2008) \bibitem{Ami} Amirat, Y., Hamdache, K.: Weak solutions to the equations of motion for compressible magnetic fuilds. \textit{J. Math. Pure. Appl} \textbf{91}, 433-467 (2009) \bibitem{Coi} Coifman, R., Lions, P.-L., Meyer, Y., Semmes, S.: Compensated compactness and Hardy spaces. \textit{J. Math. Pure. Appl} \textbf{72}, 247-286 (1993) \bibitem{Lions1} Diperna, R.-J., Lions, P.-L.: Ordinary differential equations, transport theory and Sobolev spaces. \textit{Invent. Math} \textbf{98}, 511-547 (1989) \bibitem{Dau} Dautray, R., Lions, P.L.: Analyse math\'{e}matique et calcul num\'{e}rique pour les sciences et les techniques, tome 5, Masson, 1984 \bibitem{Fei4} Ducomet, B., Feireisl, E.: The equations of Magnetohydrodynamics: on the interaction between matter and radiation in the evolution of gaseous stars. \textit{Commun. Math. Phys.} \textbf{226}, 595-629 (2006) \bibitem{Duv} Duvaut, G., Lions, J. L.: In\'{e}quation en thermo\'{e}lasticit\'{e} et magn\'{e}to-hydrodynamique. \textit{Arch. Rational Mech. Anal.} \textbf{46}, 241-279 (1972) \bibitem{Fei1} Feireisl, E., Ne\u{c}asov\'a, S.-M., Petzeltov\'a, H., Stra\u{s}kraba, I.: On the motion of a viscous compressible fluid driven by a time periodic external force. \textit{Arch. Rational Mech. Anal} \textbf{149}, 69-96 (1999) \bibitem{Fei2} Feireisl, E., Petzeltov\'a, H.: Large time behaviour of solutions to the Navier-Stokes equations of compressible flow. \textit{Arch. Rational Mech. Anal} \textbf{150}, 77-96 (1999) \bibitem{Fei3} Feireisl, E.: Dynamics of viscous compressible fluids. Oxford Lecture Series in Mathematics and its applications, vol. 26. Oxford University Press, Oxford, 2004 \bibitem{Gas} Gaspari, G. D.: Bloch equation for conduction electron spin resonance. \textit{Phys. Rev.} \textbf{131}, 215-219 (1966) \bibitem{Grief} Grief, A. D., Richardson, G.: Mathematical modelling of magnetically targeted drug delivery. \textit{J. Magnetism Magnetic Materials.} \textbf{293} (73), 455-463 (2005) \bibitem{Hu1} Hu, X., Wang, D.: Global solutions to the three dimensional full compressible magnetohydrodynamic flows. \textit{Comm. Math. Phys} \textbf{283}, 255-284 (2008) \bibitem{Hu2} Hu, X., Wang, D.: Global existence and large time behavior of solutions to the three dimensional equations of compressible magnetohydrodynamic flows. \textit{Arch. Rational Mech. Anal.} \textbf{197}, 203-238 (2010) \bibitem{Hu3} Hu, X., Wang, D.: Low mach number limit of viscous compressible magnetohydrodynamic flows. \textit{SIAM J. Math. Anal.} \textbf{41}, 1272-1294 (2009) \bibitem{Jiang} Jiang, S., Ju, Q.C., Li, F.C.: Incompressible limit of the compressible Magnetohydrodynamic equations with periodic boundary conditions. \textit{Comm. Math. Phys.} \textbf{297}, 371-400 (2010) \bibitem{Jiang1} Jiang, S., Ju, Q.C., Li, F.C.: Incompressible limit of the compressible Magnetohydrodynamic equations with vanishing viscosity coefficients. \textit{SIAM J. Math. Anal.} \textbf{42}, 2539-2553 (2010) \bibitem{Kuli} Kulikovskiy, A.G., Lifshitz, E.M.: Magnetohydrodynamics. Reading, MA: Addison-Wesley, 1965 \bibitem{Laudau} Laudau, L.D., Lifshitz, E.M.: Electrodynamics of continuous media. 2nd ed., New York: Pergamon, 1984 \bibitem{Lions3} Lions, P.-L.: Mathematical topics in fluid dynamics. vol. 2. Compressible models. Oxford University Press: Oxford, 1998 \bibitem{Lions4} Lions, P.-L.: Quelques M\'{e}thodes de r\'{e}solution des probl\`{e}mes aux limites non lin\'{e}aires, Dunod-Gauthier-Villars, 1969 \bibitem{Lions5} Lions, P.-L: Compacit\'{e} des solutions des \'{e}equation de Navier-Stokes compressible isentropiques. \textit{C. R. Acad. Sci. Paris, S\'{e}r I.} \textbf{316} 1335-1340 (1993) \bibitem{Pankhurst} Pankhurst, Q. Q. A., Connolly, J., Jones, S. K., Dobson, J.: Applications of magnetic nonoparticles in biomedicine. \textit{J. Phys. D: Appl. Phys.} \textbf{36}, R167-R181 (2003) \bibitem{Rosensweig} Rosensweig, R.E., in: S. Odenbache(Ed.), Ferrofluid: Magnetically controllable fluids and their applications, in: Lecture Notes in Physics, vol, \textbf{594}, Springer-Verlag, Heidelberg, 61-84 (2002) \bibitem{Polovin} Polovin, R. V., Demutskii, V. P.: Fundamentals of Magnetohydrodynamics. New York: Consultants, Bureau, 1990 \bibitem{Prouse} Prouse, G.: Soluzioni periodiche dell�equazione di Navier-Stokes. \textit{Atti Accad. Naz. Lincei Rend. Cl. Sci. Fis. Mat. Natur} \textbf{35}, 443-447 (1963) \bibitem{Ser} Sermange, M., Temam, R.: Some mathematical questions related to the MHD equations. \textit{Comm. Pure Appl. Math.} \textbf{36}, 635-664 (1983) \bibitem{Temam} Temam, R.: Navier-Stokes equations, third(revised)ed., Elsevier Science Publishers B.V., Amsterdam, 1984 \bibitem{Torrey} Torrey, H.C.: Bloch equations with diffusion terms. \textit{Phys. Rev} \textbf{104}: 563-565 (1956) \bibitem{Venk} Venkatasubramanian, S., Kaloni, P.: Stability and uniqueness of magnetic fluid motion. \textit{Proc. Roy. Soc. London A} \textbf{458}, 1189-1204 (2002) \bibitem{Yan} Yan, W., Li, Y.: Existence of periodic flows for compressible Magnetohydrodynamics in $\mathbb{T}^3$. Submitted. \bibitem{Zahn} Zahn, M.: Magnetic fluid and nonoparticle applications to nanotechnology. \textit{J. Nanoparticle Res.} \textbf{3} (73), 73-78 (2001) \end{thebibliography} \end{document}