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 Allaire4–Reference Chi, Potomkin, Zhang, Berlyand and Aranson10] and reviews [Reference Chisholm, Legendre, Lauga and Khair11–Reference 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 Gompper19–Reference 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 Aranson25–Reference 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 Karpeev31–Reference 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 Tornberg36–Reference 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 Aranson52–Reference 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$ :
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
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
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])
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:
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
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
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
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
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
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:
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
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]
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
Remark 2.1. Note that the term $\ell$ admits a simplified form:
Here, we used boundary conditions (6). Next, for any symmetric matrix $\boldsymbol{B}=(B_{ij})_{i,j=1}^{d}$ we have
Indeed, from properties of the Levi–Civita symbol we have
On the other hand, due to symmetry of $\boldsymbol{B}$ we have
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}$ ):
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:
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:
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:
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:
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
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
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
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
which due to the same arguments as in derivation of (31) is equivalent to
Using (32) and the divergence theorem for the first term in the equation above, we get
Therefore, an external force results in the pressure difference
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:
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
and the term $S$ given by (9) satisfies the following equality:
We also impose
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
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
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
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):
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
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
Proof. Recall the equation for $\boldsymbol{Q}$ :
By multiplying the above by $\boldsymbol{Q}$ , taking the trace of the resulting expression and using (40) and (41), we get
As $|\nabla{\boldsymbol{Q}}|^2$ is non-negative, we obtain the inequality:
Now introduce $\psi ({\boldsymbol{Q}})\;:\!=\; (|{\boldsymbol{Q}}|^2- \alpha ^2)_{+}$ ( $\alpha$ is from (49)). Note that
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
The first term in the left-hand side of the above inequality vanishes due to (25) while the second term is negative:
Hence,
Next, by (48), we have
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
Note that the domain $\Pi$ is a bounded periodic box and both Laplacian and Stokes operators have a discrete spectrum implying existence of bases:
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:
Anticipating that ${\boldsymbol{u}}={\boldsymbol{u}}_{\textrm{os}} + \hat{{\boldsymbol{u}}}$ , from (57) and (25), we have
To continue, for an appropriately large constant $M\gt 0$ , we introduce a truncated potential $\hat{\mathcal{F}}_M\geq 0$ as follows:
or, more explicitly,
The functional derivative of $\hat{\mathcal{F}}_{M}$ is given by
We have the following bound on $\hat{H}_M({\boldsymbol{Q}})$ :
We now define the function ${\boldsymbol{Q}}_m$ , corresponding to the Galerkin approximation ${\bf H}_m$ as the solution to the following system:
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
where
We define Galerkin approximations ( $\boldsymbol{u}_m, \textbf{H}_m$ ) as solutions of the two equations below for each $k=1,...,m$ :
Here, $F_{\textrm{ext},M}$ is defined as a continuous function such that:
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
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:
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
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
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:
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:
2. Here, we bound the second integral in the right-hand side of (70) by the Cauchy–Schwarz inequality:
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:
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
5. We use the Cauchy–Schwarz inequality, (59), and the a priori bound (64) to estimate the 5th line of (70):
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):
7. Finally, the last term in (70) is estimated as follows:
Collect (74)–(80) and substitute them in (70):
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$ :
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
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
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
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
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
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
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
Then $p_m\rightharpoonup p$ in $L^q(\mathcal{O})$ .
From (82) and Lemma 4.5, we get
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):
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}}}$ ):
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
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
where we define linear operator $\mathcal{L}$ and nonlinear operator $\mathcal{N}$ as
and
To handle the nonlinear and inhomogeneous boundary conditions, we represent unknown functions $\boldsymbol{u}$ and $\boldsymbol{Q}$ as
The offset function ${\boldsymbol{u}}_{\textrm{os}}$ is given by
The offset function ${\boldsymbol{Q}}_{\textrm{os}}$ is defined such that
Specifically, we define
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:
To describe the domains of the operator $\mathcal{L}$ , we introduce the following Banach spaces:
with corresponding norms
Introduce also
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
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:
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:
Proof of (i):
Proof of (ii):
Proof of (iii):
Next, estimate each term in the right-hand side of (117):
Combining (117)–(121), we obtain (106).
Proof of (iv):
Next, estimate each term in the right-hand side of (122):
Remark 5.1. It is useful to rewrite (107) with the norm of space $X_{\boldsymbol{Q}}$ :
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
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
Due to the stability of the Stokes operator (similar to [Reference Turiv, Koizumi and Thijssen58, Theorem IV.6.1])
and smooth dependence of $u_{\textrm{sq}}$ in $\alpha (t)$ , we have
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
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
(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:
Indeed,
STEP 2. Here, we establish the following inequality:
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
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
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}})$ .
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
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}})$ .
Applying (107) for the first term in the right-hand side of (143), one can get
Applying same arguments for the second term in the right-hand side of (143), one can obtain
Part 3: $\sigma _{s}^1({\boldsymbol{Q}})$ .
The first two terms in the right-hand side of (145) are bounded as follows:
Next, we bound the third (cubic) term in the right-hand side of (145). Note
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
Thus, we have
which, in view of (145), implies
Part 4: $ \sigma _{s}^2({\boldsymbol{Q}})$ . In this part, we will need the following bound:
which can be obtained by applying same arguments as in (145)–(147) for $i=1,2$
Now, we can estimate,
Next, we use the triangle inequality, (127), (148) and (150):
Therefore,
Part 5: $ \sigma _{s}^3({\boldsymbol{Q}})$ .
Next, applying the same arguments as in (147) with bounds (148) and (150) we get
STEP 3. Here, we establish the following inequality:
To this end, similar to how we treated $\sigma _s$ in STEP 2, we split $S(\nabla{\boldsymbol{u}},{\boldsymbol{Q}})$ into three parts:
where
are correspondingly the linear, bilinear and trilinear part of $S$ . First note that using (134), (135) and (136), one gets
and
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
Then
Part 2: $S^2(\nabla{\boldsymbol{u}},{\boldsymbol{Q}})$ .
Apply (107), (138) and (157) to obtain
Applying the similar arguments to the second term in the right-hand side of (161), we obtain
Part 3: $S^3(\nabla{\boldsymbol{u}},{\boldsymbol{Q}})$ .
Using same arguments as in (147) and taking into account (136) and (138), we obtain
Applying similar arguments for the other two terms in the right-hand side of (163), we obtain
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:
We apply the same arguments in Part 1 of STEP 2, that is, apply (106), along with (157) and (158), to get:
Thus, we obtained
Part 2: ${\boldsymbol{u}} \cdot \nabla{\boldsymbol{Q}}$ .
Apply the same arguments as in Part 1 of STEP 2, that is, apply (106):
Using (157), (158) and (138) one gets
Thus, we obtained
Part 3: $\hat{{\bf H}}({\boldsymbol{Q}})$ . Here, we need to show
This bound follows directly from the proof on Part 3 of STEP 2.
Part 4: $\dfrac{\textrm{d}{\boldsymbol{V}}}{\textrm{d}t}$ .
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}$ :
where coefficients $b_{klmn}$ depend on ${\boldsymbol{Q}}_{\infty }$ . Then use of the triangle inequality and (127):
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)
Using trace theorem, we get
To estimate the term with $\ell$ , recall its simplified form (19):
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
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
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
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
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
We will frequently use following identities for $f({\boldsymbol{x}},\varepsilon )=\bar{f}({\boldsymbol{x}},{\boldsymbol{y}})$ with ${\boldsymbol{y}}=\varepsilon ^{-1}{\boldsymbol{x}}$ :
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}$ :
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})$ :
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$ :
Rewrite the first integral in (179) in domain $\Omega _1$ :
Expanding analogously other terms in (179) and using that ${\boldsymbol{u}}^{(0)}=0$ , we get at level $\varepsilon ^{d-1}$ :
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
Note that the above integral relation is the weak formulation for the following boundary-value problem:
Here,
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:
To evaluate the right-hand side, we use the fact that ${\boldsymbol{Q}}^{(0)}$ is independent of $\boldsymbol{y}$ and
Hence, we have
Substituting (185) and (186) into (183), we get the equation for ${\boldsymbol{Q}}^{(0)}$ :
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
we get
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
Thus, we can rewrite (188) as
The above can be rewritten in component-wise form as
where
Next, rewrite (191) in a vectorial form:
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)}$ :
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:
The term $\overline{{\boldsymbol{u}}}_{\textrm{sq}}$ is defined as the solution of
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)}$ :
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:
(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.
(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].
-
(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].
-
(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
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:
By contradiction, we assume that there exists a sequence $\{{\boldsymbol{Q}}_n\}_{n=1}^{\infty }$ :
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 )$ :
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
Using integration by parts, we get
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:
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:
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
This shows (64). Next, due to the elliptic regularity result (see Appendix B)
we have
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:
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
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
In variable $\boldsymbol{y}$ , the problem (204) has the form:
where
We note that $L$ is a positive definite symmetric matrix with the smallest eigenvalue
Thus, $L$ is uniformly positive definite:
Lemma B.1. Let $q$ be the solution of
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:
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)$ :
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:
Integration by parts for the difference quotient allows us to rewrite the first integral in (211) as follows:
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):
Finally, we estimate the third integral in (211) using [Reference Zick and Homsy65], Theorem 3 from Section 5.8.2] as follows:
Combining (211) with estimates (213), (214) and (215), we get
Using [Reference Zick and Homsy65, Theorem 3(i) from Section 5.8.2] we get
Then [Reference Zick and Homsy65, Theorem 3(ii) from Section 5.8.2] implies
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,
Analogous estimate is valid for $\partial ^2_{x_2}q$ since
and then
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:
To obtain a bound on $\|q\|_{L^2(\Omega )}$ , we will use that $q$ from (204) minimises the energy functional
From $\mathcal{E}_0(q)\leq \mathcal{E}_0(0)$ , we get for all $\delta \gt 0$
Next, we use Poincaré estimate (199):
Take $\delta \;:\!=\;\dfrac{1}{2C}\!\left (1+\dfrac{1}{\beta }\right )^{-1}$ and we get
Finally, we conclude that
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
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
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
Here and below in this section, all spatial derivatives are taken with respect to $\tilde{{\boldsymbol{x}}}$ .
Boundary condition (25) and (28) becomes
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:
Specific values of these parameters can be found in Table C.1.
Then PDEs (220) and (221) become
with boundary conditions
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}$ .