Mixed FEM for the LL-Bar and the regularised LL-Bloch equations in micromagnetics]Mixed finite element methods for the Landau–Lifshitz–Baryakhtar and the regularised Landau–Lifshitz–Bloch equations in micromagnetics

Agus L. Soenjaya]Agus L. Soenjaya
School of Mathematics and Statistics, The University of New South Wales, Sydney 2052, Australia

65M12, 65M60, 35Q60

The Landau–Lifshitz–Baryakhtar (LLBar) and the Landau–Lifshitz–Bloch (LLBloch) equations are nonlinear vector-valued PDEs which arise in the theory of micromagnetics to describe the dynamics of magnetic spin field in a ferromagnet at elevated temperatures. We consider the LLBar and the regularised LLBloch equations in a unified manner, thus allowing us to treat the numerical approximations for both problems at once. In this paper, we propose a semi-discrete mixed finite element scheme and two fully discrete mixed finite element schemes based on a semi-implicit Euler method and a semi-implicit Crank–Nicolson method to solve the problems. These numerical schemes provide accurate approximations to both the magnetisation vector and the effective magnetic field. Moreover, they are proven to be unconditionally energy-stable and preserve energy dissipativity of the system at the discrete level. Error analysis is performed which shows optimal rates of convergence in $\bb{L}^2$, $\bb{L}^\infty$, and $\bb{H}^1$ norms. These theoretical results are further corroborated by several numerical experiments.

§ INTRODUCTION

Micromagnetics is a field of physics that deals with magnetic behaviours at sub-micrometre length scales. The length scales considered are large enough for small-scale atomic structure to be ignored, but still fine enough to resolve some local behaviour. The standard model to describe time evolution of the magnetic configuration of a ferromagnet was proposed by Landau and Lifshitz <cit.>, and is commonly known as the Landau–Lifshitz (LL) equation. According to this model, the dynamics of the magnetic spin field is influenced by a precessional term (which tends to cause precession of the magnetisation vector) and a damping term (which tends to dissipate energy and align the magnetisation vector with the effective magnetic field). The damping term is chosen purely for phenomenological reasons, namely that the damping process must lead to the state with minimum energy, and that the magnetisation magnitude must remain constant.

The Landau–Lifshitz equation has been remarkably successful in describing the magnetisation dynamics at low temperature, and was widely analysed in the physics and mathematics literature, see <cit.> and references therein. However, it is not sufficient to describe many experimental observations in modern physics, especially at high temperature <cit.>, where the magnitude of the magnetisation vector is known to be varying in time. Various magnetic recording devices, including the heat-assisted magnetic recording (HAMR) and the thermally-assisted magnetic random access memory (TA-MRAM), operate at high temperatures <cit.>, thus an accurate model and simulation of magnetisation dynamics at elevated temperatures is necessary. Indeed, the Landau–Lifshitz equation is essentially a zero-temperature equation, since it cuts off all contributions from high-frequency spin waves responsible for longitudinal magnetisation dynamics <cit.>.

To rectify these problems, several approaches are proposed in the physics literature, most prominently due to Baryakhtar <cit.> and Garanin <cit.>. Baryakhtar formulates his model based on the Onsager principles and the laws of thermodynamics, resulting in the Landau–Lifshitz–Baryakhtar (LLBar) equation, which also takes into account long-range interactions. Garanin based his model on the thermal averaging of many exchange-coupled atomistic spins, resulting in the Landau–Lifshitz–Bloch (LLBloch) equation, which can be seen as an interpolation between the standard Landau–Lifshitz equation and the Ginzburg–Landau system. In both models, the dynamics of a magnetisation vector in an effective field is influenced by three factors: precession, damping, and torque. The net effect on the magnetic spin field is a competition between these factors.

Mathematically, the dynamics of the magnetic spin field $\bff{u}(t):\mathscr{D}\to\bb{R}^3$ in a bounded magnetic domain $\mathscr{D}$ can be described as

\begin{align*}
\partial_t \bff{u} = \mathcal{P}(\bff{u}) + \mathcal{D}(\bff{u}) + \mathcal{S}(\bff{u}),
\end{align*}

where $\mathcal{P}(\bff{u})$ is the precessional term, $\mathcal{D}(\bff{u})$ is the damping term, and $\mathcal{S}(\bff{u})$ is the torque term. Quantum mechanics dictates that the precessional term has the form $\mathcal{P}(\bff{u})= -\gamma \bff{u}\times \bff{H}$, as in the usual Landau–Lifshitz equation, where $\bff{H}(t):\mathscr{D}\to \bb{R}^3$ is the effective field. Baryakhtar proposes that the damping term should be proportional to the effective field and its second spatial derivatives <cit.>. In most ferromagnetic materials, the damping term for the LLBar equation [13] can be written as

\begin{equation}\label{equ:damping}
\mathcal{D}_{\text{Bar}}(\bff{u}) = \lambda_r \bff{H} - \lambda_e \Delta\bff{H},
\end{equation}

where $\bff{H}$ itself depends on $\bff{u}$. The positive constants $\lambda_r$ and $\lambda_e$ are the phenomenological relativistic damping constant and the exchange damping constant, respectively. These constants can also be replaced by positive definite symmetric tensors without changing the arguments here substantially. The second term in (<ref>) is responsible for the longitudinal damping and long range interaction between magnetic spins [13].
For the LLBloch equation above $T_\mathrm{c}$, the damping term is simply (<ref>) with $\lambda_e=0$.

The torque term has various forms depending on the physical situations considered. In the presence of applied current, the flow of electrons move charges and spins across space, affecting the magnetic properties of the material. A commonly used torque term is given by Zhang and Li <cit.>:

\begin{equation}\label{equ:spin}
\mathcal{S}(\bff{u})= \Lambda_1 \bff{u}\times (\bff{u}\times (\bff{j}\cdot \nabla)\bff{u}) + \Lambda_2 \bff{u}\times (\bff{j}\cdot \nabla)\bff{u},
\end{equation}

where $\Lambda_1$ and $\Lambda_2$ are constants (which can be positive or negative), and $\bff{j}=\bff{j}(t,\bff{x})$ is the current density.

The effective magnetic field $\bff{H}$ consists of the usual contributions from the theory of micromagnetism, namely the exchange field, the external (Zeeman) field, and the anisotropy field. To account for elevated temperatures and phase transition, the Ginzburg–Landau internal exchange field is added <cit.>, giving

\begin{equation}\label{equ:eff field}
\bff{H}(t,\bff{x})
= \underbrace{\alpha \Delta \bff{u}}_{\text{exchange}} + \underbrace{\kappa \mu \bff{u}- \kappa |\bff{u}|^2 \bff{u}}_{\text{internal exchange}} + \underbrace{\bff{B}(t)}_{\text{applied}}
- \underbrace{\beta \bff{e}(\bff{e}\cdot \bff{u})}_{\text{anisotropy}}.
\end{equation}

Here, $\alpha$ is a positive constant depending on the material structure. The unit vector $\bff{e}$ is related to the axis of anisotropy of the material, while the vector $\bff{B}(t)$ describes an applied magnetic field. Physical considerations dictate that $\kappa$ is positive, while $\mu>0$ for temperatures above the Curie temperature $T_\mathrm{c}$ and $\mu<0$ below $T_\mathrm{c}$. The constant $\beta$ corresponds to the uniaxial anisotropy constant of the object, which can be positive or negative. To simplify presentation, we assume $\beta<0$ in (<ref>), which corresponds to the presence of an easy plane. Higher order anisotropy field of the form $\beta_1 \bff{e}(\bff{e}\cdot\bff{u})- \beta_2 \bff{e}(\bff{e}\cdot \bff{u})^3$, where $\beta_1,\beta_2>0$, could also be treated in a similar manner as the internal exchange field.

Subsequently, we take $\bff{B}(t)=\mathcal{S}(\bff{u})\equiv \bff{0}$ and set $\alpha=1$ for simplicity. The presence of applied field and the Zhang–Li torque term (<ref>) could still be handled by the methods presented here with minor modifications. However, we choose to focus on the simplified version since in this case the system dissipates energy, and we want to highlight this energy dissipation property in our numerical schemes.

With these in mind, we can state the LLBar and the LLBloch equations in a unified manner. Given a magnetic body $\mathscr{D}\subset \bb{R}^d$, the magnetisation vector $\bff{u}(t,\bff{x})$ for any time $t\geq 0$ and at any point $\bff{x}\in \mathscr{D}$ evolves according to a nonlinear vector-valued PDE, which can be written in the following mixed form:

\begin{equation}\label{equ:llbar}
\begin{cases}
\displaystyle
\partial_t \bff{u}
\lambda_r \bff{H}
- \lambda_e \Delta \bff{H}
- \gamma \bff{u} \times \bff{H},
\quad & \text{for $(t,\bff{x})\in(0,T)\times\mathscr{D}$,}
\\[1mm]
\bff{H}
\Delta \bff{u}
+ \kappa \mu \bff{u}
- \kappa |\bff{u}|^2 \bff{u}
- \beta \bff{e}(\bff{e}\cdot \bff{u}),
\quad & \text{for $(t,\bff{x})\in(0,T)\times\mathscr{D}$,}
\\[1mm]
\bff{u}(0,\bff{x})= \bff{u_0}(\bff{x})
\quad & \text{for } \bff{x}\in \mathscr{D},
\\[1mm]
\displaystyle{
\frac{\partial \bff{u}}{\partial \bff{n}}= \bff{0}},
\;\displaystyle{\frac{\partial \bff{H}}{\partial \bff{n}}= \bff{0}}
\quad & \text{for } (t,\bff{x})\in (0,T) \times \partial \mathscr{D},
\end{cases}
\end{equation}

where $\partial\mathscr{D}$ is the boundary of $\mathscr{D}$, with exterior unit normal vector denoted by $\bff{n}$. Depending on the signs of some parameters, problem (<ref>) describes several models:

* If $\lambda_e>0$ and $\mu>0$, this is the LLBar equation below $T_\mathrm{c}$.

* If $\lambda_e=0$ and $\mu<0$, this is the LLBloch equation above $T_\mathrm{c}$.

* If $\lambda_e>0$ and $\mu<0$, this is the LLBar equation above $T_\mathrm{c}$, which can also be considered as the regularised LLBloch equation above $T_\mathrm{c}$ (cf. section <ref>).

We remark that the LLBloch equation below $T_\mathrm{c}$ cannot be written in the above form and will be analysed separately in an upcoming paper, although for temperatures not far from (but below) $T_\mathrm{c}$, the LLBar equation with $\lambda_e>0$ and $\mu>0$ can be considered as an approximation to the LLBloch equation below $T_\mathrm{c}$ [13].

The energy $\mathcal{E}(\bff{u})$ of the system described by (<ref>) is

\begin{equation}\label{equ:energy}
\mathcal{E}(\bff{u})
\begin{cases}
\displaystyle
\frac{1}{2} \norm{\nabla \bff{u}}{\bb{L}^2}^2
\frac{\kappa}{4} \norm{|\bff{u}|^2-\mu}{\bb{L}^2}^2
\frac{\beta}{2} \int_\mathscr{D} (\bff{e}\cdot \bff{u})^2\, \dx,
\quad & \text{if } \mu>0,
\\[3mm]
\displaystyle
\frac{1}{2} \norm{\nabla \bff{u}}{\bb{L}^2}^2
\frac{\kappa}{4} \norm{\bff{u}}{\bb{L}^4}^4
\frac{\kappa\mu}{2} \norm{\bff{u}}{\bb{L}^2}^2
\frac{\beta}{2} \int_\mathscr{D} (\bff{e}\cdot \bff{u})^2\, \dx,
\quad & \text{if } \mu< 0,
\end{cases}
\end{equation}

where $\kappa, \mu, \beta$ are constants appearing in (<ref>). Note that the two expressions in (<ref>) differ only by a constant $\frac{1}{4} \kappa \mu^2 |\mathscr{D}|$, and in principle one could take the first expression as the energy for any values of $\mu$. However, it is known <cit.> that $\bff{u}(t)\to \bff{0}$ as $t\to\infty$ for $\mu< 0$, and so the energy in this case is modified to ensure $\mathcal{E}(\bff{0})= 0$. In the absence of the spin-torque term, this system dissipates energy.

Global well-posedness of the LLBar equation (assuming exchange-dominated field) is shown in [11], while the existence and uniqueness of strong solution to the LLBloch equation are shown in <cit.> (see also <cit.> and <cit.> for the stochastic equations). On the numerical aspect, some linear $C^1$-conforming finite element methods based on a semi-implicit Euler and a BDF schemes are proposed to solve the LLBar equation <cit.>. They are proven to be stable and convergent to the actual solution at an optimal rate, without assuming quasi-uniformity of the triangulation in most cases. However, the implementation is computationally costly, since the method requires $C^1$-continuity across element boundaries. A linear conforming finite element scheme based on the linearised Euler method has been proposed for the regularised LLBloch equation with a different regularisation term, where an optimal convergence rate in $\bb{H}^1$ is shown <cit.>. Other regularisation term is used in <cit.>, where a $C^1$-conforming method is proposed for the stochastic regularised LLBloch equation, but a strong convergence rate is only provided in $\bb{L}^2$ and the results are limited to $d=1$ or $2$. In all the above, energy dissipativity of the schemes is not addressed.

Here, we continue and build on the studies done in <cit.> and <cit.> by proposing numerical schemes which can be used to solve the LLBar and the LLBloch equations for $d=1,2,3$ in a unified manner. More precisely, we propose some mixed finite element methods to solve the LLBar and the regularised LLBloch equations, including a semi-discrete (in space) finite element scheme and two fully-discrete finite element schemes based on the Euler and the Crank–Nicolson methods. These schemes are proven to be uniquely solvable and unconditionally energy-stable, and this $\bb{H}^1$-stability is robust with respect to the parameter $\lambda_e$. Moreover, the energy dissipativity of the system is preserved at the discrete level. Since we are using the mixed formulation (<ref>), these methods provide an accurate approximation to both the magnetisation vector $\bff{u}$ and the effective magnetic field $\bff{H}$, both of which are physically significant quantities.

To summarise, the main results of this paper include proving the following:

* convergence of a semi-discrete conforming finite element scheme (Theorem <ref>),

* convergence of a fully discrete conforming finite element scheme based on the semi-implicit Euler method (Theorem <ref>),

* convergence of a fully discrete conforming finite element scheme based on the Crank–Nicolson method (Theorem <ref>),

* convergence of the strong solution of the LLBar equation to that of the LLBloch equation at a certain rate (in $L^\infty(0,T;\bb{H}^1(\mathscr{D})$), thus implying the above schemes can also be used to solve the LLBloch equation. Convergence of the corresponding effective field (in $L^2(0,T;\bb{L}^2(\mathscr{D})$) is also shown (Theorem <ref>).

In all cases, we obtain the expected rates of convergence in $\bb{L}^2$, $\bb{L}^\infty$, and $\bb{H}^1$ norms. We remark that systems of vector-valued PDEs similar to (<ref>) also appear in chemistry and biology to model long-range diffusion <cit.>, anomalous bi-flux diffusion <cit.>, and population dynamics <cit.>. As such, numerical schemes proposed in this paper would also apply to certain cases of these models.

This paper is organised as follows. Notations and various assumptions on the exact solution and the finite element space are outlined in Section <ref>. A semi-discrete finite element approximation is described in Section <ref>. Two fully discrete mixed finite element schemes based on the Euler and the Crank–Nicolson methods are proposed in Section <ref> and Section <ref>, respectively. The regularised LLBloch equation and its approximation are discussed in Section <ref>. Finally, some numerical simulations which support the theoretical results are described in Section <ref>.

§ PRELIMINARIES

§.§ Notations

We begin by defining some notations used in this paper. The function space $\bb{L}^p := \bb{L}^p(\mathscr{D}; \bb{R}^3)$ denotes the usual space of $p$-th integrable functions taking values in $\bb{R}^3$ and $\bb{W}^{k,p} := \bb{W}^{k,p}(\mathscr{D}; \bb{R}^3)$ denotes the usual Sobolev space of
functions on $\mathscr{D} \subset \bb{R}^d$, for $d=1,2,3$, taking values in $\bb{R}^3$. Also, we write $\bb{H}^k := \bb{W}^{k,2}$. Let $\Delta$ be the Neumann Laplacian operator acting on $\bb{R}^3$-valued functions with domain $\mathrm{D}(\Delta)$ given by

\begin{equation*}
\mathrm{D}(\Delta):= \left\{\bff{v}\in \bb{H}^2 : \frac{\partial \bff{v}}{\partial \bff{n}} = \bff{0} \text{ on } \partial\mathscr{D} \right\}.
\end{equation*}

If $X$ is a Banach space, the spaces $L^p(0,T;X)$ and $W^{m,r}(0,T;X)$ denote respectively the Lebesgue and Sobolev spaces of functions on $(0,T)$ taking values in $X$, where $T$ can be a finite number or $\infty$. The space $C([0,T];X)$ denotes the space of continuous function on $[0,T]$ taking values in $X$. For simplicity, we will write $L^p(\bb{W}^{m,r}) := L^p(0,T; \bb{W}^{m,r})$ and $L^p(\bb{L}^q) := L^p(0,T; \bb{L}^q)$.

Throughout this paper, we denote the scalar product in a Hilbert space $H$ by $\inpro{ \cdot}{ \cdot}_H$ and its corresponding norm by $\norm{\cdot}{H}$. We will not distinguish between the scalar product of $\bb{L}^2$ vector-valued functions taking values in $\bb{R}^3$ and the scalar product of $\bb{L}^2$ matrix-valued functions taking values in $\bb{R}^{3\times 3}$, and still denote them by $\inpro{\cdot}{\cdot}$.

Here, $\mathscr{D}$ is assumed to be a bounded domain such that the $\bb{H}^2$-regularity result holds, namely:

\begin{align*}
\norm{\bff{v}}{\bb{H}^2}^2
\norm{\bff{v}}{\bb{L}^2}^2 + \norm{\Delta \bff{v}}{\bb{L}^2}^2
\end{align*}

for all $\bff{v}\in \mathrm{D}(\Delta)$. The above is true for any domain with $C^2$-smooth boundary. For domain with less regular boundary, it is known, for instance, that the $\bb{H}^2$-regularity result is guaranteed to hold for any convex Lipschitz domains [5]. Henceforth, we will assume that $\mathscr{D}$ is a smooth, or convex polygonal or polyhedral domain. In this domain, the Gagliardo–Nirenberg inequalities and the Sobolev embedding theorems also hold.

Finally, throughout this paper, the constant $C$ in the estimate denotes a
generic constant which takes different values at different occurrences. If
the dependence of $C$ on some variable, e.g. $T$, is highlighted, we often write
$C(T)$. The notation $A \lesssim B$ means $A \le C B$ where the specific form of
the constant $C$ is not important to clarify.

§.§ Assumptions

Given $T>0$ and $\bff{u}_0\in\bb{H}^1(\mathscr{D})$, the pair of functions $(\bff{u}, \bff{H}):[0,T]\to \bb{H}^1 \times \bb{H}^1$ is a weak solution to the problem
(<ref>) if $(\bff{u}, \bff{H}) $ satisfies

\begin{align}\label{equ:weakform}
\left\{
\begin{alignedat}{1}
\inpro{ \partial_t \bff{u}}{\bff{\chi}}
\lambda_r \inpro{\bff{H}}{\bff{\chi}}
\lambda_e \inpro{\nabla \bff{H}}{\nabla \bff{\chi}}
\gamma \inpro{\bff{u} \times \bff{H}}{\bff{\chi}}
\\
\inpro{\bff{H}}{\bff{\phi}}
\inpro{\nabla \bff{u}}{\nabla \bff{\phi}}
\kappa \mu \inpro{\bff{u}}{\bff{\phi}}
\kappa \inpro{|\bff{u}|^2 \bff{u}}{\bff{\phi}}
\beta \inpro{\bff{e}(\bff{e}\cdot \bff{u})}{\bff{\phi}},
\end{alignedat}
\right.
\end{align}

for all $\bff{\chi}, \bff{\phi}\in \bb{H}^1$ and $t\in [0,T]$, with $\bff{u}(0)=\bff{u}_0$.

Throughout this paper, we assume that problem (<ref>) possesses
solution $(\bff{u}, \bff{H})$ which satisfies

\begin{equation}\label{equ:ass 2 u}
\begin{aligned}
\norm{\bff{u}}{L^\infty(\bb{H}^{r+1})}
+ \norm{\partial_t \bff{u}}{L^\infty(\bb{H}^{r+1})}
+ \norm{\partial_t^2 \bff{u}}{L^\infty(\bb{H}^2)}
+ \norm{\partial_t^3 \bff{u}}{L^\infty(\bb{L}^2)}
&\le K_0,
\\
\norm{\bff{H}}{L^\infty(\bb{H}^{r+1})}
+ \norm{\partial_t \bff{H}}{L^\infty(\bb{H}^{r+1})}
+ \norm{\partial_t^2 \bff{H}}{L^\infty(\bb{H}^2)}
+ \norm{\partial_t^3 \bff{H}}{L^\infty(\bb{L}^2)}
&\le K_0,
\end{aligned}
\end{equation}

where $r$ is the degree of piecewise continuous polynomials used as the finite element space and $K_0>0$ depends on $\bff{u}_0$. The existence of an arbitrarily smooth solution to the LLBar equation on $(0,T)\times\mathscr{D}$ is guaranteed for any initial data $\bff{u_0}\in\bb{H}^1$ (cf. <cit.>).

For simplicity of presentation, throughout this paper we assume that $\mu>0$, except in Section <ref> where we specifically discuss the regularised LLBloch equation and outline the modifications needed for the case $\mu<0$, and in Section <ref> where some numerical simulations are performed.

§.§ Finite Element Approximation

Let $\mathcal{T}_h$ be a shape-regular triangulation of $\mathscr{D}\subset \bb{R}^d$ with maximal mesh-size $h$. To discretise the LLBar equation (<ref>), we introduce the finite element space $\bb{V}_h \subset \bb{H}^1$, which is the space of all piecewise continuous polynomials on $\mathcal{T}_h$ of degree at most $r$. By the Bramble–Hilbert lemma, there exists a constant $C$ independent of $h$ such that for any $\bff{v} \in \bb{H}^{r+1}$,

\begin{align*}
\inf_{\chi \in {\bb{V}}_h} \left\{ \norm{\bff{v} - \bff{\chi}}{\bb{L}^2}
h \norm{\nabla (\bff{v}-\bff{\chi})}{\bb{L}^2}
\right\}
\leq
C h^{r+1} \norm{\bff{v}}{\bb{H}^{r+1}}.
\end{align*}

We shall use several operators in the analysis. Firstly, define the $\bb{L}^2$-projection operator $\Pi_h: \bb{L}^2 \to \bb{V}_h$ such that

\begin{align}\label{equ:orth proj}
\inpro{\Pi_h \bff{v}-\bff{v}}{\bff{\chi}}=0,
\quad
\forall \bff{\chi}\in \bb{V}_h.
\end{align}

It is well known that if $\bff{v}\in \bb{H}^{r+1}$, then

\begin{align}\label{equ:proj approx}
\norm{\bff{v}- \Pi_h\bff{v}}{\bb{L}^2}
h \norm{\nabla( \bff{v}-\Pi_h\bff{v})}{\bb{L}^2}
\leq
C h^{r+1} \norm{\bff{v}}{\bb{H}^{r+1}}.
\end{align}

Moreover, if the triangulation is globally quasi-uniform, then we have the $\bb{H}^1$-stability of the $\bb{L}^2$-projection operator <cit.>, namely

\begin{align}\label{equ:H1 stab proj}
\norm{\nabla \Pi_h \bff{v}}{\bb{L}^2} \leq C \norm{\nabla \bff{v}}{\bb{L}^2}, \quad \forall \bff{v}\in \bb{H}^1.
\end{align}

Next, we introduce the discrete Laplacian operator $\Delta_h: \bb{V}_h \to \bb{V}_h$ defined by

\begin{align}\label{equ:disc laplacian}
\inpro{\Delta_h \bff{v}_h}{\bff{\chi}}
- \inpro{\nabla \bff{v}_h}{\nabla \bff{\chi}},
\quad
\forall \bff{v}_h, \bff{\chi} \in \bb{V}_h,
\end{align}

and the Ritz projection operator $R_h: \bb{H}^1 \to \bb{V}_h$ by

\begin{align}\label{equ:Ritz}
\inpro{\nabla R_h \bff{v}- \nabla \bff{v}}{\nabla \bff{\chi}}=0
\;\text{ such that }
\inpro{R_h\bff{v}-\bff{v}}{\bff{1}}=0,
\quad
\forall \bff{\chi}\in \bb{V}_h,
\end{align}

For any $\bff{v}\in \bb{H}^{r+1}$, let $\bff{\omega}(t):= \bff{v}(t)-R_h\bff{v}(t)$. The following estimate is well known [12] for $s=0$ or $1$:

\begin{align}\label{equ:Ritz ineq}
\norm{\bff{\omega}(t)}{\bb{H}^s} + \norm{\partial_t \bff{\omega}(t)}{\bb{H}^s}
\leq
C h^{r+1-s} \norm{\bff{v}(t)}{\bb{H}^{r+1}}.
\end{align}

Moreover, if the triangulation is globally quasi-uniform, then by [10],

\begin{align}\label{equ:Ritz ineq L infty}
\norm{\bff{\omega}(t)}{\bb{L}^\infty} \leq Ch^{r+1} \abs{\ln h} \norm{\bff{v}(t)}{\bb{W}^{r+1,\infty}}.
\end{align}

It is known that the term $\abs{\ln h}$ can be removed in case $r\geq 2$ and $\mathscr{D}$ is a polygonal domain.
Throughout, we shall assume sufficient conditions on the regularity of the domain and the geometry of the mesh so that the maximum-norm stability of the Ritz projection holds, namely

\begin{equation}\label{equ:Ritz stab infty}
\norm{R_h \bff{v}}{\bb{X}}
\leq
C \norm{\bff{v}}{\bb{X}},
\quad \text{where either } \bb{X}:= \bb{L}^\infty \text{ or } \bb{W}^{1,\infty}.
\end{equation}

This holds for a globally quasi-uniform triangulation [4]. However, (<ref>) also holds under more general conditions, for instance in a convex polygonal or polyhedral domain with mildly graded mesh satisfying certain assumptions [4, 8], or in a non-convex polygonal domain with locally refined mesh <cit.>.

§.§ Auxiliary Results

In the analysis, we use the following vector identities: for any vectors $\bff{a}$, $\bff{b}\in \bb{R}^3$,

\begin{align}\label{equ:a dot ab}
2\bff{a}\cdot (\bff{a}-\bff{b})
\abs{\bff{a}}^2 - \abs{\bff{b}}^2 + \abs{\bff{a}-\bff{b}}^2,
\\
\label{equ:a2a dot ab}
4\abs{\bff{a}}^2 \bff{a} \cdot (\bff{a}-\bff{b})
\abs{\bff{a}}^4 - \abs{\bff{b}}^4 + \left(\abs{\bff{a}}^2 - \abs{\bff{b}}^2\right)^2
+ 2 \abs{\bff{a}}^2 \abs{\bff{a}-\bff{b}}^2,
\\
\label{equ:a2a b2b dot ab}
2(\abs{\bff{a}}^2 \bff{a} - \abs{\bff{b}}^2 \bff{b}) \cdot (\bff{a}-\bff{b})
\left(\abs{\bff{a}}^2 - \abs{\bff{b}}^2 \right)^2 + \left(\abs{\bff{a}}^2+\abs{\bff{b}}^2\right) \abs{\bff{a}-\bff{b}}^2.
\end{align}

The following inequalities will also be used frequently.

Let $\epsilon>0$ be given. Then there exists a positive constant $C$ depending only on $\mathscr{D}$ such that the following inequalities hold.

* For any $\bff{v} \in \bb{H}^1(\mathscr{D})$,

\begin{align}
\label{equ:L4 gal nir}
\norm{\bff{v}}{\bb{L}^4}^2
C \norm{\bff{v}}{\bb{L}^2}^{2-d/2}
\norm{\bff{v}}{\bb{H}^1}^{d/2},
\\
\label{equ:L4 young}
\norm{\bff{v}}{\bb{L}^4}^2
C \norm{\bff{v}}{\bb{L}^2}^2
\epsilon \norm{\nabla \bff{v}}{\bb{L}^2}^2.
\end{align}

* For any $\bff{v} \in \mathrm{D}(\Delta)$,

\begin{align}
\label{equ:equivnorm-nabL2 young}
\norm{\nabla \bff{v}}{\bb{L}^2}^2
\frac{1}{4\epsilon} \norm{\bff{v}}{\bb{L}^2}^2
\epsilon \norm{\Delta \bff{v}}{\bb{L}^2}^2.
\end{align}

* For any $\bff{v}_h\in \bb{V}_h$,

\begin{align}
\label{equ:disc lapl L2}
\norm{\nabla \bff{v}_h}{\bb{L}^2}^2
\leq
\norm{\bff{v}_h}{\bb{L}^2} \norm{\Delta_h \bff{v}_h}{\bb{L}^2}.
\end{align}

* Let $\mathscr{D}$ be a convex polygonal or polyhedral domain with globally quasi-uniform triangulation. For any $\bff{v}_h\in \bb{V}_h$,

\begin{align}
\label{equ:disc lapl L infty}
\norm{\bff{v}_h}{\bb{L}^\infty}
C \norm{\bff{v}_h}{\bb{L}^2}^{1-\frac{d}{4}} \left(\norm{\bff{v}_h}{\bb{L}^2}^\frac{d}{4} + \norm{\Delta_h \bff{v}_h}{\bb{L}^2}^\frac{d}{4} \right),
\\
\label{equ:disc lapl L6}
\norm{\nabla \bff{v}_h}{\bb{L}^6}
C \norm{\Delta_h \bff{v}_h}{\bb{L}^2}.
\end{align}

Inequality (<ref>) follows from the Gagliardo–Nirenberg inequality. Moreover, applying Young's inequality to (<ref>) gives

\begin{align*}
\norm{\bff{v}}{\bb{L}^4}^2
C \norm{\bff{v}}{\bb{L}^2}^{2-d/2}
\big( \norm{\bff{v}}{\bb{L}^2}^{d/2}
\norm{\nabla \bff{v}}{\bb{L}^2}^{d/2} \big)
\leq
C \norm{\bff{v}}{\bb{L}^2}^2
\epsilon \norm{\nabla \bff{v}}{\bb{L}^2}^2.
\end{align*}

Inequality (<ref>) follows from integration by parts and Young's inequality, while (<ref>) follows by taking $\bff{\chi}=\bff{v}_h$ in (<ref>) and applying Hölder's inequality.

Finally, (<ref>) and (<ref>) are proven in <cit.>). This completes the proof of the lemma.

§ A SEMI-DISCRETE GALERKIN APPROXIMATION

A semi-discrete Galerkin approximation to the problem (<ref>) is $(\bff{u}_h, \bff{H}_h): [0,T] \to \bb{V}_h \times \bb{V}_h$ such that for all $t\in [0,T]$ and $\bff{\chi}, \bff{\phi} \in \bb{V}_h$,

\begin{align}\label{equ:weaksemidisc}
\left\{
\begin{alignedat}{1}
\inpro{ \partial_t \bff{u}_h}{\bff{\chi}}
\lambda_r \inpro{\bff{H}_h}{\bff{\chi}}
\lambda_e \inpro{\nabla \bff{H}_h}{\nabla \bff{\chi}}
\gamma \inpro{\bff{u}_h \times \bff{H}_h}{\bff{\chi}},
\\
\inpro{\bff{H}_h}{\bff{\phi}}
\inpro{\nabla \bff{u}_h}{\nabla \bff{\phi}}
\kappa \mu \inpro{\bff{u}_h}{\bff{\phi}}
\kappa \inpro{|\bff{u}_h|^2 \bff{u}_h}{\bff{\phi}}
\beta \inpro{\bff{e}(\bff{e}\cdot \bff{u}_h)}{\bff{\phi}},
\end{alignedat}
\right.
\end{align}

with $\bff{u}_h(0)= \bff{u}_{0,h}$, an approximation of $\bff{u_0}$ in $\bb{V}_h$.
Note that (<ref>) can be written as

\begin{align}\label{equ:weakproj}
\left\{
\begin{alignedat}{1}
\partial_t \bff{u}_h
\lambda_r \bff{H}_h
\lambda_e \Delta_h \bff{H}_h
\gamma \Pi_h (\bff{u}_h \times \bff{H}_h),
\\
\bff{H}_h
\Delta_h \bff{u}_h
\kappa \mu \bff{u}_h
\kappa \Pi_h (|\bff{u}_h|^2 \bff{u}_h)
\beta \Pi_h \big(\bff{e}(\bff{e}\cdot \bff{u}_h)\big).
\end{alignedat}
\right.
\end{align}

Substituting the second equation into the first gives

\begin{align}\label{equ:ode sub}
\nonumber
\partial_t \bff{u}_h
(\lambda_r-\kappa\lambda_e \mu) \Delta_h \bff{u}_h
\lambda_e \Delta_h^2 \bff{u}_h
\kappa \mu \lambda_r \bff{u}_h
\kappa\lambda_r \Pi_h (|\bff{u}_h|^2 \bff{u}_h)
\kappa\lambda_e \Delta_h \Pi_h (|\bff{u}_h|^2 \bff{u}_h)
\\
\gamma \Pi_h (\bff{u}_h \times \Delta_h \bff{u}_h)
\kappa \gamma \Pi_h \big( \bff{u}_h \times \Pi_h (|\bff{u}_h|^2 \bff{u}_h) \big)
\beta\gamma \Pi_h \big(\bff{u}_h\times \Pi_h \big(\bff{e}(\bff{e}\cdot \bff{u}_h)\big)\big)
\nonumber \\
\beta \lambda_r \Pi_h \big(\bff{e}(\bff{e}\cdot \bff{u}_h)\big)
\beta\lambda_e \Delta_h \Pi_h \big(\bff{e}(\bff{e}\cdot \bff{u}_h)\big).
\end{align}

Noting that all norms on the finite-dimensional space $\bb{V}_h$ are equivalent and each term in (<ref>) is locally Lipschitz, we obtain a unique solution $\bff{u}_h\in \bb{V}_h$ for (<ref>) defined on the interval $[0, t_h]\subseteq [0,T]$ by the standard theory for ordinary differential equation (thus also giving a unique solution $\bff{H}_h \in \bb{V}_h$ by substituting $\bff{u}_h$ back to (<ref>)). We will prove several stability results, which will be used to ensure the semi-discrete solution $(\bff{u}_h, \bff{H}_h)$ can be continued globally to $[0,\infty)$ for any initial data $\bff{u}_{0,h}\in \bb{V}_h$. Note that these estimates hold uniformly in time $t\in (0,\infty)$.

Let $h>0$ and initial data $\bff{u}_0$ be given. For any $s$, $t\in (0,\infty)$ such that $s\leq t$,

\begin{align}
\label{equ:semidisc-est ener}
\cal{E}(\bff{u}_h(t))
\leq
\cal{E}(\bff{u}_h(s)),
\end{align}

where $\cal{E}$ is the energy defined in (<ref>), and

\begin{align}
\label{equ:semidisc-est1}
\norm{\bff{u}_h(t)}{\bb{L}^4}^4
\norm{\bff{u}_h(t)}{\bb{H}^1}^2
\lambda_r \int_0^t \norm{\bff{H}_h(s)}{\bb{L}^2}^2 \ds
\lambda_e \int_0^t \norm{\nabla \bff{H}_h(s)}{\bb{L}^2}^2 \ds
\leq
C_1 \left(\norm{\bff{u}_0}{\bb{H}^1}^4
|\mathscr{D}| \right),
\end{align}

where $C_1$ depends only on $\kappa$, $\mu$, and $\beta$.

Taking $\bff{\chi} = \bff{H}_h$ and $\bff{\phi}= \partial_t \bff{u}_h$ in (<ref>), we obtain

\begin{align*}
\inpro{\partial_t \bff{u}_h}{\bff{H}_h}
\lambda_r \norm{\bff{H}_h}{\bb{L}^2}^2
\lambda_e \norm{\nabla \bff{H}_h}{\bb{L}^2}^2,
\\
\inpro{\bff{H}_h}{\partial_t \bff{u}_h}
\frac{1}{2} \ddt \norm{\nabla \bff{u}_h}{\bb{L}^2}^2
\frac{\kappa \mu}{2} \ddt \norm{\bff{u}_h}{\bb{L}^2}^2
\frac{\kappa}{4} \ddt \norm{\bff{u}_h}{\bb{L}^4}^4
\frac{\beta}{2} \ddt (\bff{e}\cdot \bff{u}_h)^2,
\end{align*}

These imply

\begin{align*}
\frac{1}{2} \ddt \norm{\nabla \bff{u}_h}{\bb{L}^2}^2
\frac{\kappa}{4} \ddt \norm{\bff{u}_h}{\bb{L}^4}^4
\frac{\kappa\mu}{2} \ddt \norm{\bff{u}_h}{\bb{L}^2}^2
\lambda_r \norm{\bff{H}_h}{\bb{L}^2}^2
\lambda_e \norm{\nabla \bff{H}_h}{\bb{L}^2}^2,
\end{align*}

or equivalently the energy identity:

\begin{equation*}
\ddt\, \cal{E}(\bff{u}_h(t)) + \lambda_r \norm{\bff{H}_h(t)}{\bb{L}^2}^2 + \lambda_e \norm{\nabla \bff{H}_h(t)}{\bb{L}^2}^2 = 0,
\end{equation*}

thus proving (<ref>). Next, integrating this with respect to $t$ and rearranging, we obtain

\begin{align}\label{equ:nab uh eq}
\frac{1}{2} \norm{\nabla \bff{u}_h}{\bb{L}^2}^2
\frac{1}{2}\norm{\nabla \bff{u}_h(0)}{\bb{L}^2}^2
\frac{\kappa}{4} \norm{\bff{u}_h}{\bb{L}^4}^4
\frac{\kappa}{4} \norm{\bff{u}_h(0)}{\bb{L}^4}^4
\frac{\beta}{2} \norm{\bff{e}\cdot \bff{u}_h}{\bb{L}^2}^2
\frac{\beta}{2} \norm{\bff{e}\cdot \bff{u}_h(0)}{\bb{L}^2}^2
\nonumber \\
&+ \lambda_r \int_0^t \norm{\bff{H}_h(s)}{\bb{L}^2}^2 \ds
+ \lambda_e \int_0^t \norm{\nabla \bff{H}_h(s)}{\bb{L}^2}^2 \ds
\leq
\frac{\kappa\mu}{2} \norm{\bff{u}_h}{\bb{L}^2}^2
\frac{\kappa\mu}{2} \norm{\bff{u}_h(0)}{\bb{L}^2}^2.
\end{align}

Let $A:=\kappa\mu+\beta$. After rearranging some terms, we have

\begin{align}\label{equ:norm uh L4}
\nonumber
&\norm{\nabla \bff{u}_h}{\bb{L}^2}^2
\frac{A}{2} \norm{\bff{u}_h}{\bb{L}^2}^2
\frac{\kappa}{4} \norm{\bff{u}_h}{\bb{L}^4}^4
+ \lambda_r \int_0^t \norm{\bff{H}_h(s)}{\bb{L}^2}^2 \ds
+ \lambda_e \int_0^t \norm{\nabla \bff{H}_h(s)}{\bb{L}^2}^2 \ds
\\
\nonumber
\kappa \norm{\bff{u}_h(0)}{\bb{L}^4}^4
C \norm{\bff{u}_h(0)}{\bb{H}^1}^2
\frac{3A}{2} \norm{\bff{u}_h}{\bb{L}^2}^2
\frac{\kappa}{4} \norm{\bff{u}_h}{\bb{L}^4}^4
\\
\nonumber
\kappa \norm{\bff{u}_h(0)}{\bb{L}^4}^4
C \norm{\bff{u}_h(0)}{\bb{H}^1}^2
\frac{\kappa}{4} \int_\mathscr{D} \left( \frac{6A}{\kappa} |\bff{u}_h(t)|^2 - |\bff{u}_h(t)|^4 \right) \dx
\\
\kappa \norm{\bff{u}_h(0)}{\bb{L}^4}^4
C \norm{\bff{u}_h(0)}{\bb{H}^1}^2
C |\mathscr{D}|
\leq
C \left( \norm{\bff{u}_h(0)}{\bb{H}^1}^4 + |\mathscr{D}| \right),
\end{align}

where in the last step we used (<ref>), and $C$ depends only on $\kappa, \mu, \beta$.
This shows (<ref>), thus completing the proof of the proposition.

Let $h>0$ and initial data $\bff{u}_0$ be given. For all $t\in (0,\infty)$,

\begin{align}\label{equ:semidisc-est2}
\lambda_e \norm{\bff{H}_h(t)}{\bb{L}^2}^2
\int_0^t \norm{\partial_t \bff{u}_h(s)}{\bb{L}^2}^2 \ds
\leq
C_1 \alpha_1 \left( \norm{\bff{u}_0}{\bb{H}^1}^2 + |\mathscr{D}| \right),
\end{align}

where $\alpha_1:= \lambda_r^{-1} \big((\lambda_r+\kappa \mu\lambda_e+\beta\lambda_e)^2 + C\gamma^2 + C(\kappa \lambda_e)^2\big)$, and the constant $C$ depends only on $\kappa$, $\mu$, $\beta$, and $\mathscr{D}$ (but is independent of $t$, $h$, and $\lambda_e$).

Taking $\bff{\chi}= \partial_t \bff{u}_h$ in (<ref>) yields

\begin{align}\label{equ:norm dt uh L2}
\norm{\partial_t \bff{u}_h}{\bb{L}^2}^2
\lambda_r \inpro{\bff{H}_h}{\partial_t \bff{u}_h}
\lambda_e \inpro{\nabla \bff{H}_h}{\nabla \partial_t \bff{u}_h}
\gamma \inpro{\bff{u}_h \times \bff{H}_h}{\partial_t \bff{u}_h}.
\end{align}

Differentiating the second equation in (<ref>) with respect to $t$, then taking $\bff{\phi}=\lambda_e \bff{H}_h$ yields

\begin{align*}
\frac{\lambda_e}{2} \ddt \norm{\bff{H}_h}{\bb{L}^2}^2
- \lambda_e
\inpro{\nabla \partial_t \bff{u}_h}{\nabla \bff{H}_h}
\kappa\mu \lambda_e \inpro{\partial_t \bff{u}_h}{\bff{H}_h}
\kappa \lambda_e \inpro{\partial_t (|\bff{u}_h|^2 \bff{u}_h)}{\bff{H}_h}
\\
\beta\lambda_e \inpro{\bff{e}(\bff{e}\cdot \partial_t \bff{u}_h)}{\bff{H}_h}.
\end{align*}

Adding this to (<ref>), then applying Hölder's inequality gives

\begin{align*}
&\frac{\lambda_e}{2} \ddt \norm{\bff{H}_h}{\bb{L}^2}^2
\norm{\partial_t \bff{u}_h}{\bb{L}^2}^2
\\
(\lambda_r + \kappa\mu \lambda_e) \inpro{\partial_t \bff{u}_h}{\bff{H}_h}
\gamma \inpro{\bff{u}_h \times \bff{H}_h}{\partial_t \bff{u}_h}
\kappa \lambda_e \inpro{\partial_t (|\bff{u}_h|^2 \bff{u}_h)}{\bff{H}_h}
\beta\lambda_e \inpro{\bff{e}(\bff{e}\cdot \partial_t \bff{u}_h)}{\bff{H}_h}
\\
B \norm{\partial_t \bff{u}_h}{\bb{L}^2}
\norm{\bff{H}_h}{\bb{L}^2}
\gamma \norm{\bff{u}_h}{\bb{L}^4}
\norm{\bff{H}_h}{\bb{L}^4}
\norm{\partial_t \bff{u}_h}{\bb{L}^2}
\kappa\lambda_e \norm{\partial_t \bff{u}_h}{\bb{L}^2}
\norm{\bff{u}_h}{\bb{L}^6}^2
\norm{\bff{H}_h}{\bb{L}^6}
\\
\frac{1}{4} \norm{\partial_t \bff{u}_h}{\bb{L}^2}^2
B^2 \norm{\bff{H}_h}{\bb{L}^2}^2
\frac{1}{4} \norm{\partial_t \bff{u}_h}{\bb{L}^2}^2
C\gamma^2 \norm{\bff{H}_h}{\bb{H}^1}^2
\frac{1}{4} \norm{\partial_t \bff{u}_h}{\bb{L}^2}^2
C (\kappa \lambda_e)^2 \norm{\bff{H}_h}{\bb{H}^1}^2,
\end{align*}

where $B:=\lambda_r+\kappa\mu\lambda_e+\beta\lambda_e$, and in the last line we used Young's inequality, Sobolev embedding $\bb{H}^1\hookrightarrow \bb{L}^6$ and (<ref>).
Rearranging the inequality and integrating with respect to $t$ (and noting (<ref>) again), we obtain

\begin{align*}
\lambda_e \norm{\bff{H}_h(t)}{\bb{L}^2}^2
\int_0^t \norm{\partial_s \bff{u}_h(s)}{\bb{L}^2}^2 \ds
\norm{\bff{H}_h(0)}{\bb{L}^2}^2
\alpha_1 \int_0^t \norm{\bff{H}_h(s)}{\bb{H}^1}^2 \ds
\\
C_1 \alpha_1 \left( \norm{\bff{u}_0}{\bb{H}^1}^2 + |\mathscr{D}| \right),
\end{align*}

where the constant $C$ depends only on $\kappa$, $\mu$, and $\mathscr{D}$, as required.

Let $h>0$ and initial data $\bff{u}_0$ be given. For all $t\in (0,\infty)$,

\begin{align}\label{equ:semidisc-est lapl}
\lambda_e \norm{\Delta_h \bff{u}_h(t)}{\bb{L}^2}^2
\leq
\end{align}

Moreover, if the triangulation $\cal{T}_h$ is globally quasi-uniform, then

\begin{align}\label{equ:semidisc-est L infty}
\lambda_e \norm{\bff{u}_h(t)}{\bb{L}^\infty}^2
\leq
\end{align}

Here, the constant $C$ depends only on $\kappa$, $\mu$, $\beta$, $\gamma$, $\lambda_r$, $\mathscr{D}$, and $K_0$ (but is independent of $t$, $h$, and $\lambda_e$).

Taking $\bff{\phi}= \lambda_e \Delta_h \bff{u}_h$, we obtain

\begin{align*}
\lambda_e \inpro{\bff{H}_h}{\Delta_h \bff{u}_h}
\lambda_e \norm{\Delta_h \bff{u}_h}{\bb{L}^2}^2
\kappa\mu \lambda_e \norm{\nabla \bff{u}_h}{\bb{L}^2}^2
\kappa \lambda_e \inpro{|\bff{u}_h|^2 \bff{u}_h}{\Delta_h \bff{u}_h}
\beta \lambda_e \inpro{\bff{e}(\bff{e}\cdot \bff{u}_h)}{\Delta_h \bff{u}_h}.
\end{align*}

Therefore, after rearranging, we have

\begin{align*}
\lambda_e \norm{\Delta_h \bff{u}_h}{\bb{L}^2}^2
\kappa\mu \lambda_e \norm{\nabla \bff{u}_h}{\bb{L}^2}^2
\lambda_e \inpro{\bff{H}_h}{\Delta_h \bff{u}_h}
\kappa \lambda_e \inpro{|\bff{u}_h|^2 \bff{u}_h}{\Delta_h \bff{u}_h}
\beta \lambda_e \inpro{\bff{e}(\bff{e}\cdot \bff{u}_h)}{\Delta_h \bff{u}_h}
\\
\kappa \mu \lambda_e \norm{\nabla \bff{u}_h}{\bb{L}^2}^2
\lambda_e \norm{\bff{H}_h}{\bb{L}^2}^2
\frac{\lambda_e}{2} \norm{\Delta_h \bff{u}_h}{\bb{L}^2}^2
\beta^2 \lambda_e  \norm{\bff{u}_h}{\bb{L}^2}^2
\\
C \lambda_e \norm{\bff{u}_h}{\bb{H}^1}^2
\lambda_e \norm{\bff{H}_h}{\bb{L}^2}^2
\frac{\lambda_e}{2}  \norm{\Delta_h \bff{u}_h}{\bb{L}^2}^2
\leq
C\lambda_e + C + \frac{\lambda_e}{2}  \norm{\Delta_h \bff{u}_h}{\bb{L}^2}^2
\end{align*}

where we used Young's inequality, (<ref>), and the Sobolev embedding $\bb{H}^1\hookrightarrow \bb{L}^6$. This proves inequality (<ref>).
Finally, (<ref>) follows from (<ref>), (<ref>), and (<ref>).

The following proposition shows stability of $\bff{H}_h$ in $L^\infty(\bb{H}^1)$ norm under some assumptions.

Let $h>0$ and initial data $\bff{u}_0$ be given. Suppose that one of the following holds:

* $d=1$ or $2$,

* $d=3$ and $\kappa^2 \lambda_e \left(\norm{\bff{u}_0}{\bb{H}^1}^4 + |\mathscr{D}|\right)\lesssim \mu$,

* $d=3$ and the triangulation is globally quasi-uniform.

For all $t\in [0,\infty)$,

\begin{align}\label{equ:semidisc-est3}
\norm{\partial_t \bff{u}_h(t)}{\bb{L}^2}^2
+ \int_0^t \norm{\partial_t \bff{H}_h(s)}{\bb{L}^2}^2 \ds
+ \int_0^t \norm{\nabla \partial_t \bff{u}_h(s)}{\bb{L}^2}^2 \ds
\leq
\end{align}

where the constant $C$ depends only on the coefficients of the equation, $\mathscr{D}$, and $K_0$.

Differentiating the first equation in (<ref>) with respect to $t$, then taking $\bff{\chi}= \partial_t \bff{u}_h$ gives

\begin{align}\label{equ:ddt dt uh}
\frac{1}{2} \ddt \norm{\partial_t \bff{u}_h}{\bb{L}^2}^2
\lambda_r \inpro{\partial_t \bff{H}_h}{\partial_t \bff{u}_h}
\lambda_e \inpro{\nabla \partial_t \bff{H}_h}{\nabla \partial_t \bff{u}_h}
\gamma \inpro{\bff{u}_h \times \partial_t \bff{H}_h}{\partial_t \bff{u}_h}.
\end{align}

Differentiating the second equation in (<ref>) with respect to $t$, then taking $\bff{\phi}= \lambda_e \partial_t \bff{H}_h$ gives

\begin{align}\label{equ:lambda e dt Hh L2}
\lambda_e \norm{\partial_t \bff{H}_h}{\bb{L}^2}^2
\inpro{\nabla \partial_t \bff{u}_h}{\nabla \partial_t \bff{H}_h}
\kappa\mu \lambda_e
\inpro{\partial_t \bff{u}_h}{\partial_t \bff{H}_h}
\kappa \lambda_e
\inpro{\partial_t (|\bff{u}_h|^2 \bff{u}_h)}{\partial_t \bff{H}_h}
\nonumber \\
\beta\lambda_e \inpro{\bff{e}(\bff{e}\cdot\partial_t \bff{u}_h}{\partial_t \bff{H}_h}.
\end{align}

while taking $\bff{\phi}= C_r \partial_t \bff{u}_h$, where $C_r:=\lambda_r+\kappa\mu\lambda_e$ gives

\begin{align}\label{equ:dt Hh dt uh L2 sub}
C_r \inpro{\partial_t \bff{H}_h}{\partial_t \bff{u}_h}
C_r \left(
\kappa \mu \norm{\partial_t \bff{u}_h}{\bb{L}^2}^2
\norm{\nabla \partial_t \bff{u}_h}{\bb{L}^2}^2
2\kappa \norm{\bff{u}_h \cdot \partial_t \bff{u}_h}{\bb{L}^2}^2
\kappa \norm{|\bff{u}_h| |\partial_t \bff{u}_h|}{\bb{L}^2}^2\right)
\nonumber \\
C_r \beta \norm{\bff{e}\cdot \partial_t \bff{u}_h}{L^2}^2.
\end{align}

Adding (<ref>), (<ref>), and (<ref>) yields

\begin{align*}
&\frac{1}{2} \ddt \norm{\partial_t \bff{u}_h}{\bb{L}^2}^2
\lambda_e \norm{\partial_t \bff{H}_h}{\bb{L}^2}^2
C_r \big( \norm{\nabla \partial_t \bff{u}_h}{\bb{L}^2}^2
2\kappa \norm{\bff{u}_h \cdot \partial_t \bff{u}_h}{\bb{L}^2}^2
\kappa \norm{|\bff{u}_h| |\partial_t \bff{u}_h|}{\bb{L}^2}^2 \big)
\\
C_r \kappa \mu \norm{\partial_t \bff{u}_h}{\bb{L}^2}^2
C_r \beta \norm{\bff{e}\cdot \partial_t \bff{u}_h}{L^2}^2
\kappa \lambda_e \inpro{\partial_t (|\bff{u}_h|^2 \bff{u}_h)}{\partial_t \bff{H}_h}
\gamma \inpro{\bff{u}_h \times \partial_t \bff{H}_h}{\partial_t \bff{u}_h},
\end{align*}

We estimate each term appearing on the right-hand side as follows.

Case 1: $d=1$ or $2$. In this case, Hölder's inequality implies

\begin{align*}
&\frac{1}{2} \ddt \norm{\partial_t \bff{u}_h}{\bb{L}^2}^2
\lambda_e \norm{\partial_t \bff{H}_h}{\bb{L}^2}^2
C_r \norm{\nabla \partial_t \bff{u}_h}{\bb{L}^2}^2
C_r\kappa \norm{|\bff{u}_h| |\partial_t \bff{u}_h|}{\bb{L}^2}^2
\\
C_r \big(\kappa\mu+\beta\big) \norm{\partial_t \bff{u}_h}{\bb{L}^2}^2
\kappa \lambda_e \norm{\bff{u}_h}{\bb{L}^8}^2
\norm{\partial_t \bff{u}_h}{\bb{L}^4}
\norm{\partial_t \bff{H}_h}{\bb{L}^2}
\gamma \norm{\bff{u}_h}{\bb{L}^4}
\norm{\partial_t \bff{H}_h}{\bb{L}^2}
\norm{\partial_t \bff{u}_h}{\bb{L}^4}
\\
C \norm{\partial_t \bff{u}_h}{\bb{L}^2}^2
\frac{C_r}{2} \norm{\nabla \partial_t \bff{u}_h}{\bb{L}^2}^2
\frac{\lambda_e}{2} \norm{\partial_t \bff{H}_h}{\bb{L}^2}^2,
\end{align*}

where in the last step we used the Sobolev embedding $\bb{H}^1 \hookrightarrow \bb{L}^8$, Young's inequality, (<ref>) and (<ref>). The required inequality then follows by integrating both sides with respect to $t$ and noting (<ref>).

Case 2: $d=3$ and $\kappa^2 \lambda_e \big(\norm{\bff{u}_0}{\bb{H}^1}^4 + |\mathscr{D}|\big)\lesssim \mu$. Similarly, by Hölder's and Young's inequalities, we have

\begin{align*}
&\frac{1}{2} \ddt \norm{\partial_t \bff{u}_h}{\bb{L}^2}^2
\lambda_e \norm{\partial_t \bff{H}_h}{\bb{L}^2}^2
C_r \norm{\nabla \partial_t \bff{u}_h}{\bb{L}^2}^2
C_r\kappa \norm{|\bff{u}_h| |\partial_t \bff{u}_h|}{\bb{L}^2}^2
\\
C_r \big(\kappa\mu+\beta\big) \norm{\partial_t \bff{u}_h}{\bb{L}^2}^2
\kappa \lambda_e \norm{\bff{u}_h}{\bb{L}^6}^2
\norm{\partial_t \bff{u}_h}{\bb{L}^6}
\norm{\partial_t \bff{H}_h}{\bb{L}^2}
\gamma \norm{\bff{u}_h}{\bb{L}^4}
\norm{\partial_t \bff{H}_h}{\bb{L}^2}
\norm{\partial_t \bff{u}_h}{\bb{L}^4}
\\
C \norm{\partial_t \bff{u}_h}{\bb{L}^2}^2
\frac{C_r}{2} \norm{\nabla \partial_t \bff{u}_h}{\bb{L}^2}^2
\frac{\lambda_e}{2} \norm{\partial_t \bff{H}_h}{\bb{L}^2}^2
C_{\mathrm{S}} \kappa^2 \lambda_e \norm{\bff{u}_h}{\bb{H}^1}^4 \norm{\partial_t \bff{u}_h}{\bb{H}^1}^2,
\end{align*}

where $C_{\mathrm{S}}$ is the constant associated with the Sobolev embedding $\bb{H}^1\hookrightarrow \bb{L}^6$.
If $\norm{\bff{u}_0}{\bb{H}^1}^4 + |\mathscr{D}|$ is sufficiently small, or more precisely (with $C_1$ as given in (<ref>))

\[
C_{\mathrm{S}} C_1^4 \kappa^2 \lambda_e \left(\norm{\bff{u}_0}{\bb{H}^1}^4 + |\mathscr{D}|\right)\leq \mu,
\]

then noting (<ref>), we can absorb the term $C_{\mathrm{S}} \kappa^2 \lambda_e \norm{\bff{u}_h}{\bb{H}^1}^4 \norm{\partial_t \bff{u}_h}{\bb{H}^1}^2$ to the left-hand side. The required inequality then follows by integrating both sides with respect to $t$ and using (<ref>).

Case 3: $d=3$ and the triangulation is quasi-uniform. In this case, by Hölder's inequality and (<ref>),

\begin{align*}
&\frac{1}{2} \ddt \norm{\partial_t \bff{u}_h}{\bb{L}^2}^2
\lambda_e \norm{\partial_t \bff{H}_h}{\bb{L}^2}^2
C_r \norm{\nabla \partial_t \bff{u}_h}{\bb{L}^2}^2
C_r\kappa \norm{|\bff{u}_h| |\partial_t \bff{u}_h|}{\bb{L}^2}^2
\\
C_r \big(\kappa\mu+\beta\big) \norm{\partial_t \bff{u}_h}{\bb{L}^2}^2
\kappa \lambda_e \norm{\bff{u}_h}{\bb{L}^\infty}^2
\norm{\partial_t \bff{u}_h}{\bb{L}^2}
\norm{\partial_t \bff{H}_h}{\bb{L}^2}
\gamma \norm{\bff{u}_h}{\bb{L}^\infty}
\norm{\partial_t \bff{H}_h}{\bb{L}^2}
\norm{\partial_t \bff{u}_h}{\bb{L}^2}
\\
C \norm{\partial_t \bff{u}_h}{\bb{L}^2}^2
\frac{C_r}{2} \norm{\nabla \partial_t \bff{u}_h}{\bb{L}^2}^2
\frac{\lambda_e}{2} \norm{\partial_t \bff{H}_h}{\bb{L}^2}^2.
\end{align*}

Integrating both sides with respect to $t$ then yields the inequality in this case.
This completes the proof of the proposition.

Let $h>0$ and initial data $\bff{u}_0$ be given. Under the same assumptions as Proposition <ref>, for all $t\in [0,\infty)$,

\begin{align*}
\norm{\nabla \bff{H}_h(t)}{\bb{L}^2}
\leq
\end{align*}

where the constant $C$ depends only on the coefficients of the equation, $\mathscr{D}$, and $K_0$.

Setting $\bff{\chi}=\bff{H}_h$ in (<ref>) gives

\begin{align*}
\lambda_e \norm{\nabla \bff{H}_h}{\bb{L}^2}^2
\inpro{\partial_t \bff{u}_h}{\bff{H}_h}
\lambda_r \norm{\bff{H}_h}{\bb{L}^2}^2
\lesssim
\norm{\partial_t \bff{u}_h}{\bb{L}^2}^2
\norm{\bff{H}_h}{\bb{L}^2}^2.
\end{align*}

The required result then follows from (<ref>) and (<ref>).

To estimate the error in the semi-discrete approximation, we write it as a sum of two terms:

\begin{align}\label{equ:theta rho}
\bff{u}_h(t) - \bff{u}(t)
(\bff{u}_h(t)- R_h \bff{u}(t)) + (R_h \bff{u}(t) - \bff{u}(t))
&=: \bff{\theta}(t) + \bff{\rho}(t),
\\
\label{equ:xi eta}
\bff{H}_h(t) - \bff{H}(t)
(\bff{H}_h(t)- R_h \bff{H}(t)) + (R_h \bff{H}(t) - \bff{H}(t))
&=: \bff{\xi}(t) + \bff{\eta}(t).
\end{align}

where $R_h$ is the Ritz projection operator (see (<ref>)).

In the analysis, to bound the nonlinear terms, we will often write (by adding and substracting $R_h \bff{u}$)

\begin{align}\label{equ:nonlinear Rh u}
\nonumber
|\bff{u}_h|^2 \bff{u}_h - |\bff{u}|^2 \bff{u}
(\bff{u}_h- R_h \bff{u}) \cdot (\bff{u}_h+ R_h \bff{u}) \bff{u}_h
|R_h \bff{u}|^2 (\bff{u}_h-\bff{u})
(R_h \bff{u} - \bff{u})\cdot (R_h \bff{u}+ \bff{u}) \bff{u}
\\
\big(\bff{\theta} \cdot (\bff{\theta}+ 2 R_h \bff{u}) \big) \bff{u}_h
\abs{R_h \bff{u}}^2 \big(\bff{\theta}+\bff{\rho}\big)
\big( \bff{\rho} \cdot (R_h \bff{u} + \bff{u}) \big) \bff{u}
\end{align}

\begin{align}\label{equ:nonlinear cross}
\bff{u}_h \times \bff{H}_h - \bff{u} \times \bff{H}
\nonumber
(\bff{u}_h - R_h \bff{u})\times \bff{H}_h
R_h \bff{u} \times (\bff{H}_h-\bff{H})
(R_h \bff{u}-\bff{u}) \times \bff{H}
\\
\bff{\theta} \times \bff{H}_h
R_h\bff{u} \times (\bff{\xi}+\bff{\eta})
\bff{\rho} \times \bff{H}.
\end{align}

We derive some estimates for the nonlinear terms in the following lemmas.

Let $\epsilon>0$ and initial data $\bff{u}_0$ be given. Let $(\bff{u},\bff{H})$ be the solution of (<ref>) satisfying (<ref>). Then we have

\begin{align}
\label{equ:uh cross Hh theta}
\big| \inpro{\bff{u}_h \times \bff{H}_h- \bff{u} \times \bff{H}}{\bff{\theta}} \big|
\norm{\bff{\theta}}{\bb{L}^2}^2
\epsilon \norm{\nabla \bff{\theta}}{\bb{L}^2}^2
\epsilon \norm{\bff{\xi}}{\bb{L}^2}^2,
\\
\label{equ:uh cross Hh xi}
\big| \inpro{\bff{u}_h \times \bff{H}_h- \bff{u} \times \bff{H}}{\bff{\xi}} \big|
\norm{\bff{\xi}}{\bb{L}^2}^2
\epsilon \norm{\nabla \bff{\xi}}{\bb{L}^2}^2
\epsilon \norm{\bff{\theta}}{\bb{L}^2}^2,
\\
\label{equ:uh cross Hh dt theta}
\big| \inpro{\bff{u}_h \times \bff{H}_h- \bff{u} \times \bff{H}}{\partial_t \bff{\theta}} \big|
\norm{\bff{\xi}}{\bb{L}^2}^2
\big( \norm{\bff{\theta}}{\bb{L}^2}^2 + \epsilon \norm{\nabla \bff{\theta}}{\bb{L}^2}^2 \big) \norm{\bff{H}_h}{\bb{H}^1}^2
\epsilon \norm{\partial_t \bff{\theta}}{\bb{L}^2}^2.
\end{align}

Furthermore, if the triangulation $\mathcal{T}_h$ is globally quasi-uniform, then for any $\bff{\zeta}\in \bb{V}_h$,

\begin{align}
\label{equ:uh cross Hh Delta}
\big| \inpro{\bff{u}_h \times \bff{H}_h- \bff{u} \times \bff{H}}{\bff{\zeta}} \big|
\norm{\bff{\theta}}{\bb{L}^2}^2
\norm{\bff{\xi}}{\bb{L}^2}^2
\epsilon \norm{\bff{\zeta}}{\bb{L}^2}^2,
\\
\label{equ:nab uh cross Hh}
\big|\inpro{\nabla\Pi_h\big(\bff{u}_h \times \bff{H}_h- \bff{u} \times \bff{H}\big)}{\bff{\zeta}} \big|
\norm{\bff{\theta}}{\bb{H}^1}^2
\norm{\bff{\xi}}{\bb{H}^1}^2
\epsilon \norm{\bff{\zeta}}{\bb{L}^2}^2.
\end{align}

Here, all constants are independent of $h$.

First, we will prove (<ref>) by writing

\begin{align}\label{equ:cross split}
\bff{u}_h \times \bff{H}_h - \bff{u} \times \bff{H}
(\bff{\theta}+\bff{\rho}) \times \bff{H}
\bff{u}_h \times (\bff{\xi}+\bff{\eta}).
\end{align}

Therefore, by Hölder's and Young's inequalities, we have

\begin{align*}
\big| \inpro{\bff{u}_h \times \bff{H}_h- \bff{u} \times \bff{H}}{\bff{\theta}} \big|
\norm{\bff{\rho}}{\bb{L}^2} \norm{\bff{H}}{\bb{L}^\infty} \norm{\bff{\theta}}{\bb{L}^2}
\norm{\bff{u}_h}{\bb{L}^4} \norm{\bff{\xi}+\bff{\eta}}{\bb{L}^2} \norm{\bff{\theta}}{\bb{L}^4}
\\
\norm{\bff{\rho}}{\bb{L}^2}^2
\norm{\bff{\theta}}{\bb{L}^2}^2
\norm{\bff{\theta}}{\bb{L}^4}^2
\epsilon \norm{\bff{\xi}+\bff{\eta}}{\bb{L}^2}^2
\\
\norm{\bff{\theta}}{\bb{L}^2}^2
\epsilon \norm{\nabla \bff{\theta}}{\bb{L}^2}^2
\epsilon \norm{\bff{\xi}}{\bb{L}^2}^2,
\end{align*}

where in the last step we used (<ref>) and (<ref>), thus proving (<ref>).

\begin{align*}
\big| \inpro{\bff{u}_h \times \bff{H}_h- \bff{u} \times \bff{H}}{\bff{\xi}} \big|
\norm{\bff{\theta}+\bff{\rho}}{\bb{L}^2} \norm{\bff{H}}{\bb{L}^\infty} \norm{\bff{\xi}}{\bb{L}^2}
\norm{\bff{u}_h}{\bb{L}^4} \norm{\bff{\eta}}{\bb{L}^2} \norm{\bff{\xi}}{\bb{L}^4}
\\
\epsilon \norm{\bff{\theta}}{\bb{L}^2}^2
\norm{\bff{\rho}}{\bb{L}^2}^2
\norm{\bff{\xi}}{\bb{L}^2}^2
\norm{\bff{\xi}}{\bb{L}^4}^2
\epsilon \norm{\bff{\eta}}{\bb{L}^2}^2
\\
\norm{\bff{\xi}}{\bb{L}^2}^2
\epsilon \norm{\nabla \bff{\xi}}{\bb{L}^2}^2
\epsilon \norm{\bff{\theta}}{\bb{L}^2}^2,
\end{align*}

proving (<ref>).
Next, using (<ref>), Hölder's and Young's inequalities, we obtain

\begin{align*}
&\big| \inpro{\bff{u}_h \times \bff{H}_h- \bff{u} \times \bff{H}}{\partial_t \bff{\theta}} \big|
\\
\norm{\bff{\theta}}{\bb{L}^4}
\norm{\bff{H}_h}{\bb{L}^4}
\norm{\partial_t \bff{\theta}}{\bb{L}^2}
\norm{R_h \bff{u}}{\bb{L}^\infty}
\norm{\bff{\xi}+\bff{\eta}}{\bb{L}^2}
\norm{\partial_t \bff{\theta}}{\bb{L}^2}
\norm{\bff{\rho}}{\bb{L}^2}
\norm{\bff{H}}{\bb{L}^\infty}
\norm{\partial_t \bff{\theta}}{\bb{L}^2}
\\
\norm{\bff{\theta}}{\bb{L}^4}^2 \norm{\bff{H}_h}{\bb{L}^4}^2
\epsilon \norm{\partial_t \bff{\theta}}{\bb{L}^2}^2
\norm{\bff{\xi}}{\bb{L}^2}^2
\epsilon \norm{\partial_t \bff{\theta}}{\bb{L}^2}^2
\\
\big(\norm{\bff{\theta}}{\bb{L}^2}^2 + \epsilon \norm{\nabla \bff{\theta}}{\bb{L}^2}^2 \big)
\norm{\bff{H}_h}{\bb{H}^1}^2
\norm{\bff{\xi}}{\bb{L}^2}^2
\epsilon \norm{\partial_t \bff{\theta}}{\bb{L}^2}^2,
\end{align*}

where in the last step we used Sobolev embedding and (<ref>). This proves (<ref>).

Finally, if the triangulation is quasi-uniform, then (<ref>) holds. Using (<ref>), Hölder's and Young's inequalities, we obtain

\begin{align*}
\big| \inpro{\bff{u}_h \times \bff{H}_h- \bff{u} \times \bff{H}}{\bff{\zeta}} \big|
\left( \norm{\bff{\theta}+\bff{\rho}}{\bb{L}^2} \norm{\bff{H}}{\bb{L}^\infty}
+ \norm{\bff{u}_h}{\bb{L}^\infty} \norm{\bff{\xi}+\bff{\eta}}{\bb{L}^2} \right) \norm{\bff{\zeta}}{\bb{L}^2}
\\
\norm{\bff{\theta}}{\bb{L}^2}^2
\norm{\bff{\xi}}{\bb{L}^2}^2
\epsilon \norm{\bff{\zeta}}{\bb{L}^2}^2,
\end{align*}

thus proving (<ref>). Similarly, noting (<ref>), we obtain

\begin{align*}
\big|\inpro{\nabla\Pi_h\big(\bff{u}_h \times \bff{H}_h- \bff{u} \times \bff{H}\big)}{\bff{\zeta}} \big|
\big( \norm{\nabla \bff{\theta}+\nabla \bff{\rho}}{\bb{L}^2} \norm{\bff{H}}{\bb{L}^\infty}
+ \norm{\bff{\theta}+\bff{\rho}}{\bb{L}^4} \norm{\nabla \bff{H}}{\bb{L}^4}
\\
+ \norm{\nabla \bff{u}_h}{\bb{L}^4} \norm{\bff{\xi}+\bff{\eta}}{\bb{L}^4}
+ \norm{\bff{u}_h}{\bb{L}^\infty} \norm{\nabla \bff{\xi}+\nabla \bff{\eta}}{\bb{L}^2} \big) \norm{\bff{\zeta}}{\bb{L}^2}
\\
\norm{\bff{\theta}}{\bb{H}^1}^2
\norm{\bff{\xi}}{\bb{H}^1}^2
\epsilon \norm{\bff{\zeta}}{\bb{L}^2}^2,
\end{align*}

where in the last step we used (<ref>), Proposition <ref>, and (<ref>), completing the proof of the lemma.

Let $\epsilon>0$ and initial data $\bff{u}_0$ be given. Let $(\bff{u},\bff{H})$ be the solution of (<ref>) satisfying (<ref>). Then we have

\begin{align}
\label{equ:uh2uh theta}
\big| \inpro{|\bff{u}_h|^2 \bff{u}_h - |\bff{u}|^2 \bff{u}}{\bff{\theta}} \big|
\norm{\bff{\theta}}{\bb{L}^2}^2
\epsilon \norm{\nabla \bff{\theta}}{\bb{L}^2}^2,
\\
\label{equ:uh2uh xi}
\big| \inpro{|\bff{u}_h|^2 \bff{u}_h - |\bff{u}|^2 \bff{u}}{\bff{\xi}} \big|
\big( 1+\norm{\bff{\xi}}{\bb{H}^1}^2 \big) \norm{\bff{\theta}}{\bb{L}^2}^2
\epsilon \norm{\nabla \bff{\theta}}{\bb{L}^2}^2
\epsilon \norm{\bff{\xi}}{\bb{L}^2}^2,
\\
\label{equ:uh2uh dt theta}
\big| \inpro{|\bff{u}_h|^2 \bff{u}_h - |\bff{u}|^2 \bff{u}}{\partial_t \bff{\theta}} \big|
\big( 1+\norm{\partial_t \bff{\theta}}{\bb{L}^2}^2 \big) \norm{\bff{\theta}}{\bb{H}^1}^2
\epsilon \norm{\partial_t \bff{\theta}}{\bb{L}^2}^2,
\\
\label{equ:dt uh2uh xi}
\big| \inpro{\partial_t(|\bff{u}_h|^2 \bff{u}_h - |\bff{u}|^2 \bff{u})}{\bff{\xi}} \big|
\big( 1 + \norm{\bff{\xi}}{\bb{H}^1}^2 + \norm{\partial_t \bff{u}_h}{\bb{L}^2}^2 \big) \norm{\bff{\theta}}{\bb{H}^1}^2
\norm{\bff{\xi}}{\bb{L}^2}^2
\epsilon \norm{\partial_t \bff{\theta}}{\bb{L}^2}^2.
\end{align}

Furthermore, if the triangulation $\mathcal{T}_h$ is globally quasi-uniform, then for any $\bff{\zeta}\in \bb{V}_h$,

\begin{align}
\label{equ:uh2uh Delta theta}
\big| \inpro{|\bff{u}_h|^2 \bff{u}_h - |\bff{u}|^2 \bff{u}}{\bff{\zeta}} \big|
\norm{\bff{\theta}}{\bb{L}^2}^2
\epsilon \norm{\bff{\zeta}}{\bb{L}^2}^2,
\\
\label{equ:dt uh2uh}
\big| \inpro{\partial_t \left(|\bff{u}_h|^2 \bff{u}_h - |\bff{u}|^2 \bff{u}\right)}{\bff{\zeta}} \big|
h^{2(r+1)} + \norm{\bff{\theta}}{\bb{L}^2}^2
+ \norm{\partial_t \bff{\theta}}{\bb{L}^2}^2
\epsilon \norm{\bff{\zeta}}{\bb{H}^1}^2.
\end{align}

In the above estimates, all constants are independent of $h$.

First, we will prove (<ref>). Note that

\begin{align}\label{equ:simple uhu}
|\bff{u}_h|^2 \bff{u}_h - |\bff{u}|^2 \bff{u}
|\bff{u}_h|^2 (\bff{\theta}+\bff{\rho}) - \big( (\bff{\theta}+\bff{\rho})\cdot (\bff{u}_h+\bff{u}) \big) \bff{u}.
\end{align}

Therefore, noting (<ref>) and (<ref>), by Hölder's inequality and Sobolev embedding $\bb{H}^1 \hookrightarrow \bb{L}^6$ we have

\begin{align*}
\big| \inpro{|\bff{u}_h|^2 \bff{u}_h - |\bff{u}|^2 \bff{u}}{\bff{\theta}} \big|
\norm{\bff{u}_h}{\bb{L}^4}^2 \norm{\bff{\theta}}{\bb{L}^4}^2
\norm{\bff{u}_h}{\bb{L}^6}^2 \norm{\bff{\rho}}{\bb{L}^2} \norm{\bff{\theta}}{\bb{L}^6}
\norm{\bff{u}}{\bb{L}^\infty}^2
\big(
\norm{\bff{\theta}}{\bb{L}^2}^2
\norm{\bff{\theta}}{\bb{L}^2}
\norm{\bff{\rho}}{\bb{L}^2}
\big)
\\
\norm{\bff{u}_h}{\bb{L}^4}
\norm{\bff{u}}{\bb{L}^\infty}
\big(
\norm{\bff{\theta}}{\bb{L}^4}
\norm{\bff{\theta}}{\bb{L}^2}
\norm{\bff{\theta}}{\bb{L}^4}
\norm{\bff{\rho}}{\bb{L}^2}
\big)
\\
\norm{\bff{\theta}}{\bb{L}^2}^2
\epsilon \norm{\nabla \bff{\theta}}{\bb{L}^2}^2
\end{align*}

for any $\epsilon>0$, where we used Young's inequality and (<ref>) in the last step. Next, we will prove (<ref>). In this case, using (<ref>), Hölder's and Young's inequality, and (<ref>), we have

\begin{align*}
\nonumber
\big| \inpro{|\bff{u}_h|^2 \bff{u}_h - |\bff{u}|^2 \bff{u}}{\bff{\xi}} \big|
\norm{\bff{\theta}}{\bb{L}^2}
\norm{\bff{\theta}}{\bb{L}^6}
\norm{\bff{u}_h}{\bb{L}^6}
\norm{\bff{\xi}}{\bb{L}^6}
\norm{\bff{\theta}}{\bb{L}^4}
\norm{R_h \bff{u}}{\bb{L}^\infty}
\norm{\bff{u}_h}{\bb{L}^4}
\norm{\bff{\xi}}{\bb{L}^2}
\\
\nonumber
\norm{R_h \bff{u}}{\bb{L}^\infty}^2
\norm{\bff{\theta}+\bff{\rho}}{\bb{L}^2}
\norm{\bff{\xi}}{\bb{L}^2}
\norm{\bff{\rho}}{\bb{L}^2}
\norm{R_h \bff{u}+ \bff{u}}{\bb{L}^\infty}
\norm{\bff{u}}{\bb{L}^\infty}
\norm{\bff{\xi}}{\bb{L}^2}
\\
\nonumber
\norm{\bff{\xi}}{\bb{L}^6}^2 \norm{\bff{\theta}}{\bb{L}^2}^2
\epsilon \norm{\bff{\theta}}{\bb{L}^6}^2
\norm{\bff{\theta}}{\bb{L}^2}^2
\epsilon \norm{\nabla \bff{\theta}}{\bb{L}^2}^2
\epsilon \norm{\bff{\xi}}{\bb{L}^2}^2
\norm{\bff{\rho}}{\bb{L}^2}^2
\epsilon \norm{\bff{\xi}}{\bb{L}^2}^2
\\
\big( 1+\norm{\bff{\xi}}{\bb{H}^1}^2 \big) \norm{\bff{\theta}}{\bb{L}^2}^2
\epsilon \norm{\nabla \bff{\theta}}{\bb{L}^2}^2
\epsilon \norm{\bff{\xi}}{\bb{L}^2}^2
\end{align*}

for any $\epsilon>0$, where in the last step we used Sobolev embedding $\bb{H}^1 \hookrightarrow \bb{L}^6$ and (<ref>).
Similarly, to prove (<ref>), we use (<ref>), Hölder's and Young's inequality to obtain

\begin{align*}
\big| \inpro{|\bff{u}_h|^2 \bff{u}_h - |\bff{u}|^2 \bff{u}}{\partial_t \bff{\theta}} \big|
\norm{\bff{\theta}}{\bb{L}^6}^2
\norm{\bff{u}_h}{\bb{L}^6}
\norm{\partial_t \bff{\theta}}{\bb{L}^2}
\norm{\bff{\theta}}{\bb{L}^6}
\norm{R_h \bff{u}}{\bb{L}^6}
\norm{\bff{u}_h}{\bb{L}^6}
\norm{\partial_t \bff{\theta}}{\bb{L}^2}
\\
\norm{R_h \bff{u}}{\bb{L}^\infty}^2
\norm{\bff{\theta}+\bff{\rho}}{\bb{L}^2}
\norm{\partial_t \bff{\theta}}{\bb{L}^2}
\norm{\bff{\rho}}{\bb{L}^2}
\norm{R_h \bff{u}+ \bff{u}}{\bb{L}^\infty}
\norm{\bff{u}}{\bb{L}^\infty}
\norm{\partial_t \bff{\theta}}{\bb{L}^2}
\\
\norm{\bff{\theta}}{\bb{H}^1}^2 \norm{\partial_t \bff{\theta}}{\bb{L}^2}^2
\norm{\bff{\theta}}{\bb{H}^1}^2
\epsilon \norm{\partial_t \bff{\theta}}{\bb{L}^2}
\norm{\bff{\rho}}{\bb{L}^2}^2
\epsilon \norm{\partial_t \bff{\theta}}{\bb{L}^2}
\\
\big( 1+\norm{\partial_t \bff{\theta}}{\bb{L}^2}^2 \big) \norm{\bff{\theta}}{\bb{H}^1}^2
\epsilon \norm{\partial_t \bff{\theta}}{\bb{L}^2}^2
\end{align*}

for any $\epsilon>0$, as required. Finally, differentiating (<ref>) with respect to $t$ yields

\begin{align}\label{equ:dt uh2uh T1 to 8}
\nonumber
&\big| \inpro{\partial_t(|\bff{u}_h|^2 \bff{u}_h - |\bff{u}|^2 \bff{u})}{\bff{\xi}} \big|
\\
\nonumber
\big| \inpro{\big( \partial_t\bff{\theta} \cdot (\bff{\theta}+ 2R_h \bff{u}) \big) \bff{u}_h}{\bff{\xi}} \big|
\big| \inpro{\big( \bff{\theta}\cdot (\partial_t \bff{\theta} + 2 R_h \partial_t \bff{u}) \big) \bff{u}_h}{\bff{\xi}} \big|
\big| \inpro{\big( \bff{\theta} \cdot (\bff{\theta} + R_h \bff{u}) \big) \partial_t \bff{u}_h}{\bff{\xi}} \big|
\\
\nonumber
\big| \inpro{2(R_h \bff{u} \cdot R_h \partial_t \bff{u}) (\bff{\theta}+\bff{\rho})}{\bff{\xi}} \big|
\big| \inpro{|R_h \bff{u}|^2 (\partial_t \bff{\theta}+\partial_t \bff{\rho})}{\bff{\xi}} \big|
\\
\nonumber
\big| \inpro{\big( \partial_t \bff{\rho}\cdot (R_h \bff{u}+\bff{u}) \big) \bff{u}}{\bff{\xi}} \big|
\big| \inpro{\big(\bff{\rho} \cdot (R_h \partial_t \bff{u} + \partial_t \bff{u}) \big) \bff{u}}{\bff{\xi}} \big|
\big| \inpro{\big(\bff{\rho} \cdot (R_h \bff{u} + \bff{u}) \big) \partial_t \bff{u}}{\bff{\xi}} \big|
\\
|T_1| + |T_2| + \ldots + |T_8|.
\end{align}

We will now estimate each term in the last step. For the term $|T_1|$, noting $\bff{u}_h= \bff{\theta}+R_h \bff{u}$ and applying Hölder's inequality, we have

\begin{align*}
\big| \inpro{( \partial_t\bff{\theta} \cdot \bff{\theta} ) \bff{u}_h}{\bff{\xi}} \big|
\big| \inpro{(\partial_t \bff{\theta} \cdot R_h \bff{u}) (\bff{\theta}+R_h \bff{u})}{\bff{\xi}} \big|
\\
\norm{\partial_t \bff{\theta}}{\bb{L}^2} \norm{\bff{\theta}}{\bb{L}^6} \norm{\bff{u}_h}{\bb{L}^6} \norm{\bff{\xi}}{\bb{L}^6}
\norm{\partial_t \bff{\theta}}{\bb{L}^2} \norm{R_h \bff{u}}{\bb{L}^\infty} \norm{\bff{\theta}}{\bb{L}^4} \norm{\bff{\xi}}{\bb{L}^4}
\norm{\partial_t \bff{\theta}}{\bb{L}^2} \norm{R_h \bff{u}}{\bb{L}^\infty}^2 \norm{\bff{\xi}}{\bb{L}^2}
\\
\norm{\bff{\xi}}{\bb{H}^1}^2 \norm{\bff{\theta}}{\bb{H}^1}^2
\norm{\bff{\xi}}{\bb{L}^2}^2
\epsilon \norm{\partial_t \bff{\theta}}{\bb{L}^2}^2
\end{align*}

for any $\epsilon>0$, where we used Young's inequality and Sobolev embedding in the last step. For the term $|T_2|$, by Hölder's and Young's inequality,

\begin{align*}
\norm{\bff{\theta}}{\bb{L}^6} \norm{\partial_t \bff{\theta}}{\bb{L}^2} \norm{\bff{u}_h}{\bb{L}^6} \norm{\bff{\xi}}{\bb{L}^6}
\norm{\bff{\theta}}{\bb{L}^4} \norm{R_h \partial_t \bff{u}}{\bb{L}^\infty} \norm{\bff{u}_h}{\bb{L}^4} \norm{\bff{\xi}}{\bb{L}^2}
\\
\norm{\bff{\xi}}{\bb{H}^1}^2 \norm{\bff{\theta}}{\bb{H}^1}^2
\epsilon \norm{\partial_t \bff{\theta}}{\bb{L}^2}^2
\norm{\bff{\theta}}{\bb{H}^1}^2
\epsilon \norm{\bff{\xi}}{\bb{L}^2}^2,
\end{align*}

where we used Sobolev embedding in the last step. For the term $|T_3|$, noting $\bff{u}_h= \bff{\theta}+R_h \bff{u}$ as done before, we obtain

\begin{align*}
\big| \inpro{( \bff{\theta} \cdot \bff{\theta}) \partial_t \bff{u}_h}{\bff{\xi}} \big|
\big| \inpro{( \bff{\theta} \cdot  R_h \bff{u} ) (\partial_t \bff{\theta}+ R_h \partial_t \bff{u})}{\bff{\xi}} \big|
\\
\norm{\bff{\theta}}{\bb{L}^6}^2 \norm{\partial_t \bff{u}_h}{\bb{L}^2} \norm{\bff{\xi}}{\bb{L}^6}
\norm{\bff{\theta}}{\bb{L}^4} \norm{R_h \bff{u}}{\bb{L}^\infty} \norm{\partial_t \bff{\theta}}{\bb{L}^2} \norm{\bff{\xi}}{\bb{L}^4}
\norm{\bff{\theta}}{\bb{L}^2} \norm{R_h \bff{u}}{\bb{L}^\infty} \norm{R_h \partial_t \bff{u}}{\bb{L}^\infty} \norm{\bff{\xi}}{\bb{L}^2}
\\
\norm{\partial_t \bff{u}_h}{\bb{L}^2}^2 \norm{\bff{\theta}}{\bb{H}^1}^2
\norm{\bff{\xi}}{\bb{H}^1}^2 \norm{\bff{\theta}}{\bb{H}^1}^2
\epsilon \norm{\partial_t \bff{\theta}}{\bb{L}^2}^2
\norm{\bff{\theta}}{\bb{L}^2}^2
\epsilon \norm{\bff{\xi}}{\bb{L}^2}^2.
\end{align*}

For the term $|T_4|$, by Hölder's and Young's inequalities, we have

\begin{align*}
\norm{R_h\bff{u}}{\bb{L}^\infty} \norm{R_h \partial_t \bff{u}}{\bb{L}^\infty} \norm{\bff{\theta}+\bff{\rho}}{\bb{L}^2} \norm{\bff{\xi}}{\bb{L}^2}
\lesssim
\norm{\bff{\theta}}{\bb{L}^2}^2
\norm{\bff{\rho}}{\bb{L}^2}^2
\epsilon \norm{\bff{\xi}}{\bb{L}^2}^2.
\end{align*}

Similarly for the next term,

\begin{align*}
\norm{R_h \bff{u}}{\bb{L}^\infty}^2 \norm{\partial_t \bff{\theta}+\partial_t \bff{\rho}}{\bb{L}^2} \norm{\bff{\xi}}{\bb{L}^2}
\lesssim
\norm{\partial_t \bff{\rho}}{\bb{L}^2}^2
\norm{\bff{\xi}}{\bb{L}^2}^2
\epsilon \norm{\partial_t \bff{\theta}}{\bb{L}^2}^2
\end{align*}

The terms $|T_6|, |T_7|$ and $|T_8|$ can be bounded in a similar way leading to

\begin{align*}
\leq
\norm{\bff{\rho}}{\bb{L}^2}^2
\norm{\partial_t \bff{\rho}}{\bb{L}^2}^2
\epsilon \norm{\bff{\xi}}{\bb{L}^2}^2.
\end{align*}

Altogether, noting (<ref>), we conclude the estimate (<ref>) from (<ref>).

Now, suppose the triangulation is globally quasi-uniform. In this case, (<ref>) holds. Using (<ref>), (<ref>), and Hölder's and Young's inequality, we obtain

\begin{align*}
\big| \inpro{|\bff{u}_h|^2 \bff{u}_h - |\bff{u}|^2 \bff{u}}{\Delta_h \bff{\theta}} \big|
\left( \norm{\bff{u}_h}{\bb{L}^\infty}^2 \norm{\bff{\theta}+\bff{\rho}}{\bb{L}^2}
\norm{\bff{\theta}+\bff{\rho}}{\bb{L}^2} \norm{\bff{u}_h+\bff{u}}{\bb{L}^\infty} \norm{\bff{u}}{\bb{L}^\infty} \right) \norm{\bff{\zeta}}{\bb{L}^2}
\\
h^{2(r+1)} + \norm{\bff{\theta}}{\bb{L}^2}^2 + \epsilon \norm{\bff{\zeta}}{\bb{L}^2}^2,
\end{align*}

proving (<ref>). Finally, writing

\begin{align*}
\partial_t \left(|\bff{u}_h|^2 \bff{u}_h - |\bff{u}|^2 \bff{u}\right)
\abs{\bff{u}_h}^2 (\partial_t \bff{u}_h-\partial_t \bff{u})
(\abs{\bff{u}_h}^2-\abs{\bff{u}}^2) \partial_t \bff{u}
2\left[\bff{u}_h \cdot (\partial_t \bff{u}_h-\partial_t \bff{u})\right] \bff{u}_h
\\
2 \left[(\bff{u}_h-\bff{u})\cdot \partial_t \bff{u}\right] \bff{u}_h
2 \left(\bff{u}\cdot \partial_t \bff{u}\right) (\bff{u}_h-\bff{u}),
\end{align*}

then applying Hölder's and Young's inequality, we have

\begin{align*}
&\big| \inpro{\partial_t \left(|\bff{u}_h|^2 \bff{u}_h - |\bff{u}|^2 \bff{u}\right)}{\bff{\zeta}} \big|
\\
\norm{\bff{u}_h}{\bb{L}^6}^2 \norm{\partial_t \bff{\theta}+\partial_t \bff{\rho}}{\bb{L}^2} \norm{\bff{\zeta}}{\bb{L}^6}
\norm{\bff{\theta}+\bff{\rho}}{\bb{L}^2} \norm{\bff{u}_h+\bff{u}}{\bb{L}^6} \norm{\partial_t \bff{u}}{\bb{L}^6} \norm{\bff{\zeta}}{\bb{L}^6}
\\
2 \norm{\bff{u}_h}{\bb{L}^6} \norm{\partial_t \bff{\theta}+\partial_t \bff{\rho}}{\bb{L}^2} \norm{\bff{u}_h}{\bb{L}^6} \norm{\bff{\zeta}}{\bb{L}^6}
2 \norm{\bff{\theta}+\bff{\rho}}{\bb{L}^2} \norm{\partial_t \bff{u}}{\bb{L}^6} \norm{\bff{u}_h}{\bb{L}^6} \norm{\bff{\zeta}}{\bb{L}^6}
\\
2 \norm{\bff{u}}{\bb{L}^6} \norm{\partial_t \bff{u}}{\bb{L}^6} \norm{\bff{\theta}+\bff{\rho}}{\bb{L}^2} \norm{\bff{\zeta}}{\bb{L}^6}
\\
h^{2(r+1)} + \norm{\bff{\theta}}{\bb{L}^2}^2
+ \norm{\partial_t \bff{\theta}}{\bb{L}^2}^2
\epsilon \norm{\bff{\zeta}}{\bb{H}^1}^2,
\end{align*}

thus proving (<ref>). This completes the proof of the lemma.

We will prove that the semi-discrete scheme converges at an optimal rate. For simplicity of presentation, we will assume that $\bff{u}_{0,h}=R_h\bff{u}(0)$ is the approximation in $\bb{V}_h$ of the initial data, so that $\bff{\theta}(0)=0$. In this case, subtracting the second equation in (<ref>) from that in (<ref>) yields

\begin{align*}
\inpro{\bff{\xi}+\bff{\eta}}{\bff{\phi}}
-\inpro{\nabla \bff{\theta}}{\nabla \bff{\phi}}
\kappa\mu \inpro{\bff{\theta}+\bff{\rho}}{\bff{\phi}}
\kappa \inpro{|\bff{u}_h|^2 \bff{u}_h - |\bff{u}|^2 \bff{u}}{\bff{\phi}}.
\end{align*}

Therefore, at $t=0$, noting (<ref>) and taking $\bff{\phi}=\bff{\xi}(0)$, we have

\begin{align*}
\norm{\bff{\xi}(0)}{\bb{L}^2}^2
\inpro{\bff{\eta}(0)}{\bff{\xi}(0)}
\kappa \mu \inpro{\bff{\rho}(0)}{\bff{\xi}(0)}
\kappa \inpro{|R_h\bff{u}(0)|^2 \bff{\rho}(0) - \big(\bff{\rho}(0)\cdot (R_h\bff{u}(0)+\bff{u}(0)) \big) \bff{u}(0)}{\bff{\xi}(0)}
\\
\norm{\bff{\eta}(0)}{\bb{L}^2}^2 + \norm{\bff{\rho}(0)}{\bb{L}^2}^2
\epsilon \norm{\bff{\xi}(0)}{\bb{L}^2}^2
\end{align*}

for any $\epsilon>0$ by Young's inequality. Rearranging the above, then applying (<ref>) and (<ref>), we obtain

\begin{align}\label{equ:assum xi}
\norm{\bff{\xi}(0)}{\bb{L}^2}
\lesssim
\end{align}

Next, we prove some bounds for $\bff{\theta}(t)$ and $\bff{\xi}(t)$ for $t\in (0,T)$. In particular, (<ref>) shows a superconvergence estimate for $\bff{\theta}$, which implies the corresponding estimate (<ref>) in the $\bb{L}^\infty$ norm under an additional assumption of global quasi-uniformity of the triangulation.

For any $t\in (0,T)$,

\begin{align}\label{equ:semidisc theta xi}
\norm{\bff{\theta}(t)}{\bb{L}^2}^2
\norm{\nabla \bff{\theta}(t)}{\bb{L}^2}^2
\norm{\bff{\xi}(t)}{\bb{L}^2}^2
\int_0^t \norm{\nabla \bff{\xi}(s)}{\bb{L}^2}^2 \ds
\int_0^t \norm{\partial_t \bff{\theta}(s)}{\bb{L}^2}^2 \ds
&\leq Ch^{2(r+1)}.
\end{align}

Moreover, if the triangulation $\mathcal{T}_h$ is globally quasi-uniform then

\begin{align}\label{equ:semidisc theta xi infty}
\norm{\bff{\theta}(t)}{\bb{L}^\infty}^2
\norm{\Delta_h \bff{\theta}(t)}{\bb{L}^2}^2
\int_0^t \norm{\bff{\xi}(s)}{\bb{L}^\infty}^2 \ds
\int_0^t \norm{\Delta_h \bff{\xi}(s)}{\bb{L}^2}^2 \ds
\\
\label{equ:semidisc xi L2}
\norm{\nabla \bff{\xi}(t)}{\bb{L}^2}^2
\int_0^t \norm{\nabla \Delta_h \bff{\xi}(s)}{\bb{L}^2}^2 \ds
\end{align}

The constant $C$ depends on the coefficients of the equation, $|\mathscr{D}|$, $T$, and $K_0$ (as defined in (<ref>)), but is independent of $h$.

Subtracting (<ref>) from (<ref>), using (<ref>), (<ref>) (and noting the definition of Ritz projection), we obtain for all $\bff{\chi}$, $\bff{\phi}\in \bb{V}_h$,

\begin{align}\label{equ:dt theta + dt rho}
\inpro{\partial_t \bff{\theta}+ \partial_t \bff{\rho}}{\bff{\chi}}
\lambda_r \inpro{\bff{\xi}+\bff{\eta}}{\bff{\chi}}
\lambda_e \inpro{\nabla \bff{\xi}}{\nabla \bff{\chi}}
\gamma \inpro{\bff{u}_h \times \bff{H}_h - \bff{u} \times \bff{H}}{\bff{\chi}}
\end{align}

\begin{align}
\label{equ:xi + eta}
\inpro{\bff{\xi}+\bff{\eta}}{\bff{\phi}}
-\inpro{\nabla \bff{\theta}}{\nabla \bff{\phi}}
\kappa\mu \inpro{\bff{\theta}+\bff{\rho}}{\bff{\phi}}
\kappa \inpro{|\bff{u}_h|^2 \bff{u}_h - |\bff{u}|^2 \bff{u}}{\bff{\phi}}
\beta \inpro{\bff{e}(\bff{e}\cdot \bff{\theta})}{\bff{\phi}}
\end{align}

Taking $\bff{\chi}= \bff{\theta}$ in (<ref>), we have

\begin{align}\label{equ:dt theta chi theta}
\frac{1}{2} \ddt \norm{\bff{\theta}}{\bb{L}^2}^2
\inpro{\partial_t \bff{\rho}}{\bff{\theta}}
\lambda_r \inpro{\bff{\xi}+\bff{\eta}}{\bff{\theta}}
\lambda_e \inpro{\nabla \bff{\xi}}{\nabla \bff{\theta}}
\gamma \inpro{\bff{u}_h \times \bff{H}_h - \bff{u} \times \bff{H}}{\bff{\theta}}.
\end{align}

Taking $\bff{\phi}= \lambda_r \bff{\theta}$ in (<ref>), we obtain

\begin{align}\label{equ:xi+eta theta}
\lambda_r \inpro{\bff{\xi}+\bff{\eta}}{\bff{\theta}}
&-\lambda_r \norm{\nabla \bff{\theta}}{\bb{L}^2}^2
\kappa\mu \lambda_r \norm{\bff{\theta}}{\bb{L}^2}^2
\beta \lambda_r \norm{\bff{e}\cdot\bff{\theta}}{L^2}^2
\kappa \mu\lambda_r \inpro{\bff{\rho}}{\bff{\theta}}
\nonumber \\
\kappa \lambda_r \inpro{|\bff{u}_h|^2 \bff{u}_h - |\bff{u}|^2 \bff{u}}{\bff{\theta}}.
\end{align}

Next, taking $\bff{\phi}= \lambda_e \bff{\xi}$, we have

\begin{align}\label{equ:xi phi=xi}
\lambda_e \norm{\bff{\xi}}{\bb{L}^2}^2
\lambda_e \inpro{\bff{\eta}}{\bff{\xi}}
\lambda_e \inpro{\nabla \bff{\theta}}{\nabla \bff{\xi}}
\kappa\mu \lambda_e \inpro{\bff{\theta}+\bff{\rho}}{\bff{\xi}}
\kappa \lambda_e \inpro{|\bff{u}_h|^2 \bff{u}_h - |\bff{u}|^2 \bff{u}}{\bff{\xi}}
\nonumber \\
\beta\lambda_e \inpro{\bff{e}(\bff{e}\cdot\bff{\theta})}{\bff{\xi}}.
\end{align}

Substituting (<ref>) into (<ref>), then adding the result to (<ref>) gives

\begin{align*}
&\frac{1}{2} \ddt \norm{\bff{\theta}}{\bb{L}^2}^2
\lambda_r \norm{\nabla \bff{\theta}}{\bb{L}^2}^2
\lambda_e \norm{\bff{\xi}}{\bb{L}^2}^2
\beta\lambda_r \norm{\bff{e}\cdot\bff{\theta}}{L^2}^2
\\
\inpro{\partial_t \bff{\rho}}{\bff{\theta}}
\lambda_e \inpro{\bff{\eta}}{\bff{\xi}}
\kappa\mu \lambda_e \inpro{\bff{\theta}+\bff{\rho}}{\bff{\xi}}
\kappa \mu\lambda_r \norm{\bff{\theta}}{\bb{L}^2}^2
\kappa \mu\lambda_r \inpro{\bff{\rho}}{\bff{\theta}}
\beta\lambda_e \inpro{\bff{e}(\bff{e}\cdot\bff{\theta})}{\bff{\xi}}
\\
\gamma \inpro{\bff{u}_h \times \bff{H}_h - \bff{u} \times \bff{H}}{\bff{\theta}}
\kappa \lambda_e \inpro{|\bff{u}_h|^2 \bff{u}_h - |\bff{u}|^2 \bff{u}}{\bff{\xi}}
\kappa \lambda_r \inpro{|\bff{u}_h|^2 \bff{u}_h - |\bff{u}|^2 \bff{u}}{\bff{\theta}}.
\end{align*}

Applying Hölder's and Young's inequalities, (<ref>), (<ref>), and (<ref>) to the above equation yields

\begin{align}\label{equ:theta nab theta xi}
\nonumber
&\frac{1}{2} \ddt \norm{\bff{\theta}}{\bb{L}^2}^2
\lambda_r \norm{\nabla \bff{\theta}}{\bb{L}^2}^2
\lambda_e \norm{\bff{\xi}}{\bb{L}^2}^2
\\
\nonumber
\norm{\partial_t \bff{\rho}}{\bb{L}^2}
\norm{\bff{\theta}}{\bb{L}^2}
\lambda_e \norm{\bff{\eta}}{\bb{L}^2}
\norm{\bff{\xi}}{\bb{L}^2}
\kappa\mu\lambda_e \left( \norm{\bff{\theta}}{\bb{L}^2}
\norm{\bff{\xi}}{\bb{L}^2}
\norm{\bff{\rho}}{\bb{L}^2}
\norm{\bff{\xi}}{\bb{L}^2} \right)
\kappa\mu\lambda_r \left( \norm{\bff{\theta}}{\bb{L}^2}^2
\norm{\bff{\rho}}{\bb{L}^2}
\norm{\bff{\theta}}{\bb{L}^2} \right)
\\
\nonumber
\beta\lambda_e \norm{\bff{\theta}}{\bb{L}^2} \norm{\bff{\xi}}{\bb{L}^2}
\gamma \big| \inpro{\bff{u}_h \times \bff{H}_h - \bff{u} \times \bff{H}}{\bff{\theta}} \big|
\kappa\lambda_e \big| \inpro{|\bff{u}_h|^2 \bff{u}_h - |\bff{u}|^2 \bff{u}}{\bff{\xi}} \big|
\\
\nonumber
\kappa\lambda_r \big| \inpro{|\bff{u}_h|^2 \bff{u}_h - |\bff{u}|^2 \bff{u}}{\bff{\theta}} \big|
\\
\nonumber
\norm{\partial_t \bff{\rho}}{\bb{L}^2}^2
\norm{\bff{\theta}}{\bb{L}^2}^2
\norm{\bff{\eta}}{\bb{L}^2}^2
\norm{\bff{\rho}}{\bb{L}^2}^2
\epsilon \norm{\bff{\xi}}{\bb{L}^2}^2
\\
\nonumber
\norm{\bff{\theta}}{\bb{L}^2}^2
\epsilon \norm{\nabla \bff{\theta}}{\bb{L}^2}^2
\epsilon \norm{\bff{\xi}}{\bb{L}^2}^2
\norm{\bff{\theta}}{\bb{L}^2}^2
\epsilon \norm{\nabla \bff{\theta}}{\bb{L}^2}^2
\\
\nonumber
\big( 1+\norm{\bff{\xi}}{\bb{H}^1}^2 \big) \norm{\bff{\theta}}{\bb{L}^2}^2
\epsilon \norm{\nabla \bff{\theta}}{\bb{L}^2}^2
\epsilon \norm{\bff{\xi}}{\bb{L}^2}^2
\\
\lesssim
\big( 1+\norm{\bff{\xi}}{\bb{H}^1}^2 \big) \norm{\bff{\theta}}{\bb{L}^2}^2
\epsilon \norm{\nabla \bff{\theta}}{\bb{L}^2}^2
\epsilon \norm{\bff{\xi}}{\bb{L}^2}^2
\end{align}

for any $\epsilon>0$, where the constant is independent of $h$, $t$, and $T$. Choosing $\epsilon>0$ sufficiently small and integrating over $(0,t)$, we obtain

\begin{align*}
\int_0^t \norm{\nabla \bff{\theta}}{\bb{L}^2}^2 \ds
\int_0^t \norm{\bff{\xi}}{\bb{L}^2}^2 \ds
\lesssim
\norm{\bff{\theta}(0)}{\bb{L}^2}^2
\int_0^t \big( 1+\norm{\bff{\xi}}{\bb{H}^1}^2 \big) \norm{\bff{\theta}}{\bb{L}^2}^2 \ds.
\end{align*}

Note that by Proposition <ref> and (<ref>),

\[
\int_0^t 1+\norm{\bff{\xi}(s)}{\bb{H}^1}^2 \ds
\lesssim
1 + \int_0^t \big( \norm{\bff{H}_h(s)}{\bb{H}^1} + \norm{R_h \bff{H}_h(s)}{\bb{H}^1} \big) \ds
\lesssim
\]

and so Gronwall's inequality applied to (<ref>) yields

\begin{align}\label{equ:theta L2}
\norm{\bff{\theta}(t)}{\bb{L}^2}^2
\int_0^t \norm{\nabla \bff{\theta}(s)}{\bb{L}^2}^2 \ds
\int_0^t \norm{\bff{\xi}(s)}{\bb{L}^2}^2 \ds
\leq
\end{align}

Next, differentiating (<ref>) with respect to $t$, then taking $\bff{\phi}= \lambda_e \bff{\xi}$ yields

\begin{align}\label{equ:difft xi phi=xi}
\nonumber
\frac{\lambda_e}{2} \ddt \norm{\bff{\xi}}{\bb{L}^2}^2
\lambda_e \inpro{\partial_t \bff{\eta}}{\bff{\xi}}
\lambda_e \inpro{\nabla \partial_t \bff{\theta}}{\nabla \bff{\xi}}
\kappa\mu \lambda_e \inpro{\partial_t \bff{\theta}+ \partial_t\bff{\rho}}{\bff{\xi}}
\\
\kappa \lambda_e \inpro{\partial_t(|\bff{u}_h|^2 \bff{u}_h - |\bff{u}|^2 \bff{u})}{\bff{\xi}}
\beta\lambda_e \inpro{\bff{e}(\bff{e}\cdot\partial_t \bff{\theta})}{\bff{\xi}}.
\end{align}

Taking $\bff{\chi}=\kappa\mu\lambda_e \bff{\xi}$ and rearranging the terms, we have

\begin{align}\label{equ:dt theta xi}
\nonumber
\kappa\mu\lambda_e \lambda_r \norm{\bff{\xi}}{\bb{L}^2}^2
\kappa\mu \lambda_e^2 \norm{\nabla \bff{\xi}}{\bb{L}^2}^2
\kappa\mu\lambda_e \inpro{\partial_t \bff{\theta}+\partial_t \bff{\rho}}{\bff{\xi}}
\kappa\mu\lambda_e \lambda_r \inpro{\bff{\eta}}{\bff{\xi}}
\\
\kappa\mu\lambda_e \gamma \inpro{\bff{u}_h \times \bff{H}_h- \bff{u}\times \bff{H}}{\bff{\xi}}.
\end{align}

Adding (<ref>) and (<ref>) yields

\begin{align}\label{equ:add xi dt}
\nonumber
&\frac{\lambda_e}{2} \ddt \norm{\bff{\xi}}{\bb{L}^2}^2
\kappa\mu\lambda_e\lambda_r \norm{\bff{\xi}}{\bb{L}^2}^2
\kappa\mu\lambda_e^2 \norm{\nabla \bff{\xi}}{\bb{L}^2}^2
\\
\nonumber
\lambda_e \inpro{\nabla \partial_t \bff{\theta}}{\nabla \bff{\xi}}
\lambda_e \inpro{\partial_t \bff{\eta}}{\bff{\xi}}
\kappa\mu\lambda_e\lambda_r \inpro{\bff{\eta}}{\bff{\xi}}
\kappa\mu\lambda_e \gamma \inpro{\bff{u}_h \times \bff{H}_h- \bff{u}\times \bff{H}}{\bff{\xi}}
\\
\kappa \lambda_e \inpro{\partial_t(|\bff{u}_h|^2 \bff{u}_h - |\bff{u}|^2 \bff{u})}{\bff{\xi}}
\beta\lambda_e \inpro{\bff{e}(\bff{e}\cdot\partial_t \bff{\theta})}{\bff{\xi}}.
\end{align}

Furthermore, taking $\bff{\chi}= \partial_t \bff{\theta}$ in (<ref>) gives

\begin{align}\label{equ:dt theta chi dt theta}
\norm{\partial_t \bff{\theta}}{\bb{L}^2}^2
\inpro{\partial_t \bff{\rho}}{\partial_t \bff{\theta}}
\lambda_r \inpro{\bff{\xi}+\bff{\eta}}{\partial_t \bff{\theta}}
\lambda_e \inpro{\nabla \bff{\xi}}{\nabla \partial_t \bff{\theta}}
\gamma \inpro{\bff{u}_h \times \bff{H}_h - \bff{u} \times \bff{H}}{\partial_t \bff{\theta}}.
\end{align}

Taking $\bff{\phi}=\lambda_r \partial_t \bff{\theta}$ in (<ref>), we obtain

\begin{align}\label{equ:xi+eta dt theta}
\nonumber
\lambda_r \inpro{\bff{\xi}+\bff{\eta}}{\partial_t \bff{\theta}}
-\frac{\lambda_r}{2} \ddt \norm{\nabla \bff{\theta}}{\bb{L}^2}^2
\kappa\mu \lambda_r \inpro{\bff{\theta}}{\partial_t \bff{\theta}}
\kappa \mu\lambda_r \inpro{\bff{\rho}}{\partial_t \bff{\theta}}
\\
\kappa \lambda_r \inpro{|\bff{u}_h|^2 \bff{u}_h - |\bff{u}|^2 \bff{u}}{\partial_t \bff{\theta}}
\beta\lambda_r \inpro{\bff{e}(\bff{e}\cdot \bff{\theta})}{\partial_t \bff{\theta}}.
\end{align}

Substituting (<ref>) into (<ref>), then adding the result to (<ref>) give

\begin{align*}
&\frac{\lambda_r}{2} \ddt \norm{\nabla \bff{\theta}}{\bb{L}^2}^2
\frac{\lambda_e}{2} \ddt \norm{\bff{\xi}}{\bb{L}^2}^2
\norm{\partial_t \bff{\theta}}{\bb{L}^2}^2
\kappa\mu\lambda_e\lambda_r \norm{\bff{\xi}}{\bb{L}^2}^2
\kappa\mu\lambda_e^2 \norm{\nabla \bff{\xi}}{\bb{L}^2}^2
\\
\inpro{\partial_t \bff{\rho}}{\partial_t \bff{\theta}}
\lambda_e \inpro{\partial_t \bff{\eta}}{\bff{\xi}}
\kappa \mu\lambda_r \inpro{\bff{\theta}}{\partial_t \bff{\theta}}
\kappa \mu\lambda_r \inpro{\bff{\rho}}{\partial_t \bff{\theta}}
\kappa\mu\lambda_e\lambda_r \inpro{\bff{\eta}}{\bff{\xi}}
\\
\gamma \inpro{\bff{u}_h \times \bff{H}_h - \bff{u} \times \bff{H}}{\partial_t \bff{\theta}}
\kappa\mu\lambda_e \gamma \inpro{\bff{u}_h \times \bff{H}_h- \bff{u}\times \bff{H}}{\bff{\xi}}
\\
\kappa \lambda_e \inpro{\partial_t(|\bff{u}_h|^2 \bff{u}_h - |\bff{u}|^2 \bff{u})}{\bff{\xi}}
\kappa \lambda_r \inpro{|\bff{u}_h|^2 \bff{u}_h - |\bff{u}|^2 \bff{u}}{\partial_t \bff{\theta}}
\beta\lambda_r \inpro{\bff{e}(\bff{e}\cdot \bff{\theta})}{\partial_t \bff{\theta}}.
\end{align*}

Applying Hölder's and Young's inequalities to bound the terms on the second line, then using (<ref>), (<ref>) and (<ref>) for the remaining terms yield

\begin{align*}
&\frac{\lambda_r}{2} \ddt \norm{\nabla \bff{\theta}}{\bb{L}^2}^2
\frac{\lambda_e}{2} \ddt \norm{\bff{\xi}}{\bb{L}^2}^2
\norm{\partial_t \bff{\theta}}{\bb{L}^2}^2
\lambda_r \norm{\bff{\xi}}{\bb{L}^2}^2
\lambda_e \norm{\nabla \bff{\xi}}{\bb{L}^2}^2
\\
\norm{\partial_t \bff{\rho}}{\bb{L}^2}
\norm{\partial_t \bff{\theta}}{\bb{L}^2}
\lambda_e \norm{\partial_t \bff{\eta}}{\bb{L}^2} \norm{\bff{\xi}}{\bb{L}^2}
\kappa\mu \lambda_r \norm{\bff{\theta}}{\bb{L}^2} \norm{\partial_t \bff{\theta}}{\bb{L}^2}
\kappa\mu \lambda_r \norm{\bff{\rho}}{\bb{L}^2}
\norm{\partial_t \bff{\theta}}{\bb{L}^2}
\kappa\mu\lambda_e\lambda_r \norm{\bff{\eta}}{\bb{L}^2}
\norm{\bff{\xi}}{\bb{L}^2}
\\
\gamma \big| \inpro{\bff{u}_h \times \bff{H}_h - \bff{u} \times \bff{H}}{\partial_t \bff{\theta}} \big|
\gamma \big| \inpro{\bff{u}_h \times \bff{H}_h- \bff{u}\times \bff{H}}{\bff{\xi}} \big|
\\
\kappa \lambda_e \big| \inpro{\partial_t(|\bff{u}_h|^2 \bff{u}_h - |\bff{u}|^2 \bff{u})}{\bff{\xi}} \big|
\kappa \lambda_r \big| \inpro{|\bff{u}_h|^2 \bff{u}_h - |\bff{u}|^2 \bff{u}}{\partial_t \bff{\theta}} \big|
\beta\lambda_r \big|\inpro{\bff{e}(\bff{e}\cdot \bff{\theta})}{\partial_t \bff{\theta}} \big|
\\
\norm{\bff{\theta}}{\bb{L}^2}^2
\norm{\partial_t \bff{\rho}}{\bb{L}^2}^2
\norm{\partial_t \bff{\eta}}{\bb{L}^2}^2
\norm{\bff{\xi}}{\bb{L}^2}^2
\norm{\bff{\rho}}{\bb{L}^2}^2
\norm{\bff{\eta}}{\bb{L}^2}^2
\epsilon \norm{\partial_t \bff{\theta}}{\bb{L}^2}^2
\\
\norm{\bff{\xi}}{\bb{L}^2}^2
\big( \norm{\bff{\theta}}{\bb{L}^2}^2 + \epsilon \norm{\nabla \bff{\theta}}{\bb{L}^2}^2 \big) \norm{\bff{H}_h}{\bb{H}^1}^2
\epsilon \norm{\partial_t \bff{\theta}}{\bb{L}^2}^2
\\
\norm{\bff{\xi}}{\bb{L}^2}^2
\epsilon \norm{\nabla \bff{\xi}}{\bb{L}^2}^2
\epsilon \norm{\bff{\theta}}{\bb{L}^2}^2
\\
\big( 1 + \norm{\bff{\xi}}{\bb{H}^1}^2 + \norm{\partial_t \bff{u}_h}{\bb{L}^2}^2 \big) \norm{\bff{\theta}}{\bb{H}^1}^2
\norm{\bff{\xi}}{\bb{L}^2}^2
\epsilon \norm{\partial_t \bff{\theta}}{\bb{L}^2}^2
\\
\big( 1+\norm{\partial_t \bff{\theta}}{\bb{L}^2}^2 \big) \norm{\bff{\theta}}{\bb{H}^1}^2
\epsilon \norm{\partial_t \bff{\theta}}{\bb{L}^2}^2
\\
\norm{\bff{\theta}}{\bb{L}^2}^2
\norm{\bff{\xi}}{\bb{L}^2}^2
\epsilon \norm{\partial_t \bff{\theta}}{\bb{L}^2}^2
\big( 1+ \norm{\bff{\xi}}{\bb{H}^1}^2 + \norm{\partial_t \bff{\theta}}{\bb{L}^2}^2
+ \norm{\partial_t \bff{u}_h}{\bb{L}^2}^2 + \norm{\bff{H}_h}{\bb{H}^1}^2 \big)
\norm{\bff{\theta}}{\bb{H}^1}^2
\end{align*}

for any $\epsilon>0$, where the constant is independent of $h$, $t$, and $T$. Choosing $\epsilon>0$ sufficiently small, then integrating over $(0,t)$ and using (<ref>) (and noting (<ref>)), we obtain

\begin{align*}
\norm{\nabla \bff{\theta}(t)}{\bb{L}^2}^2
\norm{\bff{\xi}(t)}{\bb{L}^2}^2
\int_0^t \norm{\partial_t \bff{\theta}(s)}{\bb{L}^2}^2 \ds
\int_0^t \norm{\nabla \bff{\xi}(s)}{\bb{L}^2}^2 \ds
\\
\int_0^t \norm{\bff{\xi}(s)}{\bb{L}^2}^2 \ds
\int_0^t \mathcal{B}(s) \norm{\bff{\theta}(s)}{\bb{H}^1}^2 \ds,
\end{align*}

where $\mathcal{B}(s):= 1+ \norm{\bff{\xi}(s)}{\bb{H}^1}^2 + \norm{\partial_t \bff{\theta}(s)}{\bb{L}^2}^2
+ \norm{\partial_t \bff{u}_h(s)}{\bb{L}^2}^2 + \norm{\bff{H}_h(s)}{\bb{H}^1}^2$. Note that by Proposition <ref>, Proposition <ref>, and inequality (<ref>),

\begin{align*}
\int_0^t \mathcal{B}(s)\, \ds
\leq
\int_0^t \big( 1 + \norm{\bff{H}_h(s)}{\bb{H}^1}^2 + \norm{R_h\bff{H}(s)}{\bb{H}^1}^2
+ \norm{\partial_t \bff{u}_h(s)}{\bb{L}^2}^2 + \norm{R_h\partial_t \bff{u}(s)}{\bb{L}^2}^2 \big)\, \ds
\lesssim 1.
\end{align*}

Therefore, by Gronwall's inequality,

\begin{align}\label{equ:nab theta xi L2 int}
\norm{\nabla \bff{\theta}(t)}{\bb{L}^2}^2
\norm{\bff{\xi}(t)}{\bb{L}^2}^2
\int_0^t \norm{\partial_t \bff{\theta}(s)}{\bb{L}^2}^2 \ds
\int_0^t \norm{\nabla \bff{\xi}(s)}{\bb{L}^2}^2 \ds
\leq
\end{align}

Inequality (<ref>) then follows from (<ref>) and (<ref>).

Suppose now that the triangulation is globally quasi-uniform. Taking $\bff{\phi}= \Delta_h \bff{\theta}$ in (<ref>), rearranging the terms, then applying Young's inequality and (<ref>) give

\begin{align*}
\norm{\Delta_h \bff{\theta}}{\bb{L}^2}^2
\inpro{\bff{\xi}+\bff{\eta}}{\Delta_h \bff{\theta}}
\kappa \mu\inpro{\bff{\theta}+\bff{\rho}}{\Delta_h \bff{\theta}}
\kappa \inpro{|\bff{u}_h|^2 \bff{u}_h - |\bff{u}|^2 \bff{u}}{\Delta_h \bff{\theta}}
\beta \inpro{\bff{e}(\bff{e}\cdot\bff{\theta}}{\Delta_h \bff{\theta}}
\\
\norm{\bff{\xi}}{\bb{L}^2}^2 +\norm{\bff{\theta}}{\bb{L}^2}^2
+ \epsilon \norm{\Delta_h \bff{\theta}}{\bb{L}^2}^2
\lesssim
h^{2(r+1)} + \epsilon \norm{\Delta_h \bff{\theta}}{\bb{L}^2}^2,
\end{align*}

where in the last step we used (<ref>) shown above. Choosing $\epsilon>0$ sufficiently small then yields

\begin{align}\label{equ:Delta theta L2}
\norm{\Delta_h \bff{\theta}}{\bb{L}^2}^2 \lesssim h^{2(r+1)}.
\end{align}

Next, we take $\bff{\chi}=\Delta_h \bff{\xi}$ and apply the same argument (using (<ref>) and (<ref>) this time) to obtain

\begin{align}\label{equ:int Delta xi L2}
\int_0^t \norm{\Delta_h \bff{\xi}(s)}{\bb{L}^2}^2 \ds \lesssim h^{2(r+1)}.
\end{align}

Inequality (<ref>) then follows from (<ref>), (<ref>), and (<ref>).

Finally, setting $\bff{\chi}=-\Delta_h^2 \bff{\xi}$ and applying (<ref>) as necessary, we have

\begin{align}\label{equ:nab Delta xi}
\nonumber
\lambda_r \norm{\Delta_h \bff{\xi}}{\bb{L}^2}^2
\lambda_e \norm{\nabla \Delta_h \bff{\xi}}{\bb{L}^2}^2
-\inpro{\partial_t \nabla \bff{\theta}}{\nabla \Delta_h \bff{\xi}}
\inpro{\nabla \Pi_h \partial_t \bff{\rho}}{\nabla \Delta_h \bff{\xi}}
\\
\gamma \inpro{\nabla \Pi_h(\bff{u}_h\times \bff{H}_h - \bff{u}\times \bff{H})}{\nabla \Delta_h \bff{\xi}}.
\end{align}

Differentiating the second equation in (<ref>) with respect to $t$, then setting $\bff{\phi}=-\Delta_h \bff{\xi}$, we obtain

\begin{align}\label{equ:ddt nab xi}
\nonumber
\frac{1}{2}\ddt \norm{\nabla \bff{\xi}}{\bb{L}^2}^2
\inpro{\partial_t \bff{\eta}}{\Delta_h \bff{\xi}}
\inpro{\partial_t \nabla \bff{\theta}}{\nabla \Delta_h \bff{\xi}}
\kappa \mu\inpro{\partial_t \bff{\theta}+\partial_t \bff{\rho}}{\Delta_h \bff{\xi}}
\\
\kappa \inpro{\partial_t \left(\abs{\bff{u}_h}^2 \bff{u}_h - \abs{\bff{u}}^2 \bff{u}\right)}{\Delta_h \bff{\xi}}
\beta \inpro{\bff{e}(\bff{e}\cdot \partial_t \bff{\theta})}{\Delta_h \bff{\xi}}
\end{align}

Adding (<ref>) and (<ref>), we obtain

\begin{align*}
&\frac{1}{2}\ddt \norm{\nabla \bff{\xi}}{\bb{L}^2}^2
\lambda_r \norm{\Delta_h \bff{\xi}}{\bb{L}^2}^2
\lambda_e \norm{\nabla \Delta_h \bff{\xi}}{\bb{L}^2}^2
\\
\inpro{\partial_t \bff{\eta}}{\Delta_h \bff{\xi}}
\inpro{\nabla \Pi_h \partial_t \bff{\rho}}{\nabla \Delta_h \bff{\xi}}
\kappa \mu\inpro{\partial_t \bff{\theta}+\partial_t \bff{\rho}}{\Delta_h \bff{\xi}}
\beta \inpro{\bff{e}(\bff{e}\cdot \partial_t \bff{\theta})}{\Delta_h \bff{\xi}}
\\
\gamma \inpro{\nabla \Pi_h(\bff{u}_h\times \bff{H}_h - \bff{u}\times \bff{H})}{\nabla \Delta_h \bff{\xi}}
\kappa \inpro{\partial_t \left(\abs{\bff{u}_h}^2 \bff{u}_h - \abs{\bff{u}}^2 \bff{u}\right)}{\Delta_h \bff{\xi}}.
\end{align*}

The first four terms on the right-hand side can be estimated using Young's inequality (together with (<ref>)) in a straightforward manner. The last two terms can be bounded using (<ref>) and (<ref>). We then have for any $\epsilon>0$,

\begin{align*}
&\frac{1}{2}\ddt \norm{\nabla \bff{\xi}}{\bb{L}^2}^2
\lambda_r \norm{\Delta_h \bff{\xi}}{\bb{L}^2}^2
\lambda_e \norm{\nabla \Delta_h \bff{\xi}}{\bb{L}^2}^2
\lesssim
\norm{\partial_t \bff{\theta}}{\bb{L}^2}^2
\norm{\bff{\theta}}{\bb{H}^1}^2
\norm{\bff{\xi}}{\bb{H}^1}^2
\epsilon \norm{\Delta_h \bff{\xi}}{\bb{H}^1}^2
\end{align*}

Choosing $\epsilon>0$ sufficiently small, integrating both sides with respect to $t$ (and noting (<ref>)), we obtain (<ref>). This completes the proof of the proposition.

We are now ready to state the main theorem of this section on the order of convergence for the semi-discrete scheme (<ref>).

Let $(\bff{u},\bff{H})$ be the solution of (<ref>) as described in Section <ref>, and let $(\bff{u}_h,\bff{H}_h)\in \bb{V}_h\times \bb{V}_h$ be the solution of (<ref>) with initial data $\bff{u}_0$. Then

\begin{align*}
\norm{\bff{u}_h - \bff{u}}{L^\infty(\bb{L}^2)}
h \norm{\nabla \bff{u}_h - \nabla \bff{u}}{L^\infty(\bb{L}^2)}
\\
\norm{\bff{H}_h - \bff{H}}{L^\infty(\bb{L}^2)}
h \norm{\nabla \bff{H}_h - \nabla \bff{H}}{L^2(\bb{L}^2)}
\end{align*}

Moreover, if the triangulation is globally quasi-uniform, then for $s=0$ or $1$,

\begin{align*}
\norm{\bff{u}_h - \bff{u}}{L^\infty(\bb{H}^s)}
\norm{\bff{H}_h - \bff{H}}{L^\infty(\bb{H}^s)}
\\
\norm{\bff{u}_h - \bff{u}}{L^\infty(\bb{L}^\infty)}
\norm{\bff{H}_h - \bff{H}}{L^2(\bb{L}^\infty)}
Ch^{r+1} \abs{\ln h}^{\frac{1}{2}}.
\end{align*}

The constant $C$ depends on the coefficients of the equation, $|\mathscr{D}|$, $T$, and $K_0$ (as defined in (<ref>)), but is independent of $h$.

The result follows from Proposition <ref>, equations (<ref>) and (<ref>), estimates (<ref>), (<ref>), and the triangle inequality.

§ A FULLY DISCRETE SCHEME BASED ON THE SEMI-IMPLICIT EULER METHOD

We propose a time-discrete scheme for the LLBar equation using the semi-implicit Euler method. First, we fix some notations. Let $k$ be the time step and $\bff{u}_h^n$ be the approximation in $\bb{V}_h$ of $\bff{u}(t)$ at time $t=t_n:=nk$, where $n=0,1,2,\ldots, \lfloor T/k \rfloor$. For any function $\bff{v}$, we denote $\bff{v}^n:= \bff{v}(t_n)$, and define

\begin{align*}
\delta \bff{v}^{n+1} :=
\frac{\bff{v}^{n+1}-\bff{v}^{n}}{k},\quad
\text{for } n=0,1,\ldots.
\end{align*}

A fully discrete scheme can now be described as follows. We start with $\bff{u}_h^0= R_h \bff{u}(0) \in \bb{V}_h$ for simplicity. For $t_n\in [0,T]$ where $n\in \bb{N}$, given $\bff{u}_h^n \in \bb{V}_h$, define $\bff{u}_h^{n+1}$ and $\bff{H}_h^{n+1}$ by

\begin{align}\label{equ:euler}
\left\{
\begin{alignedat}{1}
\inpro{\delta \bff{u}_h^{n+1}}{\bff{\chi}}
\lambda_r \inpro{\bff{H}_h^{n+1}}{\bff{\chi}}
\lambda_e \inpro{\nabla \bff{H}_h^{n+1}}{\nabla \bff{\chi}}
\gamma \inpro{\bff{u}_h^{n} \times \bff{H}_h^{n+1}}{\bff{\chi}}
\\
\inpro{\bff{H}_h^{n+1}}{\bff{\phi}}
\inpro{\nabla \bff{u}_h^{n+1}}{\nabla \bff{\phi}}
\kappa \mu \inpro{\bff{u}_h^n}{\bff{\phi}}
\kappa \inpro{\abs{\bff{u}_h^{n+1}}^2 \bff{u}_h^{n+1}}{\bff{\phi}}
\beta \inpro{\bff{e}\big(\bff{e}\cdot \bff{u}_h^{n+1}\big)}{\bff{\phi}},
\end{alignedat}
\right.
\end{align}

for all $\bff{\chi}, \bff{\phi}\in \bb{V}_h$.

Note that under the assumptions (<ref>), for $p\in [1,\infty]$ we have

\begin{align}\label{equ:delta un Lp euler}
\norm{\delta \bff{u}^{n+1}}{\bb{L}^p}
\norm{\frac{1}{k} \int_{t_n}^{t_{n+1}} \partial_t \bff{u}(t)\, \dt}{\bb{L}^p}
\leq
\norm{\partial_t \bff{u}(t)}{\bb{L}^p},
\\
\label{equ:delta un min u euler}
\norm{\delta \bff{u}^n-\partial_t \bff{u}^{n+1}}{\bb{L}^p}
\norm{\delta \bff{u}^n-\partial_t \bff{u}^n}{\bb{L}^p}
\norm{\int_{t_n}^{t_{n+1}} \partial_{tt}\bff{u}(t)\, \dt}{\bb{L}^p}
\leq
\end{align}

where Taylor's theorem was used in the last step of (<ref>), and $C$ depends on $\norm{\partial_{tt} \bff{u}}{L^\infty(\bb{L}^p)}$.

We now show that the scheme (<ref>) is well-posed for any time step size $k$.

Let $\bff{u}_h^0 \in \bb{V}_h$ and $k>0$ be given. For each $n\in \bb{N}$, given $\bff{u}_h^n$, there exists unique $\bff{u}_h^{n+1}$ and $\bff{H}_h^{n+1}$ solving (<ref>).

Let $\bb{X}:= \bb{V}_h\times \bb{V}_h$ (which is a subspace of $\bb{H}^1\times \bb{H}^1$), equipped with norm $\norm{(\bff{u},\bff{v})}{\bb{X}} := \norm{\bff{u}}{\bb{H}^1}+ \norm{\bff{v}}{\bb{H}^1}$. Define a nonlinear form $\mathcal{A}:\bb{X}\times \bb{X} \to \bb{R}$ by

\begin{align*}
\nonumber
\mathcal{A}\big((\bff{u},\bff{v}), (\bff{\phi},\bff{\chi})\big)
\inpro{\bff{v}}{\bff{\phi}}
\inpro{\bff{u}}{\bff{\chi}}
k\lambda_r \inpro{\bff{v}}{\bff{\chi}}
k\lambda_e \inpro{\nabla \bff{v}}{\nabla \bff{\chi}}
k\gamma \inpro{\bff{u}_h^n\times \bff{v}}{\bff{\chi}}
\\
\inpro{\nabla \bff{u}}{\nabla \bff{\phi}}
\kappa \inpro{\abs{\bff{u}}^2 \bff{u}}{\bff{\phi}}
\beta\inpro{\bff{e}(\bff{e}\cdot\bff{u})}{\bff{\phi}}.
\end{align*}

and a linear form $f:\bb{X}\to \bb{R}$ by

\[
\kappa \mu \inpro{\bff{u}_h^n}{\bff{\phi}}
\inpro{\bff{u}_h^n}{\bff{\chi}}.
\]

Multiplying the first equation in (<ref>) by $-k$ and adding it to the second equation, we see that solving (<ref>) is equivalent to solving

\begin{align}\label{equ:A equal f}
\mathcal{A}\left((\bff{u}_h^{n+1}, \bff{H}_h^{n+1}), (\bff{\phi},\bff{\chi})\right)
= f\big((\bff{\phi},\bff{\chi})\big),
\quad
\forall (\bff{\phi},\bff{\chi})\in \bb{X}.
\end{align}

For each fixed $(\bff{u},\bff{v})\in \bb{X}$, the map $(\bff{\phi},\bff{\chi})\mapsto \mathcal{A}\big((\bff{u},\bff{v}), (\bff{\phi},\bff{\chi})\big)$ is a bounded linear functional. Thus, there exists a map $T:\bb{X}\to \bb{X}^\ast$ defined by

\[
(\bff{u},\bff{v})\mapsto T((\bff{u},\bff{v})), \quad \text{such that} \quad \mathcal{A}\big((\bff{u},\bff{v}), (\bff{\phi},\bff{\chi})\big)
= \inpro{T\big((\bff{u},\bff{v})\big)}{(\bff{\phi},\bff{\chi})},
\quad \forall (\bff{\phi},\bff{\chi})\in \bb{X},
\]

where $\inpro{\cdot}{\cdot}$ denotes the duality pairing which extends the usual $\bb{L}^2$-inner product.
We aim to use the Browder–Minty theorem to deduce the existence and uniqueness of (<ref>). To this end, it remains to show $T$ is bounded, continuous, strictly monotone, and coercive.

Boundedness: We have by Hölder's inequality and Sobolev embedding $\bb{H}^1\subset \bb{L}^6$,

\begin{align*}
\norm{T\big((\bff{u},\bff{v})\big)}{\bb{X}^\ast}
\sup \left\{
\left| \mathcal{A}\big((\bff{u},\bff{v}), (\bff{\phi},\bff{\chi})\big) \right| :
(\bff{\phi},\bff{\chi})\in \bb{X} \text{ and } \norm{(\bff{\phi},\bff{\chi})}{\bb{X}}\leq 1 \right\}
\\
C\norm{\bff{u}}{\bb{H}^1}+ \left(1+k\lambda_r+k\lambda_e+k\gamma \norm{\bff{u}_h^n}{\bb{H}^1}\right) \norm{\bff{v}}{\bb{H}^1} + \norm{\bff{u}}{\bb{H}^1}^3
\leq
C \norm{(\bff{u},\bff{v})}{\bb{X}} \left(1+ \norm{(\bff{u},\bff{v})}{\bb{X}}^2\right).
\end{align*}

This shows $T$ maps bounded sets in $\bb{X}$ into bounded sets in $\bb{X}^\ast$.

Continuity: It is clear that $(\bff{u}_n, \bff{v}_n)\to (\bff{u},\bff{v})$ strongly in $\bb{X}$ implies $\norm{T\big((\bff{u}_n,\bff{v}_n)\big) - T\big((\bff{u}, \bff{v})\big)}{\bb{X}^\ast} \to 0$ as $n\to\infty$. The details are omitted.

Strict monotonicity: Let $(\bff{u},\bff{v}), (\bff{u}',\bff{v}')\in \bb{X}$. Then we have

\begin{align*}
&\inpro{T\big((\bff{u},\bff{v})\big)- T\big((\bff{u}',\bff{v}')\big)}{(\bff{u},\bff{v})-(\bff{u}',\bff{v}')}
\\
\mathcal{A}\big((\bff{u},\bff{v}), (\bff{u}-\bff{u}',\bff{v}-\bff{v}')\big)
\mathcal{A}\big((\bff{u}',\bff{v}'), (\bff{u}-\bff{u}',\bff{v}-\bff{v}')\big)
\\
k\lambda_r \norm{\bff{v}-\bff{v}'}{\bb{L}^2}^2
k \lambda_e \norm{\nabla \bff{v}-\nabla \bff{v}'}{\bb{L}^2}^2
\norm{\nabla \bff{u}-\nabla \bff{u}'}{\bb{L}^2}^2
\beta \norm{\bff{e}\cdot (\bff{u}-\bff{u}')}{L^2}^2
\\
\kappa \norm{\abs{\bff{u}}^2- \abs{\bff{u}'}^2}{\bb{L}^2}^2
\kappa \norm{\abs{\bff{u}} \abs{\bff{u}-\bff{u}'}}{\bb{L}^2}^2
\kappa \norm{\abs{\bff{u}'} \abs{\bff{u}-\bff{u}'}}{\bb{L}^2}^2 \geq 0,
\end{align*}

where in the last step we used (<ref>). Moreover, equality holds if and only if $(\bff{u},\bff{v})=(\bff{u}',\bff{v}')$. This shows the strict monotonicity of $T$.

Coercivity: Let $(\bff{u},\bff{v})\in \bb{X}$. We have

\begin{align*}
\inpro{T\big((\bff{u},\bff{v})\big)}{(\bff{u},\bff{v})}
\mathcal{A}\big((\bff{u},\bff{v}), (\bff{u},\bff{v})\big)
\\
k\lambda_r \norm{\bff{v}}{\bb{L}^2}^2
k\lambda_e \norm{\nabla \bff{v}}{\bb{L}^2}^2
\norm{\nabla \bff{u}}{\bb{L}^2}^2
\kappa \norm{\bff{u}}{\bb{L}^4}^4
\beta \norm{\bff{e}\cdot\bff{u}}{L^2}^2
\\
C_1 \norm{\bff{v}}{\bb{H}^1}^2
C_2 \norm{\bff{u}}{\bb{H}^1}^2
\kappa \abs{\mathscr{D}}^{-4}
\\
\min\{C_1,C_2\} \norm{(\bff{u},\bff{v})}{\bb{X}}^2 - \kappa\abs{\mathscr{D}}^{-4},
\end{align*}

where $C_1:=\min\{k\lambda_r,k\lambda_e\}$ and $C_2:= \min\{1,\kappa\abs{\mathscr{D}}^{-4}\}$ are both positive constants, and we also used the inequality

\[
\norm{\bff{u}}{\bb{L}^4}^4 \geq \abs{\mathscr{D}}^{-4} \left(\norm{\bff{u}}{\bb{L}^2}^2-1\right),
\]

which follows from Hölder's and Young's inequalities. Therefore,

\begin{align*}
\frac{\inpro{T\big((\bff{u},\bff{v})\big)}{(\bff{u},\bff{v})}}{\norm{(\bff{u},\bff{v})}{\bb{X}}}
\to \infty
\quad \text{as } \norm{(\bff{u},\bff{v})}{\bb{X}}\to \infty,
\end{align*}

showing the coercivity of $T$.

The existence and uniqueness of  $\bff{u}_h^{n+1}$ and $\bff{H}_h^{n+1}$ solving (<ref>) then follows from the Browder–Minty theorem. This completes the proof of the proposition.

Next, we show some stability results. The following proposition shows that $\bff{u}_h^n$ is stable in $\ell^\infty(\bb{H}^1)$ norm, while $\bff{H}_h^n$ is stable in $\ell^2(\bb{H}^1)$ norm (over an arbitrary number of iterations $n$, even if $T=\infty$). Moreover, the energy dissipation property (<ref>) is satisfied unconditionally.

Let $\bff{u}_h^0 \in \bb{V}_h$ be given and let $\left(\bff{u}_h^n, \bff{H}_h^n\right)$ be defined by
(<ref>). Then for any $k>0$ and $n\in \bb{N}$,

\begin{align}\label{equ:ene decrease euler}
\mathcal{E}(\bff{u}_h^{n+1}) \leq \mathcal{E}(\bff{u}_h^n),
\end{align}

where $\mathcal{E}$ was defined in (<ref>). Moreover,

\begin{align}\label{equ:stab L4 H1 euler}
\norm{\bff{u}_h^n}{\bb{L}^4}^4
\norm{\bff{u}_h^n}{\bb{H}^1}^2
k\lambda_r \sum_{m=1}^n \norm{\bff{H}_h^{m}}{\bb{L}^2}^2
k\lambda_e \sum_{m=1}^n \norm{\nabla \bff{H}_h^{m}}{\bb{L}^2}^2
\leq
C \left( \norm{\bff{u}_h^0}{\bb{H}^1}^2 + |\mathscr{D}| \right),
\end{align}

where $C$ depends only on $\kappa$ and $\mu$.

Setting $\bff{\chi}=\bff{H}_h^{n+1}$ in (<ref>) gives

\begin{align}\label{equ:ene stab Un Hn euler}
\inpro{\delta \bff{u}_h^{n+1}}{\bff{H}_h^{n+1}}
\lambda_r \norm{\bff{H}_h^{n+1}}{\bb{L}^2}^2
\lambda_e \norm{\nabla \bff{H}_h^{n+1}}{\bb{L}^2}^2,
\end{align}

while setting $\bff{\phi}=\delta \bff{u}_h^{n+1}$ gives

\begin{align}\label{equ:ene stab Hn Un euler}
\nonumber
\inpro{\bff{H}_h^{n+1}}{\delta \bff{u}_h^{n+1}}
\frac{1}{2k} \big(\norm{\nabla \bff{u}_h^{n+1}}{\bb{L}^2}^2 - \norm{\nabla \bff{u}_h^n}{\bb{L}^2}^2 \big)
\frac{1}{2k}\norm{\nabla \bff{u}_h^{n+1}- \nabla \bff{u}_h^n}{\bb{L}^2}^2
\\
\nonumber
\frac{\kappa}{4k} \left(\norm{\abs{\bff{u}_h^{n+1}}^2 - \mu}{\bb{L}^2}^2 - \norm{\abs{\bff{u}_h^n}^2 -\mu}{\bb{L}^2}^2 \right)
\frac{\kappa}{4k} \norm{\abs{\bff{u}_h^{n+1}}^2 - \abs{\bff{u}_h^n}^2}{\bb{L}^2}^2
\\
\frac{\kappa k}{2} \norm{\abs{\bff{u}_h^{n+1}} \abs{\delta \bff{u}_h^{n+1}}}{\bb{L}^2}^2
\frac{\mu k}{2} \norm{\bff{u}_h^{n+1}-\bff{u}_h^n}{\bb{L}^2}^2
\nonumber \\
\frac{\beta}{2k} \left(\norm{\bff{e}\cdot \bff{u}_h^{n+1}}{L^2}^2 - \norm{\bff{e}\cdot \bff{u}_h^n}{L^2}^2 \right)
\frac{\beta}{2k} \norm{\bff{e}\cdot (\bff{u}_h^{n+1}-\bff{u}_h^n)}{L^2}^2.
\end{align}

Substituting (<ref>) into (<ref>), using identities (<ref>) and (<ref>), and rearranging the terms yield

\begin{align*}
&\frac{1}{2} \big(\norm{\nabla \bff{u}_h^{n+1}}{\bb{L}^2}^2 - \norm{\nabla \bff{u}_h^n}{\bb{L}^2}^2 \big)
\frac{\kappa}{4} \left(\norm{\abs{\bff{u}_h^{n+1}}^2 - \mu}{\bb{L}^2}^2 - \norm{\abs{\bff{u}_h^n}^2 -\mu}{\bb{L}^2}^2 \right)
\frac{\beta}{2} \left(\norm{\bff{e}\cdot \bff{u}_h^{n+1}}{L^2}^2 - \norm{\bff{e}\cdot \bff{u}_h^n}{L^2}^2 \right)
\\
\frac{1}{2}\norm{\nabla \bff{u}_h^{n+1}- \nabla \bff{u}_h^n}{\bb{L}^2}^2
\frac{\kappa}{4} \norm{\abs{\bff{u}_h^{n+1}}^2 - \abs{\bff{u}_h^n}^2}{\bb{L}^2}^2
\frac{\kappa k^2}{2} \norm{\abs{\bff{u}_h^{n+1}} \abs{\delta \bff{u}_h^{n+1}}}{\bb{L}^2}^2
\frac{\mu k^2}{2} \norm{\bff{u}_h^{n+1}-\bff{u}_h^n}{\bb{L}^2}^2
\\
\frac{\beta}{2} \norm{\bff{e}\cdot (\bff{u}_h^{n+1}-\bff{u}_h^n)}{L^2}^2
k\lambda_r \norm{\bff{H}_h^{n+1}}{\bb{L}^2}^2
k\lambda_e \norm{\nabla \bff{H}_h^{n+1}}{\bb{L}^2}^2
= 0,
\end{align*}

which implies (<ref>). Finally, applying similar arguments as (<ref>) and (<ref>) yield (<ref>).

We also derive the stability of $\bff{H}_h^{n}$ in $\ell^\infty(\bb{L}^2)$ norm. The following identity will be used in the proof:

\begin{align}
\label{equ:identity psi euler}
\abs{\bff{u}_h^{n+1}}^2 \bff{u}_h^{n+1} - \abs{\bff{u}_h^n}^2 \bff{u}_h^n
\abs{\bff{u}_h^{n+1}}^2 \left(\bff{u}_h^{n+1}-\bff{u}_h^n\right)
\left(\bff{u}_h^{n+1}-\bff{u}_h^n\right) \cdot \left(\bff{u}_h^{n+1}+\bff{u}_h^n\right) \bff{u}_h^n.
\end{align}

Let $\bff{u}_h^0 \in \bb{V}_h$ be given and let $\bff{u}_h^n$, $\bff{H}_h^{n}$ be defined by (<ref>). Then for any $k>0$ and $n\in \bb{N}$,

\begin{align}\label{equ:Hn L2 stab euler}
\norm{\bff{H}_h^{n}}{\bb{L}^2}^2
k \sum_{m=1}^n \norm{\delta \bff{u}_h^m}{\bb{L}^2}^2
\leq
\end{align}

where $C$ depends on the coefficients of the equation, $K_0$, and $|\mathscr{D}|$, but is independent of $n$, $k$, and $T$.

Taking $\bff{\chi}=\delta \bff{u}_h^{n+1}$ in (<ref>) gives

\begin{align}\label{equ:delta n half euler}
\norm{\delta \bff{u}_h^{n+1}}{\bb{L}^2}^2
\lambda_r \inpro{\bff{H}_h^{n+1}}{\delta \bff{u}_h^{n+1}}
\lambda_e \inpro{\nabla \bff{H}_h^{n+1}}{\nabla \delta \bff{u}_h^{n+1}}
\gamma \inpro{\bff{u}_h^{n} \times \bff{H}_h^{n+1}}{\delta \bff{u}_h^{n+1}}.
\end{align}

Subtracting the second equation in (<ref>) at time step $n-1$ from the same equation at time step $n$, then setting $\bff{\phi}= \lambda_e \bff{H}_h^{n+1}/k$ gives

\begin{align}\label{equ:H n half euler}
\nonumber
&\frac{\lambda_e}{2k} \left( \norm{\bff{H}_h^{n+1}}{\bb{L}^2}^2 - \norm{\bff{H}_h^n}{\bb{L}^2}^2 \right)
\frac{\lambda_e}{2k} \norm{\bff{H}_h^{n+1}-\bff{H}_h^n}{\bb{L}^2}^2
\\
\nonumber
\kappa\mu \lambda_e \inpro{\bff{H}_h^{n+1}}{\delta \bff{u}_h^{n+1}}
\lambda_e \inpro{\nabla \bff{H}_h^{n+1}}{\nabla \delta \bff{u}_h^{n+1}}
\beta\lambda_e \inpro{\bff{e}(\bff{e}\cdot \delta \bff{u}_h^{n+1})}{\bff{H}_h^{n+1}}
\\
\kappa \lambda_e \inpro{\abs{\bff{u}_h^{n+1}}^2 \delta\bff{u}_h^{n+1} + \delta\bff{u}_h^{n+1} \cdot (\bff{u}_h^{n+1}+\bff{u}_h^n) \bff{u}_h^n}{\bff{H}_h^{n+1}},
\end{align}

where we used (<ref>) for the left-hand side, and (<ref>) for the last inner product. Adding (<ref>) and (<ref>) gives

\begin{align*}
&\frac{\lambda_e}{2k} \left( \norm{\bff{H}_h^{n+1}}{\bb{L}^2}^2 - \norm{\bff{H}_h^n}{\bb{L}^2}^2 \right)
\frac{\lambda_e}{2k} \norm{\bff{H}_h^{n+1}-\bff{H}_h^n}{\bb{L}^2}^2
\norm{\delta \bff{u}_h^{n+1}}{\bb{L}^2}^2
\\
(\lambda_r+\kappa\mu \lambda_e) \inpro{\bff{H}_h^{n+1}}{\delta \bff{u}_h^{n+1}}
\beta\lambda_e \inpro{\bff{e}(\bff{e}\cdot \delta \bff{u}_h^{n+1})}{\bff{H}_h^{n+1}}
\gamma \inpro{\bff{u}_h^{n} \times \bff{H}_h^{n+1}}{\delta \bff{u}_h^{n+1}}
\\
\kappa \lambda_e \inpro{\abs{\bff{u}_h^{n+1}}^2 \delta\bff{u}_h^{n+1} + \delta\bff{u}_h^{n+1} \cdot (\bff{u}_h^{n+1}+\bff{u}_h^n) \bff{u}_h^n}{\bff{H}_h^{n+1}}
\\
\end{align*}

We will estimate each term on the right-hand side by applying Hölder's and Young's inequalities. Firstly, for the terms $S_1$ and $S_2$,

\begin{align}\label{equ:S1 euler}
|S_1| + |S_2|
C \norm{\bff{H}_h^{n+1}}{\bb{L}^2}^2
\frac{1}{4} \norm{\delta \bff{u}_h^{n+1}}{\bb{L}^2}^2.
\end{align}

For the third term, similarly we have

\begin{align}\label{equ:S2 euler}
\norm{\bff{u}_h^n}{\bb{L}^4} \norm{\bff{H}_h^{n+1}}{\bb{L}^4} \norm{\delta \bff{u}_h^{n+1}}{\bb{L}^2}
\leq
C \norm{\bff{H}_h^{n+1}}{\bb{H}^1}^2
\frac{1}{4} \norm{\delta \bff{u}_h^{n+1}}{\bb{L}^2}^2,
\end{align}

where in the last step we used (<ref>) and the Sobolev embedding $\bb{H}^1 \hookrightarrow \bb{L}^4$. Similarly for the last term,

\begin{align}\label{equ:S3 euler}
\nonumber
C \left( \norm{\bff{u}_h^{n+1}}{\bb{L}^6}^2 \norm{\delta \bff{u}_h^{n+1}}{\bb{L}^2} + \norm{\delta \bff{u}_h^{n+1}}{\bb{L}^2} \norm{\bff{u}_h^{n+1}+\bff{u}_h^n}{\bb{L}^6} \norm{\bff{u}_h^n}{\bb{L}^6} \right) \norm{\bff{H}_h^{n+1}}{\bb{L}^6}
\\
C \norm{\bff{H}_h^{n+1}}{\bb{H}^1}^2
\frac{1}{4} \norm{\delta \bff{u}_h^{n+1}}{\bb{L}^2}^2.
\end{align}

Altogether, upon rearranging and summing the terms, (<ref>), (<ref>), and (<ref>) imply

\begin{align*}
\lambda_e \norm{\bff{H}_h^n}{\bb{L}^2}^2
\lambda_e \sum_{m=0}^n \norm{\bff{H}_h^{m+1}-\bff{H}_h^m}{\bb{L}^2}^2
k \sum_{m=0}^n \norm{\delta \bff{u}_h^n}{\bb{L}^2}^2
\lambda_e \norm{\bff{H}_h^0}{\bb{L}^2}^2
Ck \sum_{m=0}^n \norm{\bff{H}_h^{m+1}}{\bb{H}^1}^2
\leq C,
\end{align*}

as required.

The above proposition implies the stability of $\bff{u}_h^n$ in $\ell^\infty(\bb{L}^\infty)$ norm under an additional assumption that the triangulation is globally quasi-uniform.

Let $\bff{u}_h^0 \in \bb{V}_h$ be given and let $\left(\bff{u}_h^n, \bff{H}_h^n\right)$ be defined by (<ref>). Then for any $k>0$ and $n\in \bb{N}$,

\begin{align}\label{equ:stab Delta uh euler}
\norm{\Delta_h \bff{u}_h^n}{\bb{L}^2}^2 \leq C.
\end{align}

Moreover, if the triangulation $\mathcal{T}_h$ is (globally) quasi-uniform, then

\begin{align}\label{equ:stab uh L infty euler}
\norm{\bff{u}_h^n}{\bb{L}^\infty}^2 \leq C.
\end{align}

Here, the constant $C$ depends on the coefficients of the equation, $K_0$, and $|\mathscr{D}|$, but is independent of $n$, $k$, and $T$.

Taking $\bff{\phi}=\Delta_h \bff{u}_h^{n+1}$ and applying Young's and Hölder's inequalities, we have

\begin{align*}
\norm{\Delta_h \bff{u}_h^{n+1}}{\bb{L}^2}^2
\inpro{\bff{H}_h^{n+1}}{\Delta_h \bff{u}_h^{n+1}}
\kappa \mu \inpro{\bff{u}_h^n}{\Delta_h \bff{u}_h^{n+1}}
\beta \inpro{\bff{e}(\bff{e}\cdot \bff{u}_h^{n+1})}{\Delta_h \bff{u}_h^{n+1}}
\\
\kappa \inpro{\abs{\bff{u}_h^{n+1}}^2 \bff{u}_h^{n+1}}{\Delta_h \bff{u}_h^{n+1}}
\\
\frac{1}{8} \norm{\Delta_h \bff{u}_h^{n+1}}{\bb{L}^2}^2
(\kappa^2 \mu^2+\beta^2) \norm{\bff{u}_h^n}{\bb{L}^2}^2
\frac{3}{8} \norm{\Delta_h \bff{u}_h^n}{\bb{L}^2}^2
\kappa^2 \norm{\bff{u}_h^{n+1}}{\bb{L}^6}^6.
\end{align*}

Therefore, rearranging the terms, using the Sobolev embedding $\bb{H}^1 \hookrightarrow \bb{L}^6$ (noting (<ref>) and (<ref>)), we infer (<ref>).
Inequality (<ref>) then follows from (<ref>), completing the proof of the proposition.

Furthermore, in dimensions 1 and 2 (or $d=3$ with an additional quasi-uniformity assumption on the triangulation), we can derive the stability of $\bff{H}_h^n$ in $\ell^\infty(\bb{H}^1)$.

Let $\bff{u}_h^0 \in \bb{V}_h$ be given and let $\bff{u}_h^n$, $\bff{H}_h^{n}$ be defined by (<ref>). Suppose that one of the following holds:

* $d=1$ or $2$,

* $d=3$ and the triangulation is globally quasi-uniform.

Then for any $k>0$ and $n\in \bb{N}$,

\begin{align}\label{equ:dt un stab euler}
\norm{\delta \bff{u}_h^n}{\bb{L}^2}^2
k \sum_{m=1}^n \norm{\delta \bff{H}_h^m}{\bb{L}^2}^2
k \sum_{m=1}^n \norm{\nabla \delta \bff{u}_h^m}{\bb{L}^2}^2
\leq
\end{align}

where $C$ depends on the coefficients of the equation, $K_0$, and $|\mathscr{D}|$, but is independent of $n$, $k$, and $T$.

For $n\in \bb{N}$, subtracting the first equation in (<ref>) at time step $n$ from the corresponding equation at time step $n+1$, dividing by $k$, then setting $\bff{\chi}=\delta\bff{u}_h^{n+1}$ give

\begin{align}\label{equ:delta uh n1}
\nonumber
\frac{1}{2k} \left(\norm{\delta \bff{u}_h^{n+1}}{\bb{L}^2}^2 - \norm{\delta \bff{u}_h^n}{\bb{L}^2}^2\right)
\frac{1}{2k} \norm{\delta \bff{u}_h^{n+1}-\delta \bff{u}_h^n}{\bb{L}^2}^2
\lambda_r \inpro{\delta \bff{H}_h^{n+1}}{\delta \bff{u}_h^{n+1}}
\lambda_e \inpro{\nabla \delta \bff{H}_h^{n+1}}{\nabla \delta \bff{u}_h^{n+1}}
\\
\gamma \inpro{\delta \bff{u}_h^n\times \bff{H}_h^{n+1} + \bff{u}_h^{n-1} \times \delta \bff{H}_h^{n+1}}{\delta \bff{u}_h^{n+1}}.
\end{align}

Applying the same operations to the second equation in (<ref>), then successively setting $\bff{\phi}=\lambda_e \delta \bff{H}_h^{n+1}$ and $\bff{\phi}=\Lambda \delta \bff{u}_h^{n+1}$ (where $\Lambda:= \lambda_r+\kappa \mu \lambda_e$) yield

\begin{align}
\label{equ:lambda e delta Hh n1}
\nonumber
\lambda_e \norm{\delta \bff{H}_h^{n+1}}{\bb{L}^2}^2
\lambda_e \inpro{\nabla \delta \bff{u}_h^{n+1}}{\nabla \delta \bff{H}_h^{n+1}}
\kappa \mu \lambda_e \inpro{\delta\bff{u}_h^n}{\delta \bff{H}_h^{n+1}}
\beta\lambda_e \inpro{\bff{e}(\bff{e}\cdot \delta\bff{u}_h^{n+1})}{\delta\bff{H}_h^{n+1}}
\\
\kappa \lambda_e \inpro{\abs{\bff{u}_h^{n+1}}^2 \delta \bff{u}_h^{n+1} + \left( (\bff{u}_h^{n+1}+\bff{u}_h^n)\cdot \delta \bff{u}_h^{n+1}\right) \bff{u}_h^n}{\delta \bff{H}_h^{n+1}},
\end{align}

and (noting the identity (<ref>)),

\begin{align}
\label{equ:delta Hh delta uh}
\Lambda \inpro{\delta \bff{H}_h^{n+1}}{\delta \bff{u}_h^{n+1}}
-\Lambda\left(\norm{\nabla \delta \bff{u}_h^{n+1}}{\bb{L}^2}^2 - \kappa \mu \inpro{\delta \bff{u}_h^n}{\delta \bff{u}_h^{n+1}}
+ P \right).
\end{align}

In the last step, we denoted $\bff{u}_h^{n+\frac{1}{2}}:= \frac{1}{2} (\bff{u}_h^{n+1}+\bff{u}_h^n)$ and

\[
P:= \kappa \norm{\bff{u}_h^{n+\frac{1}{2}}\cdot \delta \bff{u}_h^{n+1}}{\bb{L}^2}^2
+ \frac{\kappa}{2} \norm{\abs{\bff{u}_h^{n+1}} \abs{\delta \bff{u}_h^{n+1}}}{\bb{L}^2}^2
+ \frac{\kappa}{2} \norm{\abs{\bff{u}_h^n} \abs{\delta \bff{u}_h^{n+1}}}{\bb{L}^2}^2
+ \beta\norm{\bff{e}\cdot \delta\bff{u}_h^{n+1}}{L^2}^2.
\]

Adding (<ref>), (<ref>), and (<ref>), then applying Young's inequality, we obtain for any $\epsilon >0$,

\begin{align*}
E&:= \frac{1}{2k} \left(\norm{\delta \bff{u}_h^{n+1}}{\bb{L}^2}^2 - \norm{\delta \bff{u}_h^n}{\bb{L}^2}^2\right)
\frac{1}{2k} \norm{\delta \bff{u}_h^{n+1}-\delta \bff{u}_h^n}{\bb{L}^2}^2
\lambda_e \norm{\delta \bff{H}_h^{n+1}}{\bb{L}^2}^2
\Lambda \norm{\nabla \delta \bff{u}_h^{n+1}}{\bb{L}^2}^2
\\
\Lambda \kappa \mu \inpro{\delta \bff{u}_h^n}{\delta \bff{u}_h^{n+1}}
\kappa \lambda_e \inpro{\abs{\bff{u}_h^{n+1}}^2 \delta \bff{u}_h^{n+1} + \left( (\bff{u}_h^{n+1}+\bff{u}_h^n)\cdot \delta \bff{u}_h^{n+1}\right) \bff{u}_h^n}{\delta \bff{H}_h^{n+1}}
\\
\beta\lambda_e \inpro{\bff{e}(\bff{e}\cdot \delta\bff{u}_h^{n+1})}{\delta\bff{H}_h^{n+1}}
\gamma \inpro{\delta \bff{u}_h^n\times \bff{H}_h^{n+1} + \bff{u}_h^{n-1} \times \delta \bff{H}_h^{n+1}}{\delta \bff{u}_h^{n+1}}
\\
C \norm{\delta \bff{u}_h^n}{\bb{L}^2} \norm{\delta \bff{u}_h^{n+1}}{\bb{L}^2}
C \left(\norm{\bff{u}_h^{n+1}}{\bb{L}^8}^2 + \norm{\bff{u}_h^n}{\bb{L}^8}^2\right) \norm{\delta \bff{u}_h^{n+1}}{\bb{L}^4} \norm{\delta \bff{H}_h^{n+1}}{\bb{L}^2}
\\
C \norm{\delta \bff{u}_h^{n+1}}{\bb{L}^2} \norm{\delta \bff{H}_h^{n+1}}{\bb{L}^2}
C \norm{\delta \bff{u}_h^n}{\bb{L}^2} \norm{\bff{H}_h^{n+1}}{\bb{L}^4} \norm{\delta \bff{u}_h^{n+1}}{\bb{L}^4}
\\
C \norm{\bff{u}_h^{n-1}}{\bb{L}^4} \norm{\delta \bff{H}_h^{n+1}}{\bb{L}^2} \norm{\delta \bff{u}_h^{n+1}}{\bb{L}^4}.
\end{align*}

Note that $P\geq 0$. The expression in the last step can be estimated as follows.

Case 1: $d=1$ or $2$. By Young's inequality, Sobolev embedding $\bb{H}^1\hookrightarrow \bb{L}^8$, and (<ref>), we have

\begin{align}\label{equ:E dun}
C \norm{\delta \bff{u}_h^n}{\bb{L}^2}^2
C\norm{\delta \bff{u}_h^{n+1}}{\bb{L}^2}^2
\epsilon \norm{\nabla \delta \bff{u}_h^{n+1}}{\bb{L}^2}^2
\epsilon \norm{\delta \bff{H}_h^{n+1}}{\bb{L}^2}^2
C \norm{\bff{H}_h^{n+1}}{\bb{L}^4}^2 \norm{\delta \bff{u}_h^n}{\bb{L}^2}^2.
\end{align}

Choosing $\epsilon>0$ sufficiently small, multiplyting by $2k$, rearranging the terms, then summing over $m\in \{1,2,\ldots,n-1\}$, we infer that

\begin{align*}
\norm{\delta \bff{u}_h^n}{\bb{L}^2}^2
k \sum_{m=1}^{n-1} \left(
\norm{\delta \bff{H}_h^{m+1}}{\bb{L}^2}^2
\norm{\nabla \delta \bff{u}_h^{m+1}}{\bb{L}^2}^2 \right)
\leq
\norm{\delta \bff{u}_h^1}{\bb{L}^2}^2
Ck\sum_{m=1}^{n-1} \left(1+ C \norm{\bff{H}_h^{m+1}}{\bb{L}^4}^2\right) \norm{\delta \bff{u}_h^m}{\bb{L}^2}^2.
\end{align*}

Note that by taking $n=0$ and $\bff{\chi}=\delta \bff{u}_h^1$ in the first equation of (<ref>), we have

\[
\norm{\delta \bff{u}_h^1}{\bb{L}^2}^2
\lesssim
\norm{\bff{H}_h^1}{\bb{H}^1}^2 + \norm{\Delta_h \bff{H}_h^1}{\bb{L}^2}^2 \lesssim 1.
\]

Therefore, the discrete Gronwall inequality and (<ref>) then imply (<ref>).

Case 2: $d=3$ and the triangulation is quasi-uniform. In this case, we apply Young's inequality, (<ref>), and (<ref>) instead to infer (<ref>).
The same conclusion is then attained.

Suppose that the hypotheses of Proposition <ref> hold. Then for any $k>0$ and $n\in \bb{N}$,

\begin{align}\label{equ:Hn H1 stab euler}
\norm{\nabla \bff{H}_h^n}{\bb{L}^2}^2
\leq
\end{align}

where $C$ depends on the coefficients of the equation, $K_0$, and $|\mathscr{D}|$, but is independent of $n$, $k$, and $T$.

Setting $\bff{\chi}=\bff{H}_h^{n+1}$ in (<ref>) gives

\begin{align*}
\lambda_e \norm{\nabla \bff{H}_h^{n+1}}{\bb{L}^2}^2
\inpro{\delta \bff{u}_h^{n+1}}{\bff{H}_h^{n+1}}
\lambda_r \norm{\bff{H}_h^{n+1}}{\bb{L}^2}^2
\leq
C\norm{\delta \bff{u}_h^{n+1}}{\bb{L}^2}^2
\leq C,
\end{align*}

where in the last step we used (<ref>) and (<ref>).

Proposition <ref> and <ref> also hold under a technical smallness assumption on the initial data similar to Proposition <ref>. We will not elaborate them further for brevity.

To derive optimal error estimates for the fully discrete scheme, we begin with some preparatory results. As done previously, we write

\begin{align}
\label{equ:Un utn}
\bff{u}_h^n-\bff{u}(t_n)
\big(\bff{u}_h^n- R_h\bff{u}(t_n) \big)
\big(R_h \bff{u}(t_n)-\bff{u}(t_n) \big)
\bff{\theta}^n + \bff{\rho}^n,
\\
\label{equ:Hn Htn}
\bff{H}_h^n-\bff{H}(t_n)
\big(\bff{H}_h^n- R_h\bff{H}(t_n) \big)
\big(R_h \bff{H}(t_n)-\bff{H}(t_n) \big)
\bff{\xi}^n + \bff{\eta}^n.
\end{align}

Recall that by the definition of Ritz projection,

\begin{align}\label{equ:nab Ritz disc}
\inpro{\nabla \bff{\rho}^n}{\nabla \bff{\chi}}
\inpro{\nabla \bff{\eta}^n}{\nabla \bff{\chi}}
\quad\forall \bff{\chi} \in \bb{V}_h.
\end{align}

The following lemmas are needed to bound the nonlinear terms in the main theorem.

Let $\epsilon>0$ be arbitrary. The following inequality holds:

\begin{align}
\label{equ:cross theta est euler}
\nonumber
\left| \inpro{\bff{u}_h^{n} \times \bff{H}_h^{n+1} - \bff{u}^{n+1} \times \bff{H}^{n+1}}{\bff{\theta}^{n+1}} \right|
h^{2(r+1)} + k^2+ \norm{\bff{\theta}^{n+1}}{\bb{L}^2}^2
\\
+ \epsilon \norm{\bff{\theta}^n}{\bb{L}^2}^2
+ \epsilon \norm{\nabla \bff{\theta}^{n+1}}{\bb{L}^2}^2
+ \epsilon \norm{\bff{\xi}^{n+1}}{\bb{L}^2}^2,
\\
\nonumber
\label{equ:cross xi euler}
\big| \inpro{\bff{u}_h^n \times \bff{H}_h^{n+1}- \bff{u}^{n+1} \times \bff{H}^{n+1}}{\bff{\xi}^{n+1}} \big|
\norm{\bff{\xi}^{n+1}}{\bb{L}^2}^2
\\
\epsilon \norm{\nabla \bff{\xi}^{n+1}}{\bb{L}^2}^2
\epsilon \norm{\bff{\theta}^n}{\bb{L}^2}^2,
\\
\label{equ:cross dt theta euler}
\nonumber
\big| \inpro{\bff{u}_h^n \times \bff{H}_h^{n+1}- \bff{u}^{n+1} \times \bff{H}^{n+1}}{\delta \bff{\theta}^{n+1}} \big|
h^{2(r+1)} + k^2
\big( \norm{\bff{\theta}^n}{\bb{L}^2}^2 + \epsilon \norm{\nabla \bff{\theta}^n}{\bb{L}^2}^2 \big) \norm{\bff{H}_h^{n+1}}{\bb{H}^1}^2
\\
\norm{\bff{\xi}^{n+1}}{\bb{L}^2}^2
\epsilon \norm{\delta \bff{\theta}^{n+1}}{\bb{L}^2}^2.
\end{align}

Moreover, if the triangulation $\mathcal{T}_h$ is globally quasi-uniform, then for any $\bff{\zeta} \in \bb{V}_h$,

\begin{align}\label{equ:cross est euler}
\left| \inpro{\bff{u}_h^{n} \times \bff{H}_h^{n+1} - \bff{u}^{n+1} \times \bff{H}^{n+1}}{\bff{\zeta}} \right|
h^{2(r+1)} + k^2+ \norm{\bff{\xi}^{n+1}}{\bb{L}^2}^2
+ \epsilon \norm{\bff{\theta}^n}{\bb{L}^2}^2
+ \epsilon \norm{\bff{\zeta}}{\bb{L}^2}^2
\\
\label{equ:cross nab est euler}
\big|\inpro{\nabla\Pi_h\big(\bff{u}_h^n \times \bff{H}_h^{n+1}- \bff{u}^{n+1} \times \bff{H}^{n+1}\big)}{\bff{\zeta}} \big|
h^{2r} + k^2
\norm{\bff{\theta}^n}{\bb{H}^1}^2
\norm{\bff{\xi}^{n+1}}{\bb{H}^1}^2
\epsilon \norm{\bff{\zeta}}{\bb{L}^2}^2.
\end{align}

We write

\begin{align}\label{equ:cross un Hn}
\bff{u}_h^{n} \times \bff{H}_h^{n+1} - \bff{u}^{n+1} \times \bff{H}^{n+1}
\bff{u}_h^{n} \times \left(\bff{\xi}^{n+1}+\bff{\eta}^{n+1}\right)
\left(\bff{\theta}^{n} + \bff{\rho}^{n} - k\cdot \delta \bff{u}^{n+1}\right) \times \bff{H}^{n+1}.
\end{align}

The proof of (<ref>) then follows by arguments similar to that in (<ref>) (noting (<ref>)), without assuming $\mathcal{T}_h$ is globally quasi-uniform. Similarly, the proof of (<ref>) follows that in (<ref>).

Next, following (<ref>), we write

\begin{align*}
\bff{u}_h^{n} \times \bff{H}_h^{n+1} - \bff{u}^{n+1} \times \bff{H}^{n+1}
\bff{\theta}^n \times \bff{H}_h^{n+1}
R_h\bff{u}^n \times (\bff{\xi}^{n+1}+\bff{\eta}^{n+1})
(\bff{\rho}^n- k \cdot \delta \bff{u}^{n+1}) \times \bff{H}^{n+1}.
\end{align*}

The proof of (<ref>) then follows along the line of (<ref>) with obvious modifications.

Now, suppose the triangulation is globally quasi-uniform (and thus (<ref>) holds in this case). Noting (<ref>), by Hölder's inequality and assumptions (<ref>) on the exact solution, we have

\begin{align*}
\nonumber
&\norm{\bff{u}_h^{n} \times \bff{H}_h^{n+1} - \bff{u}^{n+1} \times \bff{H}^{n+1}}{\bb{L}^2}
\\
\nonumber
\norm{\bff{u}_h^{n}}{\bb{L}^\infty} \norm{\bff{\xi}^{n+1}+\bff{\eta}^{n+1}}{\bb{L}^2}
\left( \norm{\bff{\theta}^{n} + \bff{\rho}^{n}}{\bb{L}^2} + k \norm{\delta \bff{u}^{n+1}}{\bb{L}^2} \right) \norm{\bff{H}^{n+1}}{\bb{L}^\infty}
\\
h^{r+1} + \norm{\bff{\xi}^{n+1}}{\bb{L}^2}
+ \norm{\bff{\theta}^n}{\bb{L}^2} + k,
\end{align*}

where in the last step we used (<ref>) and (<ref>). Estimate (<ref>) follows by Young's inequality.
Similarly, noting (<ref>), we obtain (<ref>) following the proof of (<ref>). This completes the proof of the lemma.

Let $\epsilon>0$ be arbitrary. The following inequalities hold:

\begin{align}\label{equ:uh2uh theta disc euler}
\left| \inpro{\abs{\bff{u}_h^{n+1}}^2 \bff{u}_h^{n+1}  - \abs{\bff{u}^{n+1}}^2 \bff{u}^{n+1}}{\bff{\theta}^{n+1}} \right|
h^{2(r+1)} + \norm{\bff{\theta}^{n+1}}{\bb{L}^2}^2
+ \epsilon \norm{\nabla \bff{\theta}^{n+1}}{\bb{L}^2}^2,
\\
\nonumber
\label{equ:uh2uh xi disc euler}
\left| \inpro{\abs{\bff{u}_h^{n+1}}^2 \bff{u}_h^{n+1}  - \abs{\bff{u}^{n+1}}^2 \bff{u}^{n+1}}{\bff{\xi}^{n+1}} \right|
+ \left( 1+\norm{\bff{\xi}^{n+1}}{\bb{H}^1}^2 \right) \norm{\bff{\theta}^{n+1}}{\bb{L}^2}^2
\\
+ \epsilon \norm{\nabla \bff{\theta}^{n+1}}{\bb{L}^2}^2
+ \epsilon \norm{\bff{\xi}^{n+1}}{\bb{L}^2}^2,
\\
\label{equ:uh2uh dt theta euler}
\nonumber
\big| \inpro{\abs{\bff{u}_h^{n+1}}^2 \bff{u}_h^{n+1}  - \abs{\bff{u}^{n+1}}^2 \bff{u}^{n+1}}{\delta \bff{\theta}^{n+1}} \big|
\big( 1+\norm{\delta \bff{\theta}^{n+1}}{\bb{L}^2}^2 \big) \norm{\bff{\theta}^{n+1}}{\bb{H}^1}^2
\\
\epsilon \norm{\delta \bff{\theta}^{n+1}}{\bb{L}^2}^2.
\end{align}

Suppose now that the triangulation $\mathcal{T}_h$ is globally quasi-uniform. Then for any $\bff{\zeta}\in \bb{V}_h$,

\begin{align}\label{equ:psi xi disc euler}
\left| \inpro{\abs{\bff{u}_h^{n+1}}^2 \bff{u}_h^{n+1}  - \abs{\bff{u}^{n+1}}^2 \bff{u}^{n+1}}{\bff{\zeta}} \right|
h^{2(r+1)}  + \norm{\bff{\theta}^{n+1}}{\bb{L}^2}^2
+ \epsilon \norm{\bff{\zeta}}{\bb{L}^2}^2,
\\
\label{equ:psi xi Delta euler}
\left| \inpro{\abs{\bff{u}_h^{n+1}}^2 \bff{u}_h^{n+1}  - \abs{\bff{u}^{n+1}}^2 \bff{u}^{n+1}}{\Delta_h \bff{\zeta}} \right|
h^{2r} + \norm{\bff{\theta}^{n+1}}{\bb{H}^1}^2
\epsilon \norm{\nabla \bff{\zeta}}{\bb{L}^2}^2.
\end{align}

First, we write

\begin{align}\label{equ:psi u2u euler}
\nonumber
&\abs{\bff{u}_h^{n+1}}^2 \bff{u}_h^{n+1}  - \abs{\bff{u}^{n+1}}^2 \bff{u}^{n+1}
\\
\abs{\bff{u}_h^{n+1}}^2  \left(\bff{\theta}^{n+1} + \bff{\rho}^{n+1}\right)
\left(\bff{\theta}^{n+1} + \bff{\rho}^{n+1}\right) \cdot \left( \bff{u}_h^{n+1}+ \bff{u}^{n+1} \right)  \bff{u}^{n+1}.
\end{align}

The proof of (<ref>) then follows the same argument as that of (<ref>).

Next, noting $\bff{u}_h^{n+1}=\bff{\theta}^{n+1}+ R_h \bff{u}^{n+1}$ and following (<ref>), we write

\begin{align}\label{equ:uhn2 minus un2}
\nonumber
\abs{\bff{u}_h^{n+1}}^2 \bff{u}_h^{n+1}  - \abs{\bff{u}^{n+1}}^2 \bff{u}^{n+1}
\big(\bff{\theta}^{n+1} \cdot (\bff{\theta}^{n+1}+ 2 R_h \bff{u}^{n+1}) \big) \bff{u}_h^{n+1}
\abs{R_h \bff{u}^{n+1}}^2 \big(\bff{\theta}^{n+1}+\bff{\rho}^{n+1}\big)
\\
\big( \bff{\rho}^{n+1} \cdot (R_h \bff{u}^{n+1} + \bff{u}^{n+1}) \big) \bff{u}^{n+1}.
\end{align}

Applying the same argument as in the proof of (<ref>), we then obtain (<ref>). The proof of (<ref>) follows along the line of (<ref>) (by replacing $\bff{\theta}$, $\bff{\rho}$, $\bff{u}$, and $\partial_t \bff{\theta}$ with $\bff{\theta}^{n+1}$, $\bff{\rho}^{n+1}$, $\bff{u}^{n+1}$, and $\delta \bff{\theta}^{n+1}$ respectively).

Next, if the triangulation $\mathcal{T}_h$ is quasi-uniform, then (<ref>) holds. As such, we can bound the $\bb{L}^\infty$ norms of $\bff{u}_h^n$ and $\bff{u}_h^{n+1}$ appearing in (<ref>) (uniformly in $n$ and $h$), giving

\begin{align*}
\norm{\abs{\bff{u}_h^{n+1}}^2 \bff{u}_h^{n+1}  - \abs{\bff{u}^{n+1}}^2 \bff{u}^{n+1}}{\bb{L}^2}
\lesssim
\norm{\bff{\theta}^{n+1}}{\bb{L}^2}
+ \norm{\bff{\rho}^{n+1}}{\bb{L}^2}.
\end{align*}

The estimate (<ref>) then follows by Young's inequality.

It remains to prove (<ref>). Let $S$ be the right-hand side of (<ref>). Then we have

\begin{align}\label{equ:inpro S Delta euler}
\inpro{S}{\Delta_h \bff{\zeta}}
\inpro{\nabla \Pi_h S}{\nabla \bff{\zeta}}.
\end{align}

We will proceed by estimating $\norm{\nabla S}{\bb{L}^2}$ using the expression (<ref>). By the product rule for gradient,

\begin{align*}
\nabla S
2\left(\bff{u}_h^{n+1}\cdot \nabla \bff{u}_h^{n+1}\right) \left(\bff{\theta}^{n+1}+\bff{\rho}^{n+1}\right)
\abs{\bff{u}_h^{n+1}}^2 \left(\nabla \bff{\theta}^{n+1}+\nabla \bff{\rho}^{n+1}\right)
\\
\left(\nabla \bff{\theta}^{n+1}+\nabla \bff{\rho}^{n+1}\right) \cdot \left(\bff{u}_h^{n+1}+\bff{u}^{n+1}\right) \bff{u}^{n+1}
\left(\bff{\theta}^{n+1}+\bff{\rho}^{n+1}\right)\cdot \left(\nabla \bff{u}_h^{n+1}+\nabla \bff{u}^{n+1}\right) \bff{u}^{n+1}
\\
\left(\bff{\theta}^{n+1}+\bff{\rho}^{n+1}\right) \cdot \left(\bff{u}_h^{n+1}+\bff{u}^{n+1}\right) \nabla \bff{u}^{n+1}.
\end{align*}

Therefore, by Hölder's inequalities and Sobolev embedding $\bb{H}^1\hookrightarrow \bb{L}^6$ (and noting (<ref>)), it is straightforward to see that

\begin{align*}
\norm{\nabla S}{\bb{L}^2} \lesssim \norm{\bff{\theta}^{n+1}}{\bb{H}^1} + \norm{\bff{\rho}^{n+1}}{\bb{H}^1}
\lesssim
\norm{\bff{\theta}^{n+1}}{\bb{H}^1} + h^r.
\end{align*}

Inequality (<ref>) then follows by applying Young's inequality and (<ref>) to (<ref>).

Let $\epsilon>0$ be arbitrary.

\begin{align}\label{equ:psi psi dt un euler}
\nonumber
&\left| \inpro{\frac{\abs{\bff{u}_h^{n+1}}^2 \bff{u}_h^{n+1}- \abs{\bff{u}_h^n}^2 \bff{u}_h^n}{k} - \partial_t \left(\abs{\bff{u}^{n+1}}^2 \bff{u}^{n+1}\right)}{\bff{\xi}^{n+1}} \right|
\\
\nonumber
\lesssim
h^{2(r+1)} + k^2
\left(1+ \norm{\bff{\xi}^{n+1}}{\bb{H}^1}^2 + \norm{\bff{\xi}^n}{\bb{H}^1}^2\right) \left(\norm{\bff{\theta}^{n+1}}{\bb{H}^1}^2 + \norm{\bff{\theta}^n}{\bb{H}^1}^2 \right)
\\
\norm{\bff{\xi}^n}{\bb{L}^2}^2
\epsilon \norm{\bff{\xi}^{n+1}}{\bb{L}^2}^2
\epsilon \norm{\delta \bff{\theta}^{n+1}}{\bb{L}^2}^2.
\end{align}

Moreover, if the triangulation $\mathcal{T}_h$ is globally quasi-uniform, then for any $\bff{\zeta}\in \bb{V}_h$,

\begin{align}\label{equ:dt psi zeta euler}
\nonumber
\left| \inpro{\frac{\abs{\bff{u}_h^{n+1}}^2 \bff{u}_h^{n+1}- \abs{\bff{u}_h^n}^2 \bff{u}_h^n}{k} - \partial_t \left(\abs{\bff{u}^{n+1}}^2 \bff{u}^{n+1}\right)}{\bff{\zeta}} \right|
\norm{\bff{\theta}^n}{\bb{L}^2}^2
\norm{\bff{\theta}^{n+1}}{\bb{L}^2}^2
\\
\norm{\delta \bff{\theta}^{n+1}}{\bb{L}^2}^2
\epsilon \norm{\bff{\zeta}}{\bb{H}^1}^2.
\end{align}

After some tedious algebra, we can write the first component in the inner product on the left-hand side as

\begin{align}\label{equ:E id}
\nonumber
\frac{\abs{\bff{u}_h^{n+1}}^2 \bff{u}_h^{n+1}- \abs{\bff{u}_h^n}^2 \bff{u}_h^n}{k} - \partial_t \left(\abs{\bff{u}^{n+1}}^2 \bff{u}^{n+1}\right)
\\
\nonumber
\abs{\bff{\theta}^{n+1}+ R_h \bff{u}^{n+1}}^2 \left(\delta \bff{u}_h^{n+1} - \partial_t \bff{u}^{n+1} \right)
\left[\big(\bff{\theta}^{n+1}+\bff{\rho}^{n+1}\big)\cdot \big(\bff{\theta}^{n+1}+R_h\bff{u}^{n+1}+\bff{u}^{n+1}\big)\right] \partial_t \bff{u}^{n+1}
\\
\nonumber
\left[\left(\delta \bff{u}_h^{n+1}- \partial_t \bff{u}^{n+1}\right)\cdot \left(\bff{\theta}^{n+1}+\bff{\theta}^n+ R_h\bff{u}^{n+1}+R_h\bff{u}^n\right) \right]
\left(\bff{\theta}^n+R_h \bff{u}^n\right)
\\
\nonumber
\left[\partial_t \bff{u}^{n+1} \cdot \left(\bff{u}_h^{n+1}+\bff{u}_h^{n}-2\bff{u}^{n+1}\right)\right] \left(\bff{\theta}^n+R_h \bff{u}^n\right)
2 \left(\bff{u}^{n+1}\cdot \partial_t \bff{u}^{n+1}\right) \left(\bff{u}_h^n- \bff{u}^{n+1}\right)
\\
&=: E_1+E_2+E_3+E_4+E_5.
\end{align}

We want to obtain a bound for $\left|\inpro{E_i}{\bff{\xi}^{n+1}}\right|$ for $i=1,2,\cdots,5$. First, we note the frequently used inequality

\begin{align}\label{equ:D est}
\nonumber
\norm{\mathcal{D}}{\bb{L}^2}
\norm{\delta \bff{u}_h^{n+1} - \partial_t \bff{u}^{n+1}}{\bb{L}^2}
\norm{\delta \bff{\theta}^{n+1}+\delta \bff{\rho}^{n+1}+\delta \bff{u}^{n+1}-\partial_t \bff{u}^{n+1}}{\bb{L}^2}
\\
\norm{\delta \bff{\theta}^{n+1}}{\bb{L}^2} + h^{r+1} + k,
\end{align}

where in the last step (<ref>) and (<ref>) were used.
Now, we begin to estimate the term with $E_1$ by applying Hölder's and Young's inequalities (noting (<ref>) and (<ref>)), giving

\begin{align}\label{equ:E1 euler}
\nonumber
\left|\inpro{E_1}{\bff{\xi}^{n+1}}\right|
\norm{\bff{\theta}^{n+1}}{\bb{L}^6}^2 \norm{\mathcal{D}}{\bb{L}^2} \norm{\bff{\xi}^{n+1}}{\bb{L}^6}
\norm{R_h \bff{u}^{n+1}}{\bb{L}^\infty}^2 \norm{\mathcal{D}}{\bb{L}^2} \norm{\bff{\xi}^{n+1}}{\bb{L}^2}
\\
\nonumber
\norm{\bff{\xi}^{n+1}}{\bb{H}^1}^2 \norm{\bff{\theta}^{n+1}}{\bb{H}^1}^2
\epsilon \norm{\mathcal{D}}{\bb{L}^2}^2
\norm{\bff{\xi}^{n+1}}{\bb{L}^2}^2
\\
\left(1+ \norm{\bff{\xi}^{n+1}}{\bb{H}^1}^2\right) \norm{\bff{\theta}^{n+1}}{\bb{H}^1}^2
\epsilon \norm{\delta \bff{\theta}^{n+1}}{\bb{L}^2}^2.
\end{align}

Similarly, for the next term,

\begin{align}\label{equ:E2 euler}
\nonumber
\left|\inpro{E_2}{\bff{\xi}^{n+1}}\right|
\norm{\bff{\theta}^{n+1}+\bff{\rho}^{n+1}}{\bb{L}^2}
\norm{\bff{\theta}^{n+1}}{\bb{L}^6}
\norm{\partial_t \bff{u}^{n+1}}{\bb{L}^6}
\norm{\bff{\xi}^{n+1}}{\bb{L}^6}
\\
\nonumber
\norm{\bff{\theta}^{n+1}+\bff{\rho}^{n+1}}{\bb{L}^2}
\norm{R_h \bff{u}^{n+1}+\bff{u}^{n+1}}{\bb{L}^\infty}
\norm{\partial_t \bff{u}^{n+1}}{\bb{L}^\infty}
\norm{\bff{\xi}^{n+1}}{\bb{L}^2}
\\
\left(1+ \norm{\bff{\xi}^{n+1}}{\bb{H}^1}^2\right) \norm{\bff{\theta}^{n+1}}{\bb{H}^1}^2
\epsilon \norm{\bff{\xi}^{n+1}}{\bb{L}^2}^2.
\end{align}

For the term with $E_3$, by similar argument we obtain

\begin{align}\label{equ:E3 euler}
\nonumber
\left|\inpro{E_3}{\bff{\xi}^{n+1}}\right|
\norm{\mathcal{D}}{\bb{L}^2}
\norm{\bff{\theta}^{n+1}+\bff{\theta}^n}{\bb{L}^6}
\norm{\bff{\theta}^n}{\bb{L}^6}
\norm{\bff{\xi}^{n+1}}{\bb{L}^6}
\norm{\mathcal{D}}{\bb{L}^2}
\norm{\bff{\theta}^{n+1}+\bff{\theta}^n}{\bb{L}^6}
\norm{R_h \bff{u}^n}{\bb{L}^6}
\norm{\bff{\xi}^{n+1}}{\bb{L}^6}
\\
\nonumber
\norm{\mathcal{D}}{\bb{L}^2}
\norm{R_h \bff{u}^{n+1}+R_h\bff{u}^n}{\bb{L}^6}
\norm{\bff{\theta}^n}{\bb{L}^6}
\norm{\bff{\xi}^{n+1}}{\bb{L}^6}
\\
\nonumber
\norm{\mathcal{D}}{\bb{L}^2}
\norm{R_h \bff{u}^{n+1}+R_h\bff{u}^n}{\bb{L}^\infty}
\norm{R_h \bff{u}^n}{\bb{L}^\infty}
\norm{\bff{\xi}^{n+1}}{\bb{L}^2}
\\
h^{2(r+1)} + k^2
\left(1+ \norm{\bff{\xi}^{n+1}}{\bb{H}^1}^2\right) \left(\norm{\bff{\theta}^{n+1}}{\bb{H}^1}^2 + \norm{\bff{\theta}^n}{\bb{H}^1}^2 \right)
\epsilon \norm{\bff{\xi}^{n+1}}{\bb{L}^2}^2.
\end{align}

By the same argument, we also have

\begin{align}\label{equ:E4 euler}
\left|\inpro{E_4}{\bff{\xi}^{n+1}}\right|
h^{2(r+1)} + k^2
\left(1+ \norm{\bff{\xi}^{n+1}}{\bb{H}^1}^2 + \norm{\bff{\xi}^n}{\bb{H}^1}^2\right) \left(\norm{\bff{\theta}^{n+1}}{\bb{H}^1}^2 + \norm{\bff{\theta}^n}{\bb{H}^1}^2 \right)
\epsilon \norm{\bff{\xi}^{n+1}}{\bb{L}^2}^2
\end{align}

\begin{align}\label{equ:E5 euler}
\left|\inpro{E_5}{\bff{\xi}^{n+1}}\right|
h^{2(r+1)} + k^2
\norm{\bff{\theta}^n}{\bb{L}^2}^2
\norm{\bff{\xi}^n}{\bb{L}^2}^2
\epsilon\norm{\bff{\xi}^{n+1}}{\bb{L}^2}^2.
\end{align}

Altogether, (<ref>), (<ref>), (<ref>), (<ref>), and (<ref>) yield (<ref>).

Next, with $E$ as defined in (<ref>), we have the identity

\begin{align*}
\abs{\bff{u}_h^{n+1}}^2 \left(\delta \bff{u}_h^{n+1}-\partial_t \bff{u}^{n+1}\right)
\left[ \big(\bff{\theta}^n+\bff{\rho}^n\big)\cdot \big(\bff{u}_h^{n+1}+\bff{u}^{n+1}\big)\right] \partial_t \bff{u}^{n+1}
\\
\left[ \big(\delta \bff{u}_h^{n+1}-\partial_t \bff{u}^{n+1}\big)\cdot \big(\bff{u}_h^{n+1}+\bff{u}_h^n\big)\right] \bff{u}_h^n
\\
\left[\partial_t \bff{u}^{n+1} \cdot \big(\bff{u}_h^{n+1}+\bff{u}_h^n-2\bff{u}^n\big) \right] \bff{u}_h^n
2 \left(\bff{u}^{n+1}\cdot \partial_t \bff{u}^n\right) \left(\bff{u}_h^n-\bff{u}^{n+1}\right).
\end{align*}

Therefore using the assumption (<ref>), applying Young's inequality, (<ref>), and Taylor's theorem (noting that (<ref>) can be applied by the quasi-uniformity of $\mathcal{T}_h$), we can follow the proof of (<ref>) to obtain the estimate (<ref>).

We also have the following superconvergence estimates on $\bff{\theta}^n$ and $\bff{\xi}^n$ analogous to Proposition <ref>, which is an essential step in the proof of the main theorem.

Suppose $\bff{u}$ and $\bff{H}$ satisfy (<ref>). Then for $h, k>0$ and $n\in \{1,2,\ldots,\lfloor T/k \rfloor\}$,

\begin{align}\label{equ:theta n L2 euler}
\norm{\bff{\theta}^n}{\bb{L}^2}^2
k \sum_{m=0}^{n-1} \norm{\nabla \bff{\theta}^{m+1}}{\bb{L}^2}^2
k \sum_{m=0}^{n-1} \norm{\bff{\xi}^{m+1}}{\bb{L}^2}^2
C \big(h^{2(r+1)}+k^2 \big).
\end{align}

Assume further that one of the following holds:

* $d=1$ or $2$,

* $d=3$ and the triangulation is globally quasi-uniform.

Then for $h, k>0$ and $n\in \{1,2,\ldots,\lfloor T/k \rfloor\}$,

\begin{align}
\label{equ:xi nab theta n L2 euler}
\norm{\bff{\xi}^n}{\bb{L}^2}^2
\norm{\nabla \bff{\theta}^{n}}{\bb{L}^2}^2
k \sum_{m=0}^{n-1} \norm{\delta \bff{\theta}^{m+1}}{\bb{L}^2}^2
k \sum_{m=0}^{n-1} \norm{\nabla \bff{\xi}^{m+1}}{\bb{L}^2}^2
C \big(h^{2(r+1)}+k^4 \big).
\end{align}

Finally, for $d=1$, $2$, $3$, if the triangulation is globally quasi-uniform, then

\begin{align}\label{equ:Delta theta n L2 euler}
\norm{\Delta_h \bff{\theta}^{n}}{\bb{L}^2}^2
\norm{\bff{\theta}^{n}}{\bb{L}^\infty}^2
k\sum_{m=0}^{n-1} \norm{\Delta_h \bff{\xi}^{m+1}}{\bb{L}^2}^2
k\sum_{m=0}^{n-1} \norm{\bff{\xi}^{m+1}}{\bb{L}^\infty}^2
C \big(h^{2(r+1)}+k^2 \big),
\\
\label{equ:nab xi n L2 euler}
\norm{\nabla \bff{\xi}^n}{\bb{L}^2}^2
k \sum_{m=0}^{n-1} \norm{\nabla \Delta_h \bff{\xi}^{m+1}}{\bb{L}^2}^2
C \big(h^{2r}+k^2 \big).
\end{align}

The constant $C$ depends on the coefficients of the equation, $|\mathscr{D}|$, $T$, and $K_0$ (as defined in (<ref>)), but is independent of $n$, $h$, and $k$.

Subtracting (<ref>) from (<ref>) at time step $n+1$, using (<ref>), (<ref>) (and noting the definition of Ritz projection), we obtain for all $\bff{\chi}$, $\bff{\phi}\in \bb{V}_h$,

\begin{align}\label{equ:dt theta disc euler}
\nonumber
\inpro{\delta \bff{\theta}^{n+1} + \delta \bff{\rho}^{n+1} + \delta \bff{u}^{n+1}- \partial_t \bff{u}^{n+1}}{\bff{\chi}}
\lambda_r \inpro{\bff{\xi}^{n+1}+\bff{\eta}^{n+1}}{\bff{\chi}}
\lambda_e \inpro{\nabla \bff{\xi}^{n+1}}{\nabla \bff{\chi}}
\\
\gamma \inpro{\bff{u}_h^n \times \bff{H}_h^{n+1} - \bff{u}^{n+1} \times \bff{H}^{n+1}}{\bff{\chi}}
\end{align}

\begin{align}
\nonumber
\label{equ:xi eta disc euler}
\inpro{\bff{\xi}^{n+1}+\bff{\eta}^{n+1}}{\bff{\phi}}
-\inpro{\nabla \bff{\theta}^{n+1}}{\nabla \bff{\phi}}
\kappa \inpro{\bff{\theta}^n+\bff{\rho}^n+\bff{u}^n-\bff{u}^{n+1}}{\bff{\phi}}
\\
\kappa \inpro{\abs{\bff{u}_h^{n+1}}^2 \bff{u}_h^{n+1} - \abs{\bff{u}^{n+1}}^2 \bff{u}^{n+1}}{\bff{\phi}}
\beta \inpro{\bff{e}\big(\bff{e}\cdot (\bff{\theta}^{n+1}+\bff{\rho}^{n+1})\big)}{\bff{\phi}}.
\end{align}

Taking $\bff{\chi}= \bff{\theta}^{n+1}$ in (<ref>) and $\bff{\phi}=\lambda_r \bff{\theta}^{n+1}$ in (<ref>), then adding the resulting equations, we obtain

\begin{align}\label{equ:dt theta chi disc euler}
\nonumber
&\frac{1}{2k} \left( \norm{\bff{\theta}^{n+1}}{\bb{L}^2}^2 - \norm{\bff{\theta}^n}{\bb{L}^2}^2 \right)
\frac{1}{2k} \norm{\bff{\theta}^{n+1}-\bff{\theta}^n}{\bb{L}^2}^2
\inpro{\delta \bff{\rho}^{n+1} + \delta \bff{u}^{n+1} -\partial_t \bff{u}^{n+1}}{\bff{\theta}^{n+1}}
\lambda_r \norm{\nabla \bff{\theta}^{n+\frac{1}{2}}}{\bb{L}^2}^2
\\
\nonumber
\kappa \lambda_r \inpro{\bff{\theta}^n}{\bff{\theta}^{n+1}}
\kappa\lambda_r \inpro{\bff{\rho}^{n}}{\bff{\theta}^{n+1}}
\kappa\lambda_r \inpro{\bff{u}^n-\bff{u}^{n+1}}{\bff{\theta}^{n+1}}
\lambda_e \inpro{\nabla \bff{\xi}^{n+1}}{\nabla \bff{\theta}^{n+1}}
\\