\section{Introduction}
\label{sec:intro}
For more than 130 years (cf. Lehmann \cite{lehmann1889krystalle}) liquid crystals inspire many scientists in a lot of ways.
Their field of applications ranges from  photonics and optics (e.g. liquid crystal displays (LCD)) over materials science (cf. \cite{palffy2007lc}) up to biophysics, for example to describe the cytoskeleton (cf. \cite{ahmadi2006activeLC, kierfeld2008active}). In the latter, there is often an activity mechanism involved which drives the system out of equilibrium and which will be discussed in more detail later.
Liquid crystals can be considered as a system which has mesomorphic states between an ordinary liquid and a crystal. In particular, such systems consist of rodlike polymers which interact pairwise and are characterized by long-range orientational order and not or only partially by positional order.
In the 1980s, Doi developed a model for rodlike molecules which describes the properties of liquid crystal polymers in a solvent (cf. \cite[Chapter 8]{doi1986theory}, \cite{doi1981molec_dynamics, degennes1993LC}). This model consists of a Smoluchowksi equation describing the evolution of a pseudo probability density of particles $\psi(\mathbf{x}, \mathbf{m}, t)$
with center of mass in $\mathbf{x} \in \Omega$ and orientation $\mathbf{m} \in {\SMat^{d-1}}$ at time $t \geq 0$, where ${\SMat^{d-1}}$ denotes the unit sphere and $\Omega \subset \mathbb{R}^d$ is a bounded domain for $d \in \{2, 3\}$.
If the interaction of the rodlike polymers is very strong or if the concentration of them is very high, the rods line up with each other and build a nematic phase. This phase transition problem was first described by Onsager (cf. \cite{onsager1949colloidal, zhang2007review}).
In the homogeneous case, that is in the absence of effects of fluid dynamics, this model is often referred to as Doi-Onsager model whereas the model is referred to as Doi-Hess model in the inhomogeneous case.
In the past, these models have been the subject of many mathematical studies (cf. \cite{constantin2005nonlFP, constantin2007smolNast, constantin2008smolNaSt2d, liu2005classification, luo2004structure, otto2008velocityGrad, zhou2005equilibria}).\\
However, the models for inhomogeneous flow did not take distortional elasticity into account (cf. \cite{e_zhang_2006_kinetic, zhang2008existenceNewModel}). There have been many attempts to extend the theory for inhomogeneous flow to include distortional elasticity on a microscopic level (cf. \cite{marruci1991nematics, tsuji1998nematicLC, rey1998LCrheology, wang2002nematicLC}). Unfortunately, these models are all phenomenological in nature and hence, they consist of many unknown parameters which are very difficult to determine in experiments.
Moreover, they often lack consistency with other theories and among themselves (cf. \cite{e_zhang_2006_kinetic}).\\
Then, E and Zhang (cf. \cite{e_zhang_2006_kinetic}) extended the Doi model for homogeneous flow of rodlike liquid crystalline polymers to inhomogeneous flow and introduced a nonlocal intermolecular potential which results in an extra term in the form of an elastic body force to take into account distortional elasticity.\\
In this paper, we consider unflexible rodlike polymers with constant length $l > 0$, diameter $b > 0$, ratio $r:= \frac{b}{l} \ll 1$ (see Figure \ref{fig:rod}), and a mean number density $\ndens > 0$ given by
\begin{align*}
\ndens := \frac{1}{\abs{\Omega}} \int_{\OmSph} \psi \, \text{d}\mathscr{H}^{d-1} \, \text{d}\xflow.
\end{align*}
In particular, we restrict ourselves to solutions with a concentrated regime described by $\nu b l^2 > 1$ such as liquid crystals (cf. \cite[Chapter 10.1]{doi1986theory}).
\begin{figure}[!htbp]
    \begin{tikzpicture}[scale=0.75]

        \def\mx{0}
        \def\my{0}
        \def3/4{0.5}
        \def6*\b{6*3/4}
        \def1.4142{1.4142}
        \def0.70711{0.70711}

        \def\llx{\mx - 0.70711*3/4/2 - 0.70711*6*\b}
        \def\lly{\my + 0.70711*3/4/2 - 0.70711*6*\b}

        \def\ulx{\mx - 0.70711*3/4/2 + 0.70711*6*\b}
        \def\uly{\mx + 0.70711*3/4/2 + 0.70711*6*\b}

        \def\lrx{\mx + 0.70711*3/4/2 - 0.70711*6*\b}
        \def\lry{\my - 0.70711*3/4/2 - 0.70711*6*\b}

        \def\urx{\mx + 0.70711*3/4/2 + 0.70711*6*\b}
        \def\ury{\my - 0.70711*3/4/2 + 0.70711*6*\b}

        \def\radellipsex{3/4/2}
        \def\radellipsey{3/4/4}

        \def\disttorod{3/4*2/3}

        \draw (\llx, \lly) -- (\ulx, \uly);
        \draw (\lrx, \lry) -- (\urx, \ury);

        \node[circle,fill=black,inner sep=0pt,minimum size=5pt,label=below left:{$\mathbf{x}$}] (a) at (\mx, \my) {};

        \draw[rotate around={-45:(\llx/2 + \lrx/2, \lly/2 + \lry/2)}] (\llx/2 + \lrx/2, \lly/2 + \lry/2) ellipse ({\radellipsex} and {\radellipsey});
        \draw[rotate around={-45:(\urx, \ury)}] (\urx, \ury) arc(0:180: {\radellipsex} and {\radellipsey});

        \draw[arrows={angle 90 - angle 90}] (\llx - \disttorod, \lly + \disttorod) -- (\ulx - \disttorod, \uly + \disttorod);
        \draw[arrows={angle 90 - angle 90}] (\ulx + \disttorod, \uly + \disttorod) -- (\urx + \disttorod, \ury + \disttorod);

        \node (nameL) at (\mx - 0.70711*3/4 - \disttorod - 3/4/2, \my + 0.70711*3/4 + \disttorod + 3/4/2) {$l$};
        \node (nameB) at (\ulx/2 + \urx/2 + \disttorod + 3/4/2, \uly/2 + \ury/2 + \disttorod + 3/4/2) {$b$};
    \end{tikzpicture}
    \caption{Rodlike polymer with length $l$, diameter $b$, and center of mass $\mathbf{x}$.}
    \label{fig:rod}
\end{figure}
The model of E and Zhang in dimensional form (cf. \cite{e_zhang_2006_kinetic}) coupling a Smoluchowski equation to a Navier-Stokes type equation reads as follows for a given time $T > 0$:
\begin{subequations}
	\label{intro:eq:dimensional}
	\begin{align}
	\begin{split}
	\label{intro:eq:psi:dimensional:srot}
	&\partial_t \psi + \Div_\xflow(\mathbf{u} \psi) - \frac{1}{ k_B \mathcal{T}} \Div_\xflow\lrcb{ \psi\lrcb{ D_\parallel \mflow\otimes\mflow + D_\perp (\mathbb{I} - \mflow\otimes\mflow) }\nabla_\xflow \overline{\mu}_\psi } \\
	&\qquad - \frac{D_r}{ k_B \mathcal{T}} \sRot \cdot \lrcb{ \psi \sRot \overline{\mu}_\psi }  = -  \sRot \cdot \lrcb{ \psi \mathbf{m} \times (\nabla_\xflow \mathbf{u} \mathbf{m})  } \quad \text{ in } {\Omega\times\Sph}\times(0,T],
	\end{split}\\
	\begin{split}
	\label{intro:eq:u:dimensional}
	&\rho \lrcb{ \partial_t \mathbf{u} + (\mathbf{u} \cdot \nabla_\xflow) \mathbf{u} } + \nabla_\xflow p = - \int_{\Sph} \psi \nabla_\xflow \overline{\mu}_\psi \, \text{d}\mathscr{H}^{d-1} + \Div_\xflow(\overline{\mathbb{T}}) \quad \text{ in } \Omega \times (0,T].
	\end{split}
	\end{align}
\end{subequations}
At this point let us refer to Section \ref{sec:model} for details about this model.
As the degrees of freedom of the microscopic equation are quite high, there have been many attempts to derive macroscopic equations from the Doi model (a good overview over the different models from a mathematical point of view can be found in \cite{emmrich2018LC, ball2017LC}).
In particular, the interest in the Ericksen-Leslie equations (cf. \cite{ericksen1961convervationLC, ericksen1991orientationLC, leslie1968LC}), the Q-Tensor equations (cf. \cite{degennes1993LC, majumdar2010LC, beris1994thermodynamics}), and the mathematical analysis of these models (cf. \cite{lin1995dissipativeLC, lin1996dynamicLC, lin2000ericksenLeslie, cavaterra2013LC, ball2010LC, paicu2011QTensor, paicu2012QTensor, huang2015LC, abels2014QTensor}) has been growing significantly over the last years.
Kuzuu and Doi's (cf. \cite{kuzuu1983LC}) formal derivation of the Ericksen-Leslie equation from the Doi-Onsager equation in the homogeneous case has been extended to the inhomogeneous case from systematic asymptotic analysis of \eqref{intro:eq:dimensional} (cf. E and Zhang \cite{e_zhang_2006_kinetic}). Moreover, Wang et al. (cf. \cite{wang2015deborahLimit}) established a rigorous derivation of the Ericksen-Leslie equation from \eqref{intro:eq:dimensional}. Furthermore, Han et al. (cf. \cite{han2015microToMacro}) derived a new Q-Tensor model from \eqref{intro:eq:dimensional}, which obeys the energy dissipation law, by applying the Bingham-closure (cf. \cite{bingham1974}). Therefore, as the model \eqref{intro:eq:dimensional} is not just a fundamental model itself but a starting point for the derivation of many equations, which find a lot of interest these days, it is important to study it analytically.

Chen and Liu (\cite{chen_liu_2013_entropy_solution}) proved the existence of global-in-time weak solutions to \eqref{intro:eq:dimensional} including active stresses induced by active particles (see Section \ref{sec:model} for more details regarding activity) in two and three space dimensions in the case that the friction constant is zero, that is $\xi_r = 0$ (see \eqref{eq:viscous_stress_p}) and that the interaction potential equals zero, that is $\pot[\psi] \equiv 0$ (see \eqref{def:pot}).
Zhang and Zhang (cf. \cite{zhang2008existenceNewModel}) considered the model \eqref{intro:eq:dimensional}, periodic domains, and high regularity assumptions on the initial data.
They proved the existence of local solutions in two and three space dimensions (cf. \cite[Theorem 1.1]{zhang2008existenceNewModel}).
In the case of two space dimensions, they additionally showed the existence of global-in-time weak solutions (in the sense that $0 < T < \infty$ is given) (cf. \cite[Remark 1.2]{zhang2008existenceNewModel}).
Moreover, they established the existence of global weak solutions (in the sense that $[0,T) = [0, \infty)$) if the Deborah number and the Reynolds number are sufficiently small for $d \in \{2, 3\}$ (cf. \cite[Theorem 1.2]{zhang2008existenceNewModel}).
However, proving the existence of global-in-time weak solutions to \eqref{intro:eq:dimensional} (for passive or active particles in the sense that $0 < T < \infty$ is given) for three space dimensions without assumptions on the Reynolds and Deborah number was stated as an open problem by Emmrich et al. (cf. \cite[p.42]{emmrich2018LC}).

The main contribution of this paper is the proof of the latter problem. In Main Theorem \ref{thrm:main}, we establish the existence of global-in-time weak solutions to \eqref{intro:eq:dimensional} for both passive and active particles (in the sense that $0 < T < \infty$ is given) in two and three space dimensions under minimal regularity assumptions on the initial data. Moreover, our result holds for arbitrary values of the Reynolds and Deborah number. To be precise, we consider the model $(\mathcal{P})$, which is introduced in Section \ref{sec:model}.

The outline of this paper is as follows:
First, we complete Section \ref{sec:intro} by introducing the model $(\mathcal{P})$, which is an extension of \eqref{intro:eq:dimensional} for active particles.
Then, we present some notation, function spaces, and frequently used formula in Section \ref{sec:preliminares}. To get started, we derive an entropy estimate \eqref{formal:entropy_est} for the model $(\mathcal{P})$ in Section \ref{sec:formal_entropy_estimate} on a formal level.
As we want to make our computations rigorous, we present a discrete-in-time scheme $\lrcb{\mathcal{P}_{\regL}^\tau}$ (see \eqref{def:PLtau}) in Section \ref{sec:discrete_in_time_approx}, where $\tau > 0$ is a small time increment and $\text{L} > 0$ is a cut-off parameter from above. The regularization parameter $\text{L}$ is used to establish the existence of discrete-in-time solutions. Nevertheless, future bounds in energy estimates shall be independent of $\text{L}$.\\
To show the existence of a solution to $\lrcb{\mathcal{P}_{\regL}^\tau}$ in Section \ref{sec:ex_of_sols_to_PLtau}, we introduce a fixed-point scheme in Definition \ref{def:fixed_point_iteration}, prove in Lemma \ref{lmm:fixed_point_iter} that the scheme is well defined, and finally apply Schaefer's fixed-point theorem in Lemma \ref{lmm:existence_fp} to complete the proof.\\
Having established the existence of a solution to $\lrcb{\mathcal{P}_{\regL}^\tau}$, we devote Section \ref{sec:uniform_entropy_estimate} to the derivation of an entropy estimate in Lemma \ref{lmm:time_discrete_entropy} which mimics the formal estimate \eqref{formal:entropy_est}. Unfortunately, the regularity results of the entropy estimate \eqref{est:discrete_entropy} are too weak to pass to the limit in $\lrcb{\mathcal{P}_{\regL}^\tau}$ as $\tau \searrow 0$ and $\text{L} \nearrow \infty$.
Moreover, testing the Smoluchowski equation with $\psi$ is fraud with difficulties to gain more regularity. However, we shall circumvent this problem by deriving an equality for the polymer number density
\begin{align*}
\omega(\mathbf{x}, t) := \int_{\Sph} \psi(\mathbf{x}, \mathbf{m}, t) \, \text{d}\mathscr{H}^{d-1}(\mathbf{m}) \quad \text{ for } (\mathbf{x}, t) \in \Omega \times [0,T]
\end{align*}
and show regularity results for $\omega$ which turn out to be very useful to improve the regularity of $\psi$ afterwards. Such a method is often applied to microscopic macroscopic models (cf. \cite{barrett_sueli_2011_polymers_I, chen_liu_2013_entropy_solution, gruen_metzger_2015microMacro}). The derivation of regularity results for the macroscopic discrete-in-time polymer number density is realized in Section \ref{sec:particle_density}.\\
Then, we establish time regularity results for $\psi$ and $\mathbf{u}$ in Section \ref{sec:time_regularity}. To pass to the limit in Section \ref{sec:passage_limit} as $\tau \searrow 0$ and $\text{L} \nearrow \infty$, we exploit our regularity results to prove strong convergence for suitable subsequences. Instead of using common results of Aubin-Lions or Simon (cf. \cite{lions1969quelques, simon1986compact}), we apply a compactness result of Dreher and J\"ungel (cf. \cite{dreher2012compact}), which relies on hypothesis on time translation and avoids the construction of linear interpolation functions in time.

\subsection{Doi model for active liquid crystals}
\label{sec:model}
Before we extend the model for active particles, let us describe the system \eqref{intro:eq:dimensional}  in more detail.
The velocity field is denoted by $\mathbf{u}: \Omega \times [0,T] \to \mathbb{R}^d$, $p: \Omega \times (0,T] \to \mathbb{R}$ is the pressure, $k_B$ is the Boltzmann constant, $\mathcal{T}$ is the absolute temperature, $\ndens > 0$ is the mean number density of polymers, $D_r = \frac{k_B \mathcal{T}}{\xi_r}$ is the rotary diffusivity, $\xi_r > 0$ is the friction coefficient, $\rho > 0$ is the fluid density, $\mathscr{H}^{d-1}$ is the Hausdorff measure, and $D_\parallel > 0$ and $D_\perp > 0$ are the translational diffusion coefficients parallel and normal to the orientation of the rods.
Moreover, it is always assumed that $\mathbf{m} \in {\SMat^{d-1}}$.
Furthermore, $\nabla_\xflow$ and $\Div_\xflow$ denote the gradient and divergence operator with respect to $\mathbf{x} \in \Omega$. The rotational gradient operator $\sRot$ on the unit sphere is defined by $\sRot = \mathbf{m} \times \nabla_\mflow$, where $\nabla_\mflow$ is the gradient with respect to $\mathbf{m} \in {\SMat^{d-1}}$. Moreover, the chemical potential $\overline{\mu}_\psi$ is given by
\begin{align}
\label{def:muPsi}
\overline{\mu}_\psi := k_B \mathcal{T} \entrop^\prime(\psi) + k_B \mathcal{T} \pot[\psi],
\end{align}
where $\mathcal{F}: \mathbb{R}_{\geq 0} \ni s \mapsto s (\log(s) - 1) + 1 \in \mathbb{R}_{\geq 0}$ is an entropy functional.

The interaction potential $\overline{\pot}[\psi] : {\Omega\times\Sph}\times[0,T] \to \mathbb{R}$ modeling the effects of alignment is defined by
\begin{align}
\label{def:pot}
\pot[\psi](\mathbf{x}, \mathbf{m}, t) :=  \potStr \Je{ \int_{\Sph} \psi(\cdot, \widetilde{\mflow}, t) K(\mathbf{m}, \widetilde{\mflow}) \, \text{d}\mathscr{H}^{d-1}(\widetilde{\mflow})} (\mathbf{x})
\end{align}
for $(\mathbf{x}, \mathbf{m}, t) \in {\Omega\times\Sph} \times [0,T]$, where the dimensionless parameter $\potStr > 0$ denotes the strength of the interaction potential, $\mathcal{J}_{\mollEps}$ is an isotropic mollifier which is defined below, and $K(\mathbf{m}, \widetilde{\mflow})$ is an interaction kernel for $\mathbf{m}, \widetilde{\mflow} \in {\SMat^{d-1}}$ which should have small values for $\mathbf{m}$ and $\widetilde{\mflow}$ being aligned and large values otherwise. In the original Onsager model (cf. \cite{onsager1949colloidal}), the interaction kernel is given by
\begin{align*}
K(\mathbf{m}, \widetilde{\mflow}) = \abs{\mathbf{m} \times \widetilde{\mflow}} \quad \text{ for } \mathbf{m}, \widetilde{\mflow} \in {\SMat^{d-1}}.
\end{align*}
Another common choice is the Maier-Saupe kernel (cf. \cite{maier_saupe_1958})
\begin{align}
\label{def:interaction_kernel}
K(\mathbf{m}, \widetilde{\mflow}) := \abs{\mathbf{m} \times \widetilde{\mflow}}^2 = 1 - (\mathbf{m} \cdot \widetilde{\mflow})^2 \quad \text{ for } \mathbf{m}, \widetilde{\mflow} \in {\SMat^{d-1}},
\end{align}
which has also been used by E and Zhang (cf. \cite{e_zhang_2006_kinetic}).
In this paper, we restrict ourselves to the Maier-Saupe interaction kernel \eqref{def:interaction_kernel}.
Moreover, the isotropic mollifier $\mathcal{J}_{\mollEps}: \L^1(\Omega) \to \mathrm{W}^{1, \infty}(\Omega)$ in \eqref{def:pot} is defined by
\begin{align}
\label{def:mollifier}
\Je{f}(\mathbf{x}) := \int_{\Omega} \zeta_\mollEps(\mathbf{x} - \mathbf{y}) f(\mathbf{y}) \, \text{d}\mathbf{y} \quad \text{ for } f \in \L^1(\Omega),
\end{align}
where $\zeta_\mollEps (\mathbf{x}) := \mollEps^{-d} \zeta(\mollEps^{-1} \mathbf{x})$ for $\mathbf{x} \in \Omega$ and $\mollEps > 0$.
Hereby, the nonnegative, rotationally symmetric, and compactly supported function $\zeta \in \C_0^\infty(\mathbb{R}^d)$ is given by $\zeta(\mathbf{x}) := c \exp(-\frac{1}{1 - \abs{\mathbf{x}}^2})$ if $\abs{\mathbf{x}} < 1$ and $\zeta(\mathbf{x}) = 0$ otherwise. In particular, the constant $c > 0$ is chosen such that $\int_{\mathbb{R}^d} \zeta(\mathbf{x}) \, \text{d}\xflow = 1$.
Therefore, on noting a property of the mollifier (see \eqref{mollifier:fg}), we have that $\overline{\mu}_\psi$ is the first variation of the energy functional
\begin{align*}
\mathcal{E}(\psi)(t) &:= k_B \mathcal{T} \int_{\OmSph}\mathcal{F}(\psi)(\mathbf{x}, \mathbf{m}, t) \,  \text{d}\mathscr{H}^{d-1}(\mathbf{m}) \, \text{d}\xflow \\
&\quad + \frac{ k_B \mathcal{T}}{2} \int_{\OmSph}  \psi(\mathbf{x}, \mathbf{m}, t) \pot[\psi](\mathbf{x}, \mathbf{m}, t) \,  \text{d}\mathscr{H}^{d-1}(\mathbf{m}) \, \text{d}\xflow
\end{align*}
for $t \in [0,T]$.

Now, let us describe the physical interpretation of the terms in \eqref{intro:eq:dimensional} in more detail.
The term on the right-hand side of \eqref{intro:eq:psi:dimensional:srot} describes the rotation of particles induced by the macroscopic velocity gradient $\nabla_\xflow \mathbf{u}$ (note \eqref{intro:eq:psi:dimensional:nablag} for the rewritten form of this term).
As illustrated in Figure \ref{fig:rotated_rod}, a rod $\mathbf{m} \in {\SMat^{d-1}}$ is rotated by the gradient of a macroscopic flow $\nabla_\xflow \mathbf{u}$. However, as the rod cannot elongate and has to maintain its constant length, the rotated rod $\nabla_\xflow \mathbf{u} \mathbf{m}$ is projected onto the tangential space of the unit sphere.
\begin{figure}[!htbp]
    \centering
    \begin{tikzpicture}[scale=1.0]

    \def\ulx{1.374}
    \def\uly{1.624}

    \def\urx{\uly}
    \def\ury{\ulx}

    \def\llx{-\urx}
    \def\lly{-\ury}

    \def\lrx{-\ulx}
    \def\lry{-\uly}

    \def\sulx{2.002}
    \def\suly{0.720}

    \def\surx{2.094}
    \def\sury{0.378}

    \def\sllx{-\surx}
    \def\slly{-\sury}

    \def\slrx{-\sulx}
    \def\slry{-\suly}

    \def\umx{\ulx/2 + \urx/2}
    \def\umy{\uly/2 + \ury/2}

    \def\lmx{\llx/2 + \lrx/2}
    \def\lmy{\lly/2 + \lry/2}

    \def\sumx{\sulx/2 + \surx/2}
    \def\sumy{\suly/2 + \sury/2}

    \def\slmx{\sllx/2 + \slrx/2}
    \def\slmy{\slly/2 + \slry/2}

    \def\mx{\ulx/4 + \urx/4 + \llx/4 + \lrx/4}
    \def\my{\uly/4 + \ury/4 + \lly/4 + \lry/4}

    \draw (\sllx, \slly) -- (\sulx, \suly) -- (\surx, \sury) -- (\slrx, \slry) -- cycle;

    \draw[fill=white] (\llx, \lly) -- (\ulx, \uly) -- (\urx, \ury) -- (\lrx, \lry) -- cycle;

    \draw[arrows={- angle 90}] (\lmx+0.1, \lmy+0.1) -- (\umx-0.1, \umy-0.1) node[pos=0.85, left=16pt, scale=0.75, rotate=45] {$\mathbf{m}$};

    \draw[arrows={- angle 90}, line width=1pt] (\umx, \umy) -- (\sumx, \sumy) node[pos=0.85, below=12pt, scale=0.75, rotate=15] {$\mathbf{(\mathbb{I} - \mathbf{m}\otimes\mathbf{m})} \nabla_{\mathbf{x}}\mathbf{u}\, \mathbf{m}$};

    \draw[dashed, rotate around={15:(\sumx, \sumy)}] (\sumx, \sumy) -- (\sumx + 4, \sumy);

    \draw[arrows={- angle 90}, line width=1pt, rotate around={-10:(\umx, \umy)}] (\umx, \umy) -- (\umx+2.5, \umy) node[pos=1.0, below=7pt, scale=0.75, rotate=15] {$\quad \nabla_{\mathbf{x}} \mathbf{u}\, \mathbf{m}$};

    \node[circle,fill=black,inner sep=0pt,minimum size=3pt] at (\mx, \my) {};
    \end{tikzpicture}
    \caption{A rod $\mathbf{m}$ is rotated by a macroscopic velocity gradient $\nabla_\xflow \mathbf{u}$.}
    \label{fig:rotated_rod}
\end{figure}
The stress tensor on the right-hand side of \eqref{intro:eq:u:dimensional}
\begin{align*}
\overline{\mathbb{T}} := \overline{\mathbb{T}}_{\text{visc, s}} + \overline{\mathbb{T}}_{\text{visc, p}} + \overline{\mathbb{T}}_{\text{elast}}
\end{align*}
is composed of viscous stresses $\overline{\mathbb{T}}_{\text{visc}} := \overline{\mathbb{T}}_{\text{visc, s}} + \overline{\mathbb{T}}_{\text{visc, p}}$ and elastic stresses $\overline{\mathbb{T}}_{\text{elast}}$. Hereby,
\begin{align*}
\overline{\mathbb{T}}_{\text{visc, s}} := \eta_s \lrcb{\nabla_\xflow \mathbf{u} + (\nabla_\xflow \mathbf{u})^T}
\end{align*}
models the viscous stresses coming from the solvent, where $\eta_s > 0$ denotes the solvent viscosity.
The viscous stresses $\overline{\mathbb{T}}_{\text{visc, p}}$ are due to the friction caused by the motion of rodlike polymers. They are related to the hydrodynamic energy dissipation $W$ (cf. \cite[(8.117)]{doi1986theory}) by
\begin{align}
\label{intro:eq:dissip}
W = \int_0^T \int_{\Omega} \nabla_\xflow\mathbf{u} : \overline{\mathbb{T}}_{\text{visc, p}} \, \text{d}\xflow\,\text{dt}.
\end{align}
On noting Figure \ref{fig:rotated_rod}, the velocity of a rod relative to the fluid can be computed by
\begin{align*}
\mathbf{u}_{\text{rod, fluid}} = (\mathbb{I} - \mflow\otimes\mflow) \nabla_\xflow \mathbf{u} \mathbf{m} - \nabla_\xflow \mathbf{u} \mathbf{m} = - (\mflow\otimes\mflow) \nabla_\xflow \mathbf{u} \mathbf{m}.
\end{align*}
Moreover, the frictional force exerted on a rod is given by
\begin{align*}
\mathbf{f}_{\text{rod, frict}} = \frac{ \xi_r}{2} \mathbf{u}_{\text{rod, fluid}},
\end{align*}
where $\xi_r > 0$ is the friction constant. Therefore, the work done by the frictional force is calculated as
\begin{align*}
W &= \int_0^T \int_{\Omega} \int_{\Sph} \psi \mathbf{f}_{\text{rod, frict}} \cdot \mathbf{u}_{\text{rod, fluid}} \, \text{d}\mathscr{H}^{d-1} \, \text{d}\xflow \,\text{dt} \\
&= \frac{ \xi_r}{2} \int_0^T\int_{\Omega}\int_{\Sph} \psi \lrcb{ \nabla_\xflow \mathbf{u} : \mflow\otimes\mflow }^2 \, \text{d}\mathscr{H}^{d-1}(\mathbf{m})\, \text{d}\xflow\,\text{dt}.
\end{align*}
Hence, on noting \eqref{intro:eq:dissip}, we have that
\begin{align}
\label{eq:viscous_stress_p}
\overline{\mathbb{T}}_{\text{visc, p}} = \frac{ \xi_r}{2} \int_{\Sph} \psi (\nabla_\xflow \mathbf{u} : \mflow\otimes\mflow) \mflow\otimes\mflow \, \text{d}\mathscr{H}^{d-1}(\mathbf{m}).
\end{align}
The body force term (the first term on the right-hand side of \eqref{intro:eq:u:dimensional}) and the elastic stress tensor $\overline{\mathbb{T}}_{\text{elast}}$ are related to the change in energy of $\mathcal{E}(\psi)$. In particular, they are inevitable for thermodynamic consistency and they have been derived by a virtual work principle (cf. \cite[Chapter 8.6]{doi1986theory}, \cite[p. 185f]{e_zhang_2006_kinetic}) resulting in an elastic stress tensor given by
\begin{align*}
\overline{\mathbb{T}}_{\text{elast}} := - \int_{\Sph} \psi \lrcb{\mathbf{m} \times \sRot \overline{\mu}_\psi} \otimes \mathbf{m} \, \text{d}\mathscr{H}^{d-1}(\mathbf{m}).
\end{align*}

In our opinion, it is more convenient to work with the surface gradient operator $\nabla_\surfaceLtr$ and the surface divergence operator $\Div_\surfaceLtr$ on ${\SMat^{d-1}}$ instead of $\sRot$.
Noting $\sRot = \mathbf{m} \times \nabla_\mflow$ and $\nabla_\surfaceLtr = (\mathbb{I} - \mflow\otimes\mflow) \nabla_\mflow$, we see that \eqref{intro:eq:psi:dimensional:srot} can equivalently be written as
\begin{align}
\label{intro:eq:psi:dimensional:nablag}
\begin{split}
&\partial_t \psi + \Div_\xflow(\mathbf{u} \psi) - \frac{1}{ k_B \mathcal{T}} \Div_\xflow\lrcb{ \psi\lrcb{ D_\parallel \mflow\otimes\mflow + D_\perp (\mathbb{I} - \mflow\otimes\mflow) }\nabla_\xflow \overline{\mu}_\psi } \\
&\quad - \frac{D_r}{ k_B \mathcal{T}} \Div_\surfaceLtr\lrcb{ \psi \nabla_\surfaceLtr \overline{\mu}_\psi }  = -  \Div_\surfaceLtr\lrcb{ \psi (\mathbb{I} - \mflow\otimes\mflow) \nabla_\xflow \mathbf{u} \mathbf{m} } \quad \text{ in } {\Omega\times\Sph}\times(0,T]
\end{split}
\end{align}
and that the elastic stress tensor can be rewritten as
\begin{equation}
\label{intro:eq:stress_sRot_divg}
\begin{aligned}
\overline{\mathbb{T}}_{\text{elast}} &= - \int_{\Sph} \psi (\mathbf{m} \times \sRot \overline{\mu}_\psi) \otimes \mathbf{m} \, \text{d}\mathscr{H}^{d-1}(\mathbf{m})  \\
&=  \int_{\Sph} \psi \left[ (\mathbb{I} - \mflow\otimes\mflow) \nabla_\surfaceLtr \overline{\mu}_\psi \right] \otimes \mathbf{m} \, \text{d}\mathscr{H}^{d-1}(\mathbf{m}) .
\end{aligned}
\end{equation}

Now, let us come back to active liquid crystals which are systems that are far from equilibrium (cf. \cite{gao2017activeNematics}). The range of active nematic systems includes cytoskeletal filaments (cf. \cite{kierfeld2008active}), bacteria (cf. \cite{darnton2010swarming}), dense suspensions of microswimmers (cf. \cite{wensink2012livingFluids}), and microtubule bundles (cf. \cite{sanchez2012spontaneous}).
Instead of modeling passive particles, we now consider active rodlike polymers which induce active stresses described by an additional stress tensor $\overline{\mathbb{T}}_{\text{act}}$. Following the lines in \cite{gao2017analytical}, we model this stress tensor: We assume that each rod induces a symmetric surface flow (see Figure \ref{fig:active_flow}, Figure \ref{fig:rod})
\begin{align*}
\mathbf{u}_{\text{act}}(s) := \text{sgn}(s) u_\text{act} \mathbf{m}, \quad - \frac{l}{2} \leq s \leq \frac{l}{2},
\end{align*}
where $u_\text{act}$ is the signed surface flow speed, $\mathbf{m}$ is the unit orientation vector of the rod, and $\text{sgn}(\cdot)$ is the sign function.
\begin{figure}[!htbp]
    \centering
    \begin{tikzpicture}[scale=1.0]

    \def\mx{0}
    \def\my{0}
    \def3/4{3/4}
    \def6*\b{6*3/4}
    \def20pt{20pt}
    \def\eps{2pt}

    \draw [|-|] (\mx - 6*\b/2, \my) -- (\mx, \my) node[pos=0, below, scale=0.75] {$\mathbf{x} - \frac{l}{2}$} node[pos=1, below=3pt, scale=0.75] {$\mathbf{x}$};
    \draw [|-|] (\mx, \my) -- (\mx + 6*\b/2, \my) node[pos=1, below, scale=0.75] {$\mathbf{x} + \frac{l}{2}$};

    \draw [->, line width=1pt, red] (\mx - \eps, \my - 20pt) -- (\mx - 6*\b/2, \my - 20pt) node[pos=0.5, below=3pt, scale=0.75] {$- u_\text{act}\, \mathbf{m}$};
    \draw [->, line width=1pt, red] (\mx + \eps, \my - 20pt) -- (\mx + 6*\b/2, \my - 20pt) node[pos=0.5, below=3pt, scale=0.75] {$u_\text{act}\, \mathbf{m}$};

    \draw[arrows={- triangle 90}, line width=1pt, blue] (\mx - 6*\b/2, \my + 20pt*3/4) -- (\mx + 6*\b/2, \my + 20pt*3/4) node[pos=0.5, above=3pt, scale=0.75] {$\mathbf{m}$};
    \end{tikzpicture}
    \caption{Symmetric surface flow induces by an active rod.}
    \label{fig:active_flow}
\end{figure}
Hereby, a rod produces an extensional straining flow along itself for $u_\text{act} > 0$ and a compressive one otherwise. Moreover, let us describe the position of a rod by its centerline $\mathbf{X}(s) := \mathbf{x} + s \mathbf{m}$, $-\frac{l}{2} \leq s \leq \frac{l}{2}$, where $\mathbf{x}$ denotes the center of mass of a rod (see Figure \ref{fig:rod}). Then, the active stress caused by a single rod can be computed as follows on assuming that $\mathbf{x} = 0$ (cf. \cite{gao2017analytical, saintillan2013nonlinear}):
\begin{align*}
\overline{\mathbb{T}}_{\text{act, rod}}
&= - \frac{4 \pi \eta}{\abs{\log(e r^2)}} \int_{-\frac{l}{2}}^{\frac{l}{2}} \mathbf{u}_{\text{act}}(s) \otimes \mathbf{X}(s) \,\text{d}s \\
&= - \frac{\pi \eta u_\text{act} l^2}{\abs{\log(e r^2)}} \mflow\otimes\mflow,
\end{align*}
where $\eta := \eta_s + \eta_p > 0$ is the total viscosity, $\eta_p := \nu \xi_r$ is the polymer viscosity, $e := \exp(1)$, and $r = \frac{b}{l}$ (see Figure \ref{fig:rod}).
Hence, the active stress tensor due to the activity of the rodlike polymers is given by
\begin{align}
\label{intro:def:Tact_dimensional}
\overline{\mathbb{T}}_{\text{act}} = - \frac{\pi \eta u_\text{act} l^2}{\abs{\log(e r^2)}} \int_{\Sph} \psi \lrcb{ \mflow\otimes\mflow - \frac{1}{d}\mathbb{I} } \, \text{d}\mathscr{H}^{d-1}(\mathbf{m}),
\end{align}
where the subtraction \eqref{intro:def:Tact_dimensional} modifies solely the pressure and has no influence on the flow. In particular, the modification in \eqref{intro:def:Tact_dimensional} is equivalent to subtracting the trace of the tensor to make it traceless.\\
It will turn out, note for example the formal energy estimate \eqref{formal:entropy_est}, that the additional active stresses may be seen as an internal source of energy as the thermodynamic equilibrium does not hold any longer in the sense that the total energy at a time $t_2 > t_1$ equals the total energy at time $t_1$ when there is no work done by external forces. This can be explained by the ability of active liquid crystals to convert energy from the local environment to mechanical work (cf. \cite{gao2017activeNematics}).\\
Next, to derive the model in nondimensionalised form, we follow the lines in \cite{e_zhang_2006_kinetic}: Let $L_0 > 0$ be a typical size of the flow region, $V_0 > 0$ the typical velocity scale, and $T_0 := \frac{L_0}{V_0}$ a typical convective time scale.
Moreover, we scale the pseudo probability density $\psi$ by $\frac{1}{\nu}$, where $\nu > 0$ is the mean number density.
We introduce the Deborah number, which is defined by the ratio of the orientational diffusion time scale of the rods and the convective time scale of the fluid, given by
\begin{align*}
\DE := \frac{ \frac{\xi_r}{k_B \mathcal{T}} }{\frac{L_0}{V_0}} = \frac{\xi_r V_0}{k_B \mathcal{T} L_0}.
\end{align*}
Moreover, we note that $\eta_p = \ndens \xi_r$, $\eta = \eta_s + \eta_p > 0$ and set $\viscFrac := \frac{\eta_s}{\eta} \in (0,1)$. Then, the Reynolds number is defined by
\begin{align*}
\RE := \frac{\rho V_0 L_0}{\eta}.
\end{align*}
Furthermore, we set $\mollEps := \frac{l}{L_0} > 0$, where $l>0$ denotes the length of rods (see Figure \ref{fig:rod}).
We define the dimensionless translational diffusion parameters by
\begin{align*}
\hat{D}_\parallel := \frac{\xi_r}{l^2 k_B \mathcal{T}} D_\parallel \ \text{ and } \ \hat{D}_\perp := \frac{\xi_r}{l^2 k_B \mathcal{T}} D_\perp.
\end{align*}
It usually holds that $\hat{D}_\parallel > \hat{D}_\perp > 0$ (\cite[p. 296f]{doi1986theory}) and in most of the cases one has that $\hat{D}_\parallel = 2 \hat{D}_\perp$ (\cite[p. 296f]{doi1986theory}). From an analytical point of view, it is important that the tensor
\begin{align*}
\hat{D}_\parallel \mflow\otimes\mflow + \hat{D}_\perp (\mathbb{I} - \mflow\otimes\mflow)
= \hat{D}_\perp \mathbb{I} + (\hat{D}_\parallel - \hat{D}_\perp) \mflow\otimes\mflow
\end{align*}
is symmetric and positive definite, which is true if $\hat{D}_\parallel \geq \hat{D}_\perp > 0$. Hence, we set $\hat{D}_\parallel = \hat{D}_\perp = 1$ for ease of notation in this paper.
On noting \eqref{intro:def:Tact_dimensional}, we introduce the dimensionless activity parameter
\begin{align*}
\activity := - \frac{\pi \eta l^2 u_\text{act}}{\abs{\log(e r^2)} k_B \mathcal{T}}.
\end{align*}

Thus, the active model in nondimensionalised form (without relabelling $\psi$, $\mathbf{u}$, $p$, $\mathbf{x}$, and $t$ for ease of notation) reads as follows:\\
\hspace*{2ex}Given $T > 0$, find $\mathbf{u}: \Omega \times [0,T] \ni (\mathbf{x}, t) \mapsto \mathbf{u}(\mathbf{x}, t) \in \mathbb{R}^d$,\\
\hspace*{2ex}$p: \Omega \times (0,T] \ni (\mathbf{x}, t) \mapsto p(\mathbf{x}, t) \in \mathbb{R}$,\\
\hspace*{2ex}and $\psi: {\Omega\times\Sph}\times[0,T] \ni (\mathbf{x}, \mathbf{m}, t) \mapsto \psi(\mathbf{x}, \mathbf{m}, t) \in \mathbb{R}$ such that
\begin{subequations}
    \label{eq:active:stokes}
    \begin{alignat}{3}
    \partial_t \mathbf{u} + (\mathbf{u} \cdot \nabla_\xflow)\mathbf{u} + \nabla_\xflow p &= - \frac{1 - \viscFrac}{\RE\DE} \int_{\Sph} \psi \nabla_\xflow \muPsi \, \text{d}\mathscr{H}^{d-1} \label{eq:active:u} \\
    &\quad + \Div_\xflow(\mathbb{T})  \quad && \text{in } \Omega \times (0,T], \nonumber \\
    \Div_\xflow(\mathbf{u}) &= 0 \quad &&\text{in } \Omega \times (0,T], \label{eq:active:divu} \\
    \mathbf{u} &= 0 \quad &&\text{on } \partial\Omega \times (0,T], \label{eq:active:bc:u} \\
    \mathbf{u}(\mathbf{x}, 0) &= \mathbf{u}_0(\mathbf{x}) \quad &&\forall \mathbf{x} \in \Omega \label{eq:active:init:u}
    \end{alignat}
\end{subequations}
and
\begin{subequations}
    \label{eq:active:smoluchowski}
    \begin{align}
    & \partial_t \psi + \mathbf{u} \cdot \nabla_\xflow \psi - \frac{\mollEps^2}{\DE} \Div_\xflow (\psi \nabla_\xflow \muPsi) - \frac{1}{\DE} \Div_\surfaceLtr \lrcb{ \psi \nabla_\surfaceLtr \muPsi } \nonumber \\
    &\quad = - \Div_\surfaceLtr \lrcb{ (\mathbb{I} - \mflow\otimes\mflow) \nabla_\xflow \mathbf{u} \mathbf{m} \psi } \qquad \text{in } \Omega \times {\SMat^{d-1}} \times (0,T], \label{eq:active:psi} \\
    &\muPsi := \entrop^\prime(\psi) + \pot[\psi] \qquad\qquad\qquad\qquad\ \ \, \text{in } \Omega \times {\SMat^{d-1}} \times (0,T], \label{eq:active:muPsi}\\
    &\psi \nabla_\xflow \muPsi \cdot \nx = 0 \qquad\qquad\qquad\qquad\qquad\ \, \text{on } \partial \Omega \times {\SMat^{d-1}} \times (0,T], \label{eq:active:bc:psi} \\
    &\psi(\mathbf{x}, \mathbf{m}, 0) = \psi_0(\mathbf{x}, \mathbf{m}) \qquad\qquad\qquad\quad\ \ \ \forall (\mathbf{x}, \mathbf{m}) \in {\Omega\times\Sph}, \label{eq:active:init:psi}
    \end{align}
\end{subequations}
where $\nx : \partial\Omega \to \mathbb{R}^d$ denotes the outward normal to $\partial\Omega$ and the stress tensor in the Navier-Stokes type equation is given by
\begin{align*}
\mathbb{T} := \mathbb{T}_{\text{visc, s}} + \mathbb{T}_{\text{visc, p}} + \mathbb{T}_{\text{elast}} + \mathbb{T}_{\text{act}}
\end{align*}
with
\begin{align*}
\mathbb{T}_{\text{visc, s}} &:= \frac{\viscFrac}{\RE} \lrcb{\nabla_\xflow \mathbf{u} + (\nabla_\xflow \mathbf{u})^T}, \\
\mathbb{T}_{\text{visc, p}} &:= \frac{1 - \viscFrac}{2 \RE}  \int_{\Sph} \psi (\nabla_\xflow \mathbf{u} : \mflow\otimes\mflow) \mflow\otimes\mflow \, \text{d}\mathscr{H}^{d-1}(\mathbf{m}), \\
\mathbb{T}_{\text{elast}} &:= \frac{1 - \viscFrac}{\RE\DE}  \int_{\Sph} \psi \left[ (\mathbb{I} - \mflow\otimes\mflow) \nabla_\surfaceLtr \muPsi \right] \otimes \mathbf{m} \, \text{d}\mathscr{H}^{d-1}(\mathbf{m}), \\
\mathbb{T}_{\text{act}} &:= \frac{\activity(1 - \viscFrac)}{\RE\DE}  \int_{\Sph} \psi \lrcb{\mflow\otimes\mflow - \frac{1}{d}\mathbb{I}} \, \text{d}\mathscr{H}^{d-1}(\mathbf{m}).
\end{align*}

The initial data $\mathbf{u}_0$ and $\psi_0$ will be specified later.

\begin{defi}
    The coupled problem \eqref{eq:active:stokes} - \eqref{eq:active:smoluchowski} will be referred to as model $(\mathcal{P})$ throughout this paper. We call this model passive if $\alpha = 0$ and active if $\alpha \neq 0$.
\end{defi}

\begin{rmrk}
    As $\activity$ has the opposite sign of $u_\text{act}$, contractile stress is modeled by $\activity > 0$ and extensile stress by $\activity < 0$.
\end{rmrk}

\section{Preliminaries}
\label{sec:preliminares}

\subsection{Notation}
We denote the real numbers by $\mathbb{R}$, the natural numbers by $\mathbb{N} := \{1, 2, \dots \}$ and set $\mathbb{N}_0 := \mathbb{N} \cup \{0\}$. In particular, we define $\mathbb{R}_{\geq 0} := [0, \infty)$ and $\mathbb{R}_{>0} := (0, \infty)$.
Moreover, we introduce the following scalar products: For $\mathbf{v}, \mathbf{w} \in \mathbb{R}^d$ and $\mathbb{A}, \mathbb{B} \in \mathbb{R}^{d \times d}$, we define $\mathbf{v} \cdot \mathbf{w} := \sum_{i=1}^d \mathbf{v}_i \mathbf{w}_i$ and $\mathbb{A} : \mathbb{B} := \sum_{i,j=1}^d \mathbb{A}_{ij} \mathbb{B}_{ij}$.
The trace of a tensor $\mathbb{A} \in \mathbb{R}^{d \times d}$ is given by $\tr{\mathbb{A}} := \sum_{i=1}^d \mathbb{A}_{ii}$.
The notation $\abs{\cdot}$ is used as follows: For a real number $a \in \mathbb{R}$, we set $\abs{a}$ to the absolute value of $a$, whereas $\abs{\mathbf{v}} := (\mathbf{v} \cdot \mathbf{v})^{1/2}$ and $\abs{\mathbb{A}} := (\mathbb{A} : \mathbb{A})^{1/2}$ for $\mathbf{v} \in \mathbb{R}^d$ and $\mathbb{A} \in \mathbb{R}^{d \times d}$.
Moreover, we set $\abs{U}$ to the $d$-dimensional Lebesgue measure of a measurable set $U \subset \mathbb{R}^d$.
For $\mathbf{v}, \mathbf{w} \in \mathbb{R}^d$, we define $\mathbf{v} \otimes \mathbf{w} := ( \mathbf{v}_i \mathbf{w}_j )_{1 \leq i,j \leq d}$.
We introduce the positive part $\pos{a}$ and the negative part $\negat{a}$ of a real number $a \in \mathbb{R}$ by $[a]_{\pm} := \frac{1}{2}(a \pm \abs{a})$.
For ease of notation, we write $\dH$ instead of $\text{d}\mathscr{H}^{d-1}$ in integrals for the Hausdorff measure.
The isotropic mollifier $\mathcal{J}_{\mollEps}: \L^1(\Omega) \to \mathrm{W}^{1, \infty}(\Omega)$ (cf. \eqref{def:mollifier}) is naturally extended to the tensor valued mollifier $\mathbb{J}_{\mollEps}: \mathbb{L}^1(\Omega) \to \mathbb{W}^{1, \infty}(\Omega)$.
The velocity gradient is defined by $(\nabla_\xflow \mathbf{u})_{1 \leq i,j \leq d} = (\partial_{\mathbf{x}_j} \mathbf{u}_i)$ for a differentiable function $\mathbf{u} : \Omega \to \mathbb{R}^d$. Moreover, we introduce the divergence of a differentiable function $\mathbb{A} : \Omega \to \mathbb{R}^{d \times d}$ by
\begin{align*}
\Div_\xflow(\mathbb{A}) := \lrcb{ \sum_{j=1}^d \partial_{\mathbf{x}_j} \mathbb{A}_{ij}}_{i=1}^d.
\end{align*}

\subsection{Function spaces}
We use the standard notation for Lebesgue and Sobolev spaces and define $\mathbf{L}^p(\Omega) := \L^p(\Omega; \mathbb{R}^d)$, $\mathbf{W}^{k, p}(\Omega) := \mathrm{W}^{k,p}(\Omega; \mathbb{R}^d)$, $\mathbf{H}^k(\Omega) := \mathbf{W}^{k, 2}($ $\Omega)$, $\mathbb{L}^p(\Omega) := \L^p(\Omega; \mathbb{R}^{d \times d})$, $\mathbb{W}^{k,p}(\Omega)$ $:=$ $\mathrm{W}^{k,p}(\Omega;$ $ \mathbb{R}^{d \times d})$, and $\mathbb{H}^k(\Omega) := \mathbb{W}^{k, 2}(\Omega)$ for $1 \leq p \leq \infty$ and $k \in \mathbb{N}$. Moreover, we introduce the spaces
\begin{align*}
\Hflow_{0,\text{div}}^1(\Omega) &:= \left\{ \mathbf{w} \in \mathbf{H}_0^1(\Omega):\ \Div_\xflow(\mathbf{w}) = 0 \right\}, \\
\Ldivp{2} &:= \left\{ \mathbf{w} \in \mathbf{L}^2(\Omega):\ \int_{\Omega} \mathbf{w} \cdot \nabla_\xflow \varphi \, \text{d}\xflow = 0 \quad \forall \varphi \in \H^1(\Omega) \right\}, \\
\mathscr{V} &:= \left\{ \mathbf{w} \in \mathbf{C}_0^\infty(\Omega):\ \Div_\xflow(\mathbf{w}) = 0 \right\}, \\
\Vk{k} &:= \overline{\mathscr{V}}^{\norm{\cdot}_{\mathbf{H}_0^1(\Omega) \cap \mathbf{H}^k(\Omega)}},\  k \in \mathbb{N}.
\end{align*}
We note that it holds that $\Vk{1} = \Hflow_{0,\text{div}}^1(\Omega)$ (cf. Temam \cite[Chapter 1, Theorem 1.6]{temam1977navierNumerics}).
An introduction to Sobolev spaces $\mathrm{W}^{k,p}({\SMat^{d-1}})$ on the sphere can be found in \cite{kufner1977functionSpaces}. Furthermore, we define the dual spaces of $\mathrm{W}^{k,p}(\Omega \times {\SMat^{d-1}})$ and $\Vk{k}$ by $(\mathrm{W}^{k,p}(\Omega \times {\SMat^{d-1}}))^\prime$ and $(\Vk{k})^\prime$, respectively. It will turn out useful to introduce the space
\begin{align*}
\Zp{p} := \left\{ \theta \in \LpOS{p}:\ \theta \geq 0 \text{ a.e. in } {\Omega\times\Sph} \right\} \quad \forall 1 \leq p < \infty.
\end{align*}

\subsection{Frequently used formula}
Let us note some results which will be often referred to in future discussions.

\begin{lmm}
	Let the isotropic mollifier $\mathcal{J}_{\mollEps}: \L^1(\Omega) \to \mathrm{W}^{1, \infty}(\Omega)$ be defined by \eqref{def:mollifier}. Then, it holds that
	\begin{subequations}
		\begin{alignat}{2}
		\norm{ \Je{f} }_{\L^2(\Omega)} &\leq \norm{f}_{\L^2(\Omega)} \quad &&\forall f \in \L^2(\Omega), \label{mollifier:L2_stab}\\

		\int_{\Omega} \Je{f} g \, \text{d}\xflow &= \int_{\Omega} f \Je{g} \, \text{d}\xflow \quad &&\forall f,g \in \L^2(\Omega), \label{mollifier:fg}\\

		\norm{\Je{f}}_{\mathrm{W}^{1,\infty}(\Omega)} &\leq \overline{\C}_\mollEps \norm{f}_{\L^1(\Omega)} \quad &&\forall f \in \L^1(\Omega). \label{mollifier:W_2_infty}
		\end{alignat}
	\end{subequations}
	Moreover, similar results hold for the tensor valued mollifier $\mathbb{J}_{\mollEps}$.
\end{lmm}

Furthermore, we will exploit the following formula on the unit sphere (cf. \cite[Lemma 2.1, Lemma 1.1]{chen_liu_2013_entropy_solution}):

\begin{lmm}
	Let $\mathbf{v} \in \mathbf{H}^1({\SMat^{d-1}})$ and $f, g \in \H^1({\SMat^{d-1}})$. Then, it holds that
	\begin{subequations}
		\begin{align}
		\int_{\Sph} \Div_\surfaceLtr(\mathbf{v}) f \dH &= - \int_{\Sph} \mathbf{v} \cdot \nabla_\surfaceLtr f \dH + (d - 1) \int_{\Sph} (\mathbf{v} \cdot \mathbf{m}) f \dH(\mflow), \label{eq:pInt:divg} \\
		\int_{\Sph} (\Delta_\surfaceLtr f) g \dH &= - \int_{\Sph} \nabla_\surfaceLtr f \cdot \nabla_\surfaceLtr g \dH. \label{eq:pInt:nablag}
		\end{align}
	\end{subequations}
\end{lmm}

\begin{lmm}
	Let $f \in \mathrm{W}^{1,1}({\SMat^{d-1}})$ and $\mathbb{A} \in \R^{\text{d}\times \text{d}}$ be a constant matrix with $\tr{\mathbb{A}} = 0$. Then, it holds that
	\begin{align}
	\label{eq:chen_stress_formula}
	\int_{\Sph} \lrcb{ (\mathbb{I} - \mflow\otimes\mflow) \mathbb{A} \mathbf{m} } \cdot \nabla_\surfaceLtr f \dH(\mflow)
	= \int_{\Sph} f (d \mflow\otimes\mflow - \mathbb{I}) : \mathbb{A} \dH(\mflow).
	\end{align}
\end{lmm}

Moreover, the following lemma is frequently used, which follows directly from $\nabla_\surfaceLtr = (\mathbb{I} - \mflow\otimes\mflow) \nabla_\mflow$, $(\mathbf{m} \cdot \widetilde{\mflow})^2 = (\mflow\otimes\mflow) : (\widetilde{\mflow} \otimes \widetilde{\mflow})$ for $\mathbf{m}, \widetilde{\mflow} \in {\SMat^{d-1}}$, and \eqref{def:pot}.
\begin{lmm}
	\label{lmm:nablag_g_pot_formula}
	It holds that
	\begin{align}
	\label{eq:nabla_g_pot_formula}
	\nabla_\surfaceLtr \pot[\theta](\mathbf{x}, \mathbf{m}) = - 2 \potStr (\mathbb{I} - \mflow\otimes\mflow) \JeMat{ \int_{\Sph} \theta(\cdot, \widetilde{\mflow}) \widetilde{\mflow} \otimes \widetilde{\mflow} \dH(\widetilde{\mflow}) } (\mathbf{x}) \mathbf{m}
	\end{align}
	for all $\theta \in \LpOS{1}$ for $(\mathbf{x}, \mathbf{m}) \in {\Omega\times\Sph}$.
\end{lmm}

\subsection{Formal entropy estimate}
\label{sec:formal_entropy_estimate}
Now, let us derive an entropy estimate on a formal level which exploits the energetic and entropic structure of the model $(\mathcal{P})$.
First, assuming that $\psi > 0$, multiplying \eqref{eq:active:psi} by $\muPsi$, integrating over ${\Omega\times\Sph}$, integrating by parts in $\Omega$ and on ${\SMat^{d-1}}$ (see \eqref{eq:pInt:divg} and \eqref{eq:pInt:nablag}), and noting \eqref{eq:active:divu}-\eqref{eq:active:bc:u}, \eqref{eq:active:bc:psi}, and \eqref{mollifier:fg}, we have that
\begin{align}
\label{formal:psi:muPsi}
\begin{split}
&\frac{\text{d}}{\,\text{dt}} \int_{\OmSph} \left[ \mathcal{F}(\psi) + \frac{1}{2} \psi \pot[\psi] \right] \dH\, \text{d}\xflow \\
&\qquad + \frac{\mollEps^2}{\DE} \int_{\OmSph} \psi \abs{\nabla_\xflow \muPsi}^2 \dH\, \text{d}\xflow + \frac{1}{\DE} \int_{\OmSph} \psi \abs{\nabla_\surfaceLtr \muPsi}^2 \dH\, \text{d}\xflow \\
&\quad = - \int_{\OmSph} \mathbf{u} \cdot \nabla_\xflow \psi \pot[\psi] \dH\, \text{d}\xflow \\
&\qquad + \int_{\OmSph} \psi \lrcb{ (\mathbb{I} - \mflow\otimes\mflow) \nabla_\xflow\mathbf{u} \mathbf{m} } \cdot \nabla_\surfaceLtr \muPsi \dH(\mflow)\, \text{d}\xflow.
\end{split}
\end{align}

Next, multiplying \eqref{eq:active:u} by $\mathbf{u}$, integrating over $\Omega$, noting \eqref{eq:active:divu} and \eqref{eq:active:bc:u} and some integration by parts yield
\begin{align*}
&\frac{1}{2} \frac{\text{d}}{\,\text{dt}} \int_{\Omega} \abs{\mathbf{u}}^2 \, \text{d}\xflow + \frac{\viscFrac}{\RE} \int_{\Omega} \abs{\nabla_\xflow \mathbf{u}}^2 \, \text{d}\xflow \\
&\qquad + \frac{1 - \viscFrac}{2 \RE} \int_{\OmSph} \psi \lrcb{ \nabla_\xflow \mathbf{u} : \mflow\otimes\mflow }^2  \dH(\mflow)\, \text{d}\xflow \\
&\quad= - \frac{1 - \viscFrac}{\RE\DE} \int_{\OmSph} \psi \lrcb{ (\mathbb{I} - \mflow\otimes\mflow) \nabla_\xflow \mathbf{u} \mathbf{m} } \cdot \nabla_\surfaceLtr \muPsi \dH(\mflow)\, \text{d}\xflow \\
&\quad\quad + \frac{1 - \viscFrac}{\RE\DE} \int_{\OmSph} \mathbf{u} \cdot \nabla_\xflow \psi \pot[\psi] \dH\, \text{d}\xflow \numberthis[formal:u:u] \\
&\quad\quad - \frac{\alpha (1 - \viscFrac)}{\RE \DE} \int_{\OmSph} \psi \nabla_\xflow \mathbf{u} : \mflow\otimes\mflow \dH(\mflow)\, \text{d}\xflow,
\end{align*}
where we used that
\begin{align}
\label{formal:stress_formula}
\begin{split}
&\int_{\OmSph} \psi \nabla_\xflow\mathbf{u} : \lrcb{ \left[ (\mathbb{I} - \mflow\otimes\mflow) \nabla_\surfaceLtr \muPsi \right] \otimes \mathbf{m} } \dH(\mflow)\, \text{d}\xflow \\
&\quad = \int_{\OmSph} \psi \lrcb{ (\mathbb{I} - \mflow\otimes\mflow) \nabla_\xflow \mathbf{u} \mathbf{m}} \cdot \nabla_\surfaceLtr \muPsi \dH(\mflow)\, \text{d}\xflow.
\end{split}
\end{align}

Now, we infer from multiplying \eqref{formal:psi:muPsi} by $\frac{1-\viscFrac}{\RE\DE}$, combining \eqref{formal:psi:muPsi} - \eqref{formal:u:u}, and applying Young's inequality that
\begin{align*}
&\frac{1}{2} \frac{\text{d}}{\,\text{dt}} \int_{\Omega} \abs{\mathbf{u}}^2 \, \text{d}\xflow + \frac{\viscFrac}{\RE} \int_{\Omega} \abs{\nabla_\xflow \mathbf{u}}^2 \, \text{d}\xflow \\
&\qquad + \frac{1 - \viscFrac}{2 \RE} \int_{\OmSph} \psi \lrcb{\nabla_\xflow\mathbf{u} : \mflow\otimes\mflow}^2 \dH(\mflow)\, \text{d}\xflow \\
&\qquad + \frac{1 - \viscFrac}{\RE\DE} \frac{\text{d}}{\,\text{dt}} \int_{\OmSph} \left[ \mathcal{F}(\psi) + \frac{1}{2}\psi\pot[\psi] \right] \dH\, \text{d}\xflow \\
&\qquad + \frac{1 - \viscFrac}{\RE\DE} \lrcb{ \frac{\mollEps^2}{\DE} \int_{\OmSph} \psi \abs{\nabla_\xflow \muPsi}^2 \dH\, \text{d}\xflow + \frac{1}{\DE} \int_{\OmSph} \psi \abs{\nabla_\surfaceLtr \muPsi}^2 \dH\, \text{d}\xflow } \\
&\quad= - \frac{\alpha (1 - \viscFrac)}{\RE \DE} \int_{\OmSph} \psi \nabla_\xflow \mathbf{u} : \mflow\otimes\mflow \dH(\mflow)\, \text{d}\xflow \numberthis[formal:summarize_1] \\
&\quad\leq \frac{1 - \viscFrac}{4\RE} \int_{\OmSph} \psi \lrcb{\nabla_\xflow\mathbf{u} : \mflow\otimes\mflow}^2 \dH(\mflow)\, \text{d}\xflow + \frac{\alpha^2 (1 - \viscFrac)}{\RE \DE^2} \int_{\OmSph} \psi \dH\, \text{d}\xflow.
\end{align*}

Furthermore, to deal with the last term in \eqref{formal:summarize_1}, we obtain from \eqref{eq:active:divu}-\eqref{eq:active:bc:u}, \eqref{eq:active:bc:psi}, \eqref{eq:pInt:divg}, and integrating \eqref{eq:active:psi} over ${\Omega\times\Sph}$ that
\begin{align}
\label{formal:mass_conserv}
\int_{\OmSph} \psi(\mathbf{x}, \mathbf{m}, t) \dH(\mflow)\, \text{d}\xflow = \int_{\OmSph} \psi_0(\mathbf{x}, \mathbf{m}) \dH(\mflow)\, \text{d}\xflow
\end{align}
for $t \in (0, T]$.
Therefore, we conclude from \eqref{formal:summarize_1} and \eqref{formal:mass_conserv} that
\begin{align*}
&\frac{\text{d}}{\,\text{dt}} \left[ \frac{1}{2} \int_{\Omega} \abs{\mathbf{u}}^2 \, \text{d}\xflow + \frac{1 - \viscFrac}{\RE\DE} \int_{\OmSph} \left[ \mathcal{F}(\psi) + \frac{1}{2}\psi\pot[\psi] \right] \dH\, \text{d}\xflow \right] \\
&\quad + \frac{\viscFrac}{\RE} \int_{\Omega} \abs{\nabla_\xflow \mathbf{u}}^2 \, \text{d}\xflow + \frac{1 - \viscFrac}{4\RE} \int_{\OmSph} \psi \lrcb{\nabla_\xflow\mathbf{u} : \mflow\otimes\mflow}^2 \dH(\mflow)\, \text{d}\xflow \\
&\quad + \frac{1 - \viscFrac}{\RE\DE} \lrcb{ \frac{\mollEps^2}{\DE} \int_{\OmSph} \psi \abs{\nabla_\xflow \muPsi}^2 \dH\, \text{d}\xflow + \frac{1}{\DE} \int_{\OmSph} \psi \abs{\nabla_\surfaceLtr \muPsi}^2 \dH\, \text{d}\xflow } \\
&\leq \frac{\alpha^2 (1 - \viscFrac)}{\RE \DE^2} \int_{\OmSph} \psi_0 \dH\, \text{d}\xflow. \numberthis[formal:entropy_est]
\end{align*}

We note that there is no need for an absorption argument in \eqref{formal:summarize_1} in the passive case as the right-hand sides vanishes for $\activity = 0$.

\section{Discrete-in-time approximations}
\label{sec:discrete_in_time_approx}
To make the formal computations in Section \ref{sec:formal_entropy_estimate} rigorous, we introduce the following regularization functions: For any $\text{L} > 1$, we define
\begin{subequations}
	\label{def:regularizations}
	\begin{align}
	\begin{split}
	\label{def:regL}
	\entrop^\regL(s) &:=
	\begin{cases}
	\mathcal{F}(s) \equiv s (\log(s) - 1) + 1 \quad & \text{ for } 0 \leq s \leq \text{L}, \\
	\frac{s^2 - \text{L}^2}{2 \text{L}} + s(\log(\text{L}) - 1) + 1 \quad & \text{ for } \text{L}\leq s,
	\end{cases}
	\end{split} \\
	\begin{split}
	\label{def:QL}
	\text{Q}^\regL(s) &:=
	\begin{cases}
	s \quad & \text{ for } s \leq \text{L}, \\
	\text{L} \quad & \text{ for } \text{L} \leq s,
	\end{cases}
	\end{split} \\
	\begin{split}
	\label{def:QZL}
	\text{Q}_0^\regL(s) &:= \max\{ 0,\ \text{Q}^\regL(s) \} \quad \text{ for } s \in \mathbb{R}.
	\end{split}
	\end{align}
\end{subequations}

Furthermore, we note that the regularization functions satisfy the following properties (cf. Chen and Liu \cite{chen_liu_2013_entropy_solution} and Barrett and S\"uli \cite{barrett_sueli_2011_polymers_I, barrett_sueli_2011_polymers_II}):
\begin{subequations}
	\label{reg:props}
	\begin{alignat}{2}
	\text{Q}^\regL \in \C^{0,1}(\mathbb{R}),\quad \entrop^\regL &\in \C^{2,1}(\mathbb{R}_{>0}) \cap \C(\mathbb{R}_{\geq 0}), \\
	\entrop^\regL & \text{ is convex on } [0, \infty), \quad && \label{reg:F_cvx} \\
	\entrop^\regL(s) & \geq \mathcal{F}(s) \quad && \forall s \geq 0, \label{reg:FL_geq_F} \\
	\entrop^\regL\lrcb{\text{Q}^\regL(s) + \delta} & \leq \delta + \frac{\delta^2}{2} + \mathcal{F}(s + \delta) \quad && \forall \delta \in (0,1), \forall s \geq 0, \label{reg:FLQL_leq_FL} \\
	\lrcb{\entropL}^{\prime\prime}(s) &= \lrcb{\text{Q}^\regL(s)}^{-1} \geq \frac{1}{s} \quad &&\forall s > 0, \label{reg:ddFL_geq} \\
	\entrop^{\prime\prime}(s + \delta) & \leq \frac{1}{\delta} \quad && \forall \delta \in (0,1),\ \forall s \geq 0, \label{reg:ddFL_leq} \\
	\entrop^\regL(s) &\geq \frac{s^2}{4\text{L}} - C(\text{L}) \quad && \forall s \geq 0, \label{reg:FL_geq_ss} \\
	\abs{\text{Q}_0^\regL(s)} &\leq \abs{\text{Q}^\regL(s)} \leq \abs{s} \quad && \forall s \in \mathbb{R}. \label{reg:QL_leq}
	\end{alignat}
\end{subequations}

\textit{\textbf{Assumptions on the initial data and the domain.}}
We assume the domain and the initial data to satisfy
\begin{subequations}
	\label{assumptions:domain_init_data}
	\begin{align}
	\begin{split}
	\label{eq:domain}
	& \Omega \subset \mathbb{R}^d \text{ is an open bounded Lipschitz domain},\quad d \in \{ 2, 3\},
	\end{split}\\
	\begin{split}
	\label{init:psi_0}
	& \psi_0 \in \L^2(\Omega; \L^1({\SMat^{d-1}})), \quad \psi_0 \geq 0 \text{ a.e. in } {\Omega\times\Sph}  \\
	& \qquad \text{ with } \mathcal{F}(\psi_0) \in \LpOS{1},
	\end{split}\\
	\begin{split}
	\label{init:u_0}
	\mathbf{u}_0 \in \Ldivp{2}.
	\end{split}
	\end{align}
\end{subequations}

Let a time $T > 0$ and a number of time steps $N \in \mathbb{N}$ be given. Then, we define a time increment $\tau > 0$ by $\tau := \frac{T}{N}$ and $t_n := n \tau$ for $n \in \{0, \dots, N\}$. Moreover, we define the initial velocity $\uL{0} \in \Hflow_{0,\text{div}}^1(\Omega)$ for our discrete-in-time scheme as the weak solution of
\begin{align}
\label{init:def:u_L_0}
\int_{\Omega} \left[ \uL{0} \cdot \mathbf{w} + \text{L}^{-1} \nabla_\xflow \uL{0} : \nabla_\xflow \mathbf{w} \right] \, \text{d}\xflow = \int_{\Omega} \mathbf{u}_0 \cdot \mathbf{w} \, \text{d}\xflow \quad \forall \mathbf{w} \in \Hflow_{0,\text{div}}^1(\Omega).
\end{align}
We deduce from \eqref{init:def:u_L_0} that
\begin{align}
\label{init:regularity_u_L_0}
\int_{\Omega} \left[ \abs{\uL{0}}^2 + \text{L}^{-1} \abs{\nabla_\xflow \uL{0}}^2 \right] \, \text{d}\xflow \leq \int_{\Omega} \abs{\mathbf{u}_0}^2 \, \text{d}\xflow \leq C < \infty,
\end{align}
where $C > 0$ is independent of $\tau$ and $\text{L}$. Furthermore, we set
\begin{align}
\label{init:def:psi_L_0}
\psiL{0} := \text{Q}^\regL(\psi_0) \in \Zp{2}.
\end{align}
In particular, we infer from \eqref{init:psi_0}, \eqref{init:def:psi_L_0}, and \eqref{reg:QL_leq} that
\begin{align}
\label{init:regularity_psi_L_0}
\int_{\OmSph} \psiL{0} \dH\, \text{d}\xflow \leq \int_{\OmSph} \psi_0 \dH\, \text{d}\xflow \leq C < \infty,
\end{align}
where $C>0$ is independent of $\tau$ and $\text{L}$.
Moreover, we obtain from the definition of the interaction potential \eqref{def:pot}, the nonnegativity of the Maier-Saupe interaction kernel \eqref{def:interaction_kernel}, \eqref{init:psi_0}, \eqref{init:def:psi_L_0}, \eqref{reg:QL_leq}, \eqref{mollifier:L2_stab}, and H\"older's inequality that
\begin{align*}
&\int_{\OmSph} \psiL{0} \pot[\psiL{0}] \dH\, \text{d}\xflow \\
&\quad\leq C(\potStr, d) \int_{\Omega} \lrcb{ \int_{\Sph} \psi_0(\mathbf{x}, \mathbf{m}) \dH(\mflow) }^2 \, \text{d}\xflow
\leq C < \infty,	\numberthis[init:psi_0_pot_0]
\end{align*}
where $C > 0$ is independent of $\tau$ and $\text{L}$.

Our discrete-in-time problem $\lrcb{\mathcal{P}_{\regL}^\tau}$ reads as follows:\\
\hspace*{3ex} $\lrcb{\mathcal{P}_{\regL}^\tau}$ Let $\uL{0} \in \Hflow_{0,\text{div}}^1(\Omega)$ and $\psiL{0} \in \Zp{2}$ be given by \eqref{init:def:u_L_0} and \eqref{init:def:psi_L_0}, respectively. Then, for $n \in \{1, \dots, N\}$, given $(\uL{n-1}, \psiL{n-1})$ $\in$ $\Hflow_{0,\text{div}}^1(\Omega) \times \Zp{2}$, find $(\uL{n}, \psiL{n})$ $\in$ $\Hflow_{0,\text{div}}^1(\Omega) \times (\HkOS{1} \cap \Zp{2})$ such that, for all $\mathbf{w} \in \Hflow_{0,\text{div}}^1(\Omega)$ and $\theta \in \H^1(\Omega \times {\SMat^{d-1}})$,
\begin{subequations}
	\label{def:PLtau}
	\begin{align}
	\begin{split}
	\label{eq:PLtau:u}
	& \RE \DE \int_{\Omega} \frac{\uL{n} - \uL{n-1}}{\tau} \cdot \mathbf{w} \, \text{d}\xflow + \RE \DE \int_{\Omega} \lrcb{ (\uL{n-1} \cdot \nabla_\xflow) \uL{n} } \cdot \mathbf{w} \, \text{d}\xflow \\
	&\qquad + \viscFrac \DE \int_{\Omega} \nabla_\xflow \uL{n} : \nabla_\xflow \mathbf{w} \, \text{d}\xflow \\
	&\qquad + \frac{(1 - \viscFrac) \DE}{2} \int_{\OmSph} \text{Q}_0^\regL(\psiL{n}) (\nabla_\xflow \uL{n} : \mflow\otimes\mflow) (\nabla_\xflow \mathbf{w} : \mflow\otimes\mflow) \dH(\mflow) \, \text{d}\xflow \\
	&\quad = - (1 - \viscFrac) \int_{\OmSph} \psiL{n} \nabla_\xflow \potL{n} \cdot \mathbf{w} \dH\, \text{d}\xflow  \\
	&\quad \quad - (1 - \viscFrac) \int_{\OmSph} \text{Q}_0^\regL(\psiL{n}) \lrcb{ (\mathbb{I} - \mflow\otimes\mflow) \nabla_\xflow \mathbf{w} \mathbf{m} } \cdot \nabla_\surfaceLtr \potL{n} \dH(\mflow)\, \text{d}\xflow \\
	&\quad \quad - (1 - \viscFrac) \int_{\OmSph} \psiL{n} (d\mflow\otimes\mflow - \mathbb{I}) : \nabla_\xflow \mathbf{w} \dH(\mflow)\, \text{d}\xflow \\
	&\quad \quad - \activity(1 - \viscFrac) \int_{\OmSph} \text{Q}_0^\regL(\psiL{n}) \mflow\otimes\mflow : \nabla_\xflow \mathbf{w} \dH(\mflow)\, \text{d}\xflow,
	\end{split} \\
	\begin{split}
	\label{eq:PLtau:psi}
	&\int_{\OmSph} \frac{\psiL{n} - \psiL{n-1}}{\tau} \theta \dH\, \text{d}\xflow - \int_{\OmSph} \psiL{n} \uL{n} \cdot \nabla_\xflow \theta \dH\, \text{d}\xflow \\
	&\quad + \frac{1}{\DE} \int_{\OmSph} \left[ \mollEps^2 \nabla_\xflow \psiL{n} \cdot \nabla_\xflow \theta + \nabla_\surfaceLtr \psiL{n} \cdot \nabla_\surfaceLtr \theta \right] \dH\, \text{d}\xflow \\
	&\quad + \frac{1}{\DE} \int_{\OmSph} \left[ \mollEps^2 \text{Q}_0^\regL(\psiL{n}) \nabla_\xflow \potL{n} \cdot \nabla_\xflow \theta + \text{Q}_0^\regL(\psiL{n}) \nabla_\surfaceLtr \potL{n} \cdot \nabla_\surfaceLtr \theta \right] \dH\, \text{d}\xflow \\
	&\quad - \int_{\OmSph} \text{Q}_0^\regL(\psiL{n}) \lrcb{(\mathbb{I} - \mflow\otimes\mflow) \nabla_\xflow \uL{n} \mathbf{m}} \cdot \nabla_\surfaceLtr \theta \dH(\mflow)\, \text{d}\xflow = 0.
	\end{split}
	\end{align}
\end{subequations}

\begin{rmrk}
	\begin{enumerate}
		\item[i)] On noting \eqref{formal:stress_formula} and \eqref{eq:chen_stress_formula}, the elastic stress tensor in the Navier-Stokes equation can be rewritten as
		\begin{align*}
		& \int_{\OmSph} \Div_\xflow \lrcb{ \psi \left[ (\mathbb{I} - \mflow\otimes\mflow) \nabla_\surfaceLtr \muPsi \right] \otimes \mathbf{m} } \cdot \mathbf{w} \dH(\mflow)\, \text{d}\xflow \\
		&\quad = - \int_{\OmSph} \psi \lrcb{ (\mathbb{I} - \mflow\otimes\mflow) \nabla_\xflow \mathbf{w} \mathbf{m} } \cdot \nabla_\surfaceLtr \muPsi  \dH(\mflow)\, \text{d}\xflow \\
		&\quad = - \int_{\OmSph} \psi (d \mflow\otimes\mflow - \mathbb{I}) : \nabla_\xflow \mathbf{w} \dH(\mflow)\, \text{d}\xflow \numberthis[rmrk:stress_tensor] \\
		&\quad \qquad - \int_{\OmSph} \psi \lrcb{ (\mathbb{I} - \mflow\otimes\mflow) \nabla_\xflow \mathbf{w} \mathbf{m} } \cdot \nabla_\surfaceLtr \pot[\psi] \dH(\mflow)\, \text{d}\xflow
		\end{align*}
		for all $\mathbf{w} \in \Hflow_{0,\text{div}}^1(\Omega)$.
		\item[ii)] For reasons of readability, we use the abbreviation
		\begin{align}
		\label{def:potL}
		\potL{n} := \frac{1}{2} \pot\left[ \psiL{n} + \psiL{n-1} \right] \qquad \text{ for } n \in \{1, \dots, N\}.
		\end{align}
	\end{enumerate}
\end{rmrk}

\subsection{Existence of solutions to \texorpdfstring{$\lrcb{\mathcal{P}_{\regL}^\tau}$}{}}
\label{sec:ex_of_sols_to_PLtau}
In this section, our goal is to prove the existence of a solution $(\uL{n}, \psiL{n})$ $\in$ $\Hflow_{0,\text{div}}^1(\Omega) \times (\HkOS{1} \cap \Zp{2})$ to $\lrcb{\mathcal{P}_{\regL}^\tau}$ by applying Schaefer's fixed-point theorem (cf. Evans \cite{evans2010pdes}).
First of all, we prove a lemma which will be used very frequently.

\begin{lmm}
	Let $\theta \in \L^1({\Omega\times\Sph})$. Then, it holds that
	\begin{align}
	\label{est:nabla_pot_Linfty}
	\norm{\pot[\theta]}_{\mathrm{W}^{1, \infty}({\Omega\times\Sph})}
	\leq C(\overline{\C}_\mollEps, \potStr, d) \norm{\theta}_{\L^1({\Omega\times\Sph})},
	\end{align}
	where $\overline{\C}_\mollEps > 0$ is the constant in \eqref{mollifier:W_2_infty}.
\end{lmm}

\begin{proof}
	On noting the definition of the interaction potential \eqref{def:pot} and $(\mathbf{m} \cdot \widetilde{\mflow})^2 = \mflow\otimes\mflow : \widetilde{\mflow} \otimes \widetilde{\mflow}$ for $\mathbf{m}, \widetilde{\mflow} \in {\SMat^{d-1}}$, we have that
	\begin{align*}
		\pot[\theta](\mathbf{x}, \mathbf{m})
		&= \potStr \int_{\Sph} \Je{\theta(\cdot, \widetilde{\mflow})}(\mathbf{x}) \lrcb{1 - (\mathbf{m} \cdot \widetilde{\mflow})^2} \dH(\widetilde{\mflow}) \\
		&= \potStr \Je{ \int_{\Sph} \theta(\cdot, \widetilde{\mflow}) \dH(\widetilde{\mflow}) }(\mathbf{x}) \\
		&\quad - \potStr \JeMat{ \int_{\Sph} \theta(\cdot, \widetilde{\mflow}) \widetilde{\mflow} \otimes \widetilde{\mflow} \dH(\widetilde{\mflow}) }(\mathbf{x}) : (\mflow\otimes\mflow)
	\end{align*}
	for $(\mathbf{x}, \mathbf{m}) \in {\Omega\times\Sph}$. Hence, we infer from \eqref{mollifier:W_2_infty} that
	\begin{align*}
		&\norm{\pot[\theta]}_{\LpOS{\infty}} + \norm{\nabla_\xflow \pot[\theta]}_{\L^\infty({\Omega\times\Sph})} \\
		&\quad\leq C(\overline{\C}_\mollEps, \potStr) \norm{\theta}_{\LpOS{1}}  + C(\overline{\C}_\mollEps, \potStr, d) \int_{\OmSph} \abs{\theta(\mathbf{x}, \widetilde{\mflow})} \abs{\widetilde{\mflow} \otimes \widetilde{\mflow}} \dH(\widetilde{\mflow})\, \text{d}\xflow \\
		&\quad\leq C(\overline{\C}_\mollEps, \potStr, d) \norm{\theta}_{\LpOS{1}}.
	\end{align*}
	Furthermore, we infer from \eqref{mollifier:W_2_infty} and \eqref{eq:nabla_g_pot_formula} that
	\begin{align*}
	\norm{ \nabla_\surfaceLtr \pot[\theta] }_{\L^\infty({\Omega\times\Sph})}
	&\quad\leq C(d) \potStr \esssup_{\mathbf{x} \in \Omega} \abs{ \JeMat{ \int_{\Sph} \theta(\cdot, \widetilde{\mflow}) \widetilde{\mflow} \otimes \widetilde{\mflow} \dH(\widetilde{\mflow}) } (\mathbf{x}) } \\
	&\quad\leq C(\overline{\C}_\mollEps, \potStr, d) \norm{\theta}_{\L^1({\Omega\times\Sph})}.
	\end{align*}
\end{proof}

\textit{\textbf{Fixed-point iteration.}} To introduce our fixed-point scheme, we rewrite \eqref{eq:PLtau:u} as
\begin{align}
\label{fp:eq:a}
a(\psiL{n})(\uL{n}, \mathbf{w}) = k(\psiL{n})(\mathbf{w}) \quad \forall \mathbf{w} \in \Hflow_{0,\text{div}}^1(\Omega),
\end{align}
where
\begin{align}
\label{fp:def:a}
\begin{split}
&a(\overline{\psi})(\mathbf{u}, \mathbf{w})\\
&\quad:= \RE \DE \int_{\Omega} \mathbf{u} \cdot \mathbf{w} \, \text{d}\xflow + \tau \RE \DE \int_{\Omega} \lrcb{ (\uL{n-1} \cdot \nabla_\xflow) \mathbf{u} } \cdot \mathbf{w} \, \text{d}\xflow \\
&\quad\quad + \tau \viscFrac \DE \int_{\Omega} \nabla_\xflow \mathbf{u} : \nabla_\xflow \mathbf{w} \, \text{d}\xflow \\
&\quad\quad + \tau \frac{(1-\viscFrac) \DE}{2} \int_{\OmSph} \text{Q}_0^\regL(\overline{\psi}) (\nabla_\xflow \mathbf{u} : \mflow\otimes\mflow) (\nabla_\xflow \mathbf{w} : \mflow\otimes\mflow) \dH(\mflow) \, \text{d}\xflow
\end{split}
\end{align}
and
\begin{align}
\label{fp:def:k}
\begin{split}
&k(\overline{\psi})(\mathbf{w})\\
&\quad:= \RE \DE \int_{\Omega} \uL{n-1} \cdot \mathbf{w} \, \text{d}\xflow - \tau \frac{1 - \viscFrac}{2} \int_{\OmSph} \overline{\psi} \nabla_\xflow \pot\left[\overline{\psi} + \psiL{n-1}\right] \cdot \mathbf{w} \dH\, \text{d}\xflow \\
&\quad\quad - \tau \frac{1 - \viscFrac}{2} \int_{\OmSph} \text{Q}_0^\regL(\overline{\psi}) \lrcb{ (\mathbb{I} - \mflow\otimes\mflow) \nabla_\xflow \mathbf{w} \mathbf{m} } \cdot \nabla_\surfaceLtr \pot\left[\overline{\psi} + \psiL{n-1}\right] \dH(\mflow)\, \text{d}\xflow \\
&\quad\quad - \tau (1 - \viscFrac) \int_{\OmSph} \overline{\psi} (d\mflow\otimes\mflow - \mathbb{I}) : \nabla_\xflow \mathbf{w} \dH(\mflow)\, \text{d}\xflow \\
&\quad\quad - \tau \activity (1 - \viscFrac) \int_{\OmSph} \text{Q}_0^\regL(\overline{\psi}) \mflow\otimes\mflow : \nabla_\xflow \mathbf{w} \dH(\mflow)\, \text{d}\xflow
\end{split}
\end{align}
for all $\mathbf{u}, \mathbf{w} \in \Hflow_{0,\text{div}}^1(\Omega)$ and for all $\overline{\psi} \in \L^2({\Omega\times\Sph})$.
Moreover, we rewrite \eqref{eq:PLtau:psi} as
\begin{align}
\label{fp:eq:b}
\begin{split}
b(\uL{n})(\psiL{n}, \theta) = l(\uL{n}, \psiL{n})(\theta) \qquad \forall \theta \in \HkOS{1},
\end{split}
\end{align}
where
\begin{align}
\label{fp:def:b}
\begin{split}
b(\mathbf{u})(\psi, \theta)
&:= \int_{\OmSph} \psi \theta \dH\, \text{d}\xflow - \tau \int_{\OmSph} \psi \mathbf{u} \cdot \nabla_\xflow \theta \dH\, \text{d}\xflow \\
&\quad + \frac{\tau}{\DE} \int_{\OmSph} \left[ \mollEps^2 \nabla_\xflow \psi \cdot \nabla_\xflow \theta + \nabla_\surfaceLtr \psi \cdot \nabla_\surfaceLtr \theta \right] \dH\, \text{d}\xflow
\end{split}
\end{align}
and
\begin{align}
\label{fp:def:l}
\begin{split}
l(\mathbf{u}, \overline{\psi})(\theta)
&:= \int_{\OmSph} \psiL{n-1} \theta \dH\, \text{d}\xflow \\
&\quad - \frac{\tau}{2 \DE} \int_{\OmSph} \mollEps^2 \text{Q}_0^\regL(\overline{\psi}) \nabla_\xflow \pot\left[\overline{\psi} + \psiL{n-1}\right] \cdot \nabla_\xflow \theta \dH\, \text{d}\xflow \\
&\quad - \frac{\tau}{2 \DE} \int_{\OmSph} \text{Q}_0^\regL(\overline{\psi}) \nabla_\surfaceLtr \pot\left[\overline{\psi} + \psiL{n-1}\right] \cdot \nabla_\surfaceLtr \theta \dH\, \text{d}\xflow \\
&\quad + \tau \int_{\OmSph} \text{Q}_0^\regL(\overline{\psi}) \lrcb{ (\mathbb{I} - \mflow\otimes\mflow) \nabla_\xflow \mathbf{u} \mathbf{m} } \cdot \nabla_\surfaceLtr \theta \dH(\mflow)\, \text{d}\xflow
\end{split}
\end{align}
for all $\psi, \theta \in \HkOS{1}$ and for all $\overline{\psi} \in \LpOS{2}$, $\mathbf{u} \in \Hflow_{0,\text{div}}^1(\Omega)$.
\begin{lmm}
	\label{lmm:fixed_point_iter}
	Suppose that the assumptions \eqref{assumptions:domain_init_data} hold and let $\overline{\psi} \in \L^2({\Omega\times\Sph})$.
	\begin{enumerate}
		\item[i)] There exists a unique $\mathbf{u} \in \Hflow_{0,\text{div}}^1(\Omega)$ such that
		\begin{align}
		\label{fp:iteration:i}
		a(\overline{\psi})(\mathbf{u}, \mathbf{w}) = k(\overline{\psi})(\mathbf{w}) \quad \forall \mathbf{w} \in \Hflow_{0,\text{div}}^1(\Omega).
		\end{align}
		\item[ii)] Now, let $\mathbf{u} \in \Hflow_{0,\text{div}}^1(\Omega)$ be the unique solution to \eqref{fp:iteration:i}. Then, there exists a unique $\psi$ $\in$ $\HkOS{1}$ such that
		\begin{align}
		\label{fp:iteration:ii}
		b(\mathbf{u})(\psi, \theta) = l(\mathbf{u}, \overline{\psi})(\theta) \quad \forall \theta \in \HkOS{1}.
		\end{align}
		\item[iii)] Moreover, solutions to \eqref{fp:iteration:i} and \eqref{fp:iteration:ii} satisfy
		\begin{align}
		\label{fp:iteration:iii:u}
		\begin{split}
			&\DE\, \min\{ \RE, \tau\viscFrac \} \norm{\mathbf{u}}_{\H^1(\Omega)} \\
			&\quad \leq C \lrcb{ 1 + \tau \text{L} + \tau \text{L} \norm{\overline{\psi}}_{\LpOS{2}} + \tau \norm{\overline{\psi}}_{\LpOS{2}}^2 },
		\end{split}
		\end{align}
		where $C = C(\RE, \DE, \uL{n-1}, \overline{\C}_\mollEps, \potStr, d, \viscFrac, \activity, \psiL{n-1}) > 0$ and
		\begin{align}
		\label{fp:iteration:iii:psi}
		\begin{split}
			&\min\left\{ 1, \frac{\tau}{\DE}, \frac{\tau \mollEps^2}{\DE} \right\} \norm{\psi}_{\HkOS{1}} \\
			&\quad \leq C \lrcb{ 1 + \tau \text{L} + \tau \text{L} \norm{\overline{\psi}}_{\LpOS{1}} + \tau \text{L} \norm{\nabla_\xflow \mathbf{u}}_{\L^2(\Omega)} },
		\end{split}
		\end{align}
		where $C = C(\psiL{n-1}, \overline{\C}_\mollEps, \potStr, d, \mollEps, \DE^{-1}) > 0$.
	\end{enumerate}
\end{lmm}

\begin{defi}[Fixed-point iteration]
	\label{def:fixed_point_iteration}
	We define the operator $\Xi: \LpOS{2} \to \LpOS{2}$ as follows: For given $\overline{\psi} \in \LpOS{2}$, we introduce $\Xi[\overline{\psi}] \in \H^1(\Omega \times {\SMat^{d-1}})$ as the solution of the variational problem \eqref{fp:iteration:ii}, where $\mathbf{u} \in \Hflow_{0,\text{div}}^1(\Omega)$ is the solution of the variational problem \eqref{fp:iteration:i}.
\end{defi}

\begin{rmrk}
	We infer from Lemma \ref{lmm:fixed_point_iter} that the operator $\Xi$ defined in Definition \ref{def:fixed_point_iteration} is well defined.
	Clearly, if $\psi \in \HkOS{1}$ is a fixed-point of $\Xi$, then the tuple $(\uL{n}, \psiL{n})$ $:=$ $(\mathbf{u}, \psi)$ is a solution to \eqref{def:PLtau}. Hence, proving the existence of a fixed-point of $\Xi$ is equivalent to proving the existence of a solution to \eqref{def:PLtau}.
\end{rmrk}

\begin{proof}[Proof of Lemma \ref{lmm:fixed_point_iter}]
	Let $\overline{\psi} \in \L^2({\Omega\times\Sph})$.\\
	\textit{Proof of i).}
	For all $\mathbf{u} \in \Hflow_{0,\text{div}}^1(\Omega)$, on noting
	\begin{align*}
	\int_{\Omega} \lrcb{ (\uL{n-1} \cdot \nabla_\xflow) \mathbf{u} } \cdot \mathbf{u} \, \text{d}\xflow = 0,
	\end{align*}
	we have that
	\begin{align}
	\label{fp:a_coercive}
	\begin{split}
	a(\overline{\psi})(\mathbf{u}, \mathbf{u})
    \geq \DE \, \min\{ \RE, \tau \viscFrac \} \norm{\mathbf{u}}_{\H^1(\Omega)}^2.
	\end{split}
	\end{align}

	Moreover, we infer from $\H^1(\Omega) \hookrightarrow \mathbf{L}^6(\Omega)$ that, for all $\mathbf{u}, \mathbf{w} \in \Hflow_{0,\text{div}}^1(\Omega)$,
	\begin{align}
	\label{fp:a_bounded}
	\begin{split}
	\abs{a(\overline{\psi})(\mathbf{u}, \mathbf{w})}
	\leq C\lrcb{ \tau, \RE, \DE, \gamma, d, \uL{n-1}, \text{L} } \norm{\mathbf{u}}_{\H^1(\Omega)} \norm{\mathbf{w}}_{\H^1(\Omega)}.
	\end{split}
	\end{align}

	Furthermore, we obtain from \eqref{reg:QL_leq} and \eqref{est:nabla_pot_Linfty} that
	\begin{align}
	\label{fp:k_bounded}
	\begin{split}
	\abs{k(\overline{\psi})(\mathbf{w})}
    \leq C \lrcb{1 + \tau \text{L} + \tau \norm{\overline{\psi}}_{\LpOS{2}}^2 + \tau\text{L} \norm{\overline{\psi}}_{\LpOS{2}}} \norm{\mathbf{w}}_{\H^1(\Omega)},
	\end{split}
	\end{align}
	where $C = C \lrcb{ \RE, \DE, \uL{n-1}, \overline{\C}_\mollEps, \potStr, d, \viscFrac, \activity, \psiL{n-1} } > 0$ and $\overline{\C}_\mollEps > 0$ is the constant in \eqref{mollifier:W_2_infty}.
	We conclude from \eqref{fp:a_coercive}-\eqref{fp:a_bounded} that $a(\overline{\psi})(\cdot, \cdot)$ is a bounded, coercive bilinear form on $\Hflow_{0,\text{div}}^1(\Omega) \times \Hflow_{0,\text{div}}^1(\Omega)$ and from \eqref{fp:k_bounded} that $k(\overline{\psi})(\cdot) \in (\Hflow_{0,\text{div}}^1(\Omega))^\prime$. Hence, we obtain from Lax-Milgram's theorem the existence of a unique solution $\mathbf{u} \in \Hflow_{0,\text{div}}^1(\Omega)$ to \eqref{fp:iteration:i}.

	\textit{Proof of ii).} Let $\mathbf{u} \in \Hflow_{0,\text{div}}^1(\Omega)$ be the unique solution to \eqref{fp:iteration:i}. For all $\psi \in \H^1(\Omega \times {\SMat^{d-1}})$,
	on noting
	\begin{align*}
	\int_{\OmSph} \psi \mathbf{u} \cdot \nabla_\xflow \psi \dH\, \text{d}\xflow = \int_{\OmSph} \frac{1}{2} \mathbf{u} \cdot \nabla_\xflow (\psi^2) \dH\, \text{d}\xflow = 0,
	\end{align*}
	we have that
	\begin{align}
	\label{fp:b_coercive}
	\begin{split}
	b(\mathbf{u})(\psi, \psi)
	\geq \min \left\{ 1, \frac{\tau}{\DE}, \frac{\tau}{\DE}\mollEps^2 \right\} \norm{\psi}_{\HkOS{1}}^2.
	\end{split}
	\end{align}

	Furthermore, we infer from $\HkOS{1} \hookrightarrow \LpOS{3}$ that
	\begin{align}
	\label{fp:b_bounded}
	\begin{split}
	\abs{b(\mathbf{u})(\psi, \theta)}
	\leq \lrcb{ 1 + \frac{\tau}{\DE} + \frac{\tau}{\DE}\mollEps^2 + C(d) \tau \norm{\mathbf{u}}_{\L^6(\Omega)} } \norm{\psi}_{\HkOS{1}} \norm{\theta}_{\HkOS{1}}
	\end{split}
	\end{align}
	for all $\psi, \theta \in$ $\H^1(\Omega$ $\times{\SMat^{d-1}})$.
	Moreover, we obtain from \eqref{est:nabla_pot_Linfty} that
	\begin{align}
	\label{fp:l_bounded}
	\begin{split}
	\abs{l(\mathbf{u}, \overline{\psi})(\theta)}
	\leq C \Big( 1 + \tau \text{L} + \tau \text{L} \norm{\overline{\psi}}_{\LpOS{1}} + \tau \text{L} \norm{\nabla_\xflow \mathbf{u}}_{\L^2(\Omega)} \Big) \norm{\theta}_{\HkOS{1}},
	\end{split}
	\end{align}
	where $C = C(\psiL{n-1}, \overline{\C}_\mollEps, \potStr, d, \mollEps, \DE^{-1}) > 0$.
	Consequently, on noting \eqref{fp:b_coercive}-\eqref{fp:l_bounded}, it follows that $b(\mathbf{u})(\cdot, \cdot)$ is a bounded, coercive bilinear form on $\H^1(\Omega \times {\SMat^{d-1}})$ $\times$ $\HkOS{1}$ and that $l(\mathbf{u}, \overline{\psi})(\cdot)$ $\in$ $\lrcb{\HkOS{1}}^\prime$. Hence, applying Lax-Milgram's theorem completes the proof of \textit{ii)}.

	\textit{Proof of iii).}
	On noting \eqref{fp:a_coercive}, \eqref{fp:k_bounded}, \eqref{fp:b_coercive}, and \eqref{fp:l_bounded}, it follows immediately that \eqref{fp:iteration:iii:u} and \eqref{fp:iteration:iii:psi} hold.
\end{proof}

As proving the existence of a fixed-point of $\Xi$ is equivalent to proving the existence of a solution to \eqref{def:PLtau}, let us show the following lemma:
\begin{lmm}[Existence of a fixed-point]
	\label{lmm:existence_fp}
	Suppose that the assumptions \eqref{assumptions:domain_init_data} hold.
	Then, the operator $\Xi: \LpOS{2} \to \LpOS{2}$ has a fixed-point.
\end{lmm}

\begin{proof}
	The idea to prove this lemma is applying Schaefer's fixed-point theorem (cf. Evans \cite{evans2010pdes}).
	Before doing so however we need to verify the following three conditions:
	\begin{enumerate}
		\item[i)] $\Xi: \LpOS{2} \to \LpOS{2}$ is continuous.
		\item[ii)] $\Xi: \LpOS{2} \to \LpOS{2}$ is compact.
		\item[iii)] The set
		\begin{align*}
		\left\{ \psi \in \LpOS{2}:\ \psi = \lambda \Xi(\psi) \text{ for some } 0 < \lambda \leq 1 \right\}
		\end{align*}
		is bounded in $\LpOS{2}$.
	\end{enumerate}

	\textit{Proof of i).} Let $\psi, \psi_k \in \LpOS{2}$, $k \in \mathbb{N}$, such that
	\begin{subequations}
		\begin{align}
		\begin{split}
		\label{fp:ex:i:cvg:psiov_k}
		\psi_k \to \psi \quad \text{ in } \LpOS{2} \text{ as } k \to \infty.
		\end{split}
		\end{align}
	\end{subequations}

	Our goal is to show that
	\begin{align}
	\label{fp:ex:i:goal}
	\Xi[\psi_k] \to \Xi[\psi] \quad \text{ in } \LpOS{2} \text{ as } k \to \infty.
	\end{align}

	First, we infer from Lemma \ref{lmm:fixed_point_iter} and the definition of $\Xi$ that there exist unique $\mathbf{u}, \mathbf{u}_k$ $\in$ $\Hflow_{0,\text{div}}^1(\Omega)$, $k \in \mathbb{N}$, such that
	\begin{subequations}
		\begin{alignat}{2}
		a(\psi)(\mathbf{u}, \mathbf{w}) &= k(\psi)(\mathbf{w}) \quad && \forall \mathbf{w} \in \Hflow_{0,\text{div}}^1(\Omega), \label{fp:ex:i:a_u}\\
		a(\psi_k)(\mathbf{u}_k, \mathbf{w}) &= k(\psi_k)(\mathbf{w}) \quad && \forall \mathbf{w} \in \Hflow_{0,\text{div}}^1(\Omega),\ \forall k \in \mathbb{N}, \label{fp:ex:i:a_uk}\\
		b(\mathbf{u})(\Xi[\psi], \theta) &= l(\mathbf{u}, \psi)(\theta) \quad && \forall \theta \in \HkOS{1}, \label{fp:ex:i:b_psi}\\
		b(\mathbf{u}_k)(\Xi[\psi_k], \theta) &= l(\mathbf{u}_k, \psi_k)(\theta) \quad && \forall \theta \in \HkOS{1},\ \forall k \in \mathbb{N}. \label{fp:ex:i:b_psi_k}
		\end{alignat}
	\end{subequations}

	Now, let us derive uniform bounds on $\{\Xi[\psi_k]\}_{k \in \mathbb{N}}$ in $\HkOS{1}$ and on $\{ \mathbf{u}_k \}_{k \in \mathbb{N}}$ in $\Hflow_{0,\text{div}}^1(\Omega)$ to extract converging subsequences in appropriate spaces.
	Taking $\theta := \Xi[\psi_k]$ $\in$ $\HkOS{1}$ as a test function in \eqref{fp:ex:i:b_psi_k} and applying \eqref{fp:iteration:iii:psi}, we have that
	\begin{align}
	\label{fp:ex:i:psi_k_bound_2}
	\norm{\Xi[\psi_k]}_{\HkOS{1}}
	\leq C \lrcb{ 1 + \norm{\psi_k}_{\LpOS{1}} + \norm{\nabla_\xflow \mathbf{u}_k}_{\L^2(\Omega)} },
	\end{align}
	where $C = C\lrcb{ \psiL{n-1}, \overline{\C}_\mollEps, \potStr, d, \text{L}, \tau, \tau^{-1}, \DE, \DE^{-1}, \mollEps, \mollEps^{-1} } > 0$.
	So far, we do not control the right-hand side of \eqref{fp:ex:i:psi_k_bound_2}.
	Consequently, it remains to derive uniform bounds on $\{ \mathbf{u}_k \}_{k \in \mathbb{N}}$ in $\Hflow_{0,\text{div}}^1(\Omega)$.
	Taking $\mathbf{w} := \mathbf{u}_k$ $\in$ $\Hflow_{0,\text{div}}^1(\Omega)$ as a test function in \eqref{fp:ex:i:a_uk}, we deduce from \eqref{fp:iteration:iii:u} that
	\begin{align}
	\label{fp:ex:i:uk_bound_1}
	\begin{split}
	&\DE \, \min \left\{ \RE, \tau \viscFrac \right\} \norm{\mathbf{u}_k}_{\H^1(\Omega)}^2 \\
	&\quad \leq C \lrcb{1 + \norm{\psi_k}_{\LpOS{2}} + \norm{\psi_k}_{\LpOS{2}}^2} \norm{\mathbf{u}_k}_{\H^1(\Omega)},
	\end{split}
	\end{align}
	where $C = C(\RE, \DE, \uL{n-1}, \overline{\C}_\mollEps, \potStr, d, \tau, \viscFrac, \text{L}, \activity, \psiL{n-1}) > 0$.
	Moreover, we obtain from \eqref{fp:ex:i:cvg:psiov_k} that there exists a constant $C > 0$ independent of $k$ such that
	\begin{align}
	\label{fp:ex:i:psiovk_bound}
	\norm{\psi_k}_{\LpOS{2}} \leq C \quad \forall k \in \mathbb{N}.
	\end{align}
	Therefore, we deduce from \eqref{fp:ex:i:uk_bound_1} and \eqref{fp:ex:i:psiovk_bound} that
	\begin{align}
	\label{fp:ex:i:uk_bound_2}
	\norm{\mathbf{u}_k}_{\H^1(\Omega)}
	\leq C,
	\end{align}
	where $C = C\lrcb{ \RE, \RE^{-1}, \DE, \DE^{-1}, \viscFrac, \viscFrac^{-1}, \tau, \tau^{-1}, \uL{n-1}, \overline{\C}_\mollEps, \potStr, d, \text{L}, \activity, \psiL{n-1} } > 0$ is independent of $k$.
	Hence, inserting \eqref{fp:ex:i:psiovk_bound} and \eqref{fp:ex:i:uk_bound_2} in \eqref{fp:ex:i:psi_k_bound_2} yields
	\begin{align}
	\label{fp:ex:i:psi_k_bound_3}
	\norm{\Xi[\psi_k]}_{\HkOS{1}}
	\leq  C,
	\end{align}
	where $C = C\lrcb{ \uL{n-1}, \psiL{n-1}, \RE, \RE^{-1}, \DE, \DE^{-1}, \tau, \tau^{-1}, \viscFrac, \viscFrac^{-1}, \mollEps, \mollEps^{-1}, \overline{\C}_\mollEps, \text{L}, \activity, \potStr, d } > 0$ is independent of $k$.
	Therefore, on noting \eqref{fp:ex:i:uk_bound_2}, \eqref{fp:ex:i:psi_k_bound_3}, the Lipschitz continuity of $\text{Q}_0^\regL$, and the compact embeddings
	\begin{align*}
	\HkOS{1} \hookrightarrow \hookrightarrow \LpOS{2} \quad \text{ and } \quad \Hflow_{0,\text{div}}^1(\Omega) \hookrightarrow \hookrightarrow \mathbf{L}^2(\Omega),
	\end{align*}
	we deduce that there exists a subsequence $\{ \mathbf{u}_{k_i}, \Xi[\psi_{k_i}] \}_{i \in \mathbb{N}}$ of $\{ \mathbf{u}_k, \Xi[\psi_k] \}_{k \in \mathbb{N}}$ and functions $\hat{\uflow} \in \Hflow_{0,\text{div}}^1(\Omega)$ and $\hat{\psi} \in \HkOS{1}$ such that
	\begin{subequations}
		\begin{alignat}{2}
		\Xi[\psi_{k_i}] &\to \hat{\psi} \quad &&\text{ strongly in } \LpOS{2}, \label{fp:ex:i:cvg:psi_ki_strong} \\
		\Xi[\psi_{k_i}] &\rightharpoonup \hat{\psi} \quad &&\text{ weakly in } \HkOS{1}, \label{fp:ex:i:cvg:psi_ki_weak}\\

		\text{Q}_0^\regL(\psi_{k_i}) &\to \text{Q}_0^\regL(\psi) \quad &&\text{ strongly in } \LpOS{2}, \label{fp:ex:i:cvg:QZL_psiov_ki}\\
		\mathbf{u}_{k_i} & \to \hat{\uflow} \quad &&\text{ strongly in } \mathbf{L}^2(\Omega), \label{fp:ex:i:cvg:u_ki_strong}\\
		\mathbf{u}_{k_i} & \rightharpoonup \hat{\uflow} \quad &&\text{ weakly in } \Hflow_{0,\text{div}}^1(\Omega) \label{fp:ex:i:cvg:u_ki_weak}
		\end{alignat}
	\end{subequations}
	as $i \to \infty$.
	In particular, we infer from \eqref{est:nabla_pot_Linfty} and $\psi_{k_i} \to \psi$ strongly in $\LpOS{2}$ that
	\begin{align}
	\label{fp:ex:i:cvg:pot_psiov_ki}
	\pot[\psi_{k_i}] \to \pot[\psi] \quad \text{ strongly in } \mathrm{W}^{1, \infty}({\Omega\times\Sph}) \text{ as } i \to \infty.
	\end{align}
	Moreover, we obtain from \eqref{fp:ex:i:a_uk} that
	\begin{align}
	\label{fp:ex:i:a_uki}
	a(\psi_{k_i})(\mathbf{u}_{k_i}, \mathbf{w}) = k(\psi_{k_i})(\mathbf{w}) \quad \forall \mathbf{w} \in \Hflow_{0,\text{div}}^1(\Omega),\ \forall i \in \mathbb{N}.
	\end{align}
	On noting \eqref{fp:ex:i:cvg:psi_ki_weak}, \eqref{fp:ex:i:cvg:QZL_psiov_ki}, \eqref{fp:ex:i:cvg:u_ki_strong}, \eqref{fp:ex:i:cvg:u_ki_weak}, and \eqref{fp:ex:i:cvg:pot_psiov_ki}, we pass to the limit in \eqref{fp:ex:i:a_uki} as $i \to \infty$ to deduce that
	\begin{align}
	\label{fp:ex:i:a_uhat}
	a(\psi)(\hat{\uflow}, \mathbf{w}) = k(\psi)(\mathbf{w}) \quad \forall \mathbf{w} \in \Hflow_{0,\text{div}}^1(\Omega).
	\end{align}

	Furthermore, it follows from \eqref{fp:ex:i:b_psi_k} that
	\begin{align}
	\label{fp:ex:i:b_psi_ki}
	b(\mathbf{u}_{k_i})(\Xi[\psi_{k_i}], \theta) = l(\mathbf{u}_{k_i}, \psi_{k_i})(\theta) \quad \forall \theta \in \HkOS{1},\ \forall i \in \mathbb{N}.
	\end{align}
	On noting \eqref{fp:ex:i:cvg:psi_ki_weak}, \eqref{fp:ex:i:cvg:QZL_psiov_ki}, \eqref{fp:ex:i:cvg:u_ki_strong}, \eqref{fp:ex:i:cvg:u_ki_weak}, and \eqref{fp:ex:i:cvg:pot_psiov_ki}, we pass to the limit in \eqref{fp:ex:i:b_psi_ki} as $i \to \infty$ to deduce that
	\begin{align}
	\label{fp:ex:i:b_psihat}
	b(\hat{\uflow})(\hat{\psi}, \theta) = l(\hat{\uflow}, \psi)(\theta) \quad \forall \theta \in \HkOS{1}.
	\end{align}
	Now, we infer from the uniqueness of solutions to \eqref{fp:iteration:i} and \eqref{fp:iteration:ii} and from the comparison of \eqref{fp:ex:i:a_u} and \eqref{fp:ex:i:a_uhat} and of \eqref{fp:ex:i:b_psi} and \eqref{fp:ex:i:b_psihat} that
	\begin{align}
	\label{fp:ex:i:limit_id}
	\mathbf{u} = \hat{\uflow} \quad \text{ and } \quad \Xi[\psi] = \hat{\psi}.
	\end{align}
	Consequently, on noting \eqref{fp:ex:i:cvg:psi_ki_strong} and \eqref{fp:ex:i:limit_id}, we have that
	\begin{align}
	\label{fp:ex:i:psi_ki_limit}
	\Xi[\psi_{k_i}] \to \hat{\psi} = \Xi[\psi] \quad \text{ in } \LpOS{2} \text{ as } i \to \infty.
	\end{align}
	Therefore, we infer from the uniqueness of solutions to \eqref{fp:iteration:i} and \eqref{fp:iteration:ii} that we can find for any subsequence $\{ \Xi[\psi_{k_j}] \}_{j \in \mathbb{N}}$ of $\{ \Xi[\psi_k] \}_{k \in \mathbb{N}}$ another subsequence $\{ \Xi[\psi_{k_{j_i}}] \}_{i \in \mathbb{N}}$ of $\{ \Xi[\psi_{k_j]} \}_{j \in \mathbb{N}}$ such that
	\begin{align}
	\label{fp:ex:i:cvg:psi_kji}
	\Xi[\psi_{k_{j_i}}] \to \Xi[\psi] \quad \text{ in } \LpOS{2} \text{ as } i \to \infty.
	\end{align}
	Hence, we can return to the original series $\{ \Xi[\psi_k] \}_{k \in \mathbb{N}}$ in \eqref{fp:ex:i:cvg:psi_kji}, which proves \eqref{fp:ex:i:goal}.

	\textit{Proof of ii).}
	It follows directly from \eqref{fp:iteration:iii:u}, \eqref{fp:iteration:iii:psi}, and the compact embedding $\H^1(\Omega \times {\SMat^{d-1}})$ $\hookrightarrow\hookrightarrow$ $\L^2(\Omega \times {\SMat^{d-1}})$ that \textit{ii)} holds.

	\textit{Proof of iii).}
	Let $\lambda \in (0,1]$ and $\psi \in \L^2({\Omega\times\Sph})$ such that
	\begin{align}
	\label{fp:ex:iii:def_fp}
	\psi = \lambda \Xi(\psi).
	\end{align}
	We infer from Lemma \ref{lmm:fixed_point_iter} and the definition of $\Xi$ that there exists a unique $\mathbf{u} \in \Hflow_{0,\text{div}}^1(\Omega)$ such that
	\begin{subequations}
		\begin{alignat}{2}
		a(\psi)(\mathbf{u}, \mathbf{w}) &= k(\psi)(\mathbf{w}) \quad && \forall \mathbf{w} \in \Hflow_{0,\text{div}}^1(\Omega), \label{fp:ex:iii:eq:a} \\
		b(\mathbf{u})(\psi, \theta) &= \lambda l(\mathbf{u}, \psi)(\theta) \quad && \forall \theta \in \HkOS{1}. \label{fp:ex:iii:eq:b}
		\end{alignat}
	\end{subequations}
	In particular, it holds that $\psi \in \HkOS{1}$.
	Moreover, on choosing $\theta := \negat{\psi} \in \HkOS{1}$ as a test function in \eqref{fp:ex:iii:eq:b}, we obtain from $\negat{\psiL{n-1}} = 0$, $\negat{\text{Q}_0^\regL(\psi)}$ $= 0$, $\mathbf{u}$ being divergence-free in $\Omega$ and having zero trace on $\partial \Omega$ that
	\begin{align}
    \label{fp:ex:iii:nonneg}
    \begin{split}
	&\int_{\OmSph} \negat{\psi}^2 \dH\, \text{d}\xflow + \frac{\tau}{\DE} \int_{\OmSph} \left[ \mollEps^2 \abs{\nabla_\xflow \negat{\psi}}^2 + \abs{\nabla_\surfaceLtr \negat{\psi}}^2 \right] \dH\, \text{d}\xflow \\
	&\quad = \lambda \int_{\OmSph} \psiL{n-1} \negat{\psi} \dH\, \text{d}\xflow
	\leq 0.
    \end{split}
	\end{align}
	Hence, we have that $\psi \geq 0$ a.e. in ${\Omega\times\Sph}$, that is $\psi \in \Zp{2}$.
	Moreover, taking $\theta \equiv 1$ in \eqref{fp:ex:iii:eq:b} yields
	\begin{align}
	\label{fp:ex:iii:mass_consv}
		\int_{\OmSph} \psi \dH\, \text{d}\xflow = \lambda \int_{\OmSph} \psiL{n-1} \dH\, \text{d}\xflow.
	\end{align}
	Hence, we infer from \eqref{fp:ex:iii:mass_consv}, \eqref{est:nabla_pot_Linfty}, and $0 < \lambda \leq 1$ that
	\begin{align}
	\label{fp:ex:iii:est_pot_infty}
		\norm{ \pot[\psi + \psiL{n-1}] }_{\mathrm{W}^{1, \infty}({\Omega\times\Sph})} \leq C(\overline{\C}_\mollEps, \potStr, d) \norm{\psiL{n-1}}_{\LpOS{1}}.
	\end{align}

	The goal is to prove the existence of a constant $C^\star > 0$, which is independent of $\lambda$ and $\psi$, such that
	\begin{align}
	\label{fp:ex:iii:goal}
	\norm{\psi}_{\LpOS{2}} \leq C^\star.
	\end{align}

	One common ansatz to prove \eqref{fp:ex:iii:goal} would be taking $\mathbf{w} := \mathbf{u} \in \Hflow_{0,\text{div}}^1(\Omega)$ in \eqref{fp:ex:iii:eq:a} and $\theta := \psi \in \HkOS{1}$ in \eqref{fp:ex:iii:eq:b} as test functions to derive the desired bound. Then, combining \eqref{fp:iteration:iii:u} and \eqref{fp:iteration:iii:psi} yields
	\begin{align}
	\label{fp:ex:iii:fail}
	\norm{\psi}_{\LpOS{2}} \leq C(\tau, \text{L}) + C(\tau, \text{L}) \norm{\psi}_{\LpOS{2}} + C(\tau, \text{L}) \norm{\psi}_{\LpOS{2}}^2.
	\end{align}
	By choosing $\tau$ and $\text{L}$ in a clever way, we could absorb the second term on the right-hand side of \eqref{fp:ex:iii:fail} into the left-hand side. However, we could not control or absorb the last term on the right-hand side of \eqref{fp:ex:iii:fail}.
	This is the reason why we change the ansatz to mimicking the entropy estimate in Section \ref{sec:formal_entropy_estimate}. If we assume to control the $\text{L}$-regularized entropy of $\psi$ analogously to the formal estimate \eqref{formal:entropy_est}, we infer from \eqref{reg:FL_geq_ss} that
	\begin{align*}
	\frac{1}{4 \text{L}} \int_{\OmSph} \abs{\psi}^2 \dH\, \text{d}\xflow - C(\text{L}) \leq \int_{\OmSph} \entrop^\regL(\psi) \dH\, \text{d}\xflow \leq C,
	\end{align*}
	which then yields the desired bound \eqref{fp:ex:iii:goal}.
	Hence, our goal is to control the $\text{L}$-regularized entropy of $\psi$.

	As $\psi$ is nonnegative and there is no reason why $\psi$ should be strictly positive, we cannot apply $\lrcb{\entropL}^\prime(\cdot)$ to $\psi$ to define our test function as $\mu_{\psi, \text{L}} = \lrcb{\entropL}^\prime(\psi) + \pot[\psi]$. Hence, for $0 < \delta \ll 1$ let us define the $\delta$-regularized chemical potential by
	\begin{align}
	\label{fp:ex:iii:def:muPsidLl}
	\mu_{\psi, \regDelta, \regL, \lambda} := \lrcb{\entropL}^\prime(\psi + \delta) + \frac{\lambda}{2} \pot\left[ \psi + \psiL{n-1} \right] \in \HkOS{1}.
	\end{align}
	In particular, the following representation of the diffusive terms motivates the multiplication of the interaction potential in \eqref{fp:ex:iii:def:muPsidLl} by $\lambda$.
	First, we rewrite the diffusive terms in \eqref{fp:ex:iii:eq:b}. It holds that, for all $\theta \in \HkOS{1}$,
	\begin{align}
	\label{fp:ex:iii:diff_x_rewritten}
	\begin{split}
	& \int_{\OmSph} \left[ \nabla_\xflow \psi \cdot \nabla_\xflow \theta + \frac{\lambda}{2} \text{Q}^\regL(\psi) \nabla_\xflow \pot[\psi + \psiL{n-1}] \cdot \nabla_\xflow \theta \right] \dH\, \text{d}\xflow \\
	&\quad= \int_{\OmSph} \text{Q}^\regL(\psi + \delta) \nabla_\xflow \mu_{\psi, \regDelta, \regL, \lambda} \cdot \nabla_\xflow \theta \dH\, \text{d}\xflow \\
	&\quad\quad + \int_{\OmSph} \frac{\lambda}{2} \lrcb{ \text{Q}^\regL(\psi) - \text{Q}^\regL(\psi + \delta) } \nabla_\xflow \pot[\psi + \psiL{n-1}] \cdot \nabla_\xflow \theta \dH\, \text{d}\xflow.
	\end{split}
	\end{align}
	Similarly, we can rewrite the analogue term with respect to the surface gradient.
	Now, let us start deriving the entropy estimate.
	First, multiplying \eqref{fp:ex:iii:eq:a} by $\lambda$, taking $\mathbf{w} := \mathbf{u} \in \Hflow_{0,\text{div}}^1(\Omega)$ as a test function in \eqref{fp:ex:iii:eq:a}, and noting the formula
	\begin{align}
	\label{eq:aab}
	a(a-b) = \frac{1}{2} \lrcb{ a^2 + (a-b)^2 - b^2 } \qquad \text{ for all } a,b \in \mathbb{R},
	\end{align}
	it follows that
	\begin{align}
	\label{fp:ex:iii:u_tested_u}
	\begin{split}
	&\frac{\lambda \RE \DE}{2 \tau} \int_{\Omega} \left[ \abs{\mathbf{u}}^2 + \abs{\mathbf{u} - \uL{n-1}}^2 - \abs{\uL{n-1}}^2 \right] \, \text{d}\xflow + \lambda \viscFrac \DE \int_{\Omega} \abs{\nabla_\xflow \mathbf{u}}^2 \, \text{d}\xflow \\
	&\quad + \frac{\lambda (1 - \viscFrac) \DE}{2} \int_{\OmSph} \text{Q}^\regL(\psi) \lrcb{\nabla_\xflow\mathbf{u} : \mflow\otimes\mflow}^2 \dH(\mflow)\, \text{d}\xflow \\
	&= - \frac{\lambda(1-\viscFrac)}{2} \int_{\OmSph} \psi \nabla_\xflow \pot[\psi + \psiL{n-1}] \cdot \mathbf{u} \dH\, \text{d}\xflow \\
	&\quad - \frac{\lambda (1- \viscFrac)}{2} \int_{\OmSph} \text{Q}^\regL(\psi) \lrcb{ (\mathbb{I} - \mflow\otimes\mflow)\nabla_\xflow\mathbf{u}\mathbf{m} } \cdot \nabla_\surfaceLtr \pot[\psi + \psiL{n-1}] \dH(\mflow)\, \text{d}\xflow \\
	&\quad - \lambda(1 - \viscFrac) \int_{\OmSph} \psi (d\mflow\otimes\mflow - \mathbb{I}):\nabla_\xflow\mathbf{u} \dH(\mflow)\, \text{d}\xflow \\
	&\quad - \lambda\activity (1 - \viscFrac) \int_{\OmSph} \text{Q}^\regL(\psi) \mflow\otimes\mflow : \nabla_\xflow \mathbf{u} \dH(\mflow)\, \text{d}\xflow.
	\end{split}
	\end{align}

	Next, noting \eqref{fp:ex:iii:diff_x_rewritten}, the convexity of $\entrop^\regL$, and taking $\theta := \mu_{\psi, \regDelta, \regL, \lambda} \in \HkOS{1}$ as a test function in \eqref{fp:ex:iii:eq:b}, we have that
	\begin{align*}
	&\frac{1}{\tau} \int_{\OmSph} \left[ \entrop^\regL(\psi + \delta) - \entrop^\regL(\lambda \psiL{n-1} + \delta) \right] \dH\, \text{d}\xflow \\
	&\quad\quad + \frac{\lambda}{2 \tau} \int_{\OmSph} \lrcb{\psi - \lambda \psiL{n-1}} \lrcb{ \pot[\psi] + \pot[\psiL{n-1}] } \dH\, \text{d}\xflow \\
	&\quad\quad + \frac{\lambda \mollEps^2}{\DE} \int_{\OmSph} \text{Q}^\regL(\psi + \delta) \abs{\nabla_\xflow \mu_{\psi, \regDelta, \regL, \lambda}}^2 \dH\, \text{d}\xflow \\
	&\quad\quad + \frac{\lambda}{\DE} \int_{\OmSph} \text{Q}^\regL(\psi + \delta) \abs{\nabla_\surfaceLtr \mu_{\psi, \regDelta, \regL, \lambda}}^2  \dH\, \text{d}\xflow \\
	&\quad\leq \frac{\lambda}{2} \int_{\OmSph} \psi \nabla_\xflow \pot[\psi + \psiL{n-1}] \cdot \mathbf{u} \dH\, \text{d}\xflow \numberthis[fp:ex:iii:psi_tested_muPsidLl] \\
	&\quad\quad\quad + \lambda \int_{\OmSph} \text{Q}^\regL(\psi) \lrcb{(\mathbb{I} - \mflow\otimes\mflow) \nabla_\xflow \mathbf{u} \mathbf{m}} \cdot \nabla_\surfaceLtr \mu_{\psi, \regDelta, \regL, \lambda} \dH(\mflow)\, \text{d}\xflow \\
	&\quad\quad\quad - \frac{\lambda \mollEps^2}{2 \DE} \int_{\OmSph} \lrcb{ \text{Q}^\regL(\psi) - \text{Q}^\regL(\psi + \delta) } \nabla_\xflow \pot[\psi + \psiL{n-1}] \cdot \nabla_\xflow \mu_{\psi, \regDelta, \regL, \lambda} \dH\, \text{d}\xflow \\
	&\quad\quad\quad - \frac{\lambda}{2 \DE} \int_{\OmSph} \lrcb{\text{Q}^\regL(\psi) - \text{Q}^\regL(\psi + \delta)} \nabla_\surfaceLtr \pot[\psi + \psiL{n-1}] \cdot \nabla_\surfaceLtr \mu_{\psi, \regDelta, \regL, \lambda} \dH\, \text{d}\xflow.
	\end{align*}

	Before we care about the error terms in \eqref{fp:ex:iii:psi_tested_muPsidLl}, we exploit the natural structure of \eqref{fp:ex:iii:u_tested_u} and \eqref{fp:ex:iii:psi_tested_muPsidLl} as good as possible to see which terms cancel out.
	We infer from \eqref{eq:chen_stress_formula} that the second term on the right-hand side of \eqref{fp:ex:iii:psi_tested_muPsidLl} can be rewritten as
	\begin{align*}
	&\lambda \int_{\OmSph} \text{Q}^\regL(\psi) \lrcb{ (\mathbb{I} - \mflow\otimes\mflow) \nabla_\xflow \mathbf{u} \mathbf{m} } \cdot \nabla_\surfaceLtr \mu_{\psi, \regDelta, \regL, \lambda} \dH(\mflow)\, \text{d}\xflow\\
	&\quad= \frac{\lambda^2}{2} \int_{\OmSph} \text{Q}^\regL(\psi) \lrcb{ (\mathbb{I} - \mflow\otimes\mflow) \nabla_\xflow \mathbf{u} \mathbf{m} } \cdot \nabla_\surfaceLtr \pot[\psi + \psiL{n-1}] \dH(\mflow)\, \text{d}\xflow \\
	&\quad\quad + \lambda \int_{\OmSph} \text{Q}^\regL(\psi) \lrcb{ (\mathbb{I} - \mflow\otimes\mflow) \nabla_\xflow \mathbf{u} \mathbf{m} } \cdot \lrcb{ \frac{1}{\text{Q}^\regL(\psi + \delta)} \nabla_\surfaceLtr \psi } \dH(\mflow)\, \text{d}\xflow \\
	&\quad= \frac{\lambda^2}{2} \int_{\OmSph} \text{Q}^\regL(\psi) \lrcb{ (\mathbb{I} - \mflow\otimes\mflow) \nabla_\xflow \mathbf{u} \mathbf{m} } \cdot \nabla_\surfaceLtr \pot[\psi + \psiL{n-1}] \dH(\mflow)\, \text{d}\xflow \\
	&\quad\quad + \lambda \int_{\OmSph} \psi (d\mflow\otimes\mflow - \mathbb{I}) : \nabla_\xflow \mathbf{u} \dH(\mflow)\, \text{d}\xflow \numberthis[fp:ex:iii:chen_stress_rewritten] \\
	&\quad\quad + \lambda \int_{\OmSph} \frac{\text{Q}^\regL(\psi) - \text{Q}^\regL(\psi + \delta)}{\text{Q}^\regL(\psi + \delta)} \lrcb{ (\mathbb{I} - \mflow\otimes\mflow) \nabla_\xflow \mathbf{u} \mathbf{m} } \cdot \nabla_\surfaceLtr \psi \dH(\mflow)\, \text{d}\xflow.
	\end{align*}
	We note that the second term on the right-hand side of \eqref{fp:ex:iii:chen_stress_rewritten} cancels out with the third term on the right-hand side of \eqref{fp:ex:iii:u_tested_u}.
	Furthermore, we deduce from \eqref{mollifier:fg}, $0 < \lambda \leq 1$, and the nonnegativity of the Maier-Saupe interaction kernel \eqref{def:interaction_kernel} and of $\psi$ and $\psiL{n-1}$ a.e. in ${\Omega\times\Sph}$ that
	\begin{align}
	\label{fp:ex:iii:eq:pot_energy}
	\begin{split}
	&\frac{\lambda}{2\tau} \int_{\OmSph} \lrcb{\psi - \lambda \psiL{n-1}} \lrcb{ \pot[\psi] + \pot[\psiL{n-1}] } \dH\, \text{d}\xflow \\
	&\quad= \frac{\lambda}{2\tau} \int_{\OmSph} \left[ \psi\pot[\psi] - \lambda\psiL{n-1}\pot[\psiL{n-1}] \right] \dH\, \text{d}\xflow \\
	&\quad\quad + \frac{\lambda(1-\lambda)}{2 \tau} \int_{\OmSph} \psi \pot[\psiL{n-1}] \dH\, \text{d}\xflow \\
	&\quad\geq \frac{\lambda}{2\tau} \int_{\OmSph} \left[ \psi\pot[\psi] - \lambda\psiL{n-1}\pot[\psiL{n-1}] \right] \dH\, \text{d}\xflow.
	\end{split}
	\end{align}

	Now, we obtain from multiplying \eqref{fp:ex:iii:psi_tested_muPsidLl} by $(1 - \viscFrac)$ and from combining \eqref{fp:ex:iii:u_tested_u} and \eqref{fp:ex:iii:psi_tested_muPsidLl} - \eqref{fp:ex:iii:eq:pot_energy} that
	\begin{align*}
	&\frac{\lambda\RE\DE}{2 \tau} \int_{\Omega} \left[ \abs{\mathbf{u}}^2 + \abs{\mathbf{u} - \uL{n-1}}^2 - \abs{\uL{n-1}}^2 \right] \, \text{d}\xflow + \lambda \viscFrac \DE \int_{\Omega} \abs{\nabla_\xflow \mathbf{u}}^2 \, \text{d}\xflow \\
	&\quad + \frac{\lambda (1 - \viscFrac) \DE}{2} \int_{\OmSph} \text{Q}^\regL(\psi) \lrcb{\nabla_\xflow\mathbf{u} : \mflow\otimes\mflow}^2 \dH(\mflow)\, \text{d}\xflow \\
	&\quad + \frac{1 - \viscFrac}{\tau} \int_{\OmSph} \left[ \entrop^\regL(\psi + \delta) - \entrop^\regL(\lambda \psiL{n-1} + \delta) \right] \dH\, \text{d}\xflow \\
	&\quad + \frac{\lambda (1 - \viscFrac)}{2 \tau} \int_{\OmSph} \left[ \psi \pot[\psi] - \lambda \psiL{n-1} \pot[\psiL{n-1}] \right] \dH\, \text{d}\xflow \\
	&\quad + \frac{\lambda (1 - \viscFrac) \mollEps^2}{\DE} \int_{\OmSph} \text{Q}^\regL(\psi + \delta) \abs{\nabla_\xflow \mu_{\psi, \regDelta, \regL, \lambda}}^2 \dH\, \text{d}\xflow \\
	&\quad + \frac{\lambda (1 - \viscFrac)}{\DE} \int_{\OmSph} \text{Q}^\regL(\psi + \delta) \abs{\nabla_\surfaceLtr \mu_{\psi, \regDelta, \regL, \lambda}}^2 \dH\, \text{d}\xflow \\
	&\leq - \lambda \alpha (1 - \viscFrac) \int_{\OmSph} \text{Q}^\regL(\psi) \mflow\otimes\mflow : \nabla_\xflow \mathbf{u} \dH(\mflow)\, \text{d}\xflow \numberthis[fp:ex:iii:u_psi_combined] \\
	&\quad + \frac{\lambda(\lambda-1)(1-\viscFrac)}{2} \int_{\OmSph} \text{Q}^\regL(\psi) \lrcb{ (\mathbb{I} - \mflow\otimes\mflow) \nabla_\xflow \mathbf{u} \mathbf{m} } \cdot  \\
	&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad \cdot \nabla_\surfaceLtr \pot[\psi + \psiL{n-1}] \dH(\mflow)\, \text{d}\xflow \\
	&\quad + \lambda(1 - \viscFrac) \int_{\OmSph} \frac{\text{Q}^\regL(\psi) - \text{Q}^\regL(\psi + \delta)}{\text{Q}^\regL(\psi + \delta)} \lrcb{ (\mathbb{I} - \mflow\otimes\mflow) \nabla_\xflow \mathbf{u} \mathbf{m} } \cdot \\
	&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad \cdot \nabla_\surfaceLtr \psi \dH(\mflow)\, \text{d}\xflow \\
	&\quad - \frac{\lambda \mollEps^2 (1 - \viscFrac)}{2 \DE} \int_{\OmSph} \lrcb{ \text{Q}^\regL(\psi) - \text{Q}^\regL(\psi + \delta) } \nabla_\xflow \pot[\psi + \psiL{n-1}] \cdot  \\
	&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad \cdot \nabla_\xflow \mu_{\psi, \regDelta, \regL, \lambda} \dH\, \text{d}\xflow \\
	&\quad - \frac{\lambda (1 - \viscFrac)}{2 \DE} \int_{\OmSph} \lrcb{\text{Q}^\regL(\psi) - \text{Q}^\regL(\psi + \delta)} \nabla_\surfaceLtr \pot[\psi + \psiL{n-1}] \cdot  \\
	&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad \cdot \nabla_\surfaceLtr \mu_{\psi, \regDelta, \regL, \lambda} \dH\, \text{d}\xflow \\
	&=: I + II + III + IV + V.
	\end{align*}
	We note that the active term $I$ appears in the same form as in the formal estimate \eqref{formal:summarize_1}. Moreover, it is worth mentioning that in the case of the usual entropy estimate, that is $\lambda = 1$, we would have that $II = 0$.
	The error terms $III-V$ are caused by the approximate scheme and the $\delta$-regularized test function $\mu_{\psi, \regDelta, \regL, \lambda}$.\\
	Let us now continue estimating the terms $I - V$. We infer from $0 < \lambda \leq 1$, \eqref{reg:QL_leq}, \eqref{fp:ex:iii:mass_consv}, and Young's inequality that
	\begin{align*}
	\abs{I}
	&\leq \frac{\lambda (1 - \viscFrac) \DE}{4} \int_{\OmSph} \text{Q}^\regL(\psi) \lrcb{\nabla_\xflow\mathbf{u} : \mflow\otimes\mflow}^2 \dH(\mflow)\, \text{d}\xflow \\
	&\quad + \frac{(1-\viscFrac) \activity^2}{\DE} \norm{\psiL{n-1}}_{\LpOS{1}}. \numberthis[fp:ex:iii:I]
	\end{align*}
	The first term on the right-hand side of \eqref{fp:ex:iii:I} can be absorbed into the left-hand side of \eqref{fp:ex:iii:u_psi_combined}.
	Furthermore, we deduce from $\abs{\mathbb{I} - \mflow\otimes\mflow} \leq 2d$, $\abs{\lambda(\lambda-1)} \leq \lambda$, \eqref{reg:QL_leq}, \eqref{fp:ex:iii:mass_consv}, \eqref{fp:ex:iii:est_pot_infty}, Jensen's inequality, and Young's inequality that
	\begin{align}
	\label{fp:ex:iii:II}
	\begin{split}
	\abs{II}
	&\leq \frac{\lambda (1-\viscFrac) 2d }{2} \int_{\OmSph} \abs{\text{Q}^\regL(\psi)} \abs{\nabla_\xflow \mathbf{u}} \abs{\nabla_\surfaceLtr \pot[\psi + \psiL{n-1}]} \dH\, \text{d}\xflow \\
	&\leq \frac{\lambda \gamma \DE}{2} \int_{\Omega} \abs{\nabla_\xflow \mathbf{u}}^2 \, \text{d}\xflow  + C_1 \text{L} \norm{\psiL{n-1}}_{\LpOS{1}}^3,
	\end{split}
	\end{align}
	where $C_1 = C_1(\viscFrac, \viscFrac^{-1}, \DE^{-1}, d, {\Omega\times\Sph}, \overline{\C}_\mollEps, \potStr) > 0$ is independent of $\delta, \tau$, $\lambda$, $\text{L}$, and $\psi$.
	Again, we can absorb the first term on the right-hand side of \eqref{fp:ex:iii:II} into the left-hand side of \eqref{fp:ex:iii:u_psi_combined}.
	Moreover, on noting the Lipschitz continuity of $\text{Q}^\regL$, the nonnegativity of $\psi$, \eqref{reg:QL_leq}, \eqref{fp:ex:iii:mass_consv}, \eqref{fp:ex:iii:est_pot_infty}, $0 < \lambda \leq 1$, and on adapting a computation of Chen and Liu (cf. \cite[(3.32)]{chen_liu_2013_entropy_solution}), we infer from Jensen's inequality and Young's inequality that
	\begin{align*}
	\abs{III}
	&\leq \lambda (1 - \viscFrac) 2d \sqrt{\delta} \int_{\OmSph} \abs{\nabla_\xflow\mathbf{u}} \abs{\nabla_\surfaceLtr \psi} \frac{1}{\sqrt{\text{Q}^\regL(\psi + \delta)}} \dH\, \text{d}\xflow \\
	&\leq \frac{\lambda (1 - \viscFrac)}{8 \DE} \int_{\OmSph} \text{Q}^\regL(\psi + \delta) \abs{\nabla_\surfaceLtr \lrcb{\entropL}^\prime(\psi + \delta)}^2 \dH\, \text{d}\xflow \\
	&\quad + C_2(\DE, \viscFrac, d) \delta \int_{\Omega} \abs{\nabla_\xflow \mathbf{u}}^2 \, \text{d}\xflow \\
	&\leq \frac{\lambda (1 - \viscFrac)}{4 \DE} \int_{\OmSph} \text{Q}^\regL(\psi + \delta) \abs{\nabla_\surfaceLtr \mu_{\psi, \regDelta, \regL, \lambda}}^2 \dH\, \text{d}\xflow \numberthis[fp:ex:iii:III] \\
	&\quad + C_3\lrcb{\DE^{-1}, \viscFrac, \overline{\C}_\mollEps, \potStr, d, \Omega} \norm{\psiL{n-1}}_{\LpOS{1}}^2 \lrcb{\norm{\psiL{n-1}}_{\LpOS{1}} + 1} \\
	&\quad + C_2(\DE, \viscFrac, d) \delta \int_{\Omega} \abs{\nabla_\xflow \mathbf{u}}^2 \, \text{d}\xflow,
	\end{align*}
	where $C_2 = C_2(\DE, \viscFrac, d) > 0$ and $C_3 =C_3\lrcb{\DE^{-1}, \viscFrac, \overline{\C}_\mollEps, \potStr, d, \Omega} > 0$ are independent of $\delta, \tau$, $\lambda$, $\text{L}$, and $\psi$.
	We note that the first term on the right-hand side of \eqref{fp:ex:iii:III} can be absorbed into the left-hand side of \eqref{fp:ex:iii:u_psi_combined}. Moreover, the last term on the right-hand side of \eqref{fp:ex:iii:III} can be absorbed as well if we choose $\delta > 0$ sufficiently small. The exact range of $\delta$ will be specified later.
	Next, we infer from $\text{Q}^\regL(s) \leq \text{Q}^\regL(s + \delta)$ for all $s \geq 0$, the Lipschitz continuity of $\text{Q}^\regL$, \eqref{reg:QL_leq}, \eqref{fp:ex:iii:mass_consv}, \eqref{fp:ex:iii:est_pot_infty}, $0 < \lambda \leq 1$, Jensen's inequality, and Young's inequality that
	\begin{align*}
	\abs{IV}
	&\leq \frac{\lambda \mollEps^2 (1 - \viscFrac)}{2 \DE} \int_{\OmSph} \text{Q}^\regL(\psi + \delta) \abs{\nabla_\xflow \mu_{\psi, \regDelta, \regL, \lambda}}^2 \dH\, \text{d}\xflow \numberthis[fp:ex:iii:IV] \\
	&\quad + C_4 \lrcb{\mollEps, \DE^{-1}, \viscFrac, \Omega, \overline{\C}_\mollEps, \potStr, d} \delta \norm{\psiL{n-1}}_{\LpOS{1}}^2
	\end{align*}
	and
	\begin{align*}
	\abs{V}
	&\leq \frac{\lambda (1 - \viscFrac)}{4 \DE} \int_{\OmSph} \text{Q}^\regL(\psi + \delta) \abs{\nabla_\surfaceLtr \mu_{\psi, \regDelta, \regL, \lambda}}^2 \dH\, \text{d}\xflow \numberthis[fp:ex:iii:V] \\
	&\quad + C_5\lrcb{\DE^{-1}, \viscFrac, \Omega, \overline{\C}_\mollEps, \potStr, d} \delta \norm{\psiL{n-1}}_{\LpOS{1}}^2,
	\end{align*}
	where $C_4 = C_4 \lrcb{\mollEps, \DE^{-1}, \viscFrac, \Omega, \overline{\C}_\mollEps, \potStr, d} > 0$ and $C_5 = C_5\lrcb{\DE^{-1}, \viscFrac, \Omega, \overline{\C}_\mollEps, \potStr, d} > 0$ are independent of $\delta, \tau$, $\lambda$, $\text{L}$, and $\psi$.
	Now, we obtain from combining \eqref{fp:ex:iii:u_psi_combined} - \eqref{fp:ex:iii:V} that
	\begin{align*}
	&\frac{\lambda \RE \DE}{2 \tau} \int_{\Omega} \left[ \abs{\mathbf{u}}^2 + \abs{\mathbf{u} - \uL{n-1}}^2 - \abs{\uL{n-1}}^2 \right] \, \text{d}\xflow + \lrcb{\frac{\lambda\viscFrac\DE}{2} - C_2 \delta} \int_{\Omega} \abs{\nabla_\xflow \mathbf{u}}^2 \, \text{d}\xflow \\
	&\quad\quad + \frac{\lambda (1 - \viscFrac) \DE}{4} \int_{\OmSph} \text{Q}^\regL(\psi) \lrcb{\nabla_\xflow\mathbf{u} : \mflow\otimes\mflow}^2 \dH(\mflow)\, \text{d}\xflow \\
	&\quad\quad + \frac{1 - \viscFrac}{\tau} \int_{\OmSph} \left[ \entrop^\regL(\psi + \delta) - \entrop^\regL(\lambda\psiL{n-1} + \delta) \right] \dH\, \text{d}\xflow \\
	&\quad\quad + \frac{\lambda (1 - \viscFrac)}{2 \tau} \int_{\OmSph} \left[ \psi\pot[\psi] - \lambda\psiL{n-1}\pot[\psiL{n-1}] \right] \dH\, \text{d}\xflow \\
	&\quad\quad + \frac{\lambda (1 - \viscFrac) \mollEps^2}{2 \DE} \int_{\OmSph} \text{Q}^\regL(\psi + \delta) \abs{\nabla_\xflow \mu_{\psi, \regDelta, \regL, \lambda}}^2 \dH\, \text{d}\xflow \\
	&\quad\quad + \frac{\lambda (1 - \viscFrac)}{2 \DE} \int_{\OmSph} \text{Q}^\regL(\psi + \delta) \abs{\nabla_\surfaceLtr \mu_{\psi, \regDelta, \regL, \lambda}}^2 \dH\, \text{d}\xflow \\
	&\quad\leq \frac{(1-\viscFrac) \activity^2}{\DE} \norm{\psiL{n-1}}_{\LpOS{1}} +  C_1 \text{L} \norm{\psiL{n-1}}_{\LpOS{1}}^3 \numberthis[fp:ex:iii:u_psi_I_to_V] \\
	&\quad\quad + C_3 \lrcb{\norm{\psiL{n-1}}_{\LpOS{1}} + 1}\norm{\psiL{n-1}}_{\LpOS{1}}^2 \\
	&\quad\quad + \lrcb{C_4 + C_5 }\delta \norm{\psiL{n-1}}_{\LpOS{1}}^2.
	\end{align*}

    Now, we deduce from \eqref{fp:ex:iii:u_psi_I_to_V}, $0 < \lambda \leq 1$, \eqref{reg:FL_geq_ss}, the nonnegativity of $\psi$, omitting some nonnegative terms, and the Dominated Convergence Theorem that
	\begin{align}
	\label{fp:ex:iii:bound_psi_psi_1}
	\begin{split}
	\frac{1}{4 \text{L}} \int_{\OmSph} \abs{\psi}^2 \dH\, \text{d}\xflow
	 \leq C,
	\end{split}
	\end{align}
	where $C = C\lrcb{C_1, C_3, \potStr, d, \Omega,  \RE, \DE, \DE^{-1}, \activity, \viscFrac, \tau, \text{L}, \uL{n-1}, \psiL{n-1}} > 0$ is independent of $\lambda$ and $\psi$.
	Hence, it finally follows that \eqref{fp:ex:iii:goal} holds. The proof of Lemma \ref{lmm:existence_fp} is thereby complete.
\end{proof}

\begin{crllr}
	\label{crllr:nonnegativity}
	Suppose that the assumptions \eqref{assumptions:domain_init_data} hold.
	Then, there exist solutions $\{ ( \uL{k}, \psiL{k} ) \}_{k=1}^N$ $\subset$ $\Hflow_{0,\text{div}}^1(\Omega) \times \HkOS{1}$ to \eqref{def:PLtau}. In particular, it holds that $\psiL{k} \in \Zp{2}$ for $k \in \{1, \dots, N\}$. Hence, $\{ ( \uL{k}, \psiL{k} ) \}_{k=1}^N$ $\subset$ $\Hflow_{0,\text{div}}^1(\Omega) \times (\HkOS{1} \cap \Zp{2})$ is a solution to $\lrcb{\mathcal{P}_{\regL}^\tau}$.
\end{crllr}

\begin{proof}
	As a fixed-point of $\Xi$ is equivalent to a solution to \eqref{def:PLtau}, applying Lemma \ref{lmm:existence_fp} iteratively yields the existence of solutions $\{ ( \uL{k}, \psiL{k} ) \}_{k=1}^N$ $\subset$ $\Hflow_{0,\text{div}}^1(\Omega) \times \H^1(\Omega \times {\SMat^{d-1}})$ to \eqref{def:PLtau}.
	It therefore remains to show that $\{\psiL{k}\}_{k=1}^N \subset \Zp{2}$.
	As $\psiL{0} \in \Zp{2}$ (see \eqref{init:def:psi_L_0}), let $k \in \{1, \dots, N\}$ and $\psiL{k-1} \in \Zp{2}$.
    Applying the same arguments as in \eqref{fp:ex:iii:nonneg} yields $\psiL{k} \in \Zp{2}$.
\end{proof}

\subsection{A uniform entropy estimate}
\label{sec:uniform_entropy_estimate}
This section is devoted to the derivation of a uniform entropy estimate with respect to $N$ and $\text{L}$, which mimics the formal entropy estimate \eqref{formal:entropy_est}.
First, we infer from taking $\psi \equiv 1$ as a test function in \eqref{eq:PLtau:psi} that the mass of $\psiL{n}$ is conserved, that is
\begin{align}
\label{eq:mass_consv_psiLn}
\int_{\OmSph} \psiL{n} \dH\, \text{d}\xflow = \int_{\OmSph} \psiL{n-1} \dH\, \text{d}\xflow = \int_{\OmSph} \psiL{0} \dH\, \text{d}\xflow
\end{align}
for all $n \in \{1, \dots, N\}$.
Moreover, for future references, we deduce from \eqref{est:nabla_pot_Linfty}, \eqref{eq:mass_consv_psiLn}, and \eqref{init:regularity_psi_L_0} that, for all $k \in \{1, \dots, N\}$,
\begin{align}
\label{eq:nabla_potLk_bound}
\begin{split}
\norm{\potL{k}}_{\mathrm{W}^{1, \infty}({\Omega\times\Sph})}
\leq C(\overline{\C}_\mollEps, \potStr, d) \norm{\psi_0}_{\LpOS{1}}
< \infty,
\end{split}
\end{align}
where $C = C(\overline{\C}_\mollEps, \potStr, d) > 0$ is independent of $\delta, \tau$, and $\text{L}$.
The following lemma holds:

\begin{lmm}[Discrete-in-time entropy estimate]
	\label{lmm:time_discrete_entropy}
	Suppose that the assumptions \eqref{assumptions:domain_init_data} hold.
	Then, solutions $\{ ( \uL{k}, \psiL{k} ) \}_{k=1}^N$ $\subset$ $\Hflow_{0,\text{div}}^1(\Omega) \times (\HkOS{1} \cap \Zp{2})$ to \eqref{def:PLtau} satisfy
	\begin{align*}
	& \frac{\RE \DE}{2} \int_{\Omega} \abs{ \uL{n} }^2 \, \text{d}\xflow + \frac{\RE\DE}{2} \sum_{k=1}^{n} \int_{\Omega} \abs{\uL{k} - \uL{k-1}}^2 \, \text{d}\xflow + \frac{\viscFrac \DE}{2} \sum_{k=1}^{n} \tau \int_{\Omega} \abs{\nabla_\xflow \uL{k}}^2 \, \text{d}\xflow \\
	& \quad + \frac{(1 - \viscFrac) \DE}{4} \sum_{k=1}^{n} \tau \int_{\OmSph} \text{Q}^\regL(\psiL{k}) \lrcb{ \nabla_\xflow \uL{k} : \mflow\otimes\mflow }^2 \dH(\mflow) \, \text{d}\xflow \\
	& \quad + (1 - \viscFrac) \int_{\OmSph} \mathcal{F}(\psiL{n}) \dH\, \text{d}\xflow + \frac{1 - \viscFrac}{2} \int_{\OmSph} \psiL{n} \pot[\psiL{n}] \dH \, \text{d}\xflow \\
	& \quad + \frac{1 - \viscFrac}{2 \text{L}} \sum_{k=1}^n \int_{\OmSph} \abs{ \psiL{k} - \psiL{k-1} }^2 \dH\, \text{d}\xflow \\
	& \quad + \frac{1 - \viscFrac}{\DE} \sum_{k=1}^{n} \tau \int_{\OmSph} \Big[ \mollEps^2 | \nabla_\xflow \sqrt{\psiL{k}} |^2 + | \nabla_\surfaceLtr \sqrt{\psiL{k}} |^2 \Big] \dH\, \text{d}\xflow \\
	&\leq \frac{\RE \DE}{2} \int_{\Omega} \abs{\mathbf{u}_0}^2 \, \text{d}\xflow + (1 - \viscFrac) \int_{\OmSph} \mathcal{F}(\psi_0) \dH\, \text{d}\xflow \numberthis[est:discrete_entropy] \\
	&\quad + \frac{1 - \viscFrac}{2} \int_{\OmSph} \psi_0 \pot[\psi_0] \dH\, \text{d}\xflow + \frac{\activity^2 (1 - \viscFrac) T}{\DE} \norm{\psi_0}_{\LpOS{1}} + C_{*} \\
	&\leq C\lrcb{ \RE, \DE, \DE^{-1}, \viscFrac, \activity, T, \overline{\C}_\mollEps, \potStr, d, \Omega, \mollEps, \mathbf{u}_0, \psi_0 }
	< \infty \qquad \text{ for all } n \in \{ 1, \dots, N \},
	\end{align*}
	where $C_{*} > 0$ depends solely on $T, \psi_0, \DE^{-1}, \viscFrac, \overline{\C}_\mollEps, \potStr, d, \Omega$, and $\mollEps$.
\end{lmm}

\begin{proof}[Proof of Lemma \ref{lmm:time_discrete_entropy}]
	Let a regularization parameter $0 < \delta \ll 1$ be given. Then, we rewrite the diffusive terms in \eqref{eq:PLtau:psi} for all $\theta \in \HkOS{1}$ analogously to \eqref{fp:ex:iii:diff_x_rewritten} as follows:
	\begin{align}
    \label{discr_entr:diff_rewritten}
    \begin{split}
	& \frac{\mollEps^2}{\DE} \int_{\OmSph} \lrcb{\nabla_\xflow \psiL{n} + \text{Q}^\regL(\psiL{n}) \nabla_\xflow \potL{n}} \cdot \nabla_\xflow \theta \dH\, \text{d}\xflow \\
	&\quad = \frac{\mollEps^2}{\DE} \int_{\OmSph} \text{Q}^\regL(\psiL{n} + \delta) \nabla_\xflow \muPsidL{n} \cdot \nabla_\xflow \theta \dH\, \text{d}\xflow \\
	&\quad \quad + \frac{\mollEps^2}{\DE} \int_{\OmSph} \lrcb{\text{Q}^\regL(\psiL{n}) - \text{Q}^\regL(\psiL{n} + \delta)} \nabla_\xflow \potL{n} \cdot \nabla_\xflow \theta \dH\, \text{d}\xflow,
    \end{split}
	\end{align}
	where we introduced the $\delta$-regularized chemical potential
	\begin{align*}
	\muPsidL{n} := \lrcb{\entropL}^\prime(\psiL{n} + \delta) + \potL{n} \quad \in \HkOS{1}.
	\end{align*}

    Moreover, we recall the following estimate, which has been derived by Barrett and S\"uli by Taylor's series expansion (cf. \cite[p. 1238]{barrett_sueli_2011_polymers_I}):
	\begin{align}
    \label{discr_entr:taylor_est}
    \begin{split}
	&\int_{\OmSph} \frac{\psiL{n} - \psiL{n-1}}{\tau} \lrcb{\entropL}^\prime(\psiL{n} + \delta) \dH\, \text{d}\xflow \\
	&\quad\geq \frac{1}{\tau} \int_{\OmSph} \left[ \entrop^\regL(\psiL{n} + \delta) - \entrop^\regL(\psiL{n-1} + \delta) \right] \dH\, \text{d}\xflow \\
	&\quad\quad + \frac{1}{2 \tau \text{L}}\int_{\OmSph} \abs{ \psiL{n} - \psiL{n-1} } \dH\, \text{d}\xflow.
	\end{split}
    \end{align}

    Now, on noting \eqref{discr_entr:diff_rewritten} and \eqref{discr_entr:taylor_est}, we infer from taking $\mathbf{w} := \uL{n} \in \Hflow_{0,\text{div}}^1(\Omega)$ and $\theta := \muPsidL{n} \in \H^1(\Omega \times {\SMat^{d-1}})$ as test functions in \eqref{eq:PLtau:u} and \eqref{eq:PLtau:psi}, respectively, and from the analogue estimate \eqref{fp:ex:iii:u_psi_combined} with $\lambda = 1$ that the following estimate holds:
	\begin{align*}
	&\frac{\RE\DE}{2 \tau} \int_{\Omega} \left[ \abs{\uL{n}}^2 + \abs{\uL{n} - \uL{n-1}}^2 - \abs{\uL{n-1}}^2 \right] \, \text{d}\xflow + \viscFrac\DE \int_{\Omega} \abs{\nabla_\xflow \uL{n}}^2 \, \text{d}\xflow \\
	&\qquad + \frac{(1-\viscFrac) \DE}{2} \int_{\OmSph} \text{Q}^\regL(\psiL{n}) \lrcb{\nabla_\xflow \uL{n} : \mflow\otimes\mflow}^2 \dH(\mflow)\, \text{d}\xflow \\
	&\qquad + \frac{1-\viscFrac}{\tau} \int_{\OmSph} \left[ \entrop^\regL(\psiL{n} + \delta) - \entrop^\regL(\psiL{n-1} + \delta) \right] \dH\, \text{d}\xflow \\
	&\qquad + \frac{1-\viscFrac}{2 \tau} \int_{\OmSph} \left[ \psiL{n} \pot[\psiL{n}] - \psiL{n-1} \pot[\psiL{n-1}] \right]\dH\, \text{d}\xflow \\
	&\qquad + \frac{1-\viscFrac}{2 \tau \text{L}} \int_{\OmSph} \abs{ \psiL{n} - \psiL{n-1} }^2 \dH\, \text{d}\xflow \\
	&\qquad + \frac{(1-\viscFrac) \mollEps^2}{\DE} \int_{\OmSph} \text{Q}^\regL(\psiL{n} + \delta) \abs{\nabla_\xflow \muPsidL{n}}^2 \dH\, \text{d}\xflow \\
	&\qquad + \frac{1-\viscFrac}{\DE} \int_{\OmSph} \text{Q}^\regL(\psiL{n} + \delta) \abs{\nabla_\surfaceLtr \muPsidL{n}}^2 \dH\, \text{d}\xflow \\
	&\quad\leq - \activity (1 - \viscFrac) \int_{\OmSph} \text{Q}^\regL(\psiL{n}) \mflow\otimes\mflow : \nabla_\xflow \uL{n} \dH(\mflow)\, \text{d}\xflow \numberthis[discr_entr:u_psi_combined] \\
	&\quad\quad + (1-\viscFrac) \int_{\OmSph} \frac{\text{Q}^\regL(\psiL{n}) - \text{Q}^\regL(\psiL{n} + \delta)}{\text{Q}^\regL(\psiL{n} + \delta)} \lrcb{ (\mathbb{I} - \mflow\otimes\mflow) \nabla_\xflow \uL{n} \mathbf{m} } \cdot  \\
	&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad \cdot \nabla_\surfaceLtr \psiL{n} \dH(\mflow)\, \text{d}\xflow \\
	&\quad\quad - \frac{(1-\viscFrac)\mollEps^2}{\DE} \int_{\OmSph} \lrcb{ \text{Q}^\regL(\psiL{n}) - \text{Q}^\regL(\psiL{n} + \delta) } \nabla_\xflow \potL{n} \cdot \nabla_\xflow \muPsidL{n} \dH\, \text{d}\xflow \\
	&\quad\quad - \frac{1-\viscFrac}{\DE} \int_{\OmSph} \lrcb{ \text{Q}^\regL(\psiL{n}) - \text{Q}^\regL(\psiL{n} + \delta) } \nabla_\surfaceLtr \potL{n} \cdot \nabla_\surfaceLtr \muPsidL{n} \dH\, \text{d}\xflow \\
	&\quad=: I + II + III + IV.
	\end{align*}

	We infer from \eqref{reg:QL_leq}, the nonnegativity of $\psiL{n}$, and \eqref{eq:mass_consv_psiLn} that
	\begin{align}
	\label{discr_entr:I}
	\begin{split}
	\abs{I}
	&\leq \frac{(1-\viscFrac)\DE}{4} \int_{\OmSph} \text{Q}^\regL(\psiL{n}) \lrcb{\nabla_\xflow\uL{n} : \mflow\otimes\mflow}^2 \dH(\mflow)\, \text{d}\xflow \\
	&\quad + \frac{\activity^2 (1 - \viscFrac)}{\DE} \int_{\OmSph} \psiL{0} \dH\, \text{d}\xflow.
	\end{split}
	\end{align}
	The first term on the right-hand side of \eqref{discr_entr:I} can be absorbed into the left-hand side of \eqref{discr_entr:u_psi_combined}.
	Next, similarly to \eqref{fp:ex:iii:III}, we deduce from the Lipschitz continuity of $\text{Q}^\regL$, the nonnegativity of $\psiL{n}$, \eqref{eq:nabla_potLk_bound}, \eqref{reg:QL_leq}, \eqref{eq:mass_consv_psiLn}, $0 < \delta < 1$, and Young's inequality that
	\begin{align}
    \label{discr_entr:II}
    \begin{split}
	\abs{II}
	&\leq \frac{1-\viscFrac}{4 \DE} \int_{\OmSph} \text{Q}^\regL(\psiL{n} + \delta) \abs{\nabla_\surfaceLtr \muPsidL{n}}^2 \dH\, \text{d}\xflow  \\
	&\quad + C_2\lrcb{\DE^{-1}, \viscFrac, \overline{\C}_\mollEps, \potStr, d, \Omega} \norm{\psiL{0}}_{\LpOS{1}}^2 \lrcb{ 1 + \norm{\psiL{0}}_{\LpOS{1}} } \\
	&\quad + C_1(\DE, \viscFrac, d) \delta \int_{\Omega} \abs{\nabla_\xflow \uL{n}}^2 \, \text{d}\xflow,
    \end{split}
	\end{align}
	where $C_1 = C_1(\DE, \viscFrac, d) > 0$ and $C_2 = C_2\lrcb{\DE^{-1}, \viscFrac, \overline{\C}_\mollEps, \potStr, d, \Omega} > 0$ are independent of $\delta$, $\tau$, and $\text{L}$.
	Furthermore, on noting $\text{Q}^\regL(s) \leq \text{Q}^\regL(s + \delta)$ for all $s \geq 0$, the Lipschitz continuity of $\text{Q}^\regL$, the nonnegativity of $\psiL{n}$, \eqref{eq:nabla_potLk_bound}, and Young's inequality, we have that
	\begin{align}
	\label{discr_entr:III}
	\begin{split}
	\abs{III}
	&\leq \frac{\mollEps^2 (1-\viscFrac)}{2 \DE} \int_{\OmSph} \text{Q}^\regL(\psiL{n} + \delta) \abs{\nabla_\xflow \muPsidL{n}}^2 \dH\, \text{d}\xflow \\
	&\quad + C_3\lrcb{ \DE^{-1}, \mollEps, \viscFrac, \Omega, \overline{\C}_\mollEps, \potStr, d } \delta \norm{\psiL{0}}_{\LpOS{1}}^2,
	\end{split}
	\end{align}
	where $C_3 = C_3\lrcb{ \DE^{-1}, \mollEps, \viscFrac, \Omega, \overline{\C}_\mollEps, \potStr, d } > 0$ is independent of $\delta$, $\tau$, and $\text{L}$.
	Analogously to \eqref{discr_entr:III}, it follows that
	\begin{align*}
	\abs{IV}
	&\leq \frac{1-\viscFrac}{4 \DE} \int_{\OmSph} \text{Q}^\regL(\psiL{n} + \delta) \abs{\nabla_\surfaceLtr \muPsidL{n}}^2 \dH\, \text{d}\xflow \\
	&\quad + C_4\lrcb{\DE^{-1}, \viscFrac, \Omega, \overline{\C}_\mollEps, \potStr, d} \delta \norm{\psiL{0}}_{\LpOS{1}}^2, \numberthis[discr_entr:IV]
	\end{align*}
	where $C_4 = C_4\lrcb{\DE^{-1}, \viscFrac, \Omega, \overline{\C}_\mollEps, \potStr, d} > 0$ is independent of $\delta$, $\tau$, and $\text{L}$.
	Now, on combining \eqref{discr_entr:u_psi_combined}-\eqref{discr_entr:IV}, restricting the range of $\delta$ to
	\begin{align*}
	0 < \delta < \frac{\viscFrac \DE}{2 C_1},
	\end{align*}
	and applying a discrete integration in time, we have that
	\begin{align*}
	&\frac{\RE \DE}{2} \int_{\Omega} \abs{\uL{n}}^2 \, \text{d}\xflow + \frac{\RE \DE}{2} \sum_{k=1}^{n} \int_{\Omega} \abs{\uL{k} - \uL{k-1}}^2 \, \text{d}\xflow + \frac{\viscFrac \DE}{2} \sum_{k=1}^{n} \tau \int_{\Omega} \abs{\nabla_\xflow \uL{k}}^2 \, \text{d}\xflow \\
	& \qquad + \frac{(1-\viscFrac) \DE}{4} \sum_{k=1}^{n} \tau \int_{\OmSph} \text{Q}^\regL(\psiL{k}) \lrcb{ \nabla_\xflow \uL{k} : \mflow\otimes\mflow }^2 \dH(\mflow)\, \text{d}\xflow \\
	& \qquad + (1-\viscFrac) \int_{\OmSph} \entrop^\regL(\psiL{n} + \delta) \dH\, \text{d}\xflow + \frac{1-\viscFrac}{2} \int_{\OmSph} \psiL{n} \pot[\psiL{n}] \dH\, \text{d}\xflow \\
	&\qquad + \frac{1-\viscFrac}{2 \text{L}} \sum_{k=1}^n \int_{\OmSph} \abs{ \psiL{k} - \psiL{k-1} }^2 \dH\, \text{d}\xflow \\
	& \qquad + \frac{(1-\viscFrac) \mollEps^2}{2 \DE} \sum_{k=1}^{n} \tau \int_{\OmSph} \text{Q}^\regL(\psiL{k} + \delta) \abs{ \nabla_\xflow \muPsidL{k} }^2 \dH\, \text{d}\xflow \\
	& \qquad + \frac{1-\viscFrac}{2 \DE} \sum_{k=1}^{n} \tau \int_{\OmSph} \text{Q}^\regL(\psiL{k} + \delta) \abs{ \nabla_\surfaceLtr \muPsidL{k} }^2 \dH\, \text{d}\xflow \\
	&\quad\leq \frac{\RE \DE}{2} \int_{\Omega} \abs{\uL{0}}^2 \, \text{d}\xflow + (1-\viscFrac) \int_{\OmSph} \entrop^\regL(\psiL{0} + \delta) \dH\, \text{d}\xflow \numberthis[discr_entr:u_psi_I_IV_combined] \\
	&\quad\quad + \frac{1-\viscFrac}{2} \int_{\OmSph} \psiL{0} \pot[\psiL{0}] \dH\, \text{d}\xflow + \frac{\activity^2 (1 - \viscFrac) T}{\DE} \norm{\psiL{0}}_{\LpOS{1}} \\
	& \quad\quad + C_2 T \norm{\psiL{0}}_{\LpOS{1}}^2 \lrcb{ 1 + \norm{\psiL{0}}_{\LpOS{1}} } \\
	& \quad\quad + (C_3 + C_4) T \delta \norm{\psiL{0}}_{\LpOS{1}}^2 \qquad \qquad \qquad  \text{ for all } n \in \{ 1, \dots, N \}.
	\end{align*}

    Next, we infer from \eqref{eq:nabla_potLk_bound}, \eqref{reg:FLQL_leq_FL}, \eqref{eq:mass_consv_psiLn}, \eqref{init:regularity_psi_L_0}, and applying Lemma of Fatou (cf. \cite[p. 24]{chen_liu_2013_entropy_solution}, \cite[p. 1244f]{barrett_sueli_2011_polymers_I}) that
    \begin{align}
    \label{discr_entr:limit:diffusive_x}
    \begin{split}
    & \frac{(1-\viscFrac)\mollEps^2}{2 \DE} \lim_{\delta \searrow 0} \int_{\OmSph} \text{Q}^\regL(\psiL{k} + \delta) \abs{ \nabla_\xflow \muPsidL{k} }^2 \dH\, \text{d}\xflow \\
    &\quad\geq \frac{(1-\viscFrac)\mollEps^2}{4 \DE} \lim_{\delta \searrow 0} \int_{\OmSph} \frac{\abs{\nabla_\xflow \psiL{k}}^2}{\psiL{k} + \delta} \dH\, \text{d}\xflow - C(\overline{\C}_\mollEps, \potStr, d, \Omega, \viscFrac, \mollEps, \DE^{-1}, \psi_0) \\
    &\quad\geq \frac{(1-\viscFrac)\mollEps^2}{\DE} \int_{\OmSph} \abs{\nabla_\xflow \sqrt{\psiL{k}}}^2 \dH\, \text{d}\xflow - C(\overline{\C}_\mollEps, \potStr, d, \Omega, \viscFrac, \mollEps, \DE^{-1}, \psi_0).
    \end{split}
    \end{align}
    We apply similar computations for the term $\text{Q}^\regL(\psiL{k} + \delta) \abs{\nabla_\surfaceLtr \muPsidL{k}}^2$ (cf. \cite[p. 24]{chen_liu_2013_entropy_solution}). Finally, on noting \eqref{discr_entr:limit:diffusive_x}, the regularities of the initial data \eqref{init:regularity_u_L_0}, \eqref{init:regularity_psi_L_0}, \eqref{init:psi_0_pot_0}, and \eqref{reg:FLQL_leq_FL}, applying Lemma of Fatou and the Dominated Convergence Theorem in \eqref{discr_entr:u_psi_I_IV_combined} yields the entropy estimate \eqref{est:discrete_entropy} when passing to the limit as $\delta \searrow 0$.
\end{proof}

\subsection{Regularity of the polymer number density}
\label{sec:particle_density}
The discrete-in-time entropy estimate \eqref{est:discrete_entropy} is a good starting point to establish regularity results for $\uL{n}$ and $\psiL{n}$. However, the bounds in \eqref{est:discrete_entropy} are too weak to extract strongly converging subsequences. Unfortunately, testing \eqref{eq:PLtau:psi} with $\theta := \psiL{n}$ is fraud with difficulties to improve the regularity of $\psiL{n}$.
Instead, we shall establish bounds on the discrete-in-time polymer number density
\begin{align}
\label{def:omegaL}
\omegaL{n}(\mathbf{x}) := \int_{\Sph} \psiL{n}(\mathbf{x}, \mathbf{m}) \dH(\mflow) \quad \text{ for } \mathbf{x} \in \Omega,\ n \in \{0, \dots, N\}.
\end{align}

\begin{lmm}
	\label{lmm:regularity_omegaL}
	Suppose that the assumptions \eqref{assumptions:domain_init_data} hold.
	Moreover, let the time increment satisfy
	\begin{align*}
	0 < \tau < \frac{1}{4 \widetilde{C} \norm{\psi_0}_{\LpOS{1}}^2},
	\end{align*}
	where $\widetilde{C} > 0$ depends only on $\DE^{-1}, \mollEps, \overline{\C}_\mollEps, \potStr$, and $d$. Then, the following estimate holds for $n \in \{1, \dots, N\}$
	\begin{align}
	\label{eq:regularity_omegaL}
	\begin{split}
	&\frac{1}{4} \int_{\Omega} \abs{\omegaL{n}}^2 \, \text{d}\xflow + \frac{1}{2} \sum_{k=1}^{n} \int_{\Omega} \abs{\omegaL{k} - \omegaL{k-1}}^2 \, \text{d}\xflow + \frac{\mollEps^2}{2 \DE} \sum_{k=1}^{n} \tau \int_{\Omega} \abs{ \nabla_\xflow \omegaL{k} }^2 \, \text{d}\xflow \\
	&\quad\leq \lrcb{ 1 + \widetilde{C} \norm{\psi_0}_{\LpOS{1}}^2 T \exp\lrcb{ \widetilde{C} \norm{\psi_0}_{\LpOS{1}}^2 T } } \norm{\psi_0}_{\L^2(\Omega; \L^1({\SMat^{d-1}}))}^2 \\
	&\quad< \infty.
	\end{split}
	\end{align}
\end{lmm}

\begin{proof}
	We deduce from choosing $\theta := \omegaL{n} \in \H^1(\Omega)$ in \eqref{eq:PLtau:psi}, \eqref{eq:nabla_potLk_bound}, \eqref{reg:QL_leq}, \eqref{eq:mass_consv_psiLn}, \eqref{def:omegaL},
	\begin{align*}
	\int_{\Omega} \omegaL{n} \nabla_\xflow \omegaL{n} \cdot \uL{n} \, \text{d}\xflow = \int_{\Omega} \frac{1}{2} \nabla_\xflow (\omegaL{n})^2 \cdot \uL{n} \, \text{d}\xflow = 0,
	\end{align*}
	and Young's inequality that
	\begin{align*}
	&\frac{1}{2 \tau} \int_{\Omega} \left[ \abs{\omegaL{n}}^2 + \abs{\omegaL{n} - \omegaL{n-1}}^2  - \abs{\omegaL{n-1}}^2 \right] \, \text{d}\xflow + \frac{\mollEps^2}{\DE} \int_{\Omega} \abs{ \nabla_\xflow \omegaL{n} }^2 \, \text{d}\xflow \\
	&\quad= -\frac{\mollEps^2}{\DE} \int_{\OmSph} \text{Q}^\regL(\psiL{n}) \nabla_\xflow \potL{n} \cdot \nabla_\xflow \omegaL{n} \dH\, \text{d}\xflow \\
    &\quad\leq \frac{\mollEps^2}{2 \DE} \int_{\Omega} \abs{\nabla_\xflow \omegaL{n}}^2 \, \text{d}\xflow + C\lrcb{\DE^{-1}, \mollEps, \overline{\C}_\mollEps, \potStr, d} \norm{\psi_0}_{\LpOS{1}}^2 \int_{\Omega} \abs{\omegaL{n}}^2 \, \text{d}\xflow.
	\end{align*}

	Thus, a discrete integration in time and \eqref{init:psi_0} yield
	\begin{align}
	\label{eq:omegaL_combined}
	\begin{split}
	&\frac{1}{2} \int_{\Omega} \abs{\omegaL{n}}^2 \, \text{d}\xflow + \frac{1}{2} \sum_{k=1}^{n} \int_{\Omega} \abs{\omegaL{k} - \omegaL{k-1}}^2 \, \text{d}\xflow + \frac{\mollEps^2}{2 \DE} \sum_{k=1}^{n} \tau \int_{\Omega} \abs{\nabla_\xflow \omegaL{k}}^2 \, \text{d}\xflow \\
	&\quad\leq \frac{1}{2} \norm{\psi_0}_{\L^2(\Omega; \L^1({\SMat^{d-1}}))}^2 + C_1\lrcb{\DE^{-1}, \mollEps, \overline{\C}_\mollEps, \potStr, d} \norm{\psi_0}_{\LpOS{1}}^2 \sum_{k=1}^n \tau \int_{\Omega} \abs{\omegaL{k}}^2 \, \text{d}\xflow
	\end{split}
	\end{align}
	for $n \in \{1, \dots, N\}$, where $C_1 = C_1\lrcb{\DE^{-1}, \mollEps, \overline{\C}_\mollEps, \potStr, d} > 0$ is independent of $\tau$ and $\text{L}$.
    We complete the proof by restricting the range of the time increment to
    \begin{align*}
    0 < \tau < \frac{1}{4 C_1 \norm{ \psi_0 }_{\LpOS{1}}^2},
    \end{align*}
    and by applying a discrete version of Gronwall's inequality (cf. \cite[Lemma 4.3.2]{werner1986ode}).
\end{proof}

\subsection{Time regularity}
\label{sec:time_regularity}
In this section, we derive bounds on the discrete time derivatives of $\uL{n}$ and $\psiL{n}$ in appropriate spaces to apply compactness results to be able to extract strongly converging subsequences.

\begin{lmm}[Time regularity of $\mathbf{u}$]
	\label{lmm:time_reg:uL}
	Suppose that the assumptions \eqref{assumptions:domain_init_data} hold and that $\tau > 0$ satisfies the condition of Lemma \ref{lmm:regularity_omegaL}.
	Then, we have that
	\begin{align}
	\label{eq:time_reg:uL}
	\sum_{k=1}^N \tau \norm{ \frac{\uL{k} - \uL{k-1}}{\tau} }_{ (\Vk{2})^\prime}^2 \leq C,
	\end{align}
	where $C>0$ depends solely on $\RE, \RE^{-1}, \DE, \DE^{-1}, \viscFrac, \viscFrac^{-1}, \activity, T$, $\overline{\C}_\mollEps, \mollEps, \potStr, d, \Omega, \mathbf{u}_0,$ and $\psi_0$.
\end{lmm}

\begin{proof}
	Let $\mathbf{w} \in \Vk{2} \hookrightarrow \L^\infty(\Omega)$. Then, we infer from \eqref{eq:PLtau:u} and \eqref{reg:QL_leq} that, for $n \in \{1, \dots, N\}$,
	\begin{align*}
	&\RE\DE \abs{ \int_{\Omega} \frac{\uL{n} - \uL{n-1}}{\tau} \cdot \mathbf{w} \, \text{d}\xflow } \\
	&\quad \leq \RE\DE \int_{\Omega} \abs{\uL{n-1}} \abs{\nabla_\xflow \uL{n}} \abs{\mathbf{w}} \, \text{d}\xflow + \viscFrac\DE \norm{\nabla_\xflow \uL{n}}_{\L^2(\Omega)} \norm{\nabla_\xflow \mathbf{w}}_{\L^2(\Omega)} \\
	&\qquad + \frac{(1-\viscFrac) \DE}{2} \int_{\OmSph} \text{Q}^\regL(\psiL{n}) \abs{\nabla_\xflow \uL{n} : \mflow\otimes\mflow} \abs{\nabla_\xflow \mathbf{w} : \mflow\otimes\mflow} \dH(\mflow) \, \text{d}\xflow \\
	&\qquad + (1-\viscFrac) \int_{\OmSph} \abs{\psiL{n}} \abs{\nabla_\xflow \potL{n}} \abs{\mathbf{w}} \dH\, \text{d}\xflow  \\
	&\qquad + (1-\viscFrac) C(d) \int_{\OmSph} \abs{\psiL{n}} \abs{\nabla_\surfaceLtr \potL{n}} \abs{\nabla_\xflow \mathbf{w}} \dH\, \text{d}\xflow \numberthis[eq:tr:uL:tested] \\
	&\qquad + (1-\viscFrac) \lrcb{ C(d) + \activity } \int_{\OmSph} \abs{\psiL{n}} \abs{\nabla_\xflow \mathbf{w}} \dH\, \text{d}\xflow \\
	&\quad =: I + II + III + IV + V + VI.
	\end{align*}

	First, on noting \eqref{eq:mass_consv_psiLn}, \eqref{reg:QL_leq}, \eqref{init:regularity_psi_L_0}, and \eqref{eq:regularity_omegaL}, applying H\"older's inequality yields
	\begin{align*}
	III
	&\leq C(\viscFrac, \DE, d) \lrcb{ \int_{\OmSph} \text{Q}^\regL(\psiL{n}) \lrcb{ \nabla_\xflow \uL{n} : \mflow\otimes\mflow }^2 \dH(\mflow) \, \text{d}\xflow }^\frac{1}{2} \cdot \\
	&\qquad\qquad\qquad \cdot \norm{\nabla_\xflow \mathbf{w}}_{\L^4(\Omega)} \norm{\omegaL{n}}_{\L^2(\Omega)}^{\frac{1}{2}} \\
	&\leq C \lrcb{ \int_{\OmSph} \text{Q}^\regL(\psiL{n}) \lrcb{ \nabla_\xflow \uL{n} : \mflow\otimes\mflow }^2 \dH(\mflow) \, \text{d}\xflow }^\frac{1}{2} \norm{\mathbf{w}}_{\Vk{2}}, \numberthis[eq:tr:uL:III]
	\end{align*}
	where $C = C(\viscFrac, \DE, d, \DE^{-1}, \mollEps, \overline{\C}_\mollEps, \potStr, \psi_0) >0$.
	Furthermore, we deduce from \eqref{eq:nabla_potLk_bound}, \eqref{def:omegaL}, and \eqref{eq:regularity_omegaL} that
	\begin{align}
	\label{eq:tr:uL:IV_V}
	\begin{split}
	IV + V
	&\leq C\lrcb{ \viscFrac, d, \overline{\C}_\mollEps, \potStr } \norm{\psi_0}_{\LpOS{1}} \lrcb{ \max_{k \in \{1, \dots, N\}} \norm{\omegaL{k}}_{\L^2(\Omega)} } \norm{\mathbf{w}}_{\Vk{2}} \\
	&\leq C\lrcb{\viscFrac, d, \overline{\C}_\mollEps, \potStr, \DE^{-1}, \mollEps, T, \psi_0} \norm{\mathbf{w}}_{\Vk{2}}.
	\end{split}
	\end{align}

    The other terms in \eqref{eq:tr:uL:tested} are handled easily such that combining \eqref{eq:tr:uL:tested} - \eqref{eq:tr:uL:IV_V} yields
	\begin{align}
	\label{eq:tr:uL:combined}
	\begin{split}
	&\abs{ \int_{\Omega} \frac{\uL{n} - \uL{n-1}}{\tau} \cdot \mathbf{w} \, \text{d}\xflow } \\
	&\quad\leq C  \bigg( 1 + \norm{\nabla_\xflow \uL{n}}_{\L^2(\Omega)} \\
	&\qquad\qquad + \lrcb{ \int_{\OmSph} \text{Q}^\regL(\psiL{n}) \lrcb{ \nabla_\xflow \uL{n} : \mflow\otimes\mflow }^2 \dH(\mflow)\, \text{d}\xflow }^\frac{1}{2} \bigg) \norm{\mathbf{w}}_{\Vk{2}},
	\end{split}
	\end{align}
	where $C = C\lrcb{ \RE, \RE^{-1}, \DE, \DE^{-1}, \viscFrac, \activity, T, \overline{\C}_\mollEps, \mollEps, \potStr, d, \Omega, \mathbf{u}_0, \psi_0 } > 0$.
	Hence, on noting \eqref{est:discrete_entropy}, multiplying \eqref{eq:tr:uL:combined} by $\tau$ and summing up from $n = 1 \to N$, we complete the proof.
\end{proof}

\begin{lmm}[Time regularity of $\psi$]
	\label{lmm:time_reg:psiL}
	Suppose that the assumptions \eqref{assumptions:domain_init_data} hold and that $\tau > 0$ satisfies the condition of Lemma \ref{lmm:regularity_omegaL}.
	Then, we have that
	\begin{align}
	\label{eq:time_reg:psiL}
	\sum_{k=1}^N \tau \norm{\frac{\psiL{k} - \psiL{k-1}}{\tau}}_{ \lrcb{\HkOS{4}}^\prime }^2 \leq C,
	\end{align}
	where $C>0$ depends solely on $\Omega, \RE, \RE^{-1}, \DE, \DE^{-1}, \viscFrac, \viscFrac^{-1}, \activity, T, \overline{\C}_\mollEps, \mollEps, \mollEps^{-1}, \potStr, d, \mathbf{u}_0$, and $\psi_0$.
\end{lmm}

\begin{proof}
	Let $n \in \{1, \dots, N\}$ and $\theta \in \H^4({\Omega\times\Sph}) \hookrightarrow \mathrm{W}^{1, \infty}({\Omega\times\Sph})$. Then, we infer from \eqref{eq:PLtau:psi} that
	\begin{align*}
	&\abs{ \int_{\OmSph} \frac{\psiL{n} - \psiL{n-1}}{\tau} \theta \dH\, \text{d}\xflow } \\
	&\quad\leq \int_{\OmSph} \abs{\psiL{n}} \abs{\uL{n}} \abs{\nabla_\xflow \theta} \dH\, \text{d}\xflow \\
	&\quad\quad + \frac{\mollEps^2}{\DE} \int_{\OmSph} \abs{\nabla_\xflow \psiL{n}} \abs{\nabla_\xflow \theta} \dH\, \text{d}\xflow
	+ \frac{1}{\DE} \int_{\OmSph} \abs{\nabla_\surfaceLtr \psiL{n}} \abs{\nabla_\surfaceLtr \theta} \dH\, \text{d}\xflow \\
	&\quad\quad + \frac{\mollEps^2}{\DE} \int_{\OmSph} \text{Q}^\regL(\psiL{n}) \abs{\nabla_\xflow \potL{n}} \abs{\nabla_\xflow \theta} \dH\, \text{d}\xflow \\
	&\quad\quad+ \frac{1}{\DE} \int_{\OmSph} \text{Q}^\regL(\psiL{n}) \abs{\nabla_\surfaceLtr \potL{n}} \abs{\nabla_\surfaceLtr \theta} \dH\, \text{d}\xflow \\
	&\quad\quad + C(d) \int_{\OmSph} \text{Q}^\regL(\psiL{n}) \abs{\nabla_\xflow \uL{n}} \abs{\nabla_\surfaceLtr \theta} \dH\, \text{d}\xflow \\
	&\quad=: I + II + III + IV + V + VI. \numberthis[eq:tr:psiL:tested]
	\end{align*}

	First, we obtain from \eqref{reg:QL_leq}, \eqref{eq:mass_consv_psiLn}, \eqref{init:regularity_psi_L_0}, and \eqref{eq:nabla_potLk_bound} that
	\begin{align}
	\label{eq:tr:psiL:IV_V}
	\begin{split}
	IV + V
	&\leq C(\mollEps, \DE^{-1}) \Big( \norm{\nabla_\xflow \potL{n}}_{\LpOS{\infty}} \\
	&\qquad\qquad + \norm{\nabla_\surfaceLtr \potL{n}}_{\LpOS{\infty}} \Big) \norm{\psiL{n}}_{\LpOS{1}} \norm{\theta}_{\mathrm{W}^{1, \infty}({\Omega\times\Sph})} \\
	&\leq C(\mollEps, \DE^{-1}, \Omega, d) \norm{\psi_0}_{\LpOS{1}}^2 \norm{\theta}_{\HkOS{4}}.
	\end{split}
	\end{align}
    The estimates of the other terms in \eqref{eq:tr:psiL:tested} can be found in \cite[p.25]{chen2017activeLC}.
	Thus, we conclude that
	\begin{align}
	\label{eq:tr:psiL:combined}
	\begin{split}
	&\abs{ \int_{\OmSph} \frac{\psiL{n} - \psiL{n-1}}{\tau} \theta \dH\, \text{d}\xflow } \\
	&\quad\leq C \Big( 1 + \norm{ \nabla_\xflow \sqrt{\psiL{n}} }_{\LpOS{2}} \\
	&\qquad\qquad + \norm{ \nabla_\surfaceLtr \sqrt{\psiL{n}} }_{\LpOS{2}} + \norm{\nabla_\xflow \uL{n}}_{\L^2(\Omega)} \Big) \norm{\theta}_{\HkOS{4}},
	\end{split}
	\end{align}
	where $C = C\lrcb{ \Omega, \RE, \RE^{-1}, \DE, \DE^{-1}, \viscFrac, \activity, T, \overline{\C}_\mollEps, \mollEps, \potStr, d, \mathbf{u}_0, \psi_0 } > 0$.
    Now, to complete the proof, we multiply \eqref{eq:tr:psiL:combined} by $\tau$, sum up from $n = 1 \to N$, and apply \eqref{est:discrete_entropy}.
\end{proof}

\section{Passage to the limit: Existence of weak solutions}
\label{sec:passage_limit}
In this section, our goal is to pass to the limit in \eqref{def:PLtau} as $\text{L} \nearrow \infty$ and $\tau \searrow 0$.
Before doing so however we reformulate \eqref{def:PLtau} continuously in time.
For this task, let us define the following continuous-in-time functions:
\begin{equation}
\label{def:time_cont_fcts}
\begin{alignedat}{2}
\uL{{\tau,+}}(\cdot, t) &:= \uL{n}(\cdot), \qquad \uL{{\tau,-}}(\cdot, t) := \uL{n-1}(\cdot),\quad && t \in (t_{n-1}, t_n] , \quad n \in \{1, \dots, N\}, \\
\uL{{\tau,-}}(\cdot, 0) &:= \uL{0}(\cdot) \in \Hflow_{0,\text{div}}^1(\Omega).
\end{alignedat}
\end{equation}

Moreover, we define similarly to Chen and Liu (cf. \cite[p. 26]{chen_liu_2013_entropy_solution}) the difference quotient with respect to time by
\begin{align}
\label{def:time_deriv}
\partial_{t}^{\tau} \uL{{\tau,+}}(\cdot, t) := \frac{\uL{n}(\cdot) - \uL{n-1}(\cdot)}{\tau}, \quad t \in (t_{n-1}, t_n], \quad n \in \{1, \dots, N\}.
\end{align}
We introduce the same notation for $\psiL{\tau, \pm}$ and $\partial_{t}^{\tau} \psiL{{\tau,+}}$
with $\psiL{{\tau,-}}(\cdot, \tilde{\cdot}, 0) = \psiL{0}(\cdot, \tilde{\cdot})$.
From now on, let
\begin{align}
\label{tc:cond:tau}
\tau \in o(\text{L}^{-1}) \text{ satisfy the condition in Lemma } \ref{lmm:regularity_omegaL}.
\end{align}

On noting the time continuous notation \eqref{def:time_cont_fcts}, we can rewrite the weak formulation \eqref{def:PLtau} after a discrete integration with respect to time as follows:\\
Find $\uL{\tau, \pm}(\cdot, t) \in \Hflow_{0,\text{div}}^1(\Omega)$, $t \in (0,T]$, such that
\begin{align*}
&\RE \DE \int_0^T\int_{\Omega} \partial_{t}^{\tau} \uL{{\tau,+}} \cdot \mathbf{w} \, \text{d}\xflow\,\text{dt} + \RE \DE \int_0^T\int_{\Omega} \lrcb{ (\uL{{\tau,-}} \cdot \nabla_\xflow) \uL{{\tau,+}} } \cdot \mathbf{w} \, \text{d}\xflow\,\text{dt} \\
&\quad\quad + \viscFrac \DE \int_0^T\int_{\Omega} \nabla_\xflow \uL{{\tau,+}} : \nabla_\xflow \mathbf{w} \, \text{d}\xflow\,\text{dt} \\
&\quad\quad + \frac{(1-\viscFrac) \DE}{2} \int_0^T\int_{\OmSph} \text{Q}^\regL(\psiL{{\tau,+}}) \lrcb{ \nabla_\xflow \uL{{\tau,+}} : \mflow\otimes\mflow } \cdot \\
&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad \cdot \lrcb{\nabla_\xflow \mathbf{w} : \mflow\otimes\mflow} \dH(\mflow)\, \text{d}\xflow\,\text{dt} \\
&\quad= - \frac{(1-\viscFrac)}{2} \int_0^T\int_{\OmSph} \psiL{{\tau,+}} \nabla_\xflow \pot[\psiL{{\tau,+}} + \psiL{{\tau,-}}] \cdot \mathbf{w} \dH\, \text{d}\xflow\,\text{dt} \numberthis[tc:eq:u] \\
&\quad\quad - \frac{1-\viscFrac}{2} \int_0^T\int_{\OmSph} \text{Q}^\regL(\psiL{{\tau,+}}) \lrcb{ (\mathbb{I} - \mflow\otimes\mflow) \nabla_\xflow \mathbf{w} \mathbf{m} } \cdot \\
&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad \cdot \nabla_\surfaceLtr \pot[\psiL{{\tau,+}} + \psiL{{\tau,-}}] \dH(\mflow)\, \text{d}\xflow\,\text{dt} \\
&\quad\quad - (1-\viscFrac) \int_0^T\int_{\OmSph} \psiL{{\tau,+}} \lrcb{d\mflow\otimes\mflow - \mathbb{I}} : \nabla_\xflow \mathbf{w} \dH(\mflow)\, \text{d}\xflow\,\text{dt} \\
&\quad\quad - \activity (1 - \viscFrac) \int_0^T\int_{\OmSph} \text{Q}^\regL(\psiL{{\tau,+}}) \lrcb{ \mflow\otimes\mflow : \nabla_\xflow \mathbf{w} } \dH(\mflow)\, \text{d}\xflow\,\text{dt}
\end{align*}
for all $\mathbf{w} \in \L^1(0,T; \Hflow_{0,\text{div}}^1(\Omega))$.
Moreover, find $\psiL{\tau, \pm}(\cdot, \tilde{\cdot}, t) \in \HkOS{1} \cap \Zp{2}$, $t \in (0,T]$, such that
\begin{align*}
&\int_0^T\int_{\OmSph} \partial_{t}^{\tau} \psiL{{\tau,+}} \theta \dH\, \text{d}\xflow\,\text{dt} - \int_0^T\int_{\OmSph} \psiL{{\tau,+}} \uL{{\tau,+}} \cdot \nabla_\xflow \theta \dH\, \text{d}\xflow\,\text{dt} \\
&\qquad + \frac{1}{\DE} \int_0^T\int_{\OmSph} \left[ \mollEps^2 \nabla_\xflow \psiL{{\tau,+}} \cdot \nabla_\xflow \theta + \nabla_\surfaceLtr \psiL{{\tau,+}} \cdot \nabla_\surfaceLtr \theta \right] \dH\, \text{d}\xflow\,\text{dt} \\
&\qquad + \frac{\mollEps^2}{2 \DE} \int_0^T\int_{\OmSph} \text{Q}^\regL(\psiL{{\tau,+}}) \nabla_\xflow \pot[\psiL{{\tau,+}} + \psiL{{\tau,-}}] \cdot \nabla_\xflow \theta \dH\, \text{d}\xflow\,\text{dt} \numberthis[tc:eq:psi] \\
&\qquad + \frac{1}{2 \DE} \int_0^T\int_{\OmSph} \text{Q}^\regL(\psiL{{\tau,+}}) \nabla_\surfaceLtr \pot[\psiL{{\tau,+}} + \psiL{{\tau,-}}] \cdot \nabla_\surfaceLtr \theta \dH\, \text{d}\xflow\,\text{dt} \\
&\qquad - \int_0^T\int_{\OmSph} \text{Q}^\regL(\psiL{{\tau,+}}) \lrcb{ (\mathbb{I} - \mflow\otimes\mflow) \nabla_\xflow \uL{{\tau,+}} \mathbf{m} } \cdot \nabla_\surfaceLtr \theta \dH(\mflow)\, \text{d}\xflow\,\text{dt} = 0
\end{align*}
for all $\theta \in \L^1(0,T; \HkOS{1})$.
We note that the problem \eqref{tc:eq:u}-\eqref{tc:eq:psi} is equivalent to \eqref{eq:PLtau:u}-\eqref{eq:PLtau:psi}, to which existence of solutions has been proved (cf. Lemma \ref{lmm:existence_fp}).
Furthermore, let us define for future purposes
\begin{subequations}
	\begin{align}
	\omegaL{\tau, \pm}(\mathbf{x}, t) &:= \int_{\Sph} \psiL{\tau, \pm}(\mathbf{x}, \mathbf{m}, t) \dH(\mflow) \in \mathbb{R}_{\geq 0}, \label{tc:def:omegaL} \\
	\StressL{\tau, \pm}(\mathbf{x}, t) &:= \int_{\Sph} \psiL{\tau, \pm}(\mathbf{x}, \mathbf{m}, t) \mflow\otimes\mflow \dH(\mflow) \in \mathbb{R}^{d \times d} \label{tc:def:StressL}
	\end{align}
\end{subequations}
for a.e. $(\mathbf{x}, t) \in \Omega \times (0, T]$. Thus, we can rewrite the interaction potential as follows:
\begin{align}
\label{tc:pot_repres}
\pot[\psiL{\tau, \pm}](\mathbf{x}, \mathbf{m}, t) = \potStr \Je{\omegaL{\tau, \pm}(\cdot, t)}(\mathbf{x}) - \potStr \JeMat{\StressL{\tau, \pm}(\cdot, t)}(\mathbf{x}) : \mflow\otimes\mflow
\end{align}
for $(\mathbf{x}, \mathbf{m}, t) \in {\Omega\times\Sph} \times (0, T]$.
Moreover, we obtain from \eqref{init:regularity_u_L_0}, $\tau \in o(\text{L}^{-1})$, \eqref{def:time_cont_fcts}, \eqref{est:discrete_entropy}, and \eqref{eq:regularity_omegaL} that
\begin{align*}
&\esssup_{t \in [0,T]} \int_{\Omega} \abs{\uL{\tau, \pm}(\mathbf{x}, t)}^2 \, \text{d}\xflow  +  \frac{1}{\tau} \int_0^T\int_{\Omega} \abs{ \uL{{\tau,+}} - \uL{{\tau,-}} }^2 \, \text{d}\xflow\,\text{dt}  \\
&\quad\quad+  \int_0^T\int_{\Omega} \abs{\nabla_\xflow \uL{\tau, \pm}}^2 \, \text{d}\xflow\,\text{dt} \\
&\quad\quad + \esssup_{t \in [0,T]} \int_{\OmSph} \mathcal{F}(\psiL{{\tau,+}}(\mathbf{x}, \mathbf{m}, t)) \dH(\mflow)\, \text{d}\xflow \\
&\quad\quad + \frac{1}{\tau \text{L}} \int_0^T\int_{\OmSph} \abs{ \psiL{{\tau,+}} - \psiL{{\tau,-}} }^2 \dH\, \text{d}\xflow\,\text{dt} \numberthis[tc:regularity_past] \\
&\quad\quad + \int_0^T\int_{\OmSph} \Big[ | \nabla_\xflow \sqrt{\psiL{{\tau,+}}} |^2 + | \nabla_\surfaceLtr \sqrt{\psiL{{\tau,+}}} |^2 \Big] \dH\, \text{d}\xflow\,\text{dt} \\
&\quad\quad + \esssup_{t \in [0,T]} \int_{\Omega} \abs{\omegaL{{\tau,+}}(\mathbf{x}, t)}^2 \, \text{d}\xflow + \int_0^T\int_{\Omega} \abs{\nabla_\xflow \omegaL{{\tau,+}}}^2 \, \text{d}\xflow\,\text{dt} \\
&\quad\quad + \frac{1}{\tau} \int_0^T\int_{\Omega} \abs{ \omegaL{{\tau,+}} - \omegaL{{\tau,-}} }^2 \, \text{d}\xflow\,\text{dt} \\
&\quad\leq C_\star < \infty,
\end{align*}
where $C_\star = C_\star \lrcb{ \RE, \RE^{-1}, \DE, \DE^{-1}, \viscFrac, \viscFrac^{-1}, \activity, T, \overline{\C}_\mollEps, \mollEps, \mollEps^{-1}, \potStr, d, \Omega, \mathbf{u}_0, \psi_0 } > 0$ is independent of $\tau$ and $\text{L}$.

\begin{thrm}[Main Theorem: Existence of weak solutions]
	\label{thrm:main}
	Suppose that the assumptions \eqref{assumptions:domain_init_data} and the condition \eqref{tc:cond:tau}, which relates $\tau$ to $\text{L}$, hold.
	Then, there exists a subsequence of $\{ (\uL{\tau, \pm}, \psiL{\tau, \pm}) \}_{\text{L} > 1}$ (not relabelled) and a pair of functions $(\mathbf{u}, \psi)$ such that
	\begin{subequations}
		\begin{gather}
		\mathbf{u} \in \L^\infty(0, T; \mathbf{L}^2(\Omega)) \cap \L^2(0, T; \Hflow_{0,\text{div}}^1(\Omega)), \label{tc:limit:u} \\
		\partial_t \mathbf{u} \in \L^2\lrcb{ 0,T; ( \Vk{2} )^\prime }, \label{tc:limit:dt_u}
		\end{gather}
	\end{subequations}
	and
	\begin{subequations}
		\begin{gather}
		\psi \in \L^\infty(0,T; \L^2(\Omega; \L^1({\SMat^{d-1}}))), \quad \sqrt{\psi} \in \L^2(0,T; \HkOS{1}), \label{tc:limit:psi} \\
		\partial_t \psi \in \L^2( 0,T; (\HkOS{4})^\prime ), \label{tc:limit:dt_psi}
		\end{gather}
	\end{subequations}
	with $\psi \geq 0$ a.e. in ${\Omega\times\Sph} \times [0,T]$ and
	\begin{align}
	\label{tc:limit:finite_entropy}
	\mathcal{F}(\psi) \in \L^\infty(0, T; \LpOS{1}).
	\end{align}
	Moreover, it holds that, as $\text{L} \nearrow \infty$ (and thereby $\tau \searrow 0$),
	\begin{subequations}
		\begin{alignat}{3}
		\uL{\tau, \pm} &\overset{*}{\rightharpoonup} \mathbf{u} \quad &&\text{ weak}^* &&\text{ in } \L^\infty(0, T; \mathbf{L}^2(\Omega)) \label{tc:limit:u:weakstar} \\
		\uL{\tau, \pm} &\rightharpoonup \mathbf{u} \quad &&\text{ weakly} &&\text{ in } \L^2(0,T; \Hflow_{0,\text{div}}^1(\Omega)), \label{tc:limit:u:weak} \\
		\uL{\tau, \pm} &\to \mathbf{u} \quad &&\text{ strongly} &&\text{ in } \L^2(0,T; \mathbf{L}^r(\Omega)), \label{tc:limit:u:strong}\\
		\partial_{t}^{\tau} \uL{{\tau,+}} &\rightharpoonup \partial_t \mathbf{u} \quad &&\text{ weakly} &&\text{ in } \L^2( 0,T; (\Vk{2})^\prime ), \label{tc:limit:u:dt}
		\end{alignat}
	\end{subequations}
	where $r \in [1, \infty)$ if $d = 2$ and $r \in [1, 6)$ if $d = 3$, and
	\begin{subequations}
		\begin{alignat}{3}
		\nabla_\xflow \sqrt{\psiL{{\tau,+}}} &\rightharpoonup \nabla_\xflow \sqrt{\psi} \quad &&\text{ weakly} &&\text{ in } \L^2(0,T; \mathbf{L}^2({\Omega\times\Sph})), \label{tc:limit:psi:nablax_sqrt} \\
		\nabla_\surfaceLtr \sqrt{\psiL{{\tau,+}}} &\rightharpoonup \nabla_\surfaceLtr \sqrt{\psi} \quad &&\text{ weakly} &&\text{ in } \L^2(0,T; \mathbf{L}^2({\Omega\times\Sph})), \label{tc:limit:psi:nablag_sqrt} \\
		\sqrt{\psiL{{\tau,+}}} &\to \sqrt{\psi} \quad&&\text{ strongly} &&\text{ in } \L^4(0,T; \L^4(\Omega; \L^2({\SMat^{d-1}}))), \label{tc:limit:psi:sqrt_strong} \\
		\psiL{\tau, \pm} &\to \psi \quad&&\text{ strongly} &&\text{ in } \L^2(0,T; \L^2(\Omega; \L^1({\SMat^{d-1}}))), \label{tc:limit:psi:strong} \\
		\text{Q}^\regL(\psiL{{\tau,+}}) &\to \psi \quad&&\text{ strongly} &&\text{ in } \L^2(0,T; \L^2(\Omega; \L^1({\SMat^{d-1}}))), \label{tc:limit:psi:QL_strong} \\
		\partial_{t}^{\tau} \psiL{{\tau,+}} &\rightharpoonup \partial_t \psi \quad&&\text{ weakly} &&\text{ in } \L^2( 0,T; (\HkOS{4})^\prime ). \label{tc:limit:psi:dt}
		\end{alignat}
	\end{subequations}
	Furthermore, we introduce
	\begin{subequations}
		\label{tc:def:omega_Stress}
		\begin{align}
		\omega(\mathbf{x}, t) &:= \int_{\Sph} \psi(\mathbf{x}, \mathbf{m}, t) \dH(\mflow) \in \mathbb{R}_{\geq 0} \ \text{ and } \label{tc:def:omega}\\
		\Stress(\mathbf{x}, t) &:= \int_{\Sph} \psi(\mathbf{x}, \mathbf{m}, t) \mflow\otimes\mflow \dH(\mflow) \in \mathbb{R}^{d \times d} \label{tc:def:Stress}
		\end{align}
	\end{subequations}
	for a.e. $(\mathbf{x}, t) \in \Omega \times (0,T]$. It holds that
	\begin{subequations}
		\begin{align}
		\omega &\in \L^\infty(0, T; \L^2(\Omega)) \cap \L^2(0, T; \H^1(\Omega)), \label{tc:limit:omega} \\
		\Stress &\in \L^\infty(0, T; \mathbb{L}^2(\Omega)), \label{tc:limit:Stress}
		\end{align}
	\end{subequations}
	and that
	\begin{subequations}
		\begin{alignat}{3}
		\omegaL{\tau, \pm} &\to \omega \quad&&\text{ strongly} &&\text{ in } \L^2(0,T; \L^2(\Omega)), \label{tc:limit:omega:strong} \\
		\StressL{\tau, \pm} &\to \Stress \quad&&\text{ strongly} &&\text{ in } \L^2(0,T; \mathbb{L}^2(\Omega)), \label{tc:limit:Stress:strong} \\
		\pot[\psiL{\tau, \pm}] &\to \pot[\psi] \quad&&\text{ strongly} &&\text{ in } \L^2(0, T; \mathrm{W}^{1, \infty}({\Omega\times\Sph})). \label{tc:limit:pot:strong}
		\end{alignat}
	\end{subequations}

	Moreover, the tuple $(\mathbf{u}, \psi)$ is a global weak solution of $(\mathcal{P})$, in the sense that
	\begin{align}
	\label{P:weak:u}
	\begin{split}
	&-\RE\DE \int_0^T\int_{\Omega} \mathbf{u} \cdot \partial_t \mathbf{w} \, \text{d}\xflow\,\text{dt}  +  \RE\DE \int_0^T\int_{\Omega} \lrcb{ (\mathbf{u} \cdot \nabla_\xflow) \mathbf{u} } \cdot \mathbf{w} \, \text{d}\xflow\,\text{dt} \\
	&\qquad + \viscFrac \DE \int_0^T\int_{\Omega} \nabla_\xflow \mathbf{u} : \nabla_\xflow \mathbf{w} \, \text{d}\xflow\,\text{dt} \\
	&\qquad + \frac{(1-\viscFrac) \DE}{2} \int_0^T\int_{\OmSph} \psi (\nabla_\xflow \mathbf{u} : \mflow\otimes\mflow) (\nabla_\xflow \mathbf{w} : \mflow\otimes\mflow) \dH(\mflow)\, \text{d}\xflow\,\text{dt} \\
	&\quad = \RE\DE \int_{\Omega} \mathbf{u}_0(\mathbf{x}) \cdot \mathbf{w}(\mathbf{x}, 0) \, \text{d}\xflow  -  (1 - \viscFrac) \int_0^T\int_{\OmSph} \psi \nabla_\xflow \pot[\psi] \cdot \mathbf{w} \dH\, \text{d}\xflow\,\text{dt} \\
	&\quad\qquad - (1 - \viscFrac) \int_0^T\int_{\OmSph} \psi \lrcb{ (\mathbb{I} - \mflow\otimes\mflow) \nabla_\xflow \mathbf{w} \mathbf{m} } \cdot \nabla_\surfaceLtr \pot[\psi] \dH(\mflow)\, \text{d}\xflow\,\text{dt} \\
	&\quad\qquad - (1 - \viscFrac) \int_0^T\int_{\OmSph} \psi (d \mflow\otimes\mflow - \mathbb{I}) : \nabla_\xflow \mathbf{w} \dH(\mflow)\, \text{d}\xflow\,\text{dt} \\
	&\quad\qquad - \alpha (1 - \viscFrac) \int_0^T\int_{\OmSph} \psi (\mflow\otimes\mflow : \nabla_\xflow \mathbf{w}) \dH(\mflow)\, \text{d}\xflow\,\text{dt} \\
	&\qquad\qquad\qquad \forall \mathbf{w} \in \mathrm{W}^{1,1}(0,T; \Vk{2}) \text{ such that } \mathbf{w}(\cdot, T) = 0,
	\end{split}
	\end{align}
	and
	\begin{align*}
	&-\int_0^T\int_{\OmSph} \psi \partial_t \theta \dH\, \text{d}\xflow\,\text{dt}  -  \int_0^T\int_{\OmSph} \psi \mathbf{u} \cdot \nabla_\xflow \theta \dH\, \text{d}\xflow\,\text{dt} \\
	&\qquad + \frac{1}{\DE} \int_0^T\int_{\OmSph} \left[ \mollEps^2 \nabla_\xflow \psi \cdot \nabla_\xflow \theta + \nabla_\surfaceLtr \psi \cdot \nabla_\surfaceLtr \theta \right] \dH\, \text{d}\xflow\,\text{dt} \\
	&\qquad + \frac{1}{\DE} \int_0^T\int_{\OmSph} \psi \left[ \mollEps^2 \nabla_\xflow \pot[\psi] \cdot \nabla_\xflow \theta + \nabla_\surfaceLtr \pot[\psi] \cdot \nabla_\surfaceLtr \theta \right] \dH\, \text{d}\xflow\,\text{dt} \\
	&\qquad - \int_0^T\int_{\OmSph} \psi \lrcb{ (\mathbb{I} - \mflow\otimes\mflow) \nabla_\xflow \mathbf{u} \mathbf{m} } \cdot \nabla_\surfaceLtr \theta \dH(\mflow)\, \text{d}\xflow\,\text{dt} \\
	&\quad = \int_{\OmSph} \psi_0(\mathbf{x}, \mathbf{m}) \theta(\mathbf{x}, \mathbf{m}, 0) \dH(\mflow)\, \text{d}\xflow \numberthis[P:weak:psi] \\
	&\qquad\qquad\qquad \forall \theta \in \mathrm{W}^{1,1}(0,T; \HkOS{4}) \text{ such that } \theta(\cdot, \tilde{\cdot}, T) = 0.
	\end{align*}

	In particular, the following energy inequality is satisfied for a.e. $s \in [0,T]$:
	\begin{align*}
	&\frac{\RE\DE}{2} \int_{\Omega} \abs{ \mathbf{u}(\mathbf{x}, s) }^2 \, \text{d}\xflow  +  \frac{\viscFrac\DE}{2} \int_0^s \abs{ \nabla_\xflow \mathbf{u} }^2 \, \text{d}\xflow\,\text{dt} \\

	&\qquad + (1 - \viscFrac) \int_{\OmSph} \mathcal{F}(\psi(\mathbf{x}, \mathbf{m}, s)) \dH(\mflow)\, \text{d}\xflow \\
	&\qquad +  \frac{1-\viscFrac}{2} \int_{\OmSph} \psi(\mathbf{x}, \mathbf{m}, s) \pot[\psi(\mathbf{x}, \mathbf{m}, s)] \dH(\mflow)\, \text{d}\xflow \\
	&\qquad + \frac{1 - \viscFrac}{\DE} \int_0^s \int_{\OmSph} \left[ \mollEps^2 \abs{ \nabla_\xflow \sqrt{\psi} }^2 + \abs{\nabla_\surfaceLtr \sqrt{\psi}}^2 \right] \dH\, \text{d}\xflow\,\text{dt} \\
	&\quad \leq \frac{\RE\DE}{2} \int_{\Omega} \abs{\mathbf{u}_0}^2 \, \text{d}\xflow + (1 - \viscFrac) \int_{\OmSph} \mathcal{F}(\psi_0) \dH\, \text{d}\xflow \numberthis[eq:main:energy_ineq] \\
	&\quad\quad + \frac{1-\viscFrac}{2} \int_{\OmSph} \psi_0 \pot[\psi_0] \dH\, \text{d}\xflow + C_{\star\star},
	\end{align*}
	where $C_{\star\star} > 0$ depends solely on $\psi_0, \alpha, \RE, \DE, \DE^{-1}, \viscFrac, \overline{\C}_\mollEps, \potStr, \Omega, d,$ and  $T$.

\end{thrm}

\begin{proof}
    For ease of notation, we note that all generic constants $0 < C < \infty$ are independent of $\text{L}$ and $\tau$ in this proof.\\
	\textit{Convergence properties of $\mathbf{u}$.}
	First, we prove the convergence properties of the velocity field. We infer from the first bound in \eqref{tc:regularity_past} that there exist subsequences $\{ \uL{\tau, \pm} \}_{\text{L} > 1}$ (not relabelled) and functions $\mathbf{u}^{\pm} \in \L^\infty(0,T; \mathbf{L}^2(\Omega))$ such that
	\begin{align}
	\label{tc:limit:proof:u:weakstar}
	\uL{\tau, \pm} \overset{*}{\rightharpoonup} \mathbf{u}^{\pm} \quad \text{ weak}^* \text{ in } \L^\infty(0,T; \mathbf{L}^2(\Omega)).
	\end{align}
	Moreover, on noting \eqref{tc:limit:proof:u:weakstar} and $\tau \in o(\text{L}^{-1})$, the second bound in \eqref{tc:regularity_past} yields
	\begin{align*}
	\int_0^T\int_{\Omega} \lrcb{ \mathbf{u}^+ - \mathbf{u}^- } \cdot \mathbf{w} \, \text{d}\xflow\,\text{dt}
	\gets \int_0^T\int_{\Omega} \lrcb{ \uL{{\tau,+}} - \uL{{\tau,-}} } \cdot \mathbf{w} \, \text{d}\xflow\,\text{dt}
	\to 0
	\end{align*}
	for all $\mathbf{w} \in \L^2(0,T; \mathbf{L}^2(\Omega)) \subset \L^1(0,T; \mathbf{L}^2(\Omega))$. Hence, we have that $\mathbf{u}^+ = \mathbf{u}^-$ and we denote this limit by $\mathbf{u} := \mathbf{u}^+ = \mathbf{u}^- \in \L^\infty(0,T; \mathbf{L}^2(\Omega))$, which yields \eqref{tc:limit:u:weakstar}.
	Furthermore, we obtain from the uniqueness of limits of sequences in the weak topology of $\L^2(0,T; \Hflow_{0,\text{div}}^1(\Omega))$, the third bound in \eqref{tc:regularity_past}, and a further extraction of a subsequence that $\mathbf{u} \in \L^2(0,T; \Hflow_{0,\text{div}}^1(\Omega))$ and that \eqref{tc:limit:u} and \eqref{tc:limit:u:weak} hold.
	On noting \eqref{eq:time_reg:uL}, we have that
	\begin{align}
	\label{tc:limit:proof:u:timereg}
	\begin{split}
	\int_0^{T - \tau} \norm{ \uL{{\tau,+}}(\cdot, t + \tau) - \uL{{\tau,+}}(\cdot, t) }_{( \Vk{2} )^\prime}^2 dt = \sum_{k=1}^{N-1} \tau \norm{ \uL{k} - \uL{k-1} }_{( \Vk{2} )^\prime}^2
	\leq C \tau^2.
	\end{split}
	\end{align}
	On noting \eqref{tc:limit:proof:u:timereg}, $\norm{\uL{{\tau,+}}}_{\L^2(0,T; \Hflow_{0,\text{div}}^1(\Omega))} \leq C$, the compact embedding
	\begin{align*}
	\Hflow_{0,\text{div}}^1(\Omega) \hookrightarrow\hookrightarrow \mathbf{L}^r(\Omega)\cap\Ldivp{2}
	\end{align*}
	with $r \in [2, \infty)$ if $d=2$ and $r \in [2, 6)$ if $d=3$, and the continuous embedding
	\begin{align}
	\label{tc:limit:proof:cont_emb:u}
	\mathbf{L}^r(\Omega) \cap \Ldivp{2} \hookrightarrow (\Vk{2})^\prime,
	\end{align}
	it follows by Theorem 1 of Dreher and J\"ungel \cite{dreher2012compact} (for a further subsequence) that $\uL{{\tau,+}} \to \mathbf{u}$ strongly in $\L^2(0,T; \mathbf{L}^r(\Omega))$.
	As $\{ \uL{{\tau,+}} \}_{\text{L} > 1}$ is uniformly bounded in $\L^2(0,T;$ $\Hflow_{0,\text{div}}^1(\Omega))$, we deduce from the second bound in \eqref{tc:regularity_past} and the continuous embedding \eqref{tc:limit:proof:cont_emb:u} that $\uL{{\tau,-}} \to \mathbf{u}$ strongly in $\L^2(0,T; \mathbf{L}^r(\Omega))$, which yields \eqref{tc:limit:u:strong}.
    The convergence result \eqref{tc:limit:u:dt} with respect to the time derivative follows similarly to \cite[p. 29ff]{chen_liu_2013_entropy_solution}.\\
	\textit{Convergence properties of $\psi$.}
	We obtain from \eqref{eq:time_reg:psiL} that
	\begin{align}
	\label{tc:limit:proof:psi:timereg}
	\begin{split}
	&\int_0^{T - \tau} \norm{ \psiL{{\tau,+}}(\cdot, \tilde{\cdot}, t + \tau) - \psiL{{\tau,+}}(\cdot, \tilde{\cdot}, t) }_{(\HkOS{4})^\prime}^2 \,\text{dt} \\
	&\qquad\qquad\qquad\quad = \sum_{k=1}^{N-1} \tau \norm{ \psiL{k} - \psiL{k-1} }_{(\HkOS{4})^\prime}^2
	\leq C \tau^2.
	\end{split}
	\end{align}
    Now, on noting \eqref{tc:limit:proof:psi:timereg}, \eqref{tc:regularity_past}, and the embeddings
    \begin{align*}
    \mathrm{W}^{1,1}({\Omega\times\Sph}) \hookrightarrow\hookrightarrow \L^s({\Omega\times\Sph}) \hookrightarrow \lrcb{\HkOS{4}}^\prime \quad \forall 1 < s < \frac{2d - 1}{2d -2},
    \end{align*}
    following the lines pp. 27-30 in \cite{chen_liu_2013_entropy_solution} yields the convergence results \eqref{tc:limit:psi:nablax_sqrt} - \eqref{tc:limit:psi:dt} with respect to $\{ \psiL{{\tau,+}} \}_{\text{L}> 1}$.
    Furthermore, we deduce from the fifth bound in \eqref{tc:regularity_past} and $\tau \in o(\text{L}^{-1})$ that \eqref{tc:limit:psi:strong} holds for $\{ \psiL{{\tau,-}} \}_{\text{L} > 1}$ as well. The nonnegativity of $\psi$ follows from the nonnegativity of $\psiL{{\tau,+}}$.
	Moreover, on noting \eqref{tc:def:omegaL}, the seventh bound in \eqref{tc:regularity_past}, and \eqref{tc:limit:psi:strong}, we have that $\psi$ $\in$ $\L^\infty(0, T; \L^2(\Omega; \L^1({\SMat^{d-1}})))$.
	This completes the proof of \eqref{tc:limit:psi}.
	Moreover, on noting \eqref{tc:limit:psi:strong}, it follows that $\eqref{tc:limit:finite_entropy}$ holds by subtracting one more subsequence.\\
	\textit{Convergence properties of $\omega$, $\Stress$, and $\pot[\psi]$.}
    On noting \eqref{tc:regularity_past}, \eqref{tc:def:omegaL}, \eqref{tc:def:omega}, \eqref{tc:def:StressL}, \eqref{tc:def:Stress}, \eqref{tc:pot_repres}, \eqref{mollifier:W_2_infty}, we infer from \eqref{tc:limit:psi:strong} that \eqref{tc:limit:omega} - \eqref{tc:limit:pot:strong} hold.\\
	\textit{Existence of weak solutions.}
    By applying the established convergence results in Theorem \ref{thrm:main}, it follows immediately that $\mathbf{u}$ and $\psi$ satisfy the weak formulations \eqref{P:weak:u} and \eqref{P:weak:psi}.\\
	\textit{Energy inequality.}
    On noting the convergence properties \eqref{tc:limit:u:weakstar}, \eqref{tc:limit:u:weak}, \eqref{tc:limit:u:strong}, \eqref{tc:limit:psi:nablax_sqrt}, \eqref{tc:limit:psi:nablag_sqrt}, and \eqref{tc:limit:psi:strong}, we infer from the (weak) lower-semicontinuity of the terms on the left-hand side of \eqref{est:discrete_entropy} that the energy inequality \eqref{eq:main:energy_ineq} holds by subtracting one more subsequence.
\end{proof}

\textbf{Acknowledgments.} This work has been supported by the Research Training Group (Graduiertenkolleg) 2339 “Interfaces, Complex Structures,
and Singular Limits in Continuum Mechanics” of the German Science Foundation (DFG). The author is grateful to
G\"unther Gr\"un and Stefan Metzger for many valuable discussions.

\bibliographystyle{acm}