Dynamics and steady state of squirmer motion in liquid crystal

Abstract

We analyse a nonlinear partial differential equation system describing the motion of a microswimmer in a nematic liquid crystal environment. For the microswimmer’s motility, the squirmer model is used in which self-propulsion enters the model through the slip velocity on the microswimmer’s surface. The liquid crystal is described using the well-established Beris–Edwards formulation. In previous computational studies, it was shown that the squirmer, regardless of its initial configuration, eventually orients itself either parallel or perpendicular to the preferred orientation dictated by the liquid crystal. Furthermore, the corresponding solution of the coupled nonlinear system converges to a steady state. In this work, we rigorously establish the existence of steady state and also the finite-time existence for the time-dependent problem in a periodic domain. Finally, we will use a two-scale asymptotic expansion to derive a homogenised model for the collective swimming of squirmers as they reach their steady-state orientation and speed.

1. Introduction

Microswimmers are objects of micron size which are immersed in a fluid and capable of autonomous motion. They are ubiquitous in nature, as exemplified by bacteria and eukaryotic cells. Recently, synthetic microswimmers for applications in medicine and material repair have been introduced in [1], see also reviews [2, 3]. Transport of microswimmers, both living and synthetic, as well as effective properties of suspensions populated by many such microswimmers, largely depends on how they respond to surrounding environment. Modelling microswimmers has become a growing area of research. The case when microswimmers are immersed in a Newtonian fluid has been intensively studied – see [4–10] and reviews [11–14]. However, bacteria often swim in biofluids which demonstrate viscoelastic or anisotropic properties very different from those of isotropic Newtonian fluids. For example, Helicobacter pylori bacteria are present in stomach and are associated with diseases such as chronic atrophic gastritis and ulcer [15, 16]. The ‘success’ in the inflammation of stomach walls by H. Pylori depends on how the bacterium reorients itself in the mucous protective layer. Note that mucus is a viscoelastic fluid which exhibits properties of a liquid crystal for a certain range of macroscopic parameters [17, 18]. In addition to medical relevance, experimental realisation which combines bacteria with a nematic water-based and non-toxic (to bacteria) liquid crystal led to a wealth of intriguing observations such as collective phenomena for small bacterial concentrations, moving topological defects and visualisation of flagella beating [19–24].

Modelling self-propulsion of microswimmers in both Newtonian and non-Newtonian fluids has received much attention in applied mathematics. These modelling approaches range from minimal one when microswimmers are modelled as point dipoles, or rigid bodies with point motors, or rigid dumbbells, e.g., [6, 7, 25–30], to more detailed ones with, for example, explicit description of self-propelling with flagella or cilia in both viscous and viscoelastic fluids, see, e.g., reviews [31–33] and the book [34].

One of the most well-established model of microswimmer is the so-called squirmer. The model was initially introduced in [35] for Paramecium, a micro-organism which swims with the help of small elastic appendages called cilia. The main modelling assumption for squirmers is that the body is non-deformable and the swimming effect is introduced via a given slip velocity profile on the body surface that models the cilia’s activity. Analysis of squirmers immersed in a Newtonian fluid, from the well-posedness to the relation between the slip profile and the resulting velocity, has been the focus of many authors [36–42].

To describe a nematic liquid crystal, we use the well-established Beris–Edwards model [43], a highly nonlinear partial differential equation (PDE) model coupling Navier–Stokes (or Stokes) equation with a PDE for the tensor order parameter which indicates the preferred local orientation as well as the strength of the local alignment of the liquid crystal molecules. Well-posedness of the Beris–Edwards model in $\mathbb R^2$ and $\mathbb R^3$ was first studied in [44, 45]. Existence of weak and strong solutions in a bounded domain with a fixed boundary and both homogeneous and inhomogeneous boundary conditions for the tensor order parameter were established in [46, 47].

In our work, we consider a model which combines a Beris–Edwards liquid crystal with a squirmer. Such a system was, for example, used as a computational model in [48] to study orientation dynamics of the spherical squirmer with respect to the preferred orientation of the liquid crystal. In [49], we extended this model to elongated squirmers and studied how the long-term orientation dynamics of the squirmer depends on physical and geometrical parameters. To the best of our knowledge, there are no analytical results, such as well-posedness or model reductions via multi-scale limits for squirmers immersed in a liquid crystal environment. On the other hand, there are classical multi-scale results for particles in isotropic fluids, e.g., dilute, Brinkman’s and Darcy’s regimes via homogenisation limits by G. Allaire [50, 51]. We also refer to [52–56] for some recent works on analysis of fluid suspensions with solid particles.

The structure of this paper is as follows. In Section 2, we present the Beris–Edwards model coupled with a squirmer for both the time-dependent and steady-state problems. The latter corresponds to a squirmer moving with a constant velocity. In Section 3, we formulate our main results on the existence of solutions to both the steady-state and time-dependent problems as well as a two-scale homogenisation limit resulting in an effective model for a suspension of squirmers swimming parallel to each other. The last statement can be considered as a steady motion of a bacterial colony. Proof of the main results is presented in Sections 4, 5 and 6. Some calculations and non-dimensionalisation are relegated to Appendix.

2. Model

2.1. Time-dependent PDE system

Consider a rigid squirmer swimming in a liquid crystal with translational and angular velocities ${\boldsymbol{V}}(t)$ and ${\boldsymbol{\omega }}(t)$ , respectively. In the context of the Beris–Edwards model, the liquid crystal is described by a velocity field ${\boldsymbol{u}}({\boldsymbol{x}},t)\;:\;\mathbb R^d\times \mathbb R_+\mapsto \mathbb R^d$ and a tensor order field ${\boldsymbol{Q}}({\boldsymbol{x}},t)\;:\;\mathbb R^d\times \mathbb R_+\mapsto \mathbb R^{d\times d}$ taking values in symmetric traceless $d\times d$ matrices. Here, $d=2,3$ is the spatial dimension. The functions ${\boldsymbol{u}}=(u_j)_{j=1}^{d}$ and ${\boldsymbol{Q}}=(Q_{ij})_{i,j=1}^{d}$ satisfy the following system of PDEs and boundary conditions, written in the frame moving with velocity ${\boldsymbol{V}}(t)$ , so the squirmer is always centred at $0$ :

(1) \begin{equation} \rho (\partial _t +{\boldsymbol{u}}\cdot \nabla ){\boldsymbol{u}} + \rho \frac{d{\boldsymbol{V}}}{dt} = \nabla \cdot (\sigma _{\textrm{hydro}}+\sigma _{\textrm{ela}}),\textrm{ in } \Pi \setminus{\mathcal{P}}(t) \end{equation}

(2) \begin{equation} \nabla \cdot{\boldsymbol{u}} = 0,\textrm{ in }\Pi \setminus{\mathcal{P}}(t), \end{equation}

(3) \begin{equation}{\boldsymbol{u}} = u_{\textrm{sq}}(\boldsymbol{\alpha }(t),{\boldsymbol{x}})\boldsymbol{\tau } +{\boldsymbol{\omega }}(t) \times{\boldsymbol{x}},\textrm{ on }\partial{\mathcal{P}}(t), \end{equation}

(4) \begin{equation} \partial _t{\boldsymbol{Q}} + ({\boldsymbol{u}}\cdot \nabla ){\boldsymbol{Q}} -S(\nabla{\boldsymbol{u}},{\boldsymbol{Q}})= \Gamma \!\left ( K\Delta{\boldsymbol{Q}} + \hat{{\bf H}}({\boldsymbol{Q}}) \right ) + F_{\textrm{ext}}({\boldsymbol{Q}},{\boldsymbol{Q}}_\infty ),\textrm{ in }\Pi \setminus{\mathcal{P}}(t), \end{equation}

(5) \begin{equation}{\boldsymbol{Q}},{\boldsymbol{u}}, \nabla p \quad \textrm{periodic in} \, \Pi \end{equation}

(6) \begin{equation} K\partial _\nu{\boldsymbol{Q}} =W({\boldsymbol{Q}}_{\textrm{pref}}-{\boldsymbol{Q}}) \textrm{ on }\partial{\mathcal{P}}(t). \end{equation}

Here, $\Pi = (\!-\!L, L)^{d}$ is a periodic box, ${\mathcal{P}}(t)$ and $\partial{\mathcal{P}}(t)$ are the domain occupied by the squirmer and its surface in the moving frame. We assume that $\overline{\mathcal{P}(t)}\subset \Pi$ for all $t\geq 0$ . We will also use the notation $\Omega (t)\;:\!=\;\Pi \setminus \mathcal{P}(t)$ to denote the fluid region.

Equation (1) is a modified Navier–Stokes equation for the velocity ${\boldsymbol{u}}({\boldsymbol{x}},t)$ which satisfies the divergence-free condition (2). To this effect, we have $\sigma _{\textrm{hydro}}=\eta (\nabla{\boldsymbol{u}}+(\nabla{\boldsymbol{u}})^{\textrm{T}})-p\mathbb{I}$ to be the standard isotropic stress tensor where $p({\boldsymbol{x}},t)$ is the pressure of the liquid crystal with uniform density $\rho$ and viscosity $\eta$ . The internal structure of the liquid crystal, i.e., local preferred direction and order, affects the flow through an additional elastic stress $\sigma _{\textrm{ela}}$ given by

(7) \begin{eqnarray} \sigma _{\textrm{ela}}&=&K\!\left [({\boldsymbol{Q}}\,\Delta{\boldsymbol{Q}}- \Delta{\boldsymbol{Q}}\,{\boldsymbol{Q}})-\nabla{\boldsymbol{Q}}\odot \nabla{\boldsymbol{Q}}\right ]\nonumber \\[5pt] &&-\xi \!\left [{\bf H}({\boldsymbol{Q}}+\dfrac{\mathbb I}{d})+({\boldsymbol{Q}}+\dfrac{\mathbb{I}}{d}){\bf H}-2({\boldsymbol{Q}}+\frac{\mathbb{I}}{d})\textrm{Tr}({\boldsymbol{Q}}{\bf H})\right ]. \end{eqnarray}

Here, $K$ is the elastic constant and $\nabla{\boldsymbol{Q}}\odot \nabla{\boldsymbol{Q}}$ is a $d\times d$ matrix with the $(k,l)$ component $\sum \limits _{i,j}\partial _{x_k}Q_{ij}\partial _{x_l}Q_{ij}$ . The parameter $\xi$ measures the ratio between tumbling and aligning that a shear flow exerts on the liquid crystal molecules. The matrix-valued function ${\bf H}={\bf H}({\boldsymbol{Q}})$ is defined as ${\bf H}({\boldsymbol{Q}})=\hat{{\bf H}}({\boldsymbol{Q}})+K\Delta{\boldsymbol{Q}}$ , where $\hat{{\bf H}}({\boldsymbol{Q}})$ is

(8) \begin{equation} \hat{{\bf H}}({\boldsymbol{Q}})\;:\!=\;a{\boldsymbol{Q}}-c{\boldsymbol{Q}} \textrm{Tr}({\boldsymbol{Q}}^2)=-\nabla _{\boldsymbol{Q}}\!\left (\dfrac{c}{4}\|{\boldsymbol{Q}}\|^4-\dfrac{a}{2}\|{\boldsymbol{Q}}\|^2 \right )=-\nabla _{\boldsymbol{Q}}\hat{\mathcal{F}}({\boldsymbol{Q}}). \end{equation}

The scalar potential $\hat{\mathcal{F}}({\boldsymbol{Q}})$ is the polynomial part of the Landau–de Gennes free energy whose coefficients $a$ and $c$ depend on macroscopic parameters of the liquid crystal such as temperature. The potential $\hat{\mathcal{F}}({\boldsymbol{Q}})$ attains minima at ${\boldsymbol{Q}}={\bf 0}$ , corresponding to the isotropic state when the liquid crystal flows as a Newtonian fluid, and at tensor order parameters $\boldsymbol{Q}$ with $q_{\infty }\;:\!=\;\|{\boldsymbol{Q}}\|=\sqrt{\dfrac{a}{c}}$ , corresponding to the equilibrium liquid crystalline states.

Boundary conditions (3) describes how the squirmer interacts with the flow $\boldsymbol{u}$ of the liquid crystal. The orientation of the squirmer is described by a vector $\boldsymbol{\alpha }(t)\in \mathcal{S}^{d-1}$ . We also let $\boldsymbol{\tau }$ to be a tangent vector field to the surface of the squirmer which can be chosen to be $\boldsymbol{\tau }\;:\!=\;(\boldsymbol{\alpha }\times \boldsymbol{\nu })\times \boldsymbol{\nu }$ , where $\boldsymbol{\nu }$ is the inward normal vector on the squirmer’s surface $\partial \mathcal{P}(t)$ . A typical example of the slip velocity $u_{\textrm{sq}}$ is given by [35] (which is also used in the computational work [48, 49])

\begin{equation*} u_{\text {sq}}(\boldsymbol {\alpha }(t),\boldsymbol{x})=v_{\text {prop}}\sin \theta (1+\beta \cos (\theta )), \text { where }\theta =\cos ^{-1}\!\left [\dfrac {\boldsymbol{x}\cdot \boldsymbol {\alpha }}{\|\boldsymbol{x}\|}\right ]. \end{equation*}

Here, $\theta$ is the polar angle between the direction vector $\boldsymbol{\alpha }$ and the vector connecting the squirmer centre and the boundary point $\boldsymbol{x}$ , the parameter $v_{\textrm{prop}}$ is proportional to propulsion strength and $\beta$ quantifies the type of the squirmer (puller vs pusher; see [49] for details). In this work, we consider $u_{\textrm{sq}}(\boldsymbol{\alpha }(t),\boldsymbol{x})=\sin\!(\theta ) g(\theta )$ with a smooth function $g(\theta )$ . Note that such $u_{\textrm{sq}}$ vanishes at points of singularity of the vector field $\boldsymbol{\tau }$ .

The instantaneous angular velocity of the squirmer is denoted by ${\boldsymbol{\omega }}(t)$ . Then any material point $\boldsymbol{x}$ on the squirmer surface $\partial \mathcal{P}(t)$ will move with velocity ${\boldsymbol{\omega }}(t)\times{\boldsymbol{x}}$ in the moving frame for which the system (1)–(6) is written. Then the boundary condition (3) states that there is a given slip velocity $u_{\textrm{sq}}(\boldsymbol{\alpha }(t),{\boldsymbol{x}})\boldsymbol{\tau }$ with the no-penetration condition:

\begin{eqnarray*} \!\left [{\boldsymbol{u}}-{\boldsymbol{\omega }}\times{\boldsymbol{x}}\right ]\times \boldsymbol{\nu }&=&u_{\textrm{sq}}(\boldsymbol{\alpha }(t),{\boldsymbol{x}})\boldsymbol{\tau }\times \boldsymbol{\nu },\\[5pt] \!\left [{\boldsymbol{u}}-{\boldsymbol{\omega }}\times{\boldsymbol{x}}\right ]\cdot \boldsymbol{\nu }&=&0. \end{eqnarray*}

The given non-zero slip velocity models self-propulsion of the squirmer. Such a condition was originally derived for micro-organisms swimming with the help of small elastic appendages (cilia) distributed on the surface [35].

The matrix-valued equation (4) describes the dynamics of ${\boldsymbol{Q}}({\boldsymbol{x}},t)$ . While the two first terms in the left-hand side of (4) are the advection derivative, the third term $S(\nabla{\boldsymbol{u}},{\boldsymbol{Q}})$ describes how the flow gradient $\nabla{\boldsymbol{u}}$ rotates and stretches the order parameter $\boldsymbol{Q}$ and is given by

(9) \begin{eqnarray} S(\nabla{\boldsymbol{u}},{\boldsymbol{Q}})=(\xi{\boldsymbol{D}} +{\boldsymbol{A}})\!\left ({\boldsymbol{Q}}+\frac{\mathbb I}{d}\right )+ \left ({\boldsymbol{Q}}+\frac{\mathbb I}{d}\right )(\xi{\boldsymbol{D}}-{\boldsymbol{A}})-2\xi \!\left ({\boldsymbol{Q}}+\frac{\mathbb I}{d}\right )\textrm{tr}({\boldsymbol{Q}}\nabla{\boldsymbol{u}}), \end{eqnarray}

where ${\boldsymbol{D}}=\dfrac{1}{2}\!\left [\nabla{\boldsymbol{u}} +(\nabla{\boldsymbol{u}})^{\textrm{T}}\right ]$ and ${\boldsymbol{A}}=\dfrac{1}{2}\!\left [\nabla{\boldsymbol{u}} -(\nabla{\boldsymbol{u}})^{\textrm{T}}\right ]$ are symmetric and anti-symmetric parts of $\nabla{\boldsymbol{u}}$ , respectively. The right-hand side of (4) consists of the term leading to the minimisation of the total Landau–de Gennes energy

(10) \begin{equation} \mathcal{E}_{\textrm{LdG}}({\boldsymbol{Q}})=\int \limits _{\Omega (t)}\hat{\mathcal{F}}({\boldsymbol{Q}})+\frac{K}{2}|\nabla{\boldsymbol{Q}}|^2\,\textrm{d}{\boldsymbol{x}} \end{equation}

with the relaxation parameter $\Gamma \gt 0$ and the term $F_{\textrm{ext}}({\boldsymbol{Q}},{\boldsymbol{Q}}_\infty )$ describing the aligning effect with an external field. This term imposes the equilibrium condition for liquid crystal, that is, in the absence of squirmer we have ${\boldsymbol{Q}}\equiv{\boldsymbol{Q}}_\infty$ . We chose ${\boldsymbol{Q}}_{\infty }=q_{\infty }(\boldsymbol{e}_x\otimes \boldsymbol{e}_x-\frac{\mathbb I}{d})$ which means that if the liquid crystal is not perturbed by a squirmer then its molecules are oriented parallel to $\boldsymbol{e}_x$ (the unit basis vector parallel to $x$ -axis). In this work, we will use the example of $F_{\textrm{ext}}({\boldsymbol{Q}},{\boldsymbol{Q}}_\infty )$ from [21] for $d=2$ , given by

(11) \begin{equation} F_{\textrm{ext}}({\boldsymbol{Q}},{\boldsymbol{Q}}_{\infty }) = -\zeta{\boldsymbol{Q}} R_{\pi/2} \textrm{Tr}[{\boldsymbol{Q}}{\boldsymbol{Q}}_{\infty }R_{\pi/2}], \end{equation}

where $\zeta \geq 0$ is the aligning parameter and $R_{\pi/2}$ is the matrix of counterclockwise rotation by $\pi/2$ . For $d=3$ , the formula for $F_{\textrm{ext}}({\boldsymbol{Q}},{\boldsymbol{Q}}_\infty )$ is

(12) \begin{equation} F_{\textrm{ext}}({\boldsymbol{Q}},{\boldsymbol{Q}}_{\infty }) =\zeta \!\left ( \textrm{tr}({\boldsymbol{Q}}^2){\boldsymbol{Q}}_\infty - \textrm{tr}({\boldsymbol{Q}}{\boldsymbol{Q}}_{\infty }){\boldsymbol{Q}} \right ) \end{equation}

We note that if one considers dynamics $\dot{{\boldsymbol{Q}}}=F_{\textrm{ext}}({\boldsymbol{Q}},{\boldsymbol{Q}}_{\infty })$ then the Euclidean norm of $\boldsymbol{Q}$ is preserved, i.e. $\textrm{tr}({\boldsymbol{Q}}^2)\equiv \textrm{const}$ , and ${\boldsymbol{Q}}(t)$ converges to a multiple of ${\boldsymbol{Q}}_{\infty }$ as $t$ increases, so that ${\boldsymbol{Q}}\cdot{\boldsymbol{Q}}_{\infty }=\textrm{tr}({\boldsymbol{Q}}{\boldsymbol{Q}}_{\infty })\gt 0$ . One can also show that (12) is equivalent to (11) in the case of two-dimensional $\boldsymbol{Q}$ and ${\boldsymbol{Q}}_\infty$ (with zero third row and column).

We impose anchoring boundary condition (6) on $\boldsymbol{Q}$ along the squirmer surface $\partial \mathcal{P}(t)$ which forces $\boldsymbol{Q}$ to be close to a given tensor ${\boldsymbol{Q}}_{\textrm{pref}}=q_\infty (\boldsymbol{n}_{\textrm{pref}}\otimes \boldsymbol{n}_{\textrm{pref}}-\frac{{\mathbb I}}{d})$ . Here, $\boldsymbol{n}_{\textrm{pref}}=\boldsymbol{\nu }$ in the case of homeotropic anchoring when the surface orients liquid crystal molecules perpendicular to it or equivalently, parallel to the normal vector $\boldsymbol{\nu }$ . On the other hand, $\boldsymbol{n}_{\textrm{pref}}=\boldsymbol{\tau }$ in the case of the planar anchoring when molecules are aligned with the tangential vector field $\boldsymbol{\tau }$ . The boundary condition (6) indeed penalises the difference ${\boldsymbol{Q}}_{\textrm{pref}}-{\boldsymbol{Q}}$ in the sense that if we drop all terms in (4) except $\Gamma \!\left ( K\Delta{\boldsymbol{Q}} + \hat{{\bf H}}({\boldsymbol{Q}}) \right )$ , then the solution $\boldsymbol{Q}$ to this truncated version of (4) with boundary condition (6) minimises the energy

(13) \begin{equation} \mathcal{E}_{\textrm{LdG}}({\boldsymbol{Q}})+W\int \limits _{\partial \mathcal{P}(t)}|{\boldsymbol{Q}}_{\textrm{pref}}-{\boldsymbol{Q}}|^2\,\textrm{d}S_{\boldsymbol{x}}. \end{equation}

The coefficient $W$ in front of the penalisation term in the energy functional (13) and also the right-hand side of (6) measures the anchoring strength. Mathematically, depending on if $W\to \infty$ or $0$ , (6) reduces to Dirichlet or Neumann boundary condition for $\boldsymbol{Q}$ .

To determine the trajectory of the squirmer, that is, its velocity $\boldsymbol{V}(t)$ and angular velocity ${\boldsymbol{\omega }}(t)$ , we consider force and torque balances for the squirmer:

(14) \begin{equation} m\frac{\textrm{d}{\boldsymbol{V}}}{\textrm{d}t} = \int \limits _{\partial{\mathcal{P}}(t)} \sigma \boldsymbol{\nu }\, \textrm{d}S_{\boldsymbol{x}}, \end{equation}

(15) \begin{equation} \frac{\textrm{d}\!\left [I(t)\boldsymbol{\omega }\right ]}{\textrm{d}t} = \int \limits _{\partial{\mathcal{P}}(t)}{\boldsymbol{x}} \times \sigma \boldsymbol{\nu } +\boldsymbol{\ell } \,\textrm{d}S_{\boldsymbol{x}}. \end{equation}

Here, $\sigma =\sigma _{\textrm{hydro}}+\sigma _{\textrm{ela}}$ is the total stress, whereas $m$ and $I(t)= \left \{I_{ij}\right \}_{i,j=1}^{d}$ are the mass and inertia tensor of the squirmer, defined via

\begin{eqnarray*} m&=& \rho _{\mathcal{P}}|\mathcal{P}(t)|,\\[5pt] I_{ij}(t)&=&\rho _{\mathcal{P}}\int \limits _{\mathcal{P}(t)}\!\left [\boldsymbol{e}_i\times{\boldsymbol{x}} \right ]\cdot \left [\boldsymbol{e}_j\times{\boldsymbol{x}}\right ]\!\textrm{d}{\boldsymbol{x}}. \end{eqnarray*}

Here, $\rho _{\mathcal{P}}$ is the squirmer’s density. The additional torque $\boldsymbol{\ell }$ comes from the internal structure of the liquid crystal, namely from that there is a preferred direction. It translates into the non-zero asymmetric part of the stress tensor $\sigma$ . The formula for this additional torque is [49]

(16) \begin{equation} \boldsymbol{\ell }=\boldsymbol{\mu }\boldsymbol{\nu }, \textrm{ where }\boldsymbol{\mu }=(\mu _{ij})_{i,j=1}^{d} \textrm{ and } \mu _{ij}=-2K\sum \limits _{m,l,k=1}^{d}\epsilon _{ilk}Q_{lm}Q_{mk,j}. \end{equation}

Here, $\epsilon _{ilk}$ is the Levi–Civita symbol. Finally, we note that the orientation $\boldsymbol{\alpha }(t)$ and the angular velocity $\boldsymbol{\omega }$ are related via

(17) \begin{equation} \dot{\boldsymbol{\alpha }}={\boldsymbol{\omega }}\times \boldsymbol{\alpha }. \end{equation}

Remark 2.1. Note that the term $\ell$ admits a simplified form:

\begin{eqnarray*} \int \limits _{\partial \mathcal{P}}\ell _i\,\textrm{d}S_x&=&-2K\sum \limits _{l,m,k=1}^{d}\int _{\partial \mathcal{P}}\epsilon _{ilk}Q_{lm}Q_{mk,j}\nu _j\,\textrm{d}S_x\\[5pt] &=&-2W\sum \limits _{l,m,k=1}^{d}\int _{\partial \mathcal{P}}\epsilon _{ilk}Q_{lm}(Q_{\textrm{pref},mk}-Q_{mk})\,\textrm{d}S_x. \end{eqnarray*}

Here, we used boundary conditions (6). Next, for any symmetric matrix $\boldsymbol{B}=(B_{ij})_{i,j=1}^{d}$ we have

(18) \begin{equation} \sum \limits _{l,m,k}\epsilon _{ilk}B_{lm}B_{mk}=0. \end{equation}

Indeed, from properties of the Levi–Civita symbol we have

\begin{equation*} \sum \limits _{l,m,k}\epsilon _{ilk}B_{lm}B_{mk}=-\sum \limits _{l,m,k}\epsilon _{ilk}B_{km}B_{ml}. \end{equation*}

On the other hand, due to symmetry of $\boldsymbol{B}$ we have

\begin{equation*} \sum \limits _{l,m,k}\epsilon _{ilk}B_{lm}B_{mk}=\sum \limits _{l,m,k}\epsilon _{ilk}B_{km}B_{ml}. \end{equation*}

Thus, we have (18), from which we have the simplified form expression (simplified because it is linear in $\boldsymbol{Q}$ as opposed to (16) which is quadratic in $\boldsymbol{Q}$ ):

(19) \begin{eqnarray} \int \limits _{\partial \mathcal{P}}\ell _i\,\textrm{d}S_x&=&-2W\sum \limits _{l,m,k=1}^{d}\int \limits _{\partial \mathcal{P}}\epsilon _{ilk}Q_{lm}Q_{\textrm{pref},mk}\,\textrm{d}S_x. \end{eqnarray}

Remark 2.2. We end the introduction of the time-dependent problem with the energy identity satisfied by solutions of this problem. First, consider the energy functional:

(20) \begin{eqnarray} \mathcal{E}(t)&=&\dfrac{m\boldsymbol{V}^2}{2}+\dfrac{I\boldsymbol{{\boldsymbol{\omega }}}\cdot{\boldsymbol{\omega }}}{2}+\dfrac{\rho }{2}\int \limits _{\Omega (t)}|\boldsymbol{u}+\boldsymbol{V}|^2\,\textrm{d}{\boldsymbol{x}}\nonumber \\[5pt] &&+\int \limits _{\Omega (t)}\hat{\mathcal{F}}({\boldsymbol{Q}})+\frac{K}{2}|\nabla{\boldsymbol{Q}}|^2\,\textrm{d}{\boldsymbol{x}}+ \dfrac{ W}{2}\int \limits _{\partial \mathcal{P}(t)}|{\boldsymbol{Q}}_{\textrm{pref}}-{\boldsymbol{Q}}|^2\,\textrm{d}S_{\boldsymbol{x}}. \end{eqnarray}

Note that in the absence of the squirmer $\mathcal{P}(t)=\emptyset$ (or when the squirmer is passive, i.e., $u_{\textrm{sq}}=0$ ) and if the external field $F_{\textrm{ext}}$ equals zero, then the system is dissipative, that is, the energy is non-increasing:

(21) \begin{equation} \dfrac{\textrm{d}}{\textrm{d}t}\mathcal{E}(t)=-\mathcal{D}(t)\leq 0, \quad \textrm{where }\mathcal{D}(t)\;:\!=\; \eta \int \limits _{\Omega (t)}|\nabla{\boldsymbol{u}}|^2 \,\textrm{d}{\boldsymbol{x}}+\Gamma \int \limits _{\Omega (t)}|{\bf H}|^2\,\textrm{d}{\boldsymbol{x}}. \end{equation}

On the other hand, when the system experiences the energy input from the self-propulsion mechanism and external field $F_{\textrm{ext}}$ , the energy identity takes the following form:

(22) \begin{equation} \dfrac{\textrm{d}}{\textrm{d}t}\mathcal{E}(t) =-\mathcal{D}(t)+\int \limits _{\partial \mathcal{P}(t)}\sigma \boldsymbol{\nu }\cdot u_{\textrm{sq}}\boldsymbol{\tau }\,\textrm{d}S_{\boldsymbol{x}}+\int \limits _{\Omega (t)}{\bf H}\;:\;F_{\textrm{ext}}\,\textrm{d}{\boldsymbol{x}}. \end{equation}

Note that boundary integral $\int _{\partial \mathcal{P}(t)}\sigma \boldsymbol{\nu }\cdot u_{\textrm{sq}}\boldsymbol{\tau }\,\textrm{d}S_{\boldsymbol{x}}$ contains nonlinear terms in $\boldsymbol{Q}$ which in turn depends on the higher order regularity property of $\boldsymbol{Q}$ . This causes difficulty in the analysis of the time-dependent problem. Hence in this paper, we will only present a short-time existence result and leave the long-time behaviour to future work.

2.2. Steady-state PDE system

In this work, we are also interested in the steady translational motion of the squirmer in the liquid crystal. In the context of the model (1)–(6)(14)–(15), the steady motion is described by the stationary solution of this system:

(23) \begin{align} \rho{\boldsymbol{u}}\cdot \nabla{\boldsymbol{u}} = \nabla \cdot (\sigma _{\textrm{hydro}}+\sigma _{\textrm{ela}}),\textrm{ in } \Pi \setminus{\mathcal{P}}_{\textrm{st}}, \end{align}

(24) \begin{align} \nabla \cdot{\boldsymbol{u}} = 0,\textrm{ in }\Pi \setminus{\mathcal{P}}_{\textrm{st}}, \end{align}

(25) \begin{align} {\boldsymbol{u}} = u_{\textrm{sq}}(\boldsymbol{\alpha }_{\textrm{st}},{\boldsymbol{x}})\boldsymbol{\tau },\textrm{ on }\partial{\mathcal{P}}_{\textrm{st}}, \end{align}

(26) \begin{align} ({\boldsymbol{u}}\cdot \nabla ){\boldsymbol{Q}} -S(\nabla{\boldsymbol{u}},{\boldsymbol{Q}})= \Gamma \!\left ( K\Delta{\boldsymbol{Q}} + \hat{{\bf H}}({\boldsymbol{Q}}) \right ) + F_{\textrm{ext}}({\boldsymbol{Q}},{\boldsymbol{Q}}_\infty ),\textrm{ in }\Pi \setminus{\mathcal{P}}_{\textrm{st}}, \end{align}

(27) \begin{align} {\boldsymbol{Q}},{\boldsymbol{u}}, \nabla p \quad \textrm{periodic in} \, \Pi \end{align}

(28) \begin{align} K\partial _\nu{\boldsymbol{Q}} =W({\boldsymbol{Q}}_{\textrm{pref}}-{\boldsymbol{Q}}) \textrm{ on }\partial{\mathcal{P}}_{\textrm{st}}. \end{align}

Here, we assume that the squirmer moves with the velocity $V_{\textrm{st}}$ with the orientation angle $\boldsymbol{\alpha }_{\textrm{st}}$ , both of which are independent of time. As equations (23)–(28) are written in the squirmer’s frame, the domain $\mathcal{P}_{\textrm{st}}$ occupied by the squirmer will then be stationary. Similar to the time-dependent case, we use $\Omega =\Pi \setminus \mathcal{P}_{\textrm{st}}$ to denote the fluid region in the steady-state case.

In this setting, the force and torque balances (14), (15) become

(29) \begin{align} \textbf{0}=\int \limits _{\partial{\mathcal{P}}_{\textrm{st}}} \sigma \boldsymbol{\nu }\, \textrm{d}S_{\boldsymbol{x}}, \end{align}

(30) \begin{align} \textbf{0} = \int \limits _{\partial{\mathcal{P}}_{\textrm{st}}}{\boldsymbol{x}} \times \sigma \boldsymbol{\nu } +\boldsymbol{\ell } \,\textrm{d}S_{\boldsymbol{x}}. \end{align}

The force balance (29), in view of periodic boundary conditions for $\boldsymbol{u}$ and $\boldsymbol{Q}$ together with ${\boldsymbol{u}}\cdot \boldsymbol{\nu }=0$ on $\partial P_{\textrm{st}}$ (follows from (25)), leads to

(31) \begin{eqnarray} 0=\int \limits _{\partial{\mathcal{P}}_{\textrm{st}}} \sigma \boldsymbol{\nu }\, \textrm{d}S_{\boldsymbol{x}}=-\int \limits _{\partial \Pi }\sigma \boldsymbol{\nu }\, \textrm{d}S_{\boldsymbol{x}} +\rho \int \limits _{\Omega }{\boldsymbol{u}}\cdot \nabla{\boldsymbol{u}} \,\textrm{d}x=-\int \limits _{\partial \Pi } p\boldsymbol{\nu }\,\textrm{d}S_x. \end{eqnarray}

Therefore, since $\nabla p$ is periodic, as imposed in (27), it follows from (31) that $p({\boldsymbol{x}})$ is periodic in $\Pi$ . Indeed, the fact that $\nabla p$ is periodic implies that

(32) \begin{equation} p ({\boldsymbol{x}})= \boldsymbol{m} \cdot{\boldsymbol{x}}+ p^{\textrm{per}}({\boldsymbol{x}}), \end{equation}

where $p^{\textrm{per}}({\boldsymbol{x}})$ is a function which is periodic in $\Pi$ and $\boldsymbol{m}\in \mathbb R^{d}$ . Substitution of (32) into (31) implies that $\boldsymbol{m}=\textbf{0}$ and $p({\boldsymbol{x}})=p^{\textrm{per}}({\boldsymbol{x}})$ . In this case, the force balance (29) is satisfied regardless of squirmer’s velocity $V_{\textrm{st}}$ .

We note that if an external force $\textbf{F}^{(e)}= \left \{F^{(e)}_i\right \}_{i=1}^{d}$ is applied on the squirmer, then the force balance in stationary case becomes

(33) \begin{equation} \int \limits _{\partial{\mathcal{P}}_{\textrm{st}}} \sigma \boldsymbol{\nu }\, \textrm{d}S_{\boldsymbol{x}}+\textbf{F}^{(e)}=0 \end{equation}

which due to the same arguments as in derivation of (31) is equivalent to

(34) \begin{equation} -\int \limits _{\partial \Pi } p\boldsymbol{\nu }\, \textrm{d}S_{\boldsymbol{x}}+\textbf{F}^{(e)}=0. \end{equation}

Using (32) and the divergence theorem for the first term in the equation above, we get

(35) \begin{equation} |\Pi |\boldsymbol{m}=\textbf{F}^{(e)}. \end{equation}

Therefore, an external force results in the pressure difference

(36) \begin{equation} F^{(e)}_i=\dfrac{L}{2}[p]_i, \, i=1,\ldots,d, \quad \textrm{where }[p]_i=p|_{x_i=L}-p|_{x_i=-L}. \end{equation}

In terms of the force balance, the periodic problem (23)–(28) is in contrast with the analogous problem in the exterior domain $\Pi =\mathbb R^{d}$ . Namely, for the latter, we need to impose additional boundary conditions at $|{\boldsymbol{x}}| \to \infty$ : ${\boldsymbol{u}}=-\boldsymbol{V}_{\textrm{st}}$ and ${\boldsymbol{Q}}={\boldsymbol{Q}}_{\infty }$ , where $\boldsymbol{V}_{\textrm{st}}$ is the steady velocity of the squirmer. Then we would have obtained a Stokes-law-like force–velocity relation instead of the force–pressure relation (36). Another difference between the periodic problem (23)–(28) and its counterpart in $\Pi =\mathbb{R}^{d}$ is that the squirmer’s velocity $\boldsymbol{V}_{\textrm{st}}$ nowhere enters the problem neither in (23)–(28) nor in the force or torque balances (29)–(30). Physically, this is due to that an infinite periodic grid of squirmer without any connection to a reference point immovable in the inertial frame (as boundary conditions at $|{\boldsymbol{x}}|\to \infty$ would have provided) under no external influence can move with a steady velocity, of any direction and magnitude. Thus, the couple $(\boldsymbol{u}(\boldsymbol{x}+\boldsymbol{V}_{\textrm{st}}t)-\boldsymbol{V}_{\textrm{st}},{\boldsymbol{Q}}(\boldsymbol{x}-\boldsymbol{V}_{\textrm{st}}t))$ , where $(\boldsymbol{u},\boldsymbol{Q},\boldsymbol{\alpha }_{\textrm{st}})$ solve (23)–(28)(29)–(30), is a travelling wave solution of the time-dependent problem (1)–(6)(14)–(15) for any constant vector $\boldsymbol{V}_{\textrm{st}}$ . We note that applications of the model (1)–(6)(14)–(15) in [48, 49] concern with orientation squirmer dynamics for which periodic boundary conditions are sufficient. If one needs to study spatial dynamics of the squirmer, boundary conditions at $|{\boldsymbol{x}}|\to \infty$ are necessary.

In this work, the squirmer swims due to self-propulsion only, without an external force, $\textbf{F}^{(e)}=\textbf{0}$ . Thus, we impose periodicity for the pressure $p$ . Taking this into account, we define a weak solution of (23)–(28) as a couple $({\boldsymbol{u}},{\boldsymbol{Q}})\in H^1_{\textrm{per}}(\Omega ;\;\mathbb R^{d})\times H^2_{\textrm{per}}(\Omega ;\;\mathbb R^{d\times d})$ such that equations (24), (25) as well as the following two equalities hold for all $\boldsymbol{\psi }\in H^1_{\textrm{per}}(\Omega ;\;\mathbb R^d)\cap \!\left \{\boldsymbol{\psi }|_{\partial \mathcal{P}_{\textrm{st}}}=0 \textrm{ and }\nabla \cdot \boldsymbol{\psi }=0\right \}$ and $\boldsymbol{\Phi }\in H^1_{\textrm{per}}(\Omega ;\;\mathbb R^{d\times d})$ and every integral term is finite:

(37) \begin{eqnarray} &&\eta \int \limits _{\Omega }\nabla{\boldsymbol{u}}\;:\;\nabla \boldsymbol{\psi }\,\textrm{d}x + \rho \int \limits _{\Omega }({\boldsymbol{u}}\cdot \nabla ){\boldsymbol{u}}\cdot \boldsymbol{\psi }\,\textrm{d}x+\int \limits _{\Omega } \sigma _{\textrm{ela}}\;:\; \nabla \boldsymbol{\psi }\,\textrm{d}x=0. \end{eqnarray}

(38) \begin{eqnarray} &&\Gamma \!\left [\!-\!K\int \limits _{\Omega }\nabla{\boldsymbol{Q}}\;:\;\nabla \boldsymbol{\Phi } \,\textrm{d}x + W\int \limits _{\partial \mathcal{P}_{\textrm{st}}}({\boldsymbol{Q}}_{\textrm{pref}}-{\boldsymbol{Q}})\;:\;\boldsymbol{\Phi }\,\textrm{d}S_x\right ]\nonumber \\[5pt]&& \hspace{80 pt}- \int \limits _{\Omega }\!\left (({\boldsymbol{u}}\cdot \nabla ){\boldsymbol{Q}}-S(\nabla{\boldsymbol{u}},{\boldsymbol{Q}})\right )\;:\;\boldsymbol{\Phi }\,\textrm{d}x+ \int _{\Omega } F_{\textrm{ext}}\;:\;\boldsymbol{\Phi } \,\textrm{d}x=0. \end{eqnarray}

3. Main results

Here, we present our three main results.

Our first main result is the existence of a weak solution of the steady-state problem (23)–(28). Note that (23)–(28) is a boundary-value problem for the couple $(\boldsymbol{u},\boldsymbol{Q})$ . Force balance (29) results into periodic boundary conditions for pressure $p$ and torque balance (30) is not included in the system (23)–(28) so the orientation vector $\boldsymbol{\alpha }_{\textrm{st}}$ enters the problem (23)–(28) as a parameter. As discussed below in Remark 3.2, to find specific values of orientation vector $\boldsymbol{\alpha }_{st}$ with which the squirmer can swim steadily, one needs to additionally impose the torque balance (30). Provided that couples $(\boldsymbol{u},\boldsymbol{Q})$ solving (23)–(28) exist for all orientation vectors $\boldsymbol{\alpha }_{\textrm{st}}$ , the equation (30) possesses at least one solution $\boldsymbol{\alpha }_{\textrm{st}}$ , which corresponds to swimming along the LC preferred orientation $\boldsymbol{e}_{x}$ .

For the steady-state problem, we restrict ourselves to the case $\xi =0$ . In this case, we can establish the maximum principle formulated in Lemma 4.1. Under this simplification, we can represent $\sigma _{\textrm{ela}}$ as

(39) \begin{equation} \sigma _{\textrm{ela}}=-K\nabla{\boldsymbol{Q}} \odot \nabla{\boldsymbol{Q}} +\sigma _{a}, \textrm{ where }\sigma _a({\boldsymbol{Q}},{\bf H})={\boldsymbol{Q}}{\bf H}-{\bf H}{\boldsymbol{Q}}=K({\boldsymbol{Q}}\Delta{\boldsymbol{Q}}-\Delta{\boldsymbol{Q}}{\boldsymbol{Q}}) \end{equation}

and the term $S$ given by (9) satisfies the following equality:

(40) \begin{equation} S(\nabla{\boldsymbol{u}},{\boldsymbol{Q}})\;:\;{\boldsymbol{Q}} =\textrm{Tr}(S(\nabla{\boldsymbol{u}},{\boldsymbol{Q}}){\boldsymbol{Q}})=0. \end{equation}

We also impose

(41) \begin{equation} F_{\textrm{ext}}({\boldsymbol{Q}},{\boldsymbol{Q}}_{\infty })\;:\;{\boldsymbol{Q}}=0. \end{equation}

This condition holds for our specific choices of $F_{\textrm{ext}}({\boldsymbol{Q}},{\boldsymbol{Q}}_{\infty })$ given by (11) or (12).

Theorem 3.1. Suppose $\xi =0$ . There is a constant $C\gt 0$ independent of $K$ , $W$ , ${\boldsymbol{Q}}_{\textrm{pref}}$ , ${\boldsymbol{Q}}_{\infty }$ , $\eta$ , $\rho$ , $u_{\textrm{sq}}$ , $\Gamma$ , $\boldsymbol{\alpha }_{\textrm{st}}$ such that if

(42) \begin{equation} \eta \gt 2C\rho \|u_{\textrm{sq}}\|_{L^{\infty }(\Omega )}\textrm{ and }\Gamma \gt 2C\!\left (\dfrac{1}{\sqrt{K}}\|u_{\textrm{sq}}\|_{L^{\infty }(\Omega )}+\|u_{\textrm{sq}}\|_{W^{1,\infty }(\Omega )}\right ), \end{equation}

then there is a weak solution $({\boldsymbol{u}},{\boldsymbol{Q}}) \in (H^{1}(\Omega ), H^{2}(\Omega ))$ of (23)–(28).

Remark 3.1. The condition (42) holds when parameters $\eta$ and $\Gamma$ are sufficiently large, given all other parameters. The condition (42) also holds when $u_{\textrm{sq}}$ is sufficiently small which means that self-propulsion is small. In the limit $u_{\textrm{sq}}\to 0$ , we recover existence of steady state for a passive swimmer without a condition on parameters.

Remark 3.2. Theorem 3.1 states the existence of a weak solution of (23)–(28) for all orientation angles $\boldsymbol{\alpha }_{\textrm{st}}$ . As discussed in Section 4, the force balance (29) is satisfied since weak solutions of (23)–(28) have periodic pressure $p$ . To determine the steady orientation $\boldsymbol{\alpha }_{\textrm{st}}$ , one needs to consider additionally the torque balance (30) which is satisfied for $\boldsymbol{\alpha }_{\textrm{st}}=\dfrac{k\pi }{2}$ ( $k$ is an integer). We note that it follows from our numerical studies in [49] that a squirmer can swim steadily only if it is oriented parallel, $\boldsymbol{\alpha }_{\textrm{st}}=k\pi$ , or perpendicularly, $\boldsymbol{\alpha }_{\textrm{st}}=(2k-1)\frac{\pi }{2}$ , to the vector $\boldsymbol{e}_x$ , the liquid crystal orientation in the absence of the squirmer.

Our second main result is the local-in-time existence for the time-dependent problem (1)–(6) with (14) and (15). Here, we simplify the system by considering a spherical squirmer $\mathcal{P}(t)$ in its own moving frame so that $\Omega$ and $\mathcal{P}$ are independent of time. Under this assumption, the torque balance equation can be simplified into

(43) \begin{eqnarray} I\frac{\textrm{d}\boldsymbol{\omega }}{\textrm{d}t} = \int \limits _{\partial{\mathcal{P}}}{\boldsymbol{x}} \times \sigma \boldsymbol{\nu } + \boldsymbol{\ell }\,\textrm{d}S_{\boldsymbol{x}}, \end{eqnarray}

where the rotating inertia $I(t) = I\mathbb{I}$ becomes also independent of time and isotropic.

Theorem 3.2. Suppose that $({\boldsymbol{u}}_{\textrm{sq}},{\boldsymbol{Q}}_{\textrm{pref}}) \in H^{5/2}(\partial{\mathcal{P}}) \times H^{5/2}(\partial{\mathcal{P}})$ , $\xi$ is not necessarily 0, and the initial data $({\boldsymbol{u}}_0,{\boldsymbol{Q}}_0) \in H^{2}_{\sigma }(\Omega ) \times H^{3}(\Omega )$ , where $H_{\sigma }^2(\Omega )=H^2(\Omega )\cap \!\left \{\nabla \cdot{\boldsymbol{u}}=0 \right \}$ . Then there exists $T \gt 0$ and a unique solution $({\boldsymbol{u}},{\boldsymbol{Q}})$ to the system (1)–(6) with (14) and (43) such that

\begin{eqnarray*} &&{\boldsymbol{u}} \in H^{1}(0, T;\;H^{2}_{\sigma }(\Omega )) \cap H^{2}(0, T;\;L^{2}_{\sigma }(\Omega )), \\[5pt] &&{\boldsymbol{Q}} \in H^{1}(0, T;\;H^{3}(\Omega )) \cap H^{2}(0, T;\;H^{1}(\Omega )). \end{eqnarray*}

Remark 3.3. We adapt techniques from [47] to prove this result in Section 5. The main idea is to rewrite the problem in a suitable Banach space and then use the Banach’s fixed point theorem. However, the difference from [47] is an additional difficulty coming from the presence of the squirmer which requires to consider inhomogeneous boundary conditions as well as force and torque balances (14) and (43). The terms in balance equations involve boundary integrals with derivatives in integrands. It led to that the spatial regularity of the solution couple $({\boldsymbol{u}},{\boldsymbol{Q}})$ is higher than it is required by a weak solution of the PDE problem (1)–(6).

Our third main result is a formal homogenisation limit in the system (23)–(28). This result can be considered as the derivation of a simplified model describing motion of a colony with periodically distributed squirmers (e.g. bacterial colony) in the liquid crystal.

Specifically, we introduce a small parameter $ \varepsilon \;:\!=\;\frac{L}{\delta _L},$ where $L$ is the linear size of a periodic box containing a single squirmer and $\delta _L$ is the observation scale. Next, we consider the problem (23)–(28) where all the parameters are written in physical dimensions. Details of non-dimensionalisation are relegated to Appendix C. After the non-dimensionalisation, we consider the steady-state problem (23)–(28) in a periodic box $\Pi _{\varepsilon }=[\!-\!\varepsilon,\varepsilon ]^d$ . The squirmer occupies domain $\mathcal{P}_{\varepsilon }$ whose linear size is $\sim \varepsilon$ . Consider the domain $U$ which is $\mathbb R^{d}$ or a sub-domain of $\mathbb R^d$ composed of many periodic boxes $\Pi _{\varepsilon }$ such that the linear size of $U$ is of the order $1$ with respect to $\varepsilon$ . Then (23)–(28) becomes (see Appendix C for details):

(44) \begin{align} \varepsilon \gamma \Delta{\boldsymbol{Q}}+\tilde{a}\,{\boldsymbol{Q}}-\tilde{c}\,{\boldsymbol{Q}}\textrm{Tr}({\boldsymbol{Q}}^2)+S(\nabla \tilde{{\boldsymbol{u}}},{\boldsymbol{Q}}) - \tilde{{\boldsymbol{u}}} \cdot \nabla{\boldsymbol{Q}}+\tilde{\zeta }\, \tilde{F}_{\textrm{ext}}={\bf G}({\boldsymbol{x}})\textrm{ in }\Omega _{\varepsilon }, \end{align}

(45) \begin{align} \partial _{\boldsymbol{\nu }}{\boldsymbol{Q}} = \tilde{W}({\boldsymbol{Q}}_{\textrm{pref}} -{\boldsymbol{Q}}) \textrm{ on } \partial{\mathcal{P}}_{\varepsilon }, \end{align}

(46) \begin{align} \varepsilon \tilde{\rho }(\tilde{\boldsymbol{u}}\cdot \nabla )\tilde{{\boldsymbol{u}}}-\varepsilon \tilde{\eta }\Delta \tilde{{\boldsymbol{u}}} + \nabla \tilde{p} = \varepsilon ^2 \kappa \nabla \cdot (\nabla{{\boldsymbol{Q}}} \odot \nabla{{\boldsymbol{Q}}} +{{\boldsymbol{Q}}}\Delta{{\boldsymbol{Q}}} - \Delta{{\boldsymbol{Q}}}{{\boldsymbol{Q}}})+{\bf F}({\boldsymbol{x}}) \textrm{ in }\Omega _{\varepsilon }, \end{align}

(47) \begin{align} \tilde{{\boldsymbol{u}}} = \varepsilon \tilde{u}_{\textrm{sq}}\boldsymbol{\tau } \textrm{ on } \partial{\mathcal{P}}_{\varepsilon }. \end{align}

Here, $\Omega _{\varepsilon } = \Pi _{\varepsilon } \setminus \mathcal{P}_{\varepsilon }$ and $\nabla \cdot \tilde{{\boldsymbol{u}}}=0$ . Here, for simplicity, we assume that the orientation vector $\boldsymbol{\alpha }_{\textrm{st}}$ is independent of $\varepsilon$ . This physically means that the bacterial colony has reached the steady state when every squirmer swims along a stable direction which is related to the preferred direction of liquid crystal $\boldsymbol{e}_x$ and is independent of $\varepsilon$ (e.g. $\boldsymbol{\alpha }_{\textrm{st}}=\boldsymbol{e}_x$ ). $\bf G$ and $\bf F$ are given external fields, varying spatially at the scale $1$ (independent of $\varepsilon$ ). Parameters $\gamma,\tilde{a},\tilde{c},\tilde{\zeta },\tilde{W},\tilde{\rho },\tilde{\eta },\kappa$ are explained in Appendix C.

Our contribution in this regard is the identification of the homogenised limit $({\boldsymbol{u}}^{(h)},{\boldsymbol{Q}}^{(h)})$ of $(\varepsilon ^{-1}\tilde{{\boldsymbol{u}}},{\boldsymbol{Q}})$ . We relegate the presentation of the limiting equations as well as their derivation via formal two-scale asymptotic expansions to Section 6. We comment here that ${\boldsymbol{Q}}^{(h)}$ solves an algebraic equation (187), whereas ${\boldsymbol{u}}^{(h)}$ admits the representation (197) similar to that in the Darcy’s law.

4. Existence of steady state $-$ proof of Theorem 3.1

In this section, we address solvability of steady-state PDE system (23)–(28). To this end, we first show in Subsection 4.1 that if a solution of the system (23)–(28) exists (with $\xi =0$ ), then it satisfies a maximum principle for $|{\boldsymbol{Q}}|$ . Next, in Subsections 4.2, 4.3, 4.4 and 4.5, we prove the existence for the system (23)–(28) where nonlinearities $\hat{{\bf H}}$ and $F_{\textrm{ext}}$ are truncated for large values of $\boldsymbol{Q}$ . Finally, combination of the maximum principle and solvability of the truncated system implies the existence of a solution to the original system (23)–(28).

4.1. $L^{\infty }$ -bound on $\boldsymbol{Q}$

Here, we adapt the strategy from [57]. First, we introduce the number $q_*\gt 0$ such that

(48) \begin{equation} \Gamma \hat{{\bf H}}({\boldsymbol{Q}})\;:\;{\boldsymbol{Q}}\leq 0 \textrm{ for all }|{\boldsymbol{Q}}|\geq q_*. \end{equation}

Such a finite number $q_{*}$ exists since $F_{\textrm{ext}}({\boldsymbol{Q}},{\boldsymbol{Q}}_{\infty })$ is a quadratic polynomial of $\boldsymbol{Q}$ whereas $\hat{{\bf H}}$ is the third-order polynomial with a definite negative sign in front the highest power.

Lemma 4.1. Let $({\boldsymbol{u}},{\boldsymbol{Q}})$ be a solution of (23)–(28). Then $\|{\boldsymbol{Q}}\|_{L^{\infty }} \leq \alpha$ , where

(49) \begin{equation} \alpha \;:\!=\; \max \!\left \{|{\boldsymbol{Q}}_{\textrm{pref}}|,q_*\right \}. \end{equation}

Proof. Recall the equation for $\boldsymbol{Q}$ :

(50) \begin{equation} ({\boldsymbol{u}}\cdot \nabla ){\boldsymbol{Q}} - S(\nabla{\boldsymbol{u}},{\boldsymbol{Q}})-K\Gamma \Delta{\boldsymbol{Q}}-\Gamma \hat{{\bf H}}({\boldsymbol{Q}})- F_{\textrm{ext}}({\boldsymbol{Q}},{\boldsymbol{Q}}_\infty )=0. \end{equation}

By multiplying the above by $\boldsymbol{Q}$ , taking the trace of the resulting expression and using (40) and (41), we get

\begin{equation*} \dfrac {1}{2}{\boldsymbol{u}}\cdot \nabla (|{\boldsymbol{Q}}|^2)-\dfrac {\Gamma K}{2}\!\left (\Delta (|{\boldsymbol{Q}}|^2) - 2 |\nabla {\boldsymbol{Q}}|^2\right )-\Gamma \hat {{\bf H}}({\boldsymbol{Q}})\;:\;{\boldsymbol{Q}}=0. \end{equation*}

As $|\nabla{\boldsymbol{Q}}|^2$ is non-negative, we obtain the inequality:

(51) \begin{equation} {\boldsymbol{u}}\cdot \nabla (|{\boldsymbol{Q}}|^2)-\Gamma K\Delta (|{\boldsymbol{Q}}|^2) -2\Gamma \hat{{\bf H}}({\boldsymbol{Q}})\;:\;{\boldsymbol{Q}}\leq 0. \end{equation}

Now introduce $\psi ({\boldsymbol{Q}})\;:\!=\; (|{\boldsymbol{Q}}|^2- \alpha ^2)_{+}$ ( $\alpha$ is from (49)). Note that

\begin{equation*} \psi ({\boldsymbol{Q}})\,\mathcal {D}(|{\boldsymbol{Q}}|^2)=\psi ({\boldsymbol{Q}})\,\mathcal {D}(|{\boldsymbol{Q}}|^2-\alpha ^2)=\psi ({\boldsymbol{Q}})\,\mathcal {D}\psi ({\boldsymbol{Q}}), \end{equation*}

where $\mathcal{D}$ is either $\Delta$ or $\nabla$ .

Next, we multiply (51) by $\psi ({\boldsymbol{Q}})$ and integrate over $\Pi \setminus \mathcal{P}_{\textrm{st}}$ . Then, we have

(52) \begin{eqnarray} &&\dfrac{1}{2}\int \limits _{\partial{\mathcal{P}}}{\boldsymbol{u}}\cdot \nu \psi ^2\,\textrm{d}s_x-\Gamma K \int \limits _{\partial{\mathcal{P}}}\dfrac{\partial \psi ({\boldsymbol{Q}})}{\partial \nu }\psi ({\boldsymbol{Q}})\,\textrm{d}s_x\nonumber \\[5pt] &&\hspace{60 pt}+\;\Gamma K\|\nabla \psi ({\boldsymbol{Q}})\|^2_{L^2(\Omega )}-2\Gamma \int \limits _{\Omega }(\hat{{\bf H}}({\boldsymbol{Q}})\;:\;{\boldsymbol{Q}})\psi ({\boldsymbol{Q}})\,\textrm{d}x\leq 0. \end{eqnarray}

The first term in the left-hand side of the above inequality vanishes due to (25) while the second term is negative:

\begin{eqnarray*} K\int _{\partial{\mathcal{P}}} \dfrac{\partial \psi ({\boldsymbol{Q}})}{\partial \nu }\psi ({\boldsymbol{Q}})\,\textrm{d}s_x&=&2K\int _{\partial{\mathcal{P}}\cap \{|{\boldsymbol{Q}}|\gt \alpha \}} (\dfrac{\partial{\boldsymbol{Q}}}{\partial \nu }\;:\;{\boldsymbol{Q}})(|{\boldsymbol{Q}}|^2-\alpha ^2)\,\textrm{d}s_x\\[5pt] &=&2W\int _{\partial{\mathcal{P}}\cap \{|{\boldsymbol{Q}}|\gt \alpha \}} (({\boldsymbol{Q}}_{\textrm{pref}}-{\boldsymbol{Q}})\;:\;{\boldsymbol{Q}})(|{\boldsymbol{Q}}|^2-\alpha ^2)\,\textrm{d}s_x\\[5pt] &=&2W\int _{\partial{\mathcal{P}}\cap \{|{\boldsymbol{Q}}|\gt \alpha \}} (({\boldsymbol{Q}}_{\textrm{pref}}\;:\;{\boldsymbol{Q}})-|{\boldsymbol{Q}}|^2)(|{\boldsymbol{Q}}|^2-\alpha ^2)\,\textrm{d}s_x\\[5pt] &\leq &W\int _{\partial{\mathcal{P}}\cap \{|{\boldsymbol{Q}}|\gt \alpha \}} (|{\boldsymbol{Q}}_{\textrm{pref}}|^2-|{\boldsymbol{Q}}|^2)(|{\boldsymbol{Q}}|^2-\alpha ^2)\,\textrm{d}s_x\\[5pt] &\leq & 0. \end{eqnarray*}

Hence,

\begin{equation*} K\|\nabla \psi ({\boldsymbol{Q}})\|^2_{L^2(\Omega )}\leq 2 \int \limits _{\Omega }\!\left (\hat {{\bf H}}({\boldsymbol{Q}})\;:\;{\boldsymbol{Q}}\right )\psi ({\boldsymbol{Q}})\,\text {d}x. \end{equation*}

Next, by (48), we have

\begin{eqnarray*} K\|\nabla \psi ({\boldsymbol{Q}})\|^2_{L^2(\Omega )}&\leq & 2 \int \limits _{\Omega }\!\left (\hat{{\bf H}}({\boldsymbol{Q}})\;:\;{\boldsymbol{Q}}\right )\psi ({\boldsymbol{Q}})\,\textrm{d}x\leq 0. \end{eqnarray*}

so that $\|\nabla \psi ({\boldsymbol{Q}})\|^2_{L^2(\Omega )}=0$ . The lemma is thus proved.

4.2. Galerkin approximation for pair $({\boldsymbol{u}},{\bf H})$

We introduce here Galerkin approximations for the system (23)–(28). For each $m\in \mathbb N$ , we define

(53) \begin{equation}{\boldsymbol{u}}_m ={\boldsymbol{u}}_{\textrm{os}}+\hat{{\boldsymbol{u}}}_m={\boldsymbol{u}}_{\textrm{os}}+\sum \limits _{k=1}^{m}u_{km}\Psi _{k}\textrm{ and }{\bf H}_{m}=\sum \limits _{k=1}^{m}h_{km}\Phi _{k}. \end{equation}

Note that the domain $\Pi$ is a bounded periodic box and both Laplacian and Stokes operators have a discrete spectrum implying existence of bases:

(54) \begin{equation} \!\left \{\Psi _{k}\!\left |\, \Psi _k|_{\partial{\mathcal{P}}_{\textrm{st}}}=0,\, \nabla \cdot \Psi _k=0,\,\Psi _k\,\textrm{is} \Pi -\textrm{periodic} \right .\right \}_{k=1}^{\infty }\, \textrm{and}\, \!\left \{\Phi _{k}\!\left |\, \Phi _k\,\textrm{is}\, \Pi -\textrm{periodic} \right .\right \}_{k=1}^{\infty } \end{equation}

in $L^2_\sigma (\Omega ;\;\mathbb R^{d})$ and $L^2(\Omega ;\;\mathbb R^{d\times d})$ , respectively. (Recall that $\Omega = \Pi \setminus \mathcal{P}_{\textrm{st}}$ and $L^2_\sigma$ means $L^2$ -space with divergence-free condition.)

The function ${\boldsymbol{u}}_{\textrm{os}}$ above is an offset function used to take care of non-zero boundary conditions for $\boldsymbol{u}$ . It solves Stokes equation:

(55) \begin{align} \,\qquad\qquad\qquad\qquad\qquad\qquad\qquad \left\{\begin{array}{l@{\qquad\qquad\qquad\qquad\qquad\qquad\qquad}l} \eta \Delta{\boldsymbol{u}}_{\textrm{os}}+\nabla p_{\textrm{os}} = 0, \textrm{ in }\Omega & (55)\\[5pt] \nabla \cdot{\boldsymbol{u}}_{\textrm{os}}=0, \textrm{ in }\Omega, & (56)\\[5pt]{\boldsymbol{u}}_{\textrm{os}}=u_{\textrm{sq}}(\boldsymbol{\alpha }_{\textrm{st}},{\boldsymbol{x}})\boldsymbol{\tau } \textrm{ on }\partial \mathcal{P}_{\textrm{st}},& (57)\\[5pt]{\boldsymbol{u}}_{\textrm{os}}, p_{\textrm{os}} \, \textrm{periodic in} \, \Pi. & (58) \end{array}\right. \nonumber\end{align}

Anticipating that ${\boldsymbol{u}}={\boldsymbol{u}}_{\textrm{os}} + \hat{{\boldsymbol{u}}}$ , from (57) and (25), we have

(59) \begin{equation} \hat{u}=0\textrm{ on }\partial \mathcal{P}_{\textrm{st}}. \end{equation}

To continue, for an appropriately large constant $M\gt 0$ , we introduce a truncated potential $\hat{\mathcal{F}}_M\geq 0$ as follows:

(60) \begin{equation} \hat{\mathcal{F}}_{M}({\boldsymbol{Q}})= \left \{ \begin{array}{l@{\quad}l} \hat{\mathcal{F}}({\boldsymbol{Q}}), & \textrm{for $|{\boldsymbol{Q}}|\leq M$}\\[5pt] \!\left .\dfrac{\partial \hat{\mathcal{F}}}{\partial |{\boldsymbol{Q}}|}\right |_{|{\boldsymbol{Q}}|=M}|{\boldsymbol{Q}}| + \!\left (\!\left .\hat{\mathcal{F}}\right |_{|{\boldsymbol{Q}}|=M}-\left .\dfrac{\partial \hat{\mathcal{F}}}{\partial |{\boldsymbol{Q}}|}\right |_{|{\boldsymbol{Q}}|=M}M\right ) &\textrm{for }|{\boldsymbol{Q}}|\gt M \end{array} \right. . \end{equation}

or, more explicitly,

\begin{equation*} \hat {\mathcal {F}}_{M}(\boldsymbol{Q})= \left \{ \begin {array}{l@{\quad}l} \dfrac {c}{4}|\boldsymbol{Q}|^4-\dfrac {a}{2}|\boldsymbol{Q}|^2, & \text {for $|{\boldsymbol{Q}}|\leq M$}\\[5pt] M(cM^2 - a)(|\boldsymbol{Q}|-M)+\dfrac {c}{4}M^4-\dfrac {a^2}{2}M^2,\quad &\text {for }|{\boldsymbol{Q}}|\gt M \end {array} \right. . \end{equation*}

The functional derivative of $\hat{\mathcal{F}}_{M}$ is given by

(61) \begin{equation} \hat{{\bf H}}_M({\boldsymbol{Q}})=-\nabla _{{\boldsymbol{Q}}}\hat{\mathcal{F}}_M= \left \{ \begin{array}{l@{\quad}r} \hat{{\bf H}}({\boldsymbol{Q}}),&|{\boldsymbol{Q}}|\leq M,\\[5pt] \gamma _M\dfrac{{\boldsymbol{Q}}}{|{\boldsymbol{Q}}|},&|{\boldsymbol{Q}}|\gt M, \end{array}\right. \quad \textrm{where $\gamma _M = \!\left .\dfrac{\partial \hat{\mathcal{F}}}{\partial |{\boldsymbol{Q}}|}\right |_{|{\boldsymbol{Q}}|=M}$.} \end{equation}

We have the following bound on $\hat{H}_M({\boldsymbol{Q}})$ :

(62) \begin{equation} \|\hat{H}_M({\boldsymbol{Q}})\|_{L^{\infty }(\Omega )}\leq \Gamma _M, \quad \textrm{where }\Gamma _M=\max \{\|\hat{H}\|_{L^{\infty }(B_M(0))},\gamma _M\}. \end{equation}

We now define the function ${\boldsymbol{Q}}_m$ , corresponding to the Galerkin approximation ${\bf H}_m$ as the solution to the following system:

(63) \begin{equation} \!\left \{ \begin{array}{l} K\Delta{\boldsymbol{Q}}_m+\hat{{\bf H}}_M({\boldsymbol{Q}}_m)={\bf H}_m, \textrm{ in }\Omega \\[5pt] K\partial _\nu{\boldsymbol{Q}}_m= W ({\boldsymbol{Q}}_{\textrm{pref}}-{\boldsymbol{Q}}_m)\textrm{ on }\partial \mathcal{P}_{\textrm{st}},\\[5pt] \quad{\boldsymbol{Q}}_m \, \textrm{periodic in} \, \Pi \end{array} \right. \end{equation}

Below, we will need a priori estimates for the solution to the problem (63), formulated in the following lemma. Its proof is given in Appendix A.

Lemma 4.2. Let ${\bf H}_m\in L^2(\Omega )$ . Then there exists a solution ${\boldsymbol{Q}}_m$ for (63). Moreover, there exists a constant $C\gt 0$ such that

(64) \begin{align} \sqrt{K}\|\nabla{\boldsymbol{Q}}_m\|_{L^2(\Omega )}+\|{\boldsymbol{Q}}_m\|_{L^2(\Omega )}&+\sqrt{W}\|{\boldsymbol{Q}}_m\|_{L^2(\partial \mathcal{P}_{\textrm{st}})} \nonumber \\[5pt] &\leq C\!\left (\|{\bf H}_m\|_{L^2(\Omega )}+\sqrt{W}\|{\boldsymbol{Q}}_{\textrm{pref}}\|_{L^2(\partial \mathcal{P}_{\textrm{st}})}+1\right ) \end{align}

(65) \begin{align} \|{\boldsymbol{Q}}_m\|_{H^2(\Omega )}\leq C \!\left (\gamma _1\|{\bf H}_m\|_{L^2(\Omega )}+\gamma _2\|{\boldsymbol{Q}}_{\textrm{pref}}\|_{C^1}+\gamma _3\Gamma _M\right ), \end{align}

where

\begin{equation*} \gamma _1=\gamma _3=\dfrac{W+K}{W}\quad \textrm{and}\quad \gamma _2=\sqrt{\dfrac{W+K}{K}}. \end{equation*}

We define Galerkin approximations ( $\boldsymbol{u}_m, \textbf{H}_m$ ) as solutions of the two equations below for each $k=1,...,m$ :

(66) \begin{eqnarray} &&\eta \int \limits _{\Omega }\nabla \hat{{\boldsymbol{u}}}_m\;:\;\nabla \Psi _k\,\textrm{d}x +\eta \int \limits _{\Omega }\nabla{\boldsymbol{u}}_{\textrm{os}}\;:\;\nabla \Psi _k\,\textrm{d}x+ \rho \int \limits _{\Omega }({\boldsymbol{u}}_{m}\cdot \nabla ){\boldsymbol{u}}_m\cdot \Psi _k\,\textrm{d}x\nonumber \\[5pt]&& \hspace{120 pt}-\int \limits _{\Omega } \!\left (K\nabla{\boldsymbol{Q}}_m\odot \nabla{\boldsymbol{Q}}_m-\sigma _a({\boldsymbol{Q}}_m,{\bf H}_m)\right )\;:\; \nabla \Psi _{k}\,\textrm{d}x=0, \end{eqnarray}

(67) \begin{eqnarray} &&\Gamma \int \limits _{\Omega }{\bf H}_{m}\;:\;\Phi _k \,\textrm{d}x - \int \limits _{\Omega }\!\left (({\boldsymbol{u}}_m\cdot \nabla ){\boldsymbol{Q}}_{m}-S(\nabla{\boldsymbol{u}}_m,{\boldsymbol{Q}}_m)\right )\;:\;\Phi _k\,\textrm{d}x\nonumber \\[5pt] &&\hspace{220 pt}+ \int _{\Omega } F_{\textrm{ext},M} : \Phi _{k} \,\textrm{d}x=0. \end{eqnarray}

Here, $F_{\textrm{ext},M}$ is defined as a continuous function such that:

(68) \begin{equation} F_{\textrm{ext},M}= \left \{ \begin{array}{r@{\quad}l} F_{\textrm{ext}},& |F_{\textrm{ext}}|\leq M, \\[5pt] M\dfrac{F_{\textrm{ext}}}{|F_{\textrm{ext}}|},& |F_{\textrm{ext}}|\gt M. \end{array} \right. \end{equation}

Next, we will prove the existence and a priori estimates for $({\boldsymbol{u}}_m,{\bf H}_m$ ).

4.3. Energy estimate for Galerkin approximations

Lemma 4.3. There is a constant $C_0\gt 0$ independent of $K$ , $W$ , ${\boldsymbol{Q}}_{\textrm{pref}}$ , ${\boldsymbol{Q}}_{\infty }$ , $\eta$ , $\rho$ , $u_{\textrm{sq}}$ , $\Gamma$ , $\boldsymbol{\alpha }_{\textrm{st}}$ such that if $ \eta \gt 2C_0\rho \|u_{\textrm{sq}}\|_{L^{\infty }(\Omega )}\textrm{ and }\Gamma \gt 2C_0\!\left (\dfrac{1}{\sqrt{K}}\|u_{\textrm{sq}}\|_{L^{\infty }(\Omega )}+\|u_{\textrm{sq}}\|_{W^{1,\infty }(\Omega )}\right ),$ then

(69) \begin{equation} \|\nabla \hat{{\boldsymbol{u}}}_m\|_{L^2(\Omega )}^2+\|{\bf H}_m\|_{L^2(\Omega )}^2 \lt C. \end{equation}

Proof. Using test function $\hat{{\boldsymbol{u}}}_m$ and ${\bf H}_m$ instead of $\Psi _{k}$ and $\Phi _{k}$ in (66)–(67) and taking the sum of two equalities, we obtain the following energy equality:

(70) \begin{eqnarray} &&\eta \int \limits _{\Omega }|\nabla \hat{{\boldsymbol{u}}}_m|^2\,\textrm{d}x+\Gamma \int \limits _{\Omega }|{\bf H}_m|^2\,\textrm{d}x\nonumber \\[5pt] && \hspace{40 pt}= - \rho \int \limits _{\Omega }({\boldsymbol{u}}_{m}\cdot \nabla ){\boldsymbol{u}}_m\cdot \hat{{\boldsymbol{u}}}_m\,\textrm{d}x-\eta \int \limits _{\Omega }\nabla{\boldsymbol{u}}_{\textrm{os}}\;:\;\nabla \hat{{\boldsymbol{u}}}_m\,\textrm{d}x\nonumber \\[5pt] && \hspace{50 pt}-\int \limits _{\Omega }\!\left (\sigma _a({\boldsymbol{Q}}_m,{\bf H}_m)\;:\; \nabla \hat{{\boldsymbol{u}}}_m+S(\nabla \hat{{\boldsymbol{u}}}_m,{\boldsymbol{Q}}_m)\;:\;{\bf H}_m\right )\,\textrm{d}x\nonumber \\[5pt] && \hspace{50 pt}+\;K\int \limits _{\Omega }\!\left (\nabla{\boldsymbol{Q}}_m\odot \nabla{\boldsymbol{Q}}_m\;:\;\nabla \hat{{\boldsymbol{u}}}_m + (\hat{{\boldsymbol{u}}}_m\cdot \nabla ){\boldsymbol{Q}}_{m}\;:\;\Delta{{\boldsymbol{Q}}}_m\right ) \,\textrm{d}x\nonumber \\[5pt] && \hspace{50 pt} +\;K\int \limits _{\Omega } ({\boldsymbol{u}}_{\textrm{os}}\cdot \nabla ){\boldsymbol{Q}}_{m}\;:\;\Delta{{\boldsymbol{Q}}}_m \,\textrm{d}x +\int \limits _{\Omega }({\boldsymbol{u}}_m\cdot \nabla ){\boldsymbol{Q}}_m\;:\;\hat{{\bf H}}_M({\boldsymbol{Q}}_m)\,\textrm{d}x\nonumber \\[5pt] && \hspace{50 pt} - \int \limits _{\Omega } S(\nabla{\boldsymbol{u}}_{\textrm{os}},{\boldsymbol{Q}}_m)\;:\;{\bf H}_m \textrm{d}x - \int \limits _{\Omega } F_{\textrm{ext},M}\;:\;{\bf H}_m \,\textrm{d}x. \end{eqnarray}

Next, we estimate each integral. Below, $C$ denotes a generic constant independent of $K,W,{\boldsymbol{Q}}_{\textrm{pref}},\eta,\rho,{\boldsymbol{u}}_{\textrm{os}},\Gamma,m$ which may change from line to line, whereas $C^{*}$ is a generic constant which is independent of $m$ only and may also change from line to line.

1. We use the representation ${\boldsymbol{u}}_m={\boldsymbol{u}}_{\textrm{os}} + \hat{{\boldsymbol{u}}}_m$ to write

(71) \begin{eqnarray} &&-\rho \int \limits _{\Omega }({\boldsymbol{u}}_{m}\cdot \nabla ){\boldsymbol{u}}_m\cdot \hat{{\boldsymbol{u}}}_m\,\textrm{d}x=-\rho \int \limits _{\Omega }({{\boldsymbol{u}}}_m\cdot \nabla )\hat{{\boldsymbol{u}}}_m\cdot \hat{{\boldsymbol{u}}}_m\,\textrm{d}x-\rho \int \limits _{\Omega }({{\boldsymbol{u}}}_m\cdot \nabla ){{\boldsymbol{u}}}_{\textrm{os}}\cdot \hat{{\boldsymbol{u}}}_m\,\textrm{d}x.\quad \,\, \end{eqnarray}

Next, using integration by parts, the non-penetration condition $\hat{{\boldsymbol{u}}}_m\cdot \boldsymbol{\nu }=0$ on $\partial \mathcal{P}_{\textrm{st}}$ and the divergence-free condition $\nabla \cdot{\boldsymbol{u}}_m=0$ , we get

(72) \begin{equation} \int \limits _{\Omega }({{\boldsymbol{u}}}_m\cdot \nabla )\hat{{\boldsymbol{u}}}_m\cdot \hat{{\boldsymbol{u}}}_m\,\textrm{d}x=\dfrac{1}{2}\int \limits _{\Omega }{{\boldsymbol{u}}}_m\cdot \nabla |\hat{{\boldsymbol{u}}}_m|^2\,\textrm{d}x = 0. \end{equation}

To estimate the second integral in the right-hand side of (71), we integrate by parts and use the non-penetration condition again on $\partial \mathcal{P}_{\textrm{st}}$ to get:

(73) \begin{eqnarray} -\rho \int \limits _{\Omega } ({\boldsymbol{u}}_m\cdot \nabla ){\boldsymbol{u}}_{\textrm{os}}\cdot \hat{{\boldsymbol{u}}}_m\,\textrm{d}x&=&\rho \int \limits _{\Omega }({\boldsymbol{u}}_m \cdot \nabla ) \hat{{\boldsymbol{u}}}_m\cdot{\boldsymbol{u}}_{\textrm{os}}\,\textrm{d}x. \end{eqnarray}

Finally, we use the Poincaré estimate for $\hat{{\boldsymbol{u}}}_m$ (one can also use (199) with $\hat{{\boldsymbol{u}}}_m$ instead of $\boldsymbol{Q}$ ) as well as that the offset function ${\boldsymbol{u}}_{\textrm{os}}$ is a smooth function with bounded derivatives:

(74) \begin{eqnarray} &&-\rho \int \limits _{\Omega }({\boldsymbol{u}}_{m}\cdot \nabla ){\boldsymbol{u}}_m\cdot \hat{{\boldsymbol{u}}}_m\,\textrm{d}x=\rho \int \limits _{\Omega }(\hat{{\boldsymbol{u}}}_m \cdot \nabla ) \hat{{\boldsymbol{u}}}_m\cdot{\boldsymbol{u}}_{\textrm{os}}\,\textrm{d}x+\rho \int \limits _{\Omega }({\boldsymbol{u}}_{\textrm{os}} \cdot \nabla ) \hat{{\boldsymbol{u}}}_m\cdot{\boldsymbol{u}}_{\textrm{os}}\,\textrm{d}x\nonumber \\[5pt] &&\hspace{40 pt}\leq \rho \int \limits _{\Omega }|\nabla \hat{{\boldsymbol{u}}}_{m}|(|\hat{{\boldsymbol{u}}}_m||{\boldsymbol{u}}_{\textrm{os}}|+|{\boldsymbol{u}}_{\textrm{os}}|^2)\,\textrm{d}x\nonumber \\[5pt] &&\hspace{40 pt}\leq C\rho \|{{\boldsymbol{u}}}_{\textrm{os}}\|_{L^\infty (\Omega )} \!\left (\int \limits _{\Omega }|\nabla \hat{{\boldsymbol{u}}}_m|^2\,\textrm{d}x\right )+C^{*}. \end{eqnarray}

2. Here, we bound the second integral in the right-hand side of (70) by the Cauchy–Schwarz inequality:

(75) \begin{eqnarray} -\eta \int \limits _{\Omega }\nabla{\boldsymbol{u}}_{\textrm{os}}\;:\;\nabla \hat{{\boldsymbol{u}}}_m\,\textrm{d}x\leq \dfrac{\eta }{4}\int \limits _{\Omega }|\nabla \hat{{\boldsymbol{u}}}_m|^2\,\textrm{d}x+C^{*}. \end{eqnarray}

3. We have the equality $ \sigma _a (A,B)\;:\;D+S(D,A)\;:\;B=0$ which holds for all matrices $A$ , $B$ and $D$ such that $A$ and $B$ are symmetric:

(76) \begin{equation} \int \limits _{\Omega }\!\left (\sigma _a({\boldsymbol{Q}}_m,{\bf H}_m)\;:\; \nabla \hat{{\boldsymbol{u}}}_m+S(\nabla \hat{{\boldsymbol{u}}}_m,{\boldsymbol{Q}}_m)\;:\;{\bf H}_m\right )\,\textrm{d}x=0. \end{equation}

4. Note that the integral in the 4th line of (70) vanishes. Indeed, using integration by parts, $\nabla \cdot \hat{{\boldsymbol{u}}}_m=0$ and (59), we get

(77) \begin{eqnarray} &&\int \limits _{\Omega }\nabla{\boldsymbol{Q}}_m\odot \nabla{\boldsymbol{Q}}_m\;:\;\nabla \hat{{\boldsymbol{u}}}_m \,\textrm{d}x\nonumber \nonumber \\[5pt]&&\hspace{20 pt} = -\dfrac{1}{2}\int \limits _{\Omega }(\hat{{\boldsymbol{u}}}_m\cdot \nabla )|\nabla{\boldsymbol{Q}}_m|^2\,\textrm{d}x-\int \limits _{\Omega }(\hat{{\boldsymbol{u}}}_m\cdot \nabla ){\boldsymbol{Q}}_m\;:\;\Delta{\boldsymbol{Q}}_m \,\textrm{d}x +\int \limits _{\partial \mathcal{P}_{\textrm{st}} }\!\left [\!\left (\nabla{\boldsymbol{Q}}_m\odot \nabla{\boldsymbol{Q}}_m\right ) \nu \right ]\cdot \hat{{\boldsymbol{u}}}_m\,\textrm{d}S_x\nonumber \\[5pt] &&\hspace{30 pt}=-\dfrac{1}{2}\int \limits _{\partial \mathcal{P}_{\textrm{st}}}(\hat{{\boldsymbol{u}}}_m\cdot \nu )|\nabla{\boldsymbol{Q}}_m|^2\,\textrm{d}S_x+\dfrac{1}{2}\int \limits _{\Omega }(\nabla \cdot \hat{{\boldsymbol{u}}}_m)|\nabla{\boldsymbol{Q}}_m|^2\,\textrm{d}x-\int \limits _{\Omega }(\hat{{\boldsymbol{u}}}_m\cdot \nabla ){\boldsymbol{Q}}_m\;:\;\Delta{\boldsymbol{Q}}_m \,\textrm{d}x\nonumber \\[5pt] &&\hspace{30 pt}=-\int \limits _{\Omega }(\hat{{\boldsymbol{u}}}_m\cdot \nabla ){\boldsymbol{Q}}_m\;:\;\Delta{\boldsymbol{Q}}_m \,\textrm{d}x. \end{eqnarray}

5. We use the Cauchy–Schwarz inequality, (59), and the a priori bound (64) to estimate the 5th line of (70):

(78) \begin{eqnarray} &&K\int \limits _{\Omega } ({\boldsymbol{u}}_{\textrm{os}}\cdot \nabla ){\boldsymbol{Q}}_{m}\;:\;\Delta{{\boldsymbol{Q}}}_m \,\textrm{d}x +\int \limits _{\Omega }({\boldsymbol{u}}_m\cdot \nabla ){\boldsymbol{Q}}_m\;:\;\hat{{\bf H}}_M({\boldsymbol{Q}}_m)\,\textrm{d}x\nonumber \\[5pt] &&\hspace{50 pt}=\int \limits _{\Omega } ({\boldsymbol{u}}_{\textrm{os}}\cdot \nabla ){\boldsymbol{Q}}_{m}\;:\;{\bf H}_m \,\textrm{d}x+\int \limits _{\Omega }(\hat{{\boldsymbol{u}}}_m\cdot \nabla )\hat{\mathcal{F}}_M({\boldsymbol{Q}}_m)\,\textrm{d}x\nonumber \\[5pt] &&\hspace{55 pt}\leq C\|{\boldsymbol{u}}_{\textrm{os}}\|_{L^{\infty }(\Omega )}(\sqrt{K}\|\nabla{\boldsymbol{Q}}_m\|_{L^2(\Omega )}^2+\dfrac{1}{\sqrt{K}}\|{\bf H}_m\|^2_{L^2(\Omega )})+\int \limits _{\partial{\mathcal{P}}_{\textrm{st}}}(\hat{{\boldsymbol{u}}}_m\cdot \nu )\hat{\mathcal{F}}_M({\boldsymbol{Q}}_m)\,\textrm{d}S_x \nonumber \\[5pt] && \hspace{55 pt} \leq C\|{\boldsymbol{u}}_{\textrm{os}}\|_{L^{\infty }(\Omega )}(\sqrt{K}\|\nabla{\boldsymbol{Q}}_m\|_{L^2(\Omega )}^2+\dfrac{1}{\sqrt{K}}\|{\bf H}_m\|^2_{L^2(\Omega )})\nonumber \\[5pt] && \hspace{55 pt}\leq \dfrac{C}{\sqrt{K}}\|{\boldsymbol{u}}_{\textrm{os}}\|_{L^{\infty }(\Omega )}\|{\bf H}_m\|_{L^2(\Omega )}^2+C^{*}. \end{eqnarray}

6. We use again the a priori bound (64) and the Cauchy–Schwarz inequality to estimate the first term in the 6th line of (70):

(79) \begin{eqnarray} \int _{\Omega } S(\nabla{\boldsymbol{u}}_{\textrm{os}},{\boldsymbol{Q}}_m)\;:\;{\bf H}_m \textrm{d}x &\lt & C\|{\boldsymbol{u}}_{\textrm{os}}\|_{W^{1,\infty }(\Omega )}(\|{\boldsymbol{Q}}_m\|_{L^2(\Omega )}^2 + \|{\bf H}_m\|_{L^{2}(\Omega )}^2) \nonumber \\[5pt] &\lt & C \|{\boldsymbol{u}}_{\textrm{os}}\|_{W^{1,\infty }(\Omega )} \|{\bf H}_m\|_{L^2(\Omega )}^2 + C^{*}. \end{eqnarray}

7. Finally, the last term in (70) is estimated as follows:

(80) \begin{eqnarray} \int _{\Omega } F_{\textrm{ext},M}\;:\;{\bf H}_{m} \textrm{d}x \lt \dfrac{\Gamma }{4}\|{\bf H}_m\|_{L^{2}(\Omega )}^2+C^{*}. \end{eqnarray}

Collect (74)–(80) and substitute them in (70):

\begin{eqnarray*} &&\!\left (\dfrac{3\eta }{4}-C\rho \|{\boldsymbol{u}}_{\textrm{os}}\|_{L^{\infty }(\Omega )}\right )\|\hat{{\boldsymbol{u}}}_{m}\|_{L^2(\Omega )}^2\\[5pt] &&\hspace{60 pt}+ \left (\dfrac{3\Gamma }{4}-\dfrac{C}{\sqrt{K}}\|{\boldsymbol{u}}_{\textrm{os}}\|_{L^{\infty }(\Omega )}-C\|{\boldsymbol{u}}_{\textrm{os}}\|_{W^{1,\infty }(\Omega )}\right )\|{\bf H}_m\|_{L^2(\Omega )}^2\leq C^{*}. \end{eqnarray*}

Under the restrictions (42), the inequality (69) holds proving the lemma.

4.4. Existence of Galerkin approximations

We will use the following result [58, Lemma IX.3.1, p. 597]:

Theorem 4.1. Let ${\mathcal{P}}\;:\;\mathbb R^{p}\to \mathbb R^{p}$ be a continuous mapping such that for some $R\gt 0$ :

(81) \begin{equation}{\mathcal{P}}(\boldsymbol{\xi })\cdot \boldsymbol{\xi } \geq 0 \quad \textrm{ for all }\boldsymbol{\xi } \in \mathbb R^p \textrm{ with }|\boldsymbol{\xi }|=R. \end{equation}

Then there exists $\boldsymbol{\xi }_0\in \mathbb R^p$ with $|\boldsymbol{\xi }_0|\leq R$ such that ${\mathcal{P}}(\boldsymbol{\xi }_0)=0$ .

Next, we introduce the mapping $\mathcal{P}$ for problem (66)–(67) with unknowns $\hat{{\boldsymbol{u}}}_m=\sum \limits _{k=1}^{m}u_{km}\Psi _{k}$ and ${\bf H}_{m}=\sum \limits _{k=1}^{m}h_{km}\Phi _{k}$ , defined in (53). Given $m\geq 1$ , let

\begin{equation*}\boldsymbol {\xi }=(u_{1m},\ldots,u_{mm},h_{1m},\ldots,h_{mm})\in \mathbb R^{2m}\end{equation*}

and the $k$ th component mapping ${\mathcal{P}}\;:\;\mathbb R^{2m}\to \mathbb R^{2m}$ $(p=2m)$ is the left-hand side of (66) for $1\leq k\leq m$ and the left-hand side of (67) for $m+1\leq k\leq 2m$ . We obtain that

\begin{equation*} {\mathcal {P}}(\boldsymbol {\xi })\cdot \boldsymbol {\xi }=\eta \|\nabla \hat {{\boldsymbol{u}}}_m\|_{L^2(\Omega )}^2+\Gamma \|{\bf H}_m\|_{L^2(\Omega )}^2-\mathcal {R}({\boldsymbol{u}}_m,{\bf H}_m), \end{equation*}

where $\mathcal{R}({\boldsymbol{u}}_m,{\bf H}_m)$ is the right-hand side of (70). In the proof of Lemma 4.3, we showed that

\begin{eqnarray*} &&|\mathcal{R}({\boldsymbol{u}}_m,{\bf H}_m)|\leq \!\left (\dfrac{\eta }{4}+C\rho \|{\boldsymbol{u}}_{\textrm{os}}\|_{L^{\infty }(\Omega )}\right )\|\hat{{\boldsymbol{u}}}_{m}\|_{L^2(\Omega )}^2\\[5pt] &&\hspace{80 pt}+ \left (\dfrac{\Gamma }{4}+\dfrac{C}{\sqrt{K}}\|{\boldsymbol{u}}_{\textrm{os}}\|_{L^{\infty }(\Omega )}+C\|{\boldsymbol{u}}_{\textrm{os}}\|_{W^{1,\infty }(\Omega )}\right )\|{\bf H}_m\|_{L^2(\Omega )}^2+ C^{*}. \end{eqnarray*}

Therefore, using this inequality we obtain the following:

Lemma 4.4. Assume $ \eta \gt 2C_0\rho \|u_{\textrm{sq}}\|_{L^{\infty }(\Omega )}\textrm{ and }\Gamma \gt 2C_0\!\left (\dfrac{1}{\sqrt{K}}\|u_{\textrm{sq}}\|_{L^{\infty }(\Omega )}+\|u_{\textrm{sq}}\|_{W^{1,\infty }(\Omega )}\right ),$ with $C_0$ from Lemma 4.3. Then there exists constants $C_1,C_2\gt 0$ independent of $m$ such that

\begin{equation*} {\mathcal {P}}(\boldsymbol {\xi })\cdot \boldsymbol {\xi }\geq C_1 \!\left (\|\nabla \hat {{\boldsymbol{u}}}_m\|_{L^2(\Omega )}^2+\|{\bf H}_m\|_{L^2(\Omega )}^2\right )-C_2. \end{equation*}

The condition (81) is satisfied for large $R\gt 0$ and thus we have the following existence result for our Galerkin approximations:

Theorem 4.2. Assume $ \eta \gt 2C_0\rho \|u_{\textrm{sq}}\|_{L^{\infty }(\Omega )}\textrm{ and }\Gamma \gt 2C_0\!\left (\dfrac{1}{\sqrt{K}}\|u_{\textrm{sq}}\|_{L^{\infty }(\Omega )}+\|u_{\textrm{sq}}\|_{W^{1,\infty }(\Omega )}\right ),$ with $C_0$ from Lemma 4.3 . Then there exists a solution $(\hat{{\boldsymbol{u}}}_m,{\bf H}_m)$ of (66)–(67). Moreover, if ${\boldsymbol{Q}}_m$ is defined via (63), then the solution satisfies

(82) \begin{equation} \|\hat{{\boldsymbol{u}}}_m\|_{H^1(\Omega )}^2+\|{\boldsymbol{Q}}_{m}\|_{H^2(\Omega )}^2+\hat{\mathcal{F}}_M({\boldsymbol{Q}}_m)+\|{\boldsymbol{Q}}_{\textrm{pref}}-{\boldsymbol{Q}}_m\|^2_{L^2(\partial{\mathcal{P}}_{\textrm{st}})}+\|{\bf H}_m\|_{L^2(\Omega )}^2 \lt C. \end{equation}

4.5. Passing to limit $m\to \infty$

From (82), we get that there is a sub-sequence of $\!\left \{(\hat{{\boldsymbol{u}}}_m,{\bf H}_m)\right \}$ such that

(83) \begin{align} \hat{{\boldsymbol{u}}}_m \rightharpoonup \hat{{\boldsymbol{u}}}\textrm{ in }H^1(\Omega ) \end{align}

(84) \begin{align} {\bf H}_m \rightharpoonup{\bf H}\textrm{ in }L^2(\Omega ) \end{align}

(85) \begin{align} {\boldsymbol{Q}}_m \rightharpoonup{\boldsymbol{Q}} \textrm{ in }H^2(\Omega ). \end{align}

Next, we will use the following auxiliary lemma [59, Lemma 1.3]:

Lemma 4.5. Let $\mathcal{O}$ be a bounded domain. Let $p_m(x)$ and $p(x)$ be such functions from $L^q(\mathcal{O})$ , $1\lt q\lt \infty$ , such that

(86) \begin{equation} \|p_m\|_{L^q(\mathcal{O})}\leq C\textrm{ and }p_m\to p \textrm{ a.e. in }\mathcal{O}. \end{equation}

Then $p_m\rightharpoonup p$ in $L^q(\mathcal{O})$ .

From (82) and Lemma 4.5, we get

(87) \begin{equation} \hat{{\bf H}}_M({\boldsymbol{Q}}_m) \rightharpoonup \hat{{\bf H}}_M({\boldsymbol{Q}})\textrm{ in } L^2(\Omega ). \end{equation}

Using (84), (85) and (87) as well as the trace theorem, we can pass to the limits $m\to \infty$ in the weak formulation of (63):

(88) \begin{equation} K\int \limits _{\Omega }\nabla{\boldsymbol{Q}} \cdot \nabla G\,\textrm{d}x+W\int \limits _{\partial{\mathcal{P}}_{\textrm{st}}}({\boldsymbol{Q}}_{\textrm{pref}}-{\boldsymbol{Q}})\;:\;G\,\textrm{d}S_x+\int \limits _{\Omega }\hat{{\bf H}}_M({\boldsymbol{Q}})\;:\;G\,\textrm{d}x=\int \limits _{\Omega }{\bf H}\;:\;G \,\textrm{d}x \end{equation}

for all smooth test functions $G$ .

Next, we pass to the limit in (66)–(67) using (83), (85), $H^2(\Omega )\hookrightarrow H^1(\Omega )\hookrightarrow L^2(\Omega )$ , and the property that product of strongly and weakly converging sequences weakly converges to the product of corresponding limits. We get ( ${\boldsymbol{u}}={\boldsymbol{u}}_{\textrm{os}} +\hat{{\boldsymbol{u}}}$ ):

(89) \begin{eqnarray} &&\eta \int \limits _{\Omega }\nabla \hat{{\boldsymbol{u}}}\;:\;\nabla \Psi _k\,\textrm{d}x +\eta \int \limits _{\Omega }\nabla{\boldsymbol{u}}_{\textrm{os}}\;:\;\nabla \Psi _k\,\textrm{d}x+ \rho \int \limits _{\Omega }({\boldsymbol{u}}\cdot \nabla ){\boldsymbol{u}}\cdot \Psi _k\,\textrm{d}x\nonumber \\[5pt]&& \hspace{120 pt}+\int \limits _{\Omega } \sigma _a({\boldsymbol{Q}},{\bf H})\;:\; \nabla \Psi _{k}\,\textrm{d}x=K\int \limits _{\Omega } \nabla{\boldsymbol{Q}}\odot \nabla{\boldsymbol{Q}}\;:\; \nabla \Psi _{k}\,\textrm{d}x. \end{eqnarray}

(90) \begin{eqnarray} &&\Gamma \int \limits _{\Omega }{\bf H}\;:\;\Phi _k \,\textrm{d}x - \int \limits _{\Omega }\!\left (({\boldsymbol{u}}\cdot \nabla ){\boldsymbol{Q}}-S(\nabla{\boldsymbol{u}},{\boldsymbol{Q}})\right )\;:\;\Phi _k\,\textrm{d}x+ \int _{\Omega } F_{\textrm{ext},M} \;:\; \Phi _{k} \,\textrm{d}x=0. \end{eqnarray}

Finally, we can drop subscript $M$ in $\hat{{\bf H}}_M({\boldsymbol{Q}})$ and $F_{\textrm{ext},M}({\boldsymbol{Q}},{\boldsymbol{Q}}_\infty )$ due to the $L^{\infty }$ -a-priori bound on solution of (23)–(28) in Lemma 4.1.

5. Well-posedness of time-dependent problem $-$ Proof of theorem 3.2

In this section, we prove the local-in-time existence of the unique solution with additional regularity by using Banach fixed point theorem. In Section 5.1, we will write the time-dependent problem in the operator form. In Sections 5.2 and 5.3, we will address the Lipschitz properties of the nonlinear part and the solvability of the linear part of PDE system. In Section 5.3, we will prove the local-in-time existence and uniqueness by Banach fixed point theorem.

5.1. Operators and function spaces

We first define the projection operator $P_{\sigma }\;:\; H^{-1}(\Omega ) \rightarrow H^{-1}_{\sigma }(\Omega )$ onto the space of divergence-free functions so that if we apply $P_{\sigma }$ to (1), the pressure $p$ is eliminated. Specifically, the equation (1) becomes

(91) \begin{eqnarray} \partial _t{\boldsymbol{u}} + P_{\sigma }( \nabla \cdot ({\boldsymbol{u}} \otimes{\boldsymbol{u}})) + \frac{d{\boldsymbol{V}}}{dt} - \rho ^{-1}\eta P_{\sigma }(\Delta{\boldsymbol{u}})= \rho ^{-1} P_{\sigma }(\nabla \cdot \sigma _{\textrm{ela}}({\boldsymbol{Q}})). \end{eqnarray}

Now we consider the problem consisting of (91), (2)–(6) with force and torque balances (14), (43). The tuple of unknowns is $\mathcal{U}=({\boldsymbol{u}},{\boldsymbol{Q}},\boldsymbol{\omega },{\boldsymbol{V}})^{\textrm{T}}$ . We rewrite the problem as

(92) \begin{equation} \mathcal{L}\mathcal{U}=\mathcal{N}(\mathcal{U}), \end{equation}

where we define linear operator $\mathcal{L}$ and nonlinear operator $\mathcal{N}$ as

(93) \begin{align} \mathcal{L} \!\left ( \begin{matrix}{\boldsymbol{u}} \\[5pt]{\boldsymbol{Q}} \\[5pt] \boldsymbol{\omega } \\[5pt]{\boldsymbol{V}} \end{matrix} \right ) = \partial _t \!\left ( \begin{matrix}{\boldsymbol{u}} \\[5pt]{\boldsymbol{Q}} \\[5pt] \boldsymbol{\omega } \\[5pt]{\boldsymbol{V}} \end{matrix} \right ) - \!\left ( \begin{matrix} \rho ^{-1}\eta P_{\sigma } ( \Delta{\boldsymbol{u}}) \\[5pt] \Gamma K\Delta{\boldsymbol{Q}} \\[5pt] 0 \\[5pt] 0 \end{matrix} \right ) \end{align}

and

(94) \begin{align} \mathcal{N} \!\left ( \begin{matrix}{\boldsymbol{u}} \\[5pt]{\boldsymbol{Q}} \\[5pt] \boldsymbol{\omega } \\[5pt]{\boldsymbol{V}} \end{matrix} \right ) = \!\left ( \begin{matrix} \rho ^{-1} P_{\sigma } \!\left (\nabla \cdot (\sigma _{\textrm{ela}}({\boldsymbol{Q}}) - \rho{\boldsymbol{u}} \otimes{\boldsymbol{u}}) \right ) - \dfrac{d{\boldsymbol{V}}}{dt} \\[5pt] -{\boldsymbol{u}} \cdot \nabla{\boldsymbol{Q}} + \Gamma \hat{{\bf H}}({\boldsymbol{Q}}) + S(\nabla{\boldsymbol{u}},{\boldsymbol{Q}}) + F_{\textrm{ext}}({\boldsymbol{Q}},{\boldsymbol{Q}}_{\infty }) \\[5pt] \frac{1}{I} \int _{\partial{\mathcal{P}}}{\boldsymbol{x}} \times \sigma \boldsymbol{\nu } + \boldsymbol{\ell }\,\textrm{d}S_{{\boldsymbol{x}}} \\[5pt] \frac{1}{m} \int _{\partial{\mathcal{P}}} \sigma \boldsymbol{\nu } \, \textrm{d}S_{\boldsymbol{x}} \end{matrix} \right ). \end{align}

To handle the nonlinear and inhomogeneous boundary conditions, we represent unknown functions $\boldsymbol{u}$ and $\boldsymbol{Q}$ as

\begin{equation*} {\boldsymbol{u}} = {\boldsymbol{u}}_h + {\boldsymbol{u}}_{\text {os}}\text { and }{\boldsymbol{Q}} = {\boldsymbol{Q}}_h + {\boldsymbol{Q}}_{\text {os}}. \end{equation*}

The offset function ${\boldsymbol{u}}_{\textrm{os}}$ is given by

(95) \begin{align} -\eta \Delta{\boldsymbol{u}}_{\textrm{os}} + \nabla p_{\textrm{os}} = 0 \textrm{ in } \Omega \end{align}

(96) \begin{align}{\boldsymbol{u}}_{\textrm{os}} = u_{\textrm{sq}}(\boldsymbol{\alpha }(t),{\boldsymbol{x}})\boldsymbol{\tau } + \boldsymbol{\omega }(t) \times{\boldsymbol{x}} \textrm{ on } \partial{\mathcal{P}} \end{align}

(97) \begin{align} {\boldsymbol{u}}_{\textrm{os}} \textrm{ periodic in } \Pi \end{align}

The offset function ${\boldsymbol{Q}}_{\textrm{os}}$ is defined such that

(98) \begin{align} K\partial _\nu{\boldsymbol{Q}}_{\textrm{os}} = W({\boldsymbol{Q}}_{\textrm{pref}}-{\boldsymbol{Q}}_{\textrm{os}}) \textrm{ on }\partial{\mathcal{P}} \end{align}

(99) \begin{align} {\boldsymbol{Q}}_{\textrm{os}} \textrm{ periodic in } \Pi \end{align}

Specifically, we define

(100) \begin{eqnarray}{\boldsymbol{Q}}_{\textrm{os}}({\boldsymbol{x}}) ={\boldsymbol{Q}}_{\textrm{pref}}\!\left ( \dfrac{{\boldsymbol{x}}}{|{\boldsymbol{x}}|}\right ) \psi (|{\boldsymbol{x}}|), \quad{\boldsymbol{x}}\in \Pi \setminus \mathcal{P} \end{eqnarray}

Here, $\psi (|{\boldsymbol{x}}|)\geq 0$ is a smooth function such that $\psi (|{\boldsymbol{x}}|) = 1$ for ${\boldsymbol{x}} \in \!\left (\partial \mathcal{P}+B_{r_{*}}(\textbf{0})\right )\cap \Pi \setminus \mathcal{P}$ with $r_{*}=\textrm{dist}(\partial \Pi, \partial \mathcal{P})/4$ , and $\psi (|{\boldsymbol{x}}|) = 0$ when $|{\boldsymbol{x}}| \gt 2r_{*}$ . Boundary condition (98) is satisfied since $\partial _{\nu }{\boldsymbol{Q}}_{\textrm{os}}=\partial _{|{\boldsymbol{x}}|}{\boldsymbol{Q}}_{\textrm{os}}=\textbf{0}$ and $\!\left .{\boldsymbol{Q}}_{\textrm{os}}\right |_{\partial \mathcal{P}}={\boldsymbol{Q}}_{\textrm{pref}}.$ The offset function ${\boldsymbol{Q}}_{\textrm{os}}$ can be extended periodically so it satisfies (99) since ${\boldsymbol{Q}}_{\textrm{os}}\equiv 0$ on $\partial \Pi$ . We point out that ${\boldsymbol{Q}}_{\textrm{os}}$ is the solution of the Poisson problem with boundary conditions (98)–(99) and the PDE $-\Delta{\boldsymbol{Q}}_{\textrm{os}}=\boldsymbol{f}$ with $\boldsymbol{f}=-\Delta \!\left ({\boldsymbol{Q}}_{\textrm{pref}}\!\left ( \dfrac{{\boldsymbol{x}}}{|{\boldsymbol{x}}|}\right ) \psi (|{\boldsymbol{x}}|)\right )$ . Note that the offset function ${\boldsymbol{u}}_{\textrm{os}}$ depends on unknown orientation angle $\boldsymbol{\alpha }(t)$ and angular velocity $\boldsymbol{\omega }(t)$ , whereas ${\boldsymbol{Q}}_{\textrm{os}}$ does not. Therefore, ${\boldsymbol{u}}_{\textrm{os}}$ changes in time $t$ while ${\boldsymbol{Q}}_{\textrm{os}}$ is independent of time $t$ .

With the above, the functions ${\boldsymbol{u}}_h$ and ${\boldsymbol{Q}}_h$ satisfy homogeneous boundary conditions. Their equations in $\Omega$ are similar to the original (91) and (4). More precisely, these equations with force and torque balances in the form of (92) look as follows:

(101) \begin{align} \mathcal{L} \!\left ( \begin{matrix}{\boldsymbol{u}}_h \\[5pt]{\boldsymbol{Q}}_h \\[5pt] \boldsymbol{\omega } \\[5pt]{\boldsymbol{V}} \end{matrix} \right ) = \mathcal{J} \!\left ( \begin{matrix}{\boldsymbol{u}}_h \\[5pt]{\boldsymbol{Q}}_h \\[5pt] \boldsymbol{\omega } \\[5pt]{\boldsymbol{V}} \end{matrix} \right ) \;=\!:\; \mathcal{N} \!\left ( \begin{matrix}{\boldsymbol{u}}_h +{\boldsymbol{u}}_{\textrm{os}} \\[5pt]{\boldsymbol{Q}}_h +{\boldsymbol{Q}}_{\textrm{os}} \\[5pt] \boldsymbol{\omega } \\[5pt]{\boldsymbol{V}} \end{matrix} \right ) - \mathcal{L} \!\left ( \begin{matrix}{\boldsymbol{u}}_{\textrm{os}} \\[5pt]{\boldsymbol{Q}}_{\textrm{os}} \\[5pt] 0 \\[5pt] 0 \end{matrix} \right ) \end{align}

To describe the domains of the operator $\mathcal{L}$ , we introduce the following Banach spaces:

(102) \begin{align} X_{{\boldsymbol{u}}} &= \!\left \{{\boldsymbol{u}} \in H^2(0, T;\;L_{\sigma }^{2}(\Omega )) \cap H^{1}(0, T;\;H_{0, \sigma }^{2}(\Omega )) \!\left |\begin{array}{l}{\boldsymbol{u}} = 0 \textrm{ on } \partial{\mathcal{P}},\\[5pt] {\boldsymbol{u}} \textrm{ periodic in } \Pi \end{array}\right. \right \} \end{align}

(103) \begin{align} X_{\boldsymbol{Q}} &= \!\left \{{\boldsymbol{Q}} \in H^{2}(0, T;\;H^1(\Omega )) \cap H^{1}(0, T;\;H^{3}(\Omega )) \!\left | \begin{array}{l}\partial _{\nu }{\boldsymbol{Q}} = -W{\boldsymbol{Q}} \textrm{ on } \partial{\mathcal{P}},\\[5pt] {\boldsymbol{Q}} \textrm{ periodic in } \Pi \end{array}\right. \right \} \end{align}

with corresponding norms

\begin{eqnarray*} &&\|{\boldsymbol{u}}\|_{X_{\boldsymbol{u}}} = \!\left (\|{\boldsymbol{u}}\|_{H^2(0, T;L_{\sigma }^{2}(\Omega ))}^2 + \|{\boldsymbol{u}}\|_{H^{1}(0, T;H_{0, \sigma }^{2}(\Omega ))}^2 + \|{\boldsymbol{u}}\big |_{t = 0}\|_{H^{2}_{0, \sigma }(\Omega )}^2 + \|\partial _t{\boldsymbol{u}}\big |_{t = 0}\|_{H^{1}_{0, \sigma }(\Omega )}^2\right )^{\frac{1}{2}}\\[5pt] &&\|{\boldsymbol{Q}}\|_{X_{\boldsymbol{Q}}} = \!\left (\|{\boldsymbol{Q}}\|_{H^2(0, T;H^{1}(\Omega ))}^2 + \|{\boldsymbol{Q}}\|_{H^{1}(0, T;H^{3}(\Omega ))}^2 + \|{\boldsymbol{Q}}\big |_{t = 0}\|_{H^{3}(\Omega )}^2 + \|\partial _t{\boldsymbol{Q}}\big |_{t = 0}\|_{H^{2}(\Omega )}^2\right )^{\frac{1}{2}}. \end{eqnarray*}

Introduce also

\begin{align*} Y_{{\boldsymbol{u}}} = H^{1}(0, T;\;L^{2}_{\sigma }(\Omega )),\quad Y_{\boldsymbol{Q}} = H^{1}(0, T;\;H^{1}(\Omega )). \end{align*}

Then $X = X_{{\boldsymbol{u}}} \times X_{\boldsymbol{Q}} \times H^{2}(0, T) \times H^{2}(0, T)$ and $Y = Y_{{\boldsymbol{u}}} \times Y_{\boldsymbol{Q}} \times H^{1}(0, T) \times H^{1}(0, T)$ are the domain and the range of the operator $\mathcal{L}$ . The corresponding norms are

\begin{eqnarray*} &&\|({\boldsymbol{u}},{\boldsymbol{Q}},{\boldsymbol{\omega }},{\boldsymbol{V}})\|_X = (\|{\boldsymbol{u}}\|_{X_{\boldsymbol{u}}}^2 + \|{\boldsymbol{Q}}\|_{X_{\boldsymbol{Q}}}^2 + \|{\boldsymbol{\omega }}\|_{H^{2}(0, T)}^2 + \|{\boldsymbol{V}}\|_{H^{2}(0, T)}^2)^\frac{1}{2}, \\[5pt] &&\|({\boldsymbol{u}},{\boldsymbol{Q}},{\boldsymbol{\omega }},{\boldsymbol{V}})\|_Y = (\|{\boldsymbol{u}}\|_{Y_{\boldsymbol{u}}}^2 + \|{\boldsymbol{Q}}\|_{Y_{\boldsymbol{Q}}}^2 + \|{\boldsymbol{\omega }}\|_{H^{1}(0, T)}^2 + \|{\boldsymbol{V}}\|_{H^{1}(0, T)}^2)^\frac{1}{2}. \end{eqnarray*}

5.2. Lipschitz property of the nonlinear part

In this section, we show the Lipschitz property of the nonlinear operator $\mathcal{J}$ with respect to the norms of the spaces $X$ and $Y$ . Below, we will use short notations for spaces of functions depending on both $t$ and $\boldsymbol{x}$ . Namely, let $V$ and $W$ stand by either of $L^{\infty }$ , $L^2$ , $L_\sigma ^2$ , $W^{1,\infty }$ and $H^{s}$ for $s\geq 0$ . We will denote $V(0,T;\;W(\Omega ))$ by $V_tW_{{\boldsymbol{x}}}$ and the corresponding norm by $\|\cdot \|_{V_tW_{{\boldsymbol{x}}}}$ . For example, $H^{2}(0, T;\;L^{2}(\Omega ))$ will be denoted by $H^2_tL^2_{{\boldsymbol{x}}}$ and the norm by $\|\cdot \|_{H^2_tL^{2}_x}$ . We start with the following estimates with constants vanishing as $T\to 0$ .

Proposition 5.1. There exists a constant $C(T)$ such that $C(T)\to 0$ as $T\to 0$ and each of inequalities (i)-(iv) below holds for all $f$ and $g$ as long as the left-hand side of the inequality is finite:

(104) \begin{align} (\boldsymbol{i})\qquad \!\left \|f\right \|_{L^{\infty }_tL^{\infty }_{{\boldsymbol{x}}}}\leq C(T)\|f\|_{H^1_tH^2_{{\boldsymbol{x}}}}, \end{align}

(105) \begin{align} (\boldsymbol{ii})\qquad \!\left \| f \right \|_{H^{1}_tL^{\infty }_{{\boldsymbol{x}}}}\leq C(T) (\!\left \| f \right \|_{H^{1}_tH^{2}_{{\boldsymbol{x}}}} + \!\left \| f \right \|_{H^{2}_tL^{2}_{{\boldsymbol{x}}}}), \end{align}

(106) \begin{align} (\boldsymbol{iii})\qquad \|fg\|_{H^1_tH^1_{{\boldsymbol{x}}}}\leq C(T) (\!\left \| f \right \|_{H^{1}_tH^{2}_{{\boldsymbol{x}}}} + \!\left \| f \right \|_{H^{2}_tL^{2}_{{\boldsymbol{x}}}}) \cdot (\!\left \| g \right \|_{H^{1}_tH^{2}_{{\boldsymbol{x}}}} + \!\left \| g \right \|_{H^{2}_tL^{2}_{{\boldsymbol{x}}}}), \end{align}

(107) \begin{align} (\boldsymbol{iv})\qquad \|fg\|_{H^1_tH^1_{{\boldsymbol{x}}}} \leq C(T) (\!\left \| f \right \|_{H^{1}_tH^{3}_{{\boldsymbol{x}}}} + \!\left \| f \right \|_{H^{2}_tH^{1}_{{\boldsymbol{x}}}}) \cdot \left \| g \right \|_{H^{1}_tH^{1}_{{\boldsymbol{x}}}}. \end{align}

Proof. In the proof, below $C$ is independent from $T$ unless the dependence is indicated via the following notation $C(T)$ . All constants $C(T)$ vanish as $T\to 0$ . We will use the following inequalities in the proof:

(108) \begin{align} \!\left \|f\right \|_{L^{\infty }(\Omega )}\leq C\|f\|_{H^2(\Omega )} \textrm{ (General Sobolev Inequality, [60])} \end{align}

(109) \begin{align} \|f\|_{L^{\infty }(0,T)}\leq C\|f\|^{1/2}_{L^2(0,T)}\|f\|^{1/2}_{H^1(0,T)} \textrm{ (Agmon's inequality in 1D)} \end{align}

(110) \begin{align} \|f\|_{L^{\infty }(\Omega )}\leq C\|f\|^{1/4}_{L^2(\Omega )}\|f\|^{3/4}_{H^2(\Omega )} \textrm{ (Agmon's inequality in 2D and 3D)} \end{align}

(111) \begin{align} \|f\|_{L^{2}(0,T)}\leq T^{1/2} \|f\|_{L^{\infty }(0,T)} \end{align}

(112) \begin{align} \|f\|_{L^{\infty }(0,T)}\leq C \|f\|_{H^1(0,T)} \end{align}

(113) \begin{align} \|fg\|_{H^1(\Omega )}\leq C \!\left (\|f\|_{H^1(\Omega )}\|g\|_{L^\infty (\Omega )} + \|f\|_{L^\infty (\Omega )}\|g\|_{H^1(\Omega )}\right ) \end{align}

(114) \begin{align} \|fg\|_{H^1(\Omega )}\leq C\!\left ( \|f\|_{L^2(\Omega )}\|g\|_{W^{1,\infty }(\Omega )} + \|f\|_{W^{1,\infty }(\Omega )}\|g\|_{L^2(\Omega )}\right ) \end{align}

Proof of (i):

(115) \begin{eqnarray} \!\left \| f \right \|_{L^{\infty }_tL^{\infty }_{{\boldsymbol{x}}}} &&\leq \!\left \| f \right \|_{L^{\infty }_tH^{2}_{{\boldsymbol{x}}}} \quad \textrm{(use (108)}\nonumber \\[5pt] &&\leq \!\left \| f \right \|_{L^{2}_tH^{2}_{{\boldsymbol{x}}}}^{\frac{1}{2}} \!\left \| f \right \|_{H^{1}_tH^{2}_{{\boldsymbol{x}}}}^{\frac{1}{2}} \, (\textrm{use (109)})\nonumber \\[5pt] &&\leq CT^{1/4}\!\left \| f \right \|_{L^{\infty }_tH^{2}_{{\boldsymbol{x}}}}^{\frac{1}{2}} \!\left \| f \right \|_{H^{1}_tH^{2}_{{\boldsymbol{x}}}}^{\frac{1}{2}} \quad (\textrm{use (111)}) \nonumber \\[5pt]&& \leq CT^{1/4}\!\left \| f \right \|_{H^{1}_tH^{2}_{{\boldsymbol{x}}}} \quad (\textrm{use (112)}). \end{eqnarray}

Proof of (ii):

(116) \begin{align} \!\left \| f \right \|_{H^{1}_tL^{\infty }_{{\boldsymbol{x}}}} & \leq C\!\left (\!\left \| f \right \|_{L^{2}_tL^{\infty }_{{\boldsymbol{x}}}} + \!\left \| \partial _t f \right \|_{L^{2}_tL^{\infty }_{{\boldsymbol{x}}}}\right ) \nonumber \\[5pt] &\leq C(T^{1/2} \!\left \| f \right \|_{L^{\infty }_tH^{2}_{{\boldsymbol{x}}}} + \!\left \| \partial _t f \right \|_{L^{2}_tL^{2}_{{\boldsymbol{x}}}}^{1/4} \cdot \left \| \partial _t f \right \|_{L^{2}_tH^{2}_{{\boldsymbol{x}}}}^{3/4}) \textrm{ (use (108), (111), and (110))}\nonumber \\[5pt] & \leq C(T^{1/2} \!\left \| f \right \|_{L^{\infty }_tH^{2}_{{\boldsymbol{x}}}} + T^{1/8} \!\left \| \partial _t f \right \|_{L^{\infty }_tL^{2}_{{\boldsymbol{x}}}}^{1/4} \cdot \left \| \partial _t f \right \|_{L^{2}_tH^{2}_{{\boldsymbol{x}}}}^{3/4}) \textrm{ (use (111))}\nonumber \\[5pt] & \leq C(T^{1/2} \!\left \| f \right \|_{H^{1}_tH^{2}_{{\boldsymbol{x}}}} + T^{1/8} \!\left \| \partial _t f \right \|_{H^{1}_tL^{2}_{{\boldsymbol{x}}}}^{1/4} \cdot \left \| \partial _t f \right \|_{L^{2}_tH^{2}_{{\boldsymbol{x}}}}^{3/4}) \textrm{ (use (112))} \nonumber \\[5pt] & \leq C(T^{1/2} \!\left \| f \right \|_{H^{1}_tH^{2}_{{\boldsymbol{x}}}} + T^{1/8} \!\left \| f \right \|_{H^{2}_tL^{2}_{{\boldsymbol{x}}}}^{1/4} \cdot \left \| f \right \|_{H^{1}_tH^{2}_{{\boldsymbol{x}}}}^{3/4}) \nonumber \\[5pt] & \leq C(T) (\!\left \| f \right \|_{H^{1}_tH^{2}_{{\boldsymbol{x}}}} + \!\left \| f \right \|_{H^{2}_tL^{2}_{{\boldsymbol{x}}}}). \end{align}

Proof of (iii):

(117) \begin{align} \|fg\|_{H^1_tH^1_{{\boldsymbol{x}}}} & \leq \|f\|_{H^1_tL^{\infty }_{{\boldsymbol{x}}}}\|g\|_{L^{\infty }_tH^{1}_{{\boldsymbol{x}}}} + \|g\|_{H^1_tL^{\infty }_{{\boldsymbol{x}}}}\|f\|_{L^{\infty }_tH^{1}_{{\boldsymbol{x}}}} \nonumber \\[5pt] & \hspace{2cm} + \|f\|_{H^1_tH^1_{{\boldsymbol{x}}}}\|g\|_{L^{\infty }_tL^{\infty }_{{\boldsymbol{x}}}} + \|g\|_{H^1_tH^1_{{\boldsymbol{x}}}}\|f\|_{L^{\infty }_tL^{\infty }_{{\boldsymbol{x}}}} \textrm{ (use (113))} \end{align}

Next, estimate each term in the right-hand side of (117):

(118) \begin{align} & \|f\|_{H^1_tL^{\infty }_{{\boldsymbol{x}}}}\|g\|_{L^{\infty }_tH^{1}_{{\boldsymbol{x}}}} \leq C(T)(\!\left \| f \right \|_{H^{1}_tH^{2}_{{\boldsymbol{x}}}} + \!\left \| f \right \|_{H^{2}_tL^{2}_{{\boldsymbol{x}}}}) \cdot \|g\|_{H^{1}_tH^{1}_{{\boldsymbol{x}}}} \textrm{ (use (116) and (112))} \end{align}

(119) \begin{align} & \|g\|_{H^1_tL^{\infty }_{{\boldsymbol{x}}}}\|f\|_{L^{\infty }_tH^{1}_{{\boldsymbol{x}}}} \leq C(T)(\!\left \| g \right \|_{H^{1}_tH^{2}_{{\boldsymbol{x}}}} + \!\left \| g \right \|_{H^{2}_tL^{2}_{{\boldsymbol{x}}}}) \cdot \|f\|_{H^{1}_tH^{1}_{{\boldsymbol{x}}}} \textrm{ (use (116) and (112))} \end{align}

(120) \begin{align} & \|f\|_{H^1_tH^1_{{\boldsymbol{x}}}}\|g\|_{L^{\infty }_tL^{\infty }_{{\boldsymbol{x}}}} \leq C(T)\|f\|_{H^1_tH^1_{{\boldsymbol{x}}}}\|g\|_{H^1_tH^2_{{\boldsymbol{x}}}} \textrm{ (use (115))} \end{align}

(121) \begin{align} & \|g\|_{H^1_tH^1_{{\boldsymbol{x}}}}\|f\|_{L^{\infty }_tL^{\infty }_{{\boldsymbol{x}}}} \leq C(T)\|g\|_{H^1_tH^1_{{\boldsymbol{x}}}}\|f\|_{H^1_tH^2_{{\boldsymbol{x}}}} \textrm{ (use (115))}. \end{align}

Combining (117)–(121), we obtain (106).

Proof of (iv):

(122) \begin{align} \|fg\|_{H^1_tH^1_{{\boldsymbol{x}}}} & \leq \|f\|_{H^1_{t}W^{1, \infty }_{{\boldsymbol{x}}}}\|g\|_{L^{\infty }_tL^{2}_{{\boldsymbol{x}}}} + \|f\|_{H^1_{t}L^{\infty }_{{\boldsymbol{x}}}}\|g\|_{L^{\infty }_tH^{1}_{{\boldsymbol{x}}}} \nonumber \\[5pt] & \hspace{2cm} + \|f\|_{L^{\infty }_tW^{1, \infty }_{{\boldsymbol{x}}}}\|g\|_{H^{1}_tL^{2}_{{\boldsymbol{x}}}} + \|f\|_{L^{\infty }_tL^{\infty }_{{\boldsymbol{x}}}}\|g\|_{H^{1}_tH^{1}_{{\boldsymbol{x}}}} \textrm{ (use (114)).} \end{align}

Next, estimate each term in the right-hand side of (122):

(123) \begin{align} \|f\|_{H^1_tW^{1, \infty }_{{\boldsymbol{x}}}}\|g\|_{L^{\infty }_tL^{2}_{{\boldsymbol{x}}}} &\leq C(T) (\!\left \| f \right \|_{H^{1}_tH^{3}_{{\boldsymbol{x}}}} + \!\left \| f \right \|_{H^{2}_tH^{1}_{{\boldsymbol{x}}}}) \cdot \|g\|_{H^{1}_{t}L^{2}_{{\boldsymbol{x}}}} \textrm{ (use (116) and (112))} \end{align}

(124) \begin{align} \|f\|_{H^1_tL^{\infty }_{{\boldsymbol{x}}}}\|g\|_{L^{\infty }_tH^1_{{\boldsymbol{x}}}} &\leq C(T) (\!\left \| f \right \|_{H^{1}_tH^{2}_{{\boldsymbol{x}}}} + \!\left \| f \right \|_{H^{2}_tL^2_{{\boldsymbol{x}}}}) \cdot \|g\|_{H^{1}_tH^{1}_{{\boldsymbol{x}}}} \textrm{ (use (116) and (112))} \end{align}

(125) \begin{align} \|f\|_{L^{\infty }_tW^{1, \infty }_{{\boldsymbol{x}}}}\|g\|_{H^{1}_tL^{2}_{{\boldsymbol{x}}}} &\leq C(T) \|f\|_{H^1_tH^3_{{\boldsymbol{x}}}} \cdot \|g\|_{H^{1}_tL^{2}_{{\boldsymbol{x}}}} \textrm{ (use (115))} \end{align}

(126) \begin{align} \|f\|_{L^{\infty }_tL^{\infty }_{{\boldsymbol{x}}}}\|g\|_{H^{1}_tH^{1}_{{\boldsymbol{x}}}} &\leq C(T) \|f\|_{H^1_tH^2_{{\boldsymbol{x}}}} \cdot \|g\|_{H^{1}_tH^{1}_{{\boldsymbol{x}}}} \textrm{ (use (115))}. \end{align}

Combining (122)–(126), we obtain (107).

Remark 5.1. It is useful to rewrite (107) with the norm of space $X_{\boldsymbol{Q}}$ :

(127) \begin{equation} \|fg\|_{H^1_tH^1_{{\boldsymbol{x}}}} \leq C(T) \!\left \| f \right \|_{X_{\boldsymbol{Q}}} \!\left \| g \right \|_{H^{1}_tH^{1}_{{\boldsymbol{x}}}}. \end{equation}

Lemma 5.1. For all $R\gt 0$ , there exists a time $T \gt 0$ such that for all $({\boldsymbol{u}}^{(i)}_h,{\boldsymbol{Q}}^{(i)}_h, \boldsymbol{\omega }^{(i)},{\boldsymbol{V}}^{(i)}) \in B_X(0, R)=\{({\boldsymbol{u}}_h,{\boldsymbol{Q}}_h, \boldsymbol{\omega },{\boldsymbol{V}}) \in X \big | \|({\boldsymbol{u}}_h,{\boldsymbol{Q}}_h, \boldsymbol{\omega },{\boldsymbol{V}})\|_{X} \leq R\}$ , $i=1,2$ , then

(128) \begin{eqnarray} &&\|\mathcal{J}({\boldsymbol{u}}^{(1)}_h,{\boldsymbol{Q}}^{(1)}_h, \boldsymbol{\omega }^{(1)},{\boldsymbol{V}}^{(1)}) - \mathcal{J}({\boldsymbol{u}}^{(2)}_h,{\boldsymbol{Q}}^{(2)}_h, \boldsymbol{\omega }^{(2)},{\boldsymbol{V}}^{(2)})\|_{Y} \nonumber \\[5pt] && \hspace{3.0cm} \leq C(T, R)\|({\boldsymbol{u}}^{(1)}_h,{\boldsymbol{Q}}^{(1)}_h, \boldsymbol{\omega }^{(1)},{\boldsymbol{V}}^{(1)}) - ({\boldsymbol{u}}^{(2)}_h,{\boldsymbol{Q}}^{(2)}_h, \boldsymbol{\omega }^{(2)},{\boldsymbol{V}}^{(2)})\|_{X}. \end{eqnarray}

Moreover, the constant coefficient $C(T, R) \rightarrow 0$ when $T \rightarrow 0$ .

Proof. $\,$

STEP 1. We first establish the Lipschitz continuity of ${\boldsymbol{u}}_{\textrm{os}}$ , the solution of (95),(96),(97), with respect to $\boldsymbol{\alpha }$ and $\boldsymbol{\omega }$ . For given $\boldsymbol{\alpha }^{(i)}(t)$ and $\boldsymbol{\omega }^{(i)}(t)$ , $i=1,2$ , such that $\boldsymbol{\alpha }^{(1)}(0)=\boldsymbol{\alpha }^{(2)}(0)$ , one has

\begin{align*} &-\eta \Delta \!\left ({\boldsymbol{u}}_{\textrm{os}}^{(1)} -{\boldsymbol{u}}_{\textrm{os}}^{(2)} \right ) + \nabla \!\left ( p_{\textrm{os}}^{(1)} - p_{\textrm{os}}^{(2)} \right ) = 0 \textrm{ in } \Omega, \\[5pt] &{\boldsymbol{u}}_{\textrm{os}}^{(1)} -{\boldsymbol{u}}_{\textrm{os}}^{(2)} = \!\left ( u_{\textrm{sq}}(\boldsymbol{\alpha }^{(1)}(t),{\boldsymbol{x}}) - u_{\textrm{sq}}(\boldsymbol{\alpha }^{(2)}(t),{\boldsymbol{x}}) \right )\boldsymbol{\tau } + (\boldsymbol{\omega }^{(1)}(t) - \boldsymbol{\omega }^{(2)}(t)) \times{\boldsymbol{x}} \textrm{ on } \partial{\mathcal{P}},\\[5pt] &{\boldsymbol{u}}_{\textrm{os}}^{(1)},{\boldsymbol{u}}_{\textrm{os}}^{(2)} \textrm{ periodic in }\Pi. \end{align*}

Due to the stability of the Stokes operator (similar to [58, Theorem IV.6.1])

(129) \begin{equation} \|{\boldsymbol{u}}_{\textrm{os}}\|_{H^2(\Omega )}\leq C \eta ^{-1} \!\left (\|u_{\textrm{sq}}\|_{L^{2}(\partial \mathcal{P})}+|{\boldsymbol{\omega }}(t)|\right ) \end{equation}

and smooth dependence of $u_{\textrm{sq}}$ in $\alpha (t)$ , we have

(130) \begin{align} &\|{\boldsymbol{u}}_{\textrm{os}}^{(1)} -{\boldsymbol{u}}_{\textrm{os}}^{(2)}\|_{H^2_tL^2_{{\boldsymbol{x}}}} + \|{\boldsymbol{u}}_{\textrm{os}}^{(1)} -{\boldsymbol{u}}_{\textrm{os}}^{(2)}\|_{H^1_tH^2_{{\boldsymbol{x}}}} \nonumber \\[5pt] &\qquad \leq C\eta ^{-1}\!\left (\|\boldsymbol{\alpha }^{(1)} - \boldsymbol{\alpha }^{(2)}\|_{H^2(0,T)} + \|\boldsymbol{\omega }^{(1)} - \boldsymbol{\omega }^{(2)}\|_{H^2(0,T)}\right ). \end{align}

(131) \begin{align} &\|{\boldsymbol{u}}_{\textrm{os}}^{(1)} -{\boldsymbol{u}}_{\textrm{os}}^{(2)}\|_{H^1_tH^2_{{\boldsymbol{x}}}} \leq C\eta ^{-1}\!\left (\|\boldsymbol{\alpha }^{(1)} - \boldsymbol{\alpha }^{(2)}\|_{H^1(0,T)} + \|\boldsymbol{\omega }^{(1)} - \boldsymbol{\omega }^{(2)}\|_{H^1(0,T)}\right ). \end{align}

Since $\boldsymbol{\alpha }^{(i)}(t) = \boldsymbol{\alpha }^{{(i)}}(0)+\int _{0}^{t} \boldsymbol{\omega }^{(i)}(\tau )\times \boldsymbol{\alpha }^{(i)}(\tau ) \textrm{d}\tau$ (see (17)) and $|\boldsymbol{\alpha }^{(i)}(t)|=1$ , $i=1,2$ as well as $\int \limits _0^{T}|h(t)|^2\,\textrm{d}t\leq T^2\int \limits _0^{T}|h_t(0,T)|^2\,\textrm{d}t$ for all $h\in H^1(0,T)$ , one gets

(132) \begin{align} \|\boldsymbol{\alpha }^{(1)} - \boldsymbol{\alpha }^{(2)}\|_{H^2(0,T)} \leq T\|\boldsymbol{\omega }^{(1)} - \boldsymbol{\omega }^{(2)}\|_{H^2(0,T)}, \end{align}

(133) \begin{align} \|\boldsymbol{\omega }^{(1)} - \boldsymbol{\omega }^{(2)}\|_{H^1(0,T)} \leq C(T)\|\boldsymbol{\omega }^{(1)} - \boldsymbol{\omega }^{(2)}\|_{H^2(0,T)}. \end{align}

Then (130) and (131) become

(134) \begin{align} \|{\boldsymbol{u}}_{\textrm{os}}^{(1)} -{\boldsymbol{u}}_{\textrm{os}}^{(2)}\|_{H^2_tL^2_{{\boldsymbol{x}}}} + \|{\boldsymbol{u}}_{\textrm{os}}^{(1)} -{\boldsymbol{u}}_{\textrm{os}}^{(2)}\|_{H^1_tH^2_{{\boldsymbol{x}}}} \leq C\eta ^{-1}(1+T)\|\boldsymbol{\omega }^{(1)} - \boldsymbol{\omega }^{(2)}\|_{H^2(0,T)}, \end{align}

(135) \begin{align} \|{\boldsymbol{u}}_{\textrm{os}}^{(1)} -{\boldsymbol{u}}_{\textrm{os}}^{(2)}\|_{H^1_tH^2_{{\boldsymbol{x}}}} \leq C(T)\eta ^{-1}\|\boldsymbol{\omega }^{(1)} - \boldsymbol{\omega }^{(2)}\|_{H^2(0,T)}. \end{align}

Applying (129) to $\partial ^{k}_t{\boldsymbol{u}}^{(i)}$ with $k=0,1,2$ and using the definition of time-independent ${\boldsymbol{Q}}_{\textrm{os}}$ (100), there is a $C\gt 0$ depending on $\Omega$ and $q_{\infty }$ such that

(136) \begin{align} \|{\boldsymbol{u}}_{\textrm{os}}^{(i)}\|_{H^2_tL^2_{{\boldsymbol{x}}}} + \|{\boldsymbol{u}}_{\textrm{os}}^{(i)}\|_{H^1_tH^2_{{\boldsymbol{x}}}} \leq C\eta ^{-1}(\|{\boldsymbol{\omega }}^{(i)}\|_{H^2(0,T)}+1), \quad i=1,2. \end{align}

(137) \begin{align} \|{\boldsymbol{Q}}_{\textrm{os}}\|_{H^2_tH^1_{{\boldsymbol{x}}}} + \|{\boldsymbol{Q}}_{\textrm{os}}\|_{H^1_tH^3_{{\boldsymbol{x}}}} \leq C. \end{align}

(Though ${\boldsymbol{Q}}_{\textrm{os}}$ is independent of time, here we use its $H^2_tH^1_{{\boldsymbol{x}}}$ and $H^1_tH^3_{{\boldsymbol{x}}}$ norms for the clarity of arguments below.) We will also need the following inequality:

(138) \begin{equation} \|{\boldsymbol{Q}}^{(i)}\|_{H^1_tH^3_{{\boldsymbol{x}}}}+\|{\boldsymbol{Q}}^{(i)}\|_{H^2_tH^1_{{\boldsymbol{x}}}}\leq C+R, \quad i=1,2. \end{equation}

Indeed,

\begin{eqnarray*} &&\|{\boldsymbol{Q}}^{(i)}\|_{H^1_tH^3_{{\boldsymbol{x}}}}+\|{\boldsymbol{Q}}^{(i)}\|_{H^2_tH^1_{{\boldsymbol{x}}}}\leq \|{\boldsymbol{Q}}^{(i)}-{\boldsymbol{Q}}_{\textrm{os}}\|_{H^1_tH^3_{{\boldsymbol{x}}}}+\|{\boldsymbol{Q}}^{(i)}-{\boldsymbol{Q}}_{\textrm{os}}\|_{H^2_tH^1_{{\boldsymbol{x}}}}\\[5pt] &&\hspace{240 pt}+ \|{\boldsymbol{Q}}_{\textrm{os}}\|_{H^1_tH^3_{{\boldsymbol{x}}}}+\|{\boldsymbol{Q}}_{\textrm{os}}\|_{H^2_tH^1_{{\boldsymbol{x}}}}\\[5pt] &&\hspace{130 pt}\leq \|{\boldsymbol{Q}}^{(i)}_{h}\|_{X_{\boldsymbol{Q}}}+ \|{\boldsymbol{Q}}_{\textrm{os}}\|_{H^1_tH^3_{{\boldsymbol{x}}}}+\|{\boldsymbol{Q}}_{\textrm{os}}\|_{H^2_tH^1_{{\boldsymbol{x}}}} \\[5pt] && \hspace{130 pt}\leq C+R. \end{eqnarray*}

STEP 2. Here, we establish the following inequality:

(139) \begin{equation} \|\rho ^{-1}P_\sigma \nabla \cdot \left [\sigma _{\textrm{ela}}({\boldsymbol{Q}}^{(1)})-\sigma _{\textrm{ela}}({\boldsymbol{Q}}^{(2)})\right ]\|_{H^1_tL^{2}_{\sigma,{\boldsymbol{x}}}}\leq C(T)\|{\boldsymbol{Q}}^{(1)}_h-{\boldsymbol{Q}}^{(2)}_h\|_{X_{\boldsymbol{Q}}}. \end{equation}

To this end, we first note that since $P_{\sigma }\nabla \cdot \;:\; H^{1}(\Omega ) \rightarrow L^2_{\sigma }(\Omega )$ is a bounded operator [60, Lemma II.2.5.2], the inequality (139) follows from

(140) \begin{equation} \|\sigma _{\textrm{ela}}({\boldsymbol{Q}}^{(1)})-\sigma _{\textrm{ela}}({\boldsymbol{Q}}^{(2)})\|_{H^1_tH^1_{{\boldsymbol{x}}}}\leq C(T)\|{\boldsymbol{Q}}^{(1)}_h-{\boldsymbol{Q}}^{(2)}_h\|_{X_{\boldsymbol{Q}}}. \end{equation}

We decompose $\sigma _{\textrm{ela}}$ into five parts $\sigma _{\textrm{ela}} = \sigma _K + \sigma _a + \sigma _{s}^1 + \sigma _{s}^2 + \sigma _{s}^3$ , where

\begin{align*} &\sigma _{K} = -K\nabla{\boldsymbol{Q}} \odot \nabla{\boldsymbol{Q}} \\[5pt] &\sigma _{a} = K({\boldsymbol{Q}}\Delta{\boldsymbol{Q}} - \Delta{\boldsymbol{Q}}{\boldsymbol{Q}}) \\[5pt] &\sigma _{s}^1 = -\dfrac{2\xi }{d}{\bf H} \\[5pt] &\sigma _{s}^2 = -\xi \!\left [{\bf H}{\boldsymbol{Q}} +{\boldsymbol{Q}}{\bf H}\right ] + \dfrac{2\xi }{d}\textrm{Tr}({\boldsymbol{Q}}{\bf H}) \\[5pt] &\sigma _{s}^3 = 2\xi \!\left [{\boldsymbol{Q}}\textrm{Tr}({\boldsymbol{Q}}{\bf H})\right ]. \end{align*}

Here, $\sigma _{s}^{1}, \sigma _{s}^{2}, \sigma _{s}^{3}$ are the linear, bilinear and trilinear part of $\sigma _{s}\;:\!=\;\sigma ^{1}_s+\sigma ^{2}_s+\sigma ^{3}_s$ , respectively.

Part 1: $\sigma _{K}({\boldsymbol{Q}})$ .

(141) \begin{eqnarray} && \hspace{-30 pt}\| \sigma _{K}({\boldsymbol{Q}}^{(1)}) - \sigma _{K}({\boldsymbol{Q}}^{(2)}) \|_{H^1_tH^1_{{\boldsymbol{x}}}} =\| (K\nabla{\boldsymbol{Q}}^{(1)} \odot \nabla{\boldsymbol{Q}}^{(1)} - K\nabla{\boldsymbol{Q}}^{(2)}\odot \nabla{\boldsymbol{Q}}^{(2)}) \| _{H^1_tH^1_{{\boldsymbol{x}}}} \nonumber \\[5pt] &&\hspace{15 pt}\leq CK \!\left (\| \nabla{\boldsymbol{Q}}^{(1)}\odot \nabla ({\boldsymbol{Q}}^{(1)} -{\boldsymbol{Q}}^{(2)})\|_{H^1_tH^1_{{\boldsymbol{x}}}} + \| \nabla ({\boldsymbol{Q}}^{(1)} -{\boldsymbol{Q}}^{(2)})\odot \nabla{\boldsymbol{Q}}^{(2)} \| _{H^1_tH^1_{{\boldsymbol{x}}}}\right ) \nonumber \\[5pt] &&\hspace{15 pt}= CK \!\left (\| \nabla{\boldsymbol{Q}}^{(1)}\odot \nabla ({\boldsymbol{Q}}^{(1)}_h -{\boldsymbol{Q}}^{(2)}_h)\| _{H^1_tH^1_{{\boldsymbol{x}}}} + \| \nabla ({\boldsymbol{Q}}^{(1)}_h -{\boldsymbol{Q}}^{(2)}_h)\odot \nabla{\boldsymbol{Q}}^{(2)} \| _{H^1_tH^1_{{\boldsymbol{x}}}}\right ). \end{eqnarray}

Using (106), (137) and that $({\boldsymbol{u}}^{(i)}_h,{\boldsymbol{Q}}^{(i)}_h, \boldsymbol{\omega }^{(i)},{\boldsymbol{V}}^{(i)}) \in B_X(0, R)$ , we get

(142) \begin{align} & CK \!\left (\| \nabla{\boldsymbol{Q}}^{(1)}\odot \nabla ({\boldsymbol{Q}}^{(1)}_h -{\boldsymbol{Q}}^{(2)}_h)\|_{H^1_tH^1_{{\boldsymbol{x}}}} + \| \nabla ({\boldsymbol{Q}}^{(1)}_h -{\boldsymbol{Q}}^{(2)}_h)\odot \nabla{\boldsymbol{Q}}^{(2)} \|_{H^1_tH^1_{{\boldsymbol{x}}}}\right ) \nonumber \\[5pt] &\quad \leq C(T)K\!\left ( \|\nabla ({\boldsymbol{Q}}^{(1)}_h -{\boldsymbol{Q}}^{(2)}_h)\|_{H^1_tH^2_{{\boldsymbol{x}}}} + \|\nabla ({\boldsymbol{Q}}^{(1)}_h -{\boldsymbol{Q}}^{(2)}_h)\|_{H^2_tL^2_{{\boldsymbol{x}}}} \right ) \nonumber \\[5pt] &\qquad \qquad \!\left [\!\left ( \|\nabla{\boldsymbol{Q}}^{(1)}\|_{H^1_tH^2_{{\boldsymbol{x}}}} + \|\nabla{\boldsymbol{Q}}^{(1)}\|_{H^2_tL^2_{{\boldsymbol{x}}}} \right ) + \!\left ( \|\nabla{\boldsymbol{Q}}^{(2)}\|_{H^1_tH^2_{{\boldsymbol{x}}}} + \|\nabla{\boldsymbol{Q}}^{(2)}\|_{H^2_tL^2_{{\boldsymbol{x}}}} \right )\right ] \nonumber \\[5pt] &\quad \leq C(T)K(R + 1)\|{\boldsymbol{Q}}^{(1)}_h -{\boldsymbol{Q}}^{(2)}_h\|_{X_{\boldsymbol{Q}}} \nonumber \\[5pt] &\quad \leq C(T)\|{\boldsymbol{Q}}^{(1)}_h -{\boldsymbol{Q}}^{(2)}_h\|_{X_{\boldsymbol{Q}}}. \end{align}

We note that the generic constant $C(T)$ may change from line to line and may depend on, for example, $K$ , $R$ , $C$ from (137) and $T$ (but recall that notation $C(T)$ also means that $C(T)\to 0$ as $T\to 0$ ). We sometimes do not merge a parameter, for example, $K$ in the second line of the above chain of inequalities, to indicate what we used to obtain a bound.

Part 2: $\sigma _{a}({\boldsymbol{Q}})$ .

(143) \begin{eqnarray} &&\hspace{-40pt} \| \sigma _{a}({\boldsymbol{Q}}^{(1)}) - \sigma _{a}({\boldsymbol{Q}}^{(2)}) \|_{H^1_tH^1_{{\boldsymbol{x}}}} \nonumber \\[5pt] &&\hspace{10pt}\leq K\|{\boldsymbol{Q}}^{(1)} \Delta{\boldsymbol{Q}}^{(1)} -{\boldsymbol{Q}}^{(2)} \Delta{\boldsymbol{Q}}^{(2)} \|_{H^1_tH^1_{{\boldsymbol{x}}}} + K\| \Delta{\boldsymbol{Q}}^{(1)}{\boldsymbol{Q}}^{(1)} - \Delta{\boldsymbol{Q}}^{(2)}{\boldsymbol{Q}}^{(2)} \|_{H^1_tH^1_{{\boldsymbol{x}}}}. \end{eqnarray}

Applying (107) for the first term in the right-hand side of (143), one can get

(144) \begin{eqnarray} && \hspace{-30pt}\|{\boldsymbol{Q}}^{(1)}\Delta{\boldsymbol{Q}}^{(1)} -{\boldsymbol{Q}}^{(2)} \Delta{\boldsymbol{Q}}^{(2)} \|_{H^1_tH^1_{{\boldsymbol{x}}}} \nonumber \\[5pt] && \leq \| ({\boldsymbol{Q}}^{(1)}_h -{\boldsymbol{Q}}^{(2)}_h) \Delta ({\boldsymbol{Q}}^{(1)}_h +{\boldsymbol{Q}}_{\textrm{os}}) \|_{H^1_tH^1_{{\boldsymbol{x}}}} + \|({\boldsymbol{Q}}^{(2)}_h +{\boldsymbol{Q}}_{\textrm{os}}) \Delta ({\boldsymbol{Q}}^{(1)}_h -{\boldsymbol{Q}}^{(2)}_h ) \|_{H^1_tH^1_{{\boldsymbol{x}}}} \nonumber \\[5pt] && \leq C(T) \|{\boldsymbol{Q}}^{(1)}_h -{\boldsymbol{Q}}^{(2)}_h \|_{X_{\boldsymbol{Q}}} \cdot \| \Delta ({\boldsymbol{Q}}^{(1)}_h +{\boldsymbol{Q}}_{\textrm{os}})\|_{H^{1}_tH^{1}_{{\boldsymbol{x}}}} \nonumber \\[5pt] && \hspace{60 pt}+\; C(T) \|{\boldsymbol{Q}}^{(2)}_h +{\boldsymbol{Q}}_{\textrm{os}} \|_{X_{\boldsymbol{Q}}} \| \Delta ({\boldsymbol{Q}}^{(1)}_h -{\boldsymbol{Q}}^{(2)}_h) \|_{H^{1}_tH^{1}_{{\boldsymbol{x}}}} \nonumber \\[5pt] &&\leq C(T)(R + 1) \|{\boldsymbol{Q}}^{(1)}_h -{\boldsymbol{Q}}^{(2)}_h\|_{X_{\boldsymbol{Q}}}. \end{eqnarray}

Applying same arguments for the second term in the right-hand side of (143), one can obtain

\begin{equation*} \|\sigma _{a}({\boldsymbol{Q}}^{(1)}) - \sigma _{a}({\boldsymbol{Q}}^{(2)}) \|_{H^1_tH^1_{{\boldsymbol{x}}}} \leq KC(T, R) \|{\boldsymbol{Q}}^{(1)}_h - {\boldsymbol{Q}}^{(2)}_h\|_{X_{\boldsymbol{Q}}}. \end{equation*}

Part 3: $\sigma _{s}^1({\boldsymbol{Q}})$ .

(145) \begin{eqnarray} &&\hspace{-30 pt}\| \sigma _{s}^1({\boldsymbol{Q}}^{(1)}) - \sigma _{s}^1({\boldsymbol{Q}}^{(2)}) \|_{H^1_tH^1_{{\boldsymbol{x}}}} =\dfrac{2\xi }{d}\|{\bf H}({\boldsymbol{Q}}^{(1)}) -{\bf H}({\boldsymbol{Q}}^{(2)})\| _{H^1_tH^1_{{\boldsymbol{x}}}} \nonumber \\[5pt] &&\hspace{10 pt}\leq \dfrac{2\xi K}{d}\| \Delta{\boldsymbol{Q}}^{(1)} - \Delta{\boldsymbol{Q}}^{(2)} \| _{H^1_tH^1_{{\boldsymbol{x}}}}+ \dfrac{2\xi |a|}{d}\|{\boldsymbol{Q}}^{(1)} -{\boldsymbol{Q}}^{(2)}\| _{H^1_tH^1_{{\boldsymbol{x}}}} \nonumber \\[5pt] &&\hspace{50 pt} +\dfrac{2\xi |c|}{d} \|{\boldsymbol{Q}}^{(1)}\textrm{Tr}(({\boldsymbol{Q}}^{(1)})^2) -{\boldsymbol{Q}}^{(2)}\textrm{Tr}(({\boldsymbol{Q}}^{(2)})^2) \| _{H^1_tH^1_{{\boldsymbol{x}}}}. \end{eqnarray}

The first two terms in the right-hand side of (145) are bounded as follows:

\begin{equation*} \dfrac {2\xi K}{d}\| \Delta {\boldsymbol{Q}}^{(1)} - \Delta {\boldsymbol{Q}}^{(2)} \| _{H^1_tH^1_{{\boldsymbol{x}}}}+ \dfrac {2\xi |a|}{d}\| {\boldsymbol{Q}}^{(1)} - {\boldsymbol{Q}}^{(2)} \| _{H^1_tH^1_{{\boldsymbol{x}}}}\leq C \|{\boldsymbol{Q}}^{(1)} - {\boldsymbol{Q}}^{(2)}\|_{X_{\boldsymbol{Q}}}. \end{equation*}

Next, we bound the third (cubic) term in the right-hand side of (145). Note

(146) \begin{align} \|{\boldsymbol{Q}}^{(1)}\textrm{Tr}(({\boldsymbol{Q}}^{(1)})^2) -{\boldsymbol{Q}}^{(2)}\textrm{Tr}(({\boldsymbol{Q}}^{(2)})^2) \| _{H^1_tH^1_{{\boldsymbol{x}}}}& \leq \| ({\boldsymbol{Q}}^{(1)} -{\boldsymbol{Q}}^{(2)})\textrm{Tr}({\boldsymbol{Q}}^{(1)}{\boldsymbol{Q}}^{(1)})\| _{H^1_tH^1_{{\boldsymbol{x}}}} \nonumber \\[5pt] &\hspace{10 pt}+\|{\boldsymbol{Q}}^{(2)}\textrm{Tr}(({\boldsymbol{Q}}^{(1)} -{\boldsymbol{Q}}^{(2)}){\boldsymbol{Q}}^{(1)})\|_{H^1_tH^1_{{\boldsymbol{x}}}}\nonumber \\[5pt]& \hspace{10 pt}+ \|{\boldsymbol{Q}}^{(2)}\textrm{Tr}({\boldsymbol{Q}}^{(2)}({\boldsymbol{Q}}^{(1)} -{\boldsymbol{Q}}^{(2)}))\|_{H^1_tH^1_{{\boldsymbol{x}}}} \end{align}

We show how to bound the first term in the right-hand side of (146). Other terms are bounded in the same way. Apply (127) twice to obtain

(147) \begin{align} &\| ({\boldsymbol{Q}}^{(1)} -{\boldsymbol{Q}}^{(2)})\textrm{Tr}({\boldsymbol{Q}}^{(1)}{\boldsymbol{Q}}^{(1)})\| _{H^1_tH^1_{{\boldsymbol{x}}}} \nonumber \\[5pt] &\hspace{20 pt}\leq C(T)\|{\boldsymbol{Q}}^{(1)}_h -{\boldsymbol{Q}}^{(2)}_h\|_{X_{\boldsymbol{Q}}} \||{\boldsymbol{Q}}^{(1)}||{\boldsymbol{Q}}^{(1)}|\|_{H^1_tH^1_{{\boldsymbol{x}}}} \nonumber \\[5pt] &\hspace{20 pt} \leq C(T)\|{\boldsymbol{Q}}^{(1)}_h -{\boldsymbol{Q}}^{(2)}_h\|_{X_{\boldsymbol{Q}}}\!\left (\|{\boldsymbol{Q}}^{(1)}_h\|_{X_{\boldsymbol{Q}}}+\|{\boldsymbol{Q}}_{\textrm{os}}\|_{L^2_tH^3_{{\boldsymbol{x}}}}\right )\!\left (\|{\boldsymbol{Q}}^{(1)}_h\|_{H^1_tH^1_{{\boldsymbol{x}}}}+\|{\boldsymbol{Q}}_{\textrm{os}}\|_{ H^1_tH^1_{{\boldsymbol{x}}}}\right )\nonumber \\[5pt] &\hspace{20 pt} \leq C(T)(R+1)^2\|{\boldsymbol{Q}}^{(1)}_h -{\boldsymbol{Q}}^{(2)}_h\|_{X_{\boldsymbol{Q}}}. \end{align}

Thus, we have

(148) \begin{align} \|{\bf H}({\boldsymbol{Q}}^{(1)}) -{\bf H}({\boldsymbol{Q}}^{(2)}) \| _{H^1_tH^1_{{\boldsymbol{x}}}}\leq C(T,R) \|{\boldsymbol{Q}}^{(1)}_h -{\boldsymbol{Q}}^{(2)}_h\|_{X_{\boldsymbol{Q}}}, \end{align}

which, in view of (145), implies

(149) \begin{align} \| \sigma _{s}^1({\boldsymbol{Q}}^i) - \sigma _{s}^1({\boldsymbol{Q}}^j) \|_{H^1_tH^1_{{\boldsymbol{x}}}} \leq C(T,R) \|{\boldsymbol{Q}}^{(1)}_h -{\boldsymbol{Q}}^{(2)}_h\|_{X_{\boldsymbol{Q}}}. \end{align}

Part 4: $ \sigma _{s}^2({\boldsymbol{Q}})$ . In this part, we will need the following bound:

(150) \begin{align} \|{\bf H}({\boldsymbol{Q}}^{(i)})\|_{H^1_tH^1_{{\boldsymbol{x}}}} \leq C(T,R), \quad i=1,2. \end{align}

which can be obtained by applying same arguments as in (145)–(147) for $i=1,2$

\begin{align*} \|{\bf H}({\boldsymbol{Q}}^{(i)})\|_{H^1_tH^1_{{\boldsymbol{x}}}} &\leq C(\|\Delta{\boldsymbol{Q}}^{(i)}\|_{H^1_tH^1_{{\boldsymbol{x}}}} + \|{\boldsymbol{Q}}^{(i)}\|_{H^1_tH^1_{{\boldsymbol{x}}}} + \|{\boldsymbol{Q}}^{(i)}\textrm{Tr}(({\boldsymbol{Q}}^{(i)})^2)\|_{H^1_tH^1_{{\boldsymbol{x}}}})\\[5pt] & \leq C(T) (\|{\boldsymbol{Q}}^{(i)}_h\|_{X_{{\boldsymbol{Q}}}} + \|{\boldsymbol{Q}}^{(i)}_h\|_{X_{{\boldsymbol{Q}}}}^3+1)\leq C(T,R). \end{align*}

Now, we can estimate,

(151) \begin{align} &\| \sigma _{s}^2({\boldsymbol{Q}}^{(1)}) - \sigma _{s}^2({\boldsymbol{Q}}^{(2)}) \|_{H^1_tH^1_{{\boldsymbol{x}}}} \leq \xi \|{\boldsymbol{Q}}^{(1)}{\bf H}({\boldsymbol{Q}}^{(1)}) -{\boldsymbol{Q}}^{(2)}{\bf H}({\boldsymbol{Q}}^{(2)})\|_{H^1_tH^1_{{\boldsymbol{x}}}}\nonumber \\[5pt] &\qquad + \xi \|{\bf H}({\boldsymbol{Q}}^{(1)}){\boldsymbol{Q}}^{(1)} -{\bf H}({\boldsymbol{Q}}^{(2)}){\boldsymbol{Q}}^{(2)}\|_{H^1_tH^1_{{\boldsymbol{x}}}} \nonumber \\[5pt] &\qquad +\dfrac{2\xi }{d}\|\textrm{Tr}({\boldsymbol{Q}}^{(1)}{\bf H}({\boldsymbol{Q}}^{(1)}) -{\boldsymbol{Q}}^{(2)}{\bf H}({\boldsymbol{Q}}^{(2)}))\|_{H^1_tH^1_{{\boldsymbol{x}}}}\nonumber \\[5pt] &\qquad \leq C\|{\boldsymbol{Q}}^{(1)}{\bf H}({\boldsymbol{Q}}^{(1)}) -{\boldsymbol{Q}}^{(2)}{\bf H}({\boldsymbol{Q}}^{(2)})\|_{H^1_tH^1_{{\boldsymbol{x}}}}. \end{align}

Next, we use the triangle inequality, (127), (148) and (150):

\begin{align*} &\|{\boldsymbol{Q}}^{(1)}{\bf H}({\boldsymbol{Q}}^{(1)}) -{\boldsymbol{Q}}^{(2)}{\bf H}({\boldsymbol{Q}}^{(2)})\|_{H^1_tH^1_{{\boldsymbol{x}}}}\\[5pt] & \hspace{40 pt} \leq C(T) \|{\boldsymbol{Q}}^{(1)}_h -{\boldsymbol{Q}}^{(2)}_h \|_{X_{\boldsymbol{Q}}} \|{\bf H}({\boldsymbol{Q}}^{(1)}) \|_{H^{1}_tH^{1}_{{\boldsymbol{x}}}} \\[5pt] & \hspace{60 pt}+ C(T) \|{\boldsymbol{Q}}^{(2)}_h +{\boldsymbol{Q}}_{\textrm{os}} \|_{X_{\boldsymbol{Q}}}\|{\bf H}({\boldsymbol{Q}}^{(1)}) -{\bf H}({\boldsymbol{Q}}^{(2)}) \|_{H^{1}_tH^{1}_{{\boldsymbol{x}}}}\\[5pt] &\hspace{40 pt} \leq C(T,R)\|{\boldsymbol{Q}}^{(1)}_h -{\boldsymbol{Q}}^{(2)}_h \|_{X_{\boldsymbol{Q}}}. \end{align*}

Therefore,

(152) \begin{align} &\| \sigma _{s}^2({\boldsymbol{Q}}^{(1)}) - \sigma _{s}^2({\boldsymbol{Q}}^{(2)}) \|_{H^1_tH^1_{{\boldsymbol{x}}}}\leq C(T,R)\|{\boldsymbol{Q}}^{(1)}_h -{\boldsymbol{Q}}^{(2)}_h \|_{X_{\boldsymbol{Q}}}. \end{align}

Part 5: $ \sigma _{s}^3({\boldsymbol{Q}})$ .

(153) \begin{align} &\| \sigma _{s}^3({\boldsymbol{Q}}^{(1)}) - \sigma _{s}^3({\boldsymbol{Q}}^{(2)}) \|_{H^1_tH^1_{{\boldsymbol{x}}}} \leq 2\xi \|{\boldsymbol{Q}}^{(1)}\textrm{Tr}({\boldsymbol{Q}}^{(1)}{\bf H}({\boldsymbol{Q}}^{(1)})) -{\boldsymbol{Q}}^{(2)}\textrm{Tr}({\boldsymbol{Q}}^{(2)}{\bf H}({\boldsymbol{Q}}^{(2)}))\|_{H^1_tH^1_{{\boldsymbol{x}}}}\nonumber \\[5pt] & \hspace{140 pt}\leq C\| ({\boldsymbol{Q}}^{(1)} -{\boldsymbol{Q}}^{(2)})\textrm{Tr}({\boldsymbol{Q}}^{(1)}{\bf H}({\boldsymbol{Q}}^{(1)}))\| _{H^1_tH^1_{{\boldsymbol{x}}}} \nonumber \\[5pt] &\hspace{160 pt}+ C\|{\boldsymbol{Q}}^{(2)}\textrm{Tr}(({\boldsymbol{Q}}^{(1)} -{\boldsymbol{Q}}^{(2)}){\bf H}({\boldsymbol{Q}}^{(1)}))\| _{H^1_tH^1_{{\boldsymbol{x}}}}\nonumber \\[5pt]& \hspace{160 pt}+ C\|{\boldsymbol{Q}}^{(2)}\textrm{Tr}({\boldsymbol{Q}}^{(2)}({\bf H}({\boldsymbol{Q}}^{(1)}) -{\bf H}({\boldsymbol{Q}}^{(2)})))\|_{H^1_tH^1_{{\boldsymbol{x}}}}. \end{align}

Next, applying the same arguments as in (147) with bounds (148) and (150) we get

(154) \begin{align} &\| \sigma _{s}^3({\boldsymbol{Q}}^{(1)}) - \sigma _{s}^3({\boldsymbol{Q}}^{(2)}) \|_{H^1_tH^1_{{\boldsymbol{x}}}}\leq C(T,R)\|{\boldsymbol{Q}}^{(1)}_h -{\boldsymbol{Q}}^{(2)}_h \|_{X_{\boldsymbol{Q}}}. \end{align}

STEP 3. Here, we establish the following inequality:

(155) \begin{eqnarray} &&\|S^2(\nabla{\boldsymbol{u}}^{(1)},{\boldsymbol{Q}}^{(1)})-S^2(\nabla{\boldsymbol{u}}^{(2)},{\boldsymbol{Q}}^{(2)})\|_{H^1_tH^{1}_{{\boldsymbol{x}}}} \nonumber \\[5pt] && \qquad \leq C(T)\!\left (\|{\boldsymbol{Q}}^{(1)}_h-{\boldsymbol{Q}}^{(2)}_h\|_{X_{\boldsymbol{Q}}} + \|{\boldsymbol{u}}^{(1)}_h-{\boldsymbol{u}}^{(2)}_h\|_{X_{\boldsymbol{u}}} + \|\boldsymbol{\omega }^{(1)} - \boldsymbol{\omega }^{(1)}\|_{H^2(0,T)} \right ). \end{eqnarray}

To this end, similar to how we treated $\sigma _s$ in STEP 2, we split $S(\nabla{\boldsymbol{u}},{\boldsymbol{Q}})$ into three parts:

(156) \begin{eqnarray} S(\nabla{\boldsymbol{u}},{\boldsymbol{Q}}) &&=(\xi{\boldsymbol{D}} +{\boldsymbol{A}})\!\left ({\boldsymbol{Q}}+\frac{\mathbb I}{d}\right )+ \left ({\boldsymbol{Q}}+\frac{\mathbb I}{d}\right )(\xi{\boldsymbol{D}}-{\boldsymbol{A}})-2\xi \!\left ({\boldsymbol{Q}}+\frac{\mathbb I}{d}\right )\textrm{tr}({\boldsymbol{Q}}\nabla{\boldsymbol{u}}) \nonumber \\[5pt] && = S^1 + S^2 + S^3, \end{eqnarray}

where

\begin{align*} & S^1 (\nabla{\boldsymbol{u}})= \frac{2\xi }{d}{\boldsymbol{D}} \\[5pt] & S^2 (\nabla{\boldsymbol{u}},{\boldsymbol{Q}})= \xi ({\boldsymbol{D}}{\boldsymbol{Q}} +{\boldsymbol{Q}}{\boldsymbol{D}}) + ({\boldsymbol{A}}{\boldsymbol{Q}} -{\boldsymbol{Q}}{\boldsymbol{A}}) - \frac{2\xi }{d}\textrm{tr}({\boldsymbol{Q}}\nabla{\boldsymbol{u}})\\[5pt] & S^3 (\nabla{\boldsymbol{u}},{\boldsymbol{Q}})= - 2\xi{\boldsymbol{Q}}\textrm{tr}({\boldsymbol{Q}}\nabla{\boldsymbol{u}}) \end{align*}

are correspondingly the linear, bilinear and trilinear part of $S$ . First note that using (134), (135) and (136), one gets

(157) \begin{align} \|{\boldsymbol{u}}^{(1)} -{\boldsymbol{u}}^{(2)}\|_{H^1_tH^2_{{\boldsymbol{x}}} \cap H^2_tL^2_{{\boldsymbol{x}}}} &\leq \|{\boldsymbol{u}}^{(1)}_h -{\boldsymbol{u}}^{(2)}_h\|_{X_{\boldsymbol{u}}} + \|{\boldsymbol{u}}^{(1)}_{\textrm{os}} -{\boldsymbol{u}}^{(2)}_{\textrm{os}}\|_{H^1_tH^2_{{\boldsymbol{x}}} \cap H^2_tL^2_{{\boldsymbol{x}}}} \nonumber \\[5pt] & \leq \|{\boldsymbol{u}}^{(1)}_h -{\boldsymbol{u}}^{(2)}_h\|_{X_{\boldsymbol{u}}} + C\eta ^{-1}(1+T)\|\boldsymbol{\omega }^{(1)} - \boldsymbol{\omega }^{(2)}\|_{H^2(0,T)} \end{align}

and

(158) \begin{align} \|{\boldsymbol{u}}^{(i)}\|_{H^1_tH^2_{{\boldsymbol{x}}} \cap H^2_tL^2_{{\boldsymbol{x}}}} &\leq \|{\boldsymbol{u}}^{(i)}_h\|_{X_{\boldsymbol{u}}} + \|{\boldsymbol{u}}^{(i)}_{\textrm{os}}\|_{H^1_tH^2_{{\boldsymbol{x}}} \cap H^2_tL^2_{{\boldsymbol{x}}}} \nonumber \\[5pt] & \leq \|{\boldsymbol{u}}^{(i)}_h\|_{X_{\boldsymbol{u}}} + C\eta ^{-1}(\|{\boldsymbol{\omega }}^{(i)}\|_{H^2(0,T)}+1), \quad i=1,2. \end{align}

Part 1: $S^1(\nabla{\boldsymbol{u}},{\boldsymbol{Q}})$ .

Since ${\boldsymbol{D}}({\boldsymbol{u}}) = \frac{1}{2}(\nabla{\boldsymbol{u}} + (\nabla{\boldsymbol{u}})^{T})$ , using (135) one gets

(159) \begin{align} & \|{\boldsymbol{D}}({\boldsymbol{u}}^{(1)}) -{\boldsymbol{D}}({\boldsymbol{u}}^{(1)}) \|_{H^1_tH^1_{{\boldsymbol{x}}}} \leq \|{\boldsymbol{u}}^{(1)} -{\boldsymbol{u}}^{(2)}\|_{H^1_tH^2_{{\boldsymbol{x}}}}\nonumber \\[5pt] & \qquad \leq \|{\boldsymbol{u}}^{(1)}_h -{\boldsymbol{u}}^{(2)}_h\|_{H^1_tH^2_{{\boldsymbol{x}}}}+\|{\boldsymbol{u}}^{(1)}_{\textrm{os}} -{\boldsymbol{u}}^{(2)}_{\textrm{os}}\|_{H^1_tH^2_{{\boldsymbol{x}}}}\nonumber \\[5pt] & \qquad \leq \|{\boldsymbol{u}}^{(1)}_h -{\boldsymbol{u}}^{(2)}_h\|_{X_{\boldsymbol{u}}} + C(T)\eta ^{-1}\|\boldsymbol{\omega }^{(1)} - \boldsymbol{\omega }^{(2)}\|_{H^2(0,T)}. \end{align}

Then

(160) \begin{eqnarray} &&\|S^1(\nabla{\boldsymbol{u}}^{(1)},{\boldsymbol{Q}}^{(1)})-S^1(\nabla{\boldsymbol{u}}^{(2)},{\boldsymbol{Q}}^{(2)})\|_{H^1_tH^{1}_{{\boldsymbol{x}}}} \nonumber \\[5pt] && \qquad \leq C(T)\!\left (\|{\boldsymbol{Q}}^{(1)}_h-{\boldsymbol{Q}}^{(2)}_h\|_{X_{\boldsymbol{Q}}} + \|{\boldsymbol{u}}^{(1)}_h-{\boldsymbol{u}}^{(2)}_h\|_{X_{\boldsymbol{u}}} + \|\boldsymbol{\omega }^{(1)} - \boldsymbol{\omega }^{(2)}\|_{H^2(0,T)} \right ). \end{eqnarray}

Part 2: $S^2(\nabla{\boldsymbol{u}},{\boldsymbol{Q}})$ .

(161) \begin{eqnarray} &&\|S^2(\nabla{\boldsymbol{u}}^{(1)},{\boldsymbol{Q}}^{(1)})-S^2(\nabla{\boldsymbol{u}}^{(2)},{\boldsymbol{Q}}^{(2)})\|_{H^1_tH^1_{{\boldsymbol{x}}}}\leq C \|\nabla{\boldsymbol{u}}^{(1)}{\boldsymbol{Q}}^{(1)}-\nabla{\boldsymbol{u}}^{(2)}{\boldsymbol{Q}}^{(2)}\|_{H^1_tH^1_{{\boldsymbol{x}}}}\nonumber \\[5pt] &&\hspace{60 pt}\leq C\|(\nabla{\boldsymbol{u}}^{(1)}-\nabla{\boldsymbol{u}}^{(2)}){\boldsymbol{Q}}^{(1)}\|_{H^1_tH^1_{{\boldsymbol{x}}}}+C \|\nabla{\boldsymbol{u}}^{(2)}({\boldsymbol{Q}}^{(1)}-{\boldsymbol{Q}}^{(2)})\|_{H^1_tH^1_{{\boldsymbol{x}}}} \end{eqnarray}

Apply (107), (138) and (157) to obtain

\begin{eqnarray*} &&\|(\nabla{\boldsymbol{u}}^{(1)}-\nabla{\boldsymbol{u}}^{(2)}){\boldsymbol{Q}}^{(1)}\|_{H^1_tH^1_{{\boldsymbol{x}}}}\\[5pt] &&\hspace{70 pt} \leq C(T)\|\nabla{\boldsymbol{u}}^{(1)}-\nabla{\boldsymbol{u}}^{(2)}\|_{H^1_tH^1_{{\boldsymbol{x}}}} \!\left (\|{\boldsymbol{Q}}^{(1)} \|_{H^1_tH^3_{{\boldsymbol{x}}}} +\|{\boldsymbol{Q}}^{(1)} \|_{H^2_tH^1_{{\boldsymbol{x}}}} \right ) \\[5pt] &&\hspace{70 pt}\leq C(T)\|{\boldsymbol{u}}^{(1)}-{\boldsymbol{u}}^{(2)}\|_{H^1_tH^2_{{\boldsymbol{x}}}} \!\left (\|{\boldsymbol{Q}}^{(1)} \|_{H^1_tH^3_{{\boldsymbol{x}}}} + \|{\boldsymbol{Q}}^{(1)} \|_{H^2_tH^1_{{\boldsymbol{x}}}} \right )\\[5pt] && \hspace{70 pt}\leq C(T) (R + 1) \!\left ( \|{\boldsymbol{u}}^{(1)}_h -{\boldsymbol{u}}^{(2)}_h\|_{X_{\boldsymbol{u}}} + \|\boldsymbol{\omega }^{(1)} - \boldsymbol{\omega }^{(2)}\|_{H^2(0,T)} \right ). \end{eqnarray*}

Applying the similar arguments to the second term in the right-hand side of (161), we obtain

(162) \begin{eqnarray} &&\|S^2(\nabla{\boldsymbol{u}}^{(1)},{\boldsymbol{Q}}^{(1)})-S^2(\nabla{\boldsymbol{u}}^{(2)},{\boldsymbol{Q}}^{(2)})\|_{H^1_tH^{1}_{{\boldsymbol{x}}}} \nonumber \\[5pt] && \qquad \leq C(T)\!\left (\|{\boldsymbol{Q}}^{(1)}_h-{\boldsymbol{Q}}^{(2)}_h\|_{X_{\boldsymbol{Q}}} + \|{\boldsymbol{u}}^{(1)}_h-{\boldsymbol{u}}^{(2)}_h\|_{X_{\boldsymbol{u}}} + \|\boldsymbol{\omega }^{(1)} - \boldsymbol{\omega }^{(2)}\|_{H^2(0,T)} \right ). \end{eqnarray}

Part 3: $S^3(\nabla{\boldsymbol{u}},{\boldsymbol{Q}})$ .

(163) \begin{eqnarray} &&\hspace{-20 pt}\|S^3(\nabla{\boldsymbol{u}}^{(1)},{\boldsymbol{Q}}^{(1)})-S^3(\nabla{\boldsymbol{u}}^{(2)},{\boldsymbol{Q}}^{(2)})\|_{H^1_tH^1_{{\boldsymbol{x}}}}\nonumber \\[5pt] &&\hspace{60 pt}\leq 2\xi \|{\boldsymbol{Q}}^{(1)}\textrm{tr}\!\left (\nabla{\boldsymbol{u}}^{(1)}{\boldsymbol{Q}}^{(1)}\right )-{\boldsymbol{Q}}^{(2)}\textrm{tr}\!\left (\nabla{\boldsymbol{u}}^{(2)}{\boldsymbol{Q}}^{(2)}\right )\|_{H^1_tH^1_{{\boldsymbol{x}}}}\nonumber \\[5pt] &&\hspace{60 pt}\leq 2\xi \|\!\left ({\boldsymbol{Q}}^{(1)} -{\boldsymbol{Q}}^{(2)}\right )\textrm{tr}\!\left (\nabla{\boldsymbol{u}}^{(1)}{\boldsymbol{Q}}^{(1)}\right )\|_{H^1_tH^1_{{\boldsymbol{x}}}}\nonumber \\[5pt] &&\hspace{70 pt} + 2\xi \|{\boldsymbol{Q}}^{(2)}\textrm{tr}\!\left (\nabla \!\left ({\boldsymbol{u}}^{(1)} -{\boldsymbol{u}}^{(2)}\right ){\boldsymbol{Q}}^{(1)}\right )\|_{H^1_tH^1_{{\boldsymbol{x}}}} \nonumber \\[5pt] &&\hspace{70 pt}+ 2\xi \|{\boldsymbol{Q}}^{(2)}\textrm{tr}\!\left (\nabla{\boldsymbol{u}}^{(2)}\!\left ({\boldsymbol{Q}}^{(1)} -{\boldsymbol{Q}}^{(2)}\right )\right )\|_{H^1_tH^1_{{\boldsymbol{x}}}}. \end{eqnarray}

Using same arguments as in (147) and taking into account (136) and (138), we obtain

\begin{eqnarray*} &&\|\!\left ({\boldsymbol{Q}}^{(1)} -{\boldsymbol{Q}}^{(2)}\right )\textrm{tr}\!\left (\nabla{\boldsymbol{u}}^{(1)}{\boldsymbol{Q}}^{(1)}\right )\|_{H^1_tH^1_{{\boldsymbol{x}}}} \\[5pt] &&\hspace{60 pt}\leq C(T) \|{\boldsymbol{Q}}^{(1)} -{\boldsymbol{Q}}^{(2)}\|_{H^1_tH^3_{{\boldsymbol{x}}} \cap H^2_tH^1_{{\boldsymbol{x}}}} \|{\boldsymbol{Q}}^{(1)}\|_{H^1_tH^3_{{\boldsymbol{x}}} \cap H^2_tH^1_{{\boldsymbol{x}}}}\|\nabla{\boldsymbol{u}}^{(1)}\|_{H^1_tH^1_{{\boldsymbol{x}}}}\\[5pt] &&\hspace{60 pt}\leq C(T)(R^2+1)\|{\boldsymbol{Q}}^{(1)}_h -{\boldsymbol{Q}}^{(2)}_h\|_{X_{\boldsymbol{Q}}}. \end{eqnarray*}

Applying similar arguments for the other two terms in the right-hand side of (163), we obtain

(164) \begin{eqnarray} &&\hspace{-20 pt}\|S^3(\nabla{\boldsymbol{u}}^{(1)},{\boldsymbol{Q}}^{(1)})-S^3(\nabla{\boldsymbol{u}}^{(2)},{\boldsymbol{Q}}^{(2)})\|_{H^1_tH^{1}_{{\boldsymbol{x}}}} \nonumber \\[5pt] && \hspace{1cm} \leq C(T,R)\!\left (\|{\boldsymbol{Q}}^{(1)}_h-{\boldsymbol{Q}}^{(2)}_h\|_{X_{\boldsymbol{Q}}} + \|{\boldsymbol{u}}^{(1)}_h-{\boldsymbol{u}}^{(2)}_h\|_{X_{\boldsymbol{u}}} + \|\boldsymbol{\omega }^{(1)} - \boldsymbol{\omega }^{(2)}\|_{H^2(0,T)} \right ). \end{eqnarray}

STEP 4. Finally, we show the Lipschitz properties of all the remaining terms in $\mathcal{J}$ .

Part 1: $P_{\sigma }(\nabla \cdot ({\boldsymbol{u}} \otimes{\boldsymbol{u}}))$ .

We again use the fact that $P_{\sigma }\nabla \cdot : H^{1}(\Omega ) \rightarrow L^2_{\sigma }(\Omega )$ is a bounded operator:

(165) \begin{eqnarray} &&\| P_{\sigma }(\nabla \cdot ({\boldsymbol{u}}^{(1)} \otimes{\boldsymbol{u}}^{(1)})) - P_{\sigma }(\nabla \cdot ({\boldsymbol{u}}^{(2)} \otimes{\boldsymbol{u}}^{(2)})) \|_{H^1_{t}L^2_{\sigma,{\boldsymbol{x}}}} \nonumber \\[5pt] && \qquad \leq \|{\boldsymbol{u}}^{(1)} \otimes{\boldsymbol{u}}^{(1)} -{\boldsymbol{u}}^{(2)} \otimes{\boldsymbol{u}}^{(2)} \|_{H^1_tH^1_{{\boldsymbol{x}}}} \nonumber \\[5pt] && \qquad =\| ({\boldsymbol{u}}^{(1)} -{\boldsymbol{u}}^{(2)}) \otimes{\boldsymbol{u}}^{(1)} -{\boldsymbol{u}}^{(2)} \otimes ({\boldsymbol{u}}^{(1)} -{\boldsymbol{u}}^{(2)}) \|_{H^1_tH^1_{{\boldsymbol{x}}}}. \end{eqnarray}

We apply the same arguments in Part 1 of STEP 2, that is, apply (106), along with (157) and (158), to get:

\begin{align*} & \|{\boldsymbol{u}}^{(1)}\otimes ({\boldsymbol{u}}^{(1)} -{\boldsymbol{u}}^{(2)}) \| _{H^1_tH^1_{{\boldsymbol{x}}}} +\| ({\boldsymbol{u}}^{(1)} -{\boldsymbol{u}}^{(2)})\otimes{\boldsymbol{u}}^{(2)} \| _{H^1_tH^1_{{\boldsymbol{x}}}} \\[5pt] &\quad \leq C(T)\!\left ( \|{\boldsymbol{u}}^{(1)} -{\boldsymbol{u}}^{(2)}\|_{H^1_tH^2_{{\boldsymbol{x}}}} + \|{\boldsymbol{u}}^{(1)} -{\boldsymbol{u}}^{(2)}\|_{H^2_tL^2_{{\boldsymbol{x}}}} \right ) \\[5pt] &\qquad \qquad \times \!\left [\!\left ( \|{\boldsymbol{u}}^{(1)}\|_{H^1_tH^2_{{\boldsymbol{x}}}} + \|{\boldsymbol{u}}^{(1)}\|_{H^2_tL^2_{{\boldsymbol{x}}}} \right ) + \!\left ( \|{\boldsymbol{u}}^{(2)}\|_{H^1_tH^2_{{\boldsymbol{x}}}} + \|{\boldsymbol{u}}^{(2)}\|_{H^2_tL^2_{{\boldsymbol{x}}}} \right )\right ] \\[5pt] & \quad \leq C(T) \!\left (\|{\boldsymbol{u}}^{(1)}_h -{\boldsymbol{u}}^{(2)}_h\|_{X_{\boldsymbol{u}}} + C\!\left (\|\boldsymbol{\omega }^{(1)} - \boldsymbol{\omega }^{(2)}\|_{H^2(0,T)}\right ) \right ) \\[5pt] &\qquad \qquad \times \!\left [\!\left ( \|{\boldsymbol{u}}^{(1)}_h\|_{X_{\boldsymbol{u}}} + C(\|{\boldsymbol{\omega }}^{(1)}\|_{H^2(0,T)}+1) \right ) + \!\left ( \|{\boldsymbol{u}}^{(2)}_h\|_{X_{\boldsymbol{u}}} + C(\|{\boldsymbol{\omega }}^{(2)}\|_{H^2(0,T)}+1) \right ) \right ] \\[5pt] & \quad \leq C(T) (R + 1) \!\left (\|({\boldsymbol{u}}^{(1)}_h,{\boldsymbol{Q}}^{(1)}_h, \boldsymbol{\omega }^{(1)},{\boldsymbol{V}}^{(1)}) - ({\boldsymbol{u}}^{(2)}_h,{\boldsymbol{Q}}^{(2)}_h, \boldsymbol{\omega }^{(2)},{\boldsymbol{V}}^{(2)})\|_{X} \right ). \end{align*}

Thus, we obtained

\begin{eqnarray*} &&\| P_{\sigma }(\nabla \cdot ({\boldsymbol{u}}^{(1)} \otimes{\boldsymbol{u}}^{(1)})) - P_{\sigma }(\nabla \cdot ({\boldsymbol{u}}^{(2)} \otimes{\boldsymbol{u}}^{(2)})) \|_{H^1_tL^2_{\sigma,{\boldsymbol{x}}}} \\[5pt] && \qquad \leq C(T,R) \!\left (\|({\boldsymbol{u}}^{(1)}_h,{\boldsymbol{Q}}^{(1)}_h, \boldsymbol{\omega }^{(1)},{\boldsymbol{V}}^{(1)}) - ({\boldsymbol{u}}^{(2)}_h,{\boldsymbol{Q}}^{(2)}_h, \boldsymbol{\omega }^{(2)},{\boldsymbol{V}}^{(2)})\|_{X} \right ). \end{eqnarray*}

Part 2: ${\boldsymbol{u}} \cdot \nabla{\boldsymbol{Q}}$ .

(166) \begin{eqnarray} &&\|{\boldsymbol{u}}^{(1)} \cdot \nabla{\boldsymbol{Q}}^{(1)} -{\boldsymbol{u}}^{(2)} \cdot \nabla{\boldsymbol{Q}}^{(2)} \|_{H^1_tH^1_{{\boldsymbol{x}}}} \nonumber \\[5pt] && \qquad \leq \| ({\boldsymbol{u}}^{(1)} -{\boldsymbol{u}}^{(2)}) \cdot \nabla{\boldsymbol{Q}}^{(1)} \|_{H^1_tH^1_{{\boldsymbol{x}}}} + \|{\boldsymbol{u}}^{(2)} \cdot \nabla ({\boldsymbol{Q}}^{(1)} -{\boldsymbol{Q}}^{(2)}) \|_{H^1_tH^1_{{\boldsymbol{x}}}}. \end{eqnarray}

Apply the same arguments as in Part 1 of STEP 2, that is, apply (106):

(167) \begin{align} & \|{\boldsymbol{u}}^{(1)}\cdot \nabla ({\boldsymbol{Q}}^{(1)} -{\boldsymbol{Q}}^{(2)})\|_{H^1_tH^1_{{\boldsymbol{x}}}} + \| ({\boldsymbol{u}}^{(1)} -{\boldsymbol{u}}^{(2)})\cdot \nabla{\boldsymbol{Q}}^{(2)} \|_{H^1_tH^1_{{\boldsymbol{x}}}} \nonumber \\[5pt] & \quad \leq C(T)\!\left ( \|{\boldsymbol{u}}^{(1)}\|_{H^1_tH^2_{{\boldsymbol{x}}}} + \|{\boldsymbol{u}}^{(1)}\|_{H^2_tL^2_{{\boldsymbol{x}}}} \right ) \!\left ( \|\nabla{\boldsymbol{Q}}^{(1)}_h - \nabla{\boldsymbol{Q}}^{(2)}_{h}\|_{H^1_tH^2_{{\boldsymbol{x}}}} + \|\nabla{\boldsymbol{Q}}^{(1)}_h - \nabla{\boldsymbol{Q}}^{(2)}_h\|_{H^2_tL^2_{{\boldsymbol{x}}}} \right ) \nonumber \\[5pt] & \qquad + C(T)\!\left ( \|{\boldsymbol{u}}^{(1)} -{\boldsymbol{u}}^{(2)}\|_{H^1_tH^2_{{\boldsymbol{x}}}} + \|{\boldsymbol{u}}^{(1)} -{\boldsymbol{u}}^{(2)}\|_{H^2_tL^2_{{\boldsymbol{x}}}} \right )\!\left ( \|\nabla{\boldsymbol{Q}}^{(2)}\|_{H^1_tH^2_{{\boldsymbol{x}}}} + \|\nabla{\boldsymbol{Q}}^{(2)}\|_{H^2_tL^2_{{\boldsymbol{x}}}}\right ) \nonumber \\[5pt] & \quad \leq C(T)\!\left ( \|{\boldsymbol{u}}^{(1)}\|_{H^1_tH^2_{{\boldsymbol{x}}}} + \|{\boldsymbol{u}}^{(1)}\|_{H^2_tL^2_{{\boldsymbol{x}}}} \right ) \|{\boldsymbol{Q}}^{(1)}_h -{\boldsymbol{Q}}^{(2)}_{h}\|_{X_{\boldsymbol{Q}}} \nonumber \\[5pt] & \qquad + C(T)\!\left ( \|{\boldsymbol{u}}^{(1)} -{\boldsymbol{u}}^{(2)}\|_{H^1_tH^2_{{\boldsymbol{x}}}} + \|{\boldsymbol{u}}^{(1)} -{\boldsymbol{u}}^{(2)}\|_{H^2_tL^2_{{\boldsymbol{x}}}} \right )\!\left ( \|{\boldsymbol{Q}}^{(2)}\|_{H^1_tH^3_{{\boldsymbol{x}}}} + \|{\boldsymbol{Q}}^{(2)}\|_{H^2_tH^1_{{\boldsymbol{x}}}}\right ). \end{align}

Using (157), (158) and (138) one gets

(168) \begin{align} & C(T)\!\left ( \|{\boldsymbol{u}}^{(1)}\|_{H^1_tH^2_{{\boldsymbol{x}}}} + \|{\boldsymbol{u}}^{(1)}\|_{H^2_tL^2_{{\boldsymbol{x}}}} \right ) \|{\boldsymbol{Q}}^{(1)}_h -{\boldsymbol{Q}}^{(2)}_{h}\|_{X_{\boldsymbol{Q}}} \nonumber \\[5pt] & \qquad + C(T)\!\left ( \|{\boldsymbol{u}}^{(1)} -{\boldsymbol{u}}^{(2)}\|_{H^1_tH^2_{{\boldsymbol{x}}}} + \|{\boldsymbol{u}}^{(1)} -{\boldsymbol{u}}^{(2)}\|_{H^2_tL^2_{{\boldsymbol{x}}}} \right )\!\left ( \|{\boldsymbol{Q}}^{(2)}\|_{H^1_tH^3_{{\boldsymbol{x}}}} + \|{\boldsymbol{Q}}^{(2)}\|_{H^2_tH^1_{{\boldsymbol{x}}}}\right ) \nonumber \\[5pt] & \quad \leq C(T)\!\left ( \|{\boldsymbol{u}}^{(1)}_h\|_{X_{\boldsymbol{u}}} + C(\|{\boldsymbol{\omega }}^{(1)}\|_{H^2(0,T)}+1) \right ) \!\left (\|{\boldsymbol{Q}}^{(1)}_h -{\boldsymbol{Q}}^{(2)}_{h}\|_{X_{\boldsymbol{Q}}} \right )\nonumber \\[5pt] & \qquad + C(T)\!\left (\|{\boldsymbol{u}}^{(1)}_h -{\boldsymbol{u}}^{(2)}_h\|_{X_{\boldsymbol{u}}} + C\!\left (\|\boldsymbol{\omega }^{(1)} - \boldsymbol{\omega }^{(2)}\|_{H^2(0,T)}\right ) \right ) \!\left ( \|{\boldsymbol{Q}}^{(2)}\|_{H^1_tH^3_{{\boldsymbol{x}}}} + \|{\boldsymbol{Q}}^{(2)}\|_{H^2_tH^1_{{\boldsymbol{x}}}}\right )\nonumber \\[5pt] & \quad \leq C(T) (R + 1) \|({\boldsymbol{u}}^{(1)}_h,{\boldsymbol{Q}}^{(1)}_h, \boldsymbol{\omega }^{(1)},{\boldsymbol{V}}^{(1)}) - ({\boldsymbol{u}}^{(2)}_h,{\boldsymbol{Q}}^{(2)}_h, \boldsymbol{\omega }^{(2)},{\boldsymbol{V}}^{(2)})\|_{X}. \end{align}

Thus, we obtained

\begin{eqnarray*} &&\|{\boldsymbol{u}}^{(1)} \cdot \nabla{\boldsymbol{Q}}^{(1)} -{\boldsymbol{u}}^{(2)} \cdot \nabla{\boldsymbol{Q}}^{(2)} \|_{H^1_tH^1_{{\boldsymbol{x}}}} \\[5pt] && \qquad \leq C(T,R)\|({\boldsymbol{u}}^{(1)}_h,{\boldsymbol{Q}}^{(1)}_h, \boldsymbol{\omega }^{(1)},{\boldsymbol{V}}^{(1)}) - ({\boldsymbol{u}}^{(2)}_h,{\boldsymbol{Q}}^{(2)}_h, \boldsymbol{\omega }^{(2)},{\boldsymbol{V}}^{(2)})\|_{X}. \end{eqnarray*}

Part 3: $\hat{{\bf H}}({\boldsymbol{Q}})$ . Here, we need to show

\begin{eqnarray*} \|{\bf H}({\boldsymbol{Q}}^{(1)}) -{\bf H}({\boldsymbol{Q}}^{(2)}) \| _{H^1_tH^1_{{\boldsymbol{x}}}} \leq C(T,R)\|({\boldsymbol{u}}^{(1)}_h,{\boldsymbol{Q}}^{(1)}_h, \boldsymbol{\omega }^{(1)},{\boldsymbol{V}}^{(1)}) - ({\boldsymbol{u}}^{(2)}_h,{\boldsymbol{Q}}^{(2)}_h, \boldsymbol{\omega }^{(2)},{\boldsymbol{V}}^{(2)})\|_{X}. \end{eqnarray*}

This bound follows directly from the proof on Part 3 of STEP 2.

Part 4: $\dfrac{\textrm{d}{\boldsymbol{V}}}{\textrm{d}t}$ .

(169) \begin{eqnarray} & \hspace{-40 pt}\|\dfrac{\textrm{d}{\boldsymbol{V}}^{(1)}}{\textrm{d}t} - \dfrac{\textrm{d}{\boldsymbol{V}}^{(2)}}{\textrm{d}t}\|_{H^1(0,T)}& \leq \|{\boldsymbol{V}}^{(1)} -{\boldsymbol{V}}^{(2)}\|_{H^2(0,T)}\nonumber \\[5pt] &&\leq C(T) (\|({\boldsymbol{u}}^{(1)}_h,{\boldsymbol{Q}}^{(1)}_h,{\boldsymbol{V}}^{(1)}, \boldsymbol{\omega }^{(1)}) - ({\boldsymbol{u}}^{(2)}_h,{\boldsymbol{Q}}^{(2)}_h,{\boldsymbol{V}}^{(2)}, \boldsymbol{\omega }^{(2)})\|_{X}). \end{eqnarray}

Part 5: $F_{\textrm{ext}}({\boldsymbol{Q}},{\boldsymbol{Q}}_{\infty })$ . We first note that for both definitions (11) and (12), $F_{\textrm{ext}}({\boldsymbol{Q}},{\boldsymbol{Q}}_{\infty })$ is a quadratic function of $\boldsymbol{Q}$ :

\begin{equation*} F_{\text {ext}}({\boldsymbol{Q}},{\boldsymbol{Q}}_{\infty })=\mathcal {B}({\boldsymbol{Q}},{\boldsymbol{Q}})=\sum \limits _{k,l,m,n=1}^{d}b_{klmn}Q_{kl}Q_{mn}, \end{equation*}

where coefficients $b_{klmn}$ depend on ${\boldsymbol{Q}}_{\infty }$ . Then use of the triangle inequality and (127):

\begin{align*} \|F_{\textrm{ext}}({\boldsymbol{Q}}^{(1)},{\boldsymbol{Q}}_{\infty })-F_{\textrm{ext}}({\boldsymbol{Q}}^{(2)},{\boldsymbol{Q}}_{\infty })&\|_{H^1_tH^1_{{\boldsymbol{x}}}}\leq \|\mathcal{B}({\boldsymbol{Q}}^{(1)}-{\boldsymbol{Q}}^{(2)},{\boldsymbol{Q}}^{(1)})\|_{H^1_tH^1_{{\boldsymbol{x}}}}\\[5pt] &\qquad\qquad\qquad\qquad+\|\mathcal{B}({\boldsymbol{Q}}^{(2)},{\boldsymbol{Q}}^{(1)}-{\boldsymbol{Q}}^{(2)})\|_{H^1_tH^1_{{\boldsymbol{x}}}}\\[5pt] &\leq C(T)\|{\boldsymbol{Q}}^{(1)}-{\boldsymbol{Q}}^{(2)}\|_{X_{\boldsymbol{Q}}}(\|Q^{(1)}\|_{H^1_tH^1_{{\boldsymbol{x}}}}+\|Q^{(2)}\|_{H^1_tH^1_{{\boldsymbol{x}}}})\\[5pt] &\leq C(T,R)\|{\boldsymbol{Q}}^{(1)}_h-{\boldsymbol{Q}}^{(2)}_h\|_{X_{\boldsymbol{Q}}}. \end{align*}

Part 6: $\frac{1}{I} \int _{\partial{\mathcal{P}}}{\boldsymbol{x}} \times \sigma \boldsymbol{\nu } + \ell \,dS_{\boldsymbol{x}}$ .

Since $\sigma = \sigma _{\textrm{hydro}} + \sigma _{\textrm{ela}}$ we use (140) and (159)

\begin{align*} &\| \sigma (\nabla{\boldsymbol{u}}^{(1)},{\boldsymbol{Q}}^{(1)}) - \sigma (\nabla{\boldsymbol{u}}^{(1)},{\boldsymbol{Q}}^{(1)})\|_{H^1_tH^{1}_{{\boldsymbol{x}}}} \\[5pt] & \hspace{2cm} \leq C(T, R) (\|({\boldsymbol{u}}^{(1)}_h,{\boldsymbol{Q}}^{(1)}_h,{\boldsymbol{V}}^{(1)}, \boldsymbol{\omega }^{(1)}) - ({\boldsymbol{u}}^{(2)}_h,{\boldsymbol{Q}}^{(2)}_h,{\boldsymbol{V}}^{(2)}, \boldsymbol{\omega }^{(2)})\|_{X}. \end{align*}

Using trace theorem, we get

(170) \begin{align} & \!\left \|\frac{1}{I} \int _{\partial{\mathcal{P}}}{\boldsymbol{x}} \times \sigma ({\boldsymbol{Q}}^{(1)}) \boldsymbol{\nu } dS_{\boldsymbol{x}} - \frac{1}{I} \int _{\partial{\mathcal{P}}}{\boldsymbol{x}} \times \sigma ({\boldsymbol{Q}}^{(2)}) \boldsymbol{\nu } dS_{\boldsymbol{x}} \right \|_{H^{1}(0,T)} \leq C\|\sigma ({\boldsymbol{Q}}^{(1)}) - \sigma ({\boldsymbol{Q}}^{(2)})\|_{H^1_tH^1_{{\boldsymbol{x}}}} \nonumber \\[5pt] &\hspace{100 pt}\leq C(T, R) \|({\boldsymbol{u}}^{(1)}_h,{\boldsymbol{Q}}^{(1)}_h,{\boldsymbol{V}}^{(1)}, \boldsymbol{\omega }^{(1)}) - ({\boldsymbol{u}}^{(2)}_h,{\boldsymbol{Q}}^{(2)}_h,{\boldsymbol{V}}^{(2)}, \boldsymbol{\omega }^{(2)})\|_{X}. \end{align}

To estimate the term with $\ell$ , recall its simplified form (19):

(171) \begin{eqnarray} \!\left \|\int \limits _{\partial \mathcal{P}}\ell ({\boldsymbol{Q}}^{(1)})\,\textrm{d}S_x-\int \limits _{\partial \mathcal{P}}\ell ({\boldsymbol{Q}}^{(2)})\,\textrm{d}S_x\right \|_{H^1(0,T)} &\leq& C \|{\boldsymbol{Q}}^{(1)}-{\boldsymbol{Q}}^{(2)}\|_{H^1(0,T;L^1(\partial \mathcal{P}))}\nonumber \\[5pt] &\leq& C\|{\boldsymbol{Q}}^{(1)}_h-{\boldsymbol{Q}}^{(2)}_h\|_{H^1_tH^1_{{\boldsymbol{x}}}}\nonumber \\[5pt] &\leq& C(T)\|{\boldsymbol{Q}}^{(1)}_h-{\boldsymbol{Q}}^{(2)}_h\|_{H^2_tH^1_{{\boldsymbol{x}}}}\nonumber \\[5pt] &\leq& C(T) \|({\boldsymbol{u}}^{(1)}_h,{\boldsymbol{Q}}^{(1)}_h,{\boldsymbol{V}}^{(1)}, \boldsymbol{\omega }^{(1)}) - ({\boldsymbol{u}}^{(2)}_h,{\boldsymbol{Q}}^{(2)}_h,{\boldsymbol{V}}^{(2)}, \boldsymbol{\omega }^{(2)})\|_{X}. \end{eqnarray}

Part 7: $\frac{1}{m} \int _{\partial{\mathcal{P}}} \sigma \boldsymbol{\nu } dS_{\boldsymbol{x}}$ . The same argument for Part 6 also works for $\frac{1}{m} \int _{\partial{\mathcal{P}}} \sigma \boldsymbol{\nu } dS_{\boldsymbol{x}}$ .

Part 8: $\partial _t{\boldsymbol{u}}_{\textrm{os}}$ . Using (135), we have

(172) \begin{align} \|\partial _t{\boldsymbol{u}}_{\textrm{os}}^{(1)} - \partial _t{\boldsymbol{u}}_{\textrm{os}}^{(2)}\|_{H^{1}_tL^{2}_{\sigma,{\boldsymbol{x}}}} \leq C(T) \|\boldsymbol{\omega }^{(1)} - \boldsymbol{\omega }^{(2)}\|_{H^{2}(0,T)}. \end{align}

Now, collecting all bounds from STEPS 2–4, we have (128) and thus Lemma 5.1 is proved.

5.3. Proof of Theorem 3.2 (local-in-time existence)

In this section, we prove the well-posedness of the time-dependent problem. The equation (101) can be rewritten as $\mathcal{K}\mathcal{U}_h=\mathcal{U}_h$ where $\mathcal{K}\;:\!=\;\mathcal{L}^{-1}\mathcal{J}\;:\;X\to X$ and $\mathcal{U}_h=({\boldsymbol{u}}_h,{\boldsymbol{Q}}_h,{\boldsymbol{\omega }},\boldsymbol{V})$ . The inverse linear operator $\mathcal{L}^{-1}$ is bounded, as stated in the following proposition.

Proposition 5.2. For all $(\boldsymbol{f}_{\boldsymbol{u}}, \boldsymbol{f}_{\boldsymbol{Q}}, \boldsymbol{f}_{{\boldsymbol{\omega }}}, \boldsymbol{f}_{\boldsymbol{V}}) \in Y$ , and time $T \in (0, 1]$ , linear system

(173) \begin{eqnarray} \mathcal{L}({\boldsymbol{u}}_h,{\boldsymbol{Q}}_h, \boldsymbol{\omega },{\boldsymbol{V}}) = (\boldsymbol{f}_{\boldsymbol{u}}, \boldsymbol{f}_{\boldsymbol{Q}}, \boldsymbol{f}_{{\boldsymbol{\omega }}}, \boldsymbol{f}_{\boldsymbol{V}}) \end{eqnarray}

has a unique solution such that ${\boldsymbol{u}}_h|_{t=0}=0$ , ${\boldsymbol{Q}}_h|_{t=0}=0$ , ${\boldsymbol{\omega }}|_{t=0}=\boldsymbol{V}|_{t=0}=0$ and

(174) \begin{eqnarray} \|({\boldsymbol{u}}_h,{\boldsymbol{Q}}_h, \boldsymbol{\omega },{\boldsymbol{V}})\|_{X} \leq C\|(\boldsymbol{f}_{\boldsymbol{u}}, \boldsymbol{f}_{\boldsymbol{Q}}, \boldsymbol{f}_{{\boldsymbol{\omega }}}, \boldsymbol{f}_{\boldsymbol{V}})\|_{Y}, \end{eqnarray}

where the constant $C$ is independent of time $T$ and choice of $(\boldsymbol{f}_{\boldsymbol{u}}, \boldsymbol{f}_{\boldsymbol{Q}}, \boldsymbol{f}_{{\boldsymbol{\omega }}}, \boldsymbol{f}_{\boldsymbol{V}})$ .

To prove this proposition, one can follow [47]. Specifically, for the first two components, ${\boldsymbol{u}}_h$ and ${\boldsymbol{Q}}_h$ , of system (173), we adapt the proof from [47, Proposition 4.2]. For the last two components, which are not present in [47], the statement naturally follows from the classical ordinary differential equation theory.

Next, according to Propositions 5.1 and 5.2, we have that

(175) \begin{eqnarray} &&\|\mathcal{K}({\boldsymbol{u}}_h^{(1)},{\boldsymbol{Q}}_h^{(1)}, \boldsymbol{\omega }^{(1)},{\boldsymbol{V}}^{(1)}) - \mathcal{K}({\boldsymbol{u}}_h^{(2)},{\boldsymbol{Q}}_h^{(2)}, \boldsymbol{\omega }^{(2)},{\boldsymbol{V}}^{(2)})\|_{X} \nonumber \\[5pt] &&\qquad \leq C(T) \|({\boldsymbol{u}}_h^{(1)},{\boldsymbol{Q}}_h^{(1)}, \boldsymbol{\omega }^{(1)},{\boldsymbol{V}}^{(1)}) - ({\boldsymbol{u}}_h^{(2)},{\boldsymbol{Q}}_h^{(2)}, \boldsymbol{\omega }^{(2)},{\boldsymbol{V}}^{(2)})\|_{X}. \end{eqnarray}

Recall that $C(T)$ depends on $T$ in such a way that $C(T) \rightarrow 0$ when $T \rightarrow 0$ . Choose $T$ such that $C(T) \lt 1$ . Then using Banach’s fixed point theorem, we obtain that there exists a unique fixed point $({\boldsymbol{u}}_h,{\boldsymbol{Q}}_h, \boldsymbol{\omega },{\boldsymbol{V}})$ of operator $\mathcal{K}$ . Next, define ${\boldsymbol{u}}_{\textrm{os}}$ via (95)–(97). Finally, the tuple $({\boldsymbol{u}}_h +{\boldsymbol{u}}_{\textrm{os}},{\boldsymbol{Q}}_h +{\boldsymbol{Q}}_{\textrm{os}}, \boldsymbol{\omega },{\boldsymbol{V}})$ is a solution to the original system (91), (2)–(6) with force and torque balances (14), (43) for $0\leq t\leq T$ . Theorem 3.2 is proved.

6. Homogenisation: two-scale expansion

In this section, we will perform formal two-scale expansion for (44)–(47). To this end, we introduce fast variable ${\boldsymbol{y}}=\varepsilon ^{-1}{\boldsymbol{x}}$ and represent the unknowns as

(176) \begin{equation} \!\left \{\begin{array}{r@{\quad}c@{\quad}l}{\boldsymbol{u}}({\boldsymbol{x}};\;\varepsilon )&=& \bar{{\boldsymbol{u}}}({\boldsymbol{x}},{\boldsymbol{y}})={\boldsymbol{u}}^{(0)}({\boldsymbol{x}},{\boldsymbol{y}})+\varepsilon{\boldsymbol{u}}^{(1)}({\boldsymbol{x}},{\boldsymbol{y}})+\cdots \\[5pt] p({\boldsymbol{x}};\;\varepsilon )&=&\bar{p}({\boldsymbol{x}},{\boldsymbol{y}})= p^{(0)}({\boldsymbol{x}},{\boldsymbol{y}})+\varepsilon p^{(1)}({\boldsymbol{x}},{\boldsymbol{y}})+\cdots \\[5pt]{\boldsymbol{Q}}({\boldsymbol{x}};\;\varepsilon )&=&\bar{{\boldsymbol{Q}}}({\boldsymbol{x}},{\boldsymbol{y}})={\boldsymbol{Q}}^{(0)}({\boldsymbol{x}},{\boldsymbol{y}})+\varepsilon{\boldsymbol{Q}}^{(1)}({\boldsymbol{x}},{\boldsymbol{y}})+\cdots \end{array}\right. \end{equation}

We will frequently use following identities for $f({\boldsymbol{x}},\varepsilon )=\bar{f}({\boldsymbol{x}},{\boldsymbol{y}})$ with ${\boldsymbol{y}}=\varepsilon ^{-1}{\boldsymbol{x}}$ :

(177) \begin{align} \nabla f=\nabla _{{\boldsymbol{x}}} \bar{f}+\varepsilon ^{-1}\nabla _{{\boldsymbol{y}}}\bar{f}, \end{align}

(178) \begin{align} \Delta f=\Delta _{{\boldsymbol{x}}} \bar{f} +2\varepsilon ^{-1}\nabla _{\boldsymbol{y}}\cdot \nabla _{{\boldsymbol{x}}}\bar{f}+\varepsilon ^{-2}\Delta _{{\boldsymbol{y}}}\bar{f}. \end{align}

The derivation of the homogenised limit consists of the following steps.

STEP 1. Show that ${\boldsymbol{u}}^{(0)}=0$ . Substitute two-scale representations (176) for ${\boldsymbol{u}},{\boldsymbol{Q}},$ and $p$ into (46) and $\nabla \cdot{\boldsymbol{u}}=0$ . We get that at level $\varepsilon ^{-1}$ :

\begin{equation*} -\tilde {\eta }\Delta _{{\boldsymbol{y}}}{\boldsymbol{u}}^{(0)}+\nabla _{{\boldsymbol{y}}}p^{(0)}=0 \text { and }\nabla _{{\boldsymbol{y}}}\cdot {\boldsymbol{u}}^{(0)}=0, \end{equation*}

with the boundary condition ${\boldsymbol{u}}^{(0)}=0$ on $\partial \mathcal{P}_{\varepsilon }$ . Thus, we can conclude that ${\boldsymbol{u}}^{(0)}=0$ .

STEP 2. Find an equation for ${\boldsymbol{Q}}^{(0)}$ . In this step, we expand equations (44) and (45) in $\varepsilon$ . To this end, we write the weak formulation of these two equations for arbitrary test function $\Phi \in H^1(\Omega _\varepsilon ;\;\mathbb{R}^{d\times d})$ :

(179) \begin{align} - \gamma \varepsilon \int \limits _{\Omega _\varepsilon } \nabla{\boldsymbol{Q}} \cdot \nabla \Phi \,\textrm{d}{\boldsymbol{x}}&+\varepsilon \tilde{W}\int \limits _{\partial \mathcal{P}_\varepsilon }({\boldsymbol{Q}}_{\textrm{pref}}-{\boldsymbol{Q}})\;:\;\Phi \,\textrm{d}S_{{\boldsymbol{x}}}\nonumber \\[5pt] & + \int \limits _{\Omega _{\varepsilon }}\!\left [\tilde{a}{\boldsymbol{Q}}-\tilde{c}{\boldsymbol{Q}}\textrm{Tr}({\boldsymbol{Q}}^2)\right ]\;:\;\Phi \;\textrm{d}{\boldsymbol{x}}+ \int \limits _{\Omega _\varepsilon }S(\nabla \tilde{{\boldsymbol{u}}},{\boldsymbol{Q}})\;:\;\Phi \,\textrm{d}{\boldsymbol{x}}\nonumber \\[5pt] & - \int \limits _{\Omega _\varepsilon }(\tilde{{\boldsymbol{u}}}\cdot \nabla ){\boldsymbol{Q}}\;:\;\Phi \,\textrm{d}{\boldsymbol{x}}+\tilde{\zeta }\int \limits _{\Omega _\varepsilon }\tilde{F}_{\textrm{ext}}\;:\;\Phi \,\textrm{d}{\boldsymbol{x}}=\int \limits _{\Omega _\varepsilon }{\bf G}\;:\;\Phi \,\textrm{d}{\boldsymbol{x}}. \end{align}

Introduce $\Omega _1=\varepsilon ^{-1}\Omega _\varepsilon$ and $\mathcal{P}_1=\varepsilon ^{-1}\mathcal{P}_\varepsilon$ . We now consider two-scale representation for the test function $\Phi$ :

\begin{equation*} \Phi ({\boldsymbol{x}};\;\varepsilon )=\bar {\Phi }({\boldsymbol{x}},{\boldsymbol{y}}) =\Phi ^{(0)}({\boldsymbol{x}}, {\boldsymbol{y}}) + \varepsilon \Phi ^{(1)}({\boldsymbol{x}}, {\boldsymbol{y}}) + \cdots. \end{equation*}

Rewrite the first integral in (179) in domain $\Omega _1$ :

(180) \begin{eqnarray} &&\hspace{-25 pt}-\gamma \varepsilon \int \limits _{\Omega _\varepsilon }\nabla{\boldsymbol{Q}} \cdot \nabla \Phi \,\textrm{d}{\boldsymbol{x}} = -\gamma \varepsilon ^{1+d}\int \limits _{\Omega _1}\nabla{\boldsymbol{Q}}(\varepsilon{\boldsymbol{y}})\cdot \nabla \Phi (\varepsilon{\boldsymbol{y}}) \,\textrm{d}{\boldsymbol{y}}\nonumber \\[5pt] &&\hspace{10 pt}= -\gamma \varepsilon ^{1+d}\int \limits _{\Omega _1}\!\left [\nabla _{{\boldsymbol{x}}}\bar{{\boldsymbol{Q}}}(\varepsilon{\boldsymbol{y}},{\boldsymbol{y}})+\varepsilon ^{-1}\nabla _{{\boldsymbol{y}}} \bar{{\boldsymbol{Q}}}(\varepsilon{\boldsymbol{y}},{\boldsymbol{y}})\right ]\cdot \left [\nabla _{{\boldsymbol{x}}} \bar{\Phi }(\varepsilon{\boldsymbol{y}},{\boldsymbol{y}})+\varepsilon ^{-1}\nabla _{{\boldsymbol{y}}} \bar{\Phi }(\varepsilon{\boldsymbol{y}},{\boldsymbol{y}})\right ] \,\textrm{d}{\boldsymbol{y}}\nonumber \\[5pt] && \hspace{10 pt}= -\Gamma K \varepsilon ^{1+d}\!\left \{\varepsilon ^{-2}\int \limits _{\Omega _1} \nabla _{{\boldsymbol{y}}}{\boldsymbol{Q}}^{(0)}\cdot \nabla _{{\boldsymbol{y}}} \Phi ^{(0)}\,\textrm{d}{\boldsymbol{y}}\;+ \right .\nonumber \\[5pt] &&\hspace{20 pt}\!\left .+\;\varepsilon ^{-1}\int \limits _{\Omega _1}\!\left [\nabla _{{\boldsymbol{x}}}{\boldsymbol{Q}}^{(0)}+\nabla _{{\boldsymbol{y}}}{\boldsymbol{Q}}^{(1)}\right ]\cdot \nabla \Phi ^{(0)}+\nabla{\boldsymbol{Q}}^{(0)}\cdot \left [\nabla _{{\boldsymbol{x}}}\Phi ^{(0)}+\nabla _{{\boldsymbol{y}}}\Phi ^{(1)}\right ]\!\textrm{d}{\boldsymbol{y}} +\cdots \right \}\hspace{-3 pt}. \end{eqnarray}

Expanding analogously other terms in (179) and using that ${\boldsymbol{u}}^{(0)}=0$ , we get at level $\varepsilon ^{d-1}$ :

\begin{equation*} \int \limits _{{\Omega _1}} \nabla _{\boldsymbol{y}}{\boldsymbol{Q}}^{(0)}\cdot \nabla _{\boldsymbol{y}}\Phi ^{(0)}\, \textrm{d}{\boldsymbol{y}} = 0, \end{equation*}

which implies, together with periodicity in $\boldsymbol{y}$ , that ${\boldsymbol{Q}}^{(0)}({\boldsymbol{x}},{\boldsymbol{y}}) ={\boldsymbol{Q}}^{(0)}({\boldsymbol{x}})$ .

At level $\varepsilon ^{d}$ , accounting for $\nabla _{{\boldsymbol{y}}}{\boldsymbol{Q}}^{(0)}=0$ , we have

(181) \begin{eqnarray} &&-\gamma \int _{{\Omega _1}} \!\left [\nabla _{\boldsymbol{x}}{\boldsymbol{Q}}^{(0)} + \nabla _{\boldsymbol{y}}{\boldsymbol{Q}}^{(1)}\right ]\cdot \nabla _{\boldsymbol{y}} \Phi ^{(0)}\, \textrm{d}{\boldsymbol{y}} + \tilde{W} \int _{\partial{{\mathcal{P}}_1}} ({\boldsymbol{Q}}_{\textrm{pref}} -{\boldsymbol{Q}}^{(0)})\;:\;\Phi _0\, \textrm{d}S_{\boldsymbol{y}} \nonumber \\[5pt] &&\hspace{2.5cm}+ \int \limits _{\Omega _1} \!\left [\!-\!\tilde{a}{\boldsymbol{Q}}^{(0)}+\tilde{c}{\boldsymbol{Q}}^{(0)}\textrm{Tr}(({\boldsymbol{Q}}^{(0)})^2)\right ]\;:\;\Phi ^{(0)}\, \textrm{d}{\boldsymbol{y}} \nonumber \\[5pt] && \hspace{2.5cm} + \int \limits _{\Omega _1} S(\nabla _{\boldsymbol{y}}{\boldsymbol{u}}^{(1)},{\boldsymbol{Q}}^{(0)})\cdot \Phi ^{(0)} \,\textrm{d}{\boldsymbol{y}} + \tilde{\zeta }\int \limits _{\Omega _1}\tilde{F}_{\textrm{ext}}({\boldsymbol{Q}}^{(0)},{\boldsymbol{Q}}_{\infty })\;:\;\Phi ^{(0)}\,\textrm{d}{\boldsymbol{y}} \nonumber \\[5pt] && \hspace{2cm} = \int \limits _{\Omega _1}{\bf G}({\boldsymbol{x}})\;:\; \Phi ^{(0)}\, \textrm{d}{\boldsymbol{y}}. \end{eqnarray}

Note that the above integral relation is the weak formulation for the following boundary-value problem:

(182) \begin{eqnarray} \!\left \{ \begin{array}{l@{\quad}l} -\gamma \Delta _{\boldsymbol{y}}{\boldsymbol{Q}}^{(1)} = \boldsymbol{f}_1, &{\boldsymbol{y}} \textrm{ in }{\Pi _1}, \\[5pt] \gamma (\nabla _{\boldsymbol{y}}{\boldsymbol{Q}}^{(1)}) \cdot \boldsymbol{\nu }_{\boldsymbol{y}} = \boldsymbol{g}_1, &{\boldsymbol{y}} \textrm{ on } \partial{\mathcal{P}}_{1}. \end{array} \right. \end{eqnarray}

Here,

\begin{eqnarray*} \boldsymbol{f}_1&=&\gamma \Delta _{{\boldsymbol{y}}{\boldsymbol{x}}}{\boldsymbol{Q}}^{(0)} - \tilde{a}{\boldsymbol{Q}}^{(0)} + \tilde{c}{\boldsymbol{Q}}^{(0)}\textrm{Tr}(({\boldsymbol{Q}}^{(0)})^2) + S(\nabla _{\boldsymbol{y}}{\boldsymbol{u}}^{(1)},{\boldsymbol{Q}}^{(0)})+\tilde{\zeta }\tilde{F}_{\textrm{ext}}({\boldsymbol{Q}}^{(0)},{\boldsymbol{Q}}_{\infty }) -{\bf G}({\boldsymbol{x}}),\\[5pt] \boldsymbol{g}_1&=&-\gamma \nabla _{\boldsymbol{x}}{\boldsymbol{Q}}^{(0)}\cdot \boldsymbol{\nu }_{\boldsymbol{y}} + \tilde{W}({\boldsymbol{Q}}_{\textrm{pref}} -{\boldsymbol{Q}}^{(0)}), \end{eqnarray*}

and $\Delta _{{\boldsymbol{y}}{\boldsymbol{x}}}\boldsymbol{h}=\nabla _{\boldsymbol{y}}\cdot \nabla _{{\boldsymbol{x}}}\boldsymbol{h}$ for arbitrary $\boldsymbol{h}$ .

Next, we have the solvability condition for (182) given as:

(183) \begin{equation} \int \limits _{\partial \mathcal{P}_1}\boldsymbol{g}_1\,\textrm{d}S_{\boldsymbol{y}}=\int \limits _{\Omega _1}\boldsymbol{f}_1\,\textrm{d}{\boldsymbol{y}}. \end{equation}

To evaluate the right-hand side, we use the fact that ${\boldsymbol{Q}}^{(0)}$ is independent of $\boldsymbol{y}$ and

(184) \begin{equation} \int \limits _{\Omega _1}\nabla _{{\boldsymbol{y}}}{\boldsymbol{u}}\,\textrm{d}{\boldsymbol{y}}=\int \limits _{\partial \mathcal{P}_1}u_{\textrm{sq}}\boldsymbol{\tau }\boldsymbol{\nu }^{\textrm{T}}\,\textrm{d}S_{\boldsymbol{y}}. \end{equation}

Hence, we have

(185) \begin{align} \int \limits _{\Omega _1}\boldsymbol{f}_1\,\textrm{d}{\boldsymbol{y}} &=|\Omega _1|\!\left ( - \tilde{a}{\boldsymbol{Q}}^{(0)} + \tilde{c}{\boldsymbol{Q}}^{(0)}\textrm{Tr}(({\boldsymbol{Q}}^{(0)})^2)+\tilde{\zeta }\tilde{F}_{\textrm{ext}}({\boldsymbol{Q}}^{(0)},{\boldsymbol{Q}}_{\infty }) -{\bf G}({\boldsymbol{x}})\right )\nonumber \\[5pt] &\quad +S\!\left(\int \limits _{\,\partial \mathcal{P}_1}u_{\textrm{sq}}\boldsymbol{\tau }\boldsymbol{\nu }^{\textrm{T}}\,\textrm{d}S_{\boldsymbol{y}},{\boldsymbol{Q}}^{(0)}\right), \end{align}

(186) \begin{align} \int \limits _{\partial \mathcal{P}_1}\boldsymbol{g}_1\,\textrm{d}S_{\boldsymbol{y}} &=\tilde{W}\int \limits _{\partial \mathcal{P}_1}{\boldsymbol{Q}}_{\textrm{pref}}\,\textrm{d}S_{\boldsymbol{y}}-\tilde{W}|\partial \mathcal{P}_1|{\boldsymbol{Q}}^{(0)}. \end{align}

Substituting (185) and (186) into (183), we get the equation for ${\boldsymbol{Q}}^{(0)}$ :

(187) \begin{align} -\left [\tilde{a}-\dfrac{\tilde{W}|\partial \mathcal{P}_1|}{|{\Omega _1}|}\right ]{\boldsymbol{Q}}^{(0)}&+\tilde{c}{\boldsymbol{Q}}^{(0)}\textrm{Tr}(({\boldsymbol{Q}}^{(0)})^2)\nonumber \\[5pt] &+ \tilde{\zeta }\tilde{F}_{\textrm{ext}}({\boldsymbol{Q}}^{(0)},{\boldsymbol{Q}}_{\infty })+S({\bf G}_{\textrm{sq}},{\boldsymbol{Q}}^{(0)})={\bf G}(x)-\dfrac{\tilde{W}}{|\Omega _1|}\overline{{\boldsymbol{Q}}}_{\textrm{pref}}. \end{align}

Here, we denote ${\bf G}_{\textrm{sq}}=\int \limits _{\partial \mathcal{P}_1}u_{\textrm{sq}}\boldsymbol{\tau }\boldsymbol{\nu }^{\textrm{T}}\,\textrm{d}S_{\boldsymbol{y}}$ and $\overline{{\boldsymbol{Q}}}_{\textrm{pref}}=\int \limits _{\partial \mathcal{P}_1}{\boldsymbol{Q}}_{\textrm{pref}}\,\textrm{d}S_{\boldsymbol{y}}$ . The function ${\boldsymbol{Q}}^{(0)}$ is the limit of $\boldsymbol{Q}$ as $\varepsilon \to 0$ ; thus, the algebraic equation (187) determines ${\boldsymbol{Q}}^{(h)}={\boldsymbol{Q}}^{(0)}$ .

STEP 3. Find an equation for ${\boldsymbol{u}}^{(1)}$ . At level $\varepsilon ^0$ in the expansion of (46), accounting for that ${\boldsymbol{u}}^{(0)}=0$ , ${\boldsymbol{Q}}^{(0)}$ is independent of $\boldsymbol{y}$ , and

\begin{eqnarray*} \varepsilon ^2\kappa \nabla \cdot (\nabla{{\boldsymbol{Q}}} \odot \nabla{{\boldsymbol{Q}}} +{{\boldsymbol{Q}}}\Delta{{\boldsymbol{Q}}} - \Delta{{\boldsymbol{Q}}}{{\boldsymbol{Q}}})=\varepsilon ^0\kappa \nabla _{{\boldsymbol{y}}}\cdot \left ({{\boldsymbol{Q}}}^{(0)}\Delta _{{\boldsymbol{y}}}{{\boldsymbol{Q}}}^{(1)} - \Delta _{{\boldsymbol{y}}}{{\boldsymbol{Q}}}^{(1)}{{\boldsymbol{Q}}}^{(0)}\right ) + \left [\begin{array}{c}\textrm{higher}\\[5pt]\textrm{order}\\[5pt]\textrm{terms}\end{array}\right ] \end{eqnarray*}

we get

(188) \begin{align} -\tilde{\eta }\Delta _{{\boldsymbol{y}}}{\boldsymbol{u}}^{(1)}+\nabla _{{\boldsymbol{x}}}p^{(0)} +\nabla _{{\boldsymbol{y}}}p^{(1)}=\kappa \nabla _{{\boldsymbol{y}}}\cdot \left ({{\boldsymbol{Q}}}^{(0)}\Delta _{{\boldsymbol{y}}}{{\boldsymbol{Q}}}^{(1)} - \Delta _{{\boldsymbol{y}}}{{\boldsymbol{Q}}}^{(1)}{{\boldsymbol{Q}}}^{(0)}\right )+{\bf F}({\boldsymbol{x}}), \,{\boldsymbol{y}}\textrm{ in }\Omega _1, \end{align}

(189) \begin{align} {\boldsymbol{u}}^{(1)}=u_{\textrm{sq}}\boldsymbol{\tau }, \,{\boldsymbol{y}}\textrm{ on }\partial \mathcal{P}_1. \end{align}

Next, we aim to remove ${\boldsymbol{Q}}^{(1)}$ from the right-hand side of (188). To this end, we notice that due to (182) we have that $-\gamma \Delta _{{\boldsymbol{y}}}{\boldsymbol{Q}}^{(1)}=\boldsymbol{f}_1$ and

\begin{eqnarray*} \Delta _{{\boldsymbol{y}}{\boldsymbol{x}}}{\boldsymbol{Q}}^{(0)}=0\textrm{ and }{\boldsymbol{Q}}^{(0)},\,{\boldsymbol{Q}}^{(0)}\textrm{Tr}({\boldsymbol{Q}}^{(0)}),\,\tilde{\zeta }\tilde{F}_{\textrm{ext}}({\boldsymbol{Q}}^{(0)},{\boldsymbol{Q}}_{\infty }),\,{\bf G}({\boldsymbol{x}})\textrm{ are independent of }{\boldsymbol{y}}. \end{eqnarray*}

Thus, we can rewrite (188) as

(190) \begin{eqnarray} &&\hspace{-30 pt}-\tilde{\eta }\Delta _{{\boldsymbol{y}}}{\boldsymbol{u}}^{(1)}+\nabla _{{\boldsymbol{x}}}p^{(0)} +\nabla _{{\boldsymbol{y}}}p^{(1)}\nonumber \\[5pt] && \hspace{20 pt}=-\kappa \gamma ^{-1} \nabla _{{\boldsymbol{y}}}\cdot \left ({{\boldsymbol{Q}}}^{(0)}S(\nabla _{\boldsymbol{y}}{\boldsymbol{u}}^{(1)},{\boldsymbol{Q}}^{(0)}) - S(\nabla _{\boldsymbol{y}}{\boldsymbol{u}}^{(1)},{\boldsymbol{Q}}^{(0)}){{\boldsymbol{Q}}}^{(0)}\right )+{\bf F}({\boldsymbol{x}}). \end{eqnarray}

The above can be rewritten in component-wise form as

(191) \begin{equation} \sum \limits _{m,j,l=1}^{d}\eta _{klmj}u_{m,jl}^{(1)}+\partial _{x_k}p^{(0)}+\partial _{y_k}p^{(1)}=F_{k}({\boldsymbol{x}}), \quad k=1,..,d, \end{equation}

where

(192) \begin{eqnarray} u_{m,jl}^{(1)} &=&\dfrac{\partial ^{2}u_{m}^{(1)}}{\partial y_{j}\partial y_{l}}\nonumber \\[5pt] \textrm{and}\quad \eta _{klmj}&=&-\tilde{\eta }\delta _{km}\delta _{jl}+ \dfrac{\kappa \xi }{2\gamma }\sum \limits _{n=1}^{d}\!\left [Q_{kn}^{(0)}Q_{nm}^{(0)}\delta _{jl}-Q_{mn}^{(0)}Q_{nl}^{(0)}\delta _{jk}-Q_{jn}^{(0)}Q_{nl}^{(0)}\delta _{km}\right ]\nonumber \\[5pt] &&+\dfrac{\kappa \xi }{d\gamma }\!\left [Q_{km}^{(0)}\delta _{jl}-Q_{ml}^{(0)}\delta _{jk}-Q_{jl}^{(0)}\delta _{km}\right ]+\dfrac{\kappa }{\gamma }\!\left [Q_{kj}^{(0)}Q_{ml}^{(0)}-Q_{km}^{(0)}Q_{jl}^{(0)}\right ]\nonumber \\[5pt] &&+\dfrac{\kappa }{\gamma }\sum \limits _{n=1}^{d}\!\left [Q_{kn}^{(0)}Q_{nm}^{(0)}\delta _{jl}-Q_{mn}^{(0)}Q_{nl}^{(0)}\delta _{jk}+Q_{jn}^{(0)}Q_{nl}^{(0)}\delta _{km}\right ]. \end{eqnarray}

Next, rewrite (191) in a vectorial form:

(193) \begin{equation} -\eta ({\boldsymbol{Q}}^{(0)})\nabla _{\boldsymbol{y}}^2{\boldsymbol{u}}^{(1)} + \nabla _{\boldsymbol{y}} p^{(1)} ={\bf F}({\boldsymbol{x}}) - \nabla _{\boldsymbol{x}} p^{(0)}. \end{equation}

Taking into account the boundary condition ${\boldsymbol{u}}^{(1)}({\boldsymbol{x}},{\boldsymbol{y}})=\tilde{u}_{\textrm{sq}}\boldsymbol{\tau }$ on $\mathcal{P}_1$ , we obtain the following representation for ${\boldsymbol{u}}^{(1)}$ :

(194) \begin{equation}{\boldsymbol{u}}^{(1)}=\mathcal{A}_{\eta ({\boldsymbol{Q}}^{(0)})}({\boldsymbol{y}})\!\left [{\bf F}({\boldsymbol{x}})-\nabla _{{\boldsymbol{x}}}p^{(0)}\right ]+\overline{{\boldsymbol{u}}}_{\textrm{sq}}, \end{equation}

where $\mathcal{A}_{\eta ({\boldsymbol{Q}}^{(0)})}({\boldsymbol{y}})$ is a $\boldsymbol{y}$ -dependent $d\times d$ matrix such that ${\boldsymbol{u}}({\boldsymbol{y}})=\mathcal{A}_{\eta ({\boldsymbol{Q}}^{(0)})}({\boldsymbol{y}})\boldsymbol{e}_i$ ( $\boldsymbol{e}_i$ is the $i$ th basis vector) is the solution of the following cell problem:

(195) \begin{equation} \!\left \{ \begin{array}{l} -\eta ({\boldsymbol{Q}}^{(0)})\nabla _{\boldsymbol{y}}^2{\boldsymbol{u}} + \nabla _{\boldsymbol{y}} p = \boldsymbol{e}_i, \textrm{ in } \Omega _1,\\[5pt] \nabla \cdot{\boldsymbol{u}}=0, \\[5pt]{\boldsymbol{u}} = 0, \textrm{ on }\mathcal{P}_1, \\[5pt]{\boldsymbol{u}} \textrm{ is 2-periodic}. \end{array} \right. \end{equation}

The term $\overline{{\boldsymbol{u}}}_{\textrm{sq}}$ is defined as the solution of

(196) \begin{equation} \!\left \{ \begin{array}{l} -\eta ({\boldsymbol{Q}}^{(0)})\nabla _{\boldsymbol{y}}^2 \overline{{\boldsymbol{u}}}_{\textrm{sq}} + \nabla _{\boldsymbol{y}} p = 0, \textrm{ in } \Omega _1,\\[5pt] \nabla \cdot \overline{{\boldsymbol{u}}}_{\textrm{sq}}=0, \\[5pt] \overline{{\boldsymbol{u}}}_{\textrm{sq}} = \tilde{u}_{\textrm{sq}}\boldsymbol{\tau }, \textrm{ on }\mathcal{P}_1, \\[5pt] \overline{{\boldsymbol{u}}}_{\textrm{sq}} \textrm{ is 2-periodic}. \end{array} \right. \end{equation}

Finally, we define the homogenised function ${\boldsymbol{u}}^{(h)}$ by averaging ${\boldsymbol{u}}^{(1)}$ and using the fact that ${\boldsymbol{Q}}^{(h)}={\boldsymbol{Q}}^{(0)}$ :

(197) \begin{eqnarray}{\boldsymbol{u}}^{(h)} = \mathcal{B}_{\eta ({\boldsymbol{Q}}^{(h)})} \!\left [{\bf F}({{\boldsymbol{x}}}) - \nabla _{\boldsymbol{x}} p\right ] + \frac{1}{|\Omega_1|}\int _{{\Omega }_1} \bar{{\boldsymbol{u}}}_{\textrm{sq}}\,\textrm{d}{\boldsymbol{y}}, \end{eqnarray}

where $\mathcal{B}_{\eta } =\frac{1}{|\Omega_1|}\int _{\Omega _1} \mathcal{A}_{\eta }(y)\, \textrm{d}{\boldsymbol{y}}$ and the pressure $p$ can be found from the divergence-free condition $\nabla \cdot{\boldsymbol{u}}=0$ .

To conclude, we have derived a system of homogenised equations in the form of (187), an algebraic equation for $\boldsymbol{Q}$ and (197) for $\boldsymbol{u}$ in the form of Darcy’s law.

7. Concluding remarks

In this work, we initiated the theoretical justification of the active microswimmer model, also known as a squirmer, in a liquid crystal. This model has been recently developed to explore a non-trivial response of the microswimmer to surrounding environment possessing a liquid crystalline structure. As computational studies [48, 49] clearly show that the squirmer eventually converges to an equilibrium, both time-dependent solutions and steady states are important and were in the focus of our work. In investigating well-posedness of the corresponding equations, we started with combining techniques from the two fields, the squirmer in Newtonian fluids and the Beris–Edwards model of liquid crystal. However, such a combination is not straightforward. As explained in Remark 2.2, one of the main difficulties, besides that our model is complex and highly nonlinear, is that it is not dissipative: there is a permanent energy input (not necessarily constant) coming from the activity of the squirmer. It makes application of a priori energy bounds established for the Beris–Edwards model not possible here. Therefore, the considered model requires novel approaches for its analysis. For the steady-state problem, using suitable offset functions, we first proved the existence of a steady state for a truncated system via Galerkin approximations and careful energy bounds using specific properties of the Beris–Edwards system (see, e.g., (76)) and then extended the well-posedness from the truncated system to the original one using the $L^{\infty }$ result formulated in Lemma 4.1. For the time-dependent problem, in order to exploit the contraction mapping principle, we considered higher regularity solutions (instead of weak ones) which allowed us to obtain all the necessary bounds including the one for integrals where activity of squirmer enters as well as the force and torque balances for the squirmer (see, e.g., (170)(171)). Periodic settings, in which we considered our model, allowed us to pose a question of homogenisation limit which would be a model describing a colony of synchronously moving squirmers. We found a scaling which, on the one hand, is consistent with experimental data (see Appendix C) and, on the other hand, allows for a non-trivial two-scale expansions so that the homogenised limit takes the form of Darcy’s law perturbed by an algebraic expression for the liquid crystal order parameter (see equations (187) and (197)).

Natural extensions of our work are:

1. (i) Stability analysis of steady states. Namely, we would like to find conditions on parameters when a steady state corresponding to swimming either parallel or perpendicular to the liquid crystal is stable. This analytical result will be compared with the main observation from [49] on bifurcations with respect to anchoring strength parameter. It would be also important to show that there is no steady state other than corresponding to swimming parallel or perpendicular to the liquid crystal.

2. (ii) Force–velocity relation for steady swimming. Though in the squirmer’s frame and periodic conditions, squirmer’s velocity is not well-defined, we can consider the so-called superficial velocity [61] $\boldsymbol{V}=-\frac{1}{|\Omega|}\int _{\Omega }{\boldsymbol{u}}\,\textrm{d}{\boldsymbol{x}}$ , which can be understood as the velocity of the squirmer with respect to the surrounding flow, and show how it depends on propulsion force entering the problem via $u_{\textrm{sq}}$ . Specifically, given the profile of the active slip velocity $u_{\textrm{sq}}$ (with all other physical parameters fixed), what is the resulting velocity $\boldsymbol{V}$ ? This question is related to the evaluation of the squirmer efficiency in Stokes fluid as a function of $u_{\textrm{sq}}$ [38, 42].

3. (iii) Rigorous justification of the homogenisation limit. We plan to justify the two-scale limit formally derived in Section 6 and in more general stochastic settings based on techniques developed for Newtonian fluids [52-56].

4. (iv) Extension to more general boundary conditions. It is also reasonable, from the physical point of view, to describe self-propulsion by a prescribed tangential force, $\sigma \boldsymbol{\nu }\times \boldsymbol{\nu } = \boldsymbol{f}_{\textrm{prop}}$ , on the squirmer’s surface, like it was implemented in [62, 63], as an alternative to the prescribed tangential velocity as in the classical squirmer’s model [35]. We plan to extend our results to the case when we allow prescribed tangential stress on a part of surface and tangential slip (or no slip) on the remaining part.