Far Field Asymptotics of Nematic Flows Around a Small Spherical Particle

Dmitry Golovaty Department of Mathematics, The University of Akron, Akron, OH 44124, USA Nung Kwan Yip Department of Mathematics, Purdue University, West Lafayette, IN 47907, USA

Abstract

Given a small spherical particle, we consider flow of a nematic liquid crystal in the corresponding exterior domain. Our focus is on precise far field asymptotic behavior of the flow in a parameter regime when the governing equations can be reduced to a system of linear partial differential equations. We are able to analytically characterize the velocity of the flow and compare it to the classical expression for the Stokes flow. The expression for velocity away from the particle can be computed either numerically or symbolically.

1 Introduction

This paper analyzes the flow pattern around a small spherical particle in a nematic liquid crystalline environment. Although a similar phenomenon in a classical isotropic fluid has been well-studied, the understanding in the case of a complex fluid is far from complete, despite the fact that it is of utmost importance in many physical and biological systems. Here we are particularly interested in flows around solid particles immersed in a nematic liquid crystalline medium. A typical nematic liquid crystal consists of molecules which possess a degree of orientational order but no positional order so that this medium can be thought of as an anisotropic Newtonian fluid with additional elastic properties. Hence, in order to understand the flow field around a particle, one needs to analyze the behavior of solutions to a coupled system of PDEs that describes both flow and orientational elasticity.

In a broader context, we are interested in multi-particle systems such as colloidal suspensions of both passive and active particles in complex media. The overall behavior of these systems certainly depends on both the particle-fluid and particle-particle interactions. As a typical initial step, one would consider a dilute system in which particles are far from each other. Therefore, it is crucial to understand far field interactions between the particles.

As a starting point, one then can investigate a single particle in a host medium and determine the far-field asymptotics of the corresponding flow patterns. In the present case, these asymptotics should also include the far-field orientational information about the liquid crystal. This is the situation considered in our work for the parameter regime in which we are able to obtain explicit asymptotics in the same spirit as those known for the classical Stokes flow. Our work seems to be the first to give a precise characterization of the flow around a particle immersed in a nematic liquid crystal.

First, we briefly discuss prior work on liquid crystalline flows with or without particles. The theory of nematic flows has proceeded along two somewhat different but related directions that rely on distinct continuum descriptions of the nematics. The first and most widely used approach, due to Ericksen and Leslie [1, 2], is based on the $\mathbb{S}^{2}$ -valued nematic director used to model local orientation of nematic molecules, along with the velocity of the flow. The corresponding system of equations is derived by assuming local balances of mass as well as linear and angular momentum and then specifying the energy dissipation density in the form that conforms with the second law of thermodynamics. As the result, the system consists of a coupled Navier-Stokes and Ginzburg-Landau type equations. The Ericksen-Leslie model has been a subject of numerous studies [3, 4, 5, 6, 7].

The alternative approach relies on the Landau-de Gennes description of orientational elasticity that is designed to take into account certain symmetry properties of orientational distribution of nematic molecules. The Landau-de Gennes variational theory replaces the nematic director with a so-called $\mathsf{Q}$ -tensor—a symmetric traceless $3\times 3$ matrix the eigenvectors of which encode orientational properties of the nematic. The energetics of $\mathsf{Q}$ is modeled by the Landau-de Gennes functional $\mathcal{F}_{\mathrm{LdG}}(\mathsf{Q})$ – see (1). When applied to nematic flow, the evolution of $\mathsf{Q}$ and the velocity ${\bf v}$ of the nematic is given by a coupled system of the Landau-de Gennes and Navier-Stokes equations and can be derived following the procedure outlined by Sonnet and Virga (SV) on the basis of the principle of minimum constrained dissipation [8, 9, 10]. Another related but more popular model, due to Beris and Edward (BE) [11], was obtained by making use of the concept of a Poisson bracket with dissipation–in a later work, we will in fact show that the BE model is a specific example of the SV formulation. The BE model has been studied extensively in recent years, see for example, [12, 13, 14] just to name a few.

Particulate flows in a classical fluid is another widely studied subject, see [15, 16] for overview of relevant literature and mathematical theory. Motion of particles in a complex fluid, and in particular in a nematic host, has also attracted attention of many investigators, especially with the expanding interest in motion of immersed active particles [17, 18, 19, 20]. The experimental and modeling work in this field has so far outpaced the analysis effort [21, 22, 23]. The goal of the present paper is to make inroads into rigorous understanding of a nematic flow around a single particle within the framework of a $\mathsf{Q}$ -tensor-based model.

The details of the model are outlined in Section 2. Here we highlight some key physically justifiable assumptions that we make in order to obtain precise asymptotics of the flow pattern. These assumptions allow us to reduce the overall dynamics to a linear system of PDEs describing an anisotropic Stokes flow with elastic contribution. We emphasize that the distinction with the standard isotropic Stokes flow is not only due to orientability of the medium (described by $\mathsf{Q}$ ) that makes the flow anisotropic, but also because the flow is coupled to the elastic properties of the medium (described by $\nabla\mathsf{Q}$ ).

Our principal assumptions are as follows. We set both Ericksen (Er) and Reynolds (Re) numbers to be small. This means that the dynamics is dominated by the elastic effects of the nematics and the flow is highly dissipative in the sense that inertial effects are negligible. As the result, our model becomes a partially decoupled system (32)–(34) in which $\mathsf{Q}$ satisfies an equation that does not involve ${\bf v}$ . The solution of this equation can then be treated as a prescribed function which enters into an inhomogeneous linear Stokes-like equation describing the evolution of ${\bf v}$ .

To obtain an explicit solution for $\mathsf{Q}$ , we further take advantage of the small particle limit considered in [24]. In this regime, the equation for $\mathsf{Q}$ becomes a Laplace equation with an exact solution given by a harmonic function.

Our main result is the precise far-field asymptotics for ${\bf v}$ where we are able to analytically characterize the deviation of ${\bf v}$ from the solution for the classical Stokes flow. In particular, we are able to obtain an analytical expression for ${\bf v}$ that can be computed either numerically or symbolically.

The outline of our paper is as follows. In Section 2, we provide a detailed description of our model and relevant parameter regimes. In Section 3, we recall the explicit solutions for a $\mathsf{Q}$ -tensor and the classical Stokes flow in the exterior of a small particle that we use in subsequent sections. In Section 4, we analyze the structure of our anisotropic Stokes system. Then in Sections 5 and 6, we prove existence of solutions to the governing system of equations and analyze their far-field asymptotic behavior, in particular their deviation from the isotropic Stokes flow. Finally, in Section 7, we present numerical results. Given the generality of the model considered in this work, the analysis relies on some tedious but routine computations that we summarize in the Appendices A-D, along with verification of our numerical experiments.

We conclude this section with the list of notation and conventions that will be used in this paper. We remind the reader that we work with objects defined on $\mathbb{R}^{3}$ .

(a).

For any vector $u\in\mathbb{R}^{3}$ , its components are represented by $u_{i}$ , for $i=1,2,3$ and $\displaystyle\widehat{u}=\frac{u}{|u|}$ is the unit vector with the same direction as $u$ . Hence we have $\displaystyle\widehat{u}_{i}=\frac{u_{i}}{|u|}$ . Whenever there is no ambiguity, we will omit $\,\,\widehat{}\,\,$ . Euclidean inner product between vectors $u,v\in\mathbb{R}^{3}$ is denoted by $u\cdot v$ . Furthermore, the symbols $e_{1},e_{2}$ and $e_{3}$ represent the coordinate vectors of $\mathbb{R}^{3}$ .
(b).

For any $\mathsf{B},\mathsf{C}\in M^{3\times 3}$ , $3\times 3$ matrices, $\mathsf{B}\cdot\mathsf{C}=\mathrm{tr\,}\left(\mathsf{C}^{T}\mathsf{B}\right)$ and $|\mathsf{B}|^{2}=\mathsf{B}\cdot\mathsf{B}$ .
(c).

Einstein’s convention for repeated indices will be used whenever possible.
(d).

For any matrix or second order tensor field $\mathsf{T}=(\mathsf{T}_{ij})_{1\leq i,j\leq 3}$ , $\mathsf{T}_{i}$ denotes its $i$ -th row. The divergence $\mathrm{div}\,\mathsf{T}$ is to be taken row-wise, i.e. for all $i$ , $(\mathrm{div}\,\mathsf{T})_{i}:=\mathrm{div}\,\mathsf{T}_{i}=\partial_{x_{j}}\mathsf{T}_{ij}$ .
(e).

A symbol such as $x=(x_{1},x_{2},x_{3})^{T}$ denotes either a generic point or a dummy integration variable in $\mathbb{R}^{3}$ . For convenience, we will also use $r$ to denote

$r=|x|=\sqrt{x_{1}^{2}+x_{2}^{2}+x_{3}^{2}}.$
(f).

The (3-dimensional) volume integration element is denoted by

$dx=dx_{1}dx_{2}dx_{3}.$

The (2-dimensional) area integration element is denoted by

$d\sigma_{x},\,\,\,\text{or simply}\,\,\,d\sigma\,\,\,\text{(if the variable is clear from context).}$

The unit outward normal to the domain $\Omega$ occupied by the liquid crystal is given by $\nu$ .
(g).

We use the $\partial_{i}=\partial_{x_{i}}$ to denote spatial partial derivatives.
(h).

Sometimes we find it advantageous to use the symbol $\big{\langle}\text{matrix},\text{vector}\big{\rangle}$ and $\big{\langle}\text{vector},\text{vector}\big{\rangle}$ to denote the matrix-vector multiplication and the dot product, respectively.
(i).

The symbol $``:"$ will be used to denote a generic tensor contraction. It will be explicitly defined whenever necessary.
(j).

The symbol $``\lesssim"$ will be used to refer to an inequality that holds up to a multiplicative constant. This constant can change from one line to another but is not important to the ultimate conclusion.

2 Problem formulation

2.1 A model for the dynamics of nematic liquid crystals

Suppose that a nematic liquid crystal occupies a region $\Omega\subset\mathbb{R}^{3}$ and let the $\mathsf{Q}$ -tensor

\mathsf{Q}:\Omega\to\mathcal{M}:=\left\{\mathsf{Q}\in M^{3\times 3}(\mathbb{R}):\mathsf{Q}^{T}=\mathsf{Q},\ \mathrm{tr\,}{\mathsf{Q}}=0\right\}

and the velocity ${\bf v}:\Omega\to\mathbb{R}^{3}$ describe the local state of the nematic. We introduce the Landau-de Gennes energy density and the corresponding energy functional as:

\mathcal{E}_{\mathrm{LdG}}(\nabla\mathsf{Q},\mathsf{Q})=\frac{K}{2}{|\nabla\mathsf{Q}|}^{2}+f(\mathsf{Q}),\quad\mathcal{F}_{\mathrm{LdG}}(\mathsf{Q})=\int\mathcal{E}_{\mathrm{LdG}}(\nabla\mathsf{Q},\mathsf{Q})\,d^{3}x.

(1)

The parameter $K$ is the elastic constant of the liquid crystal. The nonlinear Landau-de Gennes potential $f$ is given by

f(\mathsf{Q})=-\frac{A}{2}\mathrm{tr\,}\mathsf{Q}^{2}+\frac{B}{3}\mathrm{tr\,}\mathsf{Q}^{3}+\frac{C}{4}{\left(\mathrm{tr\,}\mathsf{Q}^{2}\right)}^{2},

(2)

with $A,C>0$ . The minimum of $f$ over $\mathcal{M}$ is attained at the nematic states given by

s_{*}\left({n_{*}}\otimes{n_{*}}-\frac{1}{3}\mathsf{I}\right),\,\,\,\text{where}\,\,\,{n_{*}}\in{\mathbb{S}}^{2}\,\,\,\text{and}\,\,\,s_{*}=\left\{\begin{array}[]{ll}-\frac{B+\sqrt{B^{2}+24AC}}{4C}&\text{for $B>0$},\\ \frac{-B+\sqrt{B^{2}+24AC}}{4C}&\text{for $B<0$}.\end{array}\right.

(3)

Let

\mathring{\mathsf{Q}}=\dot{\mathsf{Q}}+\mathsf{Q}\mathsf{W}-\mathsf{W}\mathsf{Q},\quad\mathsf{A}=\frac{1}{2}\left(\nabla{\bf v}+\nabla{\bf v}^{T}\right),\quad\mathsf{W}=\frac{1}{2}\left(\nabla{\bf v}-\nabla{\bf v}^{T}\right)

(4)

where $\dot{\mathsf{Q}}$ denotes the convective derivative

\dot{\mathsf{Q}}=\partial_{t}\mathsf{Q}+{\bf v}\cdot\nabla\mathsf{Q}.

The expression $\mathring{\mathsf{Q}}$ will be used below instead of $\dot{\mathsf{Q}}$ because it is frame indifferent [10, Section 2.1.3, Eq. (2.87), Section 4.1.3]. The $3\times 3$ -matrices $\mathsf{A}$ and $\mathsf{W}$ are respectively symmetric and skew-symmetric, and the expressions such as $\mathsf{Q}\mathsf{W}$ and so forth refer to matrix multiplications. From the definition of $\mathring{\mathsf{Q}}$ , we infer that it is also a symmetric, traceless matrix.

We consider the incompressible flow of the nematic in $\Omega$ described by the following system of equations

$\displaystyle[left=\empheqlbrace]$	$\displaystyle\frac{\partial\mathcal{E}_{\mathrm{LdG}}}{\partial\mathsf{Q}}-\mathrm{div}\,\left[\frac{\partial\mathcal{E}_{\mathrm{LdG}}}{\partial\nabla\mathsf{Q}}\right]-\Lambda\mathsf{I}+\zeta_{1}\mathring{\mathsf{Q}}+\zeta_{2}\mathsf{A}+\frac{\zeta_{3}}{2}(\mathsf{A}\mathsf{Q}+\mathsf{Q}\mathsf{A})+\zeta_{9}(\mathsf{A}\cdot\mathsf{Q})\mathsf{Q}=0,$	(5)
	$\displaystyle\rho\,\dot{\bf v}+\mathrm{div}\,\left[p\mathsf{I}-\mathsf{T}^{\mathrm{v}}_{\text{SV}}-\mathsf{T}^{\mathrm{el}}\right]=0,$	(6)
	$\displaystyle\mathrm{div}\,{\bf v}=0,$	(7)

where we have

	$\displaystyle\frac{\partial\mathcal{E}_{\mathrm{LdG}}}{\partial\mathsf{Q}}$	$\displaystyle=$	$\displaystyle-A\mathsf{Q}+B\left(\mathsf{Q}^{2}-\frac{1}{3}\|\mathsf{Q}\|^{2}\mathsf{I}\right)+C\|\mathsf{Q}\|^{2}\mathsf{Q},$		(8)
	$\displaystyle\mathrm{div}\,\left[\frac{\partial\mathcal{E}_{\mathrm{LdG}}}{\partial\nabla\mathsf{Q}}\right]$	$\displaystyle=$	$\displaystyle K\triangle\mathsf{Q}.$		(9)

The parameter $\rho$ denotes the density of the nematic and $\zeta_{1}$ through $\zeta_{11}$ are the viscosity coefficients of the nematic with $\zeta_{8}$ being the isotropic viscosity. The pressure $p$ and the function $\Lambda$ are the Lagrange multipliers corresponding, respectively, to the incompressibility and tracelessness constraints for ${\bf v}$ and $\mathsf{Q}$ . The viscous and elastic stress tensors are

$\displaystyle\mathsf{T}^{\mathrm{v}}_{\text{SV}}$	$\displaystyle=$	$\displaystyle\zeta_{1}\left(\mathsf{Q}\mathring{\mathsf{Q}}-\mathring{\mathsf{Q}}\mathsf{Q}\right)+\zeta_{2}\left(\mathring{\mathsf{Q}}+\mathsf{Q}\mathsf{A}-\mathsf{A}\mathsf{Q}\right)+\frac{\zeta_{3}}{2}\left(\mathsf{Q}\mathring{\mathsf{Q}}+\mathring{\mathsf{Q}}\mathsf{Q}+\mathsf{Q}^{2}\mathsf{A}-\mathsf{A}\mathsf{Q}^{2}\right)$	(10)
		$\displaystyle+\zeta_{4}\left(\mathsf{Q}\mathsf{A}+\mathsf{A}\mathsf{Q}\right)+\zeta_{5}\left(\mathsf{Q}^{2}\mathsf{A}+\mathsf{A}\mathsf{Q}^{2}\right)+\zeta_{6}(\mathsf{A}\cdot\mathsf{Q})\mathsf{Q}+\zeta_{7}{\|\mathsf{Q}\|}^{2}\mathsf{A}+\zeta_{8}\mathsf{A}$
		$\displaystyle+\zeta_{9}(\mathring{{\mathsf{Q}}}\cdot\mathsf{Q})\mathsf{Q}+\zeta_{10}\left(\left(\mathsf{Q}^{2}\cdot\mathsf{A}\right)\mathsf{Q}+\left(\mathsf{Q}\cdot\mathsf{A}\right)\mathsf{Q}^{2}\right)+\zeta_{11}\|\mathsf{Q}\|^{2}(\mathsf{A}\cdot\mathsf{Q})\mathsf{Q}$

and

$\displaystyle\mathsf{T}^{\mathrm{el}}$	$\displaystyle=$	$\displaystyle-\frac{\partial\mathcal{E}_{\mathrm{LdG}}}{\partial\left(\partial_{j}\mathsf{Q}_{mn}\right)}\left(\partial_{i}\mathsf{Q}_{mn}\right){\mathbf{e}}_{i}\otimes{\mathbf{e}}_{j}=-K(\partial_{j}\mathsf{Q}_{mn})(\partial_{i}\mathsf{Q}_{mn}){\mathbf{e}}_{i}\otimes{\mathbf{e}}_{j}$	(11)
	$\displaystyle=$	$\displaystyle-K\left(\partial_{i}\mathsf{Q}\cdot\partial_{j}\mathsf{Q}\right){\mathbf{e}}_{i}\otimes{\mathbf{e}}_{j}$
	$\displaystyle=:$	$\displaystyle-K\left(\nabla\mathsf{Q}\odot\nabla\mathsf{Q}\right),$

respectively, so that $\left(\nabla\mathsf{Q}\odot\nabla\mathsf{Q}\right)_{ij}=\partial_{i}\mathsf{Q}\cdot\partial_{j}\mathsf{Q}$ component-wise.

The model (5)–(7) is established in [10, Section 4.1.3], where a dissipation function $R=R(\mathsf{Q};\mathsf{A},\mathring{\mathsf{Q}})$ is introduced from which the viscous stress tensor (10) can be derived. In particular, the following general form of a dissipation function is proposed in [10, Eq. (4.23)]:

R(\mathsf{Q};\mathsf{A},\mathring{\mathsf{Q}})=\frac{\zeta_{1}}{2}\mathring{\mathsf{Q}}\cdot\mathring{\mathsf{Q}}+\zeta_{2}\mathsf{A}\cdot\mathring{\mathsf{Q}}+\zeta_{3}(\mathring{\mathsf{Q}}\mathsf{Q})\cdot\mathsf{A}+\zeta_{4}\mathsf{Q}\cdot\mathsf{A}^{2}+\zeta_{5}\mathsf{Q}^{2}\cdot\mathsf{A}^{2}+\frac{\zeta_{6}}{2}(\mathsf{Q}\cdot\mathsf{A})^{2}\\ +\frac{\zeta_{7}}{2}{|\mathsf{A}|}^{2}{|\mathsf{Q}|}^{2}+\frac{\zeta_{8}}{2}{|\mathsf{A}|}^{2}+\zeta_{9}(\mathring{\mathsf{Q}}\cdot\mathsf{Q})(\mathsf{A}\cdot\mathsf{Q})+\zeta_{10}\left(\mathsf{Q}^{2}\cdot\mathsf{A}\right)\left(\mathsf{Q}\cdot\mathsf{A}\right)+\frac{\zeta_{11}}{2}|\mathsf{Q}|^{2}(\mathsf{A}\cdot\mathsf{Q})^{2}.

(12)

We remark that [10] considers only the terms involving $\zeta_{1}$ to $\zeta_{8}$ so that $R$ is exactly quadratic in the rates $\mathring{\mathsf{Q}}$ and $\mathsf{A}$ and at most quadratic in $\mathsf{Q}$ but the idea can certainly be generalized. In particular, $R$ can involve a linear combination of nineteen invariants of the tensor triple $(\mathsf{Q},\mathring{\mathsf{Q}},\mathsf{A})$ [10, p. 223]. A simpler version of (12), introduced in [25], given by $\zeta_{3}=\zeta_{5}=\zeta_{7}=0$ , also appears in [10, Eq. (4.25)]. Here, we include the terms with coefficients $\zeta_{9}$ to $\zeta_{11}$ because with these terms the model (5)-(7) subsumes the model in [11], derived by Beris and Edward. In a forthcoming work, we will show that the Beris-Edward model is, in fact, a particular version of (5)-(7), corresponding to a specific choice of dissipation constants $\zeta_{i}$ , $i=1,\ldots,11$ .

The above Sonnet-Virga model is derived using a variational framework together with the Principle of Minimum Constrained Dissipation [10, Section 2.2, Eqn. (2.172), (2.178), (2.179)]. Using $R$ , we can rewrite (5)-(7) as – see also [10, Eq. (4.21) and (4.22)]:

$\displaystyle[left=\empheqlbrace]$	$\displaystyle\frac{\partial\mathcal{E}_{\mathrm{LdG}}}{\partial\mathsf{Q}}-\mathrm{div}\,\left[\frac{\partial\mathcal{E}_{\mathrm{LdG}}}{\partial\nabla\mathsf{Q}}\right]+\frac{\partial R}{\partial\mathring{\mathsf{Q}}}=0,$	(13)
	$\displaystyle\rho\,\dot{\bf v}-\mathrm{div}\,\left[\mathsf{T}\right]=0,$	(14)
	$\displaystyle\mathrm{div}\,{\bf v}=0,$	(15)

where

\mathsf{T}=-p\mathsf{I}+\mathsf{T}_{\text{SV}}^{\text{v}}+\mathsf{T}^{\text{el}}

(16)

and

\mathsf{T}^{\mathrm{v}}_{\text{SV}}=\frac{\partial R}{\partial\mathsf{A}}+\mathsf{Q}\frac{\partial R}{\partial\mathring{\mathsf{Q}}}-\frac{\partial R}{\partial\mathring{\mathsf{Q}}}\mathsf{Q}\,\,\,\,\,\,\text{and}\,\,\,\,\,\,\mathsf{T}^{\text{el}}=-\nabla\mathsf{Q}\odot\frac{\partial\mathcal{E}_{\mathrm{LdG}}}{\partial\left(\nabla\mathsf{Q}\right)}.

(17)

To complete the description of the model, we point out that the dissipation function $R$ has to be positive semidefinite. The set of conditions satisfied by $\mathsf{Q}$ and $\zeta_{i},\ i=1,\ldots,11$ to ensure that this property holds is too complicated to list here but some of these conditions can be found in [10, pp. 221-222]. However, if $R$ only has the terms associated with $\zeta_{1},\zeta_{2}$ and $\zeta_{8}$ , i.e.

R(\mathsf{Q};\mathsf{A},\mathring{\mathsf{Q}})=\frac{\zeta_{1}}{2}\mathring{\mathsf{Q}}\cdot\mathring{\mathsf{Q}}+\zeta_{2}\mathsf{A}\cdot\mathring{\mathsf{Q}}+\frac{\zeta_{8}}{2}{|\mathsf{A}|}^{2},

then $R$ is positive semidefinite if and only if

\zeta_{1}>0\,\,\,\text{and}\,\,\,\zeta_{2}^{2}\leq\zeta_{1}\zeta_{8}.

(18)

2.2 Non-dimensionalization

Suppose $L$ is the characteristic length of the system and $v$ is the characteristic velocity. With these, we introduce

\tilde{x}=\frac{x}{L},\quad\tilde{{\bf v}}=\frac{{\bf v}}{v},\quad\tilde{t}=\frac{vt}{L},\quad\tilde{f}(\mathsf{Q})=\frac{f(\mathsf{Q})}{C},\,\,\,\text{and}\,\,\,\tilde{\mathcal{E}}_{\mathrm{LdG}}={\frac{L^{2}}{K}\mathcal{E}_{\mathrm{LdG}}},

so that

\partial_{t}=\frac{v}{L}\partial_{\tilde{t}},\quad\partial_{x}=\frac{1}{L}\partial_{\tilde{x}},\quad\mathsf{Q}=\tilde{\mathsf{Q}},\quad\mathsf{A}=\frac{v}{L}\tilde{\mathsf{A}},\quad\tilde{\mathsf{W}}=\frac{v}{L}\tilde{\mathsf{W}},\quad\dot{\mathsf{Q}}=\frac{v}{L}\dot{\widetilde{\mathsf{Q}}},\quad\mathring{\mathsf{Q}}=\frac{v}{L}\mathring{\widetilde{\mathsf{Q}}},

\tilde{f}(\mathsf{Q})=-\frac{A}{2C}\mathrm{tr\,}\mathsf{Q}^{2}+\frac{B}{3C}\mathrm{tr\,}\mathsf{Q}^{3}+\frac{1}{4}{\left(\mathrm{tr\,}\mathsf{Q}^{2}\right)}^{2},\quad\text{and}\quad\tilde{\mathcal{E}}_{\mathrm{LdG}}(\tilde{\mathsf{Q}})=\frac{1}{2}{|\tilde{\nabla}\tilde{Q}|}^{2}+\frac{CL^{2}}{K}\tilde{f}(\tilde{\mathsf{Q}}).

The forms of $\tilde{f}$ and $\tilde{\mathcal{E}}_{\mathrm{LdG}}$ are consistent with the following change of variable relation for the Ginzburg-Landau functional,

\mathcal{F}_{\mathrm{LdG}}(\mathsf{Q})=\int\frac{K}{2}|\nabla\mathsf{Q}|^{2}+f(\mathsf{Q})\,dx=KL\int\frac{1}{2}|\tilde{\nabla}\tilde{\mathsf{Q}}|^{2}+\frac{CL^{2}}{K}\tilde{f}(\tilde{\mathsf{Q}})\,d\tilde{x}=KL\int\tilde{\mathcal{E}}_{\mathrm{LdG}}(\tilde{Q})\,d\tilde{x}.

If we further introduce $\displaystyle{\varepsilon}=\frac{1}{L}\sqrt{\frac{K}{C}}$ as the ratio between the nematic correlation length $\displaystyle\sqrt{\frac{K}{C}}$ and $L$ , we can then write the energy density in the following non-dimensional form,

\tilde{\mathcal{E}}_{\mathrm{LdG}}(\tilde{\nabla}\tilde{\mathsf{Q}},\tilde{\mathsf{Q}})=\frac{1}{2}{|\tilde{\nabla}\tilde{Q}|}^{2}+\frac{1}{{\varepsilon}^{2}}\tilde{f}(\tilde{\mathsf{Q}}).

If ${\varepsilon}\ll 1$ , then the potential function imposes a heavy penalty on deviations of $\mathsf{Q}$ from the nematic states (3).

The nondimensional system of the governing equations (5)-(7) is then

	$\displaystyle\frac{\partial\tilde{\mathcal{E}}_{\mathrm{LdG}}}{\partial\tilde{\mathsf{Q}}}-\widetilde{\mathrm{div}\,}\left[\frac{\partial\tilde{\mathcal{E}}_{\mathrm{LdG}}}{\partial\tilde{\nabla}\tilde{\mathsf{Q}}}\right]-\tilde{\Lambda}\mathsf{I}+\mathrm{Er}\left(\gamma_{1}\mathring{\tilde{\mathsf{Q}}}+\gamma_{2}\tilde{\mathsf{A}}+\frac{\gamma_{3}}{2}(\tilde{\mathsf{A}}\tilde{\mathsf{Q}}+\tilde{\mathsf{Q}}\tilde{\mathsf{A}})+\gamma_{9}(\tilde{\mathsf{A}}\cdot\tilde{\mathsf{Q}})\tilde{\mathsf{Q}}\right)=0,$		(19)
	$\displaystyle\mathrm{Re}\,\dot{\tilde{{\bf v}}}+\widetilde{\mathrm{div}\,}\left[\tilde{p}\mathsf{I}-\tilde{\mathsf{T}}^{\mathrm{v}}_{\text{SV}}-\frac{1}{\mathrm{Er}}\tilde{\mathsf{T}}^{\mathrm{el}}\right]=0,$		(20)
	$\displaystyle\widetilde{\mathrm{div}\,}\tilde{\bf v}=0,$		(21)

where the nondimensional viscous and elastic stress tensors are

$\displaystyle\tilde{\mathsf{T}}^{\mathrm{v}}_{\text{SV}}$	$\displaystyle=$	$\displaystyle\gamma_{1}\left(\tilde{\mathsf{Q}}\mathring{\tilde{\mathsf{Q}}}-\mathring{\tilde{\mathsf{Q}}}\tilde{\mathsf{Q}}\right)+\gamma_{2}\left(\mathring{\tilde{\mathsf{Q}}}+\tilde{\mathsf{Q}}\tilde{\mathsf{A}}-\tilde{\mathsf{A}}\tilde{\mathsf{Q}}\right)+\frac{\gamma_{3}}{2}\left(\tilde{\mathsf{Q}}\mathring{\tilde{\mathsf{Q}}}+\mathring{\tilde{\mathsf{Q}}}\tilde{\mathsf{Q}}+\tilde{\mathsf{Q}}^{2}\tilde{\mathsf{A}}-\tilde{\mathsf{A}}\tilde{\mathsf{Q}}^{2}\right)$	(22)
		$\displaystyle+\gamma_{4}\left(\tilde{\mathsf{Q}}\tilde{\mathsf{A}}+\tilde{\mathsf{A}}\tilde{\mathsf{Q}}\right)+\gamma_{5}\left(\tilde{\mathsf{Q}}^{2}\tilde{\mathsf{A}}+\tilde{\mathsf{A}}\tilde{\mathsf{Q}}^{2}\right)+\gamma_{6}(\tilde{\mathsf{A}}\cdot\tilde{\mathsf{Q}})\mathsf{Q}+\gamma_{7}{\|\tilde{\mathsf{Q}}\|}^{2}\tilde{\mathsf{A}}+\tilde{\mathsf{A}}$
		$\displaystyle+\gamma_{9}(\mathring{\tilde{\mathsf{Q}}}\cdot\tilde{\mathsf{Q}})\tilde{\mathsf{Q}}+\gamma_{10}\left(\left(\tilde{\mathsf{Q}}^{2}\cdot\tilde{\mathsf{A}}\right)\tilde{\mathsf{Q}}+\left(\tilde{\mathsf{Q}}\cdot\tilde{\mathsf{A}}\right)\tilde{\mathsf{Q}}^{2}\right)+\gamma_{11}\|\tilde{\mathsf{Q}}\|^{2}(\tilde{\mathsf{A}}\cdot\tilde{\mathsf{Q}})\tilde{\mathsf{Q}}$

and

\tilde{\mathsf{T}}^{\mathrm{el}}=-\tilde{\nabla}\tilde{\mathsf{Q}}\odot\tilde{\nabla}\tilde{\mathsf{Q}},

(23)

respectively. In the above, the nondimensional groups

\mathrm{Er}=\frac{\zeta_{8}vL}{K}\,\,\,\mbox{ and }\,\,\,\mathrm{Re}=\frac{\rho vL}{\zeta_{8}}

(24)

are respectively, the Ericksen and Reynolds number with

\gamma_{i}=\frac{\zeta_{i}}{\zeta_{8}},\ i=1,\ldots,7,9,\ldots,11.

(25)

We note that the constant $\mathrm{Er}$ gives the ratio between the viscous and elastic forces while $\mathrm{Re}$ gives the ratio between the inertial and viscous forces.

For the rest of the paper, we will refer to our nondimensionalized system (19)–(21) but for simplicity omit the tildes in all variables.

2.3 Decoupling of the governing equations

In this section, we introduce physical regimes in which the far-field behavior of $\mathsf{Q}$ and ${\bf v}$ can be explicitly characterized. In particular, we will assume that the Ericksen number is small so that – as we will see shortly – the equation for the $\mathsf{Q}-$ tensor decouples from the equation for the velocity ${\bf v}$ . To this end, observe that

	$\displaystyle\mathrm{div}\,\mathsf{T}^{\mathrm{el}}$	$\displaystyle=$	$\displaystyle-\mathrm{div}\,\left(\left(\partial_{i}\mathsf{Q}\cdot\partial_{j}\mathsf{Q}\right){\mathbf{e}}_{i}\otimes{\mathbf{e}}_{j}\right)=-\partial_{j}(\partial_{i}\mathsf{Q}\cdot\partial_{j}\mathsf{Q}){\mathbf{e}}_{i}$		(26)
		$\displaystyle=$	$\displaystyle-\nabla\left(\frac{1}{2}\|\nabla\mathsf{Q}\|^{2}\right)-\left(\partial_{i}\mathsf{Q}\cdot\Delta\mathsf{Q}\right){\mathbf{e}}_{i},$		(26)

from which $\mathrm{div}\,\mathsf{T}^{\mathrm{el}}$ can be eliminated. More precisely, using (19), we have

\partial_{i}\mathsf{Q}\cdot\Delta\mathsf{Q}=\frac{1}{{\varepsilon}^{2}}\frac{\partial f}{\partial\mathsf{Q}}\cdot\partial_{i}\mathsf{Q}+\mathrm{Er}\left(\gamma_{1}\mathring{\mathsf{Q}}+\gamma_{2}\mathsf{A}+\frac{\gamma_{3}}{2}(\mathsf{A}\mathsf{Q}+\mathsf{Q}\mathsf{A})+\gamma_{9}(\mathsf{A}\cdot\mathsf{Q})\mathsf{Q}\right)\cdot\partial_{i}\mathsf{Q},

(27)

where the $\Lambda\mathsf{I}$ disappears as $\mathrm{tr\,}{\mathsf{Q}}=0$ : $\Lambda\mathsf{I}\cdot\partial_{i}Q=\Lambda\partial_{i}\text{tr}{\mathsf{Q}}=0$ . Combining (26) and (27), we have

			$\displaystyle\mathrm{div}\,\mathsf{T}^{\mathrm{el}}$
		$\displaystyle=$	$\displaystyle-\nabla\left(\frac{1}{2}\|\nabla\mathsf{Q}\|^{2}+\frac{1}{{\varepsilon}^{2}}f(\mathsf{Q})\right)-\mathrm{Er}\left(\gamma_{1}\mathring{\mathsf{Q}}+\gamma_{2}\mathsf{A}+\frac{\gamma_{3}}{2}(\mathsf{A}\mathsf{Q}+\mathsf{Q}\mathsf{A})+\gamma_{9}(\mathsf{A}\cdot\mathsf{Q})\mathsf{Q}\right)\cdot\partial_{i}\mathsf{Q}{\mathbf{e}}_{i}$
		$\displaystyle=$	$\displaystyle-\nabla\mathcal{E}_{\mathrm{LdG}}(\mathsf{Q})-\mathrm{Er}\left(\gamma_{1}\mathring{\mathsf{Q}}+\gamma_{2}\mathsf{A}+\frac{\gamma_{3}}{2}(\mathsf{A}\mathsf{Q}+\mathsf{Q}\mathsf{A})+\gamma_{9}(\mathsf{A}\cdot\mathsf{Q})\mathsf{Q}\right)\cdot\partial_{i}\mathsf{Q}{\mathbf{e}}_{i}.$

Hence

	$\displaystyle\mathrm{div}\,\left(-p\mathsf{I}+\mathsf{T}^{\mathrm{v}}_{\text{SV}}+\frac{1}{\mathrm{Er}}\mathsf{T}^{\mathrm{el}}\right)\,\,\,=\,\,\,$	$\displaystyle\mathrm{div}\,\left(-p\mathsf{I}-\frac{1}{\mathrm{Er}}\mathcal{E}_{\mathrm{LdG}}(\mathsf{Q})\mathsf{I}+\mathsf{T}^{\mathrm{v}}_{\text{SV}}\right)$
		$\displaystyle-\left(\gamma_{1}\mathring{\mathsf{Q}}+\gamma_{2}\mathsf{A}+\frac{\gamma_{3}}{2}(\mathsf{A}\mathsf{Q}+\mathsf{Q}\mathsf{A})+\gamma_{9}(\mathsf{A}\cdot\mathsf{Q})\mathsf{Q}\right)\cdot\partial_{i}\mathsf{Q}{\mathbf{e}}_{i}.$		(28)

Absorbing the gradient of the Landau-de Gennes energy density into the pressure field leads to the following form of our system (19)–(21):

	$\displaystyle[left=\empheqlbrace]$	$\displaystyle\frac{\partial\mathcal{E}_{\mathrm{LdG}}}{\partial\mathsf{Q}}-\mathrm{div}\,\left[\frac{\partial\mathcal{E}_{\mathrm{LdG}}}{\partial\nabla\mathsf{Q}}\right]-\Lambda\mathsf{I}+\mathrm{Er}\left(\gamma_{1}\mathring{\mathsf{Q}}+\gamma_{2}\mathsf{A}+\frac{\gamma_{3}}{2}(\mathsf{A}\mathsf{Q}+\mathsf{Q}\mathsf{A})+\gamma_{9}(\mathsf{A}\cdot\mathsf{Q})\mathsf{Q}\right)=0,$
		$\displaystyle\mathrm{Re}\,\dot{\bf v}+\mathrm{div}\,\left[p\mathsf{I}-\mathsf{T}^{\mathrm{v}}_{\text{SV}}\right]+\left(\gamma_{1}\mathring{\mathsf{Q}}+\gamma_{2}\mathsf{A}+\frac{\gamma_{3}}{2}(\mathsf{A}\mathsf{Q}+\mathsf{Q}\mathsf{A})+\gamma_{9}(\mathsf{A}\cdot\mathsf{Q})\mathsf{Q}\right)\cdot\partial_{i}\mathsf{Q}{\mathbf{e}}_{i}=0,$
		$\displaystyle\mathrm{div}\,{\bf v}=0.$

In this paper, we consider a specific regime of the above system so that the far-field spatial behavior can be revealed explicitly. This is described as follows.

(a).

The Ericksen number is small, i.e. $\mathrm{Er}\to 0$ , so that the elastic stress in the liquid crystal dominates the viscous stress. Formally, this leads to

	$\displaystyle[left=\empheqlbrace]$	$\displaystyle\frac{\partial\mathcal{E}_{\mathrm{LdG}}}{\partial\mathsf{Q}}-\mathrm{div}\,\left[\frac{\partial\mathcal{E}_{\mathrm{LdG}}}{\partial\nabla\mathsf{Q}}\right]=0,$
		$\displaystyle\mathrm{Re}\,\dot{\bf v}+\mathrm{div}\,\left[p\mathsf{I}-\mathsf{T}^{\mathrm{v}}_{\text{SV}}\right]+\left(\gamma_{1}\mathring{\mathsf{Q}}+\gamma_{2}\mathsf{A}+\frac{\gamma_{3}}{2}(\mathsf{A}\mathsf{Q}+\mathsf{Q}\mathsf{A})+\gamma_{9}(\mathsf{A}\cdot\mathsf{Q})\mathsf{Q}\right)\cdot\partial_{i}\mathsf{Q}{\mathbf{e}}_{i}=0,$
		$\displaystyle\mathrm{div}\,{\bf v}=0$

where there is no need for $\Lambda\mathsf{I}$ as the tracelessness condition for $\mathsf{Q}$ is already incorporated in (8). Note that the first equation is the Euler-Lagrange equation for the Landau-de Gennes energy and it is decoupled from the equation for the velocity. In other words, the tensor field $\mathsf{Q}$ serves as the inhomogeneous source term for the velocity field ${\bf v}$ .

(b).

The characteristic length of the problem is much smaller than the nematic correlation length, i.e., $\varepsilon\to\infty$ so that $\displaystyle\frac{\partial\mathcal{E}_{\mathrm{LdG}}}{\partial\mathsf{Q}}=\frac{1}{\varepsilon^{2}}\nabla f(\mathsf{Q})\approx 0$ . Hence we are led to the following system:

$\displaystyle[left=\empheqlbrace]$	$\displaystyle-\Delta\mathsf{Q}=0,$	(29)
	$\displaystyle\mathrm{Re}\,\dot{\bf v}+\mathrm{div}\,\left[p\mathsf{I}-\mathsf{T}^{\mathrm{v}}_{\text{SV}}\right]=-\left(\gamma_{1}\mathring{\mathsf{Q}}+\gamma_{2}\mathsf{A}+\frac{\gamma_{3}}{2}(\mathsf{A}\mathsf{Q}+\mathsf{Q}\mathsf{A})+\gamma_{9}(\mathsf{A}\cdot\mathsf{Q})\mathsf{Q}\right)\cdot\partial_{i}\mathsf{Q}{\mathbf{e}}_{i},$	(30)
	$\displaystyle\mathrm{div}\,{\bf v}=0,$	(31)

Note that now $\mathsf{Q}$ is harmonic. We remark that the positive semidefiniteness of the dissipation function is ensured if among others, the inequality (18) is satisfied. However, as $\mathsf{Q}$ is fixed or actually “prescribed” by (29) in our asymptotic regime, these inequalities can simply be replaced by the condition that all the coefficients, $\gamma_{1}$ through $\gamma_{11}$ are sufficiently small.

(c).

The Reynolds number is small and the system has reached stationarity, i.e. $\mathrm{Re}\rightarrow 0$ and $\partial_{t}\mathsf{Q}=\partial_{t}{\bf v}=0$ . Hence the system (29)–(31) becomes

$\displaystyle[left=\empheqlbrace]$	$\displaystyle-\Delta\mathsf{Q}=0,$	(32)
	$\displaystyle\mathrm{div}\,\left[p\mathsf{I}-\mathsf{T}^{\mathrm{v}}_{\text{SV}}\right]=-\left(\gamma_{1}\mathring{\mathsf{Q}}+\gamma_{2}\mathsf{A}+\frac{\gamma_{3}}{2}(\mathsf{A}\mathsf{Q}+\mathsf{Q}\mathsf{A})+\gamma_{9}(\mathsf{A}\cdot\mathsf{Q})\mathsf{Q}\right)\cdot\partial_{i}\mathsf{Q}{\mathbf{e}}_{i},$	(33)
	$\displaystyle\mathrm{div}\,{\bf v}=0,$	(34)

Note that in this case, we have $\mathring{\mathsf{Q}}={\bf v}\cdot\nabla\mathsf{Q}+\mathsf{Q}\mathsf{W}-\mathsf{W}\mathsf{Q}$ . We record ${\mathsf{T}}^{\mathrm{v}}_{\text{SV}}$ (22) here again for convenience,

$\displaystyle{\mathsf{T}}^{\mathrm{v}}_{\text{SV}}$	$\displaystyle=$	$\displaystyle\gamma_{1}\left(\mathsf{Q}\mathring{{\mathsf{Q}}}-\mathring{{\mathsf{Q}}}\mathsf{Q}\right)+\gamma_{2}\left(\mathring{{\mathsf{Q}}}+\mathsf{Q}\mathsf{A}-\mathsf{A}\mathsf{Q}\right)+\frac{\gamma_{3}}{2}\left(\mathsf{Q}\mathring{{\mathsf{Q}}}+\mathring{{\mathsf{Q}}}\mathsf{Q}+\mathsf{Q}^{2}\mathsf{A}-\mathsf{A}\mathsf{Q}^{2}\right)$	(35)
		$\displaystyle+\gamma_{4}\left(\mathsf{Q}\mathsf{A}+\mathsf{A}\mathsf{Q}\right)+\gamma_{5}\left(\mathsf{Q}^{2}\mathsf{A}+\mathsf{A}\mathsf{Q}^{2}\right)+\gamma_{6}(\mathsf{A}\cdot\mathsf{Q})\mathsf{Q}+\gamma_{7}{\|\mathsf{Q}\|}^{2}\mathsf{A}+\mathsf{A}$
		$\displaystyle+\gamma_{9}(\mathring{{\mathsf{Q}}}\cdot\mathsf{Q})\mathsf{Q}+\gamma_{10}\left(\left(\mathsf{Q}^{2}\cdot\mathsf{A}\right)\mathsf{Q}+\left(\mathsf{Q}\cdot\mathsf{A}\right)\mathsf{Q}^{2}\right)+\gamma_{11}\|\mathsf{Q}\|^{2}(\mathsf{A}\cdot\mathsf{Q})\mathsf{Q}.$

To conclude, the current paper analyzes the stationary system (32)-(34) in the domain exterior to a sphere of radius $a>0$ . We are particularly interested in the far-field spatial behavior of the flow.

To complete the description of the above system, we need to incorporate boundary conditions for ${\bf v}$ and $\mathsf{Q}$ which are discussed next.

2.4 Boundary conditions

We will solve the above system in the exterior domain $\Omega=\mathbb{R}^{3}\backslash{\bf B}_{a}(0)$ in a moving frame. The following boundary conditions will be imposed for $\mathsf{Q}$ and ${\bf v}$ :

\mathsf{Q}\longrightarrow\mathsf{Q}_{*},\quad{\bf v}\longrightarrow{\bf v}_{*}\quad\text{as $|x|\longrightarrow\infty$.}

(36)

\text{on $\partial\Omega$}:\quad\left\{\begin{array}[]{rl}\frac{1}{w}\frac{\partial\mathsf{Q}}{\partial\nu}&=\mathsf{Q}_{b}-\mathsf{Q},\\ {\bf v}&={\bf v}_{b}.\end{array}\right.

(37)

In the above, $\mathsf{Q}_{*},{\bf v}_{*}$ are the far-field states for $\mathsf{Q}$ and ${\bf v}$ , $\mathsf{Q}_{b}$ and ${\bf v}_{*}$ are to be specified, $w$ is some positive number, and $\nu$ is the outward unit normal for $\partial\Omega$ (or inward to $\mathbf{B}_{a}(0)$ ). We make the following remarks about the above boundary conditions.

Remark 1

The boundary condition for $\mathsf{Q}$ is associated with the following surface anchoring energy:

\mathcal{F}_{s}(\mathsf{Q})=\frac{w}{2}\int_{\partial\Omega}|\mathsf{Q}-\mathsf{Q}_{b}|^{2}\,d\sigma.

We will choose $\mathsf{Q}_{b}$ so as to have an explicit solution for $\mathsf{Q}$ – see Section 3.1.

Remark 2

The problem (29)-(31), (36)-(37) describes the flow of a nematic liquid crystal in the exterior of a colloidal particle under various scenarios. We emphasize that the quantity ${\bf v}_{b}$ is defined in the frame associated with the moving particle.

(a).

For a passive particle, the condition ${\bf v}_{b}=0$ describes no-slip boundary conditions on the surface of a particle that is stationary with respect to an inertial frame. The second condition in (36) imposes the constant velocity ${\bf v}_{*}$ of the flow at infinity.
(b).

Now suppose the passive particle moves in the nematic fluid with an externally imposed velocity ${\bf v}_{c}$ subject to the no-slip boundary condition, while the nematic is stationary at infinity, then the velocity on the boundary of the particle and at infinity equal ${\bf v}_{b}={\bf v}_{c}$ and ${\bf v}_{*}=0$ , respectively. In this case, if we go to a frame moving with the particle, then the velocity of the nematic liquid crystal at infinity will equal ${\bf v}_{*}=-{\bf v}_{c}$ and the velocity on the surface of the particle will vanish, that is, ${\bf v}_{b}=0$ . However, this change of frame will induce an additional forcing term in (33) due to the presence of the convective derivative ${\bf v}\cdot\nabla\mathsf{Q}$ of $\mathsf{Q}$ . This term will be described more explicitly in Section 4.1, in particular, equation (93).
(c).

For a general active particle, ${\bf v}_{b}$ is typically prescribed and nonconstant on the surface of the particle. Similar to the previous paragraph, if the particle moves with constant velocity ${\bf v}_{c}$ , then changing to a frame moving with the particle, we can replace ${\bf v}_{b}$ and ${\bf v}_{*}$ by ${\bf v}_{b}-{\bf v}_{c}$ and ${\bf v}_{*}=-{\bf v}_{c}$ . The value of ${\bf v}_{c}$ can be determined by solving the problem (29)-(31), (36)-(37) and choosing ${\bf v}_{c}$ so that the total stress on the surface of the particle vanishes.

From the mathematical point of view, in dimensons three or higher, the problem (29)-(31), (36)-(37) is uniquely solvable in exterior domains for any ${\bf v}_{b}$ and ${\bf v}_{*}$ . This is in contrast to the situation in bounded domain for which, due to the incompressibility condition, the total flux $\int{\bf v}_{b}\cdot\nu\,d\sigma$ at the boundary must vanish. In dimension two, in general, there is no solution in the exterior domain with a bounded velocity field. This is the origin of Stokes’ paradox.

As a final remark, we point out that it is certainly advantageous to consider our problem in the frame of the (moving) particle so that the domain does not change in time. Thus, in the case of passive particle which is the emphasis of the current paper, we set ${\bf v}_{b}=0$ and ${\bf v}_{*}$ to be some prescribed value. (As mentioned above, the extension to active particle in the current framework is achieved by setting ${\bf v}_{b}$ to be some general non-constant function. See [26] for examples of such a function.) Our goal is then to compute and analyze the flow pattern of the nematic fluid.

3 Preliminary information

We will make use of a known explicit stationary solution for $\mathsf{Q}$ in the exterior domain $\Omega$ and investigate its role in determining the flow pattern.

3.1 Stationary state for $\mathsf{Q}$ .

Under the physical regime and boundary conditions considered in the Sections 2.3 and 2.4, we are looking for a $\mathsf{Q}$ -tensor function satisfying

$\displaystyle\triangle\mathsf{Q}$	$\displaystyle=$	$\displaystyle 0,\quad\text{in $\Omega$}$
$\displaystyle\frac{1}{w}\frac{\partial\mathsf{Q}}{\partial\nu}$	$\displaystyle=$	$\displaystyle\mathsf{Q}_{b}-\mathsf{Q}\quad\text{on $\partial\Omega$,}$
$\displaystyle\mathsf{Q}$	$\displaystyle=$	$\displaystyle\mathsf{Q}_{*}\quad\text{at $\|x\|=\infty$.}$

The work Alama-Bronsard-Lamy [24, Theorem 1] gives the following explicit solution:

$\displaystyle\mathsf{Q}$	$\displaystyle=$	$\displaystyle\left(1-\frac{w}{1+w}\frac{1}{r}\right)\mathsf{Q}_{*}+\frac{w}{3+w}\frac{1}{r^{3}}\mathsf{Q}_{b},\quad r>1,$	(38)
$\displaystyle\text{where}\quad\mathsf{Q}_{*}$	$\displaystyle=$	$\displaystyle s_{}\left({n_{}}\otimes{n_{}}-\frac{\mathsf{I}}{3}\right),\quad\text{with a given ${n_{}}\in\mathbb{S}^{2}$}$	(39)
$\displaystyle\text{and}\quad\mathsf{Q}_{b}$	$\displaystyle=$	$\displaystyle s_{*}\left(\widehat{x}\otimes\widehat{x}-\frac{\mathsf{I}}{3}\right),\quad\widehat{x}=\frac{x}{\|x\|},$	(40)

with parameters $s_{*},w>0$ and the particle radius $a$ is taken to be one. Note that $\mathsf{Q}$ is harmonic and the boundary function $\mathsf{Q}_{b}$ is of “hedgehog” type. The function $\mathsf{Q}$ has the following far-field spatial asymptotics

\mathsf{Q}\sim\mathsf{Q}_{*}+O\left(\frac{1}{r}\right)\,\,\,\text{and}\,\,\,\nabla\mathsf{Q}\sim O\left(\frac{1}{r^{2}}\right),\,\,\,\text{for $r\gg 1$}.

(41)

We remark that [24] derives the above equation in the small particle regime $a^{2}\ll K$ , corresponding to $\varepsilon\gg 1$ . For large particle, $a^{2}\gg K$ , corresponding to $\varepsilon\ll 1$ , the solution $\mathsf{Q}$ tends to a harmonic map into ${\mathbb{S}}^{2}$ taking the form, $\mathsf{Q}(x)=s_{*}\left(n(x)\otimes n(x)-\frac{\mathsf{I}}{3}\right)$ [24, Theorem 2]. From the work [27], it is also shown that $n(x)$ has comparable spatial decay as the harmonic $\mathsf{Q}$ .

3.2 Green’s function for classical isotropic Stokes system

We introduce here the fundamental solution $(\mathsf{E},\mathsf{q})$ of the Stokes system which solves the system of equations,

\left\{\begin{array}[]{rcl}-\Delta\mathsf{E}_{ij}(x)-\frac{\partial}{\partial x_{i}}\mathsf{q}_{j}(x)&=&\delta_{ij}\delta_{0}(x),\\ \frac{\partial}{\partial x_{i}}\mathsf{E}_{ij}(x)&=&0.\end{array}\right.

(42)

Following [16, Chapter 4.2, p. 238], $(\mathsf{E},\mathsf{q})$ is given as:

	$\displaystyle\mathsf{E}(x)$	$\displaystyle=$	$\displaystyle\frac{1}{8\pi}\left[\frac{\mathsf{I}}{r}+\frac{x\otimes x}{r^{3}}\right],\quad\text{i.e.}\quad\mathsf{E}_{ij}(x)=\frac{1}{8\pi}\left[\frac{\delta_{ij}}{r}+\frac{x_{i}x_{j}}{r^{3}}\right],\,\,\,i,j=1,2,3,$		(43)
	$\displaystyle\mathsf{q}(x)$	$\displaystyle=$	$\displaystyle-\frac{1}{4\pi}\frac{x}{r^{3}},\quad\text{i.e.}\quad\mathsf{q}_{i}(x)=-\frac{x_{i}}{r^{3}},\,\,\,i=1,2,3.$		(44)

Note that

\mathsf{E}(x)\sim\frac{1}{r},\quad\text{and}\quad\mathsf{q}(x)\sim\frac{1}{r^{2}}\quad\text{for}\quad r\gg 1.

The above fundamental solution can be used to produce solutions of Stokes system on the whole $\mathbb{R}^{3}$ . More precisely, if $({\bf u},p)$ solves

$\displaystyle-\Delta{\bf u}+\nabla p$	$\displaystyle=$	$\displaystyle{\bf f},\quad\text{on $\mathbb{R}^{3}$},$
$\displaystyle\mathrm{div}\,{\bf u}$	$\displaystyle=$	$\displaystyle 0,\quad\text{on $\mathbb{R}^{3}$},$
$\displaystyle{\bf u}$	$\displaystyle=$	$\displaystyle 0,\quad\text{at $\|x\|=\infty$},$

then it is given by:

{\bf u}(x)=\int_{\mathbb{R}^{3}}\mathsf{E}(x-y){\bf f}(y)\,dy\quad\text{and}\quad p(x)=\int_{\mathbb{R}^{3}}-\mathsf{q}(x-y)\cdot{\bf f}(y)\,dy.

(45)

The above integrals are well-defined for ${\bf f}$ with sufficient spatial decay, for example, ${\bf f}(x)\lesssim 1\wedge\frac{1}{x^{3}}$ . For general existence theorems in $L^{p}(\mathbb{R}^{3})$ , we refer to [28, 16].

Similarly, in the case of a Stokes system in an exterior domain, for example,

$\displaystyle-\Delta{\bf u}+\nabla p$	$\displaystyle=$	$\displaystyle{\bf f},\quad\text{on $\Omega$},$
$\displaystyle\mathrm{div}\,{\bf u}$	$\displaystyle=$	$\displaystyle 0,\quad\text{on $\Omega$,}$
$\displaystyle{\bf u}$	$\displaystyle=$	$\displaystyle 0,\quad\text{at $\|x\|=\infty$},$

then ${\bf u}$ and $p$ can be represented by:

$\displaystyle{\bf u}(x)$	$\displaystyle=$	$\displaystyle\int_{\Omega}\mathsf{E}(x-y){\bf f}(y)\,dy$	(46)
		$\displaystyle+\int_{\partial\Omega}\big{\langle}\mathsf{E}(x-y),\mathsf{T}({\bf u},p)(y)\nu_{y}\big{\rangle}\,d\sigma_{y}-\int_{\partial\Omega}\big{\langle}\mathsf{T}(\mathsf{E},\mathsf{q})(x-y)\nu_{y},{\bf u}(y)\big{\rangle}\,d\sigma_{y},$	(46)
$\displaystyle p(x)$	$\displaystyle=$	$\displaystyle\int_{\Omega}-\mathsf{q}(x-y)\cdot{\bf f}(y)\,dy$	(47)
		$\displaystyle-\int_{\partial\Omega}\big{\langle}\mathsf{T}({\bf u},p)(y)\nu_{y},\mathsf{q}(x-y)\big{\rangle}\,d\sigma_{y}+2\int_{\partial\Omega}\big{\langle}\nabla_{x}\mathsf{q}(x-y)\nu_{y},{\bf u}(y)\big{\rangle}\,d\sigma_{y}.$	(47)

In the above, for any given vector and scalar fields $\bf w$ and $\pi$ , we define the stress tensor as

\mathsf{T}({\bf w},\pi)=\nabla{\bf w}+(\nabla{\bf w})^{T}-\pi\mathsf{I},

(48)

and for any matrix and vector fields $\mathbf{M}$ and $\mathbf{g}$ and vector $v$ , $\mathsf{T}(\mathbf{M},\mathbf{g})v$ is a matrix with its $i$ -th row given by

\big{(}\mathsf{T}(\mathbf{M},\mathbf{g})v\big{)}_{i}=\mathsf{T}(\mathbf{M}_{i},\mathbf{g}_{i})v.

We also recall the convention as stated in item (i) at the end of the Introduction.

In the above and the bulk of this paper, we are dealing with solutions ${\bf u}$ that converge to their far-field limits with rate $\displaystyle\frac{1}{r}$ . With this in mind, we expect the following estimates for the boundary integrals,

$\displaystyle\left\|\int_{\partial\Omega}\big{\langle}\mathsf{E}(x,-y),\mathsf{T}({\bf u},p)(y)\nu_{y}\big{\rangle}\,d\sigma_{y}\right\|$	$\displaystyle\lesssim$	$\displaystyle\frac{1}{r},$	(49)
$\displaystyle\left\|\int_{\partial\Omega}\big{\langle}\mathsf{T}(\mathsf{E},\mathsf{q}_{j})(x-y)\nu_{y},{\bf u}(y)\big{\rangle}\,d\sigma_{y}\right\|$	$\displaystyle\lesssim$	$\displaystyle\frac{1}{r^{2}},$	(50)
$\displaystyle\left\|\int_{\partial\Omega}\big{\langle}\mathsf{T}({\bf u},p)(y)\nu_{y},\mathsf{q}(x-y)\big{\rangle}\,d\sigma_{y}\right\|$	$\displaystyle\lesssim$	$\displaystyle\frac{1}{r^{2}},$	(51)
$\displaystyle\left\|\int_{\partial\Omega}\big{\langle}\nabla_{x}\mathsf{q}(x-y)\nu_{y},{\bf u}(y)\big{\rangle}\,d\sigma_{y}\right\|$	$\displaystyle\lesssim$	$\displaystyle\frac{1}{r^{3}}.$	(52)

Clearly (49) gives the dominating far-field behavior. It can be decomposed as:

	$\displaystyle\int_{\partial\Omega}\big{\langle}\mathsf{E}(x-y),\mathsf{T}({\bf u},p)(y)\nu_{y}\big{\rangle}\,d\sigma_{y}$	(53)
$\displaystyle=$	$\displaystyle\Big{\langle}\mathsf{E}(x),\int_{\partial\Omega}\mathsf{T}({\bf u},p)(y)\nu_{y}d\sigma_{y}\Big{\rangle}+\int_{\partial\Omega}\big{\langle}\mathsf{E}(x-y)-\mathsf{E}(x),\mathsf{T}({\bf u},p)(y)\nu_{y}\big{\rangle}d\sigma_{y}$
$\displaystyle=$	$\displaystyle\langle\mathsf{E}(x),{\mathcal{F}}\rangle+O\left(\frac{1}{r^{2}}\right)$

where

{\mathcal{F}}:=\int_{\partial\Omega}\mathsf{T}({\bf u},p)(y)\nu_{y}d\sigma_{y}

(54)

denotes the boundary stress or drag force on the particle. The decomposition (53) is the same as [29, Theorem 1, Eq. (4.2a)].

Precise far-field asymptotics of the bulk integrals and the method of solution for the Stokes system (32)–(34) will be given in Sections 6 and Appendix A.

3.3 Uniform Stokes Flow

Here we provide the solution of a uniform Stokes flow ${\bf U}$ with far-field velocity ${\bf U}_{*}$ , passing a sphere of radius $a$ . It solves the following system of equations:

$\displaystyle-\Delta{\bf U}+\nabla p$	$\displaystyle=$	$\displaystyle 0,\quad\text{for $\|x\|>a$};$	(55)
$\displaystyle\mathrm{div}\,{\bf U}$	$\displaystyle=$	$\displaystyle 0,\quad\text{for $\|x\|>a$};$	(56)
$\displaystyle{\bf U}$	$\displaystyle=$	$\displaystyle 0,\quad\text{at $\|x\|=a$};$	(57)
$\displaystyle{\bf U}$	$\displaystyle=$	$\displaystyle{\bf U}_{*},\quad\text{at $\|x\|=\infty$}.$	(58)

Using the spherical coordinates (following the physicists’ convention) with $\theta$ being the polar angle (measured from the polar axis) and $\phi$ being the azimuthal angle (measured from the meridian plane), we can write the velocity flowing along the polar axis as ${\bf U}=(u_{r},u_{\theta},u_{\phi}=0)$ . Following [30, Section 7.2], we have

u_{r}=\frac{1}{r^{2}\sin\theta}\frac{\partial\Psi}{\partial\theta},\quad\text{and}\quad u_{\theta}=-\frac{1}{r\sin\theta}\frac{\partial\Psi}{\partial r},

(59)

where

$\displaystyle\Psi(r,\theta)$	$\displaystyle=f(r)\sin^{2}\theta,\quad f(r)=\frac{V}{4}\left(2r^{2}-3ar+\frac{a^{3}}{r}\right),\quad V=\|{\bf U}_{*}\|,$	(60)
$\displaystyle u_{r}$	$\displaystyle=\frac{V}{4}\left(2-\frac{3a}{r}+\frac{a^{3}}{r^{3}}\right)2\cos\theta,$	(61)
$\displaystyle u_{\theta}$	$\displaystyle=-\frac{V}{4}\left(4-\frac{3a}{r}-\frac{a^{3}}{r^{3}}\right)\sin\theta.$	(62)

We note that the form of $f$ is found by solving $\displaystyle\left(\frac{\partial^{2}}{\partial r^{2}}-\frac{2}{r^{2}}\right)^{2}f(r)=0$ . The method of finding $f$ in a bounded (annular) domain will be presented in Appendix D.

If the flow is in the direction of the $x_{1}$ -axis, i.e., ${\bf U}_{*}=Ve_{1}$ , then we have in Cartesian coordinates that

$\displaystyle{\bf U}$	$\displaystyle=$	$\displaystyle V\left({e_{1}}-\frac{3a}{4}\left(\frac{x_{1}x+r^{2}{e_{1}}}{r^{3}}\right)-\frac{a^{3}}{4}\left(\frac{r^{2}{e_{1}}-3x_{1}x}{r^{5}}\right)\right)$	(63)
	$\displaystyle=$	$\displaystyle V\left({e_{1}}-\frac{3a}{4}\left(\frac{2x_{1}^{2}+x_{2}^{2}+x_{3}^{2},x_{1}x_{2},x_{1}x_{3}}{r^{3}}\right)\right.$
		$\displaystyle\left.-\frac{a^{3}}{4}\left(\frac{-2x_{1}^{2}+x_{2}^{2}+x_{3}^{2},-3x_{1}x_{2},-3x_{1}x_{3}}{r^{5}}\right)\right)$
$\displaystyle p$	$\displaystyle=$	$\displaystyle p_{\infty}-\frac{3}{2}\frac{Va}{r^{2}}\cos\theta=p_{\infty}-\frac{3}{2}\frac{Va}{r^{2}}\frac{x_{1}}{\sqrt{x_{1}^{2}+x_{2}^{2}}}.$	(64)

More generally, in vector form, we have

	$\displaystyle{\bf U}$	$\displaystyle=$	$\displaystyle{\bf U}_{}-\frac{3a}{4}\left(\frac{x\otimes x{\bf U}_{}+r^{2}{\bf U}_{}}{r^{3}}\right)-\frac{a^{3}}{4}\left(\frac{r^{2}{\bf U}_{}-3x\otimes x{\bf U}_{*}}{r^{5}}\right)$		(65)
		$\displaystyle=$	$\displaystyle\left[\mathsf{I}-\frac{3a}{4}\left(\frac{x\otimes x+r^{2}\mathsf{I}}{r^{3}}\right)+\frac{a^{3}}{4}\left(\frac{3x\otimes x-r^{2}\mathsf{I}}{r^{5}}\right)\right]{\bf U}_{*}$		(65)

Making use of the Green’s function $\mathsf{E}$ (43) and upon introducing

\mathsf{F}(x)=\frac{3x\otimes x-r^{2}\mathsf{I}}{r^{5}},

(66)

we have

\displaystyle{\bf U}

\displaystyle=

\displaystyle{\bf U}_{*}-6\pi a\mathsf{E}(x){\bf U}_{*}+\frac{a^{3}}{4}\mathsf{F}(x){\bf U}_{*},

or more compactly,

{\bf U}=\mathsf{E}_{S}{\bf U}_{*},\,\,\,\text{where}\,\,\,\mathsf{E}_{S}=\mathsf{I}-6\pi a\mathsf{E}(x)+\frac{a^{3}}{4}\mathsf{F}(x).

(67)

We note the following asymptotics,

\mathsf{E}(x)\sim O\left(\frac{1}{r}\right),\,\,\,\mathsf{F}(x)\sim O\left(\frac{1}{r^{3}}\right),\,\,\,

so that

{\bf U}\sim{\bf U}_{*}+O\left(\frac{1}{r}\right),\,\,\,\text{and}\,\,\,\nabla{\bf U},\,\mathsf{A},\,\mathsf{W}\sim O\left(\frac{1}{r^{2}}\right),\,\,\,\text{for $r\gg 1$.}

(68)

Next we compute the drag force ${\mathcal{F}}$ (54) on the moving particle. For this purpose, the components of $\mathsf{T}({\bf u},p)$ on $\partial{\bf B}_{a}$ , expressed in spherical coordinates are:

$\displaystyle\mathsf{T}_{rr}$	$\displaystyle=$	$\displaystyle-p+2\frac{\partial u_{r}}{\partial r}=-p_{\infty}+\frac{3}{2}\frac{V}{a}\cos\theta,$
$\displaystyle\mathsf{T}_{r\theta}$	$\displaystyle=$	$\displaystyle r\frac{\partial}{\partial r}\left(\frac{u_{\theta}}{r}\right)+\frac{1}{r}\frac{\partial u_{r}}{\partial\theta}=-\frac{3}{2}\frac{U}{a}\sin\theta,$
$\displaystyle\mathsf{T}_{r\phi}$	$\displaystyle=$	$\displaystyle 0.$

With the above, $\mathcal{F}$ , in the direction of ${\bf U}_{*}$ , is given by

\displaystyle{\mathcal{F}}=\int_{0}^{2\pi}\int_{0}^{\pi}\Big{(}T_{rr}\cos\theta-\mathsf{T}_{r\theta}\sin\theta\Big{)}a^{2}\sin\theta d\theta d\phi=6\pi aV

(69)

or in vector form, written as,

{\mathcal{F}}=6\pi a{\bf U}_{*}

(70)

which is the celebrated Stokes Law.

Remark 3

The classical isotropic Stokeslet fundamental solution given by

{\bf v}_{s}=\mathsf{E}\boldsymbol{\alpha},\quad p_{s}=\mathsf{q}\cdot\boldsymbol{\alpha},

for any $\boldsymbol{\alpha}\in\mathbb{R}^{3}$ corresponds to Stokes flow driven by a point source of strength $\boldsymbol{\alpha}$ located at the origin. Note that the far-field behavior of the Stokeslet is the same as the leading order asymptotics of classical Stokes flow in the exterior of a spherical domain. See [31] for a explanation and application of such a notion.

4 Structure of the anisotropic Stokes equation (33)

Here we will start our analysis for the stationary system (32)–(34) with boundary conditions (36)–(37). We emphasize again the feature that the tensor field $\mathsf{Q}$ given by (38)–(40) acts as an inhomogeneous term for the Stokes equation (33) for ${\bf v}$ . The key is to understand the far-field spatial behavior ( $r\gg 1$ ) of ${\bf v}$ .

4.1 Decomposition of equation (33)

Here we decompose system (29)–(31) into a form amenable for asymptotic analysis. We first substitute form (38) for $\mathsf{Q}$ into (33) for ${\bf v}$ and analyze the resulting system.

(a).

Recall the decay property (41) of $\mathsf{Q}$ : $\mathsf{Q}=\mathsf{Q}_{*}+O\left(\frac{1}{r}\right)$ and $\nabla\mathsf{Q}=O\left(\frac{1}{r^{2}}\right)$ . We look for a solution ${\bf v}$ satisfying ${\bf v}={\bf v}_{*}+O\left(\frac{1}{r}\right)$ and $\nabla{\bf v}=O\left(\frac{1}{r^{2}}\right)$ . Hence, we expect

$\displaystyle\mathsf{A},\,\,\,\mathsf{W}$	$\displaystyle\sim$	$\displaystyle O\left(\frac{1}{r^{2}}\right),$
$\displaystyle\mathring{\mathsf{Q}}$	$\displaystyle=$	$\displaystyle{\bf v}\cdot\nabla\mathsf{Q}+\mathsf{Q}\mathsf{W}-\mathsf{W}\mathsf{Q}$
	$\displaystyle=$	$\displaystyle\underset{O\left(\frac{1}{r^{2}}\right)}{\underbrace{{\bf v}_{}\cdot\nabla\mathsf{Q}+\mathsf{Q}_{}\mathsf{W}-\mathsf{W}\mathsf{Q}_{}}}+\underset{O\left(\frac{1}{r^{3}}\right)}{\underbrace{({\bf v}-{\bf v}_{})\cdot\nabla\mathsf{Q}+(\mathsf{Q}-\mathsf{Q}_{})\mathsf{W}-\mathsf{W}(\mathsf{Q}-\mathsf{Q}_{})}}$
	$\displaystyle\sim$	$\displaystyle O\left(\frac{1}{r^{2}}\right).$

(b).

With the above, the right hand side of (33) becomes

\mathcal{D}:=-\left(\gamma_{1}\mathring{\mathsf{Q}}+\gamma_{2}\mathsf{A}+\frac{\gamma_{3}}{2}(\mathsf{A}\mathsf{Q}+\mathsf{Q}\mathsf{A})+\gamma_{9}(\mathsf{A}\cdot\mathsf{Q})\mathsf{Q}\right)\cdot\partial_{i}\mathsf{Q}{\mathbf{e}}_{i}\sim O\left(\frac{1}{r^{4}}\right).

(71)

(c).

Here we analyze the terms constituting the viscous stress $\mathsf{T}^{\mathrm{v}}_{\text{SV}}$ given by (35).

(a)

The $\gamma_{1}$ -term:

	$\displaystyle\left[\mathsf{Q}\mathring{\mathsf{Q}}-\mathring{\mathsf{Q}}\mathsf{Q}\right]$	$\displaystyle=$	$\displaystyle\left[\mathsf{Q}\left({\bf v}\cdot\nabla\mathsf{Q}+\mathsf{Q}\mathsf{W}-\mathsf{W}\mathsf{Q}\right)-\left({\bf v}\cdot\nabla\mathsf{Q}+\mathsf{Q}\mathsf{W}-\mathsf{W}\mathsf{Q}\right)\mathsf{Q}\right]$		(72)
		$\displaystyle=$	$\displaystyle\left[\mathsf{Q}\left({\bf v}\cdot\nabla\mathsf{Q}\right)-\left({\bf v}\cdot\nabla\mathsf{Q}\right)\mathsf{Q}\right]+\left[\mathsf{Q}^{2}\mathsf{W}-2\mathsf{Q}\mathsf{W}\mathsf{Q}+\mathsf{W}\mathsf{Q}^{2}\right].$		(72)

For the first bracketed term in (72), we have the following decomposition

$\displaystyle\left[\mathsf{Q}\left({\bf v}\cdot\nabla\mathsf{Q}\right)-\left({\bf v}\cdot\nabla\mathsf{Q}\right)\mathsf{Q}\right]$	$\displaystyle=$	$\displaystyle\left[\mathsf{Q}_{}({\bf v}_{}\cdot\nabla\mathsf{Q})-({\bf v}_{}\cdot\nabla\mathsf{Q})\mathsf{Q}_{}\right]$	(73)
		$\displaystyle+\left[\mathsf{Q}_{}({\bf v}-{\bf v}_{})\cdot\nabla\mathsf{Q}-({\bf v}-{\bf v}_{})\cdot(\nabla\mathsf{Q})\mathsf{Q}_{}\right]$
		$\displaystyle+\left[(\mathsf{Q}-\mathsf{Q}_{}){\bf v}\cdot\nabla\mathsf{Q}-({\bf v}\cdot\nabla\mathsf{Q})(\mathsf{Q}-\mathsf{Q}_{})\right]$
	$\displaystyle=$	$\displaystyle\left[\mathsf{Q}_{}({\bf v}_{}\cdot\nabla\mathsf{Q})-({\bf v}_{}\cdot\nabla\mathsf{Q})\mathsf{Q}_{}\right]+O\left(\frac{1}{r^{3}}\right)$
	$\displaystyle=:$	$\displaystyle\mathcal{A}_{1}(x)+O\left(\frac{1}{r^{3}}\right).$

For the second bracketed term in (72), we have

	$\displaystyle\left[\mathsf{Q}^{2}\mathsf{W}-2\mathsf{Q}\mathsf{W}\mathsf{Q}+\mathsf{W}\mathsf{Q}^{2}\right]$	$\displaystyle=$	$\displaystyle\left[\mathsf{Q}_{}^{2}\mathsf{W}-2\mathsf{Q}_{}\mathsf{W}\mathsf{Q}_{}+\mathsf{W}\mathsf{Q}^{2}_{}\right]+O\left(\frac{1}{r^{3}}\right)$		(74)
		$\displaystyle=:$	$\displaystyle\mathcal{B}_{1}[\nabla{\bf v}]+O\left(\frac{1}{r^{3}}\right).$		(74)

In the above, $\mathcal{A}_{1}$ and $\mathcal{B}_{1}$ are linear in ${\bf v}_{*}\cdot\nabla\mathsf{Q}$ and $\nabla{\bf v}$ respectively, both having coefficients depending only on $\mathsf{Q}_{*}$ . Furthermore both $\mathcal{A}_{1}$ and $\mathcal{B}_{1}$ decay as $r^{-2}$ . A similar structure exists for all the remaining terms in the stress tensor as it will become clear while we proceed through the rest of this computation.

(b)

The $\gamma_{2}$ -term:

			$\displaystyle\left[\mathring{\mathsf{Q}}+\mathsf{Q}\mathsf{A}-\mathsf{A}\mathsf{Q}\right]$
		$\displaystyle=$	$\displaystyle\left[{\bf v}\cdot\nabla\mathsf{Q}+\mathsf{Q}\nabla{\bf v}-(\nabla{\bf v})\mathsf{Q}\right]$
		$\displaystyle=$	$\displaystyle\left[{\bf v}_{}\cdot\nabla\mathsf{Q}+\mathsf{Q}_{}\nabla{\bf v}-(\nabla{\bf v})\mathsf{Q}_{}\right]+\left[({\bf v}-{\bf v}_{})\cdot\nabla\mathsf{Q}+(\mathsf{Q}-\mathsf{Q}_{})\nabla{\bf v}-\nabla{\bf v}(\mathsf{Q}-\mathsf{Q}_{})\right]$
		$\displaystyle=$	$\displaystyle\mathcal{A}_{2}(x)+\mathcal{B}_{2}[\nabla{\bf v}]+O\left(\frac{1}{r^{3}}\right),$

where

\mathcal{A}_{2}(x):={\bf v}_{*}\cdot\nabla\mathsf{Q}\,\,\mbox{ and }\,\,\mathcal{B}_{2}[\nabla{\bf v}]:=\mathsf{Q}_{*}\nabla{\bf v}-(\nabla{\bf v})\mathsf{Q}_{*}.

(75)

(c)

The $\gamma_{3}$ -term:

			$\displaystyle\frac{1}{2}\left(\mathsf{Q}\mathring{\mathsf{Q}}+\mathring{\mathsf{Q}}\mathsf{Q}+\mathsf{Q}^{2}\mathsf{A}-\mathsf{A}\mathsf{Q}^{2}\right)$
		$\displaystyle=$	$\displaystyle\frac{1}{2}\Big{(}\mathsf{Q}_{}\big{(}{\bf v}_{}\cdot\nabla\mathsf{Q}+\mathsf{Q}_{}\mathsf{W}-\mathsf{W}\mathsf{Q}_{}\big{)}+\big{(}{\bf v}_{}\cdot\nabla\mathsf{Q}+\mathsf{Q}_{}\mathsf{W}-\mathsf{W}\mathsf{Q}_{}\big{)}\mathsf{Q}_{}$
			$\displaystyle+\mathsf{Q}_{}^{2}\mathsf{A}-\mathsf{A}\mathsf{Q}_{}^{2}\Big{)}+O\left(\frac{1}{r^{3}}\right)$
		$\displaystyle=$	$\displaystyle\frac{1}{2}\left(\mathsf{Q}_{}({\bf v}_{}\cdot\nabla\mathsf{Q})+({\bf v}_{}\cdot\nabla\mathsf{Q})\mathsf{Q}_{}+\mathsf{Q}_{}^{2}(\nabla{\bf v})-(\nabla{\bf v})\mathsf{Q}_{}^{2}\right)+O\left(\frac{1}{r^{3}}\right)$
		$\displaystyle=$	$\displaystyle\mathcal{A}_{3}(x)+\mathcal{B}_{3}[\nabla{\bf v}]+O\left(\frac{1}{r^{3}}\right),$

where

	$\displaystyle\mathcal{A}_{3}(x)$	$\displaystyle:=\frac{1}{2}\left(\mathsf{Q}_{}({\bf v}_{}\cdot\nabla\mathsf{Q})+({\bf v}_{}\cdot\nabla\mathsf{Q})\mathsf{Q}_{}\right),$		(76)
	$\displaystyle\mathcal{B}_{3}[\nabla{\bf v}]$	$\displaystyle:=\frac{1}{2}\left(\mathsf{Q}_{}^{2}(\nabla{\bf v})-(\nabla{\bf v})\mathsf{Q}_{}^{2}\right).$		(76)

(d)

The $\gamma_{4}$ -term:

\mathsf{Q}\mathsf{A}+\mathsf{A}\mathsf{Q}=\mathsf{Q}_{*}\mathsf{A}+\mathsf{A}\mathsf{Q}_{*}+O\left(\frac{1}{r^{3}}\right)=:\mathcal{B}_{4}[\nabla{\bf v}]+O\left(\frac{1}{r^{3}}\right).

(77)

(e)

The $\gamma_{5}$ -term:

\mathsf{Q}^{2}\mathsf{A}+\mathsf{A}\mathsf{Q}^{2}=\mathsf{Q}_{*}^{2}\mathsf{A}+\mathsf{A}\mathsf{Q}_{*}^{2}+O\left(\frac{1}{r^{3}}\right)=:\mathcal{B}_{5}[\nabla{\bf v}]+O\left(\frac{1}{r^{3}}\right).

(78)

(f)

The $\gamma_{6}$ -term:

(\mathsf{A}\cdot\mathsf{Q})\mathsf{Q}=(\mathsf{A}\cdot\mathsf{Q}_{*})\mathsf{Q}_{*}+O\left(\frac{1}{r^{3}}\right)=:\mathcal{B}_{6}[\nabla{\bf v}]+O\left(\frac{1}{r^{3}}\right).

(79)

(g)

The $\gamma_{7}$ -term:

|\mathsf{Q}^{2}|\mathsf{A}=|\mathsf{Q}^{2}_{*}|\mathsf{A}+O\left(\frac{1}{r^{3}}\right)=:\mathcal{B}_{7}[\nabla{\bf v}]+O\left(\frac{1}{r^{3}}\right).

(80)

(h)

The $\gamma_{9}$ -term:

$\displaystyle(\mathring{\mathsf{Q}}\cdot\mathsf{Q})\mathsf{Q}$	$\displaystyle=$	$\displaystyle\Big{(}({\bf v}\cdot\nabla\mathsf{Q}+\mathsf{Q}\mathsf{W}-\mathsf{W}\mathsf{Q})\cdot\mathsf{Q}\Big{)}\mathsf{Q}=\Big{(}({\bf v}\cdot\nabla\mathsf{Q})\cdot\mathsf{Q}\Big{)}\mathsf{Q}$	(81)
	$\displaystyle=$	$\displaystyle\Big{(}({\bf v}_{}\cdot\nabla\mathsf{Q})\cdot\mathsf{Q}_{}\Big{)}\mathsf{Q}_{*}+O\left(\frac{1}{r^{3}}\right)$
	$\displaystyle=$	$\displaystyle\mathcal{A}_{9}(x)+O\left(\frac{1}{r^{3}}\right),$

where we have used the fact that $(\mathsf{Q}\mathsf{W}-\mathsf{W}\mathsf{Q})\cdot\mathsf{Q}=0$ as $(\mathsf{C}\mathsf{D}-\mathsf{D}\mathsf{C})\cdot\mathsf{C}=0$ for any $\mathsf{C},\mathsf{D}\in M^{3\times 3}$ with $\mathsf{C}$ symmetric.

(i)

The $\gamma_{10}$ -term:

	$\displaystyle(\mathsf{Q}^{2}\cdot\mathsf{A})\mathsf{Q}+(\mathsf{Q}\cdot\mathsf{A})\mathsf{Q}^{2}$	$\displaystyle=$	$\displaystyle(\mathsf{Q}_{}^{2}\cdot\mathsf{A})\mathsf{Q}_{}+(\mathsf{Q}_{}\cdot\mathsf{A})\mathsf{Q}^{2}_{}+O\left(\frac{1}{r^{3}}\right)$		(82)
		$\displaystyle=:$	$\displaystyle\mathcal{B}_{10}[\nabla{\bf v}]+O\left(\frac{1}{r^{3}}\right).$		(82)

(j)

The $\gamma_{11}$ -term:

\displaystyle|\mathsf{Q}|^{2}(\mathsf{A}\cdot\mathsf{Q})\mathsf{Q}

\displaystyle=

\displaystyle|\mathsf{Q}_{*}|^{2}(\mathsf{A}\cdot\mathsf{Q}_{*})\mathsf{Q}_{*}+O\left(\frac{1}{r^{3}}\right)=:\mathcal{B}_{11}[\nabla{\bf v}]+O\left(\frac{1}{r^{3}}\right).

(83)

From the above, we note that

\mathcal{A}_{4}=\mathcal{A}_{5}=\mathcal{A}_{6}=\mathcal{A}_{7}=\mathcal{A}_{10}=\mathcal{A}_{11}=0\,\,\text{and}\,\,\mathcal{B}_{9}=0.

(84)

Taking into account the items (a)-(c) above, we can write equation (33) as

-\left(\Delta{\bf v}+\text{div}\Big{[}\mathcal{B}_{\gamma}[\nabla{\bf v}]+\mathcal{A}_{\gamma}(x)+\mathcal{C}_{\gamma}(x)\Big{]}\right)+\nabla p=\mathcal{D}_{\gamma}(x)

(85)

where

$\displaystyle\mathcal{B}_{\gamma}[\nabla{\bf v}]$	$\displaystyle=$	$\displaystyle\mathcal{B}_{\gamma,\mathsf{Q}_{*}}[\nabla{\bf v}]:=\sum_{i=1,\neq 8,9}^{11}\gamma_{i}\mathcal{B}_{i}[\nabla{\bf v}],$	(86)
$\displaystyle\mathcal{A}_{\gamma}(x)$	$\displaystyle=$	$\displaystyle\mathcal{A}_{\gamma,\mathsf{Q}_{}}[{\bf v}_{}\cdot\nabla\mathsf{Q}]:=\sum_{i=1,2,3,9}\gamma_{i}\mathcal{A}_{i}(x)=O\left(\frac{1}{r^{2}}\right),$	(87)
$\displaystyle\mathcal{C}_{\gamma}(x)$	$\displaystyle=$	$\displaystyle\mathcal{C}_{\gamma}({\bf v},x):=\mathsf{T}_{\text{SV}}^{\text{v}}-\mathsf{A}-\mathcal{B}_{\gamma}[\nabla{\bf v}]-\mathcal{A}_{\gamma}(x)=O\left(\frac{1}{r^{3}}\right),$	(88)
$\displaystyle\mathcal{D}_{\gamma}(x)$	$\displaystyle=$	$\displaystyle\mathcal{D}_{\gamma}({\bf v},x):=-\left(\gamma_{1}\mathring{\mathsf{Q}}+\gamma_{2}\mathsf{A}+\frac{\gamma_{3}}{2}(\mathsf{A}\mathsf{Q}+\mathsf{Q}\mathsf{A})+\gamma_{9}(\mathsf{A}\cdot\mathsf{Q})\mathsf{Q}\right)\cdot\partial_{i}\mathsf{Q}{\mathbf{e}}_{i}$
	$\displaystyle=$	$\displaystyle O\left(\frac{1}{r^{4}}\right).$

In the above and what follows, we will use the symbol $\gamma$ to denote an expression or quantity that genuinely depends on $\gamma_{1},\ldots,\gamma_{11}$ . Furthermore, we set

|\gamma|:=\max\{|\gamma_{1}|,\ldots,|\gamma_{11}|\}.

(90)

We might omit $\gamma$ if the dependence is clear from the context.

The advantage of representation (85)–(4.1) is highlighted as follows.

(a).

The linear form $\mathcal{B}_{\gamma}[\cdot]$ in $\nabla{\bf v}$ corresponds to the leading order far-field contribution to the diffusivity matrix originating from the stress tensor. In particular, $\mathrm{div}\,[\mathcal{B}_{\gamma}[\nabla v]]$ is linear in $D^{2}{\bf v}$ with constant coefficients depending only on $\mathsf{Q}_{*}$ .
(b).

The term $\mathrm{div}\,[\mathcal{A}_{\gamma}(x)]$ decays as $r^{-3}$ for $r\gg 1$ . It can be treated as a purely inhomogeneous forcing term involving $D^{2}\mathsf{Q}$ , ${\bf v}_{*}$ and $\mathsf{Q}_{*}$ . Note that it is linear in the expression ${\bf v}_{*}\cdot\nabla\mathsf{Q}$ so it vanishes if ${\bf v}_{*}=0$ .

(c).

The terms $\mathrm{div}\,[\mathcal{C}_{\gamma}(x)]$ and $\mathcal{D}_{\gamma}(x)$ decay as $r^{-4}$ for $r\gg 1$ . They are integrable in the exterior domain $\Omega$ ,

\int_{\Omega}\big{|}\mathrm{div}\,[\mathcal{C}_{\gamma}(x)]\big{|}\,dx,\,\,\,\int_{\Omega}\big{|}\mathcal{D}_{\gamma}(x)\big{|}\,dx<\infty.

(91)

We point out that the presence of $\mathcal{A}_{\gamma}$ and $\mathcal{D}_{\gamma}$ is due to the dependence of the interacting potential function $R$ on $\mathring{\mathsf{Q}}$ so that only $\gamma_{1},\gamma_{2},\gamma_{3}$ and $\gamma_{9}$ appear in the expressions for $\mathcal{A}_{\gamma}$ and $\mathcal{D}_{\gamma}$ . (See the form (12) of $R$ .)

For the purpose of analyzing equation (85), we will keep $\mathrm{div}\,\big{[}\mathcal{B}_{\gamma}[\nabla{\bf v}]\big{]}$ in the left hand side but move $\mathcal{A}_{\gamma},\mathcal{C}_{\gamma},\mathcal{D}_{\gamma}$ to the right of that equation and re-write it as:

-\Delta{\bf v}-\mathrm{div}\,\big{[}\mathcal{B}_{\gamma}[\nabla{\bf v}]\big{]}+\nabla p=\text{div}\big{[}\mathcal{A}_{\gamma}(x)\big{]}+\text{div}\big{[}\mathcal{C}_{\gamma}(x)\big{]}+\mathcal{D}_{\gamma}(x).

(92)

Note that $-\Delta{\bf v}-\mathrm{div}\,\big{[}\mathcal{B}_{\gamma}[\nabla{\bf v}]\big{]}$ is a second order differential operator in ${\bf v}$ with constant coefficients.

Before proceeding further, we will express $\mathcal{C}_{\gamma}$ and $\mathcal{D}_{\gamma}$ in more transparent forms.

(a).

For $\mathcal{C}_{\gamma}$ , from (88), we compute

	$\displaystyle\mathcal{C}_{\gamma}\,\,=\,\,$	$\displaystyle\mathsf{T}_{\text{SV}}^{\bf v}-\mathsf{A}-\mathcal{B}_{\gamma}[\nabla{\bf v}]]-\mathcal{A}_{\gamma}(x)$
	$\displaystyle\,\,=\,\,$	$\displaystyle\gamma_{1}\Big{[}\big{(}\mathsf{Q}\left({\bf v}\cdot\nabla\mathsf{Q}\right)-\left({\bf v}\cdot\nabla\mathsf{Q}\right)\mathsf{Q}\big{)}+\big{(}\mathsf{Q}^{2}\mathsf{W}-2\mathsf{Q}\mathsf{W}\mathsf{Q}+\mathsf{W}\mathsf{Q}^{2}\big{)}$
		$\displaystyle\hskip 15.0pt-\big{(}\mathsf{Q}_{}({\bf v}_{}\cdot\nabla\mathsf{Q})-({\bf v}_{}\cdot\nabla\mathsf{Q})\mathsf{Q}_{}\big{)}-\big{(}\mathsf{Q}_{}^{2}\mathsf{W}-2\mathsf{Q}_{}\mathsf{W}\mathsf{Q}_{}+\mathsf{W}\mathsf{Q}^{2}_{}\big{)}\Big{]}$
		$\displaystyle+\gamma_{2}\Big{[}\big{(}{\bf v}-{{\bf v}_{}}\big{)}\cdot\nabla\mathsf{Q}+\big{(}\mathsf{Q}-\mathsf{Q}_{}\big{)}\nabla{\bf v}-(\nabla{\bf v})\big{(}\mathsf{Q}-\mathsf{Q}_{*}\big{)}\Big{]}$
		$\displaystyle+\frac{\gamma_{3}}{2}\Big{[}\left(\mathsf{Q}({\bf v}\cdot\nabla\mathsf{Q})+({\bf v}\cdot\nabla\mathsf{Q})\mathsf{Q}+\mathsf{Q}^{2}\nabla{\bf v}-(\nabla{\bf v})\mathsf{Q}^{2}\right)$
		$\displaystyle\hskip 25.0pt-\left(\mathsf{Q}_{}({\bf v}_{}\cdot\nabla\mathsf{Q})+({\bf v}_{}\cdot\nabla\mathsf{Q})\mathsf{Q}_{}\right)-\left(\mathsf{Q}_{}^{2}(\nabla{\bf v})-(\nabla{\bf v})\mathsf{Q}_{}^{2}\right)\Big{]}$
		$\displaystyle+\gamma_{4}\Big{[}(\mathsf{Q}-\mathsf{Q}_{})\mathsf{A}+\mathsf{A}(\mathsf{Q}-\mathsf{Q}_{})\Big{]}$
		$\displaystyle+\gamma_{5}\Big{[}(\mathsf{Q}^{2}-\mathsf{Q}_{}^{2})\mathsf{A}+\mathsf{A}(\mathsf{Q}^{2}-\mathsf{Q}_{}^{2})\Big{]}$
		$\displaystyle+\gamma_{6}\Big{[}(\mathsf{A}\cdot\mathsf{Q})\mathsf{Q}-(\mathsf{A}\cdot\mathsf{Q}_{})\mathsf{Q}_{}\Big{]}$
		$\displaystyle+\gamma_{7}\Big{[}\left(\|\mathsf{Q}^{2}\|-\|\mathsf{Q}^{2}_{*}\|\right)\mathsf{A}\Big{]}$
		$\displaystyle+\gamma_{9}\Big{[}\big{(}({\bf v}\cdot\nabla\mathsf{Q})\cdot\mathsf{Q}\big{)}\mathsf{Q}-\Big{(}({\bf v}_{}\cdot\nabla\mathsf{Q})\cdot\mathsf{Q}_{}\Big{)}\mathsf{Q}_{*}\Big{]}$
		$\displaystyle+\gamma_{10}\Big{[}(\mathsf{Q}^{2}\cdot\mathsf{A})\mathsf{Q}+(\mathsf{Q}\cdot\mathsf{A})\mathsf{Q}^{2}-(\mathsf{Q}_{}^{2}\cdot\mathsf{A})\mathsf{Q}_{}-(\mathsf{Q}_{}\cdot\mathsf{A})\mathsf{Q}^{2}_{}\Big{]}$
		$\displaystyle+\gamma_{11}\Big{[}\|\mathsf{Q}\|^{2}(\mathsf{A}\cdot\mathsf{Q})\mathsf{Q}-\|\mathsf{Q}_{}\|^{2}(\mathsf{A}\cdot\mathsf{Q}_{})\mathsf{Q}_{*}\Big{]}$
	$\displaystyle\,\,\sim\,\,$	$\displaystyle O\left(\frac{1}{r^{3}}\right).$

(b).

For $\mathcal{D}_{\gamma}$ , from (71), we have,

$\displaystyle\mathcal{D}_{\gamma}$	$\displaystyle=$	$\displaystyle-\Big{[}\gamma_{1}\big{(}{\bf v}_{}\cdot\nabla\mathsf{Q}+\mathsf{Q}_{}\mathsf{W}-\mathsf{W}\mathsf{Q}_{}+({\bf v}-{\bf v}_{})\cdot\nabla\mathsf{Q}+(\mathsf{Q}-\mathsf{Q}_{})\mathsf{W}-\mathsf{W}(\mathsf{Q}-\mathsf{Q}_{})\big{)}$
		$\displaystyle\hskip 12.0pt+\gamma_{2}\mathsf{A}+\frac{\gamma_{3}}{2}\big{(}\mathsf{A}\mathsf{Q}_{}+\mathsf{Q}_{}\mathsf{A}+\mathsf{A}(\mathsf{Q}-\mathsf{Q}_{})+(\mathsf{Q}-\mathsf{Q}_{})\mathsf{A}\big{)}$
		$\displaystyle\hskip 12.0pt+\gamma_{9}\big{(}(\mathsf{A}\cdot\mathsf{Q}_{}+\mathsf{A}\cdot(\mathsf{Q}-\mathsf{Q}_{}))\mathsf{Q}_{}+(\mathsf{A}\cdot\mathsf{Q}_{}+\mathsf{A}\cdot(\mathsf{Q}-\mathsf{Q}_{}))(\mathsf{Q}-\mathsf{Q}_{})\big{)}\Big{]}\cdot\partial_{i}\mathsf{Q}{\mathbf{e}}_{i}$
	$\displaystyle\sim$	$\displaystyle O\left(\frac{1}{r^{4}}\right).$

Note that given $\mathsf{Q}$ , both $\mathcal{C}_{\gamma}$ and $\mathcal{D}_{\gamma}$ are linear in ${\bf v}$ .

Combining the above with the explicit expression (38) of $\mathsf{Q}$ , the terms on the right hand side of (92) take the following forms:

$\displaystyle\mathrm{div}\,\mathcal{A}_{\gamma}(x)$	$\displaystyle=$	$\displaystyle\gamma\frac{{\bf F}(\hat{x})}{r^{3}}+O\left(\frac{1}{r^{4}}\right),$	(93)
$\displaystyle\mathrm{div}\,\mathcal{C}_{\gamma}(x)+\mathcal{D}_{\gamma}(x)$	$\displaystyle=$	$\displaystyle\gamma\left(\frac{{\bf G}(\hat{x})}{r^{4}}+\frac{{\bf H}(\hat{x}):({\bf v}-{\bf v}_{*})}{r^{3}}+\frac{{\bf I}(\hat{x}):D{\bf v}}{r^{2}}+\frac{{\bf J}(\hat{x}):D^{2}{\bf v}}{r}\right)$	(94)
		$\displaystyle+O\left(\frac{1}{r^{5}}\right)$	(94)

for some spatially bounded vector or tensor fields $\bf F,G,H,I,J$ defined on $\mathbb{S}^{2}$ which depend on $\mathsf{Q}$ and $\mathsf{Q}_{*}$ but not on ${\bf v}$ and ${\bf v}_{*}$ . We again recall the convention about the symbol $``:"$ stated in item (i) at the end of Introduction. The precise forms of the contractions between tensors are not too important for our analysis. The key is the homogeneity in $r$ leading to appropriate spatial decays. The presence of $\gamma$ on the right hand side refers to the fact that the terms are multiplied by $\gamma_{i}$ ’s. In particular, $\mathcal{A}_{\gamma},\mathcal{B}_{\gamma},\mathcal{C}_{\gamma},\mathcal{D}_{\gamma}$ are all bounded in magnitude by $O(|\gamma|)$ . Hence, if $|\gamma|=0$ , then (92) simply becomes the classical isotropic Stokes equation (55).

Relating back to Section 2.4, Remark 2, note that the system (32)–(34) with boundary conditions (36)–(37) is solvable for any ${\bf v}_{b}$ and ${\bf v}_{*}$ . From the form of $\mathcal{A}_{\gamma}$ written in (93), it seems the system becomes simpler by setting ${\bf v}_{*}=0$ . This can be achieved by a change of frame or simply consider the new vector field ${\bf v}^{\prime}={\bf v}-{\bf v}_{*}$ . However, due to the presence of the convective derivative ${\bf v}\cdot\nabla\mathsf{Q}$ of $\mathsf{Q}$ , either of these procedures will necessarily give rise to the term ${\bf v}_{*}\nabla\mathsf{Q}$ which at leading order is embedded in $\mathcal{A}_{\gamma}$ .

4.2 Identification of $\mathrm{div}\,\big{[}\mathcal{B}_{\gamma}[\nabla{\bf v}]\big{]}$

In this section, we express $\mathrm{div}\,\big{[}\mathcal{B}_{\gamma}[\nabla{\bf v}]\big{]}$ in (92) in more explicit form. More precisely, we will write

\mathrm{div}\,\big{[}\mathcal{B}_{\gamma}[\nabla{\bf v}]\big{]}=\mathcal{M}_{\gamma}:D^{2}{\bf v}

(95)

for some constant fourth order tensor $\mathcal{M}_{\gamma}$ . Upon introducing the coordinates $\mathcal{M}_{\gamma}=(\mathcal{M}_{i,j;k,l})$ , the right hand side of the above is understood as

\Big{(}\mathcal{M}_{\gamma}:D^{2}{\bf v}\Big{)}_{i}=\mathcal{M}_{i,j;k,l}\partial_{kl}{\bf v}_{j},\quad\text{for $i=1,2,3$.}

(96)

The main purpose of this section is to identify $\mathcal{M}_{i,j;k,l}$ explicitly. We note the following symmetry action of $\mathcal{M}_{i,j;k,l}$ with respect to $k$ and $l$ ,

\mathcal{M}_{i,j;k,l}\partial_{kl}{\bf v}_{j}=\mathcal{M}_{i,j;l,k}\partial_{kl}{\bf v}_{j}.

(97)

Furthermore, we emphasize that $\mathcal{M}_{\gamma}$ will only act on incompressible vector fields ${\bf v}$ : $\partial_{j}{\bf v}_{j}=0$ .

Before proceeding, using (39), we record that

{\mathsf{Q}_{*}}_{ij}=s_{*}\left({n_{*}}_{i}{n_{*}}_{j}-\frac{1}{3}\delta_{ij}\right),\,\,\,{\mathsf{Q}_{*}^{2}}_{ij}=\frac{s_{*}^{2}}{3}\left({n_{*}}_{i}{n_{*}}_{j}+\frac{1}{3}\delta_{ij}\right),\,\,\,|\mathsf{Q}_{*}|^{2}=\frac{2s_{*}^{2}}{3}.

(98)

With that, we compute.

(a).

$\mathrm{div}(\mathcal{B}_{1}[\nabla{\bf v}]):$

			$\displaystyle\mathrm{div}\,\left(\left[\mathsf{Q}_{}^{2}\mathsf{W}+\mathsf{W}\mathsf{Q}_{}^{2}-2\mathsf{Q}_{}\mathsf{W}\mathsf{Q}_{}\right]_{i}\right)=\partial_{j}\Big{[}\mathsf{Q}_{}^{2}\mathsf{W}+\mathsf{W}\mathsf{Q}_{}^{2}-2\mathsf{Q}_{}\mathsf{W}\mathsf{Q}_{}\Big{]}_{ij}$
		$\displaystyle=$	$\displaystyle\partial_{j}\Big{[}\frac{1}{2}\left(\mathsf{Q}_{}^{2}\right)_{ik}(\partial_{j}{\bf v}_{k}-\partial_{k}{\bf v}_{j})+\frac{1}{2}\left(\mathsf{Q}_{}^{2}\right)_{kj}\big{(}\partial_{k}{\bf v}_{i}-\partial_{i}{\bf v}_{k}\Big{)}-{\left(\mathsf{Q}_{}\right)}_{ik}(\partial_{l}{\bf v}_{k}-\partial_{k}{\bf v}_{l}){\left(\mathsf{Q}_{}\right)}_{lj}\Big{]}$
		$\displaystyle=$	$\displaystyle\frac{1}{2}\left(\mathsf{Q}_{}^{2}\right)_{ik}\partial_{jj}{\bf v}_{k}+\frac{1}{2}\left(\mathsf{Q}_{}^{2}\right)_{kj}\Big{(}\partial_{kj}{\bf v}_{i}-\partial_{ij}{\bf v}_{k}\Big{)}-{\left(\mathsf{Q}_{}\right)}_{ik}{\left(\mathsf{Q}_{}\right)}_{lj}\big{(}\partial_{jl}{\bf v}_{k}-\partial_{jk}{\bf v}_{l}\big{)}$
		$\displaystyle=$	$\displaystyle\left(\frac{1}{2}\delta_{k,l}\left(\mathsf{Q}_{}^{2}\right)_{ij}+\frac{1}{2}\delta_{i,j}\left(\mathsf{Q}_{}^{2}\right)_{kl}-\frac{1}{2}\delta_{i,l}\left(\mathsf{Q}_{}^{2}\right)_{jk}-{\left(\mathsf{Q}_{}\right)}_{ij}{\left(\mathsf{Q}_{}\right)}_{lk}+{\left(\mathsf{Q}_{}\right)}_{ik}{\left(\mathsf{Q}_{*}\right)}_{lj}\right)\partial_{kl}{\bf v}_{j}$
		$\displaystyle=:$	$\displaystyle\mathcal{M}^{1}_{i,j;k,l}\partial_{kl}{\bf v}_{j},$

where

$\displaystyle\mathcal{M}^{1}_{i,j;k,l}$	$\displaystyle=$	$\displaystyle s_{}^{2}\left[\frac{1}{6}\delta_{k,l}\left({n_{}}_{i}{n_{}}_{j}+\frac{1}{3}\delta_{i,j}\right)+\frac{1}{6}\delta_{i,j}\left({n_{}}_{k}{n_{}}_{l}+\frac{1}{3}\delta_{k,l}\right)-\frac{1}{6}\delta_{i,l}\left({n_{}}_{k}{n_{*}}_{j}+\frac{1}{3}\delta_{k,j}\right)\right.$	(99)
		$\displaystyle\left.-{\left({n_{}}_{i}{n_{}}_{j}-\frac{1}{3}\delta_{i,j}\right)}{\left({n_{}}_{k}{n_{}}_{l}-\frac{1}{3}\delta_{k,l}\right)}+{\left({n_{}}_{k}{n_{}}_{i}-\frac{1}{3}\delta_{k,i}\right)}{\left({n_{}}_{j}{n_{}}_{l}-\frac{1}{3}\delta_{j,l}\right)}\right]$
	$\displaystyle=$	$\displaystyle\frac{s_{}^{2}}{2}\left(\delta_{kl}{n_{}}_{i}{n_{}}_{j}+\delta_{ij}{n_{}}_{k}{n_{}}_{l}-\delta_{ik}{n_{}}_{l}{n_{*}}_{j}\right).$

In the above, we have used the symmetry property (97) of $\mathcal{M}$ and the incompressibility of ${\bf v}$ :

\delta_{il}{n_{*}}_{k}{n_{*}}_{j}\partial_{kl}{\bf v}_{j}=\delta_{ik}{n_{*}}_{l}{n_{*}}_{j}\partial_{kl}{\bf v}_{j}\,\,\,\text{and}\,\,\,\delta_{il}\delta_{kj}\partial_{kl}{\bf v}_{j}=\delta_{ik}\delta_{lj}\partial_{kl}{\bf v}_{j}=\partial_{ij}{\bf v}_{j}=0.

(100)

The above property will be used in several places in what follows.

(b).

$\mathrm{div}(\mathcal{B}_{2}[\nabla{\bf v}]):$

$\displaystyle\mathrm{div}\,\left[\mathsf{Q}_{}\nabla{\bf v}-(\nabla{\bf v})\mathsf{Q}_{}\right]_{i}$	$\displaystyle=$	$\displaystyle\partial_{j}\left[{\left(\mathsf{Q}_{}\right)}_{ik}\partial_{j}{\bf v}_{k}-\partial_{k}{\bf v}_{i}{\left(\mathsf{Q}_{}\right)}_{kj}\right]$
	$\displaystyle=$	$\displaystyle{\left(\mathsf{Q}_{}\right)}_{ik}\partial_{jj}{\bf v}_{k}-{\left(\mathsf{Q}_{}\right)}_{kj}\partial_{kj}{\bf v}_{i}$
	$\displaystyle=$	$\displaystyle\left(\delta_{k,l}{\left(\mathsf{Q}_{}\right)}_{ij}-\delta_{i,j}{\left(\mathsf{Q}_{}\right)}_{kl}\right)\partial_{kl}{\bf v}_{j}$
	$\displaystyle=:$	$\displaystyle\mathcal{M}^{2}_{i,j;k,l}\partial_{kl}{\bf v}_{j}.$

Here

\mathcal{M}^{2}_{i,j;k,l}=s_{*}\left(\delta_{k,l}{n_{*}}_{i}{n_{*}}_{j}-\delta_{i,j}{n_{*}}_{k}{n_{*}}_{l}\right).

(101)

(c).

$\mathrm{div}(\mathcal{B}_{3}[\nabla{\bf v}]):$

$\displaystyle\frac{1}{2}\mathrm{div}\,\left[\mathsf{Q}_{}^{2}\nabla{\bf v}-(\nabla{\bf v})\mathsf{Q}_{}^{2}\right]_{i}$	$\displaystyle=$	$\displaystyle\frac{1}{2}\left({\left(\mathsf{Q}_{}^{2}\right)}_{ik}\partial_{jj}{\bf v}_{k}-\left({\mathsf{Q}_{}^{2}}\right)_{kj}\partial_{kj}{\bf v}_{i}\right)$
	$\displaystyle=$	$\displaystyle\frac{1}{2}\left(\delta_{k,l}{\left(\mathsf{Q}_{}^{2}\right)}_{ij}-\delta_{i,j}\left({\mathsf{Q}_{}^{2}}\right)_{kl}\right)\partial_{kl}{\bf v}_{j}$
	$\displaystyle=:$	$\displaystyle\mathcal{M}^{3}_{i,j;k,l}\partial_{kl}{\bf v}_{j}.$

Here

\mathcal{M}^{3}_{i,j;k,l}=\frac{s_{*}^{2}}{6}\left(\delta_{k,l}{n_{*}}_{i}{n_{*}}_{j}-\delta_{i,j}{n_{*}}_{k}{n_{*}}_{l}\right).

(102)

(d).

$\mathrm{div}(\mathcal{B}_{4}[\nabla{\bf v}]):$

$\displaystyle\mathrm{div}\,\left[\mathsf{Q}_{}\mathsf{A}+\mathsf{A}\mathsf{Q}_{}\right]_{i}$	$\displaystyle=$	$\displaystyle\frac{1}{2}\left(\mathsf{Q}_{}\right)_{ik}\partial_{jj}{\bf v}_{k}+\frac{1}{2}\left(\partial_{ij}{\bf v}_{k}+\partial_{kj}{\bf v}_{i}\right)\left(\mathsf{Q}_{}\right)_{kj}$
	$\displaystyle=$	$\displaystyle\frac{1}{2}\left(\delta_{k,l}\left(\mathsf{Q}_{}\right)_{ij}+\delta_{i,l}\left(\mathsf{Q}_{}\right)_{jk}+\delta_{i,j}\left(\mathsf{Q}_{*}\right)_{kl}\right)\partial_{kl}{\bf v}_{j}$
	$\displaystyle=:$	$\displaystyle\mathcal{M}^{4}_{i,j;k,l}\partial_{kl}{\bf v}_{j}.$

Here

	$\displaystyle\mathcal{M}^{4}_{i,j;k,l}$	$\displaystyle=$	$\displaystyle\frac{s_{}}{2}\left(\delta_{k,l}\left({n_{}}_{i}{n_{}}_{j}-\frac{1}{3}\delta_{i,j}\right)+\delta_{i,l}\left({n_{}}_{j}{n_{}}_{k}-\frac{1}{3}\delta_{j,k}\right)+\delta_{i,j}\left({n_{}}_{k}{n_{*}}_{l}-\frac{1}{3}\delta_{k,l}\right)\right)$		(103)
		$\displaystyle=$	$\displaystyle\frac{s_{}}{2}\left(\delta_{kl}{n_{}}_{i}{n_{}}_{j}+\delta_{ij}{n_{}}_{k}{n_{}}_{l}+\delta_{ik}{n_{}}_{l}{n_{*}}_{j}-\frac{2}{3}\delta_{kl}\delta_{ij}\right).$		(103)

(e).

$\mathrm{div}(\mathcal{B}_{5}[\nabla{\bf v}]):$

$\displaystyle\mathrm{div}\,\left[\mathsf{Q}_{}^{2}\mathsf{A}+\mathsf{A}\mathsf{Q}_{}^{2}\right]_{i}$	$\displaystyle=$	$\displaystyle\frac{1}{2}\left(\mathsf{Q}_{}^{2}\right)_{ik}\partial_{jj}{\bf v}_{k}+\frac{1}{2}\left(\partial_{ij}{\bf v}_{k}+\partial_{kj}{\bf v}_{i}\right)\left(\mathsf{Q}_{}^{2}\right)_{kj}$
	$\displaystyle=$	$\displaystyle\frac{1}{2}\left(\delta_{k,l}\left(\mathsf{Q}_{}^{2}\right)_{ij}+\delta_{i,l}\left(\mathsf{Q}_{}^{2}\right)_{jk}+\delta_{i,j}\left(\mathsf{Q}_{*}^{2}\right)_{kl}\right)\partial_{kl}{\bf v}_{j}$
	$\displaystyle=:$	$\displaystyle\mathcal{M}^{5}_{i,j;k,l}\partial_{kl}{\bf v}_{j}.$

Here

	$\displaystyle\mathcal{M}^{5}_{i,j;k,l}$	$\displaystyle=$	$\displaystyle\frac{s_{}^{2}}{6}\left(\delta_{k,l}\left({n_{}}_{i}{n_{}}_{j}+\frac{1}{3}\delta_{i,j}\right)+\delta_{i,l}\left({n_{}}_{j}{n_{}}_{k}+\frac{1}{3}\delta_{j,k}\right)+\delta_{i,j}\left({n_{}}_{k}{n_{*}}_{l}+\frac{1}{3}\delta_{k,l}\right)\right)$		(104)
		$\displaystyle=$	$\displaystyle\frac{s_{}^{2}}{6}\left(\delta_{kl}{n_{}}_{i}{n_{}}_{j}+\delta_{ij}{n_{}}_{k}{n_{}}_{l}+\delta_{ik}{n_{}}_{l}{n_{*}}_{j}+\frac{2}{3}\delta_{ij}\delta_{kl}\right).$		(104)

(f).

$\mathrm{div}(\mathcal{B}_{6}[\nabla{\bf v}]):$

$\displaystyle\mathrm{div}\,\left[(\mathsf{A}\cdot\mathsf{Q}_{})\mathsf{Q}_{}\right]_{i}$	$\displaystyle=$	$\displaystyle\frac{1}{2}\left(\mathsf{Q}_{}\right)_{ij}\left(\mathsf{Q}_{}\right)_{kl}\left(\partial_{lj}{\bf v}_{k}+\partial_{kj}{\bf v}_{l}\right)$
	$\displaystyle=$	$\displaystyle\frac{1}{2}\left(\left(\mathsf{Q}_{}\right)_{ik}\left(\mathsf{Q}_{}\right)_{jl}+\left(\mathsf{Q}_{}\right)_{il}\left(\mathsf{Q}_{}\right)_{kj}\right)\partial_{kl}{\bf v}_{j}$
	$\displaystyle=:$	$\displaystyle\mathcal{M}^{6}_{i,j;k,l}\partial_{kl}{\bf v}_{j}.$

Here

$\displaystyle\mathcal{M}^{6}_{i,j;k,l}$	$\displaystyle=$	$\displaystyle\frac{s_{}^{2}}{2}\left(\left({n_{}}_{i}{n_{}}_{k}-\frac{1}{3}\delta_{i,k}\right)\left({n_{}}_{j}{n_{*}}_{l}-\frac{1}{3}\delta_{j,l}\right)\right.$	(105)
		$\displaystyle\left.+\left({n_{}}_{i}{n_{}}_{l}-\frac{1}{3}\delta_{i,l}\right)\left({n_{}}_{j}{n_{}}_{k}-\frac{1}{3}\delta_{j,k}\right)\right)$
	$\displaystyle=$	$\displaystyle s_{}^{2}\left({n_{}}_{i}{n_{}}_{j}{n_{}}_{k}{n_{}}_{l}-\frac{1}{3}\delta_{ik}{n_{}}_{l}{n_{*}}_{j}\right).$

(g).

$\mathrm{div}(\mathcal{B}_{7}[\nabla{\bf v}]):$

\displaystyle\mathrm{div}\,\left[|\mathsf{Q}_{*}|^{2}\mathsf{A}\right]_{i}=\mathrm{div}\,\left[\frac{2s_{*}^{2}}{3}\mathsf{A}\right]_{i}=:\mathcal{M}^{7}_{i,j;k,l}\partial_{kl}{\bf v}_{j}.

Here

\mathcal{M}^{7}_{i,j;k,l}=\frac{s_{*}^{2}}{3}\delta_{k,l}\delta_{i,j}.

(106)

(h).

$\mathrm{div}(\mathcal{B}_{10}[\nabla{\bf v}]):$

$\displaystyle\mathrm{div}\,\left[(\mathsf{Q}_{}^{2}\cdot\mathsf{A})\mathsf{Q}_{}+(\mathsf{Q}_{}\cdot\mathsf{A})\mathsf{Q}^{2}_{}\right]_{i}$	$\displaystyle=$	$\displaystyle\frac{1}{2}\left(\left(\mathsf{Q}_{}^{2}\right)_{kl}\left(\mathsf{Q}_{}\right)_{ij}+\left(\mathsf{Q}_{}\right)_{kl}\left(\mathsf{Q}_{}^{2}\right)_{ij}\right)\left(\partial_{lj}{\bf v}_{k}+\partial_{kj}{\bf v}_{l}\right)$
	$\displaystyle=$	$\displaystyle\left(\left(\mathsf{Q}_{}^{2}\right)_{jk}\left(\mathsf{Q}_{}\right)_{il}+\left(\mathsf{Q}_{}\right)_{jk}\left(\mathsf{Q}_{}^{2}\right)_{il}\right)\partial_{kl}{\bf v}_{j}$
	$\displaystyle=:$	$\displaystyle\mathcal{M}^{10}_{i,j;k,l}\partial_{kl}{\bf v}_{j}.$

Here

$\displaystyle\mathcal{M}^{10}_{i,j;k,l}$	$\displaystyle=$	$\displaystyle\frac{s_{}^{3}}{3}\left(\left({n_{}}_{j}{n_{}}_{k}+\frac{1}{3}\delta_{j,k}\right)\left({n_{}}_{i}{n_{*}}_{l}-\frac{1}{3}\delta_{i,l}\right)\right.$	(107)
		$\displaystyle\left.+\left({n_{}}_{j}{n_{}}_{k}-\frac{1}{3}\delta_{j,k}\right)\left({n_{}}_{i}{n_{}}_{l}+\frac{1}{3}\delta_{i,l}\right)\right)$
	$\displaystyle=$	$\displaystyle\frac{2s_{}^{3}}{3}{n_{}}_{i}{n_{}}_{j}{n_{}}_{k}{n_{*}}_{l}.$

(i).

$\mathrm{div}(\mathcal{B}_{11}[\nabla{\bf v}]):$

\displaystyle\mathrm{div}\,\left[|\mathsf{Q}_{*}|^{2}(\mathsf{A}\cdot\mathsf{Q}_{*})\mathsf{Q}_{*}\right]_{i}=\frac{2s_{*}^{2}}{3}\left(\mathsf{Q}_{*}\right)_{il}\left(\mathsf{Q}_{*}\right)_{jk}\partial_{kl}{\bf v}_{j}=:\mathcal{M}^{11}_{i,j;k,l}\partial_{kl}{\bf v}_{j}.

Here

	$\displaystyle\mathcal{M}^{11}_{i,j;k,l}$	$\displaystyle=$	$\displaystyle\frac{2s_{}^{4}}{3}\left({n_{}}_{i}{n_{}}_{l}-\frac{1}{3}\delta_{i,l}\right)\left({n_{}}_{j}{n_{*}}_{k}-\frac{1}{3}\delta_{j,k}\right)$		(108)
		$\displaystyle=$	$\displaystyle\frac{2s_{}^{4}}{3}\left({n_{}}_{i}{n_{}}_{j}{n_{}}_{k}{n_{}}_{l}-\frac{1}{3}\delta_{ik}{n_{}}_{l}{n_{*}}_{j}\right).$		(108)

With the above, the fourth order tensor $\mathcal{M}_{\gamma}$ can be decomposed as

\mathcal{M}_{i,j;k,l}=\sum_{p=1,\ldots,7,10,11}\gamma_{p}\mathcal{M}^{p}_{i,j;k,l}.

(109)

From (99) to (108), the $\mathcal{M}^{p}_{i,j;k,l}$ ’s are effectively given by a linear combination of the following tensors:

\delta_{kl}\delta_{ij},\,\,\,\delta_{kl}{n_{*}}_{i}{n_{*}}_{j},\,\,\,\delta_{ij}{n_{*}}_{k}{n_{*}}_{l},\,\,\,\delta_{ik}{n_{*}}_{l}{n_{*}}_{j}\,\,\,(\text{or equivalently, $\delta_{il}{n_{*}}_{k}{n_{*}}_{j}$}),\,\,\,{n_{*}}_{i}{n_{*}}_{j}{n_{*}}_{k}{n_{*}}_{l}.

(110)

5 Analysis of the anisotropic Stokes equation (33)–(34)

Collecting the forms from (93) and (94), we write again here the governing system as:

$\displaystyle\mathcal{L}_{\gamma}{\bf v}+\nabla p$	$\displaystyle=$	$\displaystyle{\bf f}_{\gamma}({\bf v}),\quad\text{for}\,\,\,\|x\|>a,$	(111)
$\displaystyle\quad\mathrm{div}\,{\bf v}$	$\displaystyle=$	$\displaystyle 0,\quad\text{for}\,\,\,\|x\|>a,$	(112)
$\displaystyle{\bf v}$	$\displaystyle=$	$\displaystyle{\bf v}_{b}\,\,(=0),\quad\text{on $\partial\mathbf{B}_{a}$}$	(113)
$\displaystyle{\bf v}$	$\displaystyle=$	$\displaystyle{\bf v}_{*},\quad\text{at $\|x\|=\infty$.}$	(114)

where

$\displaystyle\mathcal{L}_{\gamma}{\bf v}$	$\displaystyle:=$	$\displaystyle-\Big{(}\Delta{\bf v}+\mathcal{M}_{\gamma}:D^{2}{\bf v}\Big{)},$	(115)
$\displaystyle{\bf f}_{\gamma}({\bf u})$	$\displaystyle:=$	$\displaystyle\mathrm{div}\,\mathcal{A}_{\gamma}(x)+\mathrm{div}\,\mathcal{C}_{\gamma}({\bf u},x)+\mathcal{D}_{\gamma}({\bf u},x),$	(116)
	$\displaystyle=$	$\displaystyle\gamma\left(\frac{{\bf F}(\hat{x})}{r^{3}}+\frac{{\bf G}(\hat{x})}{r^{4}}+\frac{{\bf H}(\hat{x}):({\bf u}-{\bf v}_{*})}{r^{3}}+\frac{{\bf I}(\hat{x}):D{\bf u}}{r^{2}}+\frac{{\bf J}(\hat{x}):D^{2}{\bf u}}{r}\right)$
		$\displaystyle+O\left(\frac{1}{r^{5}}\right).$

Note that we have emphasized the dependence of $\mathcal{C}_{\gamma}$ and $\mathcal{D}_{\gamma}$ on ${\bf u}$ . Given $\mathsf{Q}$ , this dependence is in fact linear in ${\bf u}$ so that

{\bf f}_{\gamma}({\bf u}_{1})-{\bf f}_{\gamma}({\bf u}_{2})={\bf f}_{\gamma}({\bf u}_{1}-{\bf u}_{2})

(117)

making the above a linear system which we designate as our anisotropic Stokes system.

We remark again that the above system describes the flow in the moving frame attached to the particle. The far-field velocity ${\bf v}_{*}$ is prescribed. For passive particle which is the case in this paper, we take ${\bf v}_{b}=0$ while for active particle, it is in general some prescribed, non-constant function. Note also that if ${\bf v}_{b}=0$ and ${\bf v}_{*}=0$ , then we have only the trivial solution ${\bf v}=0$ .

Upon introducing

\mathsf{T}_{\gamma}({\bf w},\pi)=\nabla{\bf w}+(\nabla{\bf w})^{T}+\mathcal{B}_{\gamma}(\nabla{\bf v})-\pi\mathsf{I},

(118)

we can write (111) as

-\mathrm{div}\,\mathsf{T}_{\gamma}({\bf v},p)={\bf f}_{\gamma}({\bf v}).

(119)

In terms of $\mathcal{M}_{i,j;k,l}$ , we have $\big{(}\mathcal{B}_{\gamma}(\nabla{\bf v})\big{)}_{ij}=\mathcal{M}_{i,m;k,j}\partial_{k}{\bf v}_{m}$ so that it is consistent with $\big{(}\mathrm{div}\,\mathcal{B}_{\gamma}(\nabla{\bf v})\big{)}_{i}=\mathcal{M}_{i,j;k,l}\partial_{kl}{\bf v}_{j}$ .

The main purpose of the next few sections is to prove the existence and uniqueness of a solution for (111)–(114) in a suitable function space when $|\gamma|\ll 1$ .

5.1 Computation of the Green’s Function for $\mathcal{L}$ .

Given a function or vector field $f$ defined on $\mathbb{R}^{3}$ , its Fourier transform and inverse are given by

\widetilde{f}(\xi)=\int_{\mathbb{R}^{n}}f(x)e^{-\mathrm{i}x\cdot\xi}\,dx,\quad f(x)=(2\pi)^{-n}\int_{\mathbb{R}^{n}}\widetilde{f}(\xi)e^{\mathrm{i}x\cdot\xi}\,d\xi,

where $\xi=(\xi_{1},\xi_{2},\xi_{3})^{T}$ . Based on (95) and (96), we introduce the $3\times 3$ matrix $[\mathcal{M}_{\gamma}](\xi)$ :

[\mathcal{M}_{\gamma}]_{ij}(\xi)=\mathcal{M}_{i,j;k,l}\xi_{k}\xi_{l}.

(120)

Then taking the Fourier transform of (92) gives

{|\xi|}^{2}\widetilde{{\bf v}}+[\mathcal{M}_{\gamma}](\xi)\widetilde{{\bf v}}+\mathrm{i}\xi\widetilde{p}=\widetilde{{\bf f}_{\gamma}}

(121)

where ${\bf f}_{\gamma}$ is the right hand side of (92). Using the incompressibility condition $\mathrm{div}\,{\bf v}=0$ written in Fourier mode $\langle\widetilde{{\bf v}},\xi\rangle=0$ , we have,

\displaystyle\langle[\mathcal{M}_{\gamma}](\xi)\widetilde{{\bf v}},\xi\rangle+\mathrm{i}|\xi|^{2}\widetilde{p}=\langle\widetilde{{\bf f}_{\gamma}},\xi\rangle,\quad\text{so that}\quad\widetilde{p}=\frac{-\big{\langle}[\mathcal{M}_{\gamma}](\xi)\widetilde{{\bf v}},\xi\big{\rangle}+\langle\widetilde{{\bf f}_{\gamma}},\xi\rangle}{\mathrm{i}|\xi|^{2}}.

Hence,

\displaystyle{|\xi|}^{2}\widetilde{{\bf v}}+[\mathcal{M}_{\gamma}](\xi)\widetilde{{\bf v}}+\frac{-\big{\langle}[\mathcal{M}_{\gamma}](\xi)\widetilde{{\bf v}},\xi\big{\rangle}\xi+\langle\widetilde{{\bf f}_{\gamma}},\xi\rangle\xi}{|\xi|^{2}}=\widetilde{{\bf f}_{\gamma}},

i.e.

\displaystyle{|\xi|}^{2}\left(\widetilde{{\bf v}}+[\mathcal{M}_{\gamma}](\hat{\xi})\widetilde{{\bf v}}-\big{\langle}[\mathcal{M}_{\gamma}](\hat{\xi})\widetilde{{\bf v}},\hat{\xi}\big{\rangle}\hat{\xi}\right)

\displaystyle=

\displaystyle\widetilde{{\bf f}_{\gamma}}-\langle\widetilde{{\bf f}_{\gamma}},\hat{\xi}\rangle\hat{\xi},\,\,\,\hat{\xi}=\frac{\xi}{|\xi|},

or equivalently,

\displaystyle{|\xi|}^{2}\left(\widetilde{{\bf v}}+[\mathcal{M}_{\gamma}](\hat{\xi})\widetilde{{\bf v}}-(\hat{\xi}\otimes\hat{\xi})[\mathcal{M}_{\gamma}](\hat{\xi})\widetilde{{\bf v}}\right)

\displaystyle=

\displaystyle(\mathsf{I}-\xi\otimes\xi)\widetilde{{\bf f}_{\gamma}}.

To conclude, the solution ${\bf v}$ and $p$ is given by

\displaystyle\widetilde{{\bf v}}

\displaystyle=

\displaystyle\frac{1}{|\xi|^{2}}\Big{(}\mathsf{I}+(\mathsf{I}-\hat{\xi}\otimes\hat{\xi})[\mathcal{M}_{\gamma}](\hat{\xi})\Big{)}^{-1}(\mathsf{I}-\hat{\xi}\otimes\hat{\xi})\widetilde{{\bf f}_{\gamma}}

(122)

and

\displaystyle\widetilde{p}=\frac{1}{\mathrm{i}|\xi|}\left\langle\left(\mathsf{I}-(\mathsf{I}-\hat{\xi}\otimes\hat{\xi})\Big{(}\mathsf{I}+(\mathsf{I}-\hat{\xi}\otimes\hat{\xi})[\mathcal{M}_{\gamma}](\hat{\xi})\Big{)}^{-T}[\mathcal{M}_{\gamma}]^{T}(\hat{\xi})\right)\hat{\xi},\widetilde{{\bf f}_{\gamma}}\right\rangle.

(123)

Note that the matrix inverse in the above is well-defined if $|\gamma|\ll 1$ .

Let ${\mathsf{G}_{\gamma}}(x)$ and ${\mathsf{h}_{\gamma}}(x)$ be the inverse Fourier transforms respectively of

\frac{1}{|\xi|^{2}}\Big{(}\mathsf{I}+(\mathsf{I}-\hat{\xi}\otimes\hat{\xi})[\mathcal{M}_{\gamma}](\hat{\xi})\Big{)}^{-1}(\mathsf{I}-\hat{\xi}\otimes\hat{\xi})

(124)

and

\frac{1}{\mathrm{i}|\xi|}\left(\mathsf{I}-(\mathsf{I}-\hat{\xi}\otimes\hat{\xi})\Big{(}\mathsf{I}+(\mathsf{I}-\hat{\xi}\otimes\hat{\xi})[\mathcal{M}_{\gamma}](\hat{\xi})\Big{)}^{-T}[\mathcal{M}_{\gamma}]^{T}(\hat{\xi})\right)\hat{\xi}.

(125)

Then we have that ${\mathsf{G}_{\gamma}}$ and ${\mathsf{h}_{\gamma}}$ are homogeneous with degrees $-1$ and $-2$ so that

{\mathsf{G}_{\gamma}}(x)=\frac{{\mathsf{G}_{\gamma}}(\hat{x})}{|x|}\,\,\,\text{and}\,\,\,{\mathsf{h}_{\gamma}}(x)=\frac{{\mathsf{h}_{\gamma}}(\hat{x})}{|x|^{2}}.

(126)

The above leads to the following properties: for $k=0,1,2\ldots$ ,

(D_{x}^{k}{\mathsf{G}_{\gamma}})(x)=\frac{(D_{x}^{k}{\mathsf{G}_{\gamma}})(\hat{x})}{|x|^{k+1}},\,\,\,\Big{|}D^{k}_{x}{\mathsf{G}_{\gamma}}(x)\Big{|}\lesssim\frac{1}{|x|^{k+1}},

(127)

and

(D_{x}^{k}{\mathsf{h}_{\gamma}})(x)=\frac{(D_{x}^{k}{\mathsf{h}_{\gamma}})(\hat{x})}{|x|^{k+2}},\,\,\,\Big{|}D^{k}_{x}{\mathsf{h}_{\gamma}}(x)\Big{|}\lesssim\frac{1}{|x|^{k+2}}.

(128)

Using the ${\mathsf{G}_{\gamma}}$ and ${\mathsf{h}_{\gamma}}$ above, we are looking for a solution ${\bf v}$ of (111) with suitable decay property at infinity such that the following representation holds for ${\bf v}$ :

	$\displaystyle{\bf v}(x)$	$\displaystyle=$	$\displaystyle\int_{\Omega}{\mathsf{G}_{\gamma}}(x-y){\bf f}_{\gamma}({\bf v}(y))\,dy+\int_{\partial\Omega}\big{\langle}{\mathsf{G}_{\gamma}}(x-y),\mathsf{T}_{\gamma}({\bf v},p)(y)\nu_{y}\big{\rangle}\,d\sigma_{y}$		(129)
			$\displaystyle-\int_{\partial\Omega}\big{\langle}\mathsf{T}_{\gamma}({\mathsf{G}_{\gamma}},{\mathsf{h}_{\gamma}})(x-y)\nu_{y},{\bf v}(y)\big{\rangle}\,d\sigma_{y}+{\bf v}_{*}.$		(129)

Similar to (49)–(52), we have

\left|\int_{\partial\Omega}\big{\langle}{\mathsf{G}_{\gamma}}(x-y),\mathsf{T}_{\gamma}({\bf v},p)(y)\nu_{y}\big{\rangle}\,d\sigma_{y}\right|\lesssim\frac{1}{r},\,\,\,\text{and}\,\,\,\left|\int_{\partial\Omega}\big{\langle}\mathsf{T}_{\gamma}({\mathsf{G}_{\gamma}},{\mathsf{h}_{\gamma}})(x-y)\nu_{y},{\bf v}(y)\big{\rangle}\,d\sigma_{y}\right|\lesssim\frac{1}{r^{2}},

(130)

so that ${\bf v}-{\bf v}_{*}\approx r^{-1}$ . More precise far-field asymptotics of ${\bf v}$ will be given in Section 6. However, in preparation for the proof of existence of ${\bf v}$ , we will first analyze the dominating term of ${\bf f}_{\gamma}({\bf v})$ given by $\mathrm{div}\,\mathcal{A}_{\gamma}(x)$ .

5.2 Far-field behavior of the inhomogeneous term

We recall form (93) for the dominating term of ${\bf f}_{\gamma}({\bf v})$ :

\mathrm{div}\,\mathcal{A}_{\gamma}(x)=\gamma\frac{{\bf F}(\hat{x})}{r^{3}}+O\left(\frac{1}{r^{4}}\right).

Hence, in order for the property ${\bf v}-{\bf v}_{*}\approx\frac{1}{r}$ to hold, from the representation formula (129), it is necessary to have $\displaystyle\int_{\Omega}{\mathsf{G}_{\gamma}}(x-y)\mathrm{div}\,\mathcal{A}_{\gamma}(y)\,dy\approx\frac{1}{r}$ . By (222), this is true only if

\int_{{\mathbb{S}}^{2}}{\bf F}(\hat{x})\,d\sigma=0.

(131)

To verify the above for our system, we write $\mathcal{A}_{\gamma}$ explicitly as:

$\displaystyle\mathcal{A}_{\gamma}$	$\displaystyle=$	$\displaystyle\gamma_{1}\mathcal{A}_{1}+\gamma_{2}\mathcal{A}_{2}+\gamma_{3}\mathcal{A}_{3}+\gamma_{9}\mathcal{A}_{9}$	(132)
	$\displaystyle=$	$\displaystyle\gamma_{1}\Big{(}\mathsf{Q}_{}({\bf v}_{}\cdot\nabla\mathsf{Q})-({\bf v}_{}\cdot\nabla\mathsf{Q})\mathsf{Q}_{}\Big{)}+\gamma_{2}\Big{(}{\bf v}_{*}\cdot\nabla\mathsf{Q}\Big{)}$
		$\displaystyle+\frac{\gamma_{3}}{2}\left(\mathsf{Q}_{}({\bf v}_{}\cdot\nabla\mathsf{Q})+({\bf v}_{}\cdot\nabla\mathsf{Q})\mathsf{Q}_{}\right)+\gamma_{9}\Big{(}({\bf v}_{}\cdot\nabla\mathsf{Q})\cdot\mathsf{Q}_{}\Big{)}\mathsf{Q}_{*}.$

By considering only the dominating term $\displaystyle\frac{1}{r}$ in the expression (38) of $\mathsf{Q}$ , we have

$\displaystyle\mathsf{Q}-\mathsf{Q}_{*}$	$\displaystyle=$	$\displaystyle\left(-\frac{w\mathsf{Q}_{*}}{1+w}\right)\frac{1}{r}+O\left(\frac{1}{r^{2}}\right)$
$\displaystyle{\bf v}_{*}\cdot\nabla\mathsf{Q}$	$\displaystyle=$	$\displaystyle\left(-\frac{w\mathsf{Q}_{}}{1+w}\right)\left\langle{\bf v}_{},\nabla\frac{1}{r}\right\rangle+O\left(\frac{1}{r^{3}}\right),$
$\displaystyle\mathsf{Q}_{}({\bf v}_{}\cdot\nabla\mathsf{Q}),\,\,\,({\bf v}_{}\cdot\nabla\mathsf{Q})\mathsf{Q}_{}$	$\displaystyle=$	$\displaystyle\left(-\frac{w\mathsf{Q}_{}^{2}}{1+w}\right)\left\langle{\bf v}_{},\nabla\frac{1}{r}\right\rangle+O\left(\frac{1}{r^{3}}\right),$
$\displaystyle\big{(}({\bf v}_{}\cdot\nabla\mathsf{Q})\cdot\mathsf{Q}_{}\big{)}\mathsf{Q}_{*}$	$\displaystyle=$	$\displaystyle\left(-\frac{w\|\mathsf{Q}_{}\|^{2}\mathsf{Q}_{}}{1+w}\right)\left\langle{\bf v}_{*},\nabla\frac{1}{r}\right\rangle+O\left(\frac{1}{r^{3}}\right).$

Hence we can write for some constant $3\times 3$ matrix ${\bf M}_{*}=(m_{ij})$ that,

\mathrm{div}\,\mathcal{A}_{\gamma}=\text{div}\left(\left\langle{\bf v}_{*},\nabla\frac{1}{r}\right\rangle{\bf M}_{*}\right)+O\left(\frac{1}{r^{4}}\right)

(133)

Note that the term in $\mathcal{A}_{\gamma}$ multiplied by $\gamma_{1}$ completely vanishes.

We next claim that

\int_{\mathbb{S}^{2}}\left.\text{div}\left(\left\langle{\bf v}_{*},\nabla\frac{1}{r}\right\rangle{\bf M}_{*}\right)\right|_{r=1}\,d\sigma=0.

(134)

To see this, let $\displaystyle g(r)=\frac{1}{r}$ , and we compute for $i=1,2,3$ ,

\text{div}\left(\left\langle{\bf v}_{*},\nabla\frac{1}{r}\right\rangle{\bf M}_{*}\right)_{i}=\partial_{j}\big{(}{{\bf v}_{*}}_{k}\partial_{k}g(r)m_{ij}\big{)}={{\bf v}_{*}}_{k}m_{ij}\partial^{2}_{kj}g(r),

(135)

and

$\displaystyle\partial^{2}_{kj}g(r)$	$\displaystyle=$	$\displaystyle g^{\prime\prime}(r)\partial_{k}r\partial_{j}r+g^{\prime}(r)\partial^{2}_{kj}r$	(136)
	$\displaystyle=$	$\displaystyle g^{\prime\prime}(r)\frac{x_{k}x_{j}}{r^{2}}+\frac{g^{\prime}(r)}{r}\left(\delta_{kj}-\frac{x_{k}x_{j}}{r^{2}}\right)$
	$\displaystyle=$	$\displaystyle\frac{g^{\prime}(r)}{r}\delta_{kj}+\left(g^{\prime\prime}(r)-\frac{g^{\prime}(r)}{r}\right)\frac{x_{k}x_{j}}{r^{2}}$
	$\displaystyle=$	$\displaystyle-\frac{1}{r^{3}}\left(\delta_{kj}-3\frac{x_{k}x_{j}}{r^{2}}\right)$

where we have used the facts that $\displaystyle g^{\prime}(r)=-\frac{1}{r^{2}}$ and $\displaystyle g^{\prime\prime}(r)=\frac{2}{r^{3}}$ . As

\int_{\mathbb{S}^{2}}3\hat{x}\otimes\hat{x}\,d\sigma=4\pi\mathsf{I},

we have for all $k,j$ that

\int_{\mathbb{S}^{2}}\partial^{2}_{kj}g(r)\,d\sigma=0.

Thus (134) holds and so does (131).

5.3 Existence of solution of (111) in Schauder spaces

We recall the form of equation (111), its inhomogeneous term (116), and the Green’s function ${\mathsf{G}_{\gamma}}$ (126). Estimate (127) for ${\mathsf{G}_{\gamma}}$ plays an important role in our analysis.

We will show the existence and uniqueness of a solution by means of the Banach Fixed Point Theorem in a suitable weighted Schauder space $\mathcal{S}$ . For this purpose, we define

\mathcal{S}=\big{\{}{\bf u}:\Omega\rightarrow\mathbb{R}^{3},\,\,\|{\bf u}\|_{\mathcal{S}}<\infty\big{\}},

(137)

where

\|{\bf u}\|_{\mathcal{S}}=\sup_{|x|\in\Omega}\big{\{}|x||{\bf u}(x)-{\bf v}_{*}|,\,\,\,|x|^{2}|D{\bf u}(x)|,\,\,\,|x|^{3}|D^{2}{\bf u}(x)|,\,\,\,|x|^{3+\alpha}[D^{2}{\bf u}]_{\alpha}(x)\big{\}}.

(138)

In the above, $0<\alpha<1$ , and

[D^{2}{\bf u}]_{\alpha}(x)=\sup_{y\in\Omega,|x-y|\leq 1}\frac{|D^{2}{\bf u}(x)-D^{2}{\bf u}(y)|}{|x-y|^{\alpha}}.

(139)

Now, given a ${\bf u}\in\mathcal{S}$ , let ${\bf v}=\mathcal{T}({\bf u})$ be the solution of (111)–(114) with ${\bf f}_{\gamma}={\bf f}_{\gamma}({\bf u})$ . We will find a fixed point ${\bf v}$ of $\mathcal{T}$ : ${\bf v}=\mathcal{T}({\bf v})$ . To achieve this, we will show the following two properties of $\mathcal{T}$ :

(a).

$\mathcal{T}$ maps $\mathcal{S}$ into $\mathcal{S}$ , in particular, there is a $C>0$ such that for any ${\bf u}\in\mathcal{S}$ ,

\|\mathcal{T}({\bf u})\|_{\mathcal{S}}\leq C\Big{[}\|{\bf u}\|_{\mathcal{S}}+\|{\bf v}_{b}\|_{C^{2,\alpha}(\partial\Omega)}+|{\bf v}_{*}|\Big{]}.

(140)

(b).

For $|\gamma|\ll 1$ , there exists a $0<c<1$ such that for any ${\bf u}_{1},{\bf u}_{2}\in\mathcal{S}$ , we have

$||\mathcal{T}({\bf u}_{1})-\mathcal{T}({\bf u}_{2})||_{\mathcal{S}}\leq c||{\bf u}_{1}-{\bf u}_{2}||_{\mathcal{S}}.$ (141)

The proof can be obtained via the following steps.

(I) Well-posedness of $\mathcal{T}$ . General existence and uniqueness theory for ${\bf f}_{\gamma}\in L^{p}(\Omega)$ can be found in [32], also described in the encyclopedic reference [16, Chapter V]. But as ${\bf f}_{\gamma}$ has specific spatial decay property and we are looking for classical solutions, we find it convenient to follow the classical approach outlined in [33, Chapter 3] using the theory of single and double layer potentials, and Schauder estimates. This theory is also outlined in [34]. The recent survey [35] covers the existence and uniqueness of solutions for the Stokes equation in Schauder spaces. For the convenience of the reader, we now outline the approach that demonstrates existence.

Given ${\bf u}\in{\mathcal{S}}$ , we write ${\bf f}(x):={\bf f}_{\gamma}({\bf u}(x))$ as

$\displaystyle{{\bf f}}(x)$	$\displaystyle=$	$\displaystyle{\bf h}_{1}(x)+{\bf h}_{2}(x),$	(142)
$\displaystyle{\bf h}_{1}(x)$	$\displaystyle:=$	$\displaystyle\mathrm{div}\,\mathcal{A}_{\gamma}(x)=\gamma\frac{{\bf F}(\hat{x})}{r^{3}}+O\left(\frac{1}{r^{4}}\right)\,\,\,\left(\text{recall}\,\,\,\int_{\mathbb{S}^{2}}{\bf F}(\hat{x})\,d\sigma=0\right),$	(143)
$\displaystyle{\bf h}_{2}(x)$	$\displaystyle:=$	$\displaystyle\mathrm{div}\,\mathcal{C}_{\gamma}({\bf u},x)+\mathcal{D}_{\gamma}({\bf u},x).$	(144)

From (93) and (94), we have $\left|{\bf h}_{1}\right|\lesssim\frac{1}{r^{3}}$ and $\left|{\bf h}_{2}\right|\lesssim\frac{1}{r^{4}}$ for $r\gg 1$ . The solution ${\bf v}$ is found by the following steps.

(a).

We extend ${\bf f}$ smoothly from $\Omega$ to $\tilde{{\bf f}}$ defined on the whole space $\mathbb{R}^{3}$ . Given the decay property of $\tilde{{\bf f}}$ inherited from $\bf f$ , we can set

\tilde{{\bf v}}_{1}(x)=\int_{\mathbb{R}^{3}}{\mathsf{G}_{\gamma}}(x-y)\tilde{{\bf f}}(y)\,dy\,\,\,\text{and}\,\,\,\tilde{p}_{1}(x)=\int_{\mathbb{R}^{3}}{\mathsf{h}_{\gamma}}(x-y)\cdot\tilde{{\bf f}}(y)\,dy.

By (142)–(144) and (222), we have that $\tilde{{\bf v}}_{1}$ and $p_{1}$ solves the following equation

-\Delta\tilde{{\bf v}}_{1}+\nabla\tilde{p}_{1}=\tilde{{\bf f}}(x),\,\,\,\mathrm{div}\,\tilde{{\bf v}}_{1}=0,\,\,\,\text{with}\,\,\,|\tilde{{\bf v}}_{1}(x)|\lesssim\frac{1}{r}\,\,\,\text{and}\,\,\,|\tilde{p}_{1}|\lesssim\frac{1}{r^{2}},

on the whole of $\mathbb{R}^{3}$ .

(b).

We find ${\bf v}_{2}$ and $p_{2}$ on $\Omega$ that solve the following equation:

$\displaystyle-\Delta{\bf v}_{2}+\nabla p_{2}$	$\displaystyle=$	$\displaystyle 0,\,\,\,\text{in}\,\,\,\Omega$
$\displaystyle\mathrm{div}\,{\bf v}_{2}$	$\displaystyle=$	$\displaystyle 0,\,\,\,\text{in}\,\,\,\Omega$
$\displaystyle{\bf v}_{2}$	$\displaystyle=$	$\displaystyle{\bf v}_{b}-\tilde{{\bf v}}_{1}-{\bf v}_{*},\,\,\,\text{on}\,\,\,\partial\Omega,$
$\displaystyle{\bf v}_{2}$	$\displaystyle=$	$\displaystyle 0,\,\,\,\text{at}\,\,\,\|x\|=\infty.$

Following [34, Theorem 6.1, eqn. (6.3)], we can have that in the class ${\bf v}_{2}=o(r)$ , the solution exists, is unique, and it satisfies $|{\bf v}_{2}(x)|\lesssim\frac{1}{r}$ .

(c).

Define

${\bf v}=\tilde{{\bf v}}_{1}+{\bf v}_{2}+{\bf v}_{*},\,\,\,p=\tilde{p}_{1}+p_{2}.$

Then ${\bf v}$ and $p$ solve (111) on $\Omega$ .

The spatial decay property of ${\bf v}$ shows that it can have the following representation,

{\bf v}(x)={\bf v}_{*}+\int_{\Omega}{\mathsf{G}_{\gamma}}(x-y){\bf f}(y)\,dy\\ +\int_{\partial\Omega}\big{\langle}{\mathsf{G}_{\gamma}}(x-y),\mathsf{T}_{\gamma}({\bf v},p)(y)\nu_{y}\big{\rangle}\,d\sigma_{y}-\int_{\partial\Omega}\big{\langle}\mathsf{T}_{\gamma}({\mathsf{G}_{\gamma}},{\mathsf{h}_{\gamma}})(x-y)\nu_{y},{\bf v}(y)\big{\rangle}\,d\sigma_{y}.

(145)

Next we use (145) to estimates ${\bf v}$ . For estimates near the boundary, we can invoke [33, Chapter 3, Theorem 5] to deduce that

||{\bf v}||_{C^{2,\alpha}(\overline{\Omega}_{a+1})},\,\,||p||_{C^{1,\alpha}(\overline{\Omega}_{a+1})}\lesssim||{\bf f}||_{C^{0,\alpha}(\Omega)}+||{\bf v}_{b}||_{C^{2,\alpha}(\partial\Omega)}+|{\bf v}_{*}|

(146)

where $\Omega_{b}=\Omega\cap\{|x|\leq b\}$ . In particular, on the boundary, we have

||{\bf v}||_{C^{0}(\partial\Omega)},\,\,||D{\bf v}||_{C^{0}(\partial\Omega)},\,\,||p||_{C^{0}(\partial\Omega)}\lesssim||{\bf f}||_{C^{0,\alpha}(\Omega)}+||{\bf v}_{b}||_{C^{2,\alpha}(\partial\Omega)}+|{\bf v}_{*}|.

(147)

Furthermore, from the form of $\bf f$ (116), we have for some constant $C_{1},C_{2}$ that

||{\bf f}||_{C^{0,\alpha}(\Omega)}\leq\gamma\big{[}C_{1}\|{\bf u}\|_{\mathcal{S}}+C_{2}\big{]}.

(148)

Now we proceed with the following two proofs.

(II) Proof of (140). With boundary estimates given by (146), we will just concentrate here on interior weighted estimates, i.e., for $x$ such that $|x|>a+1$ . Again, we will utilize the representation (145).

For the boundary integrals in (145), we have

\displaystyle\left|\int_{\partial\Omega}\big{\langle}{\mathsf{G}_{\gamma}}(x-y),\mathsf{T}_{\gamma}({\bf v},p)(y)\nu_{y}\big{\rangle}\,d\sigma_{y}\right|

\displaystyle\lesssim\frac{||D{\bf v}||_{C^{0}(\partial\Omega)}+||p||_{C^{0}(\partial\Omega)}}{r}

and more generally, for $k=1,2,\ldots$ ,

	$\displaystyle\left\|D^{k}_{x}\int_{\partial\Omega}\big{\langle}{\mathsf{G}_{\gamma}}(x-y),\mathsf{T}_{\gamma}({\bf v},p)(y)\nu_{y}\big{\rangle}\,d\sigma_{y}\right\|$	$\displaystyle\lesssim\int_{\partial\Omega}\left\|\big{\langle}D^{k}_{x}{\mathsf{G}_{\gamma}}(x-y),\mathsf{T}_{\gamma}({\bf v},p)(y)\nu_{y}\big{\rangle}\right\|\,d\sigma_{y}$
		$\displaystyle\lesssim\frac{\|\|D{\bf v}\|\|_{C^{0}(\partial\Omega)}+\|\|p\|\|_{C^{0}(\partial\Omega)}}{r^{1+k}}.$

Furthermore, for $h\in\mathbb{R}^{n}$ with $|h|\leq 1$ ,

			$\displaystyle\frac{\|x\|^{3+\alpha}}{\|h\|^{\alpha}}\left\|\int_{\partial\Omega}\big{\langle}D_{x}^{2}{\mathsf{G}_{\gamma}}(x+h-y)-D_{x}^{2}{\mathsf{G}_{\gamma}}(x-y),\mathsf{T}_{\gamma}({\bf v},p)(y)\nu_{y}\big{\rangle}\,d\sigma_{y}\right\|$
		$\displaystyle\lesssim$	$\displaystyle\big{(}\\|{\bf v}\\|_{C^{1}(\partial\Omega)}+\\|p\\|_{C^{0}(\partial\Omega)}\big{)}\frac{\|x\|^{3+\alpha}}{\|h\|^{\alpha}}\left\|\int_{\partial\Omega}\frac{\big{\|}\|x+h-y\|^{3}-\|x-y\|^{3}\big{\|}}{\|x+h-y\|^{3}\|x-y\|^{3}}\,d\sigma_{y}\right\|$
		$\displaystyle\lesssim$	$\displaystyle\big{(}\\|{\bf v}\\|_{C^{1}(\partial\Omega)}+\\|p\\|_{C^{0}(\partial\Omega)}\big{)}\left(\frac{\|x\|^{3+\alpha}}{\|h\|^{\alpha}}\right)\frac{\|x\|^{2}\|h\|}{\|x\|^{6}}$
		$\displaystyle\lesssim$	$\displaystyle\big{(}\\|{\bf v}\\|_{C^{1}(\partial\Omega)}+\\|p\\|_{C^{0}(\partial\Omega)}\big{)}\left(\frac{\|h\|^{1-\alpha}}{\|x\|^{1-\alpha}}\right)$
		$\displaystyle\lesssim$	$\displaystyle\big{(}\\|{\bf v}\\|_{C^{1}(\partial\Omega)}+\\|p\\|_{C^{0}(\partial\Omega)}\big{)}.$

Hence,

$\displaystyle\left\\|\int_{\partial\Omega}\big{\langle}{\mathsf{G}_{\gamma}}(x-y),\mathsf{T}_{\gamma}({\bf v},p)(y)\nu_{y}\big{\rangle}\,d\sigma_{y}\right\\|_{\cal S}$	$\displaystyle\lesssim$	$\displaystyle\\|{\bf v}\\|_{C^{1}(\partial\Omega)}+\\|p\\|_{C^{0}(\partial\Omega)}$
	$\displaystyle\lesssim$	$\displaystyle\|\|{\bf f}\|\|_{C^{0,\alpha}(\Omega)}+\|\|{\bf v}_{b}\|\|_{C^{2,\alpha}(\partial\Omega)}$
	$\displaystyle\lesssim$	$\displaystyle\|\gamma\|\Big{[}\\|{\bf u}\\|_{\mathcal{S}}+1\Big{]}+\|\|{\bf v}_{b}\|\|_{C^{2,\alpha}(\partial\Omega)}+\|{\bf v}_{*}\|.$

Using the same technique, we can similarly have,

$\displaystyle\left\\|\|x\|\left\|\int_{\partial\Omega}\big{\langle}\mathsf{T}_{\gamma}({\mathsf{G}_{\gamma}},{\mathsf{h}_{\gamma}})(x-y)\nu_{y},{\bf v}(y)\big{\rangle}d\sigma_{y}\right\|\right\\|_{\cal S}$	$\displaystyle\lesssim$	$\displaystyle\\|{\bf v}\\|_{C^{1}(\partial\Omega)}+\\|p\\|_{C^{0}(\partial\Omega)}$
	$\displaystyle\lesssim$	$\displaystyle\|\|{\bf f}\|\|_{C^{0,\alpha}(\Omega)}+\|\|{\bf v}_{b}\|\|_{C^{2,\alpha}(\partial\Omega)}$
	$\displaystyle\lesssim$	$\displaystyle\|\gamma\|\Big{[}\\|{\bf u}\\|_{\mathcal{S}}+1\Big{]}+\|\|{\bf v}_{b}\|\|_{C^{2,\alpha}(\partial\Omega)}+\|{\bf v}_{*}\|.$

For the bulk integrals in (145), the proof is very similar to the classical estimates for Newtonian potential – see [36, Chapter 4]. We will produce the proof for our weighted space $\cal S$ in Appendix B.

(III) Proof of (141). With (140), estimate (141) follows easily, due to the linearity of the equation. Given ${\bf u}_{1},{\bf u}_{2}\in\mathcal{S}$ , we have solutions ${\bf v}_{1}=\mathcal{T}({\bf u}_{1})$ and ${\bf v}_{2}=\mathcal{T}({\bf u}_{2})$ . Upon subtracting the corresponding equations, we deduce that

\displaystyle{\mathcal{L}}_{\gamma}({\bf v}_{1}-{\bf v}_{2})+\nabla(p_{1}-p_{2})={\bf f}_{\gamma}({\bf u}_{1})-{\bf f}_{\gamma}({\bf u}_{2}),

(149)

with ${\bf v}_{1}-{\bf v}_{2}=\bf 0$ on $\partial\Omega$ and at $|x|=\infty$ . Hence

||{\bf v}_{1}-{\bf v}_{2}||_{\mathcal{S}}\lesssim||{\bf f}_{\gamma}({\bf u}_{1})-{\bf f}_{\gamma}({\bf u}_{2})||_{C^{0,\alpha}(\Omega)}\lesssim|\gamma|\|{\bf u}_{1}-{\bf u}_{2}\|_{\mathcal{S}},

where we have used (117). Now choosing $|\gamma|\ll 1$ gives the result.

6 Properties of anisotropic Stokes flows

Here we make more precise the far-field behavior of the solution ${\bf v}$ . As an application of our analysis, we will analyze the symmetry property of the solution and give a decomposition formula for the Stokes drag. They will be validated and illustrated by numerical simulations.

We start from representation (145) for the solution ${\bf v}$ :

$\displaystyle{\bf v}(x)$	$\displaystyle=$	$\displaystyle\int_{\Omega}{\mathsf{G}_{\gamma}}(x-y){\bf f}_{\gamma}({\bf v}(y))\,dy+\int_{\partial\Omega}\big{\langle}{\mathsf{G}_{\gamma}}(x-y),\mathsf{T}_{\gamma}({\bf v},p)(y)\nu_{y}\big{\rangle}\,d\sigma_{y}$	(150)
		$\displaystyle-\int_{\partial\Omega}\big{\langle}\mathsf{T}_{\gamma}({\mathsf{G}_{\gamma}},{\mathsf{h}_{\gamma}})(x-y)\nu_{y},{\bf v}(y)\big{\rangle}\,d\sigma_{y}+{\bf v}_{*}$
	$\displaystyle=:$	$\displaystyle\mathsf{S}_{1}(x)+\mathsf{S}_{2}(x)+\mathsf{S}_{3}(x)+{\bf v}_{*}$

where

$\displaystyle\mathsf{S}_{1}(x)$	$\displaystyle:=$	$\displaystyle\int_{\Omega}{\mathsf{G}_{\gamma}}(x-y){\bf f}_{\gamma}({\bf v}(y))\,dy$
	$\displaystyle=$	$\displaystyle\int_{\Omega}{\mathsf{G}_{\gamma}}(x-y)\mathrm{div}\,\mathcal{A}_{\gamma}(y)\,dy+\int_{\Omega}{\mathsf{G}_{\gamma}}(x-y)\big{(}\mathrm{div}\,\mathcal{C}_{\gamma}({\bf v})+\mathcal{D}_{\gamma}({\bf v})\big{)}(y)\,dy$
	$\displaystyle=$	$\displaystyle\int_{\Omega}{\mathsf{G}_{\gamma}}(x-y)\mathrm{div}\,\mathcal{A}_{\gamma}(y)\,dy+{\mathsf{G}_{\gamma}}(x)\left(\int_{\Omega}\big{(}\mathrm{div}\,\mathcal{C}_{\gamma}({\bf v})+\mathcal{D}_{\gamma}({\bf v})\big{)}(y)\,dy\right)$
		$\displaystyle+\int_{\Omega}({\mathsf{G}_{\gamma}}(x-y)-{\mathsf{G}_{\gamma}}(x))\big{(}\mathrm{div}\,\mathcal{C}_{\gamma}({\bf v})+\mathcal{D}_{\gamma}({\bf v})\big{)}(y)\,dy,$
$\displaystyle\mathsf{S}_{2}(x)$	$\displaystyle:=$	$\displaystyle\int_{\partial\Omega}\big{\langle}{\mathsf{G}_{\gamma}}(x-y),\mathsf{T}_{\gamma}({\bf v},p)(y)\nu_{y}\big{\rangle}\,d\sigma_{y}$
	$\displaystyle=$	$\displaystyle\left\langle{\mathsf{G}_{\gamma}}(x),\int_{\partial\Omega}\mathsf{T}_{\gamma}({\bf v},p)(y)\nu_{y}\,d\sigma_{y}\right\rangle+\int_{\partial\Omega}\Big{\langle}{\mathsf{G}_{\gamma}}(x-y)-{\mathsf{G}_{\gamma}}(x),\mathsf{T}_{\gamma}({\bf v},p)(y)\nu_{y}\Big{\rangle}\,d\sigma_{y},$
$\displaystyle\mathsf{S}_{3}(x)$	$\displaystyle:=$	$\displaystyle-\int_{\partial\Omega}\big{\langle}\mathsf{T}_{\gamma}({\mathsf{G}_{\gamma}},{\mathsf{h}_{\gamma}})(x-y)\nu_{y},{\bf v}(y)\big{\rangle}\,d\sigma_{y}.$

We re-arrange the above terms in the following way,

$\displaystyle{\bf v}-{\bf v}_{*}$	$\displaystyle=$	$\displaystyle\left\langle{\mathsf{G}_{\gamma}}(x),\int_{\partial\Omega}\mathsf{T}_{\gamma}({\bf v},p)(y)\nu_{y}\,d\sigma_{y}\right\rangle$	(155)
		$\displaystyle+\int_{\Omega}{\mathsf{G}_{\gamma}}(x-y)\mathrm{div}\,\mathcal{A}_{\gamma}(y)\,dy+{\mathsf{G}_{\gamma}}(x)\left(\int_{\Omega}\big{(}\mathrm{div}\,\mathcal{C}_{\gamma}({\bf v})+\mathcal{D}_{\gamma}({\bf v})\big{)}(y)\,dy\right)$
		$\displaystyle+\int_{\partial\Omega}\Big{\langle}{\mathsf{G}_{\gamma}}(x-y)-{\mathsf{G}_{\gamma}}(x),\mathsf{T}_{\gamma}({\bf v},p)(y)\nu_{y}\Big{\rangle}\,d\sigma_{y}$
		$\displaystyle-\int_{\partial\Omega}\big{\langle}\mathsf{T}_{\gamma}({\mathsf{G}_{\gamma}},{\mathsf{h}_{\gamma}})(x-y)\nu_{y},{\bf v}(y)\big{\rangle}\,d\sigma_{y}$
		$\displaystyle+\int_{\Omega}({\mathsf{G}_{\gamma}}(x-y)-{\mathsf{G}_{\gamma}}(x))\big{(}\mathrm{div}\,\mathcal{C}_{\gamma}({\bf v})+\mathcal{D}_{\gamma}({\bf v})\big{)}(y)\,dy.$

Note that for $|x|\gg 1$ , we have

(155)	$\displaystyle=$	$\displaystyle\frac{O(1)}{\|x\|},$
(155)	$\displaystyle=$	$\displaystyle\frac{O(\gamma)}{\|x\|},$
(155) and (155)	$\displaystyle=$	$\displaystyle\frac{O(1)}{\|x\|^{2}},$
(155)	$\displaystyle=$	$\displaystyle\frac{O(\gamma)\log\|x\|}{\|x\|^{2}}\quad\text{(by \eqref{logx2})}$

so that up to $O(\gamma)$ , (155) and (155) are the dominant terms in the expression for the anisotropic flow.

We next describe asymptotically (155) and (155) for $\gamma\ll 1$ .

6.1 Precise asymptotics: deviation from isotropic Stokes flow

The purpose of this section is to reveal more clearly the difference between ${\bf v}$ and the classical Stokes flow ${\bf v}_{0}$ which is set to satisfy

$\displaystyle-\Delta{\bf v}_{0}+\nabla p_{0}$	$\displaystyle=$	$\displaystyle 0,\quad\text{for}\,\,\,\|x\|>a,$	(156)
$\displaystyle\mathrm{div}\,({\bf v}_{0})$	$\displaystyle=$	$\displaystyle 0,\quad\text{for}\,\,\,\|x\|>a,$
$\displaystyle{\bf v}_{0}$	$\displaystyle=$	$\displaystyle 0,\,\,\,\text{on $\|x\|=a$,}$
$\displaystyle{\bf v}_{0}$	$\displaystyle=$	$\displaystyle{\bf v}_{*},\,\,\,\text{at $\|x\|=\infty$.}$

Define $\varphi_{\gamma}:={\bf v}-{\bf v}_{0}$ . Then we have

	$\displaystyle-\Delta\varphi_{\gamma}+\nabla(p-p_{0})$	$\displaystyle=$	$\displaystyle-\mathcal{M}_{\gamma}:D^{2}{\bf v}+{\bf f}_{\gamma}({\bf v})$
		$\displaystyle=$	$\displaystyle-\mathcal{M}_{\gamma}:D^{2}{\bf v}_{0}+{\bf f}_{\gamma}({\bf v}_{0})+\big{(}{\bf f}_{\gamma}({\bf v})-{\bf f}_{\gamma}({\bf v}_{0})-\mathcal{M}_{\gamma}:(D^{2}{\bf v}-D^{2}{\bf v}_{0})\big{)}.$

Note that $|\mathcal{M}_{\gamma}|,\,\,|{\bf f}_{\gamma}({\bf v}_{0})|\lesssim|\gamma|$ and $|{\bf f}_{\gamma}({\bf v})-{\bf f}_{\gamma}({\bf v}_{0})|=|{\bf f}_{\gamma}({\bf v}-{\bf v}_{0})|\lesssim|\gamma|\|{\bf v}-{\bf v}_{0}\|_{\mathcal{S}}$ . Now let $\overline{\varphi}_{\gamma}$ solve

$\displaystyle-\Delta\overline{\varphi}_{\gamma}+\nabla(\overline{p}_{\gamma})$	$\displaystyle=$	$\displaystyle-\mathcal{M}_{\gamma}:D^{2}{\bf v}_{0}+{\bf f}_{\gamma}({\bf v}_{0}),$	(157)
$\displaystyle\text{div}(\overline{\varphi}_{\gamma})$	$\displaystyle=$	$\displaystyle 0,\quad\text{for $\|x\|>a$,}$	(158)
$\displaystyle\overline{\varphi}_{\gamma}$	$\displaystyle=$	$\displaystyle 0,\quad\text{on $\|x\|=0$ and at $\|x\|=\infty$.}$	(159)

Then the same approach in deriving estimates for ${\bf v}$ gives

\|\overline{\varphi}_{\gamma}\|_{\mathcal{S}}\lesssim O(\gamma),\quad\text{and}\quad\|\overline{\varphi}_{\gamma}-\varphi_{\gamma}\|_{\mathcal{S}}\lesssim O(\gamma^{2}).

(160)

Hence we have

{\bf v}={\bf v}_{0}+\overline{\varphi}_{\gamma}+O(\gamma^{2}).

(161)

Finally, using the Green’s function $\mathsf{E}$ (43) for the classical Stokes flow, we have the following representation of $\overline{\varphi}_{\gamma}$ ,

$\displaystyle\overline{\varphi}_{\gamma}(x)$	$\displaystyle=$	$\displaystyle\int_{\Omega}\mathsf{E}(x-y)\big{(}-\mathcal{M}_{\gamma}:D^{2}{\bf v}_{0}+{\bf f}_{\gamma}({\bf v}_{0}(y))\big{)}\,dy+\int_{\partial\Omega}\big{\langle}\mathsf{E}(x-y),\mathsf{T}(\overline{\varphi}_{\gamma},\overline{p}_{\gamma})(y)\nu_{y}\big{\rangle}\,d\sigma_{y}.$	(162)
	$\displaystyle=$	$\displaystyle\int_{\Omega}\mathsf{E}(x-y)\big{(}-\mathcal{M}_{\gamma}:D^{2}{\bf v}_{0}+{\bf f}_{\gamma}({\bf v}_{0}(y))\big{)}\,dy+\mathsf{E}(x)\int_{\partial\Omega}\mathsf{T}(\overline{\varphi}_{\gamma},\overline{p}_{\gamma})(y)\nu_{y}\,d\sigma_{y}$
		$\displaystyle+\int_{\partial\Omega}\big{\langle}\mathsf{E}(x-y)-\mathsf{E}(x),\mathsf{T}(\overline{\varphi}_{\gamma},\overline{p}_{\gamma})(y)\nu_{y}\big{\rangle}\,d\sigma_{y}$

where $\mathsf{T}$ is the Stokes stress tensor (48). Note that the last term in the above decays as $\frac{1}{|x|^{2}}$ . Hence we have,

$\displaystyle\overline{\varphi}_{\gamma}(x)$	$\displaystyle=$	$\displaystyle\mathcal{I}_{\gamma}(x)+\mathsf{E}(x)\mathcal{J}_{\gamma}+O\left(\frac{1}{\|x\|^{2}}\right),$	(163)
$\displaystyle\text{where}\quad\mathcal{I}_{\gamma}(x)$	$\displaystyle:=$	$\displaystyle\int_{\Omega}\mathsf{E}(x-y)\big{(}-\mathcal{M}_{\gamma}:D^{2}{\bf v}_{0}+{\bf f}_{\gamma}({\bf v}_{0}(y))\big{)}\,dy,$	(164)
$\displaystyle\mathcal{J}_{\gamma}$	$\displaystyle:=$	$\displaystyle\int_{\partial\Omega}\mathsf{T}(\overline{\varphi}_{\gamma},\overline{p}_{\gamma})(y)\nu_{y}d\sigma_{y}.$	(165)

The rest of this section will describe more explicitly the bulk and boundary integrals $\mathcal{I}_{\gamma}$ and $\mathcal{J}_{\gamma}$ .

6.1.1 Analysis of $\mathcal{I}_{\gamma}$

For this purpose, recalling from (67) that ${\bf v}_{0}=\mathsf{E}_{S}{\bf v}_{*}$ , we then have

$\displaystyle\mathcal{I}_{\gamma}(x)$	$\displaystyle=$	$\displaystyle\int_{\Omega}\mathsf{E}(x-y)\big{[}-\mathcal{M}_{\gamma}:D^{2}{\bf v}_{0}+{\bf f}_{\gamma}({\bf v}_{0}(y))\big{]}\,dy$
	$\displaystyle=$	$\displaystyle\int_{\Omega}\mathsf{E}(x-y)\Big{[}6\pi a\mathcal{M}_{\gamma}:D^{2}\mathsf{E}(y){\bf v}_{}-\frac{a^{3}}{4}\mathcal{M}_{\gamma}:D^{2}\mathsf{F}(y){\bf v}_{}$
		$\displaystyle\quad\quad\quad\quad\quad+\mathrm{div}\,\mathcal{A}_{\gamma}(y)+\mathrm{div}\,\mathcal{C}_{\gamma}(y)+\mathcal{D}_{\gamma}(y)\Big{]}\,dy$
	$\displaystyle=$	$\displaystyle\int_{\Omega}\mathsf{E}(x-y)\Big{[}6\pi a\mathcal{M}_{\gamma}:D^{2}\mathsf{E}(y){\bf v}_{*}+\mathrm{div}\,\mathcal{A}_{\gamma}(y)\Big{]}\,dy$
		$\displaystyle+\mathsf{E}(x)\int_{\Omega}\left[-\frac{a^{3}}{4}\mathcal{M}_{\gamma}:D^{2}\mathsf{F}(y){\bf v}_{*}+\mathrm{div}\,\mathcal{C}_{\gamma}(y)+\mathcal{D}_{\gamma}(y)\right]\,dy$
		$\displaystyle+\int_{\Omega}(\mathsf{E}(x-y)-\mathsf{E}(x))\left[-\frac{a^{3}}{4}\mathcal{M}_{\gamma}:D^{2}\mathsf{F}(y){\bf v}_{*}+\mathrm{div}\,\mathcal{C}_{\gamma}(y)+\mathcal{D}_{\gamma}(y)\right]\,dy$
	$\displaystyle=:$	$\displaystyle I_{1}+I_{2}+I_{3}.$	(167)

Note that $I_{2}$ decays as $\mathsf{E}(x)\sim\frac{1}{|x|}$ because the integrand is integrable. By (205), the integral $I_{3}$ decays as $\frac{\log|x|}{|x|^{2}}$ .

For $I_{1}$ , we will show that $D^{2}\mathsf{E}(x)$ satisfies the mean zero condition. (Such a condition is already verified for $\mathrm{div}\,\mathcal{A}_{\gamma}$ in Section 5.2.) By (222), we can then conclude that $I_{1}(x)$ decays as $\frac{1}{|x|}$ . To this end, by (232), we have

$\displaystyle\int_{\mathbb{S}^{2}}\partial_{kl}\mathsf{E}_{ij}(x)\,d\sigma$	$\displaystyle=$	$\displaystyle 4\pi\left(-\delta_{ij}\delta_{kl}+\delta_{ik}\delta_{jl}+\delta_{jk}\delta_{il}\right)-4\pi\left(-\delta_{ij}\delta_{kl}+\delta_{ik}\delta_{jl}+\delta_{jk}\delta_{il}\right)$	(168)
		$\displaystyle-4\pi\left(\delta_{il}\delta_{jk}+\delta_{jl}\delta_{ik}+\delta_{kl}\delta_{ij}\right)+15\int_{\mathbb{S}^{2}}\hat{x}_{i}\hat{x}_{j}\hat{x}_{k}\hat{x}_{l}\,d\sigma$
	$\displaystyle=$	$\displaystyle-4\pi\left(\delta_{il}\delta_{jk}+\delta_{jl}\delta_{ik}+\delta_{kl}\delta_{ij}\right)+15\int_{\mathbb{S}^{2}}\hat{x}_{i}\hat{x}_{j}\hat{x}_{k}\hat{x}_{l}\,d\sigma.$

Using spherical coordinates and considering symmetry, we check

(a).

$i\neq j,k\neq l$ , $(k,l)=(i,j)$ or $(j,i)$

$\displaystyle\int_{\mathbb{S}^{2}}\partial_{kl}\mathsf{E}_{ij}(x)\,d\sigma$	$\displaystyle=$	$\displaystyle-4\pi+15\int_{\mathbb{S}^{2}}\hat{x}_{1}^{2}\hat{x}_{2}^{2}\,d\sigma$
	$\displaystyle=$	$\displaystyle-4\pi+15\int_{0}^{\pi}\int_{0}^{2\pi}(\sin\phi\cos\theta)^{2}(\sin\phi\sin\theta)^{2}\sin\phi\,d\phi\,d\theta$
	$\displaystyle=$	$\displaystyle-4\pi+15\int_{0}^{\pi}\sin^{5}\phi\,d\phi\int_{0}^{2\pi}\cos^{2}\theta\sin^{2}\theta\,d\theta=0$

(b).

$i\neq j,k\neq l$ , $(k,l)\neq(i,j),(j,i)$

$\displaystyle\int_{\mathbb{S}^{2}}\partial_{kl}\mathsf{E}_{ij}(x)\,d\sigma=0$

(c).

$i\neq j,k=l=i$ or $k=l=j$

\displaystyle\int_{\mathbb{S}^{2}}\partial_{kl}\mathsf{E}_{ij}(x)\,d\sigma=+15\int_{\mathbb{S}^{2}}\hat{x}_{1}^{3}\hat{x}_{2}\,d\sigma=0.

(d).

$i\neq j,k=l,\neq i,j$

$\displaystyle\int_{\mathbb{S}^{2}}\partial_{kl}\mathsf{E}_{ij}(x)\,d\sigma=0$
(e).

$i=j,k\neq l$

$\displaystyle\int_{\mathbb{S}^{2}}\partial_{kl}\mathsf{E}_{ij}(x)\,d\sigma=0$

(f).

$i=j,k=l,\neq i,j$

\displaystyle\int_{\mathbb{S}^{2}}\partial_{kl}\mathsf{E}_{ij}(x)\,d\sigma=-4\pi+15\int_{\mathbb{S}^{2}}\hat{x}_{1}^{2}\hat{x}_{2}^{2}\,d\sigma=0,\quad\text{as in (a)}.

(g).

$i=j=k=l$

\displaystyle\int_{\mathbb{S}^{2}}\partial_{kl}\mathsf{E}_{ij}(x)\,d\sigma=-12\pi+15\int_{\mathbb{S}^{2}}\hat{x}_{3}^{4}\,d\sigma=-12\pi+15\int_{0}^{\pi}\int_{0}^{2\pi}\cos^{4}\phi\sin\phi\,d\theta\,d\phi=0.

Hence we can invoke (222) to conclude that

$\displaystyle\mathcal{I}_{\gamma}(x)$	$\displaystyle=$	$\displaystyle\frac{1}{\|x\|}\int_{\mathbb{S}^{2}}\mathsf{H}(\hat{x},\omega)\Big{[}6\pi a\mathcal{M}_{\gamma}:D^{2}\mathsf{E}(\omega){\bf v}_{*}+\mathrm{div}\,\mathcal{A}_{\gamma}(\omega)\Big{]}\,d\sigma_{\omega}$
		$\displaystyle+\mathsf{E}(x)\int_{\Omega}\left[-\frac{a^{3}}{4}\mathcal{M}_{\gamma}:D^{2}\mathsf{F}(y){\bf v}_{*}+\mathrm{div}\,\mathcal{C}_{\gamma}(y)+\mathcal{D}_{\gamma}(y)\right]\,dy$
		$\displaystyle+O\left(\frac{\log\|x\|}{\|x\|^{2}}\right)$

6.1.2 Analysis of $\mathcal{J}_{\gamma}$

To obtain an explicit formula for the boundary stress $\mathcal{J}_{\gamma}$ associated with $\overline{\varphi}_{\gamma}$ , we multiply (157) by a test function given by $\tilde{{\bf v}}_{0}-\tilde{{\bf v}}_{*}$ where $\tilde{{\bf v}}_{0}$ solves the following equation:

$\displaystyle-\Delta\tilde{{\bf v}}_{0}+\nabla\tilde{p}_{0}$	$\displaystyle=$	$\displaystyle 0,\quad\text{for}\,\,\,\|x\|>a,$
$\displaystyle\mathrm{div}\,(\tilde{{\bf v}}_{0})$	$\displaystyle=$	$\displaystyle 0,\quad\text{for}\,\,\,\|x\|>a,$
$\displaystyle\tilde{{\bf v}}_{0}$	$\displaystyle=$	$\displaystyle 0,\,\,\,\text{on $\|x\|=a$,}$
$\displaystyle\tilde{{\bf v}}_{0}$	$\displaystyle=$	$\displaystyle\tilde{{\bf v}}_{*},\,\,\,\text{at $\|x\|=\infty$.}$

with an arbitrary $\tilde{{\bf v}}_{*}$ . Using integration by parts – see [33, p.53 (10), (11)], we obtain

-\int_{\partial\Omega}\Big{\langle}\mathsf{T}(\overline{\varphi}_{\gamma},\overline{p}_{\gamma})\nu_{y},(\tilde{{\bf v}}_{0}-\tilde{{\bf v}}_{*})\Big{\rangle}\,d\sigma_{y}+\int_{\partial\Omega}\Big{\langle}\mathsf{T}(\tilde{{\bf v}}_{0},\tilde{p}_{0})\nu_{y},\overline{\varphi}_{\gamma}\Big{\rangle}\,d\sigma_{y}\\ =\int_{\Omega}\Big{\langle}-\mathcal{M}_{\gamma}:D^{2}{\bf v}_{0}+{\bf f}_{\gamma}({\bf v}_{0}(y)),\tilde{{\bf v}}_{0}-\tilde{{\bf v}}_{*}\Big{\rangle}\,dy.

The above is justified by the decay estimates for $\overline{\varphi}_{\gamma}$ and $\tilde{{\bf v}}_{0}-\tilde{{\bf v}}_{*}$ :

	$\displaystyle\left\|\int_{\{\|x\|=R\}}\Big{\langle}\mathsf{T}(\overline{\varphi}_{\gamma},\overline{p}_{\gamma})\nu_{y},(\tilde{{\bf v}}_{0}-\tilde{{\bf v}}_{*})\Big{\rangle}\,d\sigma_{y}\right\|\lesssim\frac{1}{R^{2}}\frac{1}{R}R^{2}=\frac{1}{R}$
	$\displaystyle\left\|\int_{\{\|x\|=R\}}\Big{\langle}\mathsf{T}(\tilde{{\bf v}}_{0},\tilde{p}_{0})\nu_{y},\overline{\varphi}_{\gamma}\Big{\rangle}\,d\sigma_{y}\right\|\lesssim\frac{1}{R^{2}}\frac{1}{R}R^{2}=\frac{1}{R}$
	$\displaystyle\int_{\{\|x\|>R\}}\Big{\|}-\mathcal{M}_{\gamma}:D^{2}{\bf v}_{0}+{\bf f}({\bf v}_{0}(y))\Big{\|}\Big{\|}\tilde{{\bf v}}_{0}-{\bf v}_{*}\Big{\|}\,dy\lesssim\int_{R}^{\infty}\frac{1}{r^{3}}\frac{1}{r}r^{2}\,dr\lesssim\frac{1}{R}$

which all vanish as $R\longrightarrow 0$ . Hence we have

			$\displaystyle\int_{\partial\Omega}\Big{\langle}\mathsf{T}(\overline{\varphi}_{\gamma},\overline{p}_{\gamma})\nu_{y},(\tilde{{\bf v}}_{0}-\tilde{{\bf v}}_{*})\Big{\rangle}\,d\sigma_{y}$
		$\displaystyle=$	$\displaystyle\int_{\partial\Omega}\Big{\langle}\mathsf{T}(\tilde{{\bf v}}_{0},\tilde{p}_{0})\nu_{y},\overline{\varphi}_{\gamma}\Big{\rangle}\,d\sigma_{y}-\int_{\Omega}\Big{\langle}-\mathcal{M}_{\gamma}:D^{2}{\bf v}_{0}+{\bf f}({\bf v}_{0}(y)),\tilde{{\bf v}}_{0}-\tilde{{\bf v}}_{*}\Big{\rangle}\,dy.$

As $\tilde{{\bf v}}_{0}=\overline{\varphi}_{\gamma}=0$ on $\partial\Omega$ , the following holds for any $\tilde{{\bf v}}_{*}$ ,

	$\displaystyle\int_{\partial\Omega}\Big{\langle}\mathsf{T}(\overline{\varphi}_{\gamma},\overline{p}_{\gamma})\nu_{y},\tilde{{\bf v}}_{*}\Big{\rangle}\,d\sigma_{y}$	$\displaystyle=$	$\displaystyle\int_{\Omega}\Big{\langle}-\mathcal{M}_{\gamma}:D^{2}{\bf v}_{0}+{\bf f}_{\gamma}({\bf v}_{0}(y)),\tilde{{\bf v}}_{0}-\tilde{{\bf v}}_{*}\Big{\rangle}\,dy$
		$\displaystyle=$	$\displaystyle\int_{\Omega}\Big{\langle}-\mathcal{M}_{\gamma}:D^{2}\mathsf{E}_{S}(y){\bf v}_{}+{\bf f}_{\gamma}(\mathsf{E}_{S}(y){\bf v}_{}),(\mathsf{E}_{S}-\mathsf{I})\tilde{{\bf v}}_{*}\Big{\rangle}\,dy.$

Hence we have

$\displaystyle\mathcal{J}_{\gamma}$	$\displaystyle=$	$\displaystyle\int_{\partial\Omega}\mathsf{T}(\overline{\varphi}_{\gamma},\overline{p}_{\gamma})\nu_{y}\,d\sigma_{y}$	(171)
	$\displaystyle=$	$\displaystyle\int_{\Omega}(\mathsf{E}_{S}-\mathsf{I})\Big{[}-\mathcal{M}_{\gamma}:D^{2}\mathsf{E}_{S}(y){\bf v}_{}+{\bf f}_{\gamma}\big{(}\mathsf{E}_{S}(y){\bf v}_{}\big{)}\Big{]}\,dy$
	$\displaystyle=$	$\displaystyle\int_{\Omega}(\mathsf{E}_{S}-\mathsf{I})\Big{[}-\mathcal{M}_{\gamma}:D^{2}\mathsf{E}_{S}(y){\bf v}_{*}$
		$\displaystyle\hskip 60.0pt+\mathrm{div}\,\mathcal{A}_{\gamma}(y)+\mathrm{div}\,\mathcal{C}_{\gamma}(\mathsf{E}_{S}(y){\bf v}_{},y)+\mathcal{D}_{\gamma}(\mathsf{E}_{S}(y){\bf v}_{},y)\Big{]}\,dy.$

6.2 Symmetry properties of solution and drag force

Here we investigate the symmetries of ${\bf v}$ with respect to our far-field data $\mathsf{Q}_{*}$ and ${\bf v}_{*}$ . Note that the interacting potential $R$ in (12) is chosen to be frame indifferent. Hence our solution is naturally invariant with respect to orthogonal transformations. For the reader’s convenience, we outline the derivation.

Let ${\mathbb{B}}=(b_{ij})$ and ${\mathbb{A}}=(a_{ij})={\mathbb{B}}^{T}(={\mathbb{B}}^{-1})$ be an orthogonal matrix and its transpose (and inverse). Upon introducing

\displaystyle y={\mathbb{B}}x,\,\,\,{\bf v}(x,t)={\mathbb{A}}\widetilde{{\bf v}}(y,t),\,\,\,p(x,t)=\widetilde{p}(y,t),\,\,\,\partial_{i}=\partial_{x_{i}},\,\,\,\widetilde{\partial}_{i}=\partial_{y_{i}}\,\,\,(\text{so that}\,\,\,\partial_{i}=a_{ij}\widetilde{\partial}_{j}),

we have,

	$\displaystyle{\bf v}\cdot\nabla=\widetilde{{\bf v}}\cdot\widetilde{\nabla},\,\,\,{\bf v}\cdot\nabla{\bf v}={\mathbb{A}}\left(\widetilde{{\bf v}}\cdot\widetilde{\nabla}\widetilde{{\bf v}}\right),\,\,\,\mathrm{div}\,{\bf v}=\widetilde{\mathrm{div}\,}\widetilde{{\bf v}},\,\,\,\triangle{\bf v}={\mathbb{A}}\,\widetilde{\triangle}\widetilde{{\bf v}},$
	$\displaystyle\nabla p={\mathbb{A}}\widetilde{\nabla}\widetilde{p},\,\,\,{\bf v}_{t}={\mathbb{A}}\widetilde{{\bf v}}_{t}.$

The above shows that the Navier-Stokes equation

{\bf v}_{t}+{\bf v}\cdot\nabla{\bf v}+\nabla p=\triangle{\bf v},\,\,\,\mathrm{div}\,{\bf v}=0,

is equivalent to

\widetilde{{\bf v}}_{t}+\widetilde{{\bf v}}\cdot\nabla\widetilde{{\bf v}}+\widetilde{\nabla}\widetilde{p}=\triangle\widetilde{{\bf v}},\,\,\,\widetilde{\mathrm{div}\,}\widetilde{{\bf v}}=0.

To take into account of our model (5)–(7), we first note that for any order two tensors (or tensor fields) $\mathsf{R}$ and $\mathsf{S}$ , suppose $\widetilde{\mathsf{R}}$ and $\widetilde{\mathsf{S}}$ are such that $\mathsf{R}(x)={\mathbb{A}}\widetilde{\mathsf{R}}(y){\mathbb{B}}$ and $\mathsf{S}(x)={\mathbb{A}}\widetilde{\mathsf{S}}(y){\mathbb{B}}$ , then it holds that

\mathsf{R}\cdot\mathsf{S}=\widetilde{\mathbb{R}}\cdot\widetilde{\mathsf{S}},\,\,\,\text{and}\,\,\,\mathrm{div}\,\mathsf{R}={\mathbb{A}}\,\widetilde{\mathrm{div}\,}\widetilde{\mathsf{R}}.

As application, consider $\mathsf{A}=\nabla{\bf v}+\nabla{\bf v}^{T}$ , $\mathsf{W}=\nabla{\bf v}-\nabla{\bf v}^{T}$ , $\widetilde{\mathsf{A}}=\widetilde{\nabla}\widetilde{{\bf v}}+\widetilde{\nabla}\widetilde{{\bf v}}^{T}$ and $\widetilde{\mathsf{W}}=\widetilde{\nabla}\widetilde{{\bf v}}-\widetilde{\nabla}\widetilde{{\bf v}}^{T}$ , then we have

\mathsf{A}={\mathbb{A}}\widetilde{\mathsf{A}}{\mathbb{B}},\,\,\,\text{and}\,\,\,\mathsf{W}={\mathbb{A}}\widetilde{\mathsf{W}}{\mathbb{B}}.

For $\mathsf{Q}$ , let $\widetilde{\mathsf{Q}}$ be such that $\mathsf{Q}(x)={\mathbb{A}}\widetilde{\mathsf{Q}}(y){\mathbb{B}}$ . In particular, if $\mathsf{Q}=s_{*}({n_{*}}\otimes{n_{*}}-\frac{\mathsf{I}}{3})$ , then $\widetilde{\mathsf{Q}}=s_{*}(({\mathbb{B}}{n_{*}})\otimes({\mathbb{B}}{n_{*}})-\frac{\mathsf{I}}{3})$ . In general, the following hold,

\mathcal{E}_{\mathrm{LdG}}(\nabla\mathsf{Q},\mathsf{Q})=\mathcal{E}_{\mathrm{LdG}}(\widetilde{\nabla}\widetilde{\mathsf{Q}},\widetilde{\mathsf{Q}}),\,\,\,\mathcal{F}_{\mathrm{LdG}}(\mathsf{Q})=\mathcal{F}_{\mathrm{LdG}}(\widetilde{\mathsf{Q}}),

and

\mathring{\mathsf{Q}}(x)=\mathsf{Q}_{t}+{\bf v}\cdot\nabla\mathsf{Q}+\mathsf{Q}\mathsf{W}-\mathsf{W}\mathsf{Q}={\mathbb{A}}\mathring{\widetilde{\mathsf{Q}}}(y){\mathbb{B}}.

Hence for $\mathsf{T}_{\text{SV}}^{\text{v}}=\mathsf{T}_{\text{SV}}^{\text{v}}(\mathring{\mathsf{Q}},\mathsf{Q},\mathsf{A})$ and $\mathsf{T}^{\text{el}}=\mathsf{T}^{\text{el}}(\nabla\mathsf{Q})$ , upon defining $\widetilde{\mathsf{T}}_{\text{SV}}^{\text{v}}=\mathsf{T}_{\text{SV}}^{\text{v}}(\mathring{\widetilde{\mathsf{Q}}},\widetilde{\mathsf{Q}},\widetilde{\mathsf{A}})$ and $\widetilde{\mathsf{T}}^{\text{el}}=\mathsf{T}^{\text{el}}(\widetilde{\nabla}\widetilde{\mathsf{Q}})$ , from (10) and (11), it holds that

\mathsf{T}_{\text{SV}}^{\text{v}}={\mathbb{A}}\widetilde{\mathsf{T}}_{\text{SV}}^{\text{v}}{\mathbb{B}},\,\,\,\,\,\,\mathsf{T}^{\mathrm{el}}={\mathbb{A}}\widetilde{\mathsf{T}}^{\mathrm{el}}{\mathbb{B}},\,\,\,\,\,\,\mathrm{div}\,\mathsf{T}_{\text{SV}}^{\text{v}}={\mathbb{A}}\,\widetilde{\mathrm{div}\,}\widetilde{\mathsf{T}}_{\text{SV}}^{\text{v}},\,\,\,\,\,\,\mathrm{div}\,\mathsf{T}^{\text{el}}={\mathbb{A}}\,\widetilde{\mathrm{div}\,}\widetilde{\mathsf{T}}^{\text{el}}.

Furthermore, let $\widetilde{\mathcal{E}}_{\mathrm{LdG}}(\widetilde{\nabla}\widetilde{\mathsf{Q}},\widetilde{\mathsf{Q}})=\mathcal{E}_{\mathrm{LdG}}(\widetilde{\nabla}\widetilde{\mathsf{Q}},\widetilde{\mathsf{Q}})$ , from (8) and (9), then

\displaystyle\frac{\partial\mathcal{E}_{\mathrm{LdG}}}{\partial\mathsf{Q}}={\mathbb{A}}\frac{\partial\widetilde{\mathcal{E}}_{\mathrm{LdG}}}{\partial\widetilde{\mathsf{Q}}}{\mathbb{B}}\,\,\,\,\,\,\text{and}\,\,\,\,\,\,\mathrm{div}\,\left[\frac{\partial\mathcal{E}_{\mathrm{LdG}}}{\partial\nabla\mathsf{Q}}\right]={\mathbb{A}}\,\mathrm{div}\,\left[\frac{\partial\widetilde{\mathcal{E}}_{\mathrm{LdG}}}{\partial\widetilde{\nabla}\widetilde{\mathsf{Q}}}\right]{\mathbb{B}}.

We then conclude that system (5)–(7) is invariant under orthogonal transformations. In particular, we say that a solution ${\bf v}$ , $\mathsf{Q}$ is invariant under an orthogonal transformation ${\mathbb{B}}$ if

	$\displaystyle\widetilde{{\bf v}}(x)={\bf v}(x),$		$\displaystyle\text{or equivalently,}\,\,\,\,\,\,{\bf v}({\mathbb{B}}x)={\mathbb{B}}{\bf v}(x),$		(172)
	$\displaystyle\widetilde{\mathsf{Q}}(x)=\mathsf{Q}(x),$		$\displaystyle\text{or equivalently,}\,\,\,\,\,\,\mathsf{Q}({\mathbb{B}}x)={\mathbb{B}}\mathsf{Q}(x){\mathbb{A}}.$		(173)

Next we investigate the symmetry property of the drag force $\mathcal{F}_{\gamma}$ on the particle which according to (6) is defined as

\mathcal{F}_{\gamma}:=\int_{\partial{\bf B}}\mathsf{T}\nu\,d\sigma_{x},\,\,\,\,\,\,\text{where}\,\,\,\mathsf{T}=-p\mathsf{I}+\mathsf{T}_{\text{SV}}^{\text{v}}+\mathsf{T}^{\text{el}}.

(174)

By the invariance property just demonstrated, we have

\displaystyle\mathcal{F}_{\gamma}=\int_{\partial{\bf B}}{\mathbb{A}}\widetilde{\mathsf{T}}{\mathbb{B}}\nu_{x}\,d\tilde{\sigma}_{y}={\mathbb{A}}\int_{\partial\widetilde{{\bf B}}}\widetilde{\mathsf{T}}\widetilde{\nu}_{y}\,d\tilde{\sigma}_{y}=:{\mathbb{A}}\widetilde{\mathcal{F}}_{\gamma},\,\,\,\,\,\,\text{where $\widetilde{\mathsf{T}}=-\widetilde{p}\mathsf{I}+\widetilde{\mathsf{T}}_{\text{SV}}^{\text{v}}+\widetilde{\mathsf{T}}^{\text{el}}$.}

(175)

Note that in the above, we have used the fact that the area element does not change under orthogonal transformation. Now if the solution is invariant under ${\mathbb{B}}$ , then $\widetilde{\mathsf{T}}(y)=\mathsf{T}(y)$ . Hence

\widetilde{\mathcal{F}}_{\gamma}=\int_{\partial\widetilde{{\bf B}}}\widetilde{\mathsf{T}}\widetilde{\nu}_{y}\,d\tilde{\sigma}_{y}=\int_{\partial{\bf B}}\mathsf{T}\nu_{x}\,d\sigma_{x}=\mathcal{F}_{\gamma}.

This leads to

\mathcal{F}_{\gamma}={\mathbb{A}}\mathcal{F}_{\gamma},\,\,\,\,\,\,\text{or equivalently,}\,\,\,\,\,\,\mathcal{F}_{\gamma}={\mathbb{B}}\mathcal{F}_{\gamma}.

(176)

Next we make use of the above to analyze the symmetries of the drag force $\mathcal{J}_{\gamma}$ . For convenience, we use $\mathcal{F}_{\gamma}({{\bf v}_{*}},\mathsf{Q}_{*})$ to denote the dependence of $\mathcal{F}_{\gamma}$ on ${{\bf v}_{*}}$ and $\mathsf{Q}_{*}$ . Note that for our linear system (111)–(116), with ${\bf v}_{b}=0$ , given $\mathsf{Q}_{*}$ , its solution ${\bf v}$ and hence $\mathcal{F}_{\gamma}$ is linear in the far-field velocity field ${\bf v}_{*}$ . To take advantage of this, for concreteness, we let ${n_{*}}=e_{3}$ in $\mathsf{Q}_{*}$ (39). Now we decompose ${\bf v}_{*}$ as,

{{\bf v}_{*}}={{\bf v}_{*}}^{\parallel}+{{\bf v}_{*}}^{\perp},\quad\text{where}\,\,\,{{\bf v}_{*}}^{\parallel}\parallel e_{3},\,\,\,{{\bf v}_{*}}^{\perp}\perp e_{3},\,\,\,\text{i.e. ${{\bf v}_{*}}^{\perp}$ lies in the $xy$-plane.}

(177)

Then the solution ${\bf v}$ and the drag force $\mathcal{F}_{\gamma}$ can be decomposed as

{\bf v}={\bf v}^{\parallel}+{\bf v}^{\perp},\quad\mathcal{F}_{\gamma}({{\bf v}_{*}},\mathsf{Q}_{*})=\mathcal{F}_{\gamma}({{\bf v}_{*}}^{\parallel},\mathsf{Q}_{*})+\mathcal{F}_{\gamma}({{\bf v}_{*}}^{\perp},\mathsf{Q}_{*})

where ${\bf v}^{\parallel},{\bf v}^{\perp}$ solve (111)–(116) with ${{\bf v}_{*}}^{\parallel},{{\bf v}_{*}}^{\perp}$ as the far-field velocity field. Furthermore, we have

\mathbb{R}\mathcal{F}_{\gamma}({{\bf v}_{*}},\mathsf{Q}_{*})=\mathcal{F}_{\gamma}(\mathbb{R}{{\bf v}_{*}},\mathbb{R}\mathsf{Q}_{*}\mathbb{R}^{-1})

(178)

for any orthogonal transformation $\mathbb{R}$ of $\mathbb{R}^{3}$ . The following arguments reveal the symmetry properties of $\mathcal{F}_{\gamma}$ with respect to ${{\bf v}_{*}}^{\parallel}$ , ${{\bf v}_{*}}^{\perp}$ .

•

Since $\mathsf{Q}_{*}$ and ${{\bf v}_{*}}^{\parallel}$ are invariant under any orthogonal transformation of $\mathbb{R}^{3}$ that leaves $e_{3}$ fixed, we can infer that ${\bf v}^{\parallel}$ satisfies the same property. Hence the drag force $\mathcal{F}_{\gamma}({{\bf v}_{*}}^{\parallel},\mathsf{Q}_{*})$ must be parallel to $e_{3}$ . By linearity, we have for some constant $\gamma_{\parallel}$ that

\mathcal{F}_{\gamma}({{\bf v}_{*}}^{\parallel},\mathsf{Q}_{*})=\gamma_{\parallel}\langle{{\bf v}_{*}}^{\parallel},e_{3}\rangle e_{3}=\gamma_{\parallel}\big{(}e_{3}\otimes e_{3}){{\bf v}_{*}}.

(179)

•

For ${{\bf v}_{*}}^{\perp}$ , we have in general that

\mathbb{R}^{\prime}\mathcal{F}_{\gamma}({{\bf v}_{*}}^{\perp},\mathsf{Q}_{*})=\mathcal{F}_{\gamma}(\mathbb{R}^{\prime}{{\bf v}_{*}}^{\perp},\mathbb{R}^{\prime}\mathsf{Q}_{*}{\mathbb{R}^{\prime}}^{-1})=\mathcal{F}_{\gamma}({{\bf v}_{*}}^{\perp},\mathsf{Q}_{*})

(180)

for any reflection $\mathbb{R}^{\prime}$ of $\mathbb{R}^{3}$ that leaves the plane spanned by $e_{3}$ and ${{\bf v}_{*}}^{\perp}$ fixed. Hence we must have $\mathcal{F}_{\gamma}({{\bf v}_{*}}^{\perp},\mathsf{Q}_{*})\in\text{Span}\{e_{3},{{\bf v}_{*}}^{\perp}\}$ . By linearity, we have for some constant $\gamma_{\perp}$ and vector ${\bf C}$ from the $xy$ -plane that

\mathcal{F}_{\gamma}({{\bf v}_{*}}^{\perp},\mathsf{Q}_{*})=\gamma_{\perp}{{\bf v}_{*}}^{\perp}+\langle{\bf C},{{\bf v}_{*}}^{\perp}\rangle e_{3}.

Now for any orthogonal transformation $\mathbb{R}^{\prime\prime}$ of $\mathbb{R}^{3}$ that leaves $e_{3}$ fixed, using

\mathbb{R}^{\prime\prime}\mathcal{F}_{\gamma}({{\bf v}_{*}}^{\perp},\mathsf{Q}_{*})=\mathcal{F}_{\gamma}(\mathbb{R}^{\prime\prime}{{\bf v}_{*}}^{\perp},\mathbb{R}^{\prime\prime}\mathsf{Q}_{*}{\mathbb{R}^{\prime\prime}}^{-1}),

we have

\gamma_{\perp}\mathbb{R}^{\prime\prime}{{\bf v}_{*}}^{\perp}+\langle{\bf C},{{\bf v}_{*}}^{\perp}\rangle e_{3}=\gamma_{\perp}\mathbb{R}^{\prime\prime}{{\bf v}_{*}}^{\perp}+\langle{\bf C},\mathbb{R}^{\prime\prime}{{\bf v}_{*}}^{\perp}\rangle e_{3}.

Hence $\langle{\bf C},{{\bf v}_{*}}^{\perp}\rangle=\langle{\bf C},\mathbb{R}^{\prime\prime}{{\bf v}_{*}}^{\perp}\rangle$ for any $\mathbb{R}^{\prime\prime}$ . Thus ${\bf C}$ must be zero. We then conclude that

\mathcal{F}_{\gamma}({{\bf v}_{*}}^{\perp},\mathsf{Q}_{*})=\gamma_{\perp}{{\bf v}_{*}}^{\perp}=\gamma_{\perp}\big{(}\mathsf{I}-e_{3}\otimes e_{3}\big{)}{{\bf v}_{*}}.

(181)

Combining the above, we finally have the following same formula as [23, (6.9)]:

\mathcal{F}_{\gamma}({{\bf v}_{*}},\mathsf{Q}_{*})=\mathcal{F}_{\gamma}({{\bf v}_{*}}^{\parallel},\mathsf{Q}_{*})+\mathcal{F}_{\gamma}({{\bf v}_{*}}^{\perp},\mathsf{Q}_{*})=\Big{[}\gamma_{\parallel}\big{(}e_{3}\otimes e_{3}\big{)}+\gamma_{\perp}\big{(}\mathsf{I}-e_{3}\otimes e_{3}\big{)}\Big{]}{{\bf v}_{*}}.

(182)

This is also consistent with the fact that if we replace ${n_{*}}$ by $-{n_{*}}$ , $\mathsf{Q}_{*}$ remains unchanged and so does the the overall system.

Even though our system of equation is a reduced version of the original model, it does not seem easy to write down an asymptotic formula for the coefficients $\gamma_{\perp}$ and $\gamma_{\parallel}$ . This is because of the presence of $\frac{1}{\text{Er}}$ in the stress tensor (20) which cannot be easily computed in the limit of vanishing Er. See [23] for a discussion of some analysis and conjectures about the drag force.

6.3 Analytical calculations

We demonstrate here analytical calculations of (6.1.1), (6.1.1) and (171). In principle, we can give analytical expressions for all the terms as they only involve homogeneous functions of negative integral degrees. More precisely, we have

	$\displaystyle\mathsf{E}(x)\,\,\text{\eqref{GreenE}}:$		$\displaystyle\frac{\delta_{ij}}{r},\,\frac{x_{i}x_{j}}{r^{3}},$
	$\displaystyle\mathsf{F}(x)\,\,\text{\eqref{GreenF}}:$		$\displaystyle\frac{\delta_{ij}}{r^{3}},\,\frac{x_{i}x_{j}}{r^{5}}$
	$\displaystyle\mathsf{Q}\,\,\text{\eqref{QABL}}:$		$\displaystyle\frac{\delta_{ij}}{r},\,\frac{\delta_{ij}}{r^{3}},\,\frac{x_{i}x_{j}}{r^{5}}$

and ${\bf v}=\mathsf{E}_{S}{\bf v}_{*}$ with $\mathsf{E}_{S}=\mathsf{I}-6\pi a\mathsf{E}(x)+\frac{a^{3}}{4}\mathsf{F}(x)$ (67). Furthermore, the terms $\mathcal{C}_{\gamma}$ and $\mathcal{D}_{\gamma}$ involve multiplications between the following matrices

\mathsf{A},\,\,\mathsf{W},\,\,\mathsf{Q}_{*},\,\,\mathsf{Q},\,\,\,{\bf v}_{*}\cdot\nabla\mathsf{Q},\,\,{\bf v}\cdot\nabla\mathsf{Q}.

For convenience, we introduce the following conventions:

$\displaystyle[A,B,C,\ldots]$	$\displaystyle=$	arbitrary linear combinations and products between $A$ , $B$ , $C$ ,	(183)
		and their powers;	(183)
$\displaystyle\{A,B,C,\ldots\}$	$\displaystyle=$	linear combinations between $A$ , $B$ , $C$ .	(184)

From Appendix C.1, we have

	$\displaystyle\mathcal{M}_{\gamma}:D^{2}\mathsf{E}(x)$	$\displaystyle\in$	$\displaystyle\frac{1}{r^{3}}\big{[}1,\langle{n_{}},\hat{x}\rangle\big{]}\Big{\{}\mathsf{I},\,\,\,{n_{}}\otimes{n_{}},\,\,\,{n_{}}\otimes\hat{x},\,\,\,\hat{x}\otimes{n_{*}},\,\,\,\hat{x}\otimes\hat{x}\Big{\}}$		(185)
	$\displaystyle\mathcal{M}_{\gamma}:D^{2}\mathsf{F}(x)$	$\displaystyle\in$	$\displaystyle\frac{1}{r^{5}}\big{[}1,\langle{n_{}},\hat{x}\rangle\big{]}\Big{\{}\mathsf{I},\,\,\,{n_{}}\otimes{n_{}},\,\,\,{n_{}}\otimes\hat{x},\,\,\,\hat{x}\otimes{n_{*}},\,\,\,\hat{x}\otimes\hat{x}\Big{\}}$		(186)

so that

	$\displaystyle\mathcal{M}_{\gamma}:D^{2}\mathsf{E}(x){\bf v}_{*}$	$\displaystyle\in$	$\displaystyle\frac{1}{r^{3}}\Big{[}1,\,\,\,\langle{n_{}},\hat{x}\rangle,\,\,\,\langle{n_{}},{\bf v}_{}\rangle,\,\,\,\langle\hat{x},{\bf v}_{}\rangle\Big{]}\Big{\{}{\bf v}_{},\,\,\,{n_{}},\,\,\,\hat{x}\Big{\}}$		(187)
	$\displaystyle\mathcal{M}_{\gamma}:D^{2}\mathsf{F}(x){\bf v}_{*}$	$\displaystyle\in$	$\displaystyle\frac{1}{r^{5}}\Big{[}1,\,\,\,\langle{n_{}},\hat{x}\rangle,\,\,\,\langle{n_{}},{\bf v}_{}\rangle,\,\,\,\langle\hat{x},{\bf v}_{}\rangle\Big{]}\Big{\{}{\bf v}_{},\,\,\,{n_{}},\,\,\,\hat{x}\Big{\}}.$		(188)

From Appendix C.2 and C.3, we have

$\displaystyle\mathrm{div}\,\mathcal{A}_{\gamma}$	$\displaystyle\in$	$\displaystyle\left\{\frac{1}{r^{3}},\frac{1}{r^{4}}\right\}\Big{[}1,\,\,\,\langle\hat{x},{n_{}}\rangle,\,\,\,\langle\hat{x},{\bf v}_{}\rangle,\,\,\,\langle{n_{}},{\bf v}_{}\rangle\Big{]}\Big{\{}\hat{x},\,\,{n_{}},\,\,{\bf v}_{}\Big{\}},$	(189)
$\displaystyle\mathcal{C}_{\gamma}$	$\displaystyle\in$	$\displaystyle\left[\frac{1}{r^{3}},\ldots,\frac{1}{r^{6}}\right]\Big{[}1,\,\,\langle\hat{x},{n_{}}\rangle,\,\,\langle\hat{x},{\bf v}_{}\rangle,\,\,\langle{n_{}},{\bf v}_{}\rangle\Big{]}\times$
		$\displaystyle\Big{\{}\mathsf{I},\,\,\,{n_{}}\otimes{n_{}},\,\,\,{n_{}}\otimes{\bf v}_{},\,\,\,{\bf v}_{}\otimes{n_{}},\,\,\,{n_{}}\otimes\hat{x},\,\,\,\hat{x}\otimes{n_{}},\,\,\,{\bf v}_{}\otimes\hat{x},\,\,\,\hat{x}\otimes{\bf v}_{},\,\,\,\hat{x}\otimes\hat{x}\Big{\}},$
$\displaystyle\mathcal{D}_{\gamma}$	$\displaystyle\in$	$\displaystyle\left[\frac{1}{r^{4}},\ldots,\frac{1}{r^{9}}\right]\Big{[}1,\,\,\,\langle\hat{x},{n_{}}\rangle,\,\,\,\langle\hat{x},{\bf v}_{}\rangle,\,\,\,\langle{n_{}},{\bf v}_{}\rangle\Big{]}\Big{\{}\hat{x},\,\,{n_{}},\,\,{\bf v}_{}\Big{\}}.$	(191)

Note that all the terms in $\mathcal{M}_{\gamma}:D^{2}\mathsf{E}$ , $\mathcal{M}_{\gamma}:D^{2}\mathsf{F}$ and $\mathrm{div}\,\mathcal{A}_{\gamma}$ have already been computed explicitly in the Appendix C.1 and C.2. From Appendix C.3, all the coefficients appearing in $\mathcal{C}_{\gamma}$ and $\mathcal{D}_{\gamma}$ are amenable for symbolic computations.

Using the above, we can give the following representations of the bulk and boundary terms appearing in (6.1.1), (6.1.1) and (171).

•

$\displaystyle\int_{\Omega}\mathsf{E}(x-y)\Big{[}6\pi a\mathcal{M}_{\gamma}:D^{2}\mathsf{E}(y){\bf v}_{*}+\mathrm{div}\,\mathcal{A}_{\gamma}(y)\Big{]}\,dy$ . By (222), we have

			$\displaystyle\int_{\Omega}\mathsf{E}(x-y)\Big{[}6\pi a\mathcal{M}_{\gamma}:D^{2}\mathsf{E}(y){\bf v}_{*}+\mathrm{div}\,\mathcal{A}_{\gamma}(y)\Big{]}\,dy$
		$\displaystyle=$	$\displaystyle\frac{1}{\|x\|}\int_{\mathbb{S}^{2}}\mathsf{H}(\hat{x},\hat{y})\Big{[}6\pi a\mathcal{M}_{\gamma}:D^{2}\mathsf{E}(\hat{y}){\bf v}_{*}+\mathrm{div}\,\mathcal{A}_{\gamma}(\hat{y})\Big{]}\,d\sigma(\hat{y})+O\left(\frac{1}{\|x\|^{2}}\right)$
		$\displaystyle=$	$\displaystyle\frac{1}{\|x\|}\int_{\mathbb{S}^{2}}\mathsf{H}(\hat{x},\hat{y})\Big{[}A_{1}(\hat{y})\hat{y}+A_{2}(\hat{y}){n_{}}+A_{3}(\hat{y}){\bf v}_{}\Big{]}\,d\sigma(\hat{y})+O\left(\frac{1}{\|x\|^{2}}\right),$

where $\mathsf{H}$ given by (220), is recorded here,

\mathsf{H}(\hat{x},\hat{y})=\lim_{\epsilon\to 0}\left[\int_{\epsilon}^{\infty}\frac{\mathsf{E}(\hat{x}-r\hat{y})}{r}\,dr+\mathsf{E}(\hat{x})\ln(\epsilon)\right],

(192)

and the $A_{i}(\hat{y})$ are appropriate polynomial functions of $\langle\hat{y},{n_{*}}\rangle,\,\langle\hat{y},{\bf v}_{*}\rangle,\,\langle{n_{*}},{\bf v}_{*}\rangle$ .

•

$\displaystyle\int_{\Omega}\left[-\frac{a^{3}}{4}\mathcal{M}_{\gamma}:D^{2}\mathsf{F}(y){\bf v}_{*}+\mathrm{div}\,\mathcal{C}_{\gamma}(y)+\mathcal{D}_{\gamma}(y)\right]\,dy$ . The terms $D^{2}\mathsf{F}$ and $\mathcal{D}_{\gamma}$ decay at least $r^{-4}$ and hence are integrable. Using their forms, we have

\displaystyle\int_{\Omega}\left[-\frac{a^{3}}{4}\mathcal{M}_{\gamma}:D^{2}\mathsf{F}(y){\bf v}_{*}+\mathcal{D}_{\gamma}(y)\right]\,dy=\int_{\mathbb{S}^{2}}\Big{[}A_{4}(\hat{y})\hat{y}+A_{5}(\hat{y}){n_{*}}+A_{6}(\hat{y}){\bf v}_{*}\Big{]}\,d\sigma(\hat{y}).

For $\mathrm{div}\,\mathcal{C}_{\gamma}=O\left(\frac{1}{r^{3}}\right)$ , it can be conveniently represented using Divergence Theorem:

\int_{\Omega}\mathrm{div}\,\mathcal{C}_{\gamma}(y)\,dy=\int_{\mathbb{S}^{2}}\Big{\langle}\mathcal{C}_{\gamma}(\hat{y}),\hat{y}\Big{\rangle}\,d\sigma(\hat{y})=\int_{\mathbb{S}^{2}}\Big{[}A_{7}(\hat{y})\hat{y}+A_{8}(\hat{y}){n_{*}}+A_{9}(\hat{y}){\bf v}_{*}\Big{]}\,d\sigma(\hat{y}).

(194)

•

$\displaystyle\int_{\Omega}(\mathsf{E}_{S}-\mathsf{I})\Big{[}-\mathcal{M}_{\gamma}:D^{2}\mathsf{E}_{S}(y){\bf v}_{*}+\mathrm{div}\,\mathcal{A}_{\gamma}(y)+\mathrm{div}\,\mathcal{C}_{\gamma}(\mathsf{E}_{S}(y){\bf v}_{*},y)+\mathcal{D}_{\gamma}(\mathsf{E}_{S}(y){\bf v}_{*},y)\Big{]}\,dy$ . As $\mathsf{E}_{S}-\mathsf{I}$ decays at least as $r^{-1}$ , the following integral is integrable:

			$\displaystyle\int_{\Omega}(\mathsf{E}_{S}-\mathsf{I})\Big{[}-\mathcal{M}_{\gamma}:D^{2}\mathsf{E}_{S}(y){\bf v}_{}+\mathrm{div}\,\mathcal{A}_{\gamma}({\bf v}_{},y)+\mathcal{D}_{\gamma}(\mathsf{E}_{S}(y){\bf v}_{*},y)\Big{]}\,dy$		(195)
		$\displaystyle=$	$\displaystyle\int_{\mathbb{S}^{2}}\Big{[}A_{10}(\hat{y})\hat{y}+A_{11}(\hat{y}){n_{}}+A_{12}(\hat{y}){\bf v}_{}\Big{]}\,d\sigma(\hat{y}),$		(196)

where we recall the form of $\mathsf{E}_{S}$ from (67). The remaining term with $\mathrm{div}\,\mathcal{C}_{\gamma}$ can also be dealt with using Divergence Theorem:

	$\displaystyle\int_{\Omega}(\mathsf{E}_{S}-\mathsf{I})\mathrm{div}\,\mathcal{C}_{\gamma}(\mathsf{E}_{S}(y){\bf v}_{*},y)\,dy$	(197)
$\displaystyle=$	$\displaystyle\int_{\mathbb{S}^{2}}\Big{\langle}(\mathsf{E}_{S}(\hat{y})-\mathsf{I}),\mathcal{C}_{\gamma}(\hat{y})\hat{y}\Big{\rangle}\,d\sigma(\hat{y})-\int_{\Omega}\Big{\langle}\nabla\mathsf{E}_{S}(y),\mathcal{C}_{\gamma}(y)\Big{\rangle}\,dy$
$\displaystyle=$	$\displaystyle\int_{\mathbb{S}^{2}}\Big{[}A_{13}(\hat{y})\hat{y}+A_{14}(\hat{y}){n_{}}+A_{15}(\hat{y}){\bf v}_{}\Big{]}\,d\sigma(\hat{y}).$

Before moving on to the next section, we consider one “simplistic model” in which $\mathsf{Q}$ is taken to be uniform in space, i.e. it equals its end state $\mathsf{Q}_{*}$ . Such an approximation was in fact used in some works, see for example [22, 37]. In this case, all the terms $\mathcal{A}_{\gamma}$ , $\mathcal{C}_{\gamma}$ and $\mathcal{D}_{\gamma}$ vanish as they involve either $\nabla\mathsf{Q}$ or $\mathsf{Q}-\mathsf{Q}_{*}$ . Then we have

$\displaystyle\mathcal{I}_{\gamma}(x)$	$\displaystyle=$	$\displaystyle\frac{1}{\|x\|}\int_{\mathbb{S}^{2}}\mathsf{H}(\hat{x},\omega)\Big{[}6\pi a\mathcal{M}_{\gamma}:D^{2}\mathsf{E}(\omega){\bf v}_{*}\Big{]}\,d\sigma_{\omega}$	(198)
		$\displaystyle+\mathsf{E}(x)\int_{\Omega}\left[-\frac{a^{3}}{4}\mathcal{M}_{\gamma}:D^{2}\mathsf{F}(y){\bf v}_{*}\right]\,dy,$	(198)
$\displaystyle\mathcal{J}_{\gamma}$	$\displaystyle=$	$\displaystyle\int_{\Omega}(\mathsf{E}_{S}-\mathsf{I})\Big{[}-\mathcal{M}_{\gamma}:D^{2}\mathsf{E}_{S}(y){\bf v}_{*}\Big{]}\,dy.$	(199)

As demonstrated numerically in Section 7.1, we see that the actual velocity flow ${\bf v}$ does depend on the overall structure of $\mathsf{Q}$ , not just its end state $\mathsf{Q}_{*}$ .

7 Numerical simulations

Here we provide numerical simulations to illustrate our analysis. The simulations are performed using a commercial finite elements software package COMSOL [38]. For validation, we used this package to compute the classical Stokes flow. The results are benchmarked against analytical solutions in a finite domain, more precisely in the annulus $a\leq r\leq R$ . For details, we refer to Appendix D.

In the following, we record our numerical results for the anisotropic Stokes system (32)-(34). Some remarks are in order.

(a).

For simplicity, we assume that only $\gamma_{1}$ and $\gamma_{2}$ are nonzero. They are fixed to be $\gamma_{1}=1$ and $\gamma_{2}=0.9$ , except in Section 7.3 where we allow them to vary.
(b).

Our analytical results show that ${\bf v}={\bf v}_{*}+O\left(\frac{1}{r}\right)$ . To better illustrate this, we will plot rescaled versions of components of ${\bf v}$ . More precisely, let $V=|{\bf v}_{*}|$ . If ${\bf v}_{*}$ is along $e_{i}$ , then we will plot the following quantity for $a<r<R:$

$r\left(\frac{{\bf v}_{i}}{V}-1\right)\,\,\,\text{and}\,\,\,r\left(\frac{{\bf v}_{j}}{V}\right)\,\,\,\text{for $j\neq i$}$ (200)

along various two-dimensional planes.
(c).

All our numerical results shows that the rescaled quantity (200) remains bounded in magnitude. Inevitably, due to boundary effect, the numerical solution is consistent with the true solution only for $a<r\ll R$ . This is made more precise in Appendix D.

7.1 ${\bf v}_{}=e_{1}$ , ${n_{}}=e_{3}$

In this section, we choose ${\bf v}_{*}$ and ${n_{*}}$ to be non-parallel to each other. Besides plotting various (rescaled) components of ${\bf v}_{i}$ ’s, we aim to illustrate the clear differences between the solution ${\bf v}$ of our anisotropic Stokes system when $\mathsf{Q}$ is set to be $\mathsf{Q}_{*}$ and given by (38).

(a) [Uncaptioned image] (b)