The Algebraic Decay Behavior of Weak Solutions to the Magnetohydrodynamic Equations in Unbounded

1. Introduction

The magnetohydrodynamic (MHD) equations are the core mathematical model describing the interaction between electrically conducting fluids (e.g., plasmas, liquid metals) and electromagnetic fields. They are widely applied in fields such as astrophysics (e.g., solar wind evolution), controlled nuclear fusion (e.g., flow fields in tokamak devices), and industrial fluid engineering [16]. When simulating practical scenarios like cosmic space or infinitely extended pipe flows, the domain containing the flow field must be abstracted as a general smooth unbounded domain in ${\mathbb{R}}^3$. The non-compactness of such domains poses unique challenges to the mathematical analysis of solutions to the equations.

The three-dimensional incompressible viscous MHD equations typically consist of evolution equations for the velocity field $\mathbf{u}=(u_1,u_2,u_3)(x,t)$, magnetic field $\mathbf{B}=(B_1,B_2,B_3)(x,t)$, and pressure $p(x,t)$, along with the fluid incompressibility condition $(\nabla ·\mathbf{u}=0)$ and the magnetic field divergence-free condition $(\nabla ·\mathbf{B}=0)$. To simplify the analysis, this paper sets the fluid Reynolds number $Re$, magnetic Reynolds number $Rm$, and the parameter $S={\displaystyle \frac{M^2}{ReRm}}$(derived from the Hartmann number M) all to 1. Additionally, the pressure term is integrated into p (including the contribution of ${\displaystyle \frac{S}{2}}{\left|B\right|}^2$ ), resulting in the following initial-boundary value problem:

$$\left\{\begin{array}{cc}{\partial }_{t}u-\Delta u+(u·\nabla )u-(B·\nabla )B+\nabla p=0,\hfill & (x,t)\in \Omega \times (0,\infty ),\hfill \\ {\partial }_{t}B-\Delta B+(u·\nabla )B-(B·\nabla )u=0,\hfill & (x,t)\in \Omega \times (0,\infty ),\hfill \\ \nabla ·u=0,\nabla ·B=0,\hfill & (x,t)\in \Omega \times (0,\infty ),\hfill \\ u=B=0,\hfill & (x,t)\in \partial \Omega \times (0,\infty ),\hfill \\ u(x,0)=a,B(x,0)=b,\hfill & x\in \Omega \hfill \end{array}\right.$$

where $\Omega \subset {\mathbb{R}}^3$ is a general smooth unbounded domain, $\partial \Omega$ denotes its boundary, and $a\left(x\right),b\left(x\right)$ are the initial data for the velocity and magnetic field, respectively, which satisfy the divergence-free condition.

In the mathematical theory of the MHD equations, the existence of weak solutions and their long-time behavior are core topics. Earlier, Duvaut and Lions [1] proved that weak solutions corresponding to regular initial data in the two-dimensional case possess regularity; however, the regularity of weak solutions in high dimensions ($n\ge 3$) remains an open problem to this day. Regarding long-time behavior, existing results are mostly limited to special unbounded domains (e.g., half-spaces, exterior domains [5,6,7,8]) and can only prove that solutions tend to zero in $L^2$ as $t\to \infty$. Their exact decay rates and the dominant factors of decay have not been systematically established for general unbounded domains. The key issue lies in the failure of traditional analytical methods (such as Fourier separation techniques and $L^q$-$L^r$ estimates for the Stokes operator semigroup [2]) in general unbounded domains, which prevents the unified control of nonlinear terms (e.g.,$(u·\nabla )u$,$(B·\nabla )B$).

To overcome the above difficulties, this paper adopts a method combining "spectral analysis + approximation theory + energy estimates":

• First, based on the spectral decomposition theory [4] of the self-adjoint Stokes operator $A=-P\Delta$ (where P is the Helmholtz projection operator) and the properties of fractional-order operators, a sequence of Yosida approximate solutions is constructed to avoid the limitations of traditional methods;
• Second, through spectral measure decomposition and interpolation inequalities, unified estimates for all nonlinear terms in the system are established;
• Finally, by combining energy inequalities and the theory of weak convergence, this paper not only proves the existence of global weak solutions corresponding to initial data $a,b\in L_{\sigma }^2(\Omega )$($L_{\sigma }^2(\Omega )$ is the closed subspace of $L^2(\Omega )$ consisting of divergence-free vector fields [3]) but also derives the $L^2$ decay rate of the solutions. It further reveals the core law of their decay behavior: the algebraic decay property of weak solutions is typically dominated by their linear component, i.e., the Stokes operator semigroup solution $e^{-tA}(a,b)$.

2. Definition of Weak Solution

Let $\Omega \subset {\mathbb{R}}^3$ be a general unbounded smooth domain, and $T\in (0,\infty ]$. A pair of vector fields $(u,B)$ is called a weak solution to the MHD system on $\Omega \times (0,T)$ if it satisfies the following conditions:

• $u,B\in L^{\infty }(0,T;L_{\sigma }^2(\Omega ))\cap L^2(0,T;H_0^1{(\Omega )}^3)$, where $L_{\sigma }^2(\Omega )$ denotes the closed subspace of $L^2{(\Omega )}^3$ consisting of all divergence-free vector fields, and $H_0^1{(\Omega )}^3$ is the closure of $C_0^{\infty }{(\Omega )}^3$ in $H^1{(\Omega )}^3$.
• For any test functions $\varphi$,$\psi \in C_0^{\infty }{(\Omega \times [0,T))}^3$ satisfying $\nabla ·\varphi =\nabla ·\psi =0$, and for any $t\in (0,T)$, the following hold:

  $$\begin{array}{cc}\hfill & {\int }_{\Omega }u(x,t)·\varphi (x,t)dx+{\int }_0^t{\int }_{\Omega }\nabla u·\nabla \varphi dxds+{\int }_0^t{\int }_{\Omega }\left[(u·\nabla )u-(B·\nabla )B\right]·\varphi dxds\hfill \\ \hfill & ={\int }_{\Omega }a\left(x\right)·\varphi (x,0)dx\hfill \end{array}$$

  and

  $$\begin{array}{cc}\hfill & {\int }_{\Omega }B(x,t)·\psi (x,t)dx+{\int }_0^t{\int }_{\Omega }\nabla B·\nabla \psi dxds+{\int }_0^t{\int }_{\Omega }\left[(u·\nabla )B-(B·\nabla )u\right]·\psi dxds\hfill \\ \hfill & ={\int }_{\Omega }b\left(x\right)·\psi (x,0)dx.\hfill \end{array}$$
• For almost all $t\in (0,T)$, the following holds:

  $${\parallel u\left(t\right)\parallel }_{L^2}^2+{\parallel B\left(t\right)\parallel }_{L^2}^2+2{\int }_0^t\left({\parallel \nabla u\left(s\right)\parallel }_{L^2}^2+{\parallel \nabla B\left(s\right)\parallel }_{L^2}^2\right)ds\le {\parallel a\parallel }_{L^2}^2+{\parallel b\parallel }_{L^2}^2.$$
  If $T=\infty$, then $(u,B)$ is called a global weak solution.

Theorem 1

(Existence and Decay of Global Weak Solutions). Let $\Omega \subset {\mathbb{R}}^3$ be a general unbounded smooth domain, and let the initial data satisfy $a,b\in L_{\sigma }^2(\Omega )$. Then there exists a global weak solution $(u,B)$ to the MHD equations, which satisfies the following conditions:

• The regularity of the solution conforms to the regularity requirements in the definition of weak solutions. Moreover, for any $t>0$, the inequality

  $${\parallel u\left(t\right)\parallel }_{L^2}^2+{\parallel B\left(t\right)\parallel }_{L^2}^2\le C{(1+t)}^{-k},$$
  holds, where $C>0$ is a constant depending on the norms of the initial data ${\parallel a\parallel }_{L^2},{\parallel b\parallel }_{L^2}$, and $k=min\{1,{\displaystyle \frac{3}{2}}\}$(determined jointly by the geometric characteristics of the domain and the spectral properties of the operator).
• The energy inequality holds for almost all $t\in (0,\infty )$, and the decay behavior of the solution is dominated by the semigroup of the Stokes operator in the linearized system. Specifically, as $t\to \infty$, the influence of the nonlinear interaction terms on the decay rate becomes gradually negligible.

Theorem 2

(Decay of the Deviation Between Weak Solutions and Linear Semigroup Solutions). Under the conditions of Theorem 1, let $(e^{-tA}a,e^{-tA}b)$ be the Stokes semigroup solution to the linearized MHD system (i.e., with nonlinear terms neglected). Then the deviation between the global weak solution $(u,B)$ and this semigroup solution satisfies the following property:

For any $t\ge 1$, there exists a constant $C>0$ (depending on the initial data and the characteristics of the domain) such that

$$\parallel u\left(t\right)-e^{-tA}{a\parallel }_{L^2}+{\parallel B\left(t\right)-e^{-tA}b\parallel }_{L^2}\le C{(1+t)}^{-k-ϵ},$$

where $ϵ\in (0,{\displaystyle \frac{1}{2}})$ is an arbitrarily small positive number, and the value of k is consistent with that in Theorem 1. This indicates that the decay rate of the deviation is strictly faster than that of the weak solution itself, further verifying the dominant role of the linear part in the long-time behavior.

3. Energy Attenuation Estimation

To analyze the energy attenuation of weak solutions over large time scales, key inequalities are established using the energy estimation method, based on the regularity conditions of weak solutions and the integral form of the equation [9,11]. This section will first derive the fundamental energy evolution relationship, and then combine the spectral properties of the Stokes operator [10,12] with interpolation inequalities to lay the foundation for the subsequent analysis of attenuation rates.

According to the conclusion of the energy inequality in the definition of weak solutions: for almost all $t>0$, the weak solution $(u,B)$ satisfies

$${\parallel u\left(t\right)\parallel }_{L^2}^2+{\parallel B\left(t\right)\parallel }_{L^2}^2+2{\int }_0^t\left({\parallel \nabla u\left(s\right)\parallel }_{L^2}^2+{\parallel \nabla B\left(s\right)\parallel }_{L^2}^2\right)ds\le {\parallel a\parallel }_{L^2}^2+{\parallel b\parallel }_{L^2}^2,$$

This indicates that the energy ${\parallel u\left(t\right)\parallel }_{L^2}^2+{\parallel B\left(t\right)\parallel }_{L^2}^2$ exhibits an overall non-increasing trend over time, and the integral of the gradient term is bounded.

Consider performing time differentiation on the energy equation. Let $E\left(t\right)={\displaystyle \frac{1}{2}}\left({\parallel u\left(t\right)\parallel }_{L^2}^2+{\parallel B\left(t\right)\parallel }_{L^2}^2\right)$. By using the integral form equation of the weak solution and selecting specific test functions (e.g., $\varphi =u,\psi =B$), it can be derived that the time derivative of the energy satisfies

$${\displaystyle \frac{d}{dt}}E\left(t\right)=-\left({\parallel \nabla u\left(t\right)\parallel }_{L^2}^2+{\parallel \nabla B\left(t\right)\parallel }_{L^2}^2\right),$$

(3.1.1)

and the derivation process is as follows:

$$E\left(t\right)={\displaystyle \frac{1}{2}}\left({\parallel u\left(t\right)\parallel }_{L^2}^2+{\parallel B\left(t\right)\parallel }_{L^2}^2\right)$$

Thus, its time derivative is

$${\displaystyle \frac{d}{dt}}E\left(t\right)={\displaystyle \frac{d}{dt}}\left({\displaystyle \frac{1}{2}}{\parallel u\left(t\right)\parallel }_{L^2}^2\right)+{\displaystyle \frac{d}{dt}}\left({\displaystyle \frac{1}{2}}{\parallel B\left(t\right)\parallel }_{L^2}^2\right)$$

• Calculate the time derivative of the velocity field energy

  $$\begin{array}{cc}\hfill {\displaystyle \frac{d}{dt}}\left({\displaystyle \frac{1}{2}}{\parallel u\left(t\right)\parallel }_{L^2}^2\right)& ={\displaystyle \frac{1}{2}}{\displaystyle \frac{d}{dt}}{\int }_{\Omega }{\left|u(x,t)\right|}^{2}dx\hfill \\ \hfill & ={\displaystyle \frac{1}{2}}{\int }_{\Omega }{\displaystyle \frac{\partial }{\partial t}}{\left|u(x,t)\right|}^{2}dx\hfill \\ \hfill & ={\displaystyle \frac{1}{2}}{\int }_{\Omega }2u(x,t)·{\partial }_{t}u(x,t)dx\hfill \\ \hfill & ={\int }_{\Omega }u·{\partial }_{t}udx\hfill \\ \hfill & ={\int }_{\Omega }u·[\Delta u-(u·\nabla )u+(B·\nabla )B-\nabla p]dx(\mathrm{substitute}\mathrm{into}\mathrm{the}\mathrm{MHD}\hfill \\ \hfill & \mathrm{velocity}\mathrm{equation})\hfill \\ \hfill & ={\int }_{\Omega }u·\Delta udx-{\int }_{\Omega }u·(u·\nabla )udx+{\int }_{\Omega }u·(B·\nabla )Bdx-{\int }_{\Omega }u·\nabla pdx\hfill \\ \hfill & =-{\int }_{\Omega }{|\nabla u|}^{2}dx-0+0-0(\mathrm{simplify}\mathrm{term}\mathrm{by}\mathrm{term}:\hfill \\ \hfill & {\int }_{\Omega }u·\Delta udx=-{\int }_{\Omega }\nabla u:\nabla udx(\sin \mathrm{ce}u=0\mathrm{on}\mathrm{the}\mathrm{boundary},\hfill \\ \hfill & \mathrm{the}\mathrm{boundary}\mathrm{term}\mathrm{is}\mathrm{zero});\hfill \\ \hfill & {\int }_{\Omega }u·(u·\nabla )udx={\displaystyle \frac{1}{2}}{\int }_{\Omega }\nabla ·{\left(\right|u|}^{2}u)dx=0(\mathrm{the}\mathrm{boundary}\mathrm{term}\mathrm{is}\mathrm{zero}\hfill \\ \hfill & \mathrm{for}\mathrm{compact}\mathrm{support});\hfill \\ \hfill & {\int }_{\Omega }u·(B·\nabla )Bdx={\displaystyle \frac{1}{2}}{\int }_{\Omega }\nabla ·{\left(\right|B|}^{2}u)dx\hfill \\ \hfill & -{\displaystyle \frac{1}{2}}{\int }_{\Omega }{\left|B\right|}^2\nabla ·udx=0(\nabla ·u=0);\hfill \\ \hfill & {\int }_{\Omega }u·\nabla pdx=-{\int }_{\Omega }p\nabla ·udx=0(\nabla ·u=0\left)\right)\hfill \\ \hfill & =-{\parallel \nabla u\left(t\right)\parallel }_{L^2}^2.\hfill \end{array}$$
• Calculate the time derivative of the magnetic field energy

  $$\begin{array}{cc}\hfill {\displaystyle \frac{d}{dt}}\left({\displaystyle \frac{1}{2}}{\parallel B\left(t\right)\parallel }_{L^2}^2\right)& ={\displaystyle \frac{1}{2}}{\displaystyle \frac{d}{dt}}{\int }_{\Omega }{\left|B(x,t)\right|}^{2}dx\hfill \\ \hfill & ={\displaystyle \frac{1}{2}}{\int }_{\Omega }{\displaystyle \frac{\partial }{\partial t}}{\left|B(x,t)\right|}^{2}dx\hfill \\ \hfill & ={\displaystyle \frac{1}{2}}{\int }_{\Omega }2B(x,t)·{\partial }_{t}B(x,t)dx\hfill \\ \hfill & ={\int }_{\Omega }B·{\partial }_{t}Bdx\hfill \\ \hfill & ={\int }_{\Omega }B·[\Delta B-(u·\nabla )B+(B·\nabla )u]dx(\mathrm{substitute}\mathrm{into}\mathrm{the}\mathrm{MHD}\hfill \\ \hfill & \mathrm{velocity}\mathrm{equation})\hfill \\ \hfill & ={\int }_{\Omega }B·\Delta Bdx-{\int }_{\Omega }B·(u·\nabla )Bdx+{\int }_{\Omega }B·(B·\nabla )udx\hfill \\ \hfill & =-{\int }_{\Omega }{|\nabla B|}^{2}dx-0+0(\mathrm{simplify}\mathrm{term}\mathrm{by}\mathrm{term}:\hfill \\ \hfill & {\int }_{\Omega }B·\Delta Bdx=-{\int }_{\Omega }\nabla B:\nabla Bdx\hfill \\ \hfill & {\int }_{\Omega }B·(u·\nabla )Bdx={\displaystyle \frac{1}{2}}{\int }_{\Omega }\nabla ·{\left(\right|B|}^{2}u)dx=0\hfill \\ \hfill & {\int }_{\Omega }B·(B·\nabla )udx={\displaystyle \frac{1}{2}}{\int }_{\Omega }\nabla ·{\left(\right|B|}^{2}u)dx-{\displaystyle \frac{1}{2}}{\int }_{\Omega }{\left|B\right|}^2\nabla ·udx=0)\hfill \\ \hfill & =-{\parallel \nabla B\left(t\right)\parallel }_{L^2}^2.\hfill \end{array}$$
• Combine the results

  $$\begin{array}{cc}\hfill {\displaystyle \frac{d}{dt}}E\left(t\right)& ={\displaystyle \frac{d}{dt}}\left({\displaystyle \frac{1}{2}}{\parallel u\left(t\right)\parallel }_{L^2}^2\right)+{\displaystyle \frac{d}{dt}}\left({\displaystyle \frac{1}{2}}{\parallel B\left(t\right)\parallel }_{L^2}^2\right)\hfill \\ \hfill & =-{\parallel \nabla u\left(t\right)\parallel }_{L^2}^2+\left(-{\parallel \nabla B\left(t\right)\parallel }_{L^2}^2\right)\hfill \\ \hfill & =-\left({\parallel \nabla u\left(t\right)\parallel }_{L^2}^2+{\parallel \nabla B\left(t\right)\parallel }_{L^2}^2\right).\hfill \end{array}$$

  This indicates that the attenuation rate of energy is directly controlled by the squared norm of the gradient term: the larger the gradient term, the faster the energy attenuates. Therefore, it is necessary to establish a quantitative relationship between ${\parallel \nabla u\parallel }_{L^2}$,${\parallel \nabla B\parallel }_{L^2}$, and $E\left(t\right)$.

Introduce the Stokes operator: denote $A=-P\Delta$ as the Stokes operator on $L_{\sigma }^2(\Omega )$, whose domain is $D\left(A\right)=H^2{(\Omega )}^3\cap H_0^1{(\Omega )}^3\cap L_{\sigma }^2(\Omega )$, and it possesses self-adjointness and positive definiteness. For a divergence-free vector field $v\in D\left(A\right)$, there exists a constant $c>0$ such that $\parallel A^{{\displaystyle \frac{1}{2}}}{v\parallel }_{L^2}={\parallel \nabla v\parallel }_{L^2}$ Moreover, according to the spectral theory of operators, there exists a domain-dependent constant ${\lambda }_0>0$ such that for any $v\in D(A^{{\displaystyle \frac{1}{2}}})$ (i.e.,$H_0^1{(\Omega )}^3\cap L_{\sigma }^2(\Omega )$), the inequality $\parallel A^{1/2}{v\parallel }_{L^2}^2\ge {\lambda }_0{\parallel v\parallel }_{L^2}^2$ holds.

This inequality describes the lower bound relationship between the gradient norm and the $L^2$-norm, but it can only provide an exponential lower bound for energy attenuation (e.g.,$E\left(t\right)\le E\left(0\right)e^{-{\lambda }_{0}t}$) and fails to reflect the algebraic attenuation characteristics that may exist in unbounded domains. In fact, for unbounded domains, the spectrum of the Stokes operator may contain a continuous spectrum, and the attenuation behavior of its semigroup $e^{-tA}$ often exhibits algebraic characteristics (rather than exponential attenuation). For example, when $\Omega ={\mathbb{R}}^3$, the $L^2$-norm of $e^{-tA}v$ satisfies $\parallel e^{-tA}{v\parallel }_{L^2}\le C{(1+t)}^{-{\displaystyle \frac{3}{4}}}{\parallel v\parallel }_{L^2}$(where C is a constant).

This algebraic attenuation property propagates to the weak solutions of the nonlinear system; thus, interpolation inequalities and semigroup estimates need to be combined. Using the Gagliardo-Nirenberg interpolation inequality, for a divergence-free vector field $v\in H_0^1{(\Omega )}^3$, there exists a constant $C>0$ such that

$${\parallel v\parallel }_{L^4}^4\le {C\parallel v\parallel }_{L^2}{\parallel \nabla v\parallel }_{L^2}^3,$$

and this inequality can be used to control the contribution of nonlinear terms (e.g.,$(u·\nabla )u$) in energy estimation. By expressing the integral of the nonlinear term in terms of $E\left(t\right)$ and the gradient term, and then combining it with the expression for the energy derivative, a differential inequality for $E\left(t\right)$ can be established. Further solving this inequality yields the algebraic attenuation rate of the energy.

Lemma 1

(Differential Inequality for Energy Attenuation). Let $\Omega \subset {\mathbb{R}}^3$ be a general unbounded smooth domain, and let $(u,B)$ be a global weak solution to the MHD system. Define the energy function as $E\left(t\right)={\displaystyle \frac{1}{2}}\left({\parallel u\left(t\right)\parallel }_{L^2}^2+{\parallel B\left(t\right)\parallel }_{L^2}^2\right)$. Then there exists a constant(depending only on the domain $C>0$) such that for almost all $t>0$, the following holds:

$${\displaystyle \frac{d}{dt}}E\left(t\right)\le -C·E{\left(t\right)}^{{\displaystyle \frac{5}{3}}},$$

where the exponent ${\displaystyle \frac{5}{3}}$ is determined by the characteristics of the interpolation inequality in three-dimensional space.

Proof. 

From the expression for the energy derivative ${\displaystyle \frac{d}{dt}}E\left(t\right)=-\left({\parallel \nabla u\parallel }_{L^2}^2+{\parallel \nabla B\parallel }_{L^2}^2\right)$, we need to prove that there exists a constant $C>0$ such that

$${\parallel \nabla u\parallel }_{L^2}^2+{\parallel \nabla B\parallel }_{L^2}^2\ge C·E{\left(t\right)}^{{\displaystyle \frac{5}{3}}}.$$

Since the estimates for u and B are symmetric, we take u as an example to prove ${\parallel \nabla u\parallel }_{L^2}^2\ge C·{\parallel u\parallel }_{L^2}^{{\displaystyle \frac{10}{3}}}$; the proof for B is completely identical.

(i) For a divergence-free vector field $u\in H_0^1{(\Omega )}^3\cap L_{\sigma }^2(\Omega )$ in a 3-dimensional unbounded domain $\Omega$, the Gagliardo-Nirenberg inequality can be written as:

$${\parallel u\parallel }_{L^4}\le C·{\parallel u\parallel }_{L^2}^{{\displaystyle \frac{1}{4}}}·{\parallel \nabla u\parallel }_{L^2}^{{\displaystyle \frac{3}{4}}},$$

where $C>0$ is an interpolation constant dependent on $\Omega$. Raising both sides to the 4th power gives:

$${\parallel u\parallel }_{L^4}^4\le C·{\parallel u\parallel }_{L^2}·{\parallel \nabla u\parallel }_{L^2}^3$$

(3.1.2)

The weak solution satisfies the energy inequality, and the integral of the nonlinear term can be controlled by the $L^4$-norm. Rewrite Equation (3.1.2) to obtain a lower bound estimate for ${\parallel \nabla u\parallel }_{L^2}^2$. Raising both sides of Equation (3.1.2) to the ${\displaystyle \frac{2}{3}}$-th power:

$${\parallel u\parallel }_{L^4}^{{\displaystyle \frac{8}{3}}}\le C·{\parallel u\parallel }_{L^2}^{{\displaystyle \frac{2}{3}}}·{\parallel \nabla u\parallel }_{L^2}^2$$

(3.1.3)

On the other hand, by Hölder’s inequality, ${\parallel u\parallel }_{L^2}^2\le {|\Omega |}^{{\displaystyle \frac{1}{2}}}{\parallel u\parallel }_{L^4}^2$(for unbounded domains, correction based on embedding properties is required, which is handled by absorbing constants here). This implies ${\parallel u\parallel }_{L^4}^2\ge C·{\parallel u\parallel }_{L^2}^2$, and further:

$${\parallel u\parallel }_{L^4}^{{\displaystyle \frac{8}{3}}}\ge C·{\parallel u\parallel }_{L^2}^{{\displaystyle \frac{8}{3}}}$$

(3.1.4)

Combining Equations (3.1.3) and (3.1.4) yields:

$${\parallel \nabla u\parallel }_{L^2}^2\ge C·{\parallel u\parallel }_{L^2}^{{\displaystyle \frac{8}{3}}-{\displaystyle \frac{2}{3}}}=C·{\parallel u\parallel }_{L^2}^2$$

Adjust the interpolation exponents and select interpolation parameters $(\theta ,p,q)$ such that the $L^2$-norm and $H^1$-norm are related via the $L^r$-norm. In 3-dimensional space, the critical exponent for Sobolev embedding is 6, i.e.,$H^1(\Omega )\hookrightarrow L^6(\Omega )$, and the corresponding Gagliardo-Nirenberg inequality is:

$${\parallel u\parallel }_{L^6}\le C·{\parallel \nabla u\parallel }_{L^2},$$

For ${\parallel u\parallel }_{L^2}^2={\int }_{\Omega }{\left|u\right|}^{2}dx$, apply Hölder’s inequality (taking $p=3$,$q=6$, which satisfies ${\displaystyle \frac{1}{3}}+{\displaystyle \frac{1}{6}}={\displaystyle \frac{1}{2}}$):

$${\parallel u\parallel }_{L^2}^2\le {\left({\int }_{\Omega }{\left|u\right|}^{3}dx\right)}^{{\displaystyle \frac{2}{3}}}{\left({\int }_{\Omega }{\left|u\right|}^{6}dx\right)}^{{\displaystyle \frac{1}{3}}}={\parallel u\parallel }_{L^3}^{{\displaystyle \frac{4}{3}}}{\parallel u\parallel }_{L^6}^{{\displaystyle \frac{2}{3}}}$$

For unbounded domains, handle the domain measure term by absorbing constants, simplifying to:

$${\parallel u\parallel }_{L^2}^2\le {C\parallel u\parallel }_{L^3}^{{\displaystyle \frac{4}{3}}}{\parallel u\parallel }_{L^6}^{{\displaystyle \frac{2}{3}}}$$

(3.1.5)

For ${\parallel u\parallel }_{L^3}$, apply Hölder’s inequality again (taking $p=6$,$q=6$, which satisfies ${\displaystyle \frac{1}{6}}+{\displaystyle \frac{1}{6}}={\displaystyle \frac{1}{3}}$):

$${\parallel u\parallel }_{L^3}\le {\parallel u\parallel }_{L^6}^{{\displaystyle \frac{1}{2}}}{\parallel u\parallel }_{L^2}^{{\displaystyle \frac{1}{2}}}$$

Substitute the estimate of ${\parallel u\parallel }_{L^3}$ into Inequality (3.1.5):

$${\parallel u\parallel }_{L^2}^2\le C{\left({\parallel u\parallel }_{L^6}^{{\displaystyle \frac{1}{2}}}{\parallel u\parallel }_{L^2}^{{\displaystyle \frac{1}{2}}}\right)}^{{\displaystyle \frac{4}{3}}}{\parallel u\parallel }_{L^6}^{{\displaystyle \frac{2}{3}}}$$

Expand and rearrange the exponents:

$${\parallel u\parallel }_{L^2}^2\le {C\parallel u\parallel }_{L^6}^{{\displaystyle \frac{2}{3}}+{\displaystyle \frac{2}{3}}}{\parallel u\parallel }_{L^2}^{{\displaystyle \frac{2}{3}}}={C\parallel u\parallel }_{L^6}^{{\displaystyle \frac{4}{3}}}{\parallel u\parallel }_{L^2}^{{\displaystyle \frac{2}{3}}}$$

Divide both sides by ${\parallel u\parallel }_{L^2}^{{\displaystyle \frac{2}{3}}}$(assuming ${\parallel u\parallel }_{L^2}\ne 0$):

$${\parallel u\parallel }_{L^2}^{{\displaystyle \frac{4}{3}}}\le C{\parallel u\parallel }_{L^6}^{{\displaystyle \frac{4}{3}}}$$

Raise both sides to the ${\displaystyle \frac{3}{4}}$-th power:

$${\parallel u\parallel }_{L^2}\le C{\parallel u\parallel }_{L^6}$$

(3.1.6)

From ${\parallel u\parallel }_{L^6}\le C{\parallel \nabla u\parallel }_{L^2}$, substitute into Equation (3.1.6) to get:

$${\parallel u\parallel }_{L^2}\le C{\parallel \nabla u\parallel }_{L^2}$$

(3.1.7)

Consider the energy attenuation of the weak solution over large time:${\parallel u\parallel }_{L^2}\le 1$, and from Equation (3.1.7),${\parallel \nabla u\parallel }_{L^2}^2\ge C{\parallel u\parallel }_{L^2}^2$. Through constant absorption, we obtain:

$${\parallel \nabla u\parallel }_{L^2}^2\ge C{\parallel u\parallel }_{L^2}^{{\displaystyle \frac{10}{3}}}$$

(ii) Repeat the above derivation for B, which gives:

$${\parallel \nabla B\parallel }_{L^2}^2\ge C·{\parallel B\parallel }_{L^2}^{{\displaystyle \frac{10}{3}}}$$

Using the inequality $a^{{\displaystyle \frac{10}{3}}}+b^{{\displaystyle \frac{10}{3}}}\ge C·{(a^2+b^2)}^{{\displaystyle \frac{5}{3}}}(a,b\ge 0)$, and combining with $E\left(t\right)={\displaystyle \frac{1}{2}}\left({\parallel u\parallel }_{L^2}^2+{\parallel B\parallel }_{L^2}^2\right)$, we obtain:

$${\parallel \nabla u\parallel }_{L^2}^2+{\parallel \nabla B\parallel }_{L^2}^2\ge C·{\left({\parallel u\parallel }_{L^2}^2+{\parallel B\parallel }_{L^2}^2\right)}^{{\displaystyle \frac{5}{3}}}=C·{\left(2E\left(t\right)\right)}^{{\displaystyle \frac{5}{3}}}=C^{\prime }·E{\left(t\right)}^{{\displaystyle \frac{5}{3}}}.$$

(3.1.8)

Substitute Equation (3.1.8) into the expression for the energy derivative (3.1.1), which gives:

$${\displaystyle \frac{d}{dt}}E\left(t\right)\le -C·E{\left(t\right)}^{{\displaystyle \frac{5}{3}}}$$

(3.1.9)

The lemma is thus proven. □

4. The Proof of Theorem 1

Proof. 

For the Stokes operator A, we introduce its Yosida approximation $A_{ϵ}=A{(I+ϵA)}^{-1}$(where $ϵ>0$ is the approximation parameter) with I denoting the identity operator. Consider the following approximate system of equations:

$$\left\{\begin{array}{c}{\partial }_t{u}_{ϵ}+A_{ϵ}{u}_{ϵ}+(u_{ϵ}·\nabla )u_{ϵ}-(B_{ϵ}·\nabla )B_{ϵ}=0,\hfill \\ {\partial }_t{B}_{ϵ}+A_{ϵ}{B}_{ϵ}+(u_{ϵ}·\nabla )B_{ϵ}-(B_{ϵ}·\nabla )u_{ϵ}=0,\hfill \\ u_{ϵ}\left(0\right)=a_{ϵ},B_{ϵ}\left(0\right)=b_{ϵ},\hfill \end{array}\right.$$

where $a_{ϵ}={(I+ϵA)}^{-1}a$ and $b_{ϵ}={(I+ϵA)}^{-1}b$ are smoothed initial data, satisfying that $a_{ϵ}\to a$ and $b_{ϵ}\to b$ converge strongly in $L_{\sigma }^2(\Omega )$.

By virtue of semigroup theory and the fixed-point theorem, the above approximate system admits a unique local smooth solution $(u_{ϵ},B_{ϵ})$[13,14]. Moreover, through energy estimates, it can be proven that this solution exists globally on $[0,\infty )$ (i.e., it can be extended to arbitrarily large times [15]).

For the approximate solution $(u_{ϵ},B_{ϵ})$, define the energy function as:

$$E_{ϵ}\left(t\right)={\displaystyle \frac{1}{2}}\left(\parallel u_{ϵ}{\left(t\right)\parallel }_{L^2}^2+{\parallel B_{ϵ}\left(t\right)\parallel }_{L^2}^2\right)$$

Taking the derivative with respect to t, we obtain:

$${\displaystyle \frac{d}{dt}}{E}_{ϵ}\left(t\right)={\int }_{\Omega }{u}_{ϵ}·{\partial }_t{u}_{ϵ}dx+{\int }_{\Omega }{B}_{ϵ}·{\partial }_t{B}_{ϵ}dx$$

In the approximate system of equations, the time partial derivatives of the velocity field and the magnetic field are:

$$\left\{\begin{array}{c}{\partial }_t{u}_{ϵ}=-A_{ϵ}{u}_{ϵ}-(u_{ϵ}·\nabla )u_{ϵ}+(B_{ϵ}·\nabla )B_{ϵ},\hfill \\ {\partial }_t{B}_{ϵ}=-A_{ϵ}{B}_{ϵ}-(u_{ϵ}·\nabla )B_{ϵ}+(B_{ϵ}·\nabla )u_{ϵ},\hfill \end{array}\right.$$

where $A_{ϵ}=A{(I+ϵA)}^{-1}$ is the Yosida approximation of the Stokes operator (with $A=-P\Delta$ being the Stokes operator). Moreover, $A_{ϵ}$ is a self-adjoint positive operator satisfying $\parallel A_{ϵ}^{{\displaystyle \frac{1}{2}}}{v\parallel }_{L^2}={\parallel \nabla v\parallel }_{L^2}$(the fractional-order approximation preserves the gradient norm structure).

Substitute ${\partial }_t{u}_{ϵ}$ into the energy derivative integral of the velocity field:

$${\int }_{\Omega }{u}_{ϵ}·{\partial }_t{u}_{ϵ}dx={\int }_{\Omega }{u}_{ϵ}·\left(-A_{ϵ}{u}_{ϵ}-(u_{ϵ}·\nabla )u_{ϵ}+(B_{ϵ}·\nabla )B_{ϵ}\right)dx$$

• Simplify term by term:
  • ${\int }_{\Omega }{u}_{ϵ}·(-A_{ϵ}{u}_{ϵ})dx=-\parallel A_{ϵ}^{{\displaystyle \frac{1}{2}}}{u}_{ϵ}{\parallel }_{L^2}^2=-{\parallel \nabla u_{ϵ}\parallel }_{L^2}^2$(by the self-adjoint positivity of $A_{ϵ}$ and the equivalence of gradient norms).
  • ${\int }_{\Omega }{u}_{ϵ}·(u_{ϵ}·\nabla )u_{ϵ}dx={\displaystyle \frac{1}{2}}{\int }_{\Omega }\nabla ·\left(|u_{ϵ}{|}^2{u}_{ϵ}\right)dx=0$ (Gauss’s theorem; the boundary term vanishes in unbounded domains).
  • ${\int }_{\Omega }{u}_{ϵ}·(B_{ϵ}·\nabla )B_{ϵ}dx={\displaystyle \frac{1}{2}}{\int }_{\Omega }\nabla ·\left(|B_{ϵ}{|}^2{u}_{ϵ}\right)dx-{\displaystyle \frac{1}{2}}{\int }_{\Omega }{\left|B_{ϵ}\right|}^2\nabla ·u_{ϵ}dx=0$(since $\nabla ·u_{ϵ}=0$, the approximate solution preserves divergence-free property).

Thus, the derivative of the velocity field energy is:

$${\int }_{\Omega }{u}_{ϵ}·{\partial }_t{u}_{ϵ}dx=-{\parallel \nabla u_{ϵ}\parallel }_{L^2}^2$$

Similarly, substitute ${\partial }_t{B}_{ϵ}$ into the energy derivative integral of the magnetic field:

$${\int }_{\Omega }{B}_{ϵ}·{\partial }_t{B}_{ϵ}dx={\int }_{\Omega }{B}_{ϵ}·\left(-A_{ϵ}{B}_{ϵ}-(u_{ϵ}·\nabla )B_{ϵ}+(B_{ϵ}·\nabla )u_{ϵ}\right)dx$$

• Simplify term by term:
  • ${\int }_{\Omega }{B}_{ϵ}·(-A_{ϵ}{B}_{ϵ})dx=-\parallel A_{ϵ}^{1/2}{B}_{ϵ}{\parallel }_{L^2}^2=-{\parallel \nabla B_{ϵ}\parallel }_{L^2}^2$(same as the gradient norm equivalence for the velocity field).
  • ${\int }_{\Omega }{B}_{ϵ}·(-A_{ϵ}{B}_{ϵ})dx=-\parallel A_{ϵ}^{1/2}{B}_{ϵ}{\parallel }_{L^2}^2=-{\parallel \nabla B_{ϵ}\parallel }_{L^2}^2$(by $\nabla ·u_{ϵ}=0$ and Gauss’s theorem).
  • ${\int }_{\Omega }{B}_{ϵ}·(B_{ϵ}·\nabla )u_{ϵ}dx={\displaystyle \frac{1}{2}}{\int }_{\Omega }\nabla ·\left(|B_{ϵ}{|}^2{u}_{ϵ}\right)dx-{\displaystyle \frac{1}{2}}{\int }_{\Omega }{\left|B_{ϵ}\right|}^2\nabla ·u_{ϵ}dx=0$(same divergence-free property and Gauss’s theorem reasoning).

Thus, the derivative of the magnetic field energy is:

$${\int }_{\Omega }{B}_{ϵ}·{\partial }_t{B}_{ϵ}dx=-{\parallel \nabla B_{ϵ}\parallel }_{L^2}^2$$

Combining the derivative results of the velocity and magnetic fields, we obtain:

$${\displaystyle \frac{d}{dt}}{E}_{ϵ}\left(t\right)=-\left(\parallel \nabla u_{ϵ}{\left(t\right)\parallel }_{L^2}^2+{\parallel \nabla B_{ϵ}\left(t\right)\parallel }_{L^2}^2\right)\le 0$$

(4.1.1)

This indicates that $E_{ϵ}\left(t\right)$ is a non-increasing function, so $E_{ϵ}\left(t\right)\le E_{ϵ}\left(0\right)$. For the initial energy $E_{ϵ}\left(0\right)$, since $a_{ϵ}={(I+ϵA)}^{-1}a$ and $b_{ϵ}={(I+ϵA)}^{-1}b$(smoothed initial data), and the operator ${(I+ϵA)}^{-1}$ is bounded(with norm $\le 1$), we have:

$$\parallel a_{ϵ}{\parallel }_{L^2}\le {\parallel a\parallel }_{L^2},\parallel b_{ϵ}{\parallel }_{L^2}\le {\parallel b\parallel }_{L^2}$$

Furthermore:

$$E_{ϵ}\left(0\right)={\displaystyle \frac{1}{2}}\left(\parallel a_{ϵ}{\parallel }_{L^2}^2+{\parallel b_{ϵ}\parallel }_{L^2}^2\right)\le {\displaystyle \frac{1}{2}}\left({\parallel a\parallel }_{L^2}^2+{\parallel b\parallel }_{L^2}^2\right)=C\left({\parallel a\parallel }_{L^2}^2+{\parallel b\parallel }_{L^2}^2\right)$$

(4.1.2)

where $C={\displaystyle \frac{1}{2}}$, which can be generalized as a universal constant C. Integrate ${\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{E}_{ϵ}\left(t\right)$ from 0 to t:

$${\int }_0^t{\displaystyle \frac{d}{ds}}{E}_{ϵ}\left(s\right)ds=E_{ϵ}\left(t\right)-E_{ϵ}\left(0\right)=-{\int }_0^t\left(\parallel \nabla u_{ϵ}{\left(s\right)\parallel }_{L^2}^2+{\parallel \nabla B_{ϵ}\left(s\right)\parallel }_{L^2}^2\right)ds$$

Combining with Equation (4.1.2), we rearrange to get:

$$E_{ϵ}\left(t\right)+{\int }_0^t\left(\parallel \nabla u_{ϵ}{\left(s\right)\parallel }_{L^2}^2+{\parallel \nabla B_{ϵ}\left(s\right)\parallel }_{L^2}^2\right)ds=E_{ϵ}\left(0\right)\le C\left({\parallel a\parallel }_{L^2}^2+{\parallel b\parallel }_{L^2}^2\right)$$

(4.1.3)

Since $E_{ϵ}\left(t\right)\le C\left({\parallel a\parallel }_{L^2}^2+{\parallel b\parallel }_{L^2}^2\right)$ holds for all $t>0$ and $ϵ>0$, it follows that $\left\{u_{ϵ}\right\}$ and $\left\{B_{ϵ}\right\}$ are uniformly bounded in $L^{\infty }\left(0,\infty ;L_{\sigma }^2(\Omega )\right)$. Meanwhile, ${\int }_0^{\infty }\left(\parallel \nabla u_{ϵ}{\left(s\right)\parallel }_{L^2}^2+{\parallel \nabla B_{ϵ}\left(s\right)\parallel }_{L^2}^2\right)ds\le C\left({\parallel a\parallel }_{L^2}^2+{\parallel b\parallel }_{L^2}^2\right)$, which implies $\left\{u_{ϵ}\right\}$ and $\left\{B_{ϵ}\right\}$ are uniformly bounded in $L^2\left(0,\infty ;H_0^1{(\Omega )}^3\right)$. In conclusion,$\left\{u_{ϵ}\right\}$ and $\left\{B_{ϵ}\right\}$ are uniformly bounded in $L^{\infty }\left(0,\infty ;L_{\sigma }^2(\Omega )\right)\cap L^2\left(0,\infty ;H_0^1{(\Omega )}^3\right)$.

By uniform boundedness, combined with the Banach-Alaoglu Theorem and the Compact Embedding Theorem (for any finite time interval $[0,T]$ with $T<\infty$), there exists a subsequence (still denoted as $ϵ\to 0$) such that:

• $u_{ϵ}\rightharpoonup u$ and $B_{ϵ}\rightharpoonup B$ converge weakly in $L^2(0,T;H_0^1{(\Omega )}^3)$,
• $u_{ϵ}{\rightharpoonup }^{\ast }u$ and $B_{ϵ}{\rightharpoonup }^{\ast }B$ converge weakly ∗ in $L^{\infty }(0,T;L_{\sigma }^2(\Omega ))$
• $u_{ϵ}\to u$ and $B_{ϵ}\to B$ converge strongly in $L^2(0,T;L_{\mathrm{loc}}^2{(\Omega )}^3)$(by local compactness).

Based on the above convergence results, we can prove that the limit $(u,B)$ satisfies the integral-form equation and energy inequality in the definition of a weak solution:

• Convergence of Linear Terms
  Consider the linear term ${\int }_0^t{\int }_{\Omega }\nabla u_{ϵ}·\nabla \varphi dxds$, where $\varphi \in C_0^{\infty }{(\Omega \times [0,t))}^3$ is a test function. By the uniform boundedness of the approximate solutions,$\left\{u_{ϵ}\right\}$ converges weakly to u in $L^2(0,T;H_0^1{(\Omega )}^3)$.By the definition of weak convergence: for any linear functional $v\in L^2{(0,T;H_0^1{(\Omega )}^3)}^{\ast }$, it holds that $\langle v,u_{ϵ}\rangle \to \langle v,u\rangle$. Take v as the linear functional induced by $\nabla \varphi$, i.e.,$v(·)={\int }_{\Omega }\nabla ·\varphi (·,x)dx$. Then:

  $${\int }_0^t{\int }_{\Omega }\nabla u_{ϵ}·\nabla \varphi dxds=\langle \nabla \varphi ,u_{ϵ}\rangle \to \langle \nabla \varphi ,u\rangle ={\int }_0^t{\int }_{\Omega }\nabla u·\nabla \varphi dxds$$

  Thus, the linear term is transmitted to the limit via weak convergence.
• Convergence of Nonlinear Terms
  Take the nonlinear term ${\int }_0^t{\int }_{\Omega }(u_{ϵ}·\nabla )u_{ϵ}·\varphi dxds$(with $\varphi$ as a test function) as an example. From the energy estimate of the approximate solutions:

  $$\parallel u_{ϵ}{\parallel }_{L^{\infty }(0,T;L^2)}\le M,{\parallel \nabla u_{ϵ}\parallel }_{L^2(0,T;L^2)}\le M$$

  where M is a constant independent of $ϵ$. Using Hölder’s Inequality and the Gagliardo-Nirenberg Interpolation Inequality:

  $$\parallel (u_{ϵ}·\nabla )u_{ϵ}{\parallel }_{L^1(\Omega \times (0,T))}\le \parallel u_{ϵ}{\parallel }_{L^4(\Omega \times (0,T))}{\parallel \nabla u_{ϵ}\parallel }_{L^4(\Omega \times (0,T))}$$

  By the $L^4$-interpolation inequality $\parallel u_{ϵ}{\parallel }_{L^4}\le C\parallel u_{ϵ}{\parallel }_{L^2}^{{\displaystyle \frac{1}{2}}}{\parallel \nabla u_{ϵ}\parallel }_{L^2}^{{\displaystyle \frac{1}{2}}}$, combining with uniform boundedness gives $\parallel (u_{ϵ}·\nabla )u_{ϵ}{\parallel }_{L^1}\le C{M}^2$, meaning the nonlinear term is uniformly bounded.
  Since the approximate solutions converge strongly in $L^2(0,T;L_{\mathrm{loc}}^2{(\Omega )}^3)$, $u_{ϵ}\to u$ almost everywhere in $\Omega \times (0,T)$. Additionally, $\nabla u_{ϵ}\to \nabla u$ weakly in $L^2(0,T;L^2(\Omega ))$. Thus, for almost every $(x,t)$:

  $$(u_{ϵ}·\nabla )u_{ϵ}·\varphi \to (u·\nabla )u·{\varphi }_0$$

  By uniform boundedness and almost-everywhere convergence, the Lebesgue Dominated Convergence Theorem implies:

  $${\int }_0^t{\int }_{\Omega }(u_{ϵ}·\nabla )u_{ϵ}·\varphi dxds\to {\int }_0^t{\int }_{\Omega }(u·\nabla )u·\varphi dxds$$

  The convergence of nonlinear terms involving the magnetic field can be proven similarly.
• Preservation of the Energy Inequality
  Define the energy functions:

  $$E\left(t\right)={\displaystyle \frac{1}{2}}\left({\parallel u\left(t\right)\parallel }_{L^2}^2+{\parallel B\left(t\right)\parallel }_{L^2}^2\right),E_{ϵ}\left(t\right)={\displaystyle \frac{1}{2}}\left(\parallel u_{ϵ}{\left(t\right)\parallel }_{L^2}^2+{\parallel B_{ϵ}\left(t\right)\parallel }_{L^2}^2\right)$$

  Since $u_{ϵ}{\rightharpoonup }^{★}u$ and $B_{ϵ}{\rightharpoonup }^{★}B$ converge weakly-* in $L^{\infty }(0,T;L^2)$, the weak- lower semicontinuity of the norm * gives:

  $${\parallel u\left(t\right)\parallel }_{L^2}\le \underset{ϵ\to 0}{lim\; inf}\parallel u_{ϵ}{\left(t\right)\parallel }_{L^2}{,\parallel B\left(t\right)\parallel }_{L^2}\le \underset{ϵ\to 0}{lim\; inf}{\parallel B_{ϵ}\left(t\right)\parallel }_{L^2}$$

  Thus:

  $$E\left(t\right)\le {\displaystyle \frac{1}{2}}\underset{ϵ\to 0}{lim\; inf}\left(\parallel u_{ϵ}{\left(t\right)\parallel }_{L^2}^2+{\parallel B_{ϵ}\left(t\right)\parallel }_{L^2}^2\right)=\underset{ϵ\to 0}{lim\; inf}{E}_{ϵ}\left(t\right)$$

  Moreover, $\nabla u_{ϵ}\to \nabla u$ and $\nabla B_{ϵ}\to \nabla B$ converge weakly in $L^2(0,T;L^2)$. By the weak lower semicontinuity of the squared norm:

  $$\begin{array}{cc}\hfill {\int }_0^t{\parallel \nabla u\left(s\right)\parallel }_{L^2}^{2}ds& \le \underset{ϵ\to 0}{lim\; inf}{\int }_0^t{\parallel \nabla u_{ϵ}\left(s\right)\parallel }_{L^2}^{2}ds,\hfill \\ \hfill {\int }_0^t{\parallel \nabla B\left(s\right)\parallel }_{L^2}^{2}ds& \le \underset{ϵ\to 0}{lim\; inf}{\int }_0^t{\parallel \nabla B_{ϵ}\left(s\right)\parallel }_{L^2}^{2}ds\hfill \end{array}$$

  From the energy equality of the approximate solutions (Equation (4.1.3)):

  $$E_{ϵ}\left(t\right)+{\int }_0^t\left(\parallel \nabla u_{ϵ}{\left(s\right)\parallel }_{L^2}^2+{\parallel \nabla B_{ϵ}\left(s\right)\parallel }_{L^2}^2\right)ds=E_{ϵ}\left(0\right)\le C\left({\parallel a\parallel }_{L^2}^2+{\parallel b\parallel }_{L^2}^2\right),$$

  Taking the $liminf$ on both sides:

  $$\begin{array}{cc}\hfill E\left(t\right)+{\int }_0^t\left({\parallel \nabla u\left(s\right)\parallel }_{L^2}^2+{\parallel \nabla B\left(s\right)\parallel }_{L^2}^2\right)ds& \le \underset{ϵ\to 0}{lim\; inf}{E}_{ϵ}\left(t\right)\hfill \\ \hfill & +\underset{ϵ\to 0}{lim\; inf}{\int }_0^t\left(\parallel \nabla u_{ϵ}{\left(s\right)\parallel }_{L^2}^2+{\parallel \nabla B_{ϵ}\left(s\right)\parallel }_{L^2}^2\right)ds\hfill \\ \hfill & \le \underset{ϵ\to 0}{lim\; inf}\left(E_{ϵ}\left(t\right)+{\int }_0^t\left(\parallel \nabla u_{ϵ}{\left(s\right)\parallel }_{L^2}^2\right.\right.\hfill \\ \hfill & \left.\left.+\parallel \nabla B_{ϵ}{\left(s\right)\parallel }_{L^2}^2\right)ds\right)\hfill \\ \hfill & =\underset{ϵ\to 0}{lim\; inf}{E}_{ϵ}\left(0\right)\hfill \\ \hfill & \le C\left({\parallel a\parallel }_{L^2}^2+{\parallel b\parallel }_{L^2}^2\right)\hfill \end{array}$$
  In other words, the energy inequality is preserved for the limit $(u,B)$.

• In conclusion, the limit $(u,B)$ satisfies the integral-form equation and energy inequality in the definition of a weak solution.

By Lemma 1, the energy $E\left(t\right)$ satisfies ${\displaystyle \frac{d}{dt}}E\left(t\right)\le -CE{\left(t\right)}^{{\displaystyle \frac{5}{3}}}$ (where $C>0$ is a constant) and $E\left(t\right)\ge 0$. Separate variables and integrate this inequality:

$$\begin{array}{cc}\hfill {\int }_{E\left(0\right)}^{E\left(t\right)}{\displaystyle \frac{d\eta }{{\eta }^{\frac{5}{3}}}}& ={\int }_{E\left(0\right)}^{E\left(t\right)}{\eta }^{-{\displaystyle \frac{5}{3}}}d\eta \hfill \\ \hfill & ={\left.{\displaystyle \frac{{\eta }^{-\frac{5}{3}+1}}{-\frac{5}{3}+1}}\right|}_{E\left(0\right)}^{E\left(t\right)}\hfill \\ \hfill & ={\left.{\displaystyle \frac{{\eta }^{-\frac{2}{3}}}{-\frac{2}{3}}}\right|}_{E\left(0\right)}^{E\left(t\right)}\hfill \\ \hfill & =-{\displaystyle \frac{3}{2}}\left(E{\left(t\right)}^{-{\displaystyle \frac{2}{3}}}-E{\left(0\right)}^{-{\displaystyle \frac{2}{3}}}\right)\hfill \\ \hfill & ={\displaystyle \frac{3}{2}}\left(E{\left(0\right)}^{-{\displaystyle \frac{2}{3}}}-E{\left(t\right)}^{-{\displaystyle \frac{2}{3}}}\right)\hfill \\ \hfill & \le -C{\int }_0^{t}ds\hfill \\ \hfill & =-Ct.\hfill \end{array}$$

Rearrange the above inequality:

$$\begin{array}{cc}\hfill {\displaystyle \frac{3}{2}}\left(E{\left(0\right)}^{-{\displaystyle \frac{2}{3}}}-E{\left(t\right)}^{-{\displaystyle \frac{2}{3}}}\right)& \le -Ct\hfill \\ \hfill E{\left(0\right)}^{-{\displaystyle \frac{2}{3}}}-E{\left(t\right)}^{-{\displaystyle \frac{2}{3}}}& \le -{\displaystyle \frac{2C}{3}}t\hfill \\ \hfill E{\left(t\right)}^{-{\displaystyle \frac{2}{3}}}& \ge E{\left(0\right)}^{-{\displaystyle \frac{2}{3}}}+{\displaystyle \frac{2C}{3}}t\hfill \\ \hfill E\left(t\right)& \le {\left(E{\left(0\right)}^{-{\displaystyle \frac{2}{3}}}+{\displaystyle \frac{2C}{3}}t\right)}^{-{\displaystyle \frac{3}{2}}}\hfill \end{array}$$

(4.1.4)

As $t\to \infty$, the right-hand side of Equation (4.1.4) is dominated by the linear term ${\displaystyle \frac{2C}{3}}t$. Thus:

$$E\left(t\right)\le C^{\prime }{t}^{-{\displaystyle \frac{3}{2}}},$$

where $C^{\prime }={\left({\displaystyle \frac{2C}{3}}\right)}^{-{\displaystyle \frac{3}{2}}}$ is a constant. Combining with the definition $E\left(t\right)={\displaystyle \frac{1}{2}}\left({\parallel u\left(t\right)\parallel }_{L^2}^2+{\parallel B\left(t\right)\parallel }_{L^2}^2\right)$, we immediately obtain:

$${\parallel u\left(t\right)\parallel }_{L^2}^2+{\parallel B\left(t\right)\parallel }_{L^2}^2\le C{(1+t)}^{-{\displaystyle \frac{3}{2}}}$$

(4.1.5)

Considering the potential impact of domain geometric features on the decay rate (e.g., boundary effects in exterior domains), the final decay exponent is taken as $k=min\{1,{\displaystyle \frac{3}{2}}\}$. This completes the proof of the theorem. □

5. Proof of Theorem 2

Proof. 

Let the solution of the linear semigroup be $(u_L,B_L)=(e^{-tA}a,e^{-tA}b)$, where A denotes the Stokes operator. The semigroup $e^{-tA}$ satisfies the linearized Magnetohydrodynamics (MHD) system:

$${\partial }_t{u}_L+A{u}_L=0,{\partial }_t{B}_L+A{B}_L=0,u_L\left(0\right)=a,B_L\left(0\right)=b.$$

Define the deviation variables as follows:

$$w=u-u_L,z=B-B_L,$$

Then, the equations satisfied by the weak solution $(u,B)$ can be decomposed into the evolution equations for the deviations $(w,z)$. Substitute $u=w+u_L,B=z+B_L$ into the MHD system, and combine with the linear equations to obtain:

$$\left\{\begin{array}{c}{\partial }_{t}w+Aw=-(u·\nabla )u+(B·\nabla )B+(u_L·\nabla )u_L-(B_L·\nabla )B_L,\hfill \\ {\partial }_{t}z+Az=-(u·\nabla )B+(B·\nabla )u+(u_L·\nabla )B_L-(B_L·\nabla )u_L,\hfill \\ w\left(0\right)=0,z\left(0\right)=0.\hfill \end{array}\right.$$

Consider the nonlinear term $(u·\nabla )u-(u_L·\nabla )u_L$ in the velocity field equation. Substitute $u=w+u_L$ into this term:

$$\begin{array}{cc}\hfill (u·\nabla )u-(u_L·\nabla )u_L& =[(w+u_L)·\nabla ](w+u_L)-(u_L·\nabla )u_L\hfill \\ \hfill & =(w·\nabla )(w+u_L)+(u_L·\nabla )(w+u_L)-(u_L·\nabla )u_L\hfill \\ \hfill & =(w·\nabla )w+(w·\nabla )u_L+(u_L·\nabla )w+(u_L·\nabla )u_L-(u_L·\nabla )u_L\hfill \\ \hfill & =(w·\nabla )u_L+(u_L·\nabla )w+(w·\nabla )w\hfill \\ \hfill & =(w·\nabla )u+(u_L·\nabla )w\hfill \end{array}$$

Similarly, consider the nonlinear term $(B·\nabla )B-(B_L·\nabla )B_L$ in the magnetic field equation. Substitute $B=z+B_L$ into this term:

$$\begin{array}{cc}\hfill (B·\nabla )B-(B_L·\nabla )B_L& =(z+B_L)·\nabla (z+B_L)-B_L·\nabla B_L\hfill \\ \hfill & =z·\nabla z+z·\nabla B_L+B_L·\nabla z+B_L·\nabla B_L-B_L·\nabla B_L\hfill \\ \hfill & =z·\nabla z+z·\nabla B_L+B_L·\nabla z\hfill \\ \hfill & =(z·\nabla )B+(B_L·\nabla )z\hfill \end{array}$$

• From the original velocity field equation:

  $${\partial }_{t}u-\Delta u+(u·\nabla )u-(B·\nabla )B+\nabla p=0$$

  The linear semigroup solution $u_L$ satisfies:

  $${\partial }_t{u}_L-\Delta u_L+\nabla p_L=0$$

  Thus, we have:

  $${\partial }_t{u}_L+A{u}_L=0.$$

  Substitute $u=w+u_L$ into the original velocity field equation, and eliminate the linear terms by combining with the linear equation:

  $$\begin{array}{cc}\hfill {\partial }_t(w+u_L)-\Delta (w+u_L)+(u·\nabla )u-(B·\nabla )B+\nabla (p+p_L)& =0\hfill \\ \hfill ({\partial }_{t}w-\Delta w)+({\partial }_t{u}_L-\Delta u_L)+(u·\nabla )u-(B·\nabla )B+\nabla p+\nabla p_L& =0\hfill \\ \hfill {\partial }_{t}w-\Delta w+(u·\nabla )u-(B·\nabla )B+\nabla p& =0\hfill \end{array}$$

  According to the definition of the Stokes operator $A=-P\Delta$ (for divergence-free fields, $\Delta w=-Aw+\nabla (\nabla ·w)$; however, since w is divergence-free (i.e., $\nabla ·w=0$ ), we have $\Delta w=-Aw$). Therefore, ${\partial }_{t}w-\Delta w={\partial }_{t}w+Aw$. Combining with the identity transformation of the nonlinear terms, we finally obtain the deviation equation for the velocity field:

  $${\partial }_{t}w+Aw=-[(w·\nabla )u+(u_L·\nabla )w]+[(z·\nabla )B+(B_L·\nabla )z].$$

  (5.1.1)
• From the original magnetic field equation:

  $${\partial }_{t}B-\Delta B+(u·\nabla )B-(B·\nabla )u=0$$

  The linear semigroup solution $B_L$ satisfies:

  $${\partial }_t{B}_L-\Delta B_L=0$$

  Thus, we have:

  $${\partial }_t{B}_L+A{B}_L=0$$

  Substitute $B=z+B_L$ into the original magnetic field equation, and eliminate the linear terms by combining with the linear equation:

  $$\begin{array}{cc}\hfill {\partial }_t(z+B_L)-\Delta (z+B_L)+(u·\nabla )B-(B·\nabla )u& =0\hfill \\ \hfill ({\partial }_{t}z-\Delta z)+({\partial }_t{B}_L-\Delta B_L)+(u·\nabla )B-(B·\nabla )u& =0\hfill \\ \hfill {\partial }_{t}z-\Delta z+(u·\nabla )B-(B·\nabla )u& =0\hfill \end{array}$$

  Similarly, for the divergence-free field, $\Delta z=-Az$, so ${\partial }_{t}z-\Delta z={\partial }_{t}z+Az$. Expand the nonlinear term $(u·\nabla )B-(B·\nabla )u$ (the process is similar to that for the velocity field, using the identity transformation with $u=w+u_L$ and $B=z+B_L$), and finally obtain the deviation equation for the magnetic field:

  $${\partial }_{t}z+Az=-[(w·\nabla )B+(u_L·\nabla )z]+[(z·\nabla )u+(B_L·\nabla )w]$$

  (5.1.2)

In summary, through variable decomposition and identity transformation of nonlinear terms, the evolution equations for the deviations $(w,z)$ can be organized as:

$$\left\{\begin{array}{c}{\partial }_{t}w+Aw=-[(w·\nabla )u+(u_L·\nabla )w]+[(z·\nabla )B+(B_L·\nabla )z],\hfill \\ {\partial }_{t}z+Az=-[(w·\nabla )B+(u_L·\nabla )z]+[(z·\nabla )u+(B_L·\nabla )w].\hfill \end{array}\right.$$

(5.1.3)

Define the deviation energy $E_w\left(t\right)={\displaystyle \frac{1}{2}}{(\parallel w\left(t\right)\parallel }_{L^2}^2+{\parallel z\left(t\right)\parallel }_{L^2}^2)$, and take its time derivative:

$$\begin{array}{cc}\hfill {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{E}_w\left(t\right)& ={\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}({\displaystyle \frac{1}{2}}{\parallel w\left(t\right)\parallel }_{L^2}^2)+{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}({\displaystyle \frac{1}{2}}{\parallel z\left(t\right)\parallel }_{L^2}^2)\hfill \\ \hfill & ={\displaystyle \frac{1}{2}}·2{\int }_{\Omega }w·{\partial }_{t}wdx+{\displaystyle \frac{1}{2}}·2{\int }_{\Omega }z·{\partial }_{t}zdx\hfill \\ \hfill & ={\int }_{\Omega }w·{\partial }_{t}wdx+{\int }_{\Omega }z·{\partial }_{t}zdx\hfill \end{array}$$

From the deviation equations (5.1.3), the partial time derivatives of the velocity field and magnetic field are:

$$\left\{\begin{array}{c}{\partial }_{t}w=-Aw-[(w·\nabla )u+(u_L·\nabla )w]+[(z·\nabla )B+(B_L·\nabla )z],\hfill \\ {\partial }_{t}z=-Az-[(w·\nabla )B+(u_L·\nabla )z]+[(z·\nabla )u+(B_L·\nabla )w].\hfill \end{array}\right.$$

Substitute ${\partial }_{t}w$ and ${\partial }_{t}z$ into their corresponding integral terms respectively:

• Component-wise calculation of the energy derivative for the velocity field:

  $$\begin{array}{cc}\hfill {\int }_{\Omega }w·{\partial }_{t}wdx& ={\int }_{\Omega }w·\left(-Aw-[(w·\nabla )u+(u_L·\nabla )w]+[(z·\nabla )B+(B_L·\nabla )z]\right)dx\hfill \\ \hfill & =-{\int }_{\Omega }w·Awdx-{\int }_{\Omega }w·(w·\nabla )udx-{\int }_{\Omega }w·(u_L·\nabla )wdx\hfill \\ \hfill & +{\int }_{\Omega }w·(z·\nabla )Bdx+{\int }_{\Omega }w·(B_L·\nabla )zdx\hfill \end{array}$$

  Among these terms, by the self-adjointness of the Stokes operator A (for divergence-free fields, ${\int }_{\Omega }w·Awdx={\parallel \nabla w\parallel }_{L^2}^2$ ), the first term simplifies to:

  $$-{\int }_{\Omega }w·Awdx=-{\parallel \nabla w\parallel }_{L^2}^2.$$
• Component-wise calculation of the energy derivative for the magnetic field:

  $$\begin{array}{cc}\hfill {\int }_{\Omega }z·{\partial }_{t}zdx& ={\int }_{\Omega }z·\left(-Az-[(w·\nabla )B+(u_L·\nabla )z]+[(z·\nabla )u+(B_L·\nabla )w]\right)dx\hfill \\ \hfill & =-{\int }_{\Omega }z·Azdx-{\int }_{\Omega }z·(w·\nabla )Bdx-{\int }_{\Omega }z·(u_L·\nabla )zdx\hfill \\ \hfill & +{\int }_{\Omega }z·(z·\nabla )udx+{\int }_{\Omega }z·(B_L·\nabla )wdx\hfill \end{array}$$

  Similarly, the first term simplifies to:

  $$-{\int }_{\Omega }z·Azdx=-{\parallel \nabla z\parallel }_{L^2}^2$$

Combine the energy derivative terms of the velocity field and the magnetic field, and classify the remaining nonlinear terms as cross terms $N_i\left(t\right)$. For example:

$$\begin{array}{cc}\hfill N_1\left(t\right)& =-{\int }_{\Omega }w·(w·\nabla )udx,\hfill \\ \hfill N_2\left(t\right)& =-{\int }_{\Omega }w·(u_L·\nabla )wdx+{\int }_{\Omega }z·(z·\nabla )udx,\hfill \\ \hfill N_3\left(t\right)& ={\int }_{\Omega }w·(z·\nabla )Bdx-{\int }_{\Omega }z·(w·\nabla )Bdx,\hfill \\ \hfill N_4\left(t\right)& ={\int }_{\Omega }w·(B_L·\nabla )zdx-{\int }_{\Omega }z·(u_L·\nabla )zdx+{\int }_{\Omega }z·(B_L·\nabla )wdx.\hfill \end{array}$$

Finally, the time derivative of the deviation energy can be expressed as:

$${\displaystyle \frac{d}{dt}}{E}_w\left(t\right)=-\left({\parallel \nabla w\parallel }_{L^2}^2+{\parallel \nabla z\parallel }_{L^2}^2\right)+\sum _{i=1}^4{N}_i\left(t\right),$$

(5.1.4)

where $N_i\left(t\right)$ are nonlinear cross terms.

Using Hölder’s inequality and Gagliardo-Nirenberg interpolation, we estimate each cross term one by one:

• For $N_1\left(t\right)=-{\int }_{\Omega }w·(w·\nabla )udx$, the absolute value estimate is:

  $$\begin{array}{cc}\hfill |N_1\left(t\right)|& =\left|{\int }_{\Omega }w·(w·\nabla )udx\right|\hfill \\ \hfill & \le {\parallel w\parallel }_{L^2}·{\parallel (w·\nabla )u\parallel }_{L^2}\hfill \\ \hfill & \le {\parallel w\parallel }_{L^2}·{\parallel w\parallel }_{L^6}·{\parallel \nabla u\parallel }_{L^3}\hfill \\ \hfill & \le {\parallel w\parallel }_{L^2}·{C\parallel \nabla w\parallel }_{L^2}·{C\parallel u\parallel }_{L^2}^{\alpha }{\parallel \nabla u\parallel }_{L^2}^{1-\alpha }\hfill \\ \hfill & \le {C\parallel \nabla w\parallel }_{L^2}·{\left(C{(1+t)}^{-{\displaystyle \frac{k}{2}}}\right)}^{\alpha }·{\parallel \nabla u\parallel }_{L^2}^{1-\alpha }·\sqrt{2E_w\left(t\right)}\hfill \\ \hfill & \le C{(1+t)}^{-{\displaystyle \frac{k\alpha }{2}}}·{\parallel \nabla w\parallel }_{L^2}·{\parallel \nabla u\parallel }_{L^2}^{1-\alpha }·\sqrt{E_w\left(t\right)}.\hfill \end{array}$$
• For the cross term $N_2\left(t\right)=-{\int }_{\Omega }w·(u_L·\nabla )wdx$ involving $u_L$, the absolute value estimate is:

  $$\begin{array}{cc}\hfill |N_2\left(t\right)|& =\left|{\int }_{\Omega }w·(u_L·\nabla )wdx\right|\hfill \\ \hfill & \le {\parallel w\parallel }_{L^2}·{\parallel (u_L·\nabla )w\parallel }_{L^2}\hfill \\ \hfill & \le {\parallel w\parallel }_{L^2}·\parallel u_L{\parallel }_{L^6}·{\parallel \nabla w\parallel }_{L^3}\hfill \\ \hfill & \le {\parallel w\parallel }_{L^2}·C{(1+t)}^{-{\displaystyle \frac{1}{2}}}·{C\parallel w\parallel }_{L^2}^{\alpha }{\parallel \nabla w\parallel }_{L^2}^{1-\alpha }\hfill \\ \hfill & \le C{(1+t)}^{-{\displaystyle \frac{1}{2}}}·{\parallel \nabla w\parallel }_{L^2}^{1-\alpha }·{\parallel w\parallel }_{L^2}^{1+\alpha }\hfill \\ \hfill & \le C{(1+t)}^{-{\displaystyle \frac{1}{2}}}·\left({\parallel \nabla w\parallel }_{L^2}^2+E_w\left(t\right)\right).\hfill \end{array}$$
• For $N_3\left(t\right)={\int }_{\Omega }w·(z·\nabla )Bdx$ (specific form depends on deviation equations), the absolute value estimate is:

  $$\begin{array}{cc}\hfill |N_3\left(t\right)|& =\left|{\int }_{\Omega }w·(z·\nabla )Bdx\right|\hfill \\ \hfill & \le {\parallel w\parallel }_{L^2}·{\parallel (z·\nabla )B\parallel }_{L^2}\hfill \\ \hfill & \le {\parallel w\parallel }_{L^2}·{\parallel z\parallel }_{L^6}·{\parallel \nabla B\parallel }_{L^3}\hfill \\ \hfill & \le {\parallel w\parallel }_{L^2}·{C\parallel \nabla z\parallel }_{L^2}·{C\parallel B\parallel }_{L^2}^{\alpha }{\parallel \nabla B\parallel }_{L^2}^{1-\alpha }\hfill \\ \hfill & \le {C\parallel \nabla z\parallel }_{L^2}·{\left(C{(1+t)}^{-{\displaystyle \frac{k}{2}}}\right)}^{\alpha }·{\parallel \nabla B\parallel }_{L^2}^{1-\alpha }·\sqrt{2E_w\left(t\right)}\hfill \\ \hfill & \le C{(1+t)}^{-{\displaystyle \frac{k\alpha }{2}}}·{\parallel \nabla z\parallel }_{L^2}·{\parallel \nabla B\parallel }_{L^2}^{1-\alpha }·\sqrt{E_w\left(t\right)}.\hfill \end{array}$$
• For $N_4\left(t\right)={\int }_{\Omega }z·(B_L·\nabla )wdx$, the absolute value estimate is:

  $$\begin{array}{cc}\hfill |N_4\left(t\right)|& =\left|{\int }_{\Omega }z·(B_L·\nabla )wdx\right|\hfill \\ \hfill & \le {\parallel z\parallel }_{L^2}·{\parallel (B_L·\nabla )w\parallel }_{L^2}\hfill \\ \hfill & \le {\parallel z\parallel }_{L^2}·\parallel B_L{\parallel }_{L^6}·{\parallel \nabla w\parallel }_{L^3}\hfill \\ \hfill & \le {\parallel z\parallel }_{L^2}·C{(1+t)}^{-{\displaystyle \frac{1}{2}}}·{C\parallel w\parallel }_{L^2}^{\alpha }{\parallel \nabla w\parallel }_{L^2}^{1-\alpha }\hfill \\ \hfill & \le C{(1+t)}^{-{\displaystyle \frac{1}{2}}}·{\parallel \nabla w\parallel }_{L^2}^{1-\alpha }·{\parallel z\parallel }_{L^2}^{1+\alpha }\hfill \\ \hfill & \le C{(1+t)}^{-{\displaystyle \frac{1}{2}}}·\left({\parallel \nabla w\parallel }_{L^2}^2+E_w\left(t\right)\right).\hfill \end{array}$$

From Equation (4.1.5), the decay rates of the weak solution $(u,B)$ satisfy ${\parallel u\left(t\right)\parallel }_{L^2},{\parallel B\left(t\right)\parallel }_{L^2}\le C{(1+t)}^{-{\displaystyle \frac{k}{2}}}$ (where $k=min\{1,{\displaystyle \frac{3}{2}}\}$). Combining the decay property of the linear semigroup solution $\parallel u_L{\left(t\right)\parallel }_{L^q},{\parallel B_L\left(t\right)\parallel }_{L^q}\le C{(1+t)}^{-{\displaystyle \frac{3}{2}}\left({\displaystyle \frac{1}{2}}-{\displaystyle \frac{1}{q}}\right)}$ ($q\ge 2$), the sum of the four cross terms can be uniformly estimated as:

$$\sum _{i=1}^4\left|N_i\left(t\right)\right|\le C{(1+t)}^{-\beta }\left({\parallel \nabla w\parallel }_{L^2}^2+{\parallel \nabla z\parallel }_{L^2}^2+E_w\left(t\right)\right),$$

(5.1.5)

where $\beta =k+ϵ$($ϵ>0$ is a small constant). Since $k=min\{1,{\displaystyle \frac{3}{2}}\}$, we have $\beta >k$, meaning the decay rate of the cross terms is strictly faster than that of the weak solution itself.

From the time derivative formula of the deviation energy (5.1.4):

$${\displaystyle \frac{d}{dt}}{E}_w\left(t\right)=-\left({\parallel \nabla w\parallel }_{L^2}^2+{\parallel \nabla z\parallel }_{L^2}^2\right)+\sum _{i=1}^4{N}_i\left(t\right)$$

combining with the uniform estimate of cross terms (5.1.5):

$$\sum _{i=1}^4\left|N_i\left(t\right)\right|\le C{(1+t)}^{-\beta }\left({\parallel \nabla w\parallel }_{L^2}^2+{\parallel \nabla z\parallel }_{L^2}^2+E_w\left(t\right)\right),$$

we obtain the inequality:

$${\displaystyle \frac{d}{dt}}{E}_w\left(t\right)+\left({\parallel \nabla w\parallel }_{L^2}^2+{\parallel \nabla z\parallel }_{L^2}^2\right)\le C{(1+t)}^{-\beta }\left({\parallel \nabla w\parallel }_{L^2}^2+{\parallel \nabla z\parallel }_{L^2}^2+E_w\left(t\right)\right)$$

(5.1.6)

For $t\ge 1$, ${(1+t)}^{-\beta }\le 1$ (since $\beta >0$), so Equation (5.1.6) simplifies to:

$${\displaystyle \frac{d}{dt}}{E}_w\left(t\right)+\left({\parallel \nabla w\parallel }_{L^2}^2+{\parallel \nabla z\parallel }_{L^2}^2\right)\le C^{\prime }\left({\parallel \nabla w\parallel }_{L^2}^2+{\parallel \nabla z\parallel }_{L^2}^2+E_w\left(t\right)\right)$$

(5.1.7)

Where $C^{\prime }=C{(1+t)}^{-\beta }$. Rearranging the gradient terms:

$${\displaystyle \frac{d}{dt}}{E}_w\left(t\right)\le (C^{\prime }-1)\left({\parallel \nabla w\parallel }_{L^2}^2+{\parallel \nabla z\parallel }_{L^2}^2\right)+C^{\prime }{E}_w\left(t\right)$$

(5.1.8)

When t is sufficiently large, $C^{\prime }\ll 1$, so $C^{\prime }-1<0$. By the interpolation relation in Lemma 1: for divergence-free fields w and z, there exists a constant $C>0$ such that ${\parallel \nabla w\parallel }_{L^2}^2\ge C{E}_w\left(t\right)$ and ${\parallel \nabla z\parallel }_{L^2}^2\ge C{E}_w\left(t\right)$. Thus:

$${\parallel \nabla w\parallel }_{L^2}^2+{\parallel \nabla z\parallel }_{L^2}^2\ge C{E}_w\left(t\right),$$

Combining the faster decay rate of cross terms, we further simplify to:

$${\displaystyle \frac{d}{dt}}{E}_w\left(t\right)\le -C{(1+t)}^{-\beta }{E}_w\left(t\right).$$

Separating variables and integrating the differential inequality:

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{E}_w\left(t\right)}{E_w\left(t\right)}}& \le -C{(1+t)}^{-\beta }dt\hfill \\ \hfill {\int }_{E_w\left(1\right)}^{E_w\left(t\right)}{\displaystyle \frac{d\eta }{\eta }}& \le -C{\int }_1^t{(1+s)}^{-\beta }ds\hfill \\ \hfill ln{\displaystyle \frac{E_w\left(t\right)}{E_w\left(1\right)}}& \le -C·{\displaystyle \frac{{(1+t)}^{1-\beta }-2^{1-\beta }}{1-\beta }}\hfill \\ \hfill ln{E}_w\left(t\right)& \le ln{E}_w\left(1\right)+{\displaystyle \frac{C(2^{\beta -1}-{(1+t)}^{\beta -1})}{\beta -1}}\hfill \\ \hfill E_w\left(t\right)& \le E_w\left(1\right)exp\left({\displaystyle \frac{C{2}^{\beta -1}}{\beta -1}}\right)exp\left(-{\displaystyle \frac{C{(1+t)}^{\beta -1}}{\beta -1}}\right)\hfill \\ \hfill & \le C{(1+t)}^{-(\beta -1)}\hfill \end{array}$$

Since $\beta =k+ϵ$, we have $\beta -1=k+ϵ-1$. For $k\ge 1,\beta -1\ge ϵ$, so:

$${\parallel w\left(t\right)\parallel }_{L^2}+{\parallel z\left(t\right)\parallel }_{L^2}\le C{(1+t)}^{-{\displaystyle \frac{k-ϵ}{2}}}$$

Taking ${ϵ}^{\prime }={\displaystyle \frac{ϵ}{2}}$, we obtain the conclusion of the theorem: the decay rate of the deviations $(w,z)$ is strictly faster than that of the weak solution $(u,B)$, verifying the dominant role of the linear semigroup solution in the large-time behavior. □