# Analysis and numerical approximation of energy-variational solutions to the
Ericksen–Leslie equations

Robert Lasarzik, Maximilian E.V. Reiter Weierstrass Institute, Mohrenstr. 39,
10117 Berlin, Germany<EMAIL_ADDRESS>University of
Berlin, Straße des 17. Juni 135, 10623 Berlin, Germany<EMAIL_ADDRESS>
###### Abstract

We define the concept of energy-variational solutions for the Ericksen–Leslie
equations in three spatial dimensions. This solution concept is finer than
dissipative solutions and satisfies the weak-strong uniqueness property. For a
certain choice of the regularity weight, the existence of energy-variational
solutions implies the existence of measure-valued solutions and for a
different choice, we construct an energy-variational solution with the help of
an implementable, structure-inheriting space-time discretization.
Computational studies are performed in order to provide some evidence of the
applicability of the proposed algorithm.
MSC(2020): 35A35, 35Q35, 65M60, 76A15
Keywords: Existence, liquid crystal, Ericksen–Leslie, energy-variational
solutions, numerical approximation, unit-norm constraint, mass-lumping, finite
element method

###### Contents

  1. 1 Introduction
  2. 2 Preliminaries, main results, and continuous system
    1. 2.1 Main results
    2. 2.2 Proofs for the continuous system
    3. 2.3 Auxiliary Lemmata
    4. 2.4 Finite element spaces
    5. 2.5 Interpolation and mass-lumping
  3. 3 Discrete system
    1. 3.1 A priori estimates
    2. 3.2 Converging subsequences
    3. 3.3 Unit-norm restriction and convergence for the director equations
    4. 3.4 Convergence to the energy-variational formulation
  4. 4 Computational studies
    1. 4.1 A smooth example in two spatial dimensions
    2. 4.2 Annihilation of two defects
    3. 4.3 Annihilation of two defects in a rotating flow

##  1 Introduction

Liquid crystals comprise the structural properties of crystals within a fluid.
The fluid flow of the nematic phase of liquid crystals can be described by the
Ericksen–Leslie system. In this model, the material behaves like a liquid,
i.e., no positional order is present, but the molecules exhibit a long-range
self-alignment along a direction (see Figure 1).

Figure 1: Alignment of the rod-like molecules in an isotropic liquid, the
nematic phase of a liquid crystal, and a solid, [24, Fig. 1.1]

In this way, liquid crystals sustain anisotropic dynamics such as the
polarization of light or the transfer of heat, but at the same time offer the
physical flexibility of a fluid, which makes these materials interesting for
engineering and sciences [23].

Ericksen [16] and Leslie [31] derived the Ericksen-Leslie equations during
their development of an instationary theory of liquid crystals in the 1960s.
Let ${\boldsymbol{v}}:\overline{\Omega}\times[0,T]\to\mathbb{R}^{3}$ denote
the velocity of the fluid,
${\boldsymbol{p}}:\overline{\Omega}\times[0,T]\to\mathbb{R}$ its pressure and
${\boldsymbol{d}}:\overline{\Omega}\times[0,T]\to\mathbb{R}^{3}$ the director.
We consider the system governed by the equations

$\displaystyle\begin{split}\partial_{t}{\boldsymbol{v}}+({\boldsymbol{v}}\cdot\nabla){\boldsymbol{v}}+\nabla{\boldsymbol{p}}+\operatorname{\nabla\cdot}\left((\nabla{\boldsymbol{d}})^{T}\nabla{\boldsymbol{d}}\right)-\nabla\cdot{\boldsymbol{T}}^{L}&=0,\end{split}$
(1a) $\displaystyle\begin{split}\nabla\cdot{\boldsymbol{v}}&=0,\end{split}$
(1b)
$\displaystyle\begin{split}\partial_{t}{\boldsymbol{d}}+({\boldsymbol{v}}\cdot\nabla){\boldsymbol{d}}-(\nabla{\boldsymbol{v}})_{skw}{\boldsymbol{d}}+\left(I-{\boldsymbol{d}}\otimes{\boldsymbol{d}}\right)\left(\lambda(\nabla{\boldsymbol{v}})_{sym}{\boldsymbol{d}}-\Delta{\boldsymbol{d}}\right)&=0,\end{split}$
(1c) $\displaystyle\begin{split}|{\boldsymbol{d}}|&=1,\end{split}$ (1d)

where we employ the initial conditions

$\displaystyle{\boldsymbol{v}}(0)={\boldsymbol{v}}_{0},\quad{\boldsymbol{d}}(0)={\boldsymbol{d}}_{0}\text{
with }\left\lvert{\boldsymbol{d}}_{0}\right\rvert=1\text{ in }\Omega$

and boundary conditions

$\displaystyle{\boldsymbol{v}}=0,\quad{\boldsymbol{d}}={\boldsymbol{d}}_{\Gamma}\text{
with }\left\lvert{\boldsymbol{d}}_{\Gamma}\right\rvert=1\text{ on
}\partial\Omega.$

The Leslie stress tensor $T^{L}$ is defined by

${\boldsymbol{T}}^{L}\coloneqq{\boldsymbol{T}}^{D}+\lambda[{\boldsymbol{d}}\otimes[{\boldsymbol{d}}]_{x}^{T}[{\boldsymbol{d}}]_{x}\cdot\Delta{\boldsymbol{d}}]_{sym}+[{\boldsymbol{d}}\otimes\Delta{\boldsymbol{d}}]_{skw},$

where $T^{D}$ collects the dissipative terms of the Leslie stress tensor, i.e.

$\displaystyle{\boldsymbol{T}}^{D}\coloneqq$
$\displaystyle(\mu_{1}+\lambda^{2})({\boldsymbol{d}}\cdot(\nabla{\boldsymbol{v}})_{sym}{\boldsymbol{d}})({\boldsymbol{d}}\otimes{\boldsymbol{d}})+\mu_{4}(\nabla{\boldsymbol{v}})_{sym}$
$\displaystyle+(\mu_{5}+\mu_{6}-\lambda^{2})({\boldsymbol{d}}\otimes(\nabla{\boldsymbol{v}})_{sym}{\boldsymbol{d}})_{sym}$

with $\mu_{4}>0,\quad\mu_{5}+\mu_{6}-\lambda^{2}\geq
0,\quad\mu_{1}+\lambda^{2}\geq 0$ in order to ensure the dissipative character
of our model.

So far, a vast majority of the mathematical work on the Ericksen–Leslie model
considers a simplified system with a relaxed unit-norm constraint that is only
enforced approximately by adding a Ginzburg–Landau penalization term
${\boldsymbol{f}}_{\epsilon}(d)=\frac{1}{4\epsilon}(\left\lvert{\boldsymbol{d}}\right\rvert^{2}-1)^{2}$
to the free energy potential $\frac{1}{2}|\nabla{\boldsymbol{d}}|^{2}$. With
simplified Leslie-stress tensor, the momentum and director equation (1) are
replaced by

$\displaystyle\begin{aligned}
\hfil\displaystyle\begin{split}\partial_{t}{\boldsymbol{v}}-\frac{1}{Re}\Delta{\boldsymbol{v}}+({\boldsymbol{v}}\cdot\nabla){\boldsymbol{v}}+\nabla{\boldsymbol{p}}+\gamma\nabla\cdot\left((\nabla{\boldsymbol{d}})^{T}\nabla{\boldsymbol{d}}\right)&=0,\quad\operatorname{\nabla\cdot}{\boldsymbol{v}}=0\,,\end{split}\\\\[8.61108pt]
\hfil\displaystyle\begin{split}\partial_{t}{\boldsymbol{d}}+({\boldsymbol{v}}\cdot\nabla){\boldsymbol{d}}-\Delta{\boldsymbol{d}}+{\boldsymbol{f}}_{\epsilon}({\boldsymbol{d}})&=0,\end{split}\end{aligned}$
(2)

respectively. In the first analysis of this system [33], the authors were able
to prove the existence of weak solutions to (2). A rather general model
including the full Leslie stress tensor is considered by [10], where again a
Ginzburg–Landau penalization approach is introduced to replace the unit-norm
constraint of the director. In this setting the authors prove that weak
solutions exist and a blow-up criterion for local strong solutions. In [15]
the existence of weak solutions is generalized to a larger class of free
energy functions. An overview of the analytical results regarding the
Ericksen–Leslie equations and its connection to other models for liquid
crystals can be found in [14]. In two spacial dimensions, for the limiting
system of (2), where $f_{\varepsilon}({\boldsymbol{d}})$ is replaced by
$-|\nabla{\boldsymbol{d}}|^{2}$, it is known that a unique weak solution
exists that is smooth except for finitely many points in time [34].

For a general model of the Ericksen–Leslie equations equipped with the
naturally arising anisotropic Oseen–Frank energy, the concept of dissipative
solutions is applied in [26], which are shown to be the local average of
measure-valued solutions [27] and inherit their weak-strong uniqueness [28].
In comparison to measure-valued solutions, dissipative solutions have the
advantage that they have less degrees of freedom and can be approximated by
numerical schemes. Nevertheless, they form only a subset of measure-valued
solutions and are thus not as precise. In the work at hand, we introduce
energy-variational solutions (cf. [30]), which have one degree of freedom more
than dissipative solutions and can be argued to be as fine as measure-valued
solutions, but as we will show, they can also be approximated by numerical
schemes and have additional advantages. This solution concept is useful where
one might either not be able to derive weak solutions for physically relevant
models or where weak solutions admit unphysical behaviour, like unphysical
non-uniqueness [22]. The solution concept was first introduced in [29]. We
will prove that this concept is finer than the so-called dissipative solutions
introduced by Lions [36] for the Euler equations and also applied to the
Ericksen–Leslie equations by one of the authors [26]. Additionally, in certain
scenarios this concept is finer than the concept of measure-valued solutions.
Energy-variational solutions do not only fulfill the standard weak-strong
uniqueness property of generalized solution concepts (cf. [28]), but they also
fulfill the semi-flow property such that prolongations and restrictions of
solutions on larger, and smaller time intervals, respectively, are energy-
variational solutions again. In [29] it was also argued that energy-
variational solutions are amenable for different selection criteria. The set
of energy-variational solutions can be seen as a convex, weakly∗-closed
superset of weak solutions but a subset of dissipative solutions and for a
special choice of the regularity weight also a subset of measure-valued
solutions. This may allow to introduce techniques from optimization theory in
order to select the physically relevant solution maximizing the dissipation in
every point-in-time [29].

In this work, we define the concept of energy-variational solutions for the
Ericksen–Leslie equations in three spatial dimensions. This definition has
some freedom, since it depends on the choice of a certain regularity weight
$\mathcal{K}$. For one choice of such a regularity weight, we prove the
equivalence to a certain class of measure-valued solutions, and for another
choice, we construct an energy-variational solution with the help of an
implementable, structure-inheriting space-time discretization based on the
finite element method.

Concerning the numerical approximation of the Ericksen–Leslie equations, a
first study for the simplified system (2) equipped with the Ginzburg–Landau
approximation is conducted in [37]. They combined an implicit Euler scheme in
time with Hermite type finite elements for the director and
$Q2$-$Q1$-Taylor–Hood elements for the velocity and pressure. In their
subsequent work [38], the authors replace the demanding Hermite finite
elements for the director by piecewise quadratic functions. Even a relaxation
from $C^{1}$ finite elements to $C^{0}$ finite elements is realized in [35]. A
different approach is proposed in [3], where the Ginzburg–Landau approximation
of the unit-norm constraint is interpreted as a saddle-point structure. For
the simplified and penalized Ericksen–Leslie system also decoupling techniques
and mixed methods are examined in [18, 9]. In [5] two numerical schemes for
the simplified model are proposed. The first one uses the Ginzburg–Landau
approximation for the unit-norm constraint and the second proposed scheme does
not depend on a regularization parameter and fulfills the unit-norm constraint
in the limit. However both schemes do not fulfill the unit-norm constraint at
the discrete level exactly. We therefore use an approach from [4] by
implementing a midpoint rule at the finite element level which solves the
sphere constraint exactly at every node of the mesh. But we have to refine
this approach by introducing a special projection (see Remark Remark), which
will allow to identify the limit as an energy-variational solution. The
proposed scheme is the first numerical scheme implementing the main properties
of the continuous system including the algebraic norm restriction at every
node of the mesh such that the approximate solutions converge to an energy-
variational solution fulfilling the physically relevant semi-flow property. We
think that the concept of energy-variational solutions is a strong tool for
identifying limits of solutions to numerical schemes and can also serve as
such in other related models.

After providing the considered system and an overview over the existing
literature, we introduce the necessary notation in Section 2, we provide the
definitions and the main results in Subsection 2.1. The weak-strong uniqueness
proof and the relation to measure-valued solutions in considered in Subsection
2.2, whereas the preliminaries are provided in the following subsections on
auxiliary lemmata 2.3, finite elements 2.4, and interpolation 2.5. The
discrete system is introduced and analysed in Section 3, we provide the
necessary a priori estimates 3.1, extract converging subsequences 3.2, and
prove convergence to the director equations 3.3 and the energy-variational
inequality 3.4. Finally, some computational studies are presented in Section 4
in order to show some evidence of the applicability of the proposed algorithm.

##  2 Preliminaries, main results, and continuous system

We denote the space of smooth solenoidal functions with compact support by
$\mathcal{C}_{c,\sigma}^{\infty}(\Omega;\operatorname{\mathbb{R}}^{3})$. By
$L^{p}_{\sigma}(\Omega)$, $H^{1}_{0,\sigma}(\Omega)$, and
$W^{1,p}_{0,\sigma}(\Omega)$, we denote the closure of
$\mathcal{C}_{c,\sigma}^{\infty}(\Omega;\operatorname{\mathbb{R}}^{3})$ with
respect to the norm of $L^{p}(\Omega)$, $H^{1}(\Omega)$, and
$W^{1,p}(\Omega)$, respectively. By $H^{1}(\Omega;\mathbb{S}^{2})$, we denote
the functions ${\boldsymbol{d}}\in H^{1}(\Omega)$ such that
$|{\boldsymbol{d}}|=1$ a.e. in $\Omega$. The $L^{2}(\Omega)$ inner product is
thereby denoted by $(.,.)$. The Dual space of a Banach space $V$ is denoted by
$V^{*}$, where the dual pairing is denoted as $\langle\cdot,\cdot\rangle$. In
order to define the cross product, we use the Levi-Cita symbol. The cross
product of two vectors $x,y\in\mathbb{R}^{3}$ is then defined as $(x\times
y)_{i}=\sum_{j,k}\epsilon_{i,j,k}x_{j}y_{k}$, the cross product of a vector
$x\in\mathbb{R}^{3}$ with a matrix $A\in\mathbb{R}^{d\times d}$ as $(x\times
A)_{i,j}=\sum_{l,m}\epsilon_{i,l,m}x_{l}A_{m,j}$. We further will make use of
the matrix notation of the cross product for three dimensions, i.e.
$([x]_{\times})_{ik}=\sum_{j}\epsilon_{i,j,k}x_{j}$. We further introduce the
discrete derivative $\boldsymbol{d}_{t}$ as
$\boldsymbol{d}_{t}f^{j}\coloneqq\frac{f^{j}-f^{j-1}}{k}$ for a constant time-
step size $k>0$ and for a sequence in a normed space $(f_{j})_{j=1,...,n}$,
where we set $\boldsymbol{d}_{t}f^{0}=0$. Usually functions in the continuous
setting are denoted by bold letters in contrast to the discrete functions.

By $\mathbb{R}^{d\times d}$ we denote $d$-dimensional quadratic matrices, by
$\mathbb{R}^{d\times d}_{\text{sym}}$ the symmetric subset, and by
$\mathbb{R}^{d\times d}_{\text{sym},+}$ the symmetric positive semi-definite
matrices. The symmetric and skew-symmetric part of a matrix
$A\in\operatorname{\mathbb{R}}^{d\times d}$ are denoted by $(A)_{\text{sym}}$
and $(A)_{\text{skw}}$, respectively. The positive and negative semi-definite
part of a matrix $A\in\operatorname{\mathbb{R}}^{d\times d}$ are denoted by
$(A)_{+}$ and $(A)_{-}$, respectively. We equip the last set with the usual
spectral norm $|A|_{2}=\max_{i\in\\{1,\ldots,d\\}}\lambda_{i}$, where
$\lambda_{i}$ are the nonnegative eigenvalues of the matrix
$A\in\mathbb{R}^{d\times d}_{\text{sym},+}$. The dual norm of the spectral
norm with respect to the Frobenius product
($A:B:=\sum_{i,j=1}^{d}A_{ji}B_{ij}$ for
$A,B\in\operatorname{\mathbb{R}}^{d\times d}$) is given by
$|A|^{\prime}_{2}=\sum_{i=1}^{d}\lambda_{i}=A:I=\operatorname{tr}(A)$ for a
matrix $A\in\mathbb{R}^{d\times d}_{\text{sym},+}$. The Radon measures taking
values in $\mathbb{R}^{d\times d}_{\text{sym}}$ are denoted by
$\mathcal{M}(\overline{{\Omega}};\mathbb{R}^{d\times d}_{\text{sym}})$, which
may be interpreted as the dual space of the continuous functions,i.e.,
$\mathcal{M}(\overline{{\Omega}};\mathbb{R}^{d\times
d}_{\text{sym}})=(\operatorname{\mathcal{C}}(\overline{{\Omega}};\mathbb{R}^{d\times
d}_{\text{sym}}))^{*}$. Note that an element
$\mu\in\mathcal{M}(\overline{{\Omega}};\mathbb{M}^{d\times d}_{\text{sym},+})$
is a Radon measure taking values in the symmetric matrices such that for any
${\boldsymbol{\xi}}\in\operatorname{\mathbb{R}}^{d}$ the measure
${\boldsymbol{\xi}}\otimes{\boldsymbol{\xi}}:\mu$ is nonnegative. By $I$, we
denote the identity matrix in $\operatorname{\mathbb{R}}^{d\times d}$.

For a given Banach space $\mathbb{X}$, the space
$\operatorname{\mathcal{C}}_{w}([0,T];\mathbb{X})$ denotes the functions on
$[0,T]$ taking values in $\mathbb{X}$ that are continuous with respect to the
weak topology of $\mathbb{X}$. The space
$L^{\infty}_{w^{*}}([0,T];\mathbb{X}^{*})$ is the space of all function on
$[0,T]$ taking values in $\mathbb{X}^{*}$ that are Bochner measurable with
respect to $\mathbb{X}^{*}$ equipped with the weak-stark topology and
essentially bounded. The total variation of a function
$E:[0,T]\operatorname{\rightarrow}\operatorname{\mathbb{R}}$ is given by

$|E|_{\text{TV}([0,T])}=\sup_{0=t_{0}<\ldots<t_{n}=T}\sum_{k=1}^{n}\lvert
E(t_{k-1})-E(t_{k})\rvert\,,$

where the supremum is taken over all finite partitions of the interval
$[0,T]$. We denote the space of all bounded functions of bounded variations on
$[0,T]$ by $\operatorname{{BV}([0,T])}$. Note that the total variation of a
monotone decreasing nonnegative function $E$ only depends on the initial
value, i.e.,

$\displaystyle|E|_{\text{TV}([0,T])}=\sup_{0=t_{0}<\ldots<t_{n}=T}\sum_{k=1}^{N}\lvert
E(t_{k-1})-E(t_{k})\rvert\leq E(0)-E(T)\leq E(0)\,.$

###  2.1 Main results

We define the total energy as

$\displaystyle\mathcal{E}({\boldsymbol{v}},{\boldsymbol{d}})\coloneqq\frac{1}{2}\lVert{\boldsymbol{v}}\rVert_{L^{2}(\Omega)}^{2}+\frac{1}{2}\left\lVert\nabla{\boldsymbol{d}}\right\rVert_{L^{2}\left(\Omega\right)}^{2}\,.$
(3)

The main new idea when defining energy-variational solutions is to introduce
an auxiliary variable $E$ as an upper bound of the total energy
$\mathcal{E}({\boldsymbol{v}},{\boldsymbol{d}})$ and add the difference of
these two variables to a variational inequality of the energy-dissipation
mechanism to close this formulation with respect to the appropriate weak
topologies. The difference between the variable $E$ and the energy
$\mathcal{E}({\boldsymbol{v}},{\boldsymbol{d}})$ can be interpreted as a
measure of the difference between weak and strong convergence of the
approximate solutions. The associated error term represents this difference in
the limit of vanishing discretization parameters.

######  Definition 2.1 (Energy-variational solutions).

We call $({\boldsymbol{v}},{\boldsymbol{d}},E)$ an energy-variational
solution, if

$\displaystyle\begin{split}&{\boldsymbol{v}}\in\operatorname{\mathcal{C}}_{w}([0,T];{L}^{2}_{\sigma}(\Omega))\cap
L^{2}(0,T;H^{1}_{0,\sigma}(\Omega))\,,\quad
E\in\operatorname{{BV}([0,T])}\,,\\\
&{\boldsymbol{d}}\in\operatorname{\mathcal{C}}_{w}([0,T];H^{1}(\Omega;\mathbb{S}^{2}))\cap
H^{1}(0,T;L^{3/2}(\Omega;\operatorname{\mathbb{R}}^{3})\text{ such that }\\\
&{\boldsymbol{d}}\times\Delta{\boldsymbol{d}}\in
L^{2}(0,T;L^{2}(\Omega;\operatorname{\mathbb{R}}^{3}))\,,\end{split}$ (4)

and $E\geq\mathcal{E}({\boldsymbol{v}},{\boldsymbol{d}})$ on $[0,T]$ as well
as $|{\boldsymbol{d}}|=1$ a.e. in $\Omega\times(0,T)$. The term
${\boldsymbol{d}}\times\Delta{\boldsymbol{d}}$ has to be understood as the
weak divergence of ${\boldsymbol{d}}\times\nabla{\boldsymbol{d}}$. The
solution fulfills the energy-variational inequality

$\left[E-\int_{\Omega}{\boldsymbol{v}}\cdot\tilde{{\boldsymbol{v}}}\mathop{}\\!\mathrm{d}x\right]\Big{|}_{s-}^{t}+\int_{s}^{t}\int_{\Omega}{\boldsymbol{v}}\cdot\partial_{t}\tilde{{\boldsymbol{v}}}-({\boldsymbol{v}}\cdot\nabla){\boldsymbol{v}}\cdot\tilde{{\boldsymbol{v}}}\mathop{}\\!\mathrm{d}x\mathop{}\\!\mathrm{d}\tau+\int_{s}^{t}\langle{\boldsymbol{T}}^{E}({\boldsymbol{d}}),\tilde{{\boldsymbol{v}}}\rangle\text{d}\tau\\\
+\int_{s}^{t}\int_{\Omega}\left((\mu_{1}+\lambda^{2})({\boldsymbol{d}}\cdot(\nabla{\boldsymbol{v}})_{sym}{\boldsymbol{d}})({\boldsymbol{d}}\otimes{\boldsymbol{d}})+\mu_{4}(\nabla{\boldsymbol{v}})_{sym}\right):(\nabla{\boldsymbol{v}}-\nabla\tilde{{\boldsymbol{v}}})\mathop{}\\!\mathrm{d}x\mathop{}\\!\mathrm{d}\tau\\\
+\int_{s}^{t}\int_{\Omega}(\mu_{5}+\mu_{6}-\lambda^{2})({\boldsymbol{d}}\otimes(\nabla{\boldsymbol{v}})_{sym}{\boldsymbol{d}})_{sym}:(\nabla{\boldsymbol{v}}-\nabla\tilde{{\boldsymbol{v}}})+|{\boldsymbol{d}}\times\Delta{\boldsymbol{d}}|^{2}\mathop{}\\!\mathrm{d}x\mathop{}\\!\mathrm{d}\tau\\\
+\int_{s}^{t}\int_{\Omega}\left(\lambda[{\boldsymbol{d}}\otimes[{\boldsymbol{d}}]_{x}^{T}[{\boldsymbol{d}}]_{x}\cdot\Delta{\boldsymbol{d}}]_{sym}+[{\boldsymbol{d}}\otimes\Delta{\boldsymbol{d}}]_{skw}\right):\nabla\tilde{{\boldsymbol{v}}}\mathop{}\\!\mathrm{d}x\mathop{}\\!\mathrm{d}\tau\\\
+\int_{s}^{t}\mathcal{K}(\tilde{{\boldsymbol{v}}})\left[\mathcal{E}({\boldsymbol{v}},{\boldsymbol{d}})-E\right]\text{d}\tau\leq
0\,.$ (5)

for all $s$, $t\in(0,T)$ and for all
$\tilde{{\boldsymbol{v}}}\in\operatorname{\mathcal{C}}^{1}([0,T];(H^{1}_{0,\sigma}(\Omega))^{*})\cap
L^{2}(0,T;H^{1}_{0,\sigma}(\Omega))$ such that
$\mathcal{K}(\tilde{{\boldsymbol{v}}})\in L^{1}(0,T)$. The initial values are
fulfilled in a weak sense ${\boldsymbol{v}}(0)={\boldsymbol{v}}_{0}$,
${\boldsymbol{d}}(0)={\boldsymbol{d}}_{0}$ and ${\boldsymbol{d}}$ fulfills the
inhomogeneuous Dirichlet boundary conditions in the sense of the trace
$\operatorname{tr}({\boldsymbol{d}}(t))=\operatorname{tr}({\boldsymbol{d}}_{\Gamma})$
for a.e. $t\in(0,T)$. Additionally, it holds

$\displaystyle\partial_{t}{\boldsymbol{d}}+({\boldsymbol{v}}\cdot\nabla){\boldsymbol{d}}-(\nabla{\boldsymbol{v}})_{\text{skw}}{\boldsymbol{d}}+(I-{\boldsymbol{d}}\otimes{\boldsymbol{d}})\left(\lambda(\nabla{\boldsymbol{v}})_{\text{sym}}{\boldsymbol{d}}-\Delta{\boldsymbol{d}}\right)=0$
(6)

a.e. in $\Omega\times(0,T)$. In this work, we choose
$\mathcal{K}_{1}(\tilde{{\boldsymbol{v}}})=\frac{1}{2}\|\tilde{{\boldsymbol{v}}}\|_{L^{\infty}(\Omega)}^{2}$
and
$\mathcal{K}_{2}(\tilde{{\boldsymbol{v}}})=2\|(\nabla\tilde{{\boldsymbol{v}}})_{\text{sym},-}\|_{\operatorname{\mathcal{C}}(\overline{{\Omega}})}$,
where $c_{L}$ is given in (22), with the associated Ericksen stresses being
defined as

$\langle{\boldsymbol{T}}_{1}^{E}({\boldsymbol{d}}),\tilde{{\boldsymbol{v}}}\rangle=\int_{\Omega}({\boldsymbol{d}}\times(\tilde{{\boldsymbol{v}}}\cdot\nabla){\boldsymbol{d}})\cdot({\boldsymbol{d}}\times\Delta{\boldsymbol{d}})\text{d}x\text{
or
}\langle{\boldsymbol{T}}_{2}^{E}({\boldsymbol{d}}),\tilde{{\boldsymbol{v}}}\rangle=-\int_{\Omega}(\nabla{\boldsymbol{d}}^{T}\nabla{\boldsymbol{d}}):\nabla\tilde{{\boldsymbol{v}}}\text{d}x\,,$

respectively.

######  Remark (Ericksen stress).

The main obstacle when passing to the limit in an approximation of the
Ericksen–Leslie equations is often the so-called Ericksen stress, which can be
seen as a kind of Korteweg stress coupling the elastic deformations of the
director field to the evolution of the velocity field. This term is the main
reason that we have to consider a generalized solution framework different
from usual weak solutions. Via Lemma 2.6 we observe that both formulations
${\boldsymbol{T}}^{E}_{1}$ and ${\boldsymbol{T}}^{E}_{2}$ are formally
equivalent via the reformulation

$\displaystyle-\int_{\Omega}(\tilde{{\boldsymbol{v}}}\cdot\nabla){\boldsymbol{d}}\cdot[{\boldsymbol{d}}]_{x}^{T}$
$\displaystyle[{\boldsymbol{d}}]_{x}\Delta{\boldsymbol{d}}\text{d}x={}-\int_{\Omega}(\tilde{{\boldsymbol{v}}}\cdot\nabla){\boldsymbol{d}}\cdot(\Delta{\boldsymbol{d}}+|\nabla{\boldsymbol{d}}|^{2}{\boldsymbol{d}})\text{d}x$
$\displaystyle={}$
$\displaystyle\int_{\Omega}(\nabla{\boldsymbol{d}}^{T}\nabla{\boldsymbol{d}}):\nabla\tilde{{\boldsymbol{v}}}\text{d}x+\frac{1}{2}\int_{\Omega}(\tilde{{\boldsymbol{v}}}\cdot\nabla)|\nabla{\boldsymbol{d}}|^{2}-|\nabla{\boldsymbol{d}}|^{2}(\tilde{{\boldsymbol{v}}}\cdot\nabla)|{\boldsymbol{d}}|^{2}\text{d}x\,,$

where the second term on the right-hand side vanish formally due to the fact
that $\tilde{{\boldsymbol{v}}}\in H^{1}_{0,\sigma}(\Omega)$ and
$|{\boldsymbol{d}}|=1$. Thus, in the case of a regular solution with
${\boldsymbol{d}}\in L^{2}(0,T;H^{2}(\Omega;\mathbb{S}^{2}))$ and
$E=\mathcal{E}({\boldsymbol{v}},{\boldsymbol{d}})$ a.e. in $(0,T)$ both
solution concepts coincide with usual weak solutions.

######  Remark (Strong continuity of the initial value).

From the monotony of $t\mapsto E(t)$ and the weak continuity of
$t\mapsto({\boldsymbol{v}}(t),{\boldsymbol{d}}(t))$, we infer by the lower
semi-continuity of the energy functional and the inequality
$E(t)\geq\mathcal{E}({\boldsymbol{v}}(t),{\boldsymbol{d}}(t))$ for all
$t\in[0,T]$ that for any sequence $t\searrow 0$ it holds

$\displaystyle\mathcal{E}({\boldsymbol{v}}_{0},{\boldsymbol{d}}_{0})=E(0)\geq\lim_{t\searrow
0}E(t)\geq\liminf_{t\searrow
0}\mathcal{E}({\boldsymbol{v}}(t),{\boldsymbol{d}}(t))\geq\mathcal{E}({\boldsymbol{v}}_{0},{\boldsymbol{d}}_{0})\,.$

This implies that all inequalities are in-fact equations such that we infer
from the uniform convexity of $L^{2}(\Omega)$ that

$\left\\{\begin{matrix}({\boldsymbol{v}}(t),{\boldsymbol{d}}(t))\rightharpoonup&({\boldsymbol{v}}_{0},{\boldsymbol{d}}_{0})\text{
as }t\searrow 0\text{ in }{L}^{2}_{\sigma}(\Omega)\times
L^{2}(\Omega;\operatorname{\mathbb{R}}^{3})\\\
\mathcal{E}({\boldsymbol{v}}(t),{\boldsymbol{d}}(t))\longrightarrow&\mathcal{E}({\boldsymbol{v}}_{0},{\boldsymbol{d}}_{0})\text{
as }t\searrow 0\text{ in
}\operatorname{\mathbb{R}}\end{matrix}\right\\}\Longrightarrow\\\
({\boldsymbol{v}}(t),{\boldsymbol{d}}(t))\longrightarrow({\boldsymbol{v}}_{0},{\boldsymbol{d}}_{0})\text{
as }t\searrow 0\text{ in }{L}^{2}_{\sigma}(\Omega)\times
L^{2}(\Omega;\operatorname{\mathbb{R}}^{3})\,.$

######  Remark (Properties of generalized solutions).

Energy-variational solutions fulfill several properties, which are desirable
for generalized solution concepts. That generalized solutions coincide with
classical solutions, if the latter exists, is the so-called weak-strong
uniqueness property and covered by Theorem 2.5.

Furthermore, energy-variational solutions fulfill the so-called semi-flow
property such that restrictions and concatinations of solutions are solutions
again, which is an important property of a reasonable solution concept (cf.
[30]).

######  Remark (Different definition).

We could also identify
$E(t)=E^{1}(t)+\frac{1}{2}\|{\boldsymbol{v}}(t)\|_{L^{2}(\Omega)}^{2}$ for
a.e. $t\in(0,T)$ such that
$E^{1}(t)\geq\frac{1}{2}\|\nabla{\boldsymbol{d}}(t)\|_{L^{2}(\Omega)}^{2}$ for
a.e. $t\in(0,T)$ in the inequality (5). This can be done due to the strong
convergence of ${\boldsymbol{v}}$ a.e. in $\Omega\times(0,T)$. We would also
replace
$\mathcal{E}({\boldsymbol{v}},{\boldsymbol{d}})-E=\frac{1}{2}\|\nabla{\boldsymbol{d}}\|_{L^{2}(\Omega)}^{2}-E^{1}$.
The disadvantage would be that the adapted inequality (5) is only fulfilled
for a.e. $s\in(0,T)$ and all $t\in(s,T]$.

######  Theorem 2.2.

Let $\Omega\subset\operatorname{\mathbb{R}}^{3}$ be a bounded convex
polyhedral domain. Let
$\bigl{(}{\boldsymbol{v}}_{0},{\boldsymbol{d}}_{0}\bigr{)}\in
H^{1}_{0,\sigma}(\Omega)\times H^{2}(\Omega;\operatorname{\mathbb{R}}^{3})$,
with $|{\boldsymbol{d}}_{0}|=1$ a.e. in $\Omega$ and
${\boldsymbol{d}}_{\Gamma}=\operatorname{tr}({\boldsymbol{d}}_{0})$. Then
there exists an energy-variational solution in the sense of Definition 2.1
with
$\mathcal{K}(\tilde{{\boldsymbol{v}}})=\frac{1}{2}\|\tilde{{\boldsymbol{v}}}\|_{L^{\infty}(\Omega)}^{2}$
such that $E(0)=\mathcal{E}({\boldsymbol{v}}_{0},{\boldsymbol{d}}_{0})$.

We are going to prove the theorem by the convergence of a fully discrete,
implementable, unconditionally solvable scheme in the following sections. The
proposed scheme in Section 3 is a fully implementable numerical scheme for the
Ericksen–Leslie system that fulfills the norm restriction
$|{\boldsymbol{d}}|=1$ for every approximate solution in every node of the
mesh. We even infer that the norm of the approximate solutions converges to
$1$ with the order $1$ (cf. Proposition 3.4). Note that the condition
${\boldsymbol{d}}_{\Gamma}=\operatorname{tr}({\boldsymbol{d}}_{0})$ is a usual
compatibility condition for parabolic problems.

Energy-variational solutions can be seen as a generalized solution concept,
which generalizes the well-known weak solutions, but is finer than the usual
concept of so-called measure-valued solutions [27].

######  Definition 2.3.

We call the tuple $({\boldsymbol{v}},{\boldsymbol{d}},{\boldsymbol{m}})$ a
measure-valued solution to the Ericksen–Leslie equations (1), if
${\boldsymbol{v}}\in\operatorname{\mathcal{C}}_{w}([0,T];{L}^{2}_{\sigma}(\Omega))\cap
L^{2}(0,T;H^{1}_{0,\sigma}(\Omega))$,
${{\boldsymbol{d}}\in\operatorname{\mathcal{C}}_{w}([0,T];H^{1}(\Omega;\mathbb{S}^{2}))}$,
as well as ${\boldsymbol{m}}\in
L^{\infty}_{w^{*}}(0,T;\mathcal{M}(\overline{{\Omega}};\operatorname{\mathbb{R}}^{3\times
3}_{\text{sym},+}))$, and if (6) is fulfilled as well as the measure-valued
formulation of the Navier–Stokes equation

$\displaystyle\begin{split}\int_{0}^{T}\Big{[}\int_{\Omega}{\boldsymbol{T}}_{L}:(\nabla\tilde{{\boldsymbol{v}}})_{\text{sym}}-({\boldsymbol{v}}\otimes{\boldsymbol{v}}+\nabla{\boldsymbol{d}}^{T}&\nabla{\boldsymbol{d}}):\nabla\tilde{{\boldsymbol{v}}}\\\
&-{\boldsymbol{v}}\cdot\partial_{t}\tilde{{\boldsymbol{v}}}\text{d}x-2\int_{\overline{{\Omega}}}(\nabla\tilde{{\boldsymbol{v}}})_{\text{sym}}:\text{d}{\boldsymbol{m}}(x)\Big{]}\text{d}s=0\end{split}$
(7)

for all
$\tilde{{\boldsymbol{v}}}\in\operatorname{\mathcal{C}}^{\infty}_{c}(\Omega\times(0,T))$
and the energy-inequality

$\displaystyle\begin{split}\frac{1}{2}\left(\|{\boldsymbol{v}}\|_{L^{2}(\Omega)}^{2}+\|\nabla{\boldsymbol{d}}\|_{L^{2}(\Omega)}^{2}+\int_{\overline{{\Omega}}}I:\text{d}{\boldsymbol{m}}(x)\right)&\Big{|}_{0}^{t}\\\
+\int_{0}^{t}\int_{\Omega}{\boldsymbol{T}}^{D}:&(\nabla{\boldsymbol{v}})_{\text{sym}}+|{\boldsymbol{d}}\times\Delta{\boldsymbol{d}}|^{2}\text{d}x\text{d}s\leq
0\end{split}$ (8)

for a.e. $t\in(0,T)$.

######  Theorem 2.4.

Let $({\boldsymbol{v}},{\boldsymbol{d}},E)$ be an energy-variational solution
with
${\mathcal{K}(\tilde{{\boldsymbol{v}}})=2\|(\nabla\tilde{{\boldsymbol{v}}})_{\text{sym},-}\|_{\operatorname{\mathcal{C}}(\overline{{\Omega}})}}$
in the sense of Definition 2.1 such that
$E(0)=\mathcal{E}({\boldsymbol{v}}_{0},{\boldsymbol{d}}_{0})$. Then there
exists a measure-valued solution
$({\boldsymbol{v}},{\boldsymbol{d}},{\boldsymbol{m}})$ in the sense of
Definition 2.3 with

$\frac{1}{2}\langle{\boldsymbol{m}},I\rangle=\frac{1}{2}\int_{\overline{{\Omega}}}I:\text{d}{\boldsymbol{m}}=\frac{1}{2}\int_{\overline{{\Omega}}}\text{d}\operatorname{tr}({\boldsymbol{m}})\leq
E-\mathcal{E}({\boldsymbol{v}},{\boldsymbol{d}})\quad\text{a.e.\leavevmode\nobreak\
in }(0,T)\,.$

######  Theorem 2.5.

Let $({\boldsymbol{v}},{\boldsymbol{d}},E)$ be an energy-variational solution
in the sense of Definition 2.1 with
$E(0)=\mathcal{E}({\boldsymbol{v}}_{0},{\boldsymbol{d}}_{0})$ and let
$(\tilde{{\boldsymbol{v}}},\tilde{{\boldsymbol{d}}})$ fulfill

$\displaystyle\begin{split}\tilde{{\boldsymbol{v}}}\in{}&L^{2}(0,T;H^{1}_{0,\sigma}(\Omega)\cap
L^{\infty}(\Omega)\cap W^{1,3}(\Omega))\cap
L^{1}(0,T;\operatorname{\mathcal{C}}^{1}(\overline{{\Omega}}))\\\
&\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\cap
H^{1}(0,T;(H^{1}_{0,\sigma}(\Omega))^{*})\end{split}$ (9a)
$\displaystyle\tilde{{\boldsymbol{d}}}\in{}$ $\displaystyle
L^{1}(0,T;W^{2,\infty}(\Omega))\cap L^{2}(0,T;W^{3,3}(\Omega))\cap
H^{1}(0,T;H^{1}(\Omega))\,$ (9b)

such that $\mathcal{K}(\tilde{{\boldsymbol{v}}})\in L^{1}(0,T)$. Then the
relative energy inequality

$\displaystyle\begin{split}\frac{1}{2}&\left[\|{\boldsymbol{v}}(t)-\tilde{{\boldsymbol{v}}}(t)\|_{L^{2}(\Omega)}^{2}+\|\nabla{\boldsymbol{d}}(t)-\nabla\tilde{{\boldsymbol{d}}}(t)\|_{L^{2}(\Omega)}^{2}\right]\\\
+&\frac{1}{2}\int_{0}^{t}\int_{\Omega}(\mu_{1}+\lambda^{2})({\boldsymbol{d}}\cdot(\nabla{\boldsymbol{v}}-\nabla\tilde{{\boldsymbol{v}}})_{sym}{\boldsymbol{d}})^{2}+\mu_{4}|(\nabla{\boldsymbol{v}}-\nabla\tilde{{\boldsymbol{v}}})_{sym}|^{2}\mathop{}\\!\mathrm{d}x\text{{e}}^{\int_{s}^{t}\mathcal{Q}(\tilde{{\boldsymbol{v}}},\tilde{{\boldsymbol{d}}})\text{d}\tau}\mathop{}\\!\mathrm{d}s\\\
+&\frac{1}{2}\int_{0}^{t}\int_{\Omega}(\mu_{5}+\mu_{6}-\lambda^{2})|(\nabla{\boldsymbol{v}}-\nabla\tilde{{\boldsymbol{v}}})_{sym}{\boldsymbol{d}}|^{2}+|{\boldsymbol{d}}\times\Delta({\boldsymbol{d}}-\tilde{{\boldsymbol{d}}})|^{2}\mathop{}\\!\mathrm{d}x\text{{e}}^{\int_{s}^{t}\mathcal{Q}(\tilde{{\boldsymbol{v}}},\tilde{{\boldsymbol{d}}})\text{d}\tau}\mathop{}\\!\mathrm{d}s\\\
+&\int_{0}^{t}\left\langle\partial_{t}\tilde{{\boldsymbol{v}}}+(\tilde{{\boldsymbol{v}}}\cdot\nabla)\tilde{{\boldsymbol{v}}}+\operatorname{\nabla\cdot}(\nabla\tilde{{\boldsymbol{d}}}^{T}\nabla\tilde{{\boldsymbol{d}}})-\operatorname{\nabla\cdot}\tilde{{\boldsymbol{T}}}^{L},{\boldsymbol{v}}-\tilde{{\boldsymbol{v}}}\right\rangle\text{{e}}^{\int_{s}^{t}\mathcal{Q}(\tilde{{\boldsymbol{v}}},\tilde{{\boldsymbol{d}}})\text{d}\tau}\mathop{}\\!\mathrm{d}s\\\
+&\int_{0}^{t}\left\langle\nabla\left(\partial_{t}\tilde{{\boldsymbol{d}}}+(\tilde{{\boldsymbol{v}}}\cdot\nabla)\tilde{{\boldsymbol{d}}}-(\nabla\tilde{{\boldsymbol{v}}})_{\text{skw}}\tilde{{\boldsymbol{d}}}\right),\nabla({\boldsymbol{d}}-\tilde{{\boldsymbol{d}}})\right\rangle\text{{e}}^{\int_{s}^{t}\mathcal{Q}(\tilde{{\boldsymbol{v}}},\tilde{{\boldsymbol{d}}})\text{d}\tau}\mathop{}\\!\mathrm{d}s\\\
+&\int_{0}^{t}\left\langle\nabla\left((I-\tilde{{\boldsymbol{d}}}\otimes\tilde{{\boldsymbol{d}}})(\lambda(\nabla\tilde{{\boldsymbol{v}}})_{\text{sym}}-\Delta\tilde{{\boldsymbol{d}}})\right),\nabla({\boldsymbol{d}}-\tilde{{\boldsymbol{d}}})\right\rangle\text{{e}}^{\int_{s}^{t}\mathcal{Q}(\tilde{{\boldsymbol{v}}},\tilde{{\boldsymbol{d}}})\text{d}\tau}\mathop{}\\!\mathrm{d}s\\\
\leq&\frac{1}{2}\left[\|{\boldsymbol{v}}_{0}-\tilde{{\boldsymbol{v}}}(0)\|_{L^{2}(\Omega)}^{2}+\|\nabla{\boldsymbol{d}}_{0}-\nabla\tilde{{\boldsymbol{d}}}(0)\|_{L^{2}(\Omega)}^{2}\right]\text{{e}}^{\int_{0}^{t}\mathcal{Q}(\tilde{{\boldsymbol{v}}},\tilde{{\boldsymbol{d}}})\text{d}s}\,,\end{split}$
(10)

holds for a.e. $t\in(0,T)$, where

$\displaystyle\mathcal{Q}(\tilde{{\boldsymbol{v}}},\tilde{{\boldsymbol{d}}}):=$
$\displaystyle\|\tilde{{\boldsymbol{v}}}\|_{L^{\infty}(\Omega)}^{2}+\|\nabla\tilde{{\boldsymbol{v}}}\|_{L^{\infty}(\Omega)}+\|\Delta\tilde{{\boldsymbol{d}}}\|_{L^{3}(\Omega)}^{2}+\|\nabla\tilde{{\boldsymbol{v}}}\|_{L^{3}(\Omega)}^{2}$
$\displaystyle+\|\Delta\tilde{{\boldsymbol{d}}}\|_{L^{\infty}(\Omega)}+\|\Delta\tilde{{\boldsymbol{d}}}\|_{W^{1,3}(\Omega)}^{2}+\mathcal{K}(\tilde{{\boldsymbol{v}}})\,.$

Additionally, if $(\tilde{{\boldsymbol{v}}},\tilde{{\boldsymbol{d}}})$ is a
classical solution to the Ericksen–Leslie system on $(0,T^{*})$, then both
solutions coincide on this time-interval, i.e.,

$\displaystyle{\boldsymbol{v}}\equiv\tilde{{\boldsymbol{v}}}\qquad{\boldsymbol{d}}\equiv\tilde{{\boldsymbol{d}}}\quad\text{
a.e.\leavevmode\nobreak\ on\leavevmode\nobreak\ }(0,T^{*})\,.$

######  Remark (Dissipative solutions).

Theorem 2.5 does not only show the weak-strong uniqueness of energy-
variational solutions, it also shows that every energy-variational solution is
a dissipative solution. A dissipative solution is defined via the relative-
energy inequality (10), i.e., $({\boldsymbol{v}},{\boldsymbol{d}})$ is called
a dissipative solution, if it fulfills the relative energy-inequality 10 for
all function $(\tilde{{\boldsymbol{v}}},\tilde{{\boldsymbol{d}}})$ fulfilling
the regularity properties (9). The concept of energy-variational solutions is
thus strictly finer than the concept of dissipative solutions.

###  2.2 Proofs for the continuous system

###### Proof of Theorem 2.4.

Let $({\boldsymbol{v}},{\boldsymbol{d}},E)$ be an energy-variational solution
according to Definition 2.1 with
$\mathcal{K}(\tilde{{\boldsymbol{v}}})=2\|\nabla\tilde{{\boldsymbol{v}}}\|_{\operatorname{\mathcal{C}}(\overline{{\Omega}})}$.
Then, we may choose $\tilde{{\boldsymbol{v}}}=0$ in order to infer

$\displaystyle
E\Big{|}_{0}^{t}+\int_{0}^{t}\int_{\Omega}(\mu_{1}+\lambda^{2})({\boldsymbol{d}}\cdot(\nabla{\boldsymbol{v}})_{\text{sym}}{\boldsymbol{d}})^{2}+\mu_{4}|(\nabla{\boldsymbol{v}})_{\text{sym}}|^{2}\mathop{}\\!\mathrm{d}x\mathop{}\\!\mathrm{d}\tau$
$\displaystyle+\int_{0}^{t}\int_{\Omega}(\mu_{5}+\mu_{6}-\lambda^{2})|(\nabla{\boldsymbol{v}})_{\text{sym}}{\boldsymbol{d}}|^{2}+|{\boldsymbol{d}}\times\Delta{\boldsymbol{d}}|^{2}\mathop{}\\!\mathrm{d}x\mathop{}\\!\mathrm{d}\tau$
$\displaystyle{}\leq 0\,.$ (11)

Furthermore, we may choose
$\tilde{{\boldsymbol{v}}}=\alpha\tilde{{\boldsymbol{u}}}$ with $\alpha\geq 0$
and multiply the energy-variational inequality (5) by $1/\alpha$. In the limit
$\alpha\operatorname{\rightarrow}\infty$, we find

$-\int_{\Omega}{\boldsymbol{v}}\cdot\tilde{{\boldsymbol{u}}}\mathop{}\\!\mathrm{d}x\Big{|}_{0}^{t}+\int_{0}^{t}\int_{\Omega}{\boldsymbol{v}}\cdot\partial_{t}\tilde{{\boldsymbol{u}}}+\left[{\boldsymbol{v}}\otimes{\boldsymbol{v}}+\nabla{\boldsymbol{d}}^{T}\nabla{\boldsymbol{d}}-{\boldsymbol{T}}_{L}\right]:\nabla\tilde{{\boldsymbol{u}}}\mathop{}\\!\mathrm{d}x\mathop{}\\!\mathrm{d}s\\\
+\int_{0}^{t}2\|(\nabla\tilde{{\boldsymbol{u}}})_{\text{sym},-}\|_{\operatorname{\mathcal{C}}(\overline{{\Omega}};\operatorname{\mathbb{R}}^{3\times
3})}\left[\mathcal{E}({\boldsymbol{v}},{\boldsymbol{d}})-E\right]\mathop{}\\!\mathrm{d}s\leq
0\,$

for all
$\tilde{{\boldsymbol{u}}}\in\operatorname{\mathcal{C}}^{1}([0,T])\otimes\operatorname{\mathcal{C}}^{1}_{0}(\Omega)$.
By the definition of the weak time derivative, we infer that

$\displaystyle-\int_{0}^{T}\langle\partial_{t}{\boldsymbol{v}},{}$
$\displaystyle\tilde{{\boldsymbol{u}}}\rangle\text{d}t=-\int_{\Omega}{\boldsymbol{v}}\cdot\tilde{{\boldsymbol{u}}}\mathop{}\\!\mathrm{d}x\Big{|}_{0}^{T}+\int_{0}^{T}\int_{\Omega}{\boldsymbol{v}}\cdot\partial_{t}\tilde{{\boldsymbol{u}}}\mathop{}\\!\mathrm{d}x\mathop{}\\!\mathrm{d}t$
$\displaystyle\leq{}$
$\displaystyle\left(\|{\boldsymbol{v}}\|_{L^{\infty}(0,T;L^{2}(\Omega))}^{2}+\|\nabla{\boldsymbol{d}}\|_{L^{\infty}(0,T;L^{2}(\Omega))}^{2}\right)\|\nabla\tilde{{\boldsymbol{u}}}\|_{L^{1}(0,T;\operatorname{\mathcal{C}}(\overline{{\Omega}}))}$
$\displaystyle+$ $\displaystyle
2\|\nabla\tilde{{\boldsymbol{u}}}\|_{L^{1}(0,T;\operatorname{\mathcal{C}}(\overline{{\Omega}}))}\left\|E-\mathcal{E}({\boldsymbol{v}},{\boldsymbol{d}})\right\|_{L^{\infty}(0,T)}+\|{\boldsymbol{T}}_{L}\|_{L^{2}(0,T;L^{2}(\Omega))}\|\nabla\tilde{{\boldsymbol{u}}}\|_{L^{2}(0,T;L^{2}(\Omega))}\,.$

First this is only well defined for
$\tilde{{\boldsymbol{u}}}\in\operatorname{\mathcal{C}}^{1}([0,T])\otimes\operatorname{\mathcal{C}}^{1}_{0}(\Omega)$.
Since this is a linear subspace of the Banach space
$L^{1}(0,T;\operatorname{\mathcal{C}}^{1}_{0}(\Omega))\cap
L^{2}(0,T;H^{1}(\Omega))$, we infer from the Hahn–Banach theorem [8, Thm. 1.1]
the existence of an element

$\partial_{t}{\boldsymbol{v}}\in\left(L^{1}(0,T;\operatorname{\mathcal{C}}^{1}_{0}(\Omega))\cap
L^{2}(0,T;H^{1}(\Omega))\right)^{*}=L^{\infty}_{w^{*}}(0,T;(\operatorname{\mathcal{C}}^{1}_{0}(\Omega))^{*})+L^{2}(0,T;W^{-1,2}(\Omega))\,.$

This allows to write the linear form

$\langle{\boldsymbol{l}}(t),\tilde{{\boldsymbol{u}}}(t)\rangle_{L^{\infty}_{w^{*}}(0,T;(\operatorname{\mathcal{C}}^{1}_{0}(\Omega))^{*}),L^{1}(0,T;\operatorname{\mathcal{C}}^{1}_{0}(\Omega))}:={}\\\
-\int_{0}^{T}\langle\partial_{t}{\boldsymbol{v}}(t),\tilde{{\boldsymbol{u}}}(t)\rangle\mathop{}\\!\mathrm{d}t+\int_{0}^{T}\int_{\Omega}\left[{\boldsymbol{v}}(t)\otimes{\boldsymbol{v}}(t)+\nabla{\boldsymbol{d}}^{T}(t)\nabla{\boldsymbol{d}}(t)-{\boldsymbol{T}}_{L}(t)\right]:\nabla\tilde{{\boldsymbol{u}}}(t)\mathop{}\\!\mathrm{d}x\mathop{}\\!\mathrm{d}t\\\
\leq{}2\int_{0}^{T}\|(\nabla\tilde{{\boldsymbol{u}}}(t))_{\text{sym},-}\|_{L^{1}(0,T;\operatorname{\mathcal{C}}(\overline{{\Omega}};\operatorname{\mathbb{R}}^{3\times
3}))}\left[E(t)-\mathcal{E}({\boldsymbol{v}}(t),{\boldsymbol{d}}(t))\right]\text{d}t\,$
(12)

for a.e. $t\in(0,T)$. This implies that the linear form ${\boldsymbol{l}}\in
L^{\infty}_{w^{*}}(0,T;(\operatorname{\mathcal{C}}^{1}_{0}(\Omega))^{*})$ can
be identified via considering the gradient operator

$\nabla:L^{1}(0,T;\operatorname{\mathcal{C}}^{1}_{0}(\Omega))\operatorname{\rightarrow}L^{1}(0,T;\operatorname{\mathcal{C}}(\overline{{\Omega}};\operatorname{\mathbb{R}}^{3\times
3}))\,,$

we set $G=\nabla(L^{1}(0,T;\operatorname{\mathcal{C}}^{1}_{0}(\Omega)))$ the
range of the gradient, the set of all gradient fields equipped with the norm
of
$L^{1}(0,T;\operatorname{\mathcal{C}}(\overline{{\Omega}};\operatorname{\mathbb{R}}^{3\times
3}))$. Note that this mapping is injective due to the prescribed Dirichlet
boundary conditions. We may set
$S=\nabla^{-1}:G\operatorname{\rightarrow}L^{1}(0,T;\operatorname{\mathcal{C}}^{1}_{0}(\Omega))$
such that for $h\in G$ the functional
$h\mapsto\langle{\boldsymbol{l}},Sh\rangle$ is a continuous linear functional
on $G$. By the Hahn–Banach theorem [8, Thm. 1.1], this linear functional can
be extended to a continuous linear functional $\Phi$ on
$(L^{1}(0,T;\operatorname{\mathcal{C}}(\overline{{\Omega}};\operatorname{\mathbb{R}}^{3\times
3})))^{*}=L^{\infty}_{w}(0,T;\mathcal{M}(\overline{{\Omega}};\operatorname{\mathbb{R}}^{3\times
3}))$ with

$\displaystyle\langle\Phi(t);{\boldsymbol{A}}\rangle\leq
2\|({\boldsymbol{A}})_{\text{sym},-}\|_{\operatorname{\mathcal{C}}(\overline{{\Omega}};\operatorname{\mathbb{R}}^{3\times
3})}[E(t)-\mathcal{E}({\boldsymbol{v}}(t),{\boldsymbol{d}}(t))]\,$ (13)

for all
${\boldsymbol{A}}\in\operatorname{\mathcal{C}}(\overline{{\Omega}};\operatorname{\mathbb{R}}^{3\times
3})$ a.e. in $(0,T)$, where the estimate follows from (12) and Hahn–Banach’s
theorem written in its point-wise a.e.-form. By the Riesz representation
theorem, we know that there exists a measure ${\boldsymbol{m}}\in
L^{\infty}_{w^{*}}(0,T;\mathcal{M}(\overline{{\Omega}};\operatorname{\mathbb{R}}^{3\times
3})$ such that

$\displaystyle\langle{\boldsymbol{l}}(t),\tilde{{\boldsymbol{u}}}(t)\rangle_{L^{\infty}_{w}(0,T:(\operatorname{\mathcal{C}}^{1}_{0}(\Omega))^{*},),L^{1}(0,T;\operatorname{\mathcal{C}}^{1}_{0}(\Omega))}$
$\displaystyle=\langle\Phi(t),\nabla\tilde{{\boldsymbol{u}}}(t)\rangle_{L^{\infty}_{w}(0,T;(\operatorname{\mathcal{C}}(\overline{{\Omega}}))^{*},),L^{1}(0,T;\operatorname{\mathcal{C}}(\overline{{\Omega}}))}$
$\displaystyle=-\int_{0}^{T}\int_{\overline{{\Omega}}}\nabla\tilde{{\boldsymbol{u}}}(t):\mathop{}\\!\mathrm{d}{\boldsymbol{m}}_{t}\text{d}t\,,$

which implies the formulation (7) by the definition of ${\boldsymbol{l}}$ in
(12). From the estimate (13), we infer that the measure ${\boldsymbol{m}}$
vanishes for all continuous functions with values in the skew-symmetric
matrices, such that ${\boldsymbol{m}}\in
L^{\infty}_{w}(0,T;\mathcal{M}(\overline{{\Omega}};\operatorname{\mathbb{R}}^{3\times
3}_{\text{sym}}))$. Additionally, (13) allows to insert
${\boldsymbol{A}}={\boldsymbol{\xi}}\otimes{\boldsymbol{\xi}}\phi$ for any
${\boldsymbol{\xi}}\in\operatorname{\mathbb{R}}^{3}$ and
$\phi\in\operatorname{\mathcal{C}}(\overline{{\Omega}};[0,\infty))$ such that
$\int_{\overline{{\Omega}}}\phi{\boldsymbol{\xi}}\otimes{\boldsymbol{\xi}}:\text{d}{\boldsymbol{m}}_{t}\geq
0$ for all
$\phi\in\operatorname{\mathcal{C}}(\overline{{\Omega}};[0,\infty))$, which is
the definition of a positive definite measure ${\boldsymbol{m}}\in
L^{\infty}_{w}(0,T;\mathcal{M}(\overline{{\Omega}};\operatorname{\mathbb{R}}^{3\times
3}_{\text{sym},+}))$. From the estimate in (13), we additionally infer by
choosing $\Phi(t)=\phi(t)I$ for
$\phi\in\operatorname{\mathcal{C}}^{1}_{0}([0,T])$ with $\phi\geq 0$ a.e. in
$(0,T)$ that

$\displaystyle\int_{0}^{T}\phi(t)\int_{\overline{{\Omega}}}I:\text{d}{\boldsymbol{m}}(x)\text{d}t=\int_{0}^{T}\phi(t)\int_{\overline{{\Omega}}}\operatorname{tr}({\boldsymbol{m}})(x)\text{d}t\leq
2\int_{0}^{T}\phi(t)[E(t)-\mathcal{E}({\boldsymbol{v}}(t),{\boldsymbol{d}}(t))]\text{d}t\,.$

This implies that $\langle{\boldsymbol{m}},I\rangle\leq
2[E-\mathcal{E}({\boldsymbol{v}},{\boldsymbol{d}})]$ for a.e. $t\in(0,T)$,
which in turn implies together with (11) the inequality (8).

Note that we use the usual spectral-norm for symmetric matrices
$|A|_{2}=\max_{i\in\\{1,2,3\\}}|\lambda_{i}|$, where $\lambda_{i}$ denote the
eigenvalues of the matrix $A\in\operatorname{\mathbb{R}}^{3\times
3}_{\text{sym}}$. The dual-norm with respect to the Frobenius product $:$ is
given by $|A|_{2}^{\prime}=\sum_{i=1}^{3}|\lambda_{i}|$ for
$A\in\operatorname{\mathbb{R}}^{3\times 3}$. For a positive definite matrix
all eigenvalues are non-negative, such that we may write
$|A|_{2}^{\prime}=\sum_{i=1}^{3}\lambda_{i}=\operatorname{tr}(A)=A:I$. This
implies that we may identify the norm $|A|_{2}^{\prime}=I:A$.

∎

###### Proof of Theorem 2.5.

We observe that (6) tested with $\Delta\tilde{{\boldsymbol{d}}}$ and
integrated-in-time gives

$\displaystyle\int_{0}^{t}\int_{\Omega}\left(\partial_{t}{\boldsymbol{d}}+({\boldsymbol{v}}\cdot\nabla){\boldsymbol{d}}-(\nabla{\boldsymbol{v}})_{\text{skw}}{\boldsymbol{d}}+(I-{\boldsymbol{d}}\otimes{\boldsymbol{d}})(\lambda(\nabla{\boldsymbol{v}})_{\text{sym}}{\boldsymbol{d}}-\Delta{\boldsymbol{d}})\right)\Delta\tilde{{\boldsymbol{d}}}\text{d}x\text{d}s=0\,.$
(14)

Additionally, the equation (1a) for the strong solution tested with
$\tilde{{\boldsymbol{v}}}-{\boldsymbol{v}}$ and (1c) for the strong solution
tested with $\Delta({\boldsymbol{d}}-\tilde{{\boldsymbol{d}}})$ implies

$\displaystyle 0{}$
$\displaystyle=\frac{1}{2}\|\tilde{{\boldsymbol{v}}}\|_{L^{2}(\Omega)}^{2}\Big{|}_{0}^{t}-\int_{0}^{t}\int_{\Omega}\partial_{t}\tilde{{\boldsymbol{v}}}\cdot{\boldsymbol{v}}+(\tilde{{\boldsymbol{v}}}\cdot\nabla)\tilde{{\boldsymbol{v}}}\cdot{\boldsymbol{v}}-\nabla\tilde{{\boldsymbol{d}}}^{T}\nabla\tilde{{\boldsymbol{d}}}:\nabla{\boldsymbol{v}}+\mu_{4}(\nabla\tilde{{\boldsymbol{v}}})_{\text{sym}}:(\nabla{\boldsymbol{v}}-\nabla\tilde{{\boldsymbol{v}}})\text{d}x\text{d}s$
$\displaystyle-\int_{0}^{t}\int_{\Omega}(\mu_{1}+\lambda^{2})(\tilde{{\boldsymbol{d}}}\cdot(\nabla\tilde{{\boldsymbol{v}}})_{\text{sym}}\tilde{{\boldsymbol{d}}})(\tilde{{\boldsymbol{d}}}\otimes\tilde{{\boldsymbol{d}}}):(\nabla{\boldsymbol{v}}-\nabla\tilde{{\boldsymbol{v}}})\text{d}x\text{d}s$
$\displaystyle-\int_{0}^{t}\int_{\Omega}(\mu_{5}+\mu_{6}-\lambda^{2})(\tilde{{\boldsymbol{d}}}\otimes(\nabla\tilde{{\boldsymbol{v}}})_{\text{sym}}\tilde{{\boldsymbol{d}}}):(\nabla{\boldsymbol{v}}-\nabla\tilde{{\boldsymbol{v}}})\text{d}x\text{d}s$
$\displaystyle-\int_{0}^{t}\int_{\Omega}\lambda\left((\nabla{\boldsymbol{v}})_{\text{sym}}\tilde{{\boldsymbol{d}}}\cdot(I-\tilde{{\boldsymbol{d}}}\otimes\tilde{{\boldsymbol{d}}})\Delta\tilde{{\boldsymbol{d}}}+\nabla\left((I-\tilde{{\boldsymbol{d}}}\otimes\tilde{{\boldsymbol{d}}})(\nabla\tilde{{\boldsymbol{v}}})_{\text{sym}}\tilde{{\boldsymbol{d}}}\right):\nabla{\boldsymbol{d}}\right)\text{d}x\text{d}s$
$\displaystyle+\frac{1}{2}\|\nabla\tilde{{\boldsymbol{d}}}\|_{L^{2}(\Omega)}^{2}\Big{|}_{0}^{t}-\int_{0}^{t}\int_{\Omega}\nabla\left((I-\tilde{{\boldsymbol{d}}}\otimes\tilde{{\boldsymbol{d}}})\Delta\tilde{{\boldsymbol{d}}}\right):(\nabla{\boldsymbol{d}}-\nabla\tilde{{\boldsymbol{d}}})\mathop{}\\!\mathrm{d}x\mathop{}\\!\mathrm{d}t$
$\displaystyle-\int_{0}^{t}\int_{\Omega}\partial_{t}\nabla\tilde{{\boldsymbol{d}}}:\nabla{\boldsymbol{d}}+\nabla\left((\tilde{{\boldsymbol{v}}}\cdot\nabla)\tilde{{\boldsymbol{d}}}\right):\nabla{\boldsymbol{d}}-(\nabla{\boldsymbol{v}})_{\text{skw}}\tilde{{\boldsymbol{d}}}\cdot\Delta\tilde{{\boldsymbol{d}}}-\nabla\left((\nabla\tilde{{\boldsymbol{v}}})_{\text{skw}}\tilde{{\boldsymbol{d}}}\right):\nabla{\boldsymbol{d}}\text{d}x\text{d}s$
$\displaystyle+\int_{0}^{t}\left\langle\partial_{t}\tilde{{\boldsymbol{v}}}+(\tilde{{\boldsymbol{v}}}\cdot\nabla)\tilde{{\boldsymbol{v}}}+\operatorname{\nabla\cdot}(\nabla\tilde{{\boldsymbol{d}}}^{T}\nabla\tilde{{\boldsymbol{d}}})-\operatorname{\nabla\cdot}\tilde{{\boldsymbol{T}}}^{L},{\boldsymbol{v}}-\tilde{{\boldsymbol{v}}}\right\rangle\mathop{}\\!\mathrm{d}s$
$\displaystyle+\int_{0}^{t}\left\langle\nabla\left(\partial_{t}\tilde{{\boldsymbol{d}}}+(\tilde{{\boldsymbol{v}}}\cdot\nabla)\tilde{{\boldsymbol{d}}}-(\nabla\tilde{{\boldsymbol{v}}})_{\text{skw}}\tilde{{\boldsymbol{d}}}+(I-\tilde{{\boldsymbol{d}}}\otimes\tilde{{\boldsymbol{d}}})(\lambda(\nabla\tilde{{\boldsymbol{v}}})_{\text{sym}}\tilde{{\boldsymbol{d}}}-\Delta\tilde{{\boldsymbol{d}}})\right),\nabla({\boldsymbol{d}}-\tilde{{\boldsymbol{d}}})\right\rangle\mathop{}\\!\mathrm{d}s\,.$
(15)

We observe that

$\displaystyle-\int_{0}^{t}\int_{\Omega}\partial_{t}\nabla\tilde{{\boldsymbol{d}}}:\nabla{\boldsymbol{d}}-\partial_{t}{\boldsymbol{d}}\cdot\Delta\tilde{{\boldsymbol{d}}}\text{d}x\text{d}s=\int_{\Omega}\nabla{\boldsymbol{d}}:\nabla\tilde{{\boldsymbol{d}}}\text{d}x\Big{|}_{0}^{t}\,.$

Adding up (5), (14), and (15) implies

$\displaystyle[E-$
$\displaystyle\mathcal{E}({\boldsymbol{v}},{\boldsymbol{d}})]\Big{|}_{0-}^{t}+\frac{1}{2}\left[\|{\boldsymbol{v}}-\tilde{{\boldsymbol{v}}}\|_{L^{2}(\Omega)}^{2}+\|\nabla{\boldsymbol{d}}-\nabla\tilde{{\boldsymbol{d}}}\|_{L^{2}(\Omega)}^{2}\right]\Big{|}_{0}^{t}$
$\displaystyle+$
$\displaystyle\int_{0}^{t}\int_{\Omega}(\mu_{1}+\lambda^{2})({\boldsymbol{d}}\cdot(\nabla{\boldsymbol{v}}-\nabla\tilde{{\boldsymbol{v}}})_{sym}{\boldsymbol{d}})^{2}+\mu_{4}|(\nabla{\boldsymbol{v}}-\nabla\tilde{{\boldsymbol{v}}})_{sym}|^{2}\mathop{}\\!\mathrm{d}x\mathop{}\\!\mathrm{d}s$
$\displaystyle+$
$\displaystyle\int_{0}^{t}\int_{\Omega}(\mu_{5}+\mu_{6}-\lambda^{2})|(\nabla{\boldsymbol{v}}-\nabla\tilde{{\boldsymbol{v}}})_{sym}{\boldsymbol{d}}|^{2}+|{\boldsymbol{d}}\times\Delta({\boldsymbol{d}}-\tilde{{\boldsymbol{d}}})|^{2}\mathop{}\\!\mathrm{d}x\mathop{}\\!\mathrm{d}s$
$\displaystyle+$
$\displaystyle\int_{0}^{t}\left\langle\partial_{t}\tilde{{\boldsymbol{v}}}+(\tilde{{\boldsymbol{v}}}\cdot\nabla)\tilde{{\boldsymbol{v}}}+\operatorname{\nabla\cdot}(\nabla\tilde{{\boldsymbol{d}}}^{T}\nabla\tilde{{\boldsymbol{d}}})-\operatorname{\nabla\cdot}\tilde{{\boldsymbol{T}}}^{L},{\boldsymbol{v}}-\tilde{{\boldsymbol{v}}}\right\rangle\mathop{}\\!\mathrm{d}s$
$\displaystyle+$
$\displaystyle\int_{0}^{t}\left\langle\nabla\left(\partial_{t}\tilde{{\boldsymbol{d}}}+(\tilde{{\boldsymbol{v}}}\cdot\nabla)\tilde{{\boldsymbol{d}}}-(\nabla\tilde{{\boldsymbol{v}}})_{\text{skw}}\tilde{{\boldsymbol{d}}}+(I-\tilde{{\boldsymbol{d}}}\otimes\tilde{{\boldsymbol{d}}})(\lambda(\nabla\tilde{{\boldsymbol{v}}})_{\text{sym}}-\Delta\tilde{{\boldsymbol{d}}})\right),\nabla({\boldsymbol{d}}-\tilde{{\boldsymbol{d}}})\right\rangle\mathop{}\\!\mathrm{d}s$
$\displaystyle\leq{}$
$\displaystyle\int_{0}^{t}\int_{\Omega}({\boldsymbol{v}}\cdot\nabla){\boldsymbol{v}}\cdot\tilde{{\boldsymbol{v}}}+(\tilde{{\boldsymbol{v}}}\cdot\nabla)\tilde{{\boldsymbol{v}}}\cdot{\boldsymbol{v}}\text{d}x\text{d}s$
$\displaystyle+\int_{0}^{t}\int_{\Omega}{\boldsymbol{T}}^{E}({\boldsymbol{d}}):\nabla\tilde{{\boldsymbol{v}}}-\nabla\tilde{{\boldsymbol{d}}}^{T}\nabla\tilde{{\boldsymbol{d}}}:\nabla{\boldsymbol{v}}-({\boldsymbol{v}}\cdot\nabla){\boldsymbol{d}}\cdot\Delta\tilde{{\boldsymbol{d}}}+\nabla\left((\tilde{{\boldsymbol{v}}}\cdot\nabla)\tilde{{\boldsymbol{d}}}\right):\nabla{\boldsymbol{d}}\text{d}x\text{d}s$
$\displaystyle-\int_{0}^{t}\int_{\Omega}(\mu_{1}+\lambda^{2})\left({\boldsymbol{d}}\cdot(\nabla\tilde{{\boldsymbol{v}}})_{\text{sym}}{\boldsymbol{d}}({\boldsymbol{d}}\otimes{\boldsymbol{d}})-\tilde{{\boldsymbol{d}}}\cdot(\nabla\tilde{{\boldsymbol{v}}})_{\text{sym}\tilde{{\boldsymbol{d}}}}(\tilde{{\boldsymbol{d}}}\otimes\tilde{{\boldsymbol{d}}})\right):(\nabla{\boldsymbol{v}}-\nabla\tilde{{\boldsymbol{v}}})_{\text{sym}}\text{d}x\text{d}s$
$\displaystyle-\int_{0}^{t}\int_{\Omega}(\mu_{5}+\mu_{6}-\lambda^{2})\left((\nabla\tilde{{\boldsymbol{v}}})_{\text{sym}}({\boldsymbol{d}}\otimes{\boldsymbol{d}}-\tilde{{\boldsymbol{d}}}\otimes\tilde{{\boldsymbol{d}}}):(\nabla{\boldsymbol{v}}-\nabla\tilde{{\boldsymbol{v}}})_{\text{sym}}\right)\text{d}x\text{d}s$
$\displaystyle-\int_{0}^{t}\int_{\Omega}\nabla\left[\left({\boldsymbol{d}}\otimes{\boldsymbol{d}}-\tilde{{\boldsymbol{d}}}\otimes\tilde{{\boldsymbol{d}}}\right)\Delta\tilde{{\boldsymbol{d}}}\right]:\left(\nabla{\boldsymbol{d}}-\nabla\tilde{{\boldsymbol{d}}}\right)\text{d}x\text{d}s$
$\displaystyle-\lambda\int_{0}^{t}\int_{\Omega}\left((I-{\boldsymbol{d}}\otimes{\boldsymbol{d}})(\nabla{\boldsymbol{v}})_{\text{sym}}{\boldsymbol{d}}-(I-\tilde{{\boldsymbol{d}}}\otimes\tilde{{\boldsymbol{d}}})(\nabla{\boldsymbol{v}})_{\text{sym}}\tilde{{\boldsymbol{d}}}\right)\cdot\Delta\tilde{{\boldsymbol{d}}}\text{d}x\text{d}s$
$\displaystyle-\lambda\int_{0}^{t}\int_{\Omega}\nabla{\boldsymbol{d}}:\nabla\left((I-{\boldsymbol{d}}\otimes{\boldsymbol{d}})(\nabla\tilde{{\boldsymbol{v}}})_{\text{sym}}{\boldsymbol{d}}-(I-\tilde{{\boldsymbol{d}}}\otimes\tilde{{\boldsymbol{d}}})(\nabla\tilde{{\boldsymbol{v}}})_{\text{sym}}\tilde{{\boldsymbol{d}}}\right)\text{d}x\text{d}s$
$\displaystyle+\int_{0}^{t}\int_{\Omega}\nabla((\nabla\tilde{{\boldsymbol{v}}})_{\text{skw}}({\boldsymbol{d}}-\tilde{{\boldsymbol{d}}})):\nabla{\boldsymbol{d}}+((\nabla{\boldsymbol{v}})_{\text{skw}}({\boldsymbol{d}}-\tilde{{\boldsymbol{d}}}))\cdot\Delta\tilde{{\boldsymbol{d}}}\text{d}x\text{d}s$
$\displaystyle\hskip
142.26378pt+\int_{0}^{t}\mathcal{K}(\tilde{{\boldsymbol{v}}})\left[E-\mathcal{E}({\boldsymbol{v}},{\boldsymbol{d}})\right]\text{d}s\,.$
(16)

We recall the identity, which is a special case of the identity (6.7) in [28],

$\displaystyle
0=\int_{\Omega}\nabla\tilde{{\boldsymbol{v}}}:(\nabla\tilde{{\boldsymbol{d}}}^{T}\nabla{\boldsymbol{d}})+\nabla\tilde{{\boldsymbol{v}}}:(\nabla{\boldsymbol{d}}^{T}\nabla\tilde{{\boldsymbol{d}}})+(\tilde{{\boldsymbol{v}}}\cdot\nabla){\boldsymbol{d}}\cdot\Delta\tilde{{\boldsymbol{d}}}-\nabla\left((\tilde{{\boldsymbol{v}}}\cdot\nabla)\tilde{{\boldsymbol{d}}}\right):\nabla{\boldsymbol{d}}\text{d}x\,.$

this identity helps to transform the second line on the right-hand side of
(16) in the case of ${\boldsymbol{T}}^{E}={\boldsymbol{T}}^{E}_{2}$ into

$\displaystyle-\int_{\Omega}$
$\displaystyle\nabla{\boldsymbol{d}}^{T}\nabla{\boldsymbol{d}}:\nabla\tilde{{\boldsymbol{v}}}+\nabla\tilde{{\boldsymbol{d}}}^{T}\nabla\tilde{{\boldsymbol{d}}}:\nabla{\boldsymbol{v}}+({\boldsymbol{v}}\cdot\nabla){\boldsymbol{d}}\cdot\Delta\tilde{{\boldsymbol{d}}}-\nabla\left((\tilde{{\boldsymbol{v}}}\cdot\nabla)\tilde{{\boldsymbol{d}}}\right):\nabla{\boldsymbol{d}}\text{d}x$
$\displaystyle=-\int_{\Omega}(\nabla{\boldsymbol{d}}-\nabla\tilde{{\boldsymbol{d}}})^{T}(\nabla{\boldsymbol{d}}-\nabla\tilde{{\boldsymbol{d}}}):\nabla\tilde{{\boldsymbol{v}}}+(\nabla\tilde{{\boldsymbol{d}}}^{T}\nabla\tilde{{\boldsymbol{d}}}):(\nabla{\boldsymbol{v}}-\nabla\tilde{{\boldsymbol{v}}})+(({\boldsymbol{v}}-\tilde{{\boldsymbol{v}}})\cdot\nabla){\boldsymbol{d}}\cdot\Delta\tilde{{\boldsymbol{d}}}\text{d}x$
$\displaystyle=-\int_{\Omega}(\nabla{\boldsymbol{d}}-\nabla\tilde{{\boldsymbol{d}}})^{T}(\nabla{\boldsymbol{d}}-\nabla\tilde{{\boldsymbol{d}}}):\nabla\tilde{{\boldsymbol{v}}}+(({\boldsymbol{v}}-\tilde{{\boldsymbol{v}}})\cdot\nabla)({\boldsymbol{d}}-\tilde{{\boldsymbol{d}}})\cdot\Delta\tilde{{\boldsymbol{d}}}\text{d}x\,.$

In case of ${\boldsymbol{T}}^{E}={\boldsymbol{T}}^{E}_{1}$ we may observe from
the uni-norm restriction and additional manipulations that

$\displaystyle\int_{\Omega}$
$\displaystyle(\tilde{{\boldsymbol{v}}}\cdot\nabla){\boldsymbol{d}}[{\boldsymbol{d}}]_{x}[{\boldsymbol{d}}]_{x}\Delta{\boldsymbol{d}}-\nabla\tilde{{\boldsymbol{d}}}^{T}\nabla\tilde{{\boldsymbol{d}}}:\nabla{\boldsymbol{v}}-({\boldsymbol{v}}\cdot\nabla){\boldsymbol{d}}\cdot\Delta\tilde{{\boldsymbol{d}}}+\nabla\left((\tilde{{\boldsymbol{v}}}\cdot\nabla)\tilde{{\boldsymbol{d}}}\right):\nabla{\boldsymbol{d}}\text{d}x\qquad\qquad\qquad$
$\displaystyle={}$
$\displaystyle\int_{\Omega}(\tilde{{\boldsymbol{v}}}\cdot\nabla){\boldsymbol{d}}[{\boldsymbol{d}}]_{x}^{T}[{\boldsymbol{d}}]_{x}\Delta{\boldsymbol{d}}+\left[[\tilde{{\boldsymbol{d}}}]_{x}^{T}[\tilde{{\boldsymbol{d}}}]_{x}\nabla\tilde{{\boldsymbol{d}}}-[{\boldsymbol{d}}]_{x}^{T}[{\boldsymbol{d}}]_{x}\nabla{\boldsymbol{d}}\right]\cdot(\Delta\tilde{{\boldsymbol{d}}}\otimes{\boldsymbol{v}})\text{d}x$
$\displaystyle+\int_{\Omega}\nabla\left((\tilde{{\boldsymbol{v}}}\cdot\nabla)\tilde{{\boldsymbol{d}}}[\tilde{{\boldsymbol{d}}}]_{x}^{T}[\tilde{{\boldsymbol{d}}}]_{x}\right):\nabla{\boldsymbol{d}}\text{d}x$
$\displaystyle={}$
$\displaystyle\int_{\Omega}(\tilde{{\boldsymbol{v}}}\cdot\nabla)({\boldsymbol{d}}-\tilde{{\boldsymbol{d}}})[{\boldsymbol{d}}]_{x}^{T}[{\boldsymbol{d}}]_{x}(\Delta{\boldsymbol{d}}-\Delta\tilde{{\boldsymbol{d}}})-(({\boldsymbol{v}}-\tilde{{\boldsymbol{v}}})\cdot\nabla)({\boldsymbol{d}}-\tilde{{\boldsymbol{d}}})\cdot[{\boldsymbol{d}}]_{x}^{T}[{\boldsymbol{d}}]_{x}\Delta\tilde{{\boldsymbol{d}}}\text{d}x$
$\displaystyle-\int_{\Omega}\nabla\left[(\tilde{{\boldsymbol{v}}}\cdot\nabla)\tilde{{\boldsymbol{d}}}([{\boldsymbol{d}}]_{x}^{T}[{\boldsymbol{d}}]_{x}-[\tilde{{\boldsymbol{d}}}]_{x}^{T}[\tilde{{\boldsymbol{d}}}]_{x})\right]:(\nabla{\boldsymbol{d}}-\nabla\tilde{{\boldsymbol{d}}})\text{d}x$
$\displaystyle+\int_{\Omega}(({\boldsymbol{v}}-\tilde{{\boldsymbol{v}}})\cdot\nabla)\tilde{{\boldsymbol{d}}}\cdot([{\boldsymbol{d}}]_{x}^{T}[{\boldsymbol{d}}]_{x}-[\tilde{{\boldsymbol{d}}}]_{x}^{T}[\tilde{{\boldsymbol{d}}}]_{x})\Delta\tilde{{\boldsymbol{d}}}\text{d}x\,.$

The resulting terms can now be estimated in terms of the relative energy and
the relative dissipation, i.e., the terms on the left-hand side. Similar
arguments have already been applied in [28] and [25]. We remark, that in the
term in the fifth line on the right-hand side of (16), we have to apply the
product rule to infer

$-\int_{\Omega}\nabla\left[\left({\boldsymbol{d}}\otimes{\boldsymbol{d}}-\tilde{{\boldsymbol{d}}}\otimes\tilde{{\boldsymbol{d}}}\right)\Delta\tilde{{\boldsymbol{d}}}\right]:\left(\nabla{\boldsymbol{d}}-\nabla\tilde{{\boldsymbol{d}}}\right)\text{d}x\\\
=-\int_{\Omega}(\nabla{\boldsymbol{d}}-\nabla\tilde{{\boldsymbol{d}}}):\nabla\left({\boldsymbol{d}}\otimes{\boldsymbol{d}}-\tilde{{\boldsymbol{d}}}\otimes\tilde{{\boldsymbol{d}}}\right)\cdot\Delta\tilde{{\boldsymbol{d}}}+\left((\nabla{\boldsymbol{d}}-\nabla\tilde{{\boldsymbol{d}}})^{T}\left({\boldsymbol{d}}\otimes{\boldsymbol{d}}-\tilde{{\boldsymbol{d}}}\otimes\tilde{{\boldsymbol{d}}}\right)\right):\nabla\Delta\tilde{{\boldsymbol{d}}}\text{d}x\,.$

In order to prove the weak-strong uniqueness result, the terms on the right-
hand side of (16) have to estimated by terms on the left-hand side and in the
end Gronwall’s inequality is applied. Since most of the estimates are rather
standard and done using Hölder’s and Young’s inequality, we refrain from
displaying it in full detail here. The interested reader may consult the
article [28] for more details, where these estimates were performed for a more
involved system.

Estimating all terms appropriately, we find that there exists a constant $C>0$
such that

$\displaystyle[E-$
$\displaystyle\mathcal{E}({\boldsymbol{v}},{\boldsymbol{d}})]\Big{|}_{0-}^{t}+\frac{1}{2}\left[\|{\boldsymbol{v}}-\tilde{{\boldsymbol{v}}}\|_{L^{2}(\Omega)}^{2}+\|\nabla{\boldsymbol{d}}-\nabla\tilde{{\boldsymbol{d}}}\|_{L^{2}(\Omega)}^{2}\right]\Big{|}_{0}^{t}$
$\displaystyle+$
$\displaystyle\frac{1}{2}\int_{0}^{t}\int_{\Omega}(\mu_{1}+\lambda^{2})({\boldsymbol{d}}\cdot(\nabla{\boldsymbol{v}}-\nabla\tilde{{\boldsymbol{v}}})_{sym}{\boldsymbol{d}})^{2}+\mu_{4}|(\nabla{\boldsymbol{v}}-\nabla\tilde{{\boldsymbol{v}}})_{sym}|^{2}\mathop{}\\!\mathrm{d}x\mathop{}\\!\mathrm{d}s$
$\displaystyle+$
$\displaystyle\frac{1}{2}\int_{0}^{t}\int_{\Omega}(\mu_{5}+\mu_{6}-\lambda^{2})|(\nabla{\boldsymbol{v}}-\nabla\tilde{{\boldsymbol{v}}})_{sym}{\boldsymbol{d}}|^{2}+|{\boldsymbol{d}}\times\Delta({\boldsymbol{d}}-\tilde{{\boldsymbol{d}}})|^{2}\mathop{}\\!\mathrm{d}x\mathop{}\\!\mathrm{d}s$
$\displaystyle+$
$\displaystyle\int_{0}^{t}\left\langle\partial_{t}\tilde{{\boldsymbol{v}}}+(\tilde{{\boldsymbol{v}}}\cdot\nabla)\tilde{{\boldsymbol{v}}}+\operatorname{\nabla\cdot}(\nabla\tilde{{\boldsymbol{d}}}^{T}\nabla\tilde{{\boldsymbol{d}}})-\operatorname{\nabla\cdot}\tilde{{\boldsymbol{T}}}^{L},{\boldsymbol{v}}-\tilde{{\boldsymbol{v}}}\right\rangle\mathop{}\\!\mathrm{d}s$
$\displaystyle+$
$\displaystyle\int_{0}^{t}\left\langle\nabla\left(\partial_{t}\tilde{{\boldsymbol{d}}}+(\tilde{{\boldsymbol{v}}}\cdot\nabla)\tilde{{\boldsymbol{d}}}-(\nabla\tilde{{\boldsymbol{v}}})_{\text{skw}}\tilde{{\boldsymbol{d}}}+(I-\tilde{{\boldsymbol{d}}}\otimes\tilde{{\boldsymbol{d}}})(\lambda(\nabla\tilde{{\boldsymbol{v}}})_{\text{sym}}-\Delta\tilde{{\boldsymbol{d}}})\right),\nabla({\boldsymbol{d}}-\tilde{{\boldsymbol{d}}})\right\rangle\mathop{}\\!\mathrm{d}s$
$\displaystyle\leq{}$ $\displaystyle
C\int_{0}^{t}\left(\|\tilde{{\boldsymbol{v}}}\|_{L^{\infty}(\Omega)}^{2}+\|\nabla\tilde{{\boldsymbol{v}}}\|_{L^{\infty}(\Omega)}+\|\Delta\tilde{{\boldsymbol{d}}}\|_{L^{3}(\Omega)}^{2}\right)\left[\|{\boldsymbol{v}}-\tilde{{\boldsymbol{v}}}\|_{L^{2}(\Omega)}^{2}+\|\nabla{\boldsymbol{d}}-\nabla\tilde{{\boldsymbol{d}}}\|_{L^{2}(\Omega)}^{2}\right]\text{d}s$
$\displaystyle+C\int_{0}^{t}\left(\|\nabla\tilde{{\boldsymbol{v}}}\|_{L^{3}(\Omega)}^{2}+\|\Delta\tilde{{\boldsymbol{d}}}\|_{L^{\infty}(\Omega)}+\|\Delta\tilde{{\boldsymbol{d}}}\|_{W^{1,3}(\Omega)}^{2}\right)\left[\|{\boldsymbol{v}}-\tilde{{\boldsymbol{v}}}\|_{L^{2}(\Omega)}^{2}+\|\nabla{\boldsymbol{d}}-\nabla\tilde{{\boldsymbol{d}}}\|_{L^{2}(\Omega)}^{2}\right]\text{d}s$
$\displaystyle+\int_{0}^{t}\mathcal{K}(\tilde{{\boldsymbol{v}}})\left[E-\mathcal{E}({\boldsymbol{v}},{\boldsymbol{d}})\right]\text{d}s\,.$
(17)

Gronwall’s inequality implies the relative energy inequality (10). Choosing
$(\tilde{{\boldsymbol{v}}},\tilde{{\boldsymbol{d}}})$ to be a solution,
implies that the fourth and fifth line of (10) vanish. Since all terms on the
left-hand side of (10) are non-negative, if
${\boldsymbol{v}}_{0}=\tilde{{\boldsymbol{v}}}(0)$ and
${\boldsymbol{d}}_{0}=\tilde{{\boldsymbol{d}}}(0)$ the right-hand side of (10)
is zero such that ${\boldsymbol{v}}(t)=\tilde{{\boldsymbol{v}}}(t)$ and
${\boldsymbol{d}}(t)=\tilde{{\boldsymbol{d}}}(t)$ for a.e. $t\in(0,T^{*})$,
which implies the assertion of Theorem 2.5.

∎

######  Remark (Weak-strong uniqueness extended).

It may seem strange that all calculations in the above proof are conducted on
the time interval $(0,t)$, when the definition of energy-variational solutions
2.1 is formulated on every interval $(t,s)$ for all $s<t\in(0,T)$. Indeed, the
calculations of the above proof can also be carried out for all intervals
$(s,t)$ such that we may infer inequality (17) with $0$ replaced by $s$. Thus,
the weak-strong uniqueness result can be sharpened: If there exists an
$s\in(0,T)$ such that there exists a solution
$(\tilde{{\boldsymbol{v}}},\tilde{{\boldsymbol{d}}})$ fulfilling the
regularity assumptions (9) on $(s,T)$ with
${\boldsymbol{v}}(s)=\tilde{{\boldsymbol{v}}}(s)$,
${\boldsymbol{d}}(s)=\tilde{{\boldsymbol{d}}}(s)$ and additionally
$\mathcal{E}(\tilde{{\boldsymbol{v}}}(s),\tilde{{\boldsymbol{d}}}(s))=E(s)$,
then it holds

$\displaystyle{\boldsymbol{v}}\equiv\tilde{{\boldsymbol{v}}}\,,\qquad{\boldsymbol{d}}\equiv\tilde{{\boldsymbol{d}}}\qquad\text{on
}(s,T)\,.$

###  2.3 Auxiliary Lemmata

######  Lemma 2.6.

Let ${\boldsymbol{d}}\in L^{\infty}(\Omega)\cap H^{1}(\Omega)$ with
${\boldsymbol{d}}\times\Delta{\boldsymbol{d}}\in L^{2}(\Omega)$ and
$|{\boldsymbol{d}}(x,t)|=1$ a.e. in $\Omega\times(0,T)$ such that for the
trace $|\operatorname{tr}({\boldsymbol{d}})|=1$ a.e. on $\partial\Omega$. Then
we can Identify the terms
$-{\boldsymbol{d}}\times({\boldsymbol{d}}\times\Delta{\boldsymbol{d}})=[{\boldsymbol{d}}]_{x}^{T}[{\boldsymbol{d}}]_{x}\Delta{\boldsymbol{d}}=(I-{\boldsymbol{d}}\otimes{\boldsymbol{d}})\Delta{\boldsymbol{d}}=\Delta{\boldsymbol{d}}+|\nabla{\boldsymbol{d}}|^{2}{\boldsymbol{d}}$
in $L^{2}(\Omega)$, which implies

$\displaystyle\int_{\Omega}|{\boldsymbol{d}}\times\Delta{\boldsymbol{d}}|^{2}\text{d}x=\int_{\Omega}|\Delta{\boldsymbol{d}}+|\nabla{\boldsymbol{d}}|^{2}{\boldsymbol{d}}|^{2}\text{d}x\,.$

######  Remark.

We note that the notation
$\Delta{\boldsymbol{d}}+|\nabla{\boldsymbol{d}}|^{2}{\boldsymbol{d}}$ is the
usual formulation in the Euler–Lagrange equation of harmonic maps into the
sphere [19]. The norm of this term appears as a dissipative term in the
energy-dissipation mechanism of the system (see (5) below). Formally this
implies that the solution to such a system is forced to approximately fulfill
the Euler–Lagrange equation of the Dirichlet energy minimized over all
functions with values in the sphere, in the large-time limit.

###### Proof.

For all ${\boldsymbol{\varphi}}\in H^{1}(\Omega;\mathbb{S}^{2})$, we find by
an integration-by-parts and $\nabla|{\boldsymbol{d}}|^{2}=0$ a.e. in $\Omega$
as well as a.e. on $\partial\Omega$ that

$\displaystyle\int_{\Omega}{\boldsymbol{\varphi}}\cdot\left(I-{\boldsymbol{d}}\otimes{\boldsymbol{d}}\right)\Delta{\boldsymbol{d}}\text{d}x={}$
$\displaystyle-\int_{\Omega}\nabla{\boldsymbol{\varphi}}:(I-{\boldsymbol{d}}\otimes{\boldsymbol{d}})\nabla{\boldsymbol{d}}-{\boldsymbol{\varphi}}\cdot{\boldsymbol{d}}|\nabla{\boldsymbol{d}}|^{2}-(\nabla{\boldsymbol{d}}^{T}{\boldsymbol{\varphi}})\cdot(\nabla{\boldsymbol{d}}^{T}{\boldsymbol{d}})\text{d}x$
$\displaystyle+\int_{\partial\Omega}{\boldsymbol{\varphi}}\cdot(I-{\boldsymbol{d}}\otimes{\boldsymbol{d}})\nabla{\boldsymbol{d}}\cdot{\boldsymbol{n}}\text{d}S$
$\displaystyle={}$
$\displaystyle-\int_{\Omega}\nabla{\boldsymbol{\varphi}}:\nabla{\boldsymbol{d}}-{\boldsymbol{\varphi}}\cdot|\nabla{\boldsymbol{d}}|^{2}{\boldsymbol{d}}\text{d}x-\frac{1}{2}\int_{\Omega}\nabla({\boldsymbol{\varphi}}\cdot{\boldsymbol{d}})\cdot\nabla|{\boldsymbol{d}}|^{2}\text{d}x$
$\displaystyle+\int_{\partial\Omega}{\boldsymbol{\varphi}}\cdot\nabla{\boldsymbol{d}}\cdot{\boldsymbol{n}}-\frac{1}{2}{\boldsymbol{\varphi}}\cdot{\boldsymbol{d}}{\boldsymbol{n}}\cdot\nabla|{\boldsymbol{d}}|^{2}\text{d}S$
$\displaystyle={}$
$\displaystyle\int_{\Omega}{\boldsymbol{\varphi}}\cdot(\Delta{\boldsymbol{d}}+|\nabla{\boldsymbol{d}}|^{2}{\boldsymbol{d}})\text{d}x\,,$

which implies the equivalence of the terms in $L^{2}(\Omega)$. The unit-norm
constraint implies
$(I-{\boldsymbol{d}}\otimes{\boldsymbol{d}})=(I-{\boldsymbol{d}}\otimes{\boldsymbol{d}})(I-{\boldsymbol{d}}\otimes{\boldsymbol{d}})$
such that we infer the second equality

$\int_{\Omega}|{\boldsymbol{d}}\times\Delta{\boldsymbol{d}}|^{2}\text{d}x=\int_{\Omega}\Delta{\boldsymbol{d}}\cdot(I-{\boldsymbol{d}}\otimes{\boldsymbol{d}})\Delta{\boldsymbol{d}}\text{d}x=\int_{\Omega}\left((I-{\boldsymbol{d}}\otimes{\boldsymbol{d}})\Delta{\boldsymbol{d}}\right)\cdot(I-{\boldsymbol{d}}\otimes{\boldsymbol{d}})\Delta{\boldsymbol{d}}\text{d}x\\\
=\int_{\Omega}|\Delta{\boldsymbol{d}}+|\nabla{\boldsymbol{d}}|^{2}{\boldsymbol{d}}|^{2}\text{d}x\,.$

This implies the assertion. ∎

######  Lemma 2.7 ([30, Lemma 2.4]).

Let $f\in L^{1}(0,T)$ and $g\in L^{\infty}(0,T)$ with $g\geq 0$ a.e. in
$(0,T)$. Then the two inequalities

$\displaystyle-\int_{0}^{T}\phi^{\prime}(t)g(t)\text{d}t+\int_{0}^{T}\phi(t)f(t)\text{d}t\leq
0$

for all $\phi\in{\operatorname{\mathcal{C}}}^{1}_{c}((0,T))$ with $\phi\geq 0$
for all $t\in(0,T)$ and

$\displaystyle g(t)-g(s)+\int_{s}^{t}f(\tau)\text{d}\tau\leq 0\quad\text{for
a.e.\leavevmode\nobreak\ }t,\,s\in(0,T),s\leq t\,$ (18)

are equivalent.

Moreover, from (18) we can infer that there exists a function
$h\in\operatorname{{BV}([0,T])}$ such that $h=g$ a.e. in $(0,T)$ and the
inequality (18) holds for every $s,t\in(0,T)$ with $g$ replaced by $h$.

###  2.4 Finite element spaces

For the rest of this paper we presume the following assumptions.

  1. 1.

The domain $\Omega\subset\operatorname{\mathbb{R}}^{3}$ is a convex
polyhedron.

  2. 2.

The subdivision $\mathcal{T}^{h}$ of tetrahedra of our domain
$\Omega\subset\operatorname{\mathbb{R}}^{3}$ is quasi-uniform (in the sense of
[7, 4.4.15]).

  3. 3.

All $T\in\mathcal{T}^{h}$ have at least one node that is not on the boundary
$\partial\Omega$.

  4. 4.

Every $T\in\mathcal{T}^{h}$ is equipped with a Lagrangian finite element.

In order to apply well known results for the numerical approximation of
incompressible flows, we choose ($P2$-$P1$) Taylor–Hood elements (cf. [21])
for the velocity. By $X_{h}$ and $M_{h}$ we denote the finite element spaces
for the velocity and pressure,

$\displaystyle M_{h}$ $\displaystyle\coloneqq\\{q\in
C(\overline{\Omega}):\int_{\Omega}q\mathop{}\\!\mathrm{d}x=0,\quad
q\lvert_{K}\in\mathbb{P}_{1}(K)\quad\forall K\in\mathcal{T}_{h}\\},$
$\displaystyle X_{h}$ $\displaystyle\coloneqq\\{v\in
C(\overline{\Omega}):v=0\text{ on }\partial\Omega,\quad
v\lvert_{K}\in\mathbb{P}_{2}(K)\quad\forall K\in\mathcal{T}_{h}\\},$

where $\mathbb{P}_{k}(\Omega)$ is the set of polynomials of degree $k$ or less
on the domain $\Omega$. For better readability, we denote approximately
divergence-free functions in $X_{h}$ as the space space $V_{h}$, i.e.

$\displaystyle V_{h}\coloneqq\\{v\in[X_{h}]^{3},\quad(\nabla\cdot v,q)=0\text{
for all }q\in M_{h}\\}.$ (19)

It is a well-known result (cf. [49]) that this choice of spaces fulfills the
inf-sup condition (sometimes also called LBB or Ladyzhenskaya–Babuska–Brezzi
condition) under our given assumptions (1\- 4), that is

$\displaystyle\sup_{v_{h}\in[X_{h}]^{3}}\frac{(\nabla\cdot
v_{h},q_{h})}{\left\lVert v_{h}\right\rVert_{L^{2}\left(\Omega\right)}}\geq
C\left\lVert q_{h}\right\rVert_{L^{2}\left(\Omega\right)}$

for all $q_{h}\in M_{h}$ and for some constant $C>0$. The director and
discrete Laplacian are approximated using piecewise linear functions,

$\displaystyle Z_{h}\coloneqq\\{v\in
C(\overline{\Omega}):v\lvert_{K}\in\mathbb{P}_{1}(K)\quad\forall
K\in\mathcal{T}_{h}\\}.$

The homogeneous Dirichlet boundary conditions can be implied in the space by

$\displaystyle Y_{h}\coloneqq\\{v\in C(\overline{\Omega}):v=0\text{ on
}\partial\Omega,\quad v\lvert_{K}\in\mathbb{P}_{1}(K)\quad\forall
K\in\mathcal{T}_{h}\\}.$

Our final mixed finite element space will be denoted by,

$\displaystyle U_{h}\coloneqq V_{h}\times[Y_{h}]^{3}$

######  Lemma 2.8 (Inverse estimate, cf. [7, Thm. 4.5.11]).

Under the assumptions (1-4), let $p,q\in[1,\infty]$ and $0\leq m\leq l\leq 1$.
Then there exists a constant $C$ independent of $h$, such that

$\displaystyle\left\lVert v\right\rVert_{W^{l,p}(\Omega)}\leq
Ch^{m-l+\min\\{0,3/p-3/q\\}}\left\lVert v\right\rVert_{W^{m,q}(\Omega)}$

holds for all $v\in[Z_{h}]^{3},V_{h}$.

For the finite element spaces $[Y_{h}]^{3}$ and $V_{h}$ and a function $f\in
L^{2}\left(\Omega\right)$ the standard $L^{2}$-projection, is denoted by
$\mathcal{P}_{h}{:}L^{2}(\Omega)^{3}\operatorname{\rightarrow}[Y_{h}]^{3}$ and
$\mathcal{Q}_{h}{:}{L}^{2}_{\sigma}(\Omega)\operatorname{\rightarrow}V_{h}$,
respectively.

######  Lemma 2.9 (Projection operator, cf. [50], [17, Lem. 1.131, Prop.
1.134]).

Under the previously given assumptions (1 \- 4), the $L^{2}$ projection onto
$[Y_{h}]^{3}$ is stable in $L^{p}\left(\Omega\right)$ for all $1\leq
p\leq\infty$ and $H^{1}\left(\Omega\right)$, i.e. we have

$\displaystyle\left\lVert\mathcal{P}_{h}({\boldsymbol{u}}_{1})\right\rVert_{L^{p}\left(\Omega\right)}$
$\displaystyle\leq
C\left\lVert{\boldsymbol{u}}_{1}\right\rVert_{L^{p}\left(\Omega\right)},$
$\displaystyle\left\lVert\mathcal{P}_{h}{\boldsymbol{u}}_{2}-{\boldsymbol{u}}_{2}\right\rVert_{L^{2}\left(\Omega\right)}+h\left\lVert\mathcal{P}_{h}({\boldsymbol{u}}_{2})\right\rVert_{H^{1}\left(\Omega\right)}$
$\displaystyle\leq
Ch\left\lVert{\boldsymbol{u}}_{2}\right\rVert_{H^{1}\left(\Omega\right)}$

for all ${\boldsymbol{u}}_{1}\in L^{2}\left(\Omega\right)\cap
L^{p}\left(\Omega\right)$, ${\boldsymbol{u}}_{2}\in H^{1}\left(\Omega\right)$.

For the projection onto the finite element space for the velocity $V_{h}$,
$\mathcal{Q}_{h}$ admits the error and stability estimates

$\displaystyle\begin{split}\left\lVert\mathcal{Q}_{h}({\boldsymbol{v}})-{\boldsymbol{v}}\right\rVert_{L^{2}\left(\Omega\right)}&\leq
Ch\left\lVert\nabla{\boldsymbol{v}}\right\rVert_{L^{2}\left(\Omega\right)}\\\
\|\mathcal{Q}_{h}\tilde{{\boldsymbol{v}}}-\tilde{{\boldsymbol{v}}}\|_{H^{1}(\Omega)}+h\|\mathcal{Q}_{h}\tilde{{\boldsymbol{v}}}\|_{H^{2}(\Omega)}&\leq
Ch\|\tilde{{\boldsymbol{v}}}\|_{H^{2}(\Omega)}\end{split}$ (20)

for all ${\boldsymbol{v}}\in H^{1}_{0,\sigma}(\Omega)$ and
$\tilde{{\boldsymbol{v}}}\in H^{2}(\Omega)\cap H^{1}_{0,\sigma}(\Omega)$. This
follows from [20, Lem. 4.3] and the regularity of the Stokes problem on convex
polygons (see [13], see also [42, Proposition 2] or [40, Cor. 1.8]).

###  2.5 Interpolation and mass-lumping

We define the nodal interpolation operator
$\mathcal{I}_{h}:\operatorname{\mathcal{C}}(\Omega)\to[Z_{h}]^{3}$ by

$\displaystyle\mathcal{I}_{h}(f)\coloneqq\sum_{z\in\mathcal{N}_{h}}f(z)\phi_{z},$

for all functions $f\in C(\Omega)$, where $\mathcal{N}_{h}$ is the set of all
nodes in our mesh. For continuous functions $y_{1},y_{2}\in
C(\Omega;\mathbb{R}^{3})$ we introduce the lumped inner product (cf. [44, 11])
and the according norm as

$\displaystyle\begin{aligned}
(y_{1},y_{2})_{h}&\coloneqq\int_{\Omega}\mathcal{I}_{h}(y_{1}\cdot
y_{2})\mathop{}\\!\mathrm{d}x=\sum_{z\in\mathcal{N}_{h}}y_{1}(z)\cdot
y_{2}(z)\int_{\Omega}\phi_{z}\mathop{}\\!\mathrm{d}x,\\\ \left\lVert
y_{1}\right\rVert_{h}^{2}&\coloneqq(y_{1},y_{1})_{h}.\end{aligned}$ (21)

For the interpolation operator we observe the following standard error
estimate:

######  Lemma 2.10.

Under the assumptions (1-4), let $3/2<p\leq\infty$. Then there exists a
constant $C>0$ independent of $h$ and the function $f$, such that

$\displaystyle\left\lVert
f-\mathcal{I}_{h}(f)\right\rVert_{W^{k,p}(\Omega)}\leq Ch^{2-k}\left\lvert
f\right\rvert_{W^{2,p}(\Omega)}$

holds for all $0\leq k\leq 2$ and $f\in W^{2,p}(\Omega)$.

The above theorem is standard and can for example be found in [7, Thm.
4.4.20]. There exists a generic constant $c_{L}>0$ independent of $h>0$ such
that for functions $y_{1},y_{2}\in[Z_{h}]^{3}$ the estimates

$\displaystyle\left\lVert
y_{1}\right\rVert_{L^{2}\left(\Omega\right)}\leq\left\lVert
y_{1}\right\rVert_{h}$ $\displaystyle\leq c_{L}\left\lVert
y_{1}\right\rVert_{L^{2}\left(\Omega\right)}$ (22)

hold.

######  Lemma 2.11.

Let $\frac{1}{p}+\frac{1}{q}=1$. There exists a generic constant $C>0$
independent of $h>0$ such that

$\displaystyle\begin{aligned} |(y_{1},y_{2})_{h}-(y_{1},y_{2})|&\leq
Ch\left\lVert y_{1}\right\rVert_{L^{p}\left(\Omega\right)}\left\lVert\nabla
y_{2}\right\rVert_{L^{q}\left(\Omega\right)}\\\ |(y_{h},g)_{h}-(y_{h},g)|&\leq
Ch\left\lVert y_{h}\right\rVert_{L^{p}\left(\Omega\right)}\left\lVert
g\right\rVert_{W^{2,q}\left(\Omega\right)},\,\,\end{aligned}$ (23)

for all $y_{1},y_{2},y_{h}\in[Z_{h}]^{3}$, $g\in W^{2,q}(\Omega)$ with
$q>3/2$.

###### Proof.

We proceed as in [11] and find for $y_{1}$, $y_{2}\in Z_{h}$

$\displaystyle\begin{aligned}
|(y_{1},y_{2})_{h}-(y_{1},y_{2})|&\leq\int_{\Omega}\left|\mathcal{I}_{h}(y_{1}\cdot
y_{2})-y_{1}\cdot y_{2}\right|\mathop{}\\!\mathrm{d}x\\\
&\leq\sum_{K\in\mathcal{T}^{h}}\left\lVert\mathcal{I}_{h}(y_{1}\cdot
y_{2})-y_{1}\cdot y_{2}\right\rVert_{L^{2}(K)}\left\lVert
1\right\rVert_{L^{2}(K)}\\\ &\leq
Ch^{2}\sum_{K\in\mathcal{T}^{h}}\left\lVert\nabla^{2}(y_{1}\cdot
y_{2})\right\rVert_{L^{2}(K)}\left\lVert 1\right\rVert_{L^{2}(K)}\\\ &\leq
Ch^{2}\sum_{K\in\mathcal{T}^{h}}\left\lVert\nabla y_{1}:\nabla
y_{2}\right\rVert_{L^{2}(K)}\left\lVert 1\right\rVert_{L^{2}(K)}\\\ &\leq
Ch^{2}\sum_{K\in\mathcal{T}^{h}}\left\lVert\nabla
y_{1}\right\rVert_{L^{p}(K)}\left\lVert\nabla y_{2}\right\rVert_{L^{q}(K)}\\\
&\leq Ch^{2}\left\lVert\nabla
y_{1}\right\rVert_{L^{p}\left(\Omega\right)}\left\lVert\nabla
y_{2}\right\rVert_{L^{q}\left(\Omega\right)}\\\ &\leq Ch\left\lVert
y_{1}\right\rVert_{L^{p}\left(\Omega\right)}\left\lVert\nabla
y_{2}\right\rVert_{L^{q}\left(\Omega\right)}.\end{aligned}$ (24)

The first inequality is a triangle inequality, the second one uses Hölder’s
inequality, the third bounds the interpolation error by Lemma 2.10. Since
$[Z_{h}]^{3}$ consists of piecewise affine-linear functions the second
derivative of a function in $[Z_{h}]^{3}$ vanishes. This together with
Hölder’s inequality, local and global inverse estimates (cf. 2.8) allow to
conclude. Note that the choice of exponents in the third step is necessary in
order to make sure that we can indeed apply Lemma 2.10.

In case that one function is not a finite element functions, we have to
estimate the difference with the interpolated function, i.e.,

$\displaystyle\begin{aligned}
|(\mathcal{I}_{h}y_{h},\mathcal{I}_{h}g)-(y_{h},g)|&=|(y_{h},\mathcal{I}_{h}g)-(y_{h},g)|\\\
&\leq C\left\lVert
y_{h}\right\rVert_{L^{p}\left(\Omega\right)}\left\lVert(\mathcal{I}_{h}-I)g\right\rVert_{L^{q}\left(\Omega\right)}\leq
Ch^{2}\left\lVert
y_{h}\right\rVert_{L^{p}\left(\Omega\right)}\left\lVert\nabla^{2}g\right\rVert_{L^{q}\left(\Omega\right)}\,,\end{aligned}$
(25)

and estimate the last factor on the right-hand side of (24) by

$\displaystyle\left\lVert\nabla\mathcal{I}_{h}g\right\rVert_{L^{q}\left(\Omega\right)}$
$\displaystyle\leq\left\lVert\nabla(\mathcal{I}_{h}-I)g\right\rVert_{L^{q}\left(\Omega\right)}+\left\lVert\nabla
g\right\rVert_{L^{q}\left(\Omega\right)}$ $\displaystyle\leq
Ch\left\lVert\nabla^{2}g\right\rVert_{L^{q}\left(\Omega\right)}+\left\lVert\nabla
g\right\rVert_{L^{q}\left(\Omega\right)}\leq C\left\lVert\nabla
g\right\rVert_{W^{1,q}\left(\Omega\right)}.$

∎

We introduce another interpolation operator
$\operatorname{\widetilde{\mathcal{I}}_{h}}:L^{1}(\Omega;\operatorname{\mathbb{R}}^{3})\operatorname{\rightarrow}[Z_{h}]^{3}$
by

$\operatorname{\widetilde{\mathcal{I}}_{h}}(g):=\sum_{z\in\mathcal{N}_{h}}\frac{\int_{\Omega}g(x)\phi_{z}(x)\mathop{}\\!\mathrm{d}x}{\int_{\Omega}\phi_{z}(x)\mathop{}\\!\mathrm{d}x}\phi_{z}\,.$
(26)

We observe the identity

$\displaystyle\begin{split}\left(\operatorname{\mathcal{I}_{h}}{f},g\right)&=\int_{\Omega}\sum_{z\in\mathcal{N}_{h}}f(z)\phi_{z}(x)g(x)\mathop{}\\!\mathrm{d}x\\\
&=\sum_{z\in\mathcal{N}_{h}}f(z)\int_{\Omega}\phi_{z}(x)g(x)\mathop{}\\!\mathrm{d}x=\left(f,\operatorname{\widetilde{\mathcal{I}}_{h}}g\right)_{h}\,.\end{split}$
(27)

######  Lemma 2.12.

Let $f\in\operatorname{\mathcal{C}}(\overline{{\Omega}})$. Then it holds

$\displaystyle\|f-\operatorname{\widetilde{\mathcal{I}}_{h}}f\|_{\operatorname{\mathcal{C}}(\overline{{\Omega}})}\to
0\text{ as }h\operatorname{\rightarrow}0\,.$

###### Proof.

Consider any $y\in\overline{{\Omega}}$. We obtain the pointwise estimate

$\displaystyle f(y)-\operatorname{\widetilde{\mathcal{I}}_{h}}f(y)$
$\displaystyle=\sum_{z}f(y)\phi_{z}(y)-\sum_{z}\int_{\Omega}f(x)\phi_{z}(x)\mathop{}\\!\mathrm{d}x\phi_{z}(y)\frac{1}{\int_{\Omega}\phi_{z}(x)\mathop{}\\!\mathrm{d}x}$
$\displaystyle=\sum_{z}\phi_{z}(y)\frac{\int_{\Omega}[f(y)-f(x)]\phi_{z}(x)\mathop{}\\!\mathrm{d}x}{\int_{\Omega}\phi_{z}(x)\mathop{}\\!\mathrm{d}x}$
$\displaystyle\leq\sum_{z}\phi_{z}(y)\max_{x\in supp(\phi_{z})}|f(x)-f(y)|$
$\displaystyle=\max_{z\in\mathcal{N}_{h}}\max_{x,s\in
supp(\phi_{z})}|f(x)-f(s)|,$

where the last term converges to $0$ as $h\to 0$, since $g$ is continuous.
Above, we used that $\sum_{z}\phi_{z}(y)=1$ and that support of the shape
functions $\phi_{z}$ vanishes as $h\to 0$. Since $y$ was arbitrary the
assertion follows. ∎

Further, we introduce the discrete Laplacian in a slightly different way than
usual, since we further equip it with mass lumping. Then the discrete
Laplacian $\Delta_{h}:H^{1}\left(\Omega\right)\to[Y_{h}]^{3}$ of a function
$f\in H^{1}\left(\Omega\right)$ is defined via

$\displaystyle(-\Delta_{h}f,b)_{h}$ $\displaystyle=(\nabla f,\nabla
b)\quad\text{for all }b\in Y_{h}.$ (28)

Similarly, we define an adapted projection via
$\mathcal{R}_{h}:L^{2}(\Omega;\operatorname{\mathbb{R}}^{3})\operatorname{\rightarrow}[Z_{h}]^{3}$
via

$(\mathcal{R}_{h}f,b)_{h}=(f,b)\quad\text{for all }b\in[Z_{h}]^{3}\,.$ (29)

This choice allows to infer with (22) that

$\left\lVert\mathcal{R}_{h}f\right\rVert_{h}^{2}=(\mathcal{R}_{h}f,\mathcal{R}_{h}f)_{h}=(f,\mathcal{R}_{h}f)\leq\|f\|_{L^{2}(\Omega)}\|\mathcal{R}_{h}f\|_{L^{2}(\Omega)}\leq\|f\|_{L^{2}(\Omega)}\left\lVert\mathcal{R}_{h}f\right\rVert_{h}\,,$
(30)

which implies together with (22) that
$\left\lVert\mathcal{R}_{h}f\right\rVert_{L^{2}\left(\Omega\right)}\leq\left\lVert\mathcal{R}_{h}f\right\rVert_{h}\leq\|f\|_{L^{2}(\Omega)}$.

######  Remark (Adapted discretization of convection term).

In comparison to the discretization in [4], we use the above projection
$\mathcal{R}_{h}$ in the definition of the convection term in comparison to
the $L^{2}$-projection $\mathcal{P}_{h}$. This change is needed in order to
show the convergence to the finer concept of energy-variational solutions. In
[4] only the convergence to dissipative solutions was shown, which is a bigger
class of solutions. One key ingredient for this improvement are the above
norm-equivalences.

##  3 Discrete system

In order to restrict ourselves to finite element spaces with a zero trace, we
decompose the spatial approximation of our director into its interior and its
non-homogeneous Dirichlet boundary condition, i.e.

$\displaystyle d^{j}={d^{j}_{\circ}}+\mathcal{I}_{h}{{\boldsymbol{d}}_{0}},$
(31)

with $\operatorname{tr}(d^{j}_{\circ})=0$ and ${\boldsymbol{d}}_{0}\in
H^{2}(\Omega)$.

######  Scheme 32.

Let $(v^{0},d^{0})=(\mathcal{Q}_{h}v_{0},\mathcal{I}_{h}d_{0})$. For $1\leq
j\leq J$ and $(v^{j-1},d^{j-1}_{\circ})\in U_{h}$, we want to find
$(v^{j},d^{j}_{\circ})\in U_{h}$, such that

$\displaystyle\begin{aligned}
(\boldsymbol{d}_{t}v^{j},a)+((&v^{j-1}\cdot\nabla)v^{j},a)+\frac{1}{2}([\nabla\cdot
v^{j-1}]v^{j},a)\\\ &-([\mathcal{R}_{h}{\nabla
d^{j-1}}]^{T}[d^{j-1/2}\times(d^{j-1/2}\times\Delta_{h}d^{j-1/2})],a)_{h}+(T_{L}^{j},\nabla
a)=0\,,\end{aligned}$ (32a) $\displaystyle\begin{aligned}
&(\boldsymbol{d}_{t}d^{j},c)_{h}+(d^{j-1/2}\times[\mathcal{R}_{h}{\nabla
d^{j-1}}v^{j}],d^{j-1/2}\times c)_{h}-((\nabla v^{j})_{skw}d^{j-1/2},c)_{h}\\\
&{+}\lambda(d^{j-1/2}\times[(\nabla v^{j})_{sym}d^{j-1/2}],d^{j-1/2}\times
c)_{h}{-}(d^{j-1/2}\times\Delta_{h}d^{j-1/2},d^{j-1/2}\times
c)_{h}=0,\end{aligned}$ (32b)

for all $(a,c)\in U_{h}$. Hereby we discretize the Leslie stress tensor as

$\displaystyle\begin{aligned} (T_{L}^{j},\nabla a)\coloneqq&(T_{D}^{j},\nabla
a)+\lambda(d^{j-1/2}\times[(\nabla
a)_{sym}d^{j-1/2}],d^{j-1/2}\times\Delta_{h}d^{j-1/2})_{h}\\\ &+((\nabla
a)_{skw}\Delta_{h}d^{j-1/2},d^{j-1/2})_{h},\end{aligned}$ (33)

where $T^{j}_{D}$ collects the dissipative terms of the discrete Leslie stress
tensor again, i.e.

$\displaystyle(T_{D}^{j},\nabla a)\coloneqq$
$\displaystyle(\mu_{1}+\lambda^{2})((d^{j}\cdot(\nabla
v^{j})_{sym}d^{j}),(d^{j}\cdot(\nabla a)_{sym}d^{j}))$
$\displaystyle+\mu_{4}\left((\nabla v^{j})_{sym},(\nabla
a)_{sym}\right)+(\mu_{5}+\mu_{6}-\lambda^{2})\left((\nabla
v^{j})_{sym}d^{j},(\nabla a)_{sym}d^{j}\right)\,.$

######  Remark (Features of the discrete system).

The mass lumping is applied on the director equation (LABEL:discrete_scheme_c)
to guarantee the unit-norm restriction of $d^{j}$ in every node. In order to
find a discrete energy law, the associated terms in the momentum equation
(32a) and the discrete Laplacian are mass-lumped as well. Since the gradient
of the director is piecewise constant and therefore not well-defined at the
nodes of the mesh, we apply the mass-lumped $L^{2}$ projection
$\mathcal{R}_{h}$ in the convection term. Note that this is computationally
not too expensive, since it only effects the gradient of the past iterate.

Additionally, we note that solving the discrete system amounts in introducing
two additional coupled equations, one to calculate the discrete Laplacian (28)
and one to assure that the solution $v^{j}$ fulfills the discrete vanishing
divergence condition (19).

###  3.1 A priori estimates

The proposed scheme is unconditionally solvable.

######  Proposition 3.1 (Existence of discrete solutions).

Let $k,h>0$, $j\in\mathbb{N}$ and $1\leq j\leq J=\lfloor T/k\rfloor$. Then
there exists a solution $u^{j}=(v^{j},d^{j}_{\circ})\in U_{h}$ solving Scheme
(32).

###### Proof.

We show that the map $\Pi:U_{h}\to U_{h}$ implicitly defined by the above
scheme (32), whose zero is the next iterate in terms of
$u=(v^{j},-\Delta_{h}d^{j-1/2})$, fulfills

$\displaystyle(\Pi(u),u)_{U_{h}}\geq 0$ (34)

for all $\left\lVert u\right\rVert_{U_{h}}\geq R_{j-1}$ for some constant
$R_{j-1}>0$. Where the dependency on $d^{j}$ has to be interpreted via the
discrete Greens function $\mathcal{G}_{h}$ associated to the discrete
Laplacian and given by

$\displaystyle\left(\nabla\mathcal{G}_{h}(w),\nabla
b\right)=(w,b)_{h}\quad\text{for all }b\in Y_{h}\,.$

Noting that this operator is a continuous bijection due to Lax-Milgram, we
observe that the associated mapping $\Pi$ corresponding to our discrete scheme
is actually well defined and continuous. Then the existence of solutions to
our discrete scheme (32) follows from a non-linear version of Brouwer’s fix-
point theorem (cf. [47]). In particular, we evaluate

$\displaystyle\begin{aligned}
(\Pi(v^{j},&-\Delta_{h}d^{j-1/2}),(v^{j},-\Delta_{h}d^{j-1/2}))_{U_{h}}\\\
&=\frac{1}{2}\boldsymbol{d}_{t}\left\lVert
v^{j}\right\rVert_{L^{2}\left(\Omega\right)}^{2}+\frac{1}{2}\boldsymbol{d}_{t}\left\lVert\nabla
d^{j}\right\rVert_{L^{2}\left(\Omega\right)}^{2}+\frac{k}{2}\left\lVert\boldsymbol{d}_{t}v^{j}\right\rVert_{L^{2}\left(\Omega\right)}^{2}\\\
&\quad+(\mu_{1}+\lambda^{2})\left\lVert d^{j}\cdot(\nabla
v^{j})_{sym}d^{j}\right\rVert_{L^{2}\left(\Omega\right)}+\mu_{4}\left\lVert(\nabla
v^{j})_{sym}\right\rVert_{L^{2}\left(\Omega\right)}\\\
&\quad+(\mu_{5}+\mu_{6}-\lambda^{2})\left\lVert(\nabla
v^{j})_{sym}d^{j}\right\rVert_{L^{2}\left(\Omega\right)}+\left\lVert\mathcal{I}_{h}(d^{j-1/2}\times\Delta_{h}d^{j-1/2})\right\rVert_{h}^{2},\end{aligned}$
(35)

where we used the discrete identities

$\displaystyle\begin{aligned}
&(\boldsymbol{d}_{t}v^{j},v^{j})=\frac{1}{2}\boldsymbol{d}_{t}\left\lVert
v^{j}\right\rVert_{L^{2}\left(\Omega\right)}^{2}+\frac{k}{2}\left\lVert\boldsymbol{d}_{t}v^{j}\right\rVert_{L^{2}\left(\Omega\right)}^{2}\,,\quad(-\Delta_{h}d^{j-1/2},\boldsymbol{d}_{t}d^{j})_{h}=\frac{1}{2}\boldsymbol{d}_{t}\left\lVert\nabla
d^{j}\right\rVert_{L^{2}\left(\Omega\right)}^{2}.\end{aligned}$ (36)

Equation (35) is indeed positive, if

$\displaystyle\left\lVert
v^{j}\right\rVert_{L^{2}\left(\Omega\right)}^{2}+\left\lVert\nabla
d^{j}\right\rVert_{L^{2}\left(\Omega\right)}^{2}\geq\left\lVert
v^{j-1}\right\rVert_{L^{2}\left(\Omega\right)}^{2}+\left\lVert\nabla
d^{j-1}\right\rVert_{L^{2}\left(\Omega\right)}^{2}$

holds, which only depends on $u^{j-1},k,h$. This is implied by choosing an
appropriate norm, since we have

$\displaystyle\left\lVert
v^{j}\right\rVert_{L^{2}\left(\Omega\right)}^{2}+\left\lVert\nabla
d^{j}\right\rVert_{L^{2}\left(\Omega\right)}^{2}+\left\lVert\nabla
d^{j-1}\right\rVert_{L^{2}\left(\Omega\right)}^{2}\geq\left\lVert
v^{j}\right\rVert_{L^{2}\left(\Omega\right)}^{2}+\frac{2(c_{l}h)^{2}}{3}\left\lVert\Delta_{h}d^{j-1/2}\right\rVert_{L^{2}\left(\Omega\right)}^{2}$

Note that the second term stems from a simple claculation obtained by applying
the inverse estimate:

$\displaystyle\left\lVert\Delta_{h}d^{j-1/2}\right\rVert_{L^{2}\left(\Omega\right)}^{2}\leq
c\left\lVert\Delta_{h}d^{j-1/2}\right\rVert_{h}^{2}$ $\displaystyle=c(\nabla
d^{j-1/2},\nabla(-\Delta_{h}d^{j-1/2}))$
$\displaystyle\leq(c_{l}h)^{-1}\left\lVert\nabla
d^{j-1/2}\right\rVert_{L^{2}\left(\Omega\right)}\left\lVert\Delta_{h}d^{j-1/2}\right\rVert_{L^{2}\left(\Omega\right)}.$

∎

######  Proposition 3.2 (Properties of the discrete system).

Let $u^{j}=(v^{j},d^{j}_{\circ})\in U_{h}$ be a solution of the discrete
scheme (32) for all $1\leq j\leq n$. Then the follwing energy equality holds,

$\displaystyle\small\begin{split}\frac{1}{2}\left\lVert
v^{n}\right\rVert_{L^{2}\left(\Omega\right)}^{2}&+\frac{1}{2}\left\lVert\nabla
d^{n}\right\rVert_{L^{2}\left(\Omega\right)}^{2}+\frac{k^{2}}{2}\sum_{j=1}^{n}\left\lVert\boldsymbol{d}_{t}v^{j}\right\rVert_{L^{2}\left(\Omega\right)}^{2}\\\
&+k\sum_{j=1}^{n}(\mu_{1}+\lambda^{2})\left\lVert d^{j}\cdot(\nabla
v^{j})_{sym}d^{j}\right\rVert_{L^{2}\left(\Omega\right)}^{2}+k\sum_{j=1}^{n}\mu_{4}\left\lVert(\nabla
v^{j})_{sym}\right\rVert_{L^{2}\left(\Omega\right)}^{2}\\\
&+k\sum_{j=1}^{n}(\mu_{5}+\mu_{6}-\lambda^{2})\left\lVert(\nabla
v^{j})_{sym}d^{j}\right\rVert_{L^{2}\left(\Omega\right)}^{2}+k\sum_{j=1}^{n}\left\lVert\mathcal{I}_{h}(d^{j-1/2}\times\Delta_{h}d^{j-1/2})\right\rVert_{h}^{2}\\\
&\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad=\frac{1}{2}\left\lVert
v^{0}\right\rVert_{L^{2}\left(\Omega\right)}^{2}+\frac{1}{2}\left\lVert\nabla
d^{0}\right\rVert_{L^{2}\left(\Omega\right)}^{2}\end{split}$ (37)

for all $1\leq n\leq J$. The algebraic restriction is fulfilled in every node

$\displaystyle|d^{j}(z)|=1\text{ for all }z\in\mathcal{N}_{h},$ (38)

where $\mathcal{N}_{h}$ is the set of all nodes of the mesh. Additionally, it
holds for some constant $C>0$ that

$\displaystyle
k\sum_{l=1}^{j}\left[\left\lVert\boldsymbol{d}_{t}d^{l}\right\rVert_{L^{3/2}(\Omega)}^{2}+\left\lVert\boldsymbol{d}_{t}v^{l}\right\rVert_{(H^{2}(\Omega)\cap
H^{1}_{0,\sigma}(\Omega))^{*}}^{2}+\left\lVert\nabla
v^{l}\right\rVert_{L^{2}\left(\Omega\right)}^{2}\right]\leq
C\left(\mathcal{E}(v_{0},d_{0})+1\right)\,.$ (39)

###### Proof.

The discrete energy estimate (37) follows simply from iterating (35). In order
to show (38), we adapt the approach in [4]. Therefore, the director equation
(LABEL:discrete_scheme_c) is tested with $d^{j-1/2}(\hat{z})\phi_{\hat{z}}$,
where ${\hat{z}}\in\mathcal{N}_{h}$ is a node of the mesh. All terms except
for the one including the time-derivative of the director vanish due to the
mass lumping and the skew-symmetry. Then, we are left with

$\displaystyle
0=(\boldsymbol{d}_{t}d^{j},d^{j-1/2}({\hat{z}})\phi_{\hat{z}})_{h}=$
$\displaystyle\sum_{z\in\mathcal{N}_{h}}\boldsymbol{d}_{t}d^{j-1/2}(z)\cdot
d^{j-1/2}({\hat{z}})\phi_{\hat{z}}\int_{\Omega}\phi_{z}\mathop{}\\!\mathrm{d}x$
$\displaystyle=$
$\displaystyle\frac{1}{2k}(|d^{j}({\hat{z}})|^{2}-|d^{j-1}({\hat{z}})|^{2})\int_{\Omega}\phi_{\hat{z}}\mathop{}\\!\mathrm{d}x,$

which yields the result for all interior nodes $\hat{z}$. At the boundary the
unit-norm constraint is fulfilled by the assumption on ${\boldsymbol{d}}_{0}$,
our inhomogeneous Dirichlet boundary condition. Korn’s inequality yields the
estimate for the last addend in (39).

Since $\boldsymbol{d}_{t}d^{j}\in Y_{h}$, the definition of the
$L^{2}$-projection yields for a smooth test function
${\boldsymbol{\phi}}\in\operatorname{\mathcal{C}}^{\infty}_{c}(\Omega)$ that

$\displaystyle(\boldsymbol{d}_{t}d^{j},{\boldsymbol{\phi}})=(\boldsymbol{d}_{t}d^{j},\mathcal{P}_{h}({\boldsymbol{\phi}})).$

Using a duality argument, we are able to estimate the time derivative of the
director in a $L^{p}$-norm by

$\displaystyle\begin{aligned}
\left\lVert\boldsymbol{d}_{t}d^{j}\right\rVert_{L^{3/2}\left(\Omega\right)}&=\sup_{{\boldsymbol{\phi}}\in
L^{3}\left(\Omega\right)}\left\lvert\frac{(\boldsymbol{d}_{t}d^{j},\mathcal{P}_{h}({\boldsymbol{\phi}}))}{\left\lVert\phi\right\rVert_{L^{3}\left(\Omega\right)}}\right\rvert\\\
&\leq\sup_{{\boldsymbol{\phi}}\in
L^{3}\left(\Omega\right)}\left\lvert\frac{(\boldsymbol{d}_{t}d^{j},\mathcal{P}_{h}({\boldsymbol{\phi}}))_{h}}{\left\lVert\phi\right\rVert_{L^{3}\left(\Omega\right)}}\right\rvert+\sup_{{\boldsymbol{\phi}}\in
L^{3}\left(\Omega\right)}\left\lvert\frac{(\boldsymbol{d}_{t}d^{j},\mathcal{P}_{h}({\boldsymbol{\phi}}))-(\boldsymbol{d}_{t}d^{j},\mathcal{P}_{h}({\boldsymbol{\phi}}))_{h}}{\left\lVert\phi\right\rVert_{L^{3}\left(\Omega\right)}}\right\rvert.\end{aligned}$
(40)

The first term can easily be estimated using our discrete director equation
(LABEL:discrete_scheme_c) for $c\in Y_{h}$:

$\displaystyle\begin{aligned}
\left\lvert(\boldsymbol{d}_{t}d^{j},c)_{h}\right\rvert\leq&C\left\lVert
d^{j-1/2}\right\rVert_{L^{\infty}\left(\Omega\right)}^{2}\left\lVert
v^{j}\right\rVert_{L^{6}\left(\Omega\right)}\left\lVert\mathcal{R}_{h}(\nabla
d^{j-1})\right\rVert_{L^{2}\left(\Omega\right)}\left\lVert
c\right\rVert_{L^{3}\left(\Omega\right)}\\\ &+C\left\lVert
d^{j-1/2}\times\Delta_{h}d^{j-1/2}\right\rVert_{h}\left\lVert
d^{j-1/2}\right\rVert_{L^{\infty}\left(\Omega\right)}\left\lVert
c\right\rVert_{L^{2}\left(\Omega\right)}\\\ &+C\left\lVert
d^{j-1/2}\right\rVert_{L^{\infty}\left(\Omega\right)}^{3}\left\lVert\nabla
v^{j}\right\rVert_{L^{2}\left(\Omega\right)}\left\lVert
c\right\rVert_{L^{2}\left(\Omega\right)}\\\ \leq&C\left(\left\lVert\nabla
v^{j}\right\rVert_{L^{2}\left(\Omega\right)}+\left\lVert\mathcal{I}_{h}\left(d^{j-1/2}\times\Delta_{h}d^{j-1/2}\right)\right\rVert_{h}\right)\left\lVert
c\right\rVert_{L^{3}\left(\Omega\right)}.\end{aligned}$ (41)

Choosing $c=\mathcal{P}_{h}({\boldsymbol{\phi}})$ for ${\boldsymbol{\phi}}\in
H^{1}(\Omega)$ and using the embedding
$H^{1}\left(\Omega\right)\hookrightarrow L^{6}\left(\Omega\right)$, the
stability of the projection, Lemma 2.9, and the a priori estimates of
Proposition 3.2 allow to obtain a bound for the right-hand side. The second
term in (40) describes the error introduced by the mass-lumping. It must be
estimated locally as in (24), i.e. by additionally applying the inverse
estimate one obtains

$\displaystyle\begin{aligned}
{\int_{\Omega}(I-\mathcal{I}_{h})(\boldsymbol{d}_{t}d^{j}\cdot\mathcal{P}_{h}({\boldsymbol{\phi}}))\mathop{}\\!\mathrm{d}x}&\leq\sum_{K\in\mathcal{T}_{h}}{\int_{K}(I-\mathcal{I}_{h})(\boldsymbol{d}_{t}d^{j}\cdot\mathcal{P}_{h}({\boldsymbol{\phi}}))\mathop{}\\!\mathrm{d}x}\\\
&\leq
C\sum_{K\in\mathcal{T}_{h}}\left\lVert\boldsymbol{d}_{t}d^{j}\right\rVert_{L^{2}(K)}\left\lVert\mathcal{P}_{h}({\boldsymbol{\phi}})\right\rVert_{L^{3}(K)}\left\lVert
1\right\rVert_{L^{6}(K)}\\\ &\leq
Ch^{1/2}\sum_{K\in\mathcal{T}_{h}}\left\lVert\boldsymbol{d}_{t}d^{j}\right\rVert_{L^{2}(K)}\left\lVert\mathcal{P}_{h}({\boldsymbol{\phi}})\right\rVert_{L^{3}(K)}\,.\end{aligned}$
(42)

For the time derivative of the director $\boldsymbol{d}_{t}d^{j}$ in the
$L^{2}$-norm, we observe locally by choosing
$c=\sum_{{\boldsymbol{z}}\in\mathcal{N}_{h}\cap
K}\boldsymbol{d}_{t}d_{j}({\boldsymbol{z}})\phi_{z}$ denotes the in a local
version of (41) that

$\displaystyle\left\lVert\boldsymbol{d}_{t}d^{j}\right\rVert_{L^{2}(K)}^{2}$
$\displaystyle\leq
C\int_{K}\mathcal{I}_{h}\left(|\boldsymbol{d}_{t}d^{j}|^{2}\right)\mathop{}\\!\mathrm{d}x$
$\displaystyle\leq C\Big{(}\left\lVert
v^{j}\right\rVert_{L^{6}(K)}\left\lVert\mathcal{R}_{h}(\nabla
d^{j-1})\right\rVert_{L^{2}(K)}\left\lVert\boldsymbol{d}_{t}d^{j}\right\rVert_{L^{3}(K)}$
$\displaystyle\quad+\left(\int_{K}\mathcal{I}_{h}\left[d^{j-1/2}\times\Delta_{h}d^{j-1/2}\right]^{2}\mathop{}\\!\mathrm{d}x+\left\lVert\nabla
v^{j}\right\rVert_{L^{2}(K)}\right)\left\lVert\boldsymbol{d}_{t}d^{j}\right\rVert_{L^{2}(K)}\Big{)}\,.$

With the inverse estimate we can also bound the time derivative in the
$L^{3}$-norm by

$\displaystyle\left\lVert\boldsymbol{d}_{t}d^{j}\right\rVert_{L^{3}(K)}\leq
Ch^{-1/2}\left\lVert\boldsymbol{d}_{t}d^{j}\right\rVert_{L^{2}(K)}.$

Putting the above together and reinserting it into (42) yields

$\displaystyle\int_{\Omega}$
$\displaystyle(I-\mathcal{I}_{h})(\boldsymbol{d}_{t}d^{j}\cdot\mathcal{P}_{h}({\boldsymbol{\phi}}))\mathop{}\\!\mathrm{d}x$
$\displaystyle\leq C\left\lVert
v^{j}\right\rVert_{L^{6}(\Omega)}\left\lVert\mathcal{R}_{h}(\nabla
d^{j-1})\right\rVert_{L^{2}(\Omega)}\left\lVert\mathcal{P}_{h}({\boldsymbol{\phi}})\right\rVert_{L^{3}(\Omega)}$
$\displaystyle\quad+C\left(\left\lVert\nabla
v^{j}\right\rVert_{L^{2}\left(\Omega\right)}+\left\lVert
d^{j-1/2}\times\Delta_{h}d^{j-1/2}\right\rVert_{h}\right)\left\lVert\mathcal{P}_{h}({\boldsymbol{\phi}})\right\rVert_{L^{2}(\Omega)}\,.$

Multiplying by $k$ and summing up leads to the estimate of
$\boldsymbol{d}_{t}d^{j}$, due to the a priori estimates (37).
For the velocity, we proceed as in [15]: We test the time derivative of the
velocity with a smooth solenoidal test function ${\boldsymbol{\phi}}$, and
apply the projection $\mathcal{Q}_{h}$ such that
$\langle\boldsymbol{d}_{t}v^{j},{\boldsymbol{\phi}}\rangle=(\boldsymbol{d}_{t}v^{j},\mathcal{Q}_{h}{\boldsymbol{\phi}})$
and equation (32a),

$\displaystyle\left\lvert\langle\boldsymbol{d}_{t}v^{j},{\boldsymbol{\phi}}\rangle\right\rvert$
$\displaystyle=|((v^{j-1}\cdot\nabla)v^{j},\mathcal{Q}_{h}({\boldsymbol{\phi}}))+\frac{1}{2}([\nabla\cdot
v^{j-1}]v^{j},\mathcal{Q}_{h}({\boldsymbol{\phi}}))$
$\displaystyle\quad-([\mathcal{R}_{h}{\nabla
d^{j-1}}]^{T}[d^{j-1/2}\times(d^{j-1/2}\times\Delta_{h}d^{j-1/2})],\mathcal{Q}_{h}({\boldsymbol{\phi}}))_{h}+(T_{L}^{j},\nabla\mathcal{Q}_{h}{\boldsymbol{\phi}})|.$

The convective terms can be estimated by a constant using (37) and

$\left\lvert((v^{j-1}\cdot\nabla)v^{j},\mathcal{Q}_{h}({\boldsymbol{\phi}}))+\frac{1}{2}([\nabla\cdot
v^{j-1}]v^{j},\mathcal{Q}_{h}({\boldsymbol{\phi}}))\right\rvert\\\
\leq(\left\lVert
v^{j-1}\right\rVert_{L^{2}\left(\Omega\right)}\left\lVert\nabla
v^{j}\right\rVert_{L^{2}\left(\Omega\right)}+\left\lVert
v^{j}\right\rVert_{L^{2}\left(\Omega\right)}\left\lVert\nabla
v^{j-1}\right\rVert_{L^{2}\left(\Omega\right)})\left\lVert\mathcal{Q}_{h}({\boldsymbol{\phi}})\right\rVert_{L^{\infty}\left(\Omega\right)}.$

For the Leslie stress tensor we infer the estimate

$\displaystyle\left\lvert(T_{L}^{j}:\nabla\mathcal{Q}_{h}{\boldsymbol{\phi}})\right\rvert\leq$
$\displaystyle\left\lvert(T_{D}^{j}:\nabla\mathcal{Q}_{h}{\boldsymbol{\phi}})\right\rvert+|(d^{j-1/2},(\nabla\mathcal{Q}_{h}{\boldsymbol{\phi}})_{skw}\Delta_{h}d^{j-1/2})_{h}|$
$\displaystyle\quad\quad+\lambda|(d^{j-1/2}\times[(\nabla\mathcal{Q}_{h}{\boldsymbol{\phi}})_{sym}d^{j-1/2}],d^{j-1/2}\times\Delta_{h}d^{j-1/2})_{h}|$
$\displaystyle\leq$ $\displaystyle C\Big{[}\left\lVert(\nabla
v^{j})_{sym}\right\rVert_{L^{2}\left(\Omega\right)}\left\lVert(\nabla\mathcal{Q}_{h}{\boldsymbol{\phi}})_{sym}\right\rVert_{L^{2}\left(\Omega\right)}$
$\displaystyle\quad\quad+\left\lVert\mathcal{I}_{h}\left(d^{j-1/2}\times\Delta_{h}d^{j-1/2}\right)\right\rVert_{L^{2}\left(\Omega\right)}\left\lVert\nabla\mathcal{Q}_{h}{\boldsymbol{\phi}}\right\rVert_{L^{2}\left(\Omega\right)}\Big{]}.$

Note that we hereby reformulated the term including the skew-symmetric part of
the velocity’s gradient on the algebraic level, i.e.

$\displaystyle(\Delta_{h}d^{j-1/2},(\nabla\mathcal{Q}_{h}{\boldsymbol{\phi}})_{skw}d^{j-1/2})_{h}=(d^{j-1/2}\times\Delta_{h}d^{j-1/2},d^{j-1/2}\times[(\nabla\mathcal{Q}_{h}{\boldsymbol{\phi}})_{skw}]d^{j-1/2})_{h},$

since $[d]_{x}^{T}[d]_{x}(A)_{\text{skw}}d=(I-d\otimes
d)(A)_{skw}d=(A)_{skw}d$ holds. The Ericksen stress tensor can be estimated in
a similar fashion

$\displaystyle\left\lvert([\mathcal{R}_{h}{\nabla
d^{j-1}}]^{T}[d^{j-1/2}\times(d^{j-1/2}\times\Delta_{h}d^{j-1/2})],\mathcal{Q}_{h}{\boldsymbol{\phi}})_{h}\right\rvert\quad\quad\quad\quad\quad$
$\displaystyle\leq C\left\lVert\mathcal{R}_{h}{\nabla
d^{j-1}}\right\rVert_{L^{2}\left(\Omega\right)}\left\lVert\mathcal{I}_{h}\left(d^{j-1/2}\times\Delta_{h}d^{j-1/2}\right)\right\rVert_{L^{2}\left(\Omega\right)}\left\lVert\mathcal{Q}_{h}({\boldsymbol{\phi}})\right\rVert_{L^{\infty}\left(\Omega\right)}.$

Summing over $j$, using the embedding $H^{2}(\Omega)\hookrightarrow
L^{\infty}(\Omega)$ and the stability of the projection onto $V_{h}$ due to
(20) yields the assertion. ∎

######  Proposition 3.3.

Let $u^{j}=(v^{j},d^{j}_{\circ})\in U_{h}$ be a solution to scheme (32). Then
the discrete energy-variational inequality

$\displaystyle\begin{split}\boldsymbol{d}_{t}E^{j}&+\frac{k}{2}\left\lVert\boldsymbol{d}_{t}v^{j}\right\rVert_{L^{2}\left(\Omega\right)}^{2}\\\
&+\left(T_{D}^{j};(\nabla
v^{j})_{\text{sym}}-(\nabla\mathcal{Q}_{h}\tilde{{\boldsymbol{v}}})_{\text{sym}}\right)+\left\lVert
d^{j-1/2}\times\Delta_{h}d^{j-1/2}\right\rVert_{h}^{2}\\\
&-(\boldsymbol{d}_{t}v^{j},\mathcal{Q}_{h}\tilde{{\boldsymbol{v}}})-\left((v^{j-1}\cdot\nabla)v^{j},\mathcal{Q}_{h}\tilde{{\boldsymbol{v}}}\right)-\frac{1}{2}\left((\operatorname{\nabla\cdot}v^{j-1})v^{j},\mathcal{Q}_{h}\tilde{{\boldsymbol{v}}}\right)\\\
&-\left([\mathcal{R}_{h}{\nabla
d^{j-1}}]^{T}[d^{j-1/2}\times(d^{j-1/2}\times\Delta_{h}d^{j-1/2})],\mathcal{Q}_{h}\tilde{{\boldsymbol{v}}}\right)_{h}\\\
&-\lambda(d^{j-1/2}\times[(\nabla\mathcal{Q}_{h}\tilde{{\boldsymbol{v}}})_{sym}d^{j-1/2}],d^{j-1/2}\times\Delta_{h}d^{j-1/2})_{h}\\\
&+((\nabla\mathcal{Q}_{h}\tilde{{\boldsymbol{v}}})_{skw}\Delta_{h}d^{j-1/2},d^{j-1/2})_{h}\\\
&+\mathcal{K}(\mathcal{Q}_{h}\tilde{{\boldsymbol{v}}})\left(\frac{1}{2}\|v^{j-1}\|_{L^{2}(\Omega)}^{2}+\frac{1}{2}\|\nabla
d^{j-1}\|_{L^{2}(\Omega)}^{2}-E^{j-1}\right)=0\end{split}$ (43)

with the energy $E^{j}$ and the regularity weight $\mathcal{K}$ given by

$\displaystyle E^{j}:=\frac{1}{2}\left\lVert
v^{j}\right\rVert_{L^{2}\left(\Omega\right)}^{2}+\frac{1}{2}\left\lVert\nabla
d^{j}\right\rVert_{L^{2}\left(\Omega\right)}^{2}\qquad\mathcal{K}(\tilde{{\boldsymbol{v}}}):={}$
$\displaystyle\frac{1}{2}\|\tilde{{\boldsymbol{v}}}\|_{L^{\infty}(\Omega)}^{2}\,.$

holds for all $\tilde{{\boldsymbol{v}}}\in
H^{1}(0,T;(H^{1}_{0,\sigma}(\Omega))^{*})\cap
L^{2}(0,T;H^{1}_{0,\sigma}(\Omega))\cap\operatorname{\mathcal{C}}(\overline{{\Omega}}\times[0,T])$.

###### Proof.

Adding (32a) tested with $-\mathcal{Q}_{h}\tilde{{\boldsymbol{v}}}^{j}$ to
(37), we find the variational-energy inequality

$\displaystyle\begin{split}\frac{1}{2k}\|v^{j}&\|_{L^{2}(\Omega)}^{2}+\frac{1}{2k}\|\nabla
d^{j}\|_{L^{2}(\Omega)}^{2}+\frac{k}{2}\left\lVert\boldsymbol{d}_{t}v^{j}\right\rVert_{L^{2}\left(\Omega\right)}^{2}\\\
&+\left(T_{D}^{j};(\nabla
v^{j})_{\text{sym}}-(\nabla\mathcal{Q}_{h}\tilde{{\boldsymbol{v}}})_{\text{sym}}\right)+\left\lVert
d^{j-1/2}\times\Delta_{h}d^{j-1/2}\right\rVert_{h}^{2}\\\
&-(\boldsymbol{d}_{t}v^{j},\mathcal{Q}_{h}\tilde{{\boldsymbol{v}}})-\left((v^{j-1}\cdot\nabla)v^{j},\mathcal{Q}_{h}\tilde{{\boldsymbol{v}}}\right)-\frac{1}{2}\left((\operatorname{\nabla\cdot}v^{j-1})v^{j},\mathcal{Q}_{h}\tilde{{\boldsymbol{v}}}\right)\\\
&-\left([\mathcal{R}_{h}{\nabla
d^{j-1}}]^{T}[d^{j-1/2}\times(d^{j-1/2}\times\Delta_{h}d^{j-1/2})],\mathcal{Q}_{h}\tilde{{\boldsymbol{v}}}\right)_{h}\\\
&-\lambda(d^{j-1/2}\times[(\nabla\mathcal{Q}_{h}\tilde{{\boldsymbol{v}}})_{sym}d^{j-1/2}],d^{j-1/2}\times\Delta_{h}d^{j-1/2})_{h}\\\
&+((\nabla\mathcal{Q}_{h}\tilde{{\boldsymbol{v}}})_{skw}\Delta_{h}d^{j-1/2},d^{j-1/2})_{h}=\frac{1}{2k}\left\lVert
v^{j-1}\right\rVert_{L^{2}\left(\Omega\right)}^{2}+\frac{1}{2k}\left\lVert\nabla
d^{j-1}\right\rVert_{L^{2}\left(\Omega\right)}^{2}.\end{split}$

By defining the variable $E^{j}$ and adding
$\mathcal{K}(\mathcal{Q}_{h}\tilde{{\boldsymbol{v}}})\left(\frac{1}{2}\|v^{j-1}\|_{L^{2}(\Omega)}^{2}+\frac{1}{2}\|\nabla
d^{j-1}\|_{L^{2}(\Omega)}^{2}-E^{j-1}\right)=0$ implies (43). ∎

###  3.2 Converging subsequences

Let $\\{u^{j}\\}_{j=0,...,J}$ be a sequence of measurable functions in space.
We define their constant and linear interpolates in time as follows,

$\displaystyle\overline{u}^{k}_{h}(t)$ $\displaystyle\coloneqq u^{j},$
$\displaystyle\quad\underline{\overline{u}}^{k}_{h}(t)$
$\displaystyle\coloneqq\frac{1}{2}(\overline{u}^{k}_{h}(t)+\underline{u}^{k}_{h}(t)),$
$\displaystyle\underline{u}^{k}_{h}(t)$ $\displaystyle\coloneqq u^{j-1},$
$\displaystyle\quad{u}^{k}_{h}(t)$
$\displaystyle\coloneqq\frac{u^{j}-u^{j-1}}{k}(t-jk)+u^{j},$

for $(j-1)k<t\leq jk$. This yields the standard relation between the
continuous and discrete time derivative,
$\partial_{t}u^{k}_{h}(t)=\frac{u^{j}-u^{j-1}}{k}=\boldsymbol{d}_{t}u^{j}$.
Analogously we define the interpolates for a smooth test function
$\tilde{u}\in C([0,T];\mathbb{Y})$ for $(j-1)k<t\leq jk$ by,

$\displaystyle\overline{\tilde{u}}^{k}(t)\coloneqq\tilde{u}(jk),\,\,\,\underline{\tilde{u}}^{k}(t)\coloneqq\tilde{u}((j-1)k),\,\,\,\hat{\tilde{u}}^{k}(t)\coloneqq\frac{\tilde{u}(jk)-\tilde{u}((j-1)k)}{k}(t-jk)+\tilde{u}(jk)\,.$

The a priori estimates 37 allow us to extract converging subsequences which we
do not relabel and $({\boldsymbol{v}},{\boldsymbol{d}},E)$ fulfilling (4) such
that

$\displaystyle\begin{split}\overline{v}^{k}_{h},\underline{v}^{k}_{h},v_{h}^{k}&\overset{\ast}{\rightharpoonup}{\boldsymbol{v}}\text{
in }L^{\infty}(0,T;{L}^{2}_{\sigma}(\Omega)),\\\
\overline{v}^{k}_{h},\underline{v}^{k}_{h},v_{h}^{k}&\rightharpoonup{\boldsymbol{v}}\text{
in }L^{2}(0,T;H^{1}_{0}(\Omega)),\\\
\partial_{t}v_{k}^{h}&\rightharpoonup\partial_{t}{\boldsymbol{v}}\text{ in
}L^{2}(0,T;(H^{2}(\Omega)\cap H^{1}_{0,\sigma}(\Omega))^{*}),\\\
\overline{d}^{k}_{h},\underline{d}^{k}_{h},d_{h}^{k}&\overset{\ast}{\rightharpoonup}{\boldsymbol{d}}\text{
in }L^{\infty}(0,T;L^{\infty}(\Omega))\cap L^{\infty}(0,T;H^{1}(\Omega)),\\\
\partial_{t}d_{k}^{h}&\rightharpoonup\partial_{t}{\boldsymbol{d}}\text{ in
}L^{2}(0,T;L^{3/2}(\Omega)),\end{split}$ (44)
$\displaystyle\begin{split}\overline{d}^{k}_{h}\cdot(\nabla\overline{v}^{k}_{h})_{sym}\overline{d}^{k}_{h}&\rightharpoonup{\boldsymbol{d}}\cdot(\nabla{\boldsymbol{v}})_{sym}{\boldsymbol{d}}\text{
in }L^{2}(0,T;L^{2}\left(\Omega\right)),\\\
(\nabla\overline{v}^{k}_{h})_{sym}\overline{d}^{k}_{h}&\rightharpoonup(\nabla{\boldsymbol{v}})_{sym}{\boldsymbol{d}}\text{
in }L^{2}(0,T;L^{2}\left(\Omega\right)),\\\
\mathcal{I}_{h}(\underline{\overline{d}}^{k}_{h}\times\Delta_{h}\overline{\underline{d}}_{k}^{h})&\rightharpoonup{\boldsymbol{d}}\times\Delta{\boldsymbol{d}}\text{
in }L^{2}(0,T;L^{2}\left(\Omega\right)),\\\ \underline{E}^{k}_{h}(t)&\to
E(t)\text{ for all }t\in[0,T]\,,\end{split}$ (45)

where the last convergence follows from Helly’s selection theorem (cf. [8, Ex.
8.3]). We can refine the convergence for the director by using the following
compact embeddings, referred to as Aubin–Lions–Simon Lemma [48],

$\displaystyle\begin{aligned} \\{{\boldsymbol{d}}\in
L^{2}(0,T;H^{1}\left(\Omega\right)):\partial_{t}{\boldsymbol{d}}\in
L^{2}(0,T;L^{3/2}\left(\Omega\right))\\}&\stackrel{{\scriptstyle
c}}{{\hookrightarrow}}L^{2}(0,T;L^{2}\left(\Omega\right)),\\\
\\{{\boldsymbol{d}}\in
L^{\infty}(0,T;H^{1}\left(\Omega\right)):\partial_{t}{\boldsymbol{d}}\in
L^{2}(0,T;L^{3/2}\left(\Omega\right))\\}&\stackrel{{\scriptstyle
c}}{{\hookrightarrow}}C^{0}(0,T;L^{2}\left(\Omega\right)),\end{aligned}$ (46)

which yield strong convergence for the director

$\displaystyle d_{k}^{h}\to{\boldsymbol{d}}\text{ in
}L^{2}(0,T;L^{2}(\Omega)).$

A standard interpolation inequality (cf. [6, p. 192 ff.]) and the uniform
bound in $L^{\infty}(0,T;L^{\infty}(\Omega))$ yields strong convergence,

$d_{k}^{h}\to{\boldsymbol{d}}\text{ in }L^{p}(0,T;L^{p}(\Omega))\text{ for all
}p\in[1,\infty).$ (47)

For the term ${\boldsymbol{d}}\times\Delta{\boldsymbol{d}}$, note that due to
the above convergences and (23), we have

$\displaystyle\int_{0}^{T}$
$\displaystyle(\mathcal{I}_{h}(\underline{\overline{d}}^{k}_{h}\times\Delta_{h}\overline{\underline{d}}_{k}^{h}),\phi)\mathop{}\\!\mathrm{d}t$
$\displaystyle=\int_{0}^{T}(\mathcal{I}_{h}(\underline{\overline{d}}^{k}_{h}\times\Delta_{h}\overline{\underline{d}}_{k}^{h}),\phi)-(\underline{\overline{d}}^{k}_{h}\times\Delta_{h}\overline{\underline{d}}_{k}^{h},\phi)_{h}\mathop{}\\!\mathrm{d}t+\int_{0}^{T}(\underline{\overline{d}}^{k}_{h}\times\Delta_{h}\overline{\underline{d}}_{k}^{h},\phi)_{h}\mathop{}\\!\mathrm{d}t$
$\displaystyle\leq
Ch\left\lVert\mathcal{I}_{h}(\underline{\overline{d}}^{k}_{h}\times\Delta_{h}\overline{\underline{d}}_{k}^{h})\right\rVert_{L^{2}(0,T;L^{2}(\Omega))}\left\lVert\nabla\phi\right\rVert_{L^{2}(0,T;H^{1}(\Omega))}-\int_{0}^{T}(\nabla\overline{\underline{d}}_{k}^{h},\nabla\phi\times\underline{\overline{d}}^{k}_{h})\mathop{}\\!\mathrm{d}t$

for $\phi\in C_{c}^{\infty}([0,T]\times\Omega;\mathbb{R}^{3})$. Now note that
the first term vanishes as $h\to 0$ and the last term converges to

$\displaystyle-\int_{0}^{T}(\nabla{\boldsymbol{d}},\nabla\phi\times{\boldsymbol{d}})\mathop{}\\!\mathrm{d}t=-\int_{0}^{T}({\boldsymbol{d}}\times\nabla{\boldsymbol{d}},\nabla\phi)\mathop{}\\!\mathrm{d}t=\int_{0}^{T}({\boldsymbol{d}}\times\Delta{\boldsymbol{d}},\phi)\mathop{}\\!\mathrm{d}t$

where the last equality,
$\nabla\cdot({\boldsymbol{d}}\times\nabla{\boldsymbol{d}})={\boldsymbol{d}}\times\Delta{\boldsymbol{d}}$
, has to be understood in the distributional sense. The fact that all
interpolates in (44) converge against the same limit can be confirmed as in
[43]. Consider examplary for the velocity,

$\displaystyle\begin{aligned} \left\lVert
v^{k}_{h}-\underline{v}^{k}_{h}\right\rVert_{L^{2}(0,T;L^{2}\left(\Omega\right))}^{2}&\leq\sum_{j=1}^{J}\left\lVert
v^{j}-v^{j-1}\right\rVert_{L^{2}\left(\Omega\right)}^{2}\int_{j(k-1)}^{jk}\left(\frac{t-jk}{k}\right)^{2}\mathop{}\\!\mathrm{d}t\\\
&\leq
k\sum_{j=1}^{J}\frac{k^{2}}{3}\left\lVert\boldsymbol{d}_{t}v^{j}\right\rVert_{L^{2}\left(\Omega\right)}^{2},\end{aligned}$

which vanishes for $k\to 0$, since the term on the right-hand side is bounded
by our energy estimate (37). In order to deduce that the velocity is
solenoidal, we use standard arguments (e.g. as in [12]). We test the
divergence of the velocity with a test function ${\boldsymbol{\phi}}\in
C_{c}^{\infty}(\Omega\times(0,T])$ with
$\int_{\Omega}{\boldsymbol{\phi}}\mathop{}\\!\mathrm{d}x=0$ and use our
discrete divergence-zero condition to obtain

$\displaystyle\int_{0}^{T}(\nabla\cdot{\boldsymbol{v}}(s),{\boldsymbol{\phi}})\mathop{}\\!\mathrm{d}s=\int_{0}^{T}(\nabla\cdot({\boldsymbol{v}}(s)-v^{h}_{k}(s)),{\boldsymbol{\phi}})+(\nabla\cdot
v^{h}_{k}(s),{\boldsymbol{\phi}}-\Pi_{h}{\boldsymbol{\phi}})\mathop{}\\!\mathrm{d}s,$

where $\Pi_{h}$ is the $L^{2}$-projection onto $M_{h}$. The first term
vanishes due to the weak convergence of the velocity in (44). The second
summand vanishes for $h\to 0$ simply due to the uniform boundedness
$v^{h}_{k}$ and a standard approximation property of the space $M_{h}$ (see
e.g. [7, Lemma 12.4.3]), by

$\displaystyle\int_{0}^{T}(\nabla\cdot
v^{h}_{k},{\boldsymbol{\phi}}-\Pi_{h}({\boldsymbol{\phi}}))\mathop{}\\!\mathrm{d}s$
$\displaystyle\leq\left\lVert
v^{h}_{k}\right\rVert_{L^{2}(0,T;H^{1}_{0}(\Omega))}\left\lVert{\boldsymbol{\phi}}-\Pi_{h}{\boldsymbol{\phi}}\right\rVert_{L^{2}(0,T;L^{2}(\Omega))}$
$\displaystyle\leq
Ch^{2}\left\lVert\nabla^{2}{\boldsymbol{\phi}}\right\rVert_{L^{2}(0,T;L^{2}(\Omega))}.$

The convergence for our non-homogeneous Dirichlet boundary conditions of the
director, follows from the trace theorem (cf. [2, Thm. 6.3.10]) and the
interpolation error estimate (2.10),

$\displaystyle\begin{aligned}
\left\lVert\operatorname{tr}({\boldsymbol{d}}_{0})-\operatorname{tr}(d^{h}_{k}(t))\right\rVert_{L^{2}(\partial\Omega)}=&\left\lVert\operatorname{tr}({\boldsymbol{d}}_{0})-\operatorname{tr}(\mathcal{I}_{h}{\boldsymbol{d}}_{0})-\operatorname{tr}(d^{h}_{k,\circ}(t))\right\rVert_{L^{2}(\partial\Omega)}\\\
=&\left\lVert\operatorname{tr}((I-\mathcal{I}_{h}){\boldsymbol{d}}_{0})\right\rVert_{L^{2}(\partial\Omega)}\\\
\leq&C\left\lVert(I-\mathcal{I}_{h}){\boldsymbol{d}}_{0}\right\rVert_{H^{1}(\Omega)}\leq
Ch\left\lVert{\boldsymbol{d}}_{0}\right\rVert_{H^{2}\left(\Omega\right)}\,,\end{aligned}$
(48)

for every $t\in[0,T]$. Analogously, we observe for the initial condition of
the director

$\displaystyle\left\lVert
d^{h}_{k}(0)-{\boldsymbol{d}}_{0}\right\rVert_{H^{1}\left(\Omega\right)}=\left\lVert\mathcal{I}_{h}({\boldsymbol{d}}_{0})-{\boldsymbol{d}}_{0}\right\rVert_{H^{1}\left(\Omega\right)}\leq
Ch\left\lVert{\boldsymbol{d}}_{0}\right\rVert_{H^{2}\left(\Omega\right)}.$
(49)

For the initial condition of the velocity, we can derive an analogous estimate
– however in the $L^{2}$ norm instead – by making use of (20), since we
approximated the initial condition with the respective projection operator,
$v^{0}=\mathcal{Q}_{h}(v_{0})$.

###  3.3 Unit-norm restriction and convergence for the director equations

From (38), we can infer that the unit-norm restriction $|d(t,x)|=1$ holds at
every node $z\in\mathcal{N}_{h}$. In the limit as $h\to 0$, this restriction
will be fulfilled almost everywhere on $\Omega$, since
$\mathcal{I}_{h}(|d^{j}|^{2}-1)=0$ holds. By subsequently applying the
interpolation error bound (cf.Lemma 2.10) and the inverse estimate (cf. Lemma
2.8), one obtains

$\displaystyle\left\lVert|d^{j}|^{2}-1\right\rVert_{L^{2}\left(\Omega\right)}$
$\displaystyle=\left\lVert(I-\mathcal{I}_{h})(|d^{j}|^{2}-1)\right\rVert_{L^{2}\left(\Omega\right)}$
$\displaystyle=\sqrt{\sum_{K\in\mathcal{T}_{h}}\left\lVert(I-\mathcal{I}_{h})(|d^{j}|^{2}-1)\right\rVert_{L^{2}(K)}^{2}}$
$\displaystyle\leq
Ch^{2}\sqrt{\sum_{K\in\mathcal{T}_{h}}\left\lVert\nabla^{2}|d^{j}|^{2}\right\rVert_{L^{2}(K)}^{2}}$
$\displaystyle\leq
Ch\left\lVert\nabla|d^{j}|^{2}\right\rVert_{L^{2}\left(\Omega\right)}$
$\displaystyle\leq Ch\left\lVert\nabla
d^{j}\right\rVert_{L^{2}\left(\Omega\right)}\left\lVert
d^{j}\right\rVert_{L^{\infty}\left(\Omega\right)},$

where the right-hand side is bounded uniformly due to our a priori energy
estimate (37). This gives us the result of Proposition 3.4.

######  Proposition 3.4.

Let $d^{k}_{h}$ be the approximate solution in the sense of the discrete
scheme (32) that converges to a limit as in (44). Then the norm of $d^{k}_{h}$
converges to $1$ in a linear order, i.e.

$\displaystyle\left\lVert|d^{k}_{h}(t)|^{2}-1\right\rVert_{L^{2}\left(\Omega\right)}\in\mathcal{O}(h).$

In order to prove the convergence of the approximate solution to
(LABEL:discrete_scheme_c) to a solution of (6), we consider equation
(LABEL:discrete_scheme_c) tested with the projection
$\mathcal{P}_{h}{\boldsymbol{c}}$ of a smooth test function
${\boldsymbol{c}}\in C^{\infty}_{c}(\Omega\times(0,T))$ and observe, for
instance by Lemma 2.11 that

$\displaystyle\int_{0}^{T}(\boldsymbol{d}_{t}d^{k}_{h},\mathcal{P}_{h}{{\boldsymbol{c}}})_{h}-(\boldsymbol{d}_{t}d^{k}_{h},\mathcal{P}_{h}{{\boldsymbol{c}}})\mathop{}\\!\mathrm{d}t$
$\displaystyle\leq
Ch\left\lVert\boldsymbol{d}_{t}d^{k}_{h}\right\rVert_{L^{2}(0,T;L^{2}(\Omega))}\left\lVert\mathcal{P}_{h}{\boldsymbol{c}}\right\rVert_{L^{2}(0,T;H^{1}(\Omega))}$
$\displaystyle\leq
Ch^{1/2}\left\lVert\boldsymbol{d}_{t}d^{k}_{h}\right\rVert_{L^{2}(0,T;L^{3/2}(\Omega))}\left\lVert\mathcal{P}_{h}{\boldsymbol{c}}\right\rVert_{L^{2}(0,T;H^{1}(\Omega))}.$

This term vanishes as $h\operatorname{\rightarrow}0$ since the last factor is
bounded for smooth ${\boldsymbol{c}}$ due to 2.9. For the convection term, we
infer due to the definition 29

$\displaystyle\left(\overline{\underline{d}}^{k}_{h}\times[\mathcal{R}_{h}{\nabla\underline{d}^{k}_{h}}\overline{v}^{k}_{h}],\overline{\underline{d}}^{k}_{h}\times\mathcal{P}_{h}{{\boldsymbol{c}}}\right)_{h}$
$\displaystyle=-\left(\mathcal{R}_{h}{\nabla\underline{d}^{k}_{h}},\mathcal{I}_{h}\left[\overline{\underline{d}}^{k}_{h}\times\left(\overline{\underline{d}}^{k}_{h}\times\mathcal{P}_{h}{{\boldsymbol{c}}}\right)\otimes\overline{v}^{k}_{h}\right]\right)_{h}$
$\displaystyle=-\left({\nabla\underline{d}^{k}_{h}},\left[\overline{\underline{d}}^{k}_{h}\times\left(\overline{\underline{d}}^{k}_{h}\times\mathcal{P}_{h}{{\boldsymbol{c}}}\right)\otimes\overline{v}^{k}_{h}\right]\right)-\left({\nabla\underline{d}^{k}_{h}},(\mathcal{I}_{h}-I)\left[\overline{\underline{d}}^{k}_{h}\times\left(\overline{\underline{d}}^{k}_{h}\times\mathcal{P}_{h}{{\boldsymbol{c}}}\right)\otimes\overline{v}^{k}_{h}\right]\right)\,.$

For the last term, we infer

$\displaystyle\left|\left({\nabla\underline{d}^{k}_{h}},(\mathcal{I}_{h}-I)\left[\overline{\underline{d}}^{k}_{h}\times\left(\overline{\underline{d}}^{k}_{h}\times\mathcal{P}_{h}{{\boldsymbol{c}}}\right)\otimes\overline{v}^{k}_{h}\right]\right)\right|$
$\displaystyle\leq{}\left\lVert{\nabla\underline{d}^{k}_{h}}\right\rVert_{L^{2}\left(\Omega\right)}\left\lVert(\mathcal{I}_{h}-I)\left[\overline{\underline{d}}^{k}_{h}\times\left(\overline{\underline{d}}^{k}_{h}\times\mathcal{P}_{h}{{\boldsymbol{c}}}\right)\otimes\overline{v}^{k}_{h}\right]\right\rVert_{L^{2}\left(\Omega\right)}$

such that we find

$\displaystyle\left\lVert(\mathcal{I}_{h}-I)\left[\overline{\underline{d}}^{k}_{h}\times\left(\overline{\underline{d}}^{k}_{h}\times\mathcal{P}_{h}{{\boldsymbol{c}}}\right)\otimes\overline{v}^{k}_{h}\right]\right\rVert_{L^{2}\left(\Omega\right)}^{2}$
$\displaystyle\leq{}Ch^{4}\sum_{K\in\mathcal{T}_{h}}\left\lVert\nabla^{2}\left[\overline{\underline{d}}^{k}_{h}\times\left(\overline{\underline{d}}^{k}_{h}\times\mathcal{P}_{h}{{\boldsymbol{c}}}\right)\otimes\overline{v}^{k}_{h}\right]\right\rVert_{L^{2}(K)}^{2}$
$\displaystyle\leq{}Ch^{4}\sum_{K\in\mathcal{T}_{h}}\left(\left\lVert\nabla\overline{\underline{d}}^{k}_{h}\times\left(\nabla\overline{\underline{d}}^{k}_{h}\times\mathcal{P}_{h}{{\boldsymbol{c}}}\right)\otimes\overline{v}^{k}_{h}+\nabla\overline{\underline{d}}^{k}_{h}\times\left(\overline{\underline{d}}^{k}_{h}\times\nabla\mathcal{P}_{h}{{\boldsymbol{c}}}\right)\otimes\overline{v}^{k}_{h}\right\rVert_{L^{2}(K)}\right)^{2}$
$\displaystyle\quad+Ch^{4}\sum_{K\in\mathcal{T}_{h}}\left(\left\lVert\nabla\overline{\underline{d}}^{k}_{h}\times\left(\overline{\underline{d}}^{k}_{h}\times\mathcal{P}_{h}{{\boldsymbol{c}}}\right)\otimes\nabla\overline{v}^{k}_{h}+\overline{\underline{d}}^{k}_{h}\times\left(\nabla\overline{\underline{d}}^{k}_{h}\times\nabla\mathcal{P}_{h}{{\boldsymbol{c}}}\right)\otimes\overline{v}^{k}_{h}\right\rVert_{L^{2}(K)}\right)^{2}$
$\displaystyle\quad+Ch^{4}\sum_{K\in\mathcal{T}_{h}}\left(\left\lVert\overline{\underline{d}}^{k}_{h}\times\left(\nabla\overline{\underline{d}}^{k}_{h}\times\mathcal{P}_{h}{{\boldsymbol{c}}}\right)\otimes\nabla\overline{v}^{k}_{h}+\overline{\underline{d}}^{k}_{h}\times\left(\overline{\underline{d}}^{k}_{h}\times\nabla\mathcal{P}_{h}{{\boldsymbol{c}}}\right)\otimes\nabla\overline{v}^{k}_{h}\right\rVert_{L^{2}(K)}\right)^{2}$
$\displaystyle\quad+Ch^{4}\sum_{K\in\mathcal{T}_{h}}\left\lVert\overline{\underline{d}}^{k}_{h}\times\left(\overline{\underline{d}}^{k}_{h}\times\mathcal{P}_{h}{{\boldsymbol{c}}}\right)\otimes\nabla^{2}\overline{v}^{k}_{h}\right\rVert_{L^{2}(K)}^{2}$
$\displaystyle\leq
C\left[h\left\lVert\nabla\overline{\underline{d}}^{k}_{h}\right\rVert_{L^{2}\left(\Omega\right)}^{2}\left\lVert\overline{v}^{k}_{h}\right\rVert_{L^{6}\left(\Omega\right)}^{2}{+}h^{2}\left\lVert\overline{\underline{d}}^{k}_{h}\right\rVert_{L^{\infty}\left(\Omega\right)}^{2}\left\lVert\nabla\overline{v}^{k}_{h}\right\rVert_{L^{2}\left(\Omega\right)}^{2}\right]\left\lVert\overline{\underline{d}}^{k}_{h}\right\rVert_{L^{\infty}\left(\Omega\right)}^{2}\left\lVert\mathcal{P}_{h}{\boldsymbol{c}}\right\rVert_{L^{\infty}\left(\Omega\right)}^{2}\,.$

The first inequality is due to Lemma 2.10, in the second inequality, we
calculate the derivative, where we used that the second derivative of $P1$
finite elements vanishes on every tetrahedron. In order to infer the last
inequality, we use Hölder’s inequality and inverse estimates.

In order to pass from the mass-lumped inner product to the standard
$L^{2}(\Omega)$ inner product, other occurring terms have to be estimated in
the same fashion, but locally. Regarding the second term for example, we
observe

$\displaystyle\int_{\Omega}$
$\displaystyle(I-\mathcal{I}_{h})(\overline{\underline{d}}^{k}_{h}\times[(\nabla\overline{v}^{k}_{h})_{\text{sym}}\overline{\underline{d}}^{k}_{h}]\cdot\overline{\underline{d}}^{k}_{h}\times\mathcal{P}_{h}{{\boldsymbol{c}}})\mathop{}\\!\mathrm{d}x$
$\displaystyle\leq
Ch^{2}\sum_{K\in\mathcal{T}_{h}}\left\lVert\nabla^{2}(\overline{\underline{d}}^{k}_{h}\times[(\nabla\overline{v}^{k}_{h})_{\text{sym}}\overline{\underline{d}}^{k}_{h}]\cdot\overline{\underline{d}}^{k}_{h}\times\mathcal{P}_{h}{{\boldsymbol{c}}})\right\rVert_{L^{2}(K)}\left\lVert
1\right\rVert_{L^{2}(K)}$ $\displaystyle\leq
Ch^{2+3/2}\sum_{K\in\mathcal{T}_{h}}\sum_{i,j=1}^{3}\left\lVert\partial_{x_{i}}\overline{\underline{d}}^{k}_{h}\times[\partial_{x_{j}}(\nabla\overline{v}^{k}_{h})_{\text{sym}}\overline{\underline{d}}^{k}_{h}]\cdot(\overline{\underline{d}}^{k}_{h}\times\mathcal{P}_{h}{{\boldsymbol{c}}})\right\rVert_{L^{2}(K)}$
$\displaystyle\quad+Ch^{2+3/2}\sum_{K\in\mathcal{T}_{h}}\sum_{i,j=1}^{3}\left\lVert\partial_{x_{i}}\overline{\underline{d}}^{k}_{h}\times[(\nabla\overline{v}^{k}_{h})_{\text{sym}}\partial_{x_{j}}\overline{\underline{d}}^{k}_{h}]\cdot(\overline{\underline{d}}^{k}_{h}\times\mathcal{P}_{h}{{\boldsymbol{c}}})\right\rVert_{L^{2}(K)}$
$\displaystyle\quad+Ch^{2+3/2}\sum_{K\in\mathcal{T}_{h}}\sum_{i,j=1}^{3}\left\lVert\partial_{x_{i}}\overline{\underline{d}}^{k}_{h}\times[(\nabla\overline{v}^{k}_{h})_{\text{sym}}\overline{\underline{d}}^{k}_{h}]\cdot\partial_{x_{j}}(\overline{\underline{d}}^{k}_{h}\times\mathcal{P}_{h}{{\boldsymbol{c}}})\right\rVert_{L^{2}(K)}$
$\displaystyle\quad+Ch^{2+3/2}\sum_{K\in\mathcal{T}_{h}}\sum_{i,j=1}^{3}\left\lVert\overline{\underline{d}}^{k}_{h}\times[\partial_{x_{i}}(\nabla\overline{v}^{k}_{h})_{\text{sym}}\partial_{x_{j}}\overline{\underline{d}}^{k}_{h}]\cdot(\overline{\underline{d}}^{k}_{h}\times\mathcal{P}_{h}{{\boldsymbol{c}}})\right\rVert_{L^{2}(K)}$
$\displaystyle\quad+Ch^{2+3/2}\sum_{K\in\mathcal{T}_{h}}\sum_{i,j=1}^{3}\left\lVert\overline{\underline{d}}^{k}_{h}\times[\partial_{x_{i}}(\nabla\overline{v}^{k}_{h})_{\text{sym}}\overline{\underline{d}}^{k}_{h}]\cdot\partial_{x_{j}}(\overline{\underline{d}}^{k}_{h}\times\mathcal{P}_{h}{{\boldsymbol{c}}})\right\rVert_{L^{2}(K)}$
$\displaystyle\quad+Ch^{2+3/2}\sum_{K\in\mathcal{T}_{h}}\sum_{i,j=1}^{3}\left\lVert\overline{\underline{d}}^{k}_{h}\times[(\nabla\overline{v}^{k}_{h})_{\text{sym}}\partial_{x_{i}}\overline{\underline{d}}^{k}_{h}]\cdot\partial_{x_{j}}(\overline{\underline{d}}^{k}_{h}\times\mathcal{P}_{h}{{\boldsymbol{c}}})\right\rVert_{L^{2}(K)}$
$\displaystyle\quad+Ch^{2+3/2}\sum_{K\in\mathcal{T}_{h}}\sum_{i,j=1}^{3}\left\lVert\overline{\underline{d}}^{k}_{h}\times[(\nabla\overline{v}^{k}_{h})_{\text{sym}}\overline{\underline{d}}^{k}_{h}]\cdot\partial^{2}_{x_{i}x_{j}}(\overline{\underline{d}}^{k}_{h}\times\mathcal{P}_{h}{{\boldsymbol{c}}})\right\rVert_{L^{2}(K)}\,,$
where we used that the second derivatives of the Lagrangian finite element
functions vanish on the set $K$. Using Hölder’s inequality and inverse
estimates for all appearing terms, we infer $\displaystyle\int_{\Omega}$
$\displaystyle(I-\mathcal{I}_{h})(\overline{\underline{d}}^{k}_{h}\times[(\nabla\overline{v}^{k}_{h})_{\text{sym}}\overline{\underline{d}}^{k}_{h}]\cdot\overline{\underline{d}}^{k}_{h}\times\mathcal{P}_{h}{{\boldsymbol{c}}})\mathop{}\\!\mathrm{d}x$
$\displaystyle\leq
Ch\left\lVert(\nabla\overline{v}^{k}_{h})_{\text{sym}}\right\rVert_{L^{2}\left(\Omega\right)}\left\lVert\overline{\underline{d}}^{k}_{h}\right\rVert_{L^{\infty}\left(\Omega\right)}^{2}\left\lVert\nabla\overline{\underline{d}}^{k}_{h}\right\rVert_{L^{2}\left(\Omega\right)}\left\lVert\mathcal{P}_{h}{\boldsymbol{c}}\right\rVert_{L^{\infty}\left(\Omega\right)}$
$\displaystyle\quad+Ch\left\lVert(\nabla\overline{v}^{k}_{h})_{\text{sym}}\right\rVert_{L^{2}\left(\Omega\right)}\left\lVert\overline{\underline{d}}^{k}_{h}\right\rVert_{L^{\infty}\left(\Omega\right)}^{2}\left\lVert\nabla\mathcal{P}_{h}{\boldsymbol{c}}\right\rVert_{L^{2}\left(\Omega\right)}\left\lVert\overline{\underline{d}}^{k}_{h}\right\rVert_{L^{\,}\left(\Omega\right)}.$

The right-hand side is bounded due to the $H^{1}$-stability of the $L^{2}$
projection (cf. Lemma 2.9). The estimates for the difference of the mass-
lumped and $L^{2}(\Omega)$ inner product for the other terms follow
analogously. Thus, it remains to pass to the limit in
(LABEL:discrete_scheme_c) with $L^{2}(\Omega)$\- inner products instead of the
lumped-inner products. Note that by standard approximation arguments
$\mathcal{P}_{h}{{\boldsymbol{c}}}\to{\boldsymbol{c}}$ in
$L^{3}(\Omega\times(0,T))$ for
${\boldsymbol{c}}\in\operatorname{\mathcal{C}}^{\infty}_{c}(\Omega\times(0,T))$,
cf. Lemma 2.9. The weak convergences of (44) and (45) allow to pass to the
limit in all occurring terms, and we infer

$\displaystyle\int_{0}^{T}\int_{\Omega}\partial_{t}{\boldsymbol{d}}\cdot{\boldsymbol{c}}-(\nabla{\boldsymbol{v}})_{\text{skw}}{\boldsymbol{d}}\cdot{\boldsymbol{c}}+({\boldsymbol{d}}\times(\nabla{\boldsymbol{d}}{\boldsymbol{v}}+\lambda(\nabla{\boldsymbol{v}})_{\text{sym}}{\boldsymbol{d}}-\Delta{\boldsymbol{d}}))\cdot({\boldsymbol{d}}\times{\boldsymbol{c}})\mathop{}\\!\mathrm{d}x\mathop{}\\!\mathrm{d}t=0$

for all ${\boldsymbol{c}}\in L^{2}(0,T;L^{3}(\Omega))$ by density arguments.
The regularity of ${\boldsymbol{d}}$ suffices to infer
$2{\boldsymbol{d}}\cdot(\nabla{\boldsymbol{d}}{\boldsymbol{v}})=({\boldsymbol{v}}\cdot\nabla)|{\boldsymbol{d}}|^{2}$
such that
$({\boldsymbol{d}}\times\nabla{\boldsymbol{d}}{\boldsymbol{v}})\cdot({\boldsymbol{d}}\times{\boldsymbol{c}})=({\boldsymbol{v}}\cdot\nabla){\boldsymbol{d}}\cdot{\boldsymbol{c}}$.
Since the relation holds for all ${\boldsymbol{c}}\in
L^{2}(0,T;L^{3}(\Omega))$ it also holds a.e. in $\Omega\times(0,T)$, which is
nothing else than (6).

###  3.4 Convergence to the energy-variational formulation

Multiplying the discrete energy-variational formulation (43) by
${{\boldsymbol{\phi}}}^{j}$ for a ${\boldsymbol{\phi}}\in C_{c}^{\infty}(0,T)$
applying an discrete integration as well as a discrete integration-by-parts,
we find

$\displaystyle\begin{split}-{}&\int_{0}^{T}(\partial_{t}\hat{{\boldsymbol{\phi}}}^{k})\underline{{E}}^{k}\text{d}t+\int_{0}^{T}\overline{{{\boldsymbol{\phi}}}}^{k}\frac{1}{2}\left\lVert\partial_{t}\hat{v}^{k}_{h}\right\rVert_{L^{2}\left(\Omega\right)}^{2}\text{d}t\\\
&+\int_{0}^{T}\overline{{{\boldsymbol{\phi}}}}^{k}\left[\left(\overline{{T_{D}}}^{k}_{h};(\nabla\overline{{v}}^{k}_{h})_{\text{sym}}-(\nabla\mathcal{Q}_{h}\overline{{\tilde{{\boldsymbol{v}}}}}^{k})_{\text{sym}}\right)+\left\lVert\overline{{\underline{{d}}}}^{k}_{h}\times\Delta_{h}\overline{{\underline{{d}}}}^{k}_{h}\right\rVert^{2}_{h}\right]\text{d}t\\\
&+\int_{0}^{T}\left[(\underline{{v}}^{k}_{h},\partial_{t}\mathcal{Q}_{h}\overline{{\tilde{{\boldsymbol{v}}}}}^{k}\overline{{{\boldsymbol{\phi}}}}^{k}+\mathcal{Q}_{h}\underline{{\tilde{{\boldsymbol{v}}}}}\partial_{t}\overline{{{\boldsymbol{\phi}}}}^{k})-\overline{{{\boldsymbol{\phi}}}}^{k}\left((\underline{{v}}^{k}_{h}\cdot\nabla)\overline{{v}}^{k}_{h},\mathcal{Q}_{h}\overline{{\tilde{{\boldsymbol{v}}}}}^{k}\right)\right]\text{d}t\\\
&-\int_{0}^{T}\overline{{{\boldsymbol{\phi}}}}^{k}\left[\frac{1}{2}\left((\operatorname{\nabla\cdot}\underline{{v}}^{k}_{h})\overline{{v}}^{k}_{h},\mathcal{Q}_{h}\overline{{\tilde{{\boldsymbol{v}}}}}^{k}\right)+\left([\mathcal{R}_{h}{\nabla\underline{{d}}^{k}}]^{T}[\overline{{\underline{{d}}}}^{k}_{h}\times(\overline{{\underline{{d}}}}^{k}_{h}\times\Delta_{h}\overline{{\underline{{d}}}}^{k}_{h})],\mathcal{Q}_{h}\overline{{\tilde{{\boldsymbol{v}}}}}^{k}\right)_{h}\right]\text{d}t\\\
&-\int_{0}^{T}\overline{{{\boldsymbol{\phi}}}}^{k}\left[\lambda(\overline{{\underline{{d}}}}^{k}_{h}\times[(\nabla\mathcal{Q}_{h}\overline{{\tilde{{\boldsymbol{v}}}}}^{k})_{sym}\overline{{\underline{{d}}}}^{k}_{h}],\overline{{\underline{{d}}}}^{k}_{h}\times\Delta_{h}\overline{{\underline{{d}}}}^{k}_{h})_{h}-((\nabla\mathcal{Q}_{h}\overline{{\tilde{{\boldsymbol{v}}}}}^{k})_{skw}\Delta_{h}\overline{{\underline{{d}}}}^{k}_{h},\overline{{\underline{{d}}}}^{k}_{h})_{h}\right]\text{d}t\\\