Hostname: page-component-586b7cd67f-l7hp2 Total loading time: 0 Render date: 2024-11-22T06:32:24.397Z Has data issue: false hasContentIssue false

Dynamics and steady state of squirmer motion in liquid crystal

Published online by Cambridge University Press:  10 July 2023

Leonid Berlyand
Affiliation:
Department of Mathematics, Pennsylvania State University, University Park, PA 16802, USA
Hai Chi
Affiliation:
Department of Mathematics, Pennsylvania State University, University Park, PA 16802, USA
Mykhailo Potomkin*
Affiliation:
Department of Mathematics, University of California at Riverside, Riverside, CA 92521, USA
Nung Kwan Yip
Affiliation:
Department of Mathematics, Purdue University, West Lafayette, IN 47907, USA
*
*Corresponding author: Mykhailo Potomkin; Email: [email protected]
Rights & Permissions [Opens in a new window]

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.

Type
Papers
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited.
Copyright
© The Author(s), 2023. Published by Cambridge University Press

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 [Reference Abels, Dolzmann and Liu1], see also reviews [Reference Abels, Dolzmann and Liu2, Reference Alexander, Pooley and Yeomans3]. 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 [Reference Allaire4Reference Chi, Potomkin, Zhang, Berlyand and Aranson10] and reviews [Reference Chisholm, Legendre, Lauga and Khair11Reference Duerinckx and Gloria14]. 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 [Reference Duerinckx and Gloria15, Reference Duerinckx and Gloria16]. 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 [Reference Duerinckx and Gloria17, Reference Elgeti, Winkler and Gompper18]. 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 [Reference Elgeti, Winkler and Gompper19Reference Galdi and Silvestre24].

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., [Reference Aranson6, Reference Aranson7, Reference Genkin, Sokolov, Lavrentovich and Aranson25Reference Guo, Zhu, Liu, Bonnet and Veerapaneni30], 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 [Reference Haines, Aranson, Berlyand and Karpeev31Reference Haines, Sokolov, Aranson, Berlyand and Karpeev33] and the book [Reference Hernandez-Ortiz, Stoltz and Graham34].

One of the most well-established model of microswimmer is the so-called squirmer. The model was initially introduced in [Reference Hernandez-Ortiz, Underhill and Graham35] 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 [Reference Kanevsky, Shelley and Tornberg36Reference Lions42].

To describe a nematic liquid crystal, we use the well-established Beris–Edwards model [Reference López, Gachelin, Douarche, Auradou and Clément43], 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 [Reference Marchetti, Joanny and Ramaswamy44, Reference Michelin and Lauga45]. 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 [Reference Michelin and Lauga46, Reference Najafi and Golestanian47].

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 [Reference Paicu and Zarnescu48] to study orientation dynamics of the spherical squirmer with respect to the preferred orientation of the liquid crystal. In [Reference Paicu and Zarnescu49], 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 [Reference Parsonnet, Friedman and Vandersteen50, Reference Paxton, Kistler and Olmeda51]. We also refer to [Reference Potomkin, Kaiser, Berlyand and Aranson52Reference Sohr56] 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 [Reference Hernandez-Ortiz, Underhill and Graham35] (which is also used in the computational work [Reference Paicu and Zarnescu48, Reference Paicu and Zarnescu49])

\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 [Reference Paicu and Zarnescu49] 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 [Reference Hernandez-Ortiz, Underhill and Graham35].

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 [Reference Figueroa-Morales, Dominguez-Rubio, Ott and Aranson21] 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 [Reference Paicu and Zarnescu49]

(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 [Reference Paicu and Zarnescu48, Reference Paicu and Zarnescu49] 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 [Reference Paicu and Zarnescu49] 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 [Reference Najafi and Golestanian47] 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 [Reference Najafi and Golestanian47] 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 [Reference Sokolov and Aranson57]. 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 [Reference Turiv, Koizumi and Thijssen58, 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 [Reference Viney59, 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 [Reference Turiv, Koizumi and Thijssen58, 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 [Reference Wang60, 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 [Reference Najafi and Golestanian47]. Specifically, for the first two components, ${\boldsymbol{u}}_h$ and ${\boldsymbol{Q}}_h$ , of system (173), we adapt the proof from [Reference Najafi and Golestanian47, Proposition 4.2]. For the last two components, which are not present in [Reference Najafi and Golestanian47], 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 [Reference Paicu and Zarnescu48, Reference Paicu and Zarnescu49] 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 [Reference Paicu and Zarnescu49] 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 [Reference Wong, Dey and Sen61] $\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}}$ [Reference Lauga and Powers38, Reference Lions42].

  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 [Reference Potomkin, Kaiser, Berlyand and Aranson52-Reference Sohr56].

  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 [Reference Zhou62, Reference Zhou, Sokolov, Lavrentovich and Aranson63], as an alternative to the prescribed tangential velocity as in the classical squirmer’s model [Reference Hernandez-Ortiz, Underhill and Graham35]. 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.

Acknowledgements

Financial support

The work of L. Berlyand was supported by NSF grant DMS-2005262 and the work of H. Chi was partially supported by NSF grant DMS-2005262.

Competing interests

The authors declare none.

Appendix

A. Proof of Lemma 4.2

Proposition A.1 (A Poincaré-type estimate). There exists $c_P\gt 0$ such that

(198) \begin{equation} \|{\boldsymbol{Q}}\|^{2}_{L^2(\Omega )}\leq c_P \!\left (\dfrac{K}{2}\|\nabla{\boldsymbol{Q}}\|^2_{L^2(\Omega )} + \dfrac{W}{2}\|{\boldsymbol{Q}}\|^2_{L^2(\partial{\mathcal{P}}_{\textrm{st}})}\right ). \end{equation}

for all ${\boldsymbol{Q}}\in H^1_{\textrm{per}}(\Omega )$ .

Proof. Here, we first show that there exists $\hat{c}_P\gt 0$ such that for arbitrary ${\boldsymbol{Q}}\in H^1_{\textrm{per}}(\Omega )$ the following inequality holds:

(199) \begin{equation} \|{\boldsymbol{Q}}\|^2_{L^2(\Omega )}\leq \hat{c}_P(\|\nabla{\boldsymbol{Q}}\|_{L^2(\Omega )}^2+\|{\boldsymbol{Q}}\|^2_{L^2(\mathcal{P}_{\textrm{st}})}). \end{equation}

By contradiction, we assume that there exists a sequence $\{{\boldsymbol{Q}}_n\}_{n=1}^{\infty }$ :

(200) \begin{equation} \|{\boldsymbol{Q}}_n\|^2_{L^2(\Omega )}=1 \textrm{ and }\|\nabla{\boldsymbol{Q}}_n\|_{L^2(\Omega )}^2+\|{\boldsymbol{Q}}_{n}\|^2_{L^2(\mathcal{P}_{\textrm{st}})}\leq \dfrac{1}{n}. \end{equation}

From boundedness of ${\boldsymbol{Q}}_n$ in $L^2(\Omega )$ , it follows that the sequence $\{{\boldsymbol{Q}}_n\}_{n=1}^{\infty }$ possesses a weakly converging sub-sequence in $L^2(\Omega )$ . Consider any such weakly converging sub-sequence $\{{\boldsymbol{Q}}_{n_k}\}_{k=1}^{\infty }$ , ${\boldsymbol{Q}}_{n_k}\rightharpoonup{\boldsymbol{Q}}^{*}$ and a function $\psi \in H^1_{\textrm{per}}(\Omega )$ :

\begin{equation*} \int \limits _{\Omega }{\boldsymbol{Q}}_{n_k}\partial _{x_i}\psi \,\text {d}x= \int \limits _{\partial \mathcal {P}_{\text {st}}}{\boldsymbol{Q}}_{n_k}\psi \nu _i\,\text {d}S_x-\int \limits _{\Omega }\partial _{x_i}{\boldsymbol{Q}}_{n_k}\psi \,\text {d}x \to 0,\quad 1\leq i\leq d. \end{equation*}

The convergence to $0$ follows from strong convergences of ${\boldsymbol{Q}}_{n_k}$ in $L^2(\partial \mathcal{P}_{\textrm{st}})$ and $\nabla{\boldsymbol{Q}}_{n_k}$ in $L^2(\Omega )$ which in turn follow from (200). Then

\begin{equation*} \int \limits _{\Omega }{\boldsymbol{Q}}^{*}\partial _{x_i}\psi \,\text {d}x=0 \quad \forall \psi \in H^1_{\text {per}}(\Omega ). \end{equation*}

Using integration by parts, we get

\begin{equation*} \int \limits _{\partial \mathcal {P}_{\text {st}}}{\boldsymbol{Q}}^{*}\psi \nu _i\,\text {d}S_x-\int \limits _{\Omega }\partial _{x_i}{\boldsymbol{Q}}^{*}\psi \,\text {d}x=0\quad \forall \psi \in H^1_{\text {per}}(\Omega ). \end{equation*}

By taking first $\psi \in C^{\infty }_{c}(\Omega )$ , we get $\partial _{x_i}{\boldsymbol{Q}}^{*}\equiv 0$ so the second integral in the equality above vanishes. Next, we get that the first integral is zero as well by taking $\psi$ with various non-zero traces. We conclude ${\boldsymbol{Q}}^{*}\equiv 0$ . Moreover, since for any weakly converging sub-sequence the limit is $0$ , then entire sequence $\{{\boldsymbol{Q}}_n\}_{n=1}^{\infty }$ weakly converges to $0$ .

Note that $H^1_{\textrm{per}}(\Omega )\subset H^1(\Omega )$ and thus $H^1_{\textrm{per}}(\Omega )$ is compactly embedded in $L^2(\Omega )$ . Hence, ${\boldsymbol{Q}}_{n}$ is strongly converging in $L^2(\Omega )$ and since the weak limit is $0$ , we conclude ${\boldsymbol{Q}}_n\to 0$ strongly in $L^2(\Omega )$ . However, it contradicts to (200) since it implies that if ${\boldsymbol{Q}}_n\to{\boldsymbol{Q}}^*$ strongly in $L^2(\Omega )$ , then $\|{\boldsymbol{Q}}^*\|_{L^2(\Omega )}^2=1$ . Thus, inequality (199) is shown.

Finally, to prove (198) take $c_P=\dfrac{2\hat{c}_P}{\min \{K,W\}}$ .

Next we turn to the proof of Lemma 4.2. We note that an equivalent way to define ${\boldsymbol{Q}}_m$ from (63) is via minimisation of the following energy functional:

(201) \begin{equation} \mathcal{E}({\boldsymbol{Q}})=\int \limits _{\Omega }\dfrac{K}{2}|\nabla{\boldsymbol{Q}}|^2 + \hat{\mathcal{F}}_M({\boldsymbol{Q}})\,\textrm{d}x+\dfrac{W}{2}\int \limits _{\partial \mathcal{P}_{\textrm{st}}}|{\boldsymbol{Q}}_{\textrm{pref}}-{\boldsymbol{Q}}|^2 \,\textrm{d}S_x + \int \limits _{\Omega }{\bf H}_m\;:\;{\boldsymbol{Q}} \,\textrm{d}x \end{equation}

among functions ${\boldsymbol{Q}}\in H^1_{\textrm{per}}(\Omega )$ . The minimiser $\boldsymbol{Q}$ of the energy functional $\mathcal{E}({\boldsymbol{Q}})$ exists.

From (198), the Cauchy inequality, and $\|{\boldsymbol{Q}}-{\boldsymbol{Q}}_{\textrm{pref}}\|^2\geq \dfrac{1}{2}\|{\boldsymbol{Q}}\|^2-\|{\boldsymbol{Q}}_{\textrm{pref}}\|^2$ we obtain that $\mathcal{E}({\boldsymbol{Q}})$ is bounded from below:

\begin{eqnarray*} \mathcal{E}({\boldsymbol{Q}})&\geq & \dfrac{K}{2}\|\nabla{\boldsymbol{Q}}\|^2_{L^2(\Omega )} + \dfrac{W}{2}\|{\boldsymbol{Q}}_{\textrm{pref}}-{\boldsymbol{Q}}\|^2_{L^2(\partial \mathcal{P}_{\textrm{st}})}-\dfrac{1}{2c_P}\|{\boldsymbol{Q}}\|^2_{L^2(\Omega )}-\dfrac{c_P}{2}\|{\bf H}_m\|^2_{L^2(\Omega )}\\[5pt] &\geq & -\dfrac{W}{2}\|{\boldsymbol{Q}}_{\textrm{pref}}\|_{L^2(\partial \mathcal{P}_{\textrm{st}})}-\dfrac{c_P}{2}\|{\bf H}_m\|^2_{L^2(\Omega )}. \end{eqnarray*}

Thus, there exists a minimising sequence ${\boldsymbol{Q}}^{(\ell )}$ weakly converging in $H^1_{\textrm{per}}(\Omega )$ . Then using the $\liminf$ property of a weakly converging sequences, we get the existence of minimiser ${\boldsymbol{Q}}={\boldsymbol{Q}}_m$ .

From $\mathcal{E}({\boldsymbol{Q}}_m)\leq \mathcal{E}({\boldsymbol{Q}}_\infty )$ and (198), we get that there exists $C\gt 0$ such that

(202) \begin{eqnarray} \dfrac{K}{8}\|\nabla{\boldsymbol{Q}}_m\|^2_{L^2(\Omega )} +\dfrac{1}{8c_P}\|{\boldsymbol{Q}}_m\|_{L^2(\Omega )}^2+\dfrac{W}{4}\|{\boldsymbol{Q}}_{\textrm{pref}}-{\boldsymbol{Q}}_m\|^2_{L^2(\partial{\mathcal{P}}_{\textrm{st}})}\leq 3c_P\|{\bf H}_m\|^2_{L^2(\Omega )}+C. \end{eqnarray}

This shows (64). Next, due to the elliptic regularity result (see Appendix B)

\begin{equation*} \!\left \{ \begin{array}{l} \Delta{\boldsymbol{Q}}= K^{-1}({\bf H}_m-\hat{{\bf H}}_M({\boldsymbol{Q}}_m))\textrm{ in }\Omega, \\[5pt] \partial _\nu{\boldsymbol{Q}}|_{\partial{\mathcal{P}}_{\textrm{st}}}=WK^{-1}({\boldsymbol{Q}}_{\textrm{pref}}-{\boldsymbol{Q}}) \end{array} \right. \end{equation*}

we have

(203) \begin{equation} \quad \|{\boldsymbol{Q}}_m\|_{H^2(\Omega )}^2\leq C \!\left (\dfrac{(W+K)^2}{W^2}\|{\bf H}_m-\hat{{\bf H}}_M({\boldsymbol{Q}}_m)\|_{L^2(\Omega )}^2+\dfrac{W+K}{K}\|{\boldsymbol{Q}}_{\textrm{pref}}\|_{C^1}^2\right ). \end{equation}

Next, using (62) we get (65) and it completes the proof of Lemma 4.2.

B. Elliptic regularity for the squirmer boundary conditions

Here, we consider the following auxiliary elliptic problem:

(204) \begin{equation} \!\left \{ \begin{array}{l} \Delta q= F\textrm{ in }\Omega, \\[5pt] \partial _\nu q|_{\partial \mathcal{P}_{\textrm{st}}}=\beta (\gamma ({\boldsymbol{x}}) -q),\\[5pt] q \textrm{ is }\Pi \textrm{-periodic.} \end{array} \right. \end{equation}

Here, $F=F({\boldsymbol{x}})$ , $\gamma =\gamma ({\boldsymbol{x}})$ and $\beta \gt 0$ are given. We aim to prove the elliptic regularity for (204):

Theorem B.1. Let $q$ be the solution of (204). Then

(205) \begin{equation} \|q\|_{H^2(\Omega )}^2\leq C \!\left (\frac{(1+\beta )^2}{\beta ^2}\|F\|_{L^2(\Omega )}^2+ (1+\beta )\|\gamma \|_{C^1}^2\right ) \end{equation}

for some constant $C$ independent of $\beta$ and $\gamma$ .

This result is well known from PDE textbooks [Reference Zhou, Tovkach, Golovaty, Sokolov, Aranson and Lavrentovich64], Theorem 8.12] and [Reference Zick and Homsy65], Theorem 4 in Section 6.3.2] for the Dirichlet boundary conditions. However, we need to re-visit this result due to our specific boundary conditions for which the afore-mentioned results are not applicable. For the sake of clarity, the proof below is written for two-dimensional case, $d=2$ .

Proof. We first address a priori estimates for regions near the boundary $\partial \mathcal{P}_{\textrm{st}}$ of the squirmer. Choose any point ${\boldsymbol{x}}_0$ on $\partial \mathcal{P}_{\textrm{st}}$ . Suppose its vicinity on the boundary can be described by equation $x_2=\varphi (x_1)$ so the domain $x_2\gt \varphi (x_1)$ is the interior of domain $\Omega$ . Then introduce change of variables

(206) \begin{equation}{\boldsymbol{y}}=\Phi ({\boldsymbol{x}}) \!\Leftrightarrow \!\left \{ \begin{array}{l} y_1=x_1,\\[5pt] y_2=x_2-\varphi (x_1). \end{array}\right. \end{equation}

In variable $\boldsymbol{y}$ , the problem (204) has the form:

\begin{equation*} \!\left \{ \begin {array}{l} \nabla _{{\boldsymbol{y}}}\cdot (L({\boldsymbol{y}})\nabla _{{\boldsymbol{y}}}q)=F,\quad y_2\gt 0\\[5pt] L({\boldsymbol{y}})\nabla _{{\boldsymbol{y}}}q\cdot \nu _{\boldsymbol{y}}=\beta \sqrt {1+(\varphi ')^2}\,(\gamma - q),\quad y_2=0, \end {array} \right. \end{equation*}

where

\begin{equation*} L= \!\left [ \begin {array}{c@{\quad}c} 1 & -\varphi ' \\[5pt] -\varphi ' & (1+(\varphi ')^2) \end {array} \right ], \quad \nu _{\boldsymbol{y}}= \!\left [ \begin {array}{c} 0\\[5pt]1 \end {array} \right ]. \end{equation*}

We note that $L$ is a positive definite symmetric matrix with the smallest eigenvalue

(207) \begin{equation} \lambda _{\textrm{min}}({\boldsymbol{y}})=\dfrac{1}{2}\!\left (2+(\varphi ')^2-\sqrt{(2+(\varphi ')^2)^2-4}\right )\geq 1. \end{equation}

Thus, $L$ is uniformly positive definite:

(208) \begin{equation} \textrm{$(L({\boldsymbol{y}}){\boldsymbol{u}}\cdot{\boldsymbol{u}})\geq |{\boldsymbol{u}}|^2$ for all ${\boldsymbol{y}},{\boldsymbol{u}}\in \mathbb R^2$.} \end{equation}

Lemma B.1. Let $q$ be the solution of

(209) \begin{equation} \!\left \{ \begin{array}{l} \nabla _{\boldsymbol{y}} \cdot (L({\boldsymbol{y}})\nabla _{\boldsymbol{y}} q) = \hat{F}, \quad{\boldsymbol{y}}\in \mathbb R_+^2= \left \{(y_1,y_2)\;:\;y_2\gt 0\right \}\\[5pt] L({\boldsymbol{y}})\partial _{{\boldsymbol{y}}}q\cdot \nu _{\boldsymbol{y}}=\hat{\beta }(\hat{\gamma }- q), \quad y_2=0, \end{array} \right. \end{equation}

for some $f\in L^2_{loc}(\mathbb R_{+}^2)$ , $\hat{\gamma }= \hat{\gamma }(y_1)\in H^1(\mathbb R)$ and $\hat{\beta }(y_1)\geq \beta _0\gt 0$ and matrix $L({\boldsymbol{y}})$ satisfying uniform positivity condition (208). Denote also $U=B_1(0)\cap \!\left \{y_2\gt 0\right \}$ and $V=B_{1/2}(0)\cap \!\left \{y_2\gt 0\right \}$ .

Then we have the following bound:

(210) \begin{equation} \|q\|_{H^2(V)}^2\leq C\!\left (\|\hat{F}\|_{L^2(U)}^2+\|q\|_{H^1(U)}^2+\|q\|_{L^2(U_0)}^2+\beta ^2 \|\gamma \|_{C^1}^2\right ). \end{equation}

Here $U_0=\{y_2=0\}\cap U$ .

Proof. We adapt arguments from [Reference Zick and Homsy65]. All gradients $\nabla$ in the proof of this lemma are taken with respect to variable $\boldsymbol{y}$ . First, we write the weak formulation of (209) for all $v\in H^1(U)$ :

(211) \begin{equation} \int \limits _{\tilde{\Pi }}L\nabla q \cdot \nabla v \,\textrm{d}y-\int \limits _{\{y_2=0\}\cap \tilde{\Pi }}\hat{\beta }(\hat{\gamma } - q)v\,\textrm{d}S_y=-\int \limits _{\tilde{\Pi }} \hat{F}v\,\textrm{d}y. \end{equation}

Here, $\tilde{\Pi }$ is the image of $\Pi \setminus \mathcal{P}_{\textrm{st}}$ under transformation (206). Next, we introduce $\zeta ({\boldsymbol{y}})$ such that $\zeta \in C^{\infty }$ and $\zeta \equiv 1$ in $V$ and $\zeta \equiv 0$ outside of $W=B_{{3/4}}(0)\cap \!\left \{y_2\gt 0\right \}$ . Take test function $v=D_{1}^{-h}(\zeta ^2D_1^{h}q)$ , where $D_1^{h}$ is the difference quotient operator:

(212) \begin{equation} D_1^{h}g=\dfrac{g(y_1+h,y_2)-g(y_1,y_2)}{h}. \end{equation}

Integration by parts for the difference quotient allows us to rewrite the first integral in (211) as follows:

(213) \begin{eqnarray} \int \limits _{\tilde{\Pi }}L\nabla q \cdot \nabla v \,\textrm{d}y&=&-\int \limits _{U}\zeta ^2\,(LD^h_1(\nabla q)\cdot D^h_1(\nabla q)) \,\textrm{d}y-2\int \limits _{U}\zeta D_1^h q (\nabla \zeta \cdot LD^h_1(\nabla q))\,\textrm{d}y\nonumber \\[5pt] &&\nonumber -\int \limits _{U}\zeta ^2 (\!\left [D^h_1L\right ](\nabla q) \cdot D^h_1(\nabla q))\,\textrm{d}y -2\int \limits _{U}\zeta D_1^h q (\nabla \zeta \cdot \left [D^h_1L\right ](\nabla q))\,\textrm{d}y\nonumber \\[5pt] &&\leq -\dfrac{1}{2}\int \limits _{U}\zeta ^2|D^h_1(\nabla q)|^2 \,\textrm{d}y+C\int \limits _{W}|D^h_1q|^2\,\textrm{d}y+C\int \limits _{U}|\nabla q|^2\,\textrm{d}y. \end{eqnarray}

Here, to obtain the estimate we used (208), uniform boundedness of $\nabla \zeta$ and $[D^h_1L]$ and Cauchy–Schwarz inequality.

Similarly, we rewrite the second integral in (211):

(214) \begin{eqnarray} -\int \limits _{\{y_2=0\}\cap \tilde{\Pi }}\hat{\beta }(\hat{\gamma } - q)v\,\textrm{d}S_y &=& -\int \limits _{\mathbb R}\hat{\beta }(\hat{\gamma } - q)D_{1}^{-h}(\zeta ^2D_1^{h}q)\,\textrm{d}y_1\nonumber \\[5pt] &=&\int \limits _{\mathbb R}\,D_1^{h}(\hat{\beta }\hat{\gamma } - \hat{\beta } q)\,(D_1^{h}q) \zeta ^2\,\textrm{d}y_1\nonumber \\[5pt] &=&\int \limits _{\mathbb R}\,\zeta ^2 \!\left [D_1^{h}(\hat{\beta }\hat{\gamma })\right ] D_1^{h}q\,\textrm{d}y_1-\int \limits _{\mathbb R}\zeta ^2 q (D^h_1\hat{\beta }) D^h_1 q\,\textrm{d}y_1 -\int \limits _{\mathbb R}\hat{\beta }\zeta ^2|D^h_1q|^2\,\textrm{d}y_1 \nonumber \\[5pt] &\leq & -\dfrac{\beta _0}{2}\int \limits _{\mathbb R}\zeta ^2|D^h_1q|^2\,\textrm{d}y_1+C\int \limits _{\mathbb R}\zeta ^2 q^2\,\textrm{d}y_1+ C\int \limits _{\mathbb R}\zeta ^2|D_1^h(\hat{\gamma }\hat{\beta })|^2\,\textrm{d}y_1. \end{eqnarray}

Finally, we estimate the third integral in (211) using [Reference Zick and Homsy65], Theorem 3 from Section 5.8.2] as follows:

(215) \begin{eqnarray} -\int \limits _{\tilde{\Pi }}\hat{F}v\,\textrm{d}y&&\leq C\!\left (\int \limits _{U}|\hat{F}|^2 \textrm{d}y\right )^{1/2} \!\left (\int \limits _{U}|\nabla (\zeta ^2 D^h_1q)|^2\,\textrm{d}y\right )^{1/2}\nonumber \\[5pt]&&\leq C\int \limits _{U}|\hat{F}|^2 \,\textrm{d}y+\dfrac{1}{8}\int \limits _{U}|\nabla (\zeta ^2D_1^hq)|^2\,\textrm{d}y\nonumber \\[5pt] &&\leq C\int \limits _{U}|\hat{F}|^2 \,\textrm{d}y+\dfrac{1}{4}\int \limits _{U}\zeta ^2|\nabla (D_1^hq)|^2\,\textrm{d}y+\dfrac{1}{4}\int \limits _{U}|\nabla (\zeta )^2\cdot D_1^hq|^2\,\textrm{d}y\nonumber \\[5pt] &&\leq C\int \limits _{U}|\hat{F}|^2 \,\textrm{d}y+\dfrac{1}{4}\int \limits _{U}\zeta ^2|\nabla (D_1^hq)|^2\,\textrm{d}y+C\int \limits _{W}|D_1^hq|^2\,\textrm{d}y. \end{eqnarray}

Combining (211) with estimates (213), (214) and (215), we get

\begin{eqnarray*} &&\dfrac{1}{4}\int \limits _{U}\zeta ^2|D^h_1(\nabla q)|^2 \,\textrm{d}y+\dfrac{\beta _0}{2}\int \limits _{\mathbb R}\zeta ^2|D^h_1q|^2\,\textrm{d}y_1 \\[5pt] &&\hspace{60 pt} \leq C\!\left [\int \limits _{U}|\hat{F}|^2 \,\textrm{d}y+\int \limits _{W}|D^h_1q|^2\,\textrm{d}y+\int \limits _{U}|\nabla q|^2\,\textrm{d}y\right .\\[5pt] &&\hspace{120 pt}\!\left .+\int \limits _{\mathbb R}\zeta ^2 q^2\,\textrm{d}y_1+ \int \limits _{\mathbb R}\zeta ^2|D_1^h(\hat{\gamma }\hat{\beta })|^2\,\textrm{d}y_1\right ]. \end{eqnarray*}

Using [Reference Zick and Homsy65, Theorem 3(i) from Section 5.8.2] we get

\begin{eqnarray*} &&\dfrac{1}{4}\int \limits _{V}|D^h_1(\nabla q)|^2 \,\textrm{d}y+\dfrac{\beta _0}{2}\int \limits _{\mathbb R \cap V}|D^h_1 q|^2\,\textrm{d}y_1 \\[5pt] &&\hspace{60 pt} \leq C\!\left [\int \limits _{U}|\hat{F}|^2 \,\textrm{d}y+\int \limits _{U}|\nabla q|^2\,\textrm{d}y+\int \limits _{\mathbb R}\zeta ^2 q^2\,\textrm{d}y_1+ \int \limits _{\mathbb R}\zeta ^2|\partial _{y_1} (\hat{\gamma }\hat{\beta })|^2\,\textrm{d}y_1\right ]. \end{eqnarray*}

Then [Reference Zick and Homsy65, Theorem 3(ii) from Section 5.8.2] implies

(216) \begin{eqnarray} &&\dfrac{1}{4}\int \limits _{V}|\partial _{y_1}(\nabla q)|^2 \,\textrm{d}y+\dfrac{\beta _0}{2}\int \limits _{\mathbb R \cap V}|\partial _{y_1} q|^2\,\textrm{d}y_1 \nonumber \\[5pt] &&\hspace{60 pt} \leq C\!\left [\int \limits _{U}|\hat{F}|^2 \,\textrm{d}y+\int \limits _{U}|\nabla q|^2\,\textrm{d}y+\int \limits _{\mathbb R}\zeta ^2 q^2\,\textrm{d}y_1+ \int \limits _{\mathbb R}\zeta ^2|\partial _{y_1} (\hat{\gamma }\hat{\beta })|^2\,\textrm{d}y_1\right ]. \end{eqnarray}

Now write $\partial _{y_1} (\hat{\gamma }\hat{\beta })=\beta \!\left (\dfrac{\varphi ''\varphi '}{\sqrt{1+(\varphi ')^2}}\gamma +\sqrt{1+(\varphi ')^2}(\gamma _{x_1}+\gamma _{x_2}\varphi ')\right )$ . Thus,

\begin{equation*} \int \limits _{\mathbb R}\zeta ^2|\partial _{y_1} (\hat {\gamma }\hat {\beta })|^2\,\text {d}y_1\leq C\beta ^2 \|\gamma \|_{C^1}^2. \end{equation*}

Analogous estimate is valid for $\partial ^2_{x_2}q$ since

\begin{equation*} \partial ^2_{y_2}q=-\dfrac {1}{1+(\varphi ')^2}\!\left [\hat {F}-\partial _{y_1}^2q+2\varphi '\partial ^2_{y_1y_2}q-\varphi ''\partial _{y_2}q+2\varphi '\varphi ''\partial _{y_2}q\right ] \end{equation*}

and then

\begin{equation*} |\partial _{y_2}^2q|\leq C\!\left (|\hat {F}|+|\partial _{y_1}(\nabla q)|+|\nabla q|\right ) \end{equation*}

and thus (210) is proved.

Next from (210) and interior regularity [Reference Zhou, Tovkach, Golovaty, Sokolov, Aranson and Lavrentovich64], Theorem 8.8] we have for $\tilde{\Omega } \subset \subset \Omega$ that:

(217) \begin{equation} \|q\|_{H^2(\tilde{\Omega })}\leq C\!\left (\|F\|_{L^2(\Omega )}+\|q\|_{L^2(\Omega )}\right ). \end{equation}

To obtain a bound on $\|q\|_{L^2(\Omega )}$ , we will use that $q$ from (204) minimises the energy functional

(218) \begin{equation} \mathcal{E}_0(q)=\dfrac{1}{2}\int \limits _{\Omega }|\nabla q|^2 \,\textrm{d}x+\dfrac{\beta }{2}\int \limits _{\Omega }|\gamma -q|^2\,\textrm{d}S_x+\int \limits _{\Omega }Fq\,\textrm{d}x. \end{equation}

From $\mathcal{E}_0(q)\leq \mathcal{E}_0(0)$ , we get for all $\delta \gt 0$

\begin{equation*} \int \limits _{\Omega }|\nabla q|^2 \,\text {d}x+\beta \int \limits _{\Omega }|q|^2\,\text {d}S_x\leq C\beta \|\gamma \|_C^2+C\delta ^{-1}\|F\|^{2}_{L^2(\Omega )}+\delta \|q\|_{L^2(\Omega )}^2. \end{equation*}

Next, we use Poincaré estimate (199):

\begin{eqnarray*} \|q\|_{L^2(\Omega )}^2\leq C\!\left (1+\dfrac{1}{\beta }\right )\!\left (\beta \|\gamma \|_C^2+C\delta ^{-1}\|F\|^{2}_{L^2(\Omega )}+\delta \|q\|_{L^2(\Omega )}^2\right ). \end{eqnarray*}

Take $\delta \;:\!=\;\dfrac{1}{2C}\!\left (1+\dfrac{1}{\beta }\right )^{-1}$ and we get

\begin{eqnarray*} \|q\|_{L^2(\Omega )}^2\leq C(\beta +1)\|\gamma \|_C^2+C\!\left (1+\dfrac{1}{\beta }\right )^2\|F\|^{2}_{L^2(\Omega )}. \end{eqnarray*}

Finally, we conclude that

\begin{eqnarray*} \|q\|^2_{H^2(\Omega )}\leq C(\beta +1)\|\gamma \|_C^2+C\!\left (1+\dfrac{1}{\beta }\right )^2\|F\|^{2}_{L^2(\Omega )}. \end{eqnarray*}

C. Rescaling

In this Appendix, we present non-dimensionalisation of the steady-state problem, showing how the scalings in (44)–(47) arise. We will assume that all quantities are in their physical dimensions. Representative values of physical parameters can be found in Table B.1.

Introduce characteristic length $\delta _{L} = 10^{-2}\,\textrm{m}$ , time $\delta _T = 1\,\textrm{s}$ and force $\delta _f = 1\,\textrm{N}$ . Non-dimensional flow velocity and pressure are

(219) \begin{eqnarray}{\boldsymbol{u}} = \tilde{{\boldsymbol{u}}}\, \frac{\delta _L}{\delta _T}\textrm{ and } p = \tilde{p}\, \frac{\delta _f}{\delta _L^2}. \end{eqnarray}

Note that tensor order parameter $\boldsymbol{Q}$ is non-dimensional and does not require a non-dimensionalisation. We also represent the external alignment field $F_{\textrm{ext}}$ as

\begin{equation*} F_{\text {ext}}=\zeta \tilde {F}_{\text {ext}}, \end{equation*}

where $\tilde{F}_{\textrm{ext}}$ is non-dimensional.

Using $\nabla _{\boldsymbol{x}}=\delta _L^{-1}\nabla _{\tilde{{\boldsymbol{x}}}}$ , PDEs (23) and (26) reduce to

(220) \begin{align} \dfrac{\delta _T\Gamma K}{\delta _L^2}\Delta{\boldsymbol{Q}}+\delta _T\Gamma a\,{\boldsymbol{Q}}-\delta _T\Gamma c\,{\boldsymbol{Q}}\textrm{Tr}({\boldsymbol{Q}}^2)+S(\nabla \tilde{{\boldsymbol{u}}},{\boldsymbol{Q}}) - \tilde{{\boldsymbol{u}}} \cdot \nabla{\boldsymbol{Q}}+\delta _T\zeta \, \tilde{F}_{\textrm{ext}}=0, \end{align}
(221) \begin{align} \dfrac{\rho \delta _L^4}{\delta _T^2\delta _f}(\tilde{\boldsymbol{u}}\cdot \nabla )\tilde{{\boldsymbol{u}}}-\dfrac{\eta \delta _L^2}{\delta _T\delta _f}\Delta \tilde{{\boldsymbol{u}}} + \nabla \tilde{p} = \dfrac{K}{\delta _f} \nabla \cdot (\nabla{{\boldsymbol{Q}}} \odot \nabla{{\boldsymbol{Q}}} +{{\boldsymbol{Q}}}\Delta{{\boldsymbol{Q}}} - \Delta{{\boldsymbol{Q}}}{{\boldsymbol{Q}}}). \end{align}

Here and below in this section, all spatial derivatives are taken with respect to $\tilde{{\boldsymbol{x}}}$ .

Boundary condition (25) and (28) becomes

(222) \begin{equation} \tilde{{\boldsymbol{u}}}=\dfrac{\delta _T}{\delta _L}v_{\textrm{prop}}{u}^{(p)}_{\textrm{sq}}\boldsymbol{\tau }\quad \textrm{and}\quad \partial _\nu{\boldsymbol{Q}}=\dfrac{W\delta _L}{K}({\boldsymbol{Q}}_{\textrm{pref}}-{\boldsymbol{Q}}). \end{equation}

Here, we represented $u_{\textrm{sq}}=v_{\textrm{prop}}{u}^{(p)}_{\textrm{sq}}$ where ${u}^{(p)}_{\textrm{sq}}$ is the profile of the propulsion such that $\max{{u}^{(p)}_{\textrm{sq}}}=1$ and $v_{\textrm{prop}}$ is the propulsion strength. Introduce rescaled parameters:

(223) \begin{eqnarray} &\varepsilon =\dfrac{L}{\delta _L},\quad \gamma =\dfrac{\delta _T\Gamma K}{L\delta _L}, \quad \tilde{a}=\delta _T\Gamma a,\quad \tilde{c}=\delta _T\Gamma c,\quad \tilde{\zeta }=\delta _T\zeta,&\nonumber \\[5pt] &\tilde{\rho }=\dfrac{\rho \delta _L^5}{L\delta _T^2\delta _f}, \quad \tilde{\eta }=\dfrac{\eta \delta _L^3}{L\delta _T\delta _f},\quad \kappa =\dfrac{K\delta _L^2}{L^2\delta _f}, \quad \tilde{W}=\dfrac{W}{\delta _L}{K},\quad \tilde{v}_{\textrm{prop}}=\dfrac{\delta _T}{L}v_{\textrm{prop}}.& \end{eqnarray}

Specific values of these parameters can be found in Table C.1.

Table C.1. Values of non-dimensional parameters introduced in (223) corresponding to values of physical parameters from Table 1

Then PDEs (220) and (221) become

\begin{eqnarray*} &&\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}}=0\textrm{ in }\Omega _{\varepsilon },\\[5pt] &&\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}}}) \textrm{ in }\Omega _{\varepsilon }, \end{eqnarray*}

with boundary conditions

\begin{eqnarray*} &&\tilde{{\boldsymbol{u}}}=\varepsilon \tilde{u}_{\textrm{sq}}\,\boldsymbol{\tau } \textrm{ and }\partial _\nu{\boldsymbol{Q}}=\tilde{W}({\boldsymbol{Q}}_{\textrm{pref}}-{\boldsymbol{Q}}) \textrm{ in } \partial \mathcal{P}_\varepsilon,\\[5pt] && \tilde{{\boldsymbol{u}}} \textrm{ and }{\boldsymbol{Q}} \textrm{ are }2\varepsilon -\textrm{periodic}. \end{eqnarray*}

Here, $\tilde{u}_{\textrm{sq}}=\tilde{v}_{\textrm{prop}}{u}^{(p)}_{\textrm{sq}}$ , $\Omega _{\varepsilon }=\delta _L^{-1}\Omega$ and $\mathcal{P}_\varepsilon =\delta _L^{-1}\mathcal{P}$ .

References

Abels, H., Dolzmann, G. & Liu, Y. (2014) Well-posedness of a fully-coupled Navier-Stokes/q-tensor system with inhomogeneous boundary data. SIAM J. Math. Anal. 46 (4), 30503077.CrossRefGoogle Scholar
Abels, H., Dolzmann, G. & Liu, Y. (2016) Strong solutions for the Beris-Edwards model for nematic liquid crystals with homogeneous Dirichlet boundary conditions. Adv. Differ. Equations 21 (1/2), 109152.CrossRefGoogle Scholar
Alexander, G., Pooley, C. & Yeomans, J. (2009) Hydrodynamics of linked sphere model swimmers. J. Phys. Condens. Matter 21 (29), 204108.CrossRefGoogle ScholarPubMed
Allaire, G. (1991) Homogenization of the Navier-Stokes equations in open sets perforated with tiny holes I. Abstract framework, a volume distribution of holes. Arch. Ration. Mech. Anal. 113 (3), 209259.CrossRefGoogle Scholar
Allaire, G. (1991) Homogenization of the Navier-Stokes equations in open sets perforated with tiny holes II: Non-critical sizes of the holes for a volume distribution and a surface distribution of holes. Arch. Ration. Mech. Anal. 113 (3), 261298.CrossRefGoogle Scholar
Aranson, I. (2022). Bacterial active matter, Rep. Progr. Phys. (in Press),Google ScholarPubMed
Aranson, I. S. (2018) Harnessing medium anisotropy to control active matter. Acc. Chem. Res. 51 (12), 30233030.Google ScholarPubMed
Beris, A. N. & Edwards, B. J. (1994). Thermodynamics of Flowing Systems: With Internal Microstructure, Oxford Engineering Sciences Series, Vol. 36, Oxford University Press, New York. Google Scholar
Bravo, D., Hoare, A., Soto, C., Valenzuela, M. & Quest, A. (2018) Helicobacter pylori in human health and disease: Mechanisms for local gastric and systemic effects. World J. Gastroenterol. 24 (28), 30713089.CrossRefGoogle ScholarPubMed
Chi, H., Potomkin, M., Zhang, L., Berlyand, L. & Aranson, I. (2020) Surface anchoring controls orientation of a microswimmer in nematic liquid crystal. Commun. Phys. 3 (1), 19.CrossRefGoogle Scholar
Chisholm, N., Legendre, D., Lauga, E. & Khair, A. (2016) A squirmer across reynolds numbers. J. Fluid Mech. 796, 233256.CrossRefGoogle Scholar
Dey, K. & Sen, A. (2017) Chemically propelled molecules and machines. J. Am. Chem. Soc. 139 (23), 76667676.CrossRefGoogle ScholarPubMed
Duerinckx, M. (2022) Effective viscosity of random suspensions without uniform separation. Ann. Inst. Henri Poincaré C.CrossRefGoogle Scholar
Duerinckx, M. & Gloria, A. (2020) On Einstein’s effective viscosity formula. Transplant. Rev. (Orlando, Fla.) 33 (1), 1728, arXiv preprint arXiv: 2008.03837.Google Scholar
Duerinckx, M. & Gloria, A. (2021) Corrector equations in fluid mechanics: Effective viscosity of colloidal suspensions. Arch. Ration. Mech. Anal. 239 (2), 10251060.CrossRefGoogle Scholar
Duerinckx, M. & Gloria, A. (2022) Sedimentation of random suspensions and the effect of hyperuniformity. Ann. PDE 8 (1), 166.CrossRefGoogle Scholar
Duerinckx, M. & Gloria, A. (2023) Continuum percolation in stochastic homogenization and the effective viscosity problem. Arch. Ration. Mech. Anal. 247 (2). arXiv preprint arXiv: 2108.09654, 2021.CrossRefGoogle Scholar
Elgeti, J., Winkler, R. & Gompper, G. (2015) Physics of microswimmers - single particle motion and collective behavior: A review. Rep. Prog. Phys. 78 (5), 056601.CrossRefGoogle ScholarPubMed
Elgeti, J., Winkler, R. & Gompper, G. (2015) Physics of microswimmers – single particle motion and collective behavior: A review. Rep. Prog. Phys. 78 (5), 056601.CrossRefGoogle ScholarPubMed
Evans, L. (1998). Partial Differential Equations, American Mathematical Society, Providence, RI.Google Scholar
Figueroa-Morales, N., Dominguez-Rubio, L., Ott, T. & Aranson, I. (2019) Mechanical shear controls bacterial penetration in mucus. Sci. Rep-UK 9 (1), 110.Google ScholarPubMed
Galdi, G. (1999) On the steady self-propelled motion of a body in a viscous incompressible fluid. Arch. Ration. Mech. Anal. 148 (1), 5388.CrossRefGoogle Scholar
Galdi, G. (2011). An Introduction to the Mathematical Theory of the Navier-Stokes Equations: Steady-state Problems, Springer Science Business Media.CrossRefGoogle Scholar
Galdi, G. & Silvestre, A. (2007) The steady motion of a Navier-Stokes liquid around a rigid body. Arch. Ration. Mech. Anal. 184 (3), 371400.CrossRefGoogle Scholar
Genkin, M., Sokolov, A., Lavrentovich, O. & Aranson, I. (2017) Topological defects in a living nematic ensnare swimming bacteria. Phys Rev X 7 (1), 011029.Google Scholar
Gilbarg, D. & Trudinger, N. S. (1983). Elliptic Partial Differential Equations of Second Order, Springer.Google Scholar
Gompper, G., Winkler, R. G., Speck, T., et al. (2020) The 2020 motile active matter roadmap. J.Phys. Condens. Matter 32 (19), 193001.CrossRefGoogle ScholarPubMed
Guillen-Gonzalez, F. & Rodriguez-Bellido, M. (2015) Weak solutions for an initial-boundary $Q$ -tensor problem related to liquid crystals. Nonlinear Anal. Theory Methods Appl. 112, 84104.CrossRefGoogle Scholar
Guo, H., Nawroth, J., Ding, Y. & Kanso, E. (2014) Cilia beating patterns are not hydrodynamically optimal. Phys. Fluids 26 (9), 091901.CrossRefGoogle Scholar
Guo, H., Zhu, H., Liu, R., Bonnet, M. & Veerapaneni, S. (2021) Optimal slip velocities of micro-swimmers with arbitrary axisymmetric shapes. J. Fluid Mech. 910, A26.CrossRefGoogle Scholar
Haines, B., Aranson, I., Berlyand, L. & Karpeev, D. (2008) Effective viscosity of dilute bacterial suspensions: A two-dimensional model. Phys. Biol. 5 (4), 046003.CrossRefGoogle ScholarPubMed
Haines, B. M., Aranson, I. S., Berlyand, L. & Karpeev, D. A. (2010) Effective viscosity of bacterial suspensions: A three-dimensional pde model with stochastic torque. Comm. Pure Appl. Anal. 11 (1), 1946.CrossRefGoogle Scholar
Haines, B. M., Sokolov, A., Aranson, I. S., Berlyand, L. & Karpeev, D. A. (2009) A three-dimensional model for the effective viscosity of bacterial suspensions. Phys. Rev. E 80 (4), 041922.Google ScholarPubMed
Hernandez-Ortiz, J., Stoltz, C. & Graham, M. (2005) Transport and collective dynamics in suspensions of confined swimming particles. Phys. Rev. Lett. 95 (20), 204501.CrossRefGoogle ScholarPubMed
Hernandez-Ortiz, J., Underhill, P. & Graham, M. (2009) Dynamics of confined suspensions of swimming particles. J. Phys. Condensed Matter 21 (20), 204107.CrossRefGoogle ScholarPubMed
Kanevsky, A., Shelley, M. & Tornberg, A. (2010) Modeling simple locomotors in Stokes flow. Conf. Proc. Lec. N Phys. 229 (4), 958977.Google Scholar
Lauga, E. (2020). The Fluid Dynamics of Cell Motility, Cambridge University Press.CrossRefGoogle Scholar
Lauga, E. & Powers, T. (2009) The hydrodynamics of swimming microorganisms. Rep. Prog. Phys. 72 (9), 096601.CrossRefGoogle Scholar
Li, G., Lauga, E. & Ardekani, A. M. (2021) Microswimming in viscoelastic fluids. J. Non-NEWTON Fluid Mech. 297, 104655.Google Scholar
Lighthill, M. (1952) On the squirming motion of nearly spherical deformable bodies through liquids at very small Reynolds numbers. Commun. Pure Appl. Math. 5 (2), 109118.CrossRefGoogle Scholar
Lintuvuori, J. S., Würger, A. & Stratford, K. (2017) Hydrodynamics defines the stable swimming direction of spherical squirmers in a nematic liquid crystal. Phys. Rev. Lett. 119 (6), 068001.CrossRefGoogle Scholar
Lions, J. (1969). Quelques méthodes de résolution de Problemes Aux Limites Non Linéaires, Dunod, Gauthier Villars, Paris.Google Scholar
López, H., Gachelin, J., Douarche, C., Auradou, H. & Clément, E. (2015) Turning bacteria suspensions into superfluids. Phys. Rev. Lett. 115 (2), 028301.CrossRefGoogle ScholarPubMed
Marchetti, M., Joanny, J., Ramaswamy, S., et al. (2013) Hydrodynamics of soft active matter. Rev. Mod. Phys. 85 (3), 11431189.CrossRefGoogle Scholar
Michelin, S. & Lauga, E. (2010) Efficiency optimization and symmetry-breaking in a model of ciliary locomotion. Phys. Fluids 22 (11), 111901.CrossRefGoogle Scholar
Michelin, S. & Lauga, E. (2013) Unsteady feeding and optimal strokes of model ciliates. J. Fluid Mech. 715, 131.CrossRefGoogle Scholar
Najafi, A. & Golestanian, R. (2004) Simple swimmer at low Reynolds number: Three linked spheres. Phys. Rev. E 69 (6), 062901.CrossRefGoogle ScholarPubMed
Paicu, M. & Zarnescu, A. (2011) Global existence and regularity for the full coupled Naiver-Stokes and Q-tensor system. SIAM J. Math. Anal. 43 (5), 20092049.CrossRefGoogle Scholar
Paicu, M. & Zarnescu, A. (2012) Energy dissipation and regularity for a coupled Navier-Stokes and Q-tensor system. Arch. Ration. Mech. Anal. 203 (1), 4567.CrossRefGoogle Scholar
Parsonnet, J., Friedman, G., Vandersteen, D., et al. (1991) Helicobacter pylori infection and the risk of gastric carcinoma. New Engl. J. Med. 325 (16), 11271131.CrossRefGoogle ScholarPubMed
Paxton, W. F., Kistler, K. C., Olmeda, C. C., et al. (2004) Catalytic nanomotors: Autonomous movement of striped nanorods. J. Am. Chem. Soc. 126 (41), 1342413431.CrossRefGoogle ScholarPubMed
Potomkin, M., Kaiser, A., Berlyand, L. & Aranson, I. (2017) Focusing of active particles in a converging flow. New J. Phys. 19 (11), 115005.CrossRefGoogle Scholar
Rubio, L., Potomkin, M., Baker, R., Sen, A., Berlyand, L. & Aranson, I. (2021) Self-propulsion and shear flow align active particles in nozzles and channels. Adv. Intell. Syst. 3 (2), 2000178.CrossRefGoogle Scholar
Ryan, S., Sokolov, A., Berlyand, L. & Aranson, I. (2013) Correlation properties of collective motion in bacterial suspension. New J. Phys. 15 (10), 105021.CrossRefGoogle Scholar
Saintillan, D. & Shelley, M. (2013) Active suspensions and their nonlinear models. C R Phys. 14 (6), 497517.CrossRefGoogle Scholar
Sohr, H. (2012). The Navier-Stokes Equations: An Elementary Functional Analytic Approach, Springer Science Business Media.Google Scholar
Sokolov, A. & Aranson, I. (2009) Reduction of viscosity in suspension of swimming bacteria. Phys. Rev. Lett. 103 (14), 148101.CrossRefGoogle ScholarPubMed
Turiv, T., Koizumi, R., Thijssen, K., et al. (2020) Polar jets of swimming bacteria condensed by a patterned liquid crystal. Nat. Phys. 16 (4), 481487.CrossRefGoogle Scholar
Viney, C. (1999) Mucus liquid crystallinity: Is function related to microstructural domain size? Biorheology 36 (4), 319323.Google ScholarPubMed
Wang, Q. (2019, 2019) Optimal strokes of low Reynolds number linked-sphere swimmers. Appl. Sci. 9 (19), 4023.CrossRefGoogle Scholar
Wong, F., Dey, K. & Sen, A. (2016) Synthetic micro/nanomotors and pumps: Fabrication and applications. Ann. Rev. Mater. Res. 46 (1), 407432.CrossRefGoogle Scholar
Zhou, S. (2018) Recent progresses in lyotropic chromonic liquid crystal research: Elasticity, viscosity, defect structures, and living liquid crystals. Liquid Cryst. Today 27 (4), 91108.CrossRefGoogle Scholar
Zhou, S., Sokolov, A., Lavrentovich, O. D. & Aranson, I. S. (2014) Living liquid crystals. Proc. Natl Acad. Sci. 111 (4), 12651270.CrossRefGoogle ScholarPubMed
Zhou, S., Tovkach, O., Golovaty, D., Sokolov, A., Aranson, I. S. & Lavrentovich, O. D. (2017) Dynamic states of swimming bacteria in a nematic liquid crystal cell with homeotropic alignment. New J. Phys. 19 (5).CrossRefGoogle Scholar
Zick, A. & Homsy, G. (1982) Stokes flow through periodic arrays of spheres. J. Fluid Mech. 115 (-1), 1326.CrossRefGoogle Scholar
Figure 0

Table B.1. Values of physical parameters, taken from [21, 49]

Figure 1

Table C.1. Values of non-dimensional parameters introduced in (223) corresponding to values of physical parameters from Table 1