Asymptotic theory of a uniform flow of a rarefied gas past a sphere at low Mach numbers

Abstract

A slow uniform flow of a rarefied gas past a sphere with a uniform temperature is considered. The steady behaviour of the gas is investigated on the basis of the Boltzmann equation by a systematic asymptotic analysis for small Mach numbers in the case where the Knudsen number is finite. Introducing a slowly varying solution whose length scale of variation is much larger than the sphere dimension, the fluid-dynamic-type equations describing the overall behaviour of the gas in the far region are derived. Then, the solution in the near region which varies on the scale of the sphere size, described by the linearised Boltzmann equation, and the solution in the far region, described by the fluid-dynamic-type equations, are sought in the form of a Mach number expansion up to the second order, in a way that they are joined in the intermediate overlapping region. As a result, the drag is derived up to the second order of the Mach number, which formally extends the linear drag obtained by Takata et al. (Phys. Fluids A, vol. 5, 1993, pp. 716–737) to a weakly nonlinear case. Numerical results for the drag on the basis of the Bhatnagar–Gross–Krook (BGK) model are also presented.

1. Introduction

The problem of slow uniform flow (or low Reynolds number flow) past a sphere is a textbook example in fluid mechanics (e.g. Landau & Lifshitz 1987). It is also an important research topic in rarefied gas dynamics (Knudsen & Weber 1911; Epstein 1924; Willis 1966; Cercignani, Pagani & Bassanini 1968; Sone & Aoki 1977a,b; Law & Loyalka 1986; Aoki & Sone 1987; Beresnev, Chernyak & Fomyagin 1990; Loyalka 1992; Takata, Sone & Aoki 1993; Torrilhon 2010), not only because of its fundamental importance but also because of its relevance to aerosol sciences as well as to small-scale engineering. The success is also highlighted by the agreement between the theoretical predictions (e.g. Cercignani et al. 1968; Loyalka 1992; Takata et al. 1993) and Millikan’s experimental results (Millikan 1923).

On the other hand, up to now, most of the results have been based on the linearised Boltzmann equation (or on linearised kinetic models). As a consequence, very little is known about the nonlinear effect caused by the slow motion of a rarefied gas around a sphere. In particular, the existing formula for the drag is linear with respect to the flow speed (or the Mach number). It would seem to be useful if we could derive a higher-order (nonlinear) correction to the linear drag on the basis of the full Boltzmann equation, thereby improving the accuracy of the existing formula.

If we consider the problem of a slow flow (small Reynolds number flow) past a sphere in the framework of conventional fluid mechanics (the Navier–Stokes equation and the no-slip boundary condition), a perturbation theory was developed around the late 1950s (Kaplun & Lagerstrom 1957; Proudman & Pearson 1957; Van Dyke 1975). According to this theory, another length scale appears, besides the sphere dimension, which is much larger than the sphere, in the far region when the Reynolds number is small. The whole domain is then divided into two regions, the near and far regions, and the solutions in the respective regions are obtained in the form of a Reynolds number expansion by joining them in the intermediate overlapping region. As a result, an asymptotic formula for the drag is derived in the form of a Reynolds number expansion with the Stokes’ drag being the leading-order term. In view of this, it is worth trying a similar analysis for the Boltzmann system using the Mach number as a small parameter. In this paper, we carry out such an analysis.

It should be mentioned that the Boltzmann equation has an intrinsic length scale known as the mean free path of gas molecules. Therefore, the ratio between the mean free path and the size of the sphere (i.e. the Knudsen number) also enters the system as a parameter. In this paper, we restrict our consideration to the case where the Knudsen number is of the order of unity, leaving the case of small Knudsen numbers to a subsequent study.

2. Formulation

2.1. Problem and basic assumptions

Consider a slow uniform flow of a monatomic rarefied gas past a sphere with radius $L$. At a great distance from the sphere, the gas is in the equilibrium state with velocity $(v_{\infty },0,0)$, pressure $p_{\infty }$ and temperature $T_{\infty }$. The sphere is kept at the same temperature as the uniform flow and is assumed to be at rest. We investigate the steady behaviour of the gas around the sphere under the following assumptions:

1. (i) the behaviour of the gas is described by the Boltzmann equation;

2. (ii) the gas molecules are reflected from the sphere according to the rule described below;

3. (iii) the Mach number of the flow is small, i.e. $\mathit{Ma}=v_{\infty }/(5RT_{\infty }/3)^{1/2}\ll 1$;

4. (iv) the Knudsen number is finite, i.e. $\mathit{Kn}=\ell _{\infty }/L=O(1)$.

Here, $R$ is the specific gas constant ($R=k_{B}/m$ with $k_{B}$ and $m$ being the Boltzmann constant and the mass of a molecule, respectively) and $\ell _{\infty }$ is the mean free path of the gas molecules at the equilibrium state at rest with pressure $p_{\infty }$ and temperature $T_{\infty }$.

2.2. Basic equations

We first introduce some basic notation. Let $Lx_{i}$ (or $L\boldsymbol{x}$) be the space coordinates and $(2RT_{\infty })^{1/2}{\it\zeta}_{i}$ (or $(2RT_{\infty })^{1/2}{\bf\zeta}$) be the molecular velocity (see figure 1). We also use the spherical coordinate system $(Lr,{\it\theta},{\it\varphi})$ throughout the paper whenever it is convenient. The $r$, ${\it\theta}$ and ${\it\varphi}$ components of ${\bf\zeta}$ are denoted by ${\it\zeta}_{r}$, ${\it\zeta}_{{\it\theta}}$ and ${\it\zeta}_{{\it\varphi}}$, respectively. The relation ${\it\zeta}_{1}={\it\zeta}_{r}\cos {\it\theta}-{\it\zeta}_{{\it\theta}}\sin {\it\theta}$ is frequently used. We denote by ${\it\rho}_{\infty }(2RT_{\infty })^{-3/2}(1+{\it\phi}(\boldsymbol{x},{\bf\zeta}))E$ the velocity distribution function, where ${\it\rho}_{\infty }=p_{\infty }/RT_{\infty }$ is the density of the uniform flow and $E={\rm\pi}^{-3/2}\exp (-|{\bf\zeta}|^{2})$. Further, we introduce the density ${\it\rho}_{\infty }(1+{\it\omega}(\boldsymbol{x}))$, the flow velocity $(2RT_{\infty })^{1/2}u_{i}(\boldsymbol{x})$, the temperature $T_{\infty }(1+{\it\tau}(\boldsymbol{x}))$, the pressure $p_{\infty }(1+P(\boldsymbol{x}))$, the stress tensor $p_{\infty }({\it\delta}_{ij}+\unicode[STIX]{x1D617}_{ij}(\boldsymbol{x}))$ (${\it\delta}_{ij}$ is the Kronecker delta) and the heat-flux vector $p_{\infty }(2RT_{\infty })^{1/2}Q_{i}(\boldsymbol{x})$ of the gas. The components of $\boldsymbol{u}=(u_{i})_{i=1,2,3}$, $\unicode[STIX]{x1D64B}=(\unicode[STIX]{x1D617}_{ij})_{i,j=1,2,3}$ and $\boldsymbol{Q}=(Q_{i})_{i=1,2,3}$ in the spherical coordinate system are denoted by using $r$, ${\it\theta}$ and ${\it\varphi}$ as the subscripts, e.g. $u_{r}$, $u_{{\it\theta}}$, $\unicode[STIX]{x1D617}_{rr}$, etc.

Figure 1. Problem and coordinate system.

The Boltzmann equation for the present steady problem is written as

(2.1)$$\begin{eqnarray}{\it\zeta}_{i}\frac{\partial {\it\phi}}{\partial x_{i}}=\frac{1}{k}(\mathscr{L}[{\it\phi}]+\mathscr{J}[{\it\phi},{\it\phi}]).\end{eqnarray}$$

Here, $\mathscr{L}[{\it\phi}]$ and $\mathscr{J}[{\it\phi},{\it\phi}]$ are the linear and nonlinear parts of the collision operator (appendix A) and $k=(\sqrt{{\rm\pi}}/2)\mathit{Kn}$. The boundary condition on the sphere is written in the form

(2.2)$$\begin{eqnarray}{\it\phi}=\mathscr{K}[{\it\phi}]\quad \text{for}~{\it\zeta}_{i}n_{i}>0~(r=1),\end{eqnarray}$$

where $\mathscr{K}[{\it\phi}]$ is the (linear) scattering operator whose explicit form depends on the prescribed rule for the molecule–surface interaction and $n_{i}$ is the unit normal vector on the surface pointing to the gas. $\mathscr{K}[{\it\phi}]$ is further expressed by a scattering kernel $\hat{K}$ as

(2.3)$$\begin{eqnarray}\mathscr{K}[{\it\phi}]=\int _{{\it\zeta}_{i\ast }n_{i}<0}\hat{K}({\bf\zeta},{\bf\zeta}_{\ast }){\it\phi}_{\ast }E_{\ast }\,\text{d}{\bf\zeta}_{\ast }\quad ({\it\zeta}_{i}n_{i}>0),\end{eqnarray}$$

where ${\it\phi}_{\ast }={\it\phi}({\bf\zeta}_{\ast })$ and $E_{\ast }={\rm\pi}^{-3/2}\exp (-|{\bf\zeta}_{\ast }|^{2})$. Since the sphere is at rest and its surface temperature is uniform, $\hat{K}$ is independent of the local position on the surface (appendix A). The derivation of (2.3) from a more general scattering kernel is given in appendix A. For the diffuse reflection condition, $\hat{K}=-2\sqrt{{\rm\pi}}{\it\zeta}_{i\ast }n_{i}$. We assume that $\hat{K}$ satisfies the following properties.

1. (i) Positivity: $\hat{K}({\bf\zeta},{\bf\zeta}_{\ast })\geqslant 0$ for ${\it\zeta}_{i}n_{i}>0$ and ${\it\zeta}_{i\ast }n_{i}<0$.

2. (ii) Impermeability: $-\!\int _{{\it\zeta}_{i}n_{i}>0}({\it\zeta}_{k}n_{k}/{\it\zeta}_{j\ast }n_{j})\hat{K}({\bf\zeta},{\bf\zeta}_{\ast })E\,\text{d}{\bf\zeta}=1$ for ${\it\zeta}_{i\ast }n_{i}<0$.

3. (iii) Let

   (2.4)$$\begin{eqnarray}(1+{\it\phi})E=\frac{1+{\hat{c}}_{0}}{{\rm\pi}^{3/2}(1+{\hat{c}}_{4})^{3/2}}\exp \left(-\frac{({\it\zeta}_{i}-{\hat{c}}_{i})^{2}}{1+{\hat{c}}_{4}}\right)\end{eqnarray}$$
   be an arbitrary Maxwellian, where ${\hat{c}}_{0}$, ${\hat{c}}_{i}$, and ${\hat{c}}_{4}$ are independent of ${\bf\zeta}$. Among such ${\it\phi}$, only ${\it\phi}={\hat{c}}_{0}$ satisfies the relation ${\it\phi}=\mathscr{K}[{\it\phi}]$ for ${\it\zeta}_{i}n_{i}>0$.

4. (iv) Local isotropy:

   (2.5)$$\begin{eqnarray}\mathscr{K}[f(\unicode[STIX]{x1D62D}_{ij}{\it\zeta}_{j})]({\bf\zeta})=\mathscr{K}[f({\bf\zeta})](\unicode[STIX]{x1D62D}_{ij}{\it\zeta}_{j}),\end{eqnarray}$$
   where $\unicode[STIX]{x1D62D}_{ij}$ is any orthogonal transformation matrix ($\unicode[STIX]{x1D62D}_{ik}\unicode[STIX]{x1D62D}_{jk}={\it\delta}_{ij}$) that satisfies the condition $\unicode[STIX]{x1D62D}_{ij}n_{j}=n_{i}$.

Finally, the boundary condition at infinity is given by

(2.6)$$\begin{eqnarray}(1+{\it\phi})E\rightarrow (1+{\it\phi}_{\infty })E=\frac{1}{{\rm\pi}^{3/2}}\exp \left(-({\it\zeta}_{i}-\hat{v}_{\infty }{\it\delta}_{i1})^{2}\right)\quad \text{as}~r\rightarrow \infty ,\end{eqnarray}$$

where $\hat{v}_{\infty }=v_{\infty }/(2RT_{\infty })^{1/2}=\sqrt{5/6}\mathit{Ma}$.

The macroscopic quantities of interest, ${\it\omega}$, $u_{i}$, ${\it\tau}$, $P$, $\unicode[STIX]{x1D617}_{ij}$ and $Q_{i}$, are expressed in terms of ${\it\phi}$ as

(2.7a)$$\begin{eqnarray}\displaystyle & {\it\omega}=\langle {\it\phi}\rangle ,\quad (1+{\it\omega})u_{i}=\langle {\it\zeta}_{i}{\it\phi}\rangle , & \displaystyle\end{eqnarray}$$

(2.7b)$$\begin{eqnarray}\displaystyle & \displaystyle {\textstyle \frac{3}{2}}(1+{\it\omega}){\it\tau}=\left\langle \left(|{\bf\zeta}|^{2}-{\textstyle \frac{3}{2}}\right){\it\phi}\right\rangle -(1+{\it\omega})u_{j}^{2}, & \displaystyle\end{eqnarray}$$

(2.7c)$$\begin{eqnarray}\displaystyle & P={\it\omega}+{\it\tau}+{\it\omega}{\it\tau},\quad \unicode[STIX]{x1D617}_{ij}=2\langle {\it\zeta}_{i}{\it\zeta}_{j}{\it\phi}\rangle -2(1+{\it\omega})u_{i}u_{j}, & \displaystyle\end{eqnarray}$$

(2.7d)$$\begin{eqnarray}\displaystyle & \displaystyle Q_{i}=\langle {\it\zeta}_{i}|{\bf\zeta}|^{2}{\it\phi}\rangle -{\textstyle \frac{5}{2}}u_{i}-u_{j}\unicode[STIX]{x1D617}_{ij}-{\textstyle \frac{3}{2}}Pu_{i}-(1+{\it\omega})u_{i}u_{j}^{2}. & \displaystyle\end{eqnarray}$$

Here, the angle brackets represent the following moment with respect to ${\bf\zeta}$:

(2.8)$$\begin{eqnarray}\langle f\rangle =\int f({\bf\zeta})E\,\text{d}{\bf\zeta},\end{eqnarray}$$

where the integration is carried out over the whole space of ${\bf\zeta}$.

To summarise, the present boundary-value problem contains two basic parameters, i.e. $k$ (or $\mathit{Kn}$) and $\hat{v}_{\infty }$ (or $\mathit{Ma}$). Because of the assumptions (iii) and (iv), we assume that

(2.9)$$\begin{eqnarray}\hat{v}_{\infty }\ll 1,\quad k=O(1).\end{eqnarray}$$

The aim of the present study is to investigate the behaviour of the gas in the situation (2.9) by a systematic asymptotic analysis for small $\hat{v}_{\infty }\ll 1$. Finally, we note that the existence of a solution in the present situation was mathematically proved by Ukai & Asano (1983).

3. Asymptotic analysis: Whitehead’s paradox in kinetic theory

In the following sections, we carry out an asymptotic analysis of the boundary-value problem (2.1), (2.2) and (2.6) for small $\hat{v}_{\infty }\ll 1$ in the case where $k=O(1)$. Hereafter, we denote the small parameter $\hat{v}_{\infty }$ by

(3.1)$$\begin{eqnarray}{\it\epsilon}=\hat{v}_{\infty }(\ll 1),\end{eqnarray}$$

and use ${\it\epsilon}$ rather than $\hat{v}_{\infty }$.

First, we try to look for a solution of the problem (2.1), (2.2) and (2.6) whose length scale of variation is of the order of unity. The solution with this length scale is denoted by attaching the subscript $S$, i.e. $\partial {\it\phi}_{S}/\partial x_{i}=O({\it\phi}_{S})$, and will be called the inner solution for later convenience. We seek ${\it\phi}_{S}$ in the form of a simple power series of ${\it\epsilon}$, i.e.

(3.2)$$\begin{eqnarray}{\it\phi}_{S}={\it\phi}_{S(1)}{\it\epsilon}+{\it\phi}_{S(2)}{\it\epsilon}^{2}+\cdots \,.\end{eqnarray}$$

Here, the expansion starts from ${\it\epsilon}$ order because ${\it\phi}_{S}$ is already a perturbation. Corresponding to the expansion of ${\it\phi}_{S}$, the macroscopic quantities are also expanded in ${\it\epsilon}$ as

(3.3)$$\begin{eqnarray}h_{S}=h_{S(1)}{\it\epsilon}+h_{S(2)}{\it\epsilon}^{2}+\cdots \,,\end{eqnarray}$$

where $h={\it\omega}$, $u_{i}$, ${\it\tau}$, etc. The relations between $h_{S(m)}$ and ${\it\phi}_{S(m)}$ are obtained by substituting the expansion (3.2) into the definitions of the macroscopic quantities (2.7) with $h=h_{S}$ and ${\it\phi}={\it\phi}_{S}$ and gathering the terms with the same power of ${\it\epsilon}$.

If we substitute the expansion (3.2) into (2.1) with ${\it\phi}={\it\phi}_{S}$ and take into account the property $\partial {\it\phi}_{S}/\partial x_{i}=O({\it\phi}_{S})$, we obtain a sequence of linearised Boltzmann equations for ${\it\phi}_{S(m)}$, i.e.

(3.4)$$\begin{eqnarray}\displaystyle & \displaystyle {\it\zeta}_{i}\frac{\partial {\it\phi}_{S(1)}}{\partial x_{i}}=\frac{1}{k}\mathscr{L}[{\it\phi}_{S(1)}], & \displaystyle\end{eqnarray}$$

(3.5)$$\begin{eqnarray}\displaystyle & \displaystyle {\it\zeta}_{i}\frac{\partial {\it\phi}_{S(m)}}{\partial x_{i}}=\frac{1}{k}\mathscr{L}[{\it\phi}_{S(m)}]+\frac{1}{k}\mathop{\sum }_{r=1}^{m-1}\mathscr{J}[{\it\phi}_{S(r)},{\it\phi}_{S(m-r)}]\quad (m=2,3,\dots ). & \displaystyle\end{eqnarray}$$

Similarly, substituting the expansion into the boundary condition (2.2) and taking into account the linearity of the operator $\mathscr{K}$, we obtain the boundary condition for ${\it\phi}_{S(m)}~(m=1,2,\dots )$ as follows:

(3.6)$$\begin{eqnarray}{\it\phi}_{S(m)}=\mathscr{K}[{\it\phi}_{S(m)}]\quad \text{for}~{\it\zeta}_{i}n_{i}>0~\text{at}~r=1.\end{eqnarray}$$

For the boundary condition at infinity (cf. (2.6)), we first consider the following expansion of the uniform Maxwellian $(1+{\it\phi}_{\infty })E$ in ${\it\epsilon}$:

(3.7)$$\begin{eqnarray}{\it\phi}_{\infty }=\frac{1}{{\rm\pi}^{3/2}}\exp (-({\it\zeta}_{i}-{\it\epsilon}{\it\delta}_{i1})^{2})E^{-1}-1=2{\it\zeta}_{1}{\it\epsilon}+(2{\it\zeta}_{1}^{2}-1){\it\epsilon}^{2}+\cdots \,.\end{eqnarray}$$

Then, we require that ${\it\phi}_{S(m)}$ approaches the corresponding term of the expansion (3.7) as $r\rightarrow \infty$.

3.1. Problem for order ${\it\epsilon}$

At the leading order, the problem is reduced to

(3.8)$$\begin{eqnarray}\displaystyle & \displaystyle {\it\zeta}_{i}\frac{\partial {\it\phi}_{S(1)}}{\partial x_{i}}=\frac{1}{k}\mathscr{L}[{\it\phi}_{S(1)}], & \displaystyle\end{eqnarray}$$

(3.9)$$\begin{eqnarray}\displaystyle & {\it\phi}_{S(1)}=\mathscr{K}[{\it\phi}_{S(1)}]\quad \text{for}~{\it\zeta}_{i}n_{i}>0~(r=1), & \displaystyle\end{eqnarray}$$

(3.10)$$\begin{eqnarray}\displaystyle & {\it\phi}_{S(1)}\rightarrow 2{\it\zeta}_{1}\quad \text{as}~r\rightarrow \infty . & \displaystyle\end{eqnarray}$$

The macroscopic quantities $h_{S(1)}$ ($h={\it\omega}$, $u_{i}$, ${\it\tau}$, etc.) are expressed in terms of ${\it\phi}_{S(1)}$ as follows:

(3.11a)$$\begin{eqnarray}\displaystyle & {\it\omega}_{S(1)}=\langle {\it\phi}_{S(1)}\rangle ,\quad u_{iS(1)}=\langle {\it\zeta}_{i}{\it\phi}_{S(1)}\rangle , & \displaystyle\end{eqnarray}$$

(3.11b)$$\begin{eqnarray}\displaystyle & \displaystyle {\it\tau}_{S(1)}={\textstyle \frac{2}{3}}\left\langle \left(|{\bf\zeta}|^{2}-{\textstyle \frac{3}{2}}\right){\it\phi}_{S(1)}\right\rangle ,\quad P_{S(1)}={\it\omega}_{S(1)}+{\it\tau}_{S(1)}, & \displaystyle\end{eqnarray}$$

(3.11c)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D617}_{ijS(1)}=2\langle {\it\zeta}_{i}{\it\zeta}_{j}{\it\phi}_{S(1)}\rangle ,\quad Q_{iS(1)}=\langle {\it\zeta}_{i}|{\bf\zeta}|^{2}{\it\phi}_{S(1)}\rangle -{\textstyle \frac{5}{2}}u_{iS(1)}. & \displaystyle\end{eqnarray}$$

The problem (3.8)–(3.10) is a boundary-value problem for the linearised Boltzmann equation, and is nothing but the problem that has been investigated by many authors in the past (see e.g. Takata et al. 1993, and the references therein). Therefore, the property of the solution is well known. Hereafter, we call this problem the linearised problem for short.

3.2. Problem for order ${\it\epsilon}^{2}$

Suppose that ${\it\phi}_{S(1)}$ is known. Then the next-order problem takes the following form:

(3.12)$$\begin{eqnarray}\displaystyle & \displaystyle {\it\zeta}_{i}\frac{\partial {\it\phi}_{S(2)}}{\partial x_{i}}=\frac{1}{k}\mathscr{L}[{\it\phi}_{S(2)}]+\frac{1}{k}\mathscr{J}[{\it\phi}_{S(1)},{\it\phi}_{S(1)}], & \displaystyle\end{eqnarray}$$

(3.13)$$\begin{eqnarray}\displaystyle & {\it\phi}_{S(2)}=\mathscr{K}[{\it\phi}_{S(2)}]\quad \text{for}~{\it\zeta}_{i}n_{i}>0~(r=1), & \displaystyle\end{eqnarray}$$

(3.14)$$\begin{eqnarray}\displaystyle & {\it\phi}_{S(2)}\rightarrow 2{\it\zeta}_{1}^{2}-1\quad \text{as}~r\rightarrow \infty . & \displaystyle\end{eqnarray}$$

The relations between the macroscopic quantities $h_{S(2)}$ ($h={\it\omega}$, $u_{i}$, ${\it\tau}$, etc.) and ${\it\phi}_{S(2)}$ are summarised as follows:

(3.15a)$$\begin{eqnarray}\displaystyle & {\it\omega}_{S(2)}=\langle {\it\phi}_{S(2)}\rangle , & \displaystyle\end{eqnarray}$$

(3.15b)$$\begin{eqnarray}\displaystyle & u_{iS(2)}=\langle {\it\zeta}_{i}{\it\phi}_{S(2)}\rangle -{\it\omega}_{S(1)}u_{iS(1)}, & \displaystyle\end{eqnarray}$$

(3.15c)$$\begin{eqnarray}\displaystyle & \displaystyle {\it\tau}_{S(2)}={\textstyle \frac{2}{3}}\left\langle \left(|{\bf\zeta}|^{2}-{\textstyle \frac{3}{2}}\right){\it\phi}_{S(2)}\right\rangle -{\textstyle \frac{2}{3}}u_{jS(1)}^{2}-{\it\omega}_{S(1)}{\it\tau}_{S(1)}, & \displaystyle\end{eqnarray}$$

(3.15d)$$\begin{eqnarray}\displaystyle & P_{S(2)}={\it\omega}_{S(2)}+{\it\tau}_{S(2)}+{\it\omega}_{S(1)}{\it\tau}_{S(1)}, & \displaystyle\end{eqnarray}$$

(3.15e)$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D617}_{ijS(2)}=2\langle {\it\zeta}_{i}{\it\zeta}_{j}{\it\phi}_{S(2)}\rangle -2u_{iS(1)}u_{jS(1)}, & \displaystyle\end{eqnarray}$$

(3.15f)$$\begin{eqnarray}\displaystyle & \displaystyle Q_{iS(2)}=\langle {\it\zeta}_{i}|{\bf\zeta}|^{2}{\it\phi}_{S(2)}\rangle -{\textstyle \frac{5}{2}}u_{iS(2)}-u_{jS(1)}\unicode[STIX]{x1D617}_{ijS(1)}-{\textstyle \frac{3}{2}}P_{S(1)}u_{iS(1)}. & \displaystyle\end{eqnarray}$$

Equations (3.12)–(3.14) again form a boundary-value problem of the linearised Boltzmann equation for ${\it\phi}_{S(2)}$, but with an inhomogeneous term. However, it can be shown that ${\it\phi}_{S(2)}$ does not satisfy the boundary condition at infinity, except the case $k=\infty$, and therefore the problem admits no solution for $k<\infty$. This situation corresponds to what is known as Whitehead’s paradox (Van Dyke 1975) in conventional fluid mechanics.

The reason for this ill-posedness can be understood as follows. We are seeking a solution of the Boltzmann equation (2.1) whose magnitude is of the order of ${\it\epsilon}$. If the length scale of variation of the solution is of the order of unity, the nonlinear term of the Boltzmann equation remains much smaller than the other terms. On the other hand, if the length scale of variation is much longer and is of the order of $1/{\it\epsilon}$ (i.e. $\partial {\it\phi}/\partial x_{i}=O({\it\epsilon}{\it\phi})$), the streaming term ${\it\zeta}_{i}\partial {\it\phi}/\partial x_{i}$ becomes comparable to the nonlinear term. In this case, the cumulative effect of the nonlinear term cannot be neglected but plays a decisive role in describing the flow behaviour in the far region, irrespective of the smallness of ${\it\phi}$. In the next section, by taking this into account, we will investigate the behaviour of the gas in the far region as described by the Boltzmann equation.

We end this section by stating the following property of the leading-order inner solution ${\it\phi}_{S(1)}$ concerning its asymptotic behaviour at large $r$.

Proposition 1. Let ${\it\phi}_{S(1)}$ be a solution of the boundary-value problem (3.8)–(3.10), and let ${\it\omega}_{S(1)}$, $u_{iS(1)}$, ${\it\tau}_{S(1)}$ and $P_{S(1)}$ be the associated macroscopic quantities given by (3.11). Then, ${\it\phi}_{S(1)}$, ${\it\omega}_{S(1)}$, $u_{iS(1)}$, ${\it\tau}_{S(1)}$ and $P_{S(1)}$ approach the corresponding values at infinity at the following rate:

(3.16)$$\begin{eqnarray}\displaystyle & {\it\phi}_{S(1)}=2{\it\zeta}_{1}+O(r^{-1}), & \displaystyle\end{eqnarray}$$

(3.17)$$\begin{eqnarray}\displaystyle & u_{iS(1)}={\it\delta}_{i1}+O(r^{-1}), & \displaystyle\end{eqnarray}$$

(3.18a−c)$$\begin{eqnarray}{\it\omega}_{S(1)}=O(r^{-2}),\quad {\it\tau}_{S(1)}=O(r^{-2}),\quad P_{S(1)}=O(r^{-2}).\end{eqnarray}$$

Proof. The estimates are the direct consequences of the results in § 5.3 below. More specifically, (3.16) follows from (5.14) and (5.24), and (3.17) and (3.18) from (5.19) and (5.23) there.

4. Slowly varying solution in the far region

Motivated by the consideration in the preceding section, we now introduce a slowly varying solution with the length scale of variation much larger than unity (or the sphere dimension in dimensional space) in the region far from the sphere. We call this solution the outer solution and denote it by attaching the subscript $O$. The length scale of variation of ${\it\phi}_{O}$ is assumed to be of the order of $1/{\it\epsilon}$, i.e. $\partial {\it\phi}_{O}/\partial x_{i}=O({\it\epsilon}{\it\phi}_{O})$. The validity of this assumption is verified if such a solution is obtained.

In order to investigate the behaviour of the outer solution, it is convenient to introduce a new variable $y_{i}$ (the outer variable) related to $x_{i}$ (the inner variable) by the relation

(4.1)$$\begin{eqnarray}y_{i}={\it\epsilon}x_{i}.\end{eqnarray}$$

Assuming that the outer solution is a function of $y_{i}$ in space (i.e. ${\it\phi}_{O}={\it\phi}_{O}(y_{i},{\it\zeta}_{i})$), we can transform the equation into the form

(4.2)$$\begin{eqnarray}{\it\zeta}_{i}\frac{\partial {\it\phi}_{O}}{\partial y_{i}}=\frac{1}{k{\it\epsilon}}(\mathscr{L}[{\it\phi}_{O}]+\mathscr{J}[{\it\phi}_{O},{\it\phi}_{O}]).\end{eqnarray}$$

Leaving aside the boundary condition, we seek the outer solution in the form of a power series of ${\it\epsilon}$, i.e.

(4.3)$$\begin{eqnarray}{\it\phi}_{O}={\it\phi}_{O(1)}{\it\epsilon}+{\it\phi}_{O(2)}{\it\epsilon}^{2}+\cdots \,.\end{eqnarray}$$

Correspondingly, the macroscopic quantities are also expanded in ${\it\epsilon}$ as

(4.4)$$\begin{eqnarray}h_{O}=h_{O(1)}{\it\epsilon}+h_{O(2)}{\it\epsilon}^{2}+\cdots \,,\end{eqnarray}$$

where $h={\it\omega}$, $u_{i}$, ${\it\tau}$, etc. The relations between $h_{O(m)}$ and ${\it\phi}_{O(m)}$ are obtained by substituting the expansion (4.3) into the definitions of macroscopic quantities (2.7) with $h=h_{O}$ and ${\it\phi}={\it\phi}_{O}$ and collecting the terms of the same power of ${\it\epsilon}$. The results for $m=1$ and 2 are given by (3.11) and (3.15) with the subscript $S$ replaced by $O$, respectively.

If we substitute the expansion (4.3) into the Boltzmann equation (4.2) and take into account the property $\partial {\it\phi}_{O}/\partial y_{i}=O({\it\phi}_{O})$, we obtain a sequence of linear integral equations for ${\it\phi}_{O(m)}$:

(4.5a)$$\begin{eqnarray}\displaystyle & \mathscr{L}[{\it\phi}_{O(1)}]=0, & \displaystyle\end{eqnarray}$$

(4.5b)$$\begin{eqnarray}\displaystyle & \displaystyle \mathscr{L}[{\it\phi}_{O(m)}]=k{\it\zeta}_{i}\frac{\partial {\it\phi}_{O(m-1)}}{\partial y_{i}}-\mathop{\sum }_{r=1}^{m-1}\mathscr{J}[{\it\phi}_{O(m-r)},{\it\phi}_{O(r)}]\quad (m=2,3,\dots ). & \displaystyle\end{eqnarray}$$

The equations can be solved successively from the lowest order, provided that certain solvability conditions are satisfied.

The solution of (4.5a) is given by a linear combination of the collision invariants $(1,{\it\zeta}_{i},|{\bf\zeta}|^{2})$ or, equivalently,

(4.6)$$\begin{eqnarray}{\it\phi}_{O(1)}={\it\omega}_{O(1)}+2u_{iO(1)}{\it\zeta}_{i}+\left(|{\bf\zeta}|^{2}-{\textstyle \frac{3}{2}}\right){\it\tau}_{O(1)},\end{eqnarray}$$

if we take into account the relations between ${\it\phi}_{O(1)}$ and ${\it\omega}_{O(1)}$, $u_{iO(1)}$ and ${\it\tau}_{O(1)}$. On the other hand, (4.5b) is an inhomogeneous linear integral equation, whose homogeneous counterpart has a non-trivial solution. Therefore, in order for the equation to have a solution, the inhomogeneous term should satisfy the following solvability condition (Sone 2002):

(4.7)$$\begin{eqnarray}\frac{\partial }{\partial y_{j}}\langle {\it\psi}_{i}{\it\zeta}_{j}{\it\phi}_{O(m-1)}\rangle =0\quad (m=2,3,\dots ),\end{eqnarray}$$

where $({\it\psi}_{0},{\it\psi}_{i},{\it\psi}_{4})=(1,{\it\zeta}_{i},|{\bf\zeta}|^{2})$ are the collision invariants. Substituting ${\it\phi}_{O(r)}$ ($r=1,2,\dots$) into the condition, the fluid-dynamic-type equations describing the overall behaviour of the gas in the far region (i.e. the outer solution) are derived.

4.0.1. Fluid-dynamic-type equations and velocity distribution function for the outer solution

Since the derivation of the fluid-dynamic-type equations is well established and is expatiated (e.g. Sone 2002, 2007), we only summarise the results of the analysis here. First, the fluid-dynamic-type equations describing the behaviour of the gas in the far region are summarised as follows:

order ${\it\epsilon}$

(4.8)$$\begin{eqnarray}\frac{\partial P_{O(1)}}{\partial y_{i}}=0,\end{eqnarray}$$

(4.9a)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial u_{iO(1)}}{\partial y_{i}}=0, & \displaystyle\end{eqnarray}$$

(4.9b)$$\begin{eqnarray}\displaystyle & \displaystyle u_{jO(1)}\frac{\partial u_{iO(1)}}{\partial y_{j}}=-\frac{1}{2}\frac{\partial P_{O(2)}}{\partial y_{i}}+\frac{{\it\gamma}_{1}k}{2}\frac{\partial ^{2}u_{iO(1)}}{\partial y_{j}^{2}}, & \displaystyle\end{eqnarray}$$

(4.9c)$$\begin{eqnarray}\displaystyle & \displaystyle u_{jO(1)}\frac{\partial {\it\tau}_{O(1)}}{\partial y_{j}}=\frac{{\it\gamma}_{2}k}{2}\frac{\partial ^{2}{\it\tau}_{O(1)}}{\partial y_{j}^{2}}, & \displaystyle\end{eqnarray}$$

(4.9d)$$\begin{eqnarray}\displaystyle & P_{O(1)}={\it\omega}_{O(1)}+{\it\tau}_{O(1)}; & \displaystyle\end{eqnarray}$$

order ${\it\epsilon}^{2}$

(4.10a)$$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\partial u_{iO(2)}}{\partial y_{i}}=-u_{jO(1)}\frac{\partial {\it\omega}_{O(1)}}{\partial y_{j}},\end{eqnarray}$$

(4.10b)$$\begin{eqnarray}\displaystyle & & \displaystyle u_{jO(1)}\frac{\partial u_{iO(2)}}{\partial y_{j}}+({\it\omega}_{O(1)}u_{jO(1)}+u_{jO(2)})\frac{\partial u_{iO(1)}}{\partial y_{j}}\nonumber\\ \displaystyle & & \displaystyle \quad =-\frac{1}{2}\frac{\partial }{\partial y_{i}}\left(P_{O(3)}-\frac{{\it\gamma}_{1}{\it\gamma}_{2}-4{\it\gamma}_{3}}{6}k^{2}\frac{\partial ^{2}{\it\tau}_{O(1)}}{\partial y_{j}^{2}}\right)\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\frac{{\it\gamma}_{1}k}{2}\frac{\partial ^{2}u_{iO(2)}}{\partial y_{j}^{2}}+\frac{{\it\gamma}_{4}k}{2}\frac{\partial }{\partial y_{j}}\left[{\it\tau}_{O(1)}\left(\frac{\partial u_{iO(1)}}{\partial y_{j}}+\frac{\partial u_{jO(1)}}{\partial y_{i}}\right)\right],\end{eqnarray}$$

(4.10c)$$\begin{eqnarray}\displaystyle & & \displaystyle u_{jO(1)}\frac{\partial {\it\tau}_{O(2)}}{\partial y_{j}}+({\it\omega}_{O(1)}u_{jO(1)}+u_{jO(2)})\frac{\partial {\it\tau}_{O(1)}}{\partial y_{j}}-\frac{2}{5}u_{jO(1)}\frac{\partial P_{O(2)}}{\partial y_{j}}\nonumber\\ \displaystyle & & \displaystyle \quad =\frac{{\it\gamma}_{1}k}{5}\left(\frac{\partial u_{iO(1)}}{\partial y_{j}}+\frac{\partial u_{jO(1)}}{\partial y_{i}}\right)^{2}+\frac{k}{2}\frac{\partial ^{2}}{\partial y_{j}^{2}}\left({\it\gamma}_{2}{\it\tau}_{O(2)}+\frac{{\it\gamma}_{5}}{2}{\it\tau}_{O(1)}^{2}\right),\end{eqnarray}$$

(4.10d)$$\begin{eqnarray}\displaystyle & & \displaystyle P_{O(2)}={\it\omega}_{O(2)}+{\it\tau}_{O(2)}+{\it\omega}_{O(1)}{\it\tau}_{O(1)},\end{eqnarray}$$

where ${\it\gamma}_{i}$ ($i=1,\dots ,5$) are constants (dimensionless transport coefficients) whose definitions are given in appendix B.

Equation (4.8) shows that the leading-order pressure is constant. On the other hand, (4.9a) and (4.9b) are the continuity equation and the Navier–Stokes equation for an incompressible fluid, respectively, while (4.9c) is an energy equation. Equation (4.9d) is the linearised equation of state. Equations (4.10a–d) are similar to (4.9a–d) but contain some extra terms corresponding to the effect of compressibility and to a non-Navier–Stokes effect (i.e. the thermal stress, the term with ${\it\gamma}_{3}$). They become important when ${\it\tau}_{O(1)}$ (${\it\omega}_{O(1)}$) or $P_{O(2)}$ is not uniform. However, as we will see, they are constant in the present situation, and, therefore, these effects become of the higher order (see § 5.1). It should also be mentioned that the derived equations are essentially the same as the fluid-dynamic-type equations obtained by the S expansion (Sone 1971, 2002, 2007). The interested reader is referred to Sone (2002, 2007) for further discussions on the features of the fluid-dynamic-type equations.

Once the macroscopic quantities are obtained (by solving the fluid-dynamic-type equations), the velocity distribution functions ${\it\phi}_{O(m)}$ are readily obtained from the expressions

(4.11a)$$\begin{eqnarray}\displaystyle & & \displaystyle {\it\phi}_{O(1)}={\it\phi}_{eO(1)},\end{eqnarray}$$

(4.11b)$$\begin{eqnarray}\displaystyle & & \displaystyle {\it\phi}_{O(2)}={\it\phi}_{eO(2)}-k{\it\zeta}_{i}{\it\zeta}_{j}B(|{\bf\zeta}|)\displaystyle \frac{\partial u_{iO(1)}}{\partial y_{j}}-k{\it\zeta}_{i}A(|{\bf\zeta}|)\displaystyle \frac{\partial {\it\tau}_{O(1)}}{\partial y_{i}},\end{eqnarray}$$

$$\begin{eqnarray}\displaystyle \cdots \,. & & \displaystyle \nonumber\end{eqnarray}$$

The functions $A(|{\bf\zeta}|)$ and $B(|{\bf\zeta}|)$ are defined in appendix B and ${\it\phi}_{eO(m)}$ ($m=1,2,\dots$) are the expansion coefficients of the local Maxwellian in ${\it\epsilon}$, i.e.

(4.12)$$\begin{eqnarray}(1+{\it\phi}_{eO})E=\frac{1+{\it\omega}_{O}}{{\rm\pi}^{3/2}(1+{\it\tau}_{O})^{3/2}}\exp \left(-\frac{({\it\zeta}_{i}-u_{iO})^{2}}{1+{\it\tau}_{O}}\right)=E(1+{\it\phi}_{eO(1)}{\it\epsilon}+{\it\phi}_{eO(2)}{\it\epsilon}^{2}+\cdots \!).\end{eqnarray}$$

More specifically,

(4.13a)$$\begin{eqnarray}\displaystyle {\it\phi}_{eO(1)} & = & \displaystyle P_{O(1)}+2{\it\zeta}_{i}u_{iO(1)}+\left(|{\bf\zeta}|^{2}-{\textstyle \frac{5}{2}}\right){\it\tau}_{O(1)},\end{eqnarray}$$

(4.13b)$$\begin{eqnarray}\displaystyle {\it\phi}_{eO(2)} & = & \displaystyle P_{O(2)}+2{\it\zeta}_{i}u_{iO(2)}+\left(|{\bf\zeta}|^{2}-\frac{5}{2}\right){\it\tau}_{O(2)}+2{\it\zeta}_{i}{\it\omega}_{O(1)}u_{iO(1)}+\left(|{\bf\zeta}|^{2}-\frac{3}{2}\right){\it\omega}_{O(1)}{\it\tau}_{O(1)}\nonumber\\ \displaystyle & & \displaystyle +\,2\left({\it\zeta}_{i}{\it\zeta}_{j}-\frac{{\it\delta}_{ij}}{2}\right)u_{iO(1)}u_{jO(1)}+2{\it\zeta}_{i}\left(|{\bf\zeta}|^{2}-\frac{5}{2}\right)u_{iO(1)}{\it\tau}_{O(1)}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{1}{2}\left(|{\bf\zeta}|^{4}-5|{\bf\zeta}|^{2}+\frac{15}{4}\right){\it\tau}_{O(1)}^{2}.\end{eqnarray}$$

4.0.2. Boundary conditions for the outer solution at infinity

Up to now, we have not considered the boundary condition for the outer solution ${\it\phi}_{O}$. At infinity, ${\it\phi}_{O}$ is required to satisfy the condition

(4.14)$$\begin{eqnarray}{\it\phi}_{O}\rightarrow {\it\phi}_{\infty }=\frac{1}{{\rm\pi}^{3/2}}\exp (-({\it\zeta}_{i}-{\it\epsilon}{\it\delta}_{i1})^{2})E^{-1}-1\quad \text{as}~{\it\eta}\rightarrow \infty ,\end{eqnarray}$$

where ${\it\eta}=(y_{i}^{2})^{1/2}(={\it\epsilon}r)$ (cf. (2.6)). If we expand ${\it\phi}_{\infty }$ as ${\it\phi}_{\infty }\!={\it\phi}_{\infty (1)}{\it\epsilon}+{\it\phi}_{\infty (2)}{\it\epsilon}^{2}+\cdots \,$, the boundary condition is satisfied by demanding

(4.15)$$\begin{eqnarray}{\it\phi}_{O(m)}\rightarrow {\it\phi}_{\infty (m)}\quad \text{as}~{\it\eta}\rightarrow \infty\end{eqnarray}$$

($m=1,2,\dots$), where

(4.16a−c)$$\begin{eqnarray}{\it\phi}_{\infty (1)}=2{\it\zeta}_{1},\quad {\it\phi}_{\infty (2)}=2{\it\zeta}_{1}^{2}-1,\quad {\it\phi}_{\infty (3)}={\textstyle \frac{4}{3}}{\it\zeta}_{1}\left({\it\zeta}_{1}^{2}-{\textstyle \frac{3}{2}}\right),\quad \text{etc}.\end{eqnarray}$$

On the other hand, since the functional form of ${\it\phi}_{O(m)}$ with respect to the molecular velocity is specified (cf. (4.11a) and (4.11b)), the above condition is satisfied if the macroscopic quantities contained in ${\it\phi}_{O(m)}$ take special values at ${\it\eta}\rightarrow \infty$. For example for ${\it\phi}_{O(1)}$ and ${\it\phi}_{O(2)}$, the conditions take the following forms:

(4.17a−c)$$\begin{eqnarray}P_{O(1)}\rightarrow 0,\quad u_{iO(1)}\rightarrow {\it\delta}_{i1},\quad {\it\tau}_{O(1)}\rightarrow 0\quad \text{as}~{\it\eta}\rightarrow \infty ,\end{eqnarray}$$

(4.18a−c)$$\begin{eqnarray}P_{O(2)}\rightarrow 0,\quad u_{iO(2)}\rightarrow 0,\quad {\it\tau}_{O(2)}\rightarrow 0\quad \text{as}~{\it\eta}\rightarrow \infty .\end{eqnarray}$$

These conditions will provide natural constraints for the fluid-dynamic-type equations at infinity. The remaining boundary conditions (at ${\it\eta}\rightarrow 0$) will be derived in the course of the analysis by considering the matching with the inner solution in § 5.

5. Asymptotic matching

We now come back to the original problem. In the preceding section, we have introduced a slowly varying solution and derived the sets of fluid-dynamic-type equations that describe the overall behaviour of the gas in the region far from the sphere. On the other hand, we can still expect that the solution varies on the scale of the sphere size in the vicinity of the sphere (i.e. the near region). If the slowly varying outer solution ${\it\phi}_{O}$, valid in the far region, and the inner solution ${\it\phi}_{S}$, valid in the near region, can be made to match in the intermediate region, they together form a meaningful solution of the original problem.

Let us recall that the inner solution ${\it\phi}_{S}$ is described by the linearised Boltzmann equation (see e.g. (3.4)). For the linearisation to be valid, the nonlinear term $\mathscr{J}({\it\phi}_{S},{\it\phi}_{S})$ should remain much smaller than the term ${\it\zeta}_{i}\partial {\it\phi}_{S}/\partial x_{i}$. From Proposition 1 in § 3, the magnitude of $\partial {\it\phi}_{S}/\partial x_{i}$ at large $r$ is estimated as $|\partial {\it\phi}_{S}/\partial x_{i}|\sim |{\it\phi}_{S}|/r\sim {\it\epsilon}/r$. Therefore, the range of $r$ where the inner solution is valid is restricted to

(5.1)$$\begin{eqnarray}r\ll \frac{1}{{\it\epsilon}}\end{eqnarray}$$

(see figure 2). On the other hand, since the fluid-dynamic-type equations describing the outer solution have been derived under the condition of slowly varying property ($\partial {\it\phi}_{O}/\partial y_{i}=O({\it\phi}_{O})$), the condition $\partial {\it\phi}_{O}/\partial y_{i}\ll O({\it\phi}_{O}/{\it\epsilon})$ is required for its validity. If we assume that $|{\it\phi}_{O}|$ grows polynomially as ${\it\eta}\rightarrow 0$, the gradient of ${\it\phi}_{O}$ is estimated as $\partial {\it\phi}_{O}/\partial y_{i}=O({\it\phi}_{O}/{\it\eta})$ for ${\it\eta}\ll 1$. In other words, the range of validity for ${\it\phi}_{O}$ is restricted to ${\it\eta}\gg {\it\epsilon}$ or, equivalently,

(5.2)$$\begin{eqnarray}r\gg 1.\end{eqnarray}$$

Therefore, the near and far regions, where the inner and outer solutions are valid respectively, overlap in the intermediate region characterised by

(5.3)$$\begin{eqnarray}1\ll r\ll \frac{1}{{\it\epsilon}},\end{eqnarray}$$

and we may expect that the two solutions can be joined there, as in the case of the classical theory for the Navier–Stokes equation (Van Dyke 1975). In this section, assuming such a structure, we carry out an asymptotic analysis of the original problem. The analysis will be carried out up to the second order in ${\it\epsilon}$.

Figure 2. Schematic illustration of the near and far regions. The near and far regions, where the inner solution ${\it\phi}_{S}$ and the outer solution ${\it\phi}_{O}$ are valid, respectively, overlap in the intermediate region characterised by $1\ll r\ll 1/{\it\epsilon}$.

5.1. Leading-order problem

In view of the fact that the boundary-value problem (3.8)–(3.10) in § 3 has a solution, we shall start with the following inner problem:

(5.4)$$\begin{eqnarray}\displaystyle & \displaystyle {\it\zeta}_{i}\frac{\partial {\it\phi}_{S(1)}}{\partial x_{i}}=\frac{1}{k}\mathscr{L}[{\it\phi}_{S(1)}], & \displaystyle\end{eqnarray}$$

(5.5)$$\begin{eqnarray}\displaystyle & {\it\phi}_{S(1)}=\mathscr{K}[{\it\phi}_{S(1)}]\quad \text{for}~{\it\zeta}_{i}n_{i}>0~\text{at}~r=1, & \displaystyle\end{eqnarray}$$

supplemented by the boundary condition

(5.6)$$\begin{eqnarray}{\it\phi}_{S(1)}\rightarrow 2{\it\zeta}_{1}\quad \text{as}~r\rightarrow \infty .\end{eqnarray}$$

The problem is the same as (3.8)–(3.10) and therefore has a solution.

Assuming that ${\it\phi}_{S(1)}$ has been obtained, we now consider the leading-order outer solution ${\it\phi}_{O(1)}$. First, we note that the macroscopic quantities $P_{O(1)}$, $u_{iO(1)}$ and ${\it\tau}_{O(1)}$ contained in ${\it\phi}_{O(1)}$ (cf. (4.11a) with (4.13a)) are described by the fluid-dynamic-type equations (4.8) and (4.9a–d). The corresponding boundary condition at infinity is given by (4.17). The remaining condition is provided by the connection (or matching) condition for ${\it\phi}_{O(1)}$ with ${\it\phi}_{S(1)}$ as follows.

To derive the connection condition, we consider a point where $r$ is large but ${\it\epsilon}r$ is small. According to Proposition 1, ${\it\phi}_{S(1)}$ behaves like

(5.7)$$\begin{eqnarray}{\it\phi}_{S(1)}=2{\it\zeta}_{1}+O(r^{-1})\end{eqnarray}$$

at such a point. This means that the inner solution ${\it\phi}_{S(1)}$, if written in terms of the outer variable ${\it\eta}(={\it\epsilon}r)$, has the form

(5.8)$$\begin{eqnarray}{\it\phi}_{S(1)}^{\ast }=2{\it\zeta}_{1}+O({\it\epsilon})\end{eqnarray}$$

at the point under consideration. Here, the asterisk has been added to emphasise that the independent variable of the inner solution has been changed from the original $r$ (or $x_{i}$) to ${\it\eta}$ (or $y_{i}$), i.e. ${\it\phi}_{S(1)}^{\ast }={\it\phi}_{S(1)}|_{r={\it\eta}/{\it\epsilon}}$. From (5.8), we immediately see that ${\it\phi}_{O(1)}$ coincides with ${\it\phi}_{S(1)}^{\ast }$ within an error of $O({\it\epsilon})$, if we impose the following connection condition:

(5.9)$$\begin{eqnarray}{\it\phi}_{O(1)}\rightarrow 2{\it\zeta}_{1}\quad \text{as}~{\it\eta}\rightarrow 0.\end{eqnarray}$$

To proceed further, we note that the functional form of ${\it\phi}_{O(1)}$ is specified as given by (4.11a) (with (4.13a)). Therefore, in order for the above condition to be satisfied, the macroscopic quantities contained in ${\it\phi}_{O(1)}$ should take the following values as ${\it\eta}\rightarrow 0$:

(5.10a−c)$$\begin{eqnarray}P_{O(1)}\rightarrow 0,\quad u_{iO(1)}\rightarrow {\it\delta}_{i1},\quad {\it\tau}_{O(1)}\rightarrow 0\quad \text{as}~{\it\eta}\rightarrow 0.\end{eqnarray}$$

In summary, the matching condition (5.10) and the boundary condition at infinity (4.17) provide appropriate boundary conditions for the leading-order fluid-dynamic-type equations (4.8) and (4.9a–d).

It is easily seen that

(5.11)$$\begin{eqnarray}P_{O(1)}={\it\tau}_{O(1)}={\it\omega}_{O(1)}=0,\end{eqnarray}$$

(5.12a,b)$$\begin{eqnarray}u_{iO(1)}={\it\delta}_{i1},\quad P_{O(2)}=\text{const.,}\end{eqnarray}$$

(trivially) satisfy the equations and the boundary condition (4.17) as well as the matching condition (5.10). Consequently, the leading-order outer solution (4.11a) is reduced to

(5.13)$$\begin{eqnarray}{\it\phi}_{O(1)}=2{\it\zeta}_{1}.\end{eqnarray}$$

That is, ${\it\phi}_{O(1)}$ is just a linearised uniform flow.

In summary, the leading-order solution is given by ${\it\phi}_{O(1)}=2{\it\zeta}_{1}$ (far from the sphere) and by ${\it\phi}_{S(1)}$ (in the vicinity of the sphere), where ${\it\phi}_{S(1)}$ solves the problem (5.4)–(5.6). Incidentally, it is worth noting that the constant pressure $P_{O(2)}$ (cf. (5.12)) is left undetermined at this stage. This constant will be determined later when we consider the next-order outer problem ((5.31) below).

5.2. Further transformation of the problem for ${\it\phi}_{S(1)}$

Before going into the next-order approximation, here and in the subsequent subsection we shall recall some basic properties of the leading-order inner problem (i.e. the problem (5.4)–(5.6)). Thanks to the axial symmetry of $\mathscr{L}$ and $\mathscr{K}$, we can employ a similarity solution of the form (Sone & Aoki 1983; Sone 2007)

(5.14)$$\begin{eqnarray}{\it\phi}_{S(1)}={\it\Phi}_{a}^{(1)}(r,{\it\zeta}_{r},|{\bf\zeta}|)\cos {\it\theta}+{\it\zeta}_{{\it\theta}}{\it\Phi}_{b}^{(1)}(r,{\it\zeta}_{r},|{\bf\zeta}|)\sin {\it\theta},\end{eqnarray}$$

where ${\it\Phi}_{a}^{(1)}$ and ${\it\Phi}_{b}^{(1)}$ are functions of $r$, ${\it\zeta}_{r}$ and $|{\bf\zeta}|$. Then, the problem is reduced to a (spatially one-dimensional) boundary-value problem for ${\it\Phi}_{a}^{(1)}$ and ${\it\Phi}_{b}^{(1)}$ as follows:

(5.15)$$\begin{eqnarray}\displaystyle & \displaystyle {\it\zeta}_{r}\frac{\partial {\it\Phi}_{a}^{(1)}}{\partial r}+\frac{|{\bf\zeta}|^{2}-{\it\zeta}_{r}^{2}}{r}\frac{\partial {\it\Phi}_{a}^{(1)}}{\partial {\it\zeta}_{r}}+\frac{|{\bf\zeta}|^{2}-{\it\zeta}_{r}^{2}}{r}{\it\Phi}_{b}^{(1)}=\frac{1}{k}\mathscr{L}[{\it\Phi}_{a}^{(1)}], & \displaystyle\end{eqnarray}$$

(5.16)$$\begin{eqnarray}\displaystyle & \displaystyle {\it\zeta}_{r}\frac{\partial {\it\Phi}_{b}^{(1)}}{\partial r}+\frac{|{\bf\zeta}|^{2}-{\it\zeta}_{r}^{2}}{r}\frac{\partial {\it\Phi}_{b}^{(1)}}{\partial {\it\zeta}_{r}}-\frac{{\it\zeta}_{r}}{r}{\it\Phi}_{b}^{(1)}-\frac{{\it\Phi}_{a}^{(1)}}{r}=\frac{1}{k}\mathscr{L}_{1}[{\it\Phi}_{b}^{(1)}], & \displaystyle\end{eqnarray}$$

(5.17)$$\begin{eqnarray}\displaystyle & \displaystyle \left[\begin{array}{@{}c@{}}{\it\Phi}_{a}^{(1)}\\ {\it\Phi}_{b}^{(1)}\end{array}\right]=\left[\begin{array}{@{}c@{}}\mathscr{K}[{\it\Phi}_{a}^{(1)}]\\ \mathscr{K}_{1}[{\it\Phi}_{b}^{(1)}]\end{array}\right]\quad \text{for}~{\it\zeta}_{r}>0~(r=1), & \displaystyle\end{eqnarray}$$

(5.18)$$\begin{eqnarray}\displaystyle & \displaystyle \left[\begin{array}{@{}c@{}}{\it\Phi}_{a}^{(1)}\\ {\it\Phi}_{b}^{(1)}\end{array}\right]\rightarrow \left[\begin{array}{@{}c@{}}2{\it\zeta}_{r}\\ -2\end{array}\right]\quad \text{as}~r\rightarrow \infty , & \displaystyle\end{eqnarray}$$

where $\mathscr{L}_{1}$ and $\mathscr{K}_{1}$ are operators defined in appendix C. Further, substituting (5.14) into (3.11), we find that the (leading-order) macroscopic variables take the following simple forms:

(5.19a)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{{\it\omega}_{S(1)}}{\cos {\it\theta}}=\widetilde{{\it\omega}}_{S(1)}(r), & \displaystyle\end{eqnarray}$$

(5.19b)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{u_{rS(1)}}{\cos {\it\theta}}=\widetilde{u}_{rS(1)}(r),\quad \frac{u_{{\it\theta}S(1)}}{\sin {\it\theta}}=\widetilde{u}_{{\it\theta}S(1)}(r),\quad u_{{\it\varphi}S(1)}=0, & \displaystyle\end{eqnarray}$$

(5.19c)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{{\it\tau}_{S(1)}}{\cos {\it\theta}}=\widetilde{{\it\tau}}_{S(1)}(r),\quad \frac{P_{S(1)}}{\cos {\it\theta}}=\widetilde{P}_{S(1)}(r), & \displaystyle\end{eqnarray}$$

(5.19d)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x1D617}_{rrS(1)}}{\cos {\it\theta}}=\widetilde{P}_{rrS(1)}(r),\quad \frac{\unicode[STIX]{x1D617}_{r{\it\theta}S(1)}}{\sin {\it\theta}}=\widetilde{P}_{r{\it\theta}S(1)}(r), & \displaystyle\end{eqnarray}$$

(5.19e)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x1D617}_{{\it\theta}{\it\theta}S(1)}}{\cos {\it\theta}}=\widetilde{P}_{{\it\theta}{\it\theta}S(1)}(r),\quad \frac{\unicode[STIX]{x1D617}_{{\it\varphi}{\it\varphi}S(1)}}{\cos {\it\theta}}=\widetilde{P}_{{\it\varphi}{\it\varphi}S(1)}(r),\quad \unicode[STIX]{x1D617}_{r{\it\varphi}S(1)}=\unicode[STIX]{x1D617}_{{\it\theta}{\it\varphi}S(1)}=0, & \displaystyle\end{eqnarray}$$

(5.19f)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{Q_{rS(1)}}{\cos {\it\theta}}=\widetilde{Q}_{rS(1)}(r),\quad \frac{Q_{{\it\theta}S(1)}}{\sin {\it\theta}}=\widetilde{Q}_{{\it\theta}S(1)}(r),\quad Q_{{\it\varphi}S(1)}=0, & \displaystyle\end{eqnarray}$$

where the quantities on the right-hand side, depending only on $r$, are given by

(5.20a)$$\begin{eqnarray}\displaystyle & \displaystyle \widetilde{{\it\omega}}_{S(1)}=2{\rm\pi}\int _{0}^{\infty }\int _{0}^{{\rm\pi}}{\it\zeta}^{2}\sin {\it\theta}_{{\it\zeta}}{\it\Phi}_{a}^{(1)}E\,\text{d}{\it\theta}_{{\it\zeta}}\,\text{d}{\it\zeta}, & \displaystyle\end{eqnarray}$$

(5.20b)$$\begin{eqnarray}\displaystyle & \displaystyle \widetilde{u}_{rS(1)}={\rm\pi}\int _{0}^{\infty }\int _{0}^{{\rm\pi}}{\it\zeta}^{3}\sin 2{\it\theta}_{{\it\zeta}}{\it\Phi}_{a}^{(1)}E\,\text{d}{\it\theta}_{{\it\zeta}}\,\text{d}{\it\zeta}, & \displaystyle\end{eqnarray}$$

(5.20c)$$\begin{eqnarray}\displaystyle & \displaystyle \widetilde{u}_{{\it\theta}S(1)}={\rm\pi}\int _{0}^{\infty }\int _{0}^{{\rm\pi}}{\it\zeta}^{4}\sin ^{3}{\it\theta}_{{\it\zeta}}{\it\Phi}_{b}^{(1)}E\,\text{d}{\it\theta}_{{\it\zeta}}\,\text{d}{\it\zeta}, & \displaystyle\end{eqnarray}$$

(5.20d)$$\begin{eqnarray}\displaystyle & \displaystyle \widetilde{{\it\tau}}_{S(1)}=\frac{4}{3}{\rm\pi}\int _{0}^{\infty }\int _{0}^{{\rm\pi}}{\it\zeta}^{2}\left({\it\zeta}^{2}-\frac{3}{2}\right)\sin {\it\theta}_{{\it\zeta}}{\it\Phi}_{a}^{(1)}E\,\text{d}{\it\theta}_{{\it\zeta}}\,\text{d}{\it\zeta}, & \displaystyle\end{eqnarray}$$

(5.20e)$$\begin{eqnarray}\displaystyle & \displaystyle \widetilde{P}_{S(1)}=\frac{4}{3}{\rm\pi}\int _{0}^{\infty }\int _{0}^{{\rm\pi}}{\it\zeta}^{4}\sin {\it\theta}_{{\it\zeta}}{\it\Phi}_{a}^{(1)}E\,\text{d}{\it\theta}_{{\it\zeta}}\,\text{d}{\it\zeta}=\widetilde{{\it\omega}}_{S(1)}+\widetilde{{\it\tau}}_{S(1)}, & \displaystyle\end{eqnarray}$$

(5.20f)$$\begin{eqnarray}\displaystyle & \displaystyle \widetilde{P}_{rrS(1)}=2{\rm\pi}\int _{0}^{\infty }\int _{0}^{{\rm\pi}}{\it\zeta}^{4}\cos {\it\theta}_{{\it\zeta}}\sin 2{\it\theta}_{{\it\zeta}}{\it\Phi}_{a}^{(1)}E\,\text{d}{\it\theta}_{{\it\zeta}}\,\text{d}{\it\zeta}, & \displaystyle\end{eqnarray}$$

(5.20g)$$\begin{eqnarray}\displaystyle & \displaystyle \widetilde{P}_{r{\it\theta}S(1)}=2{\rm\pi}\int _{0}^{\infty }\int _{0}^{{\rm\pi}}{\it\zeta}^{5}\cos {\it\theta}_{{\it\zeta}}\sin ^{3}{\it\theta}_{{\it\zeta}}{\it\Phi}_{b}^{(1)}E\,\text{d}{\it\theta}_{{\it\zeta}}\,\text{d}{\it\zeta}, & \displaystyle\end{eqnarray}$$

(5.20h)$$\begin{eqnarray}\displaystyle & \displaystyle \widetilde{P}_{{\it\theta}{\it\theta}S(1)}=\widetilde{P}_{{\it\varphi}{\it\varphi}S(1)}=2{\rm\pi}\int _{0}^{\infty }\int _{0}^{{\rm\pi}}{\it\zeta}^{4}\sin ^{3}{\it\theta}_{{\it\zeta}}{\it\Phi}_{a}^{(1)}E\,\text{d}{\it\theta}_{{\it\zeta}}\,\text{d}{\it\zeta}, & \displaystyle\end{eqnarray}$$

(5.20i)$$\begin{eqnarray}\displaystyle & \displaystyle \widetilde{Q}_{rS(1)}={\rm\pi}\int _{0}^{\infty }\int _{0}^{{\rm\pi}}{\it\zeta}^{3}\left({\it\zeta}^{2}-\frac{5}{2}\right)\sin 2{\it\theta}_{{\it\zeta}}{\it\Phi}_{a}^{(1)}E\,\text{d}{\it\theta}_{{\it\zeta}}\,\text{d}{\it\zeta}, & \displaystyle\end{eqnarray}$$

(5.20j)$$\begin{eqnarray}\displaystyle & \displaystyle \widetilde{Q}_{{\it\theta}S(1)}={\rm\pi}\int _{0}^{\infty }\int _{0}^{{\rm\pi}}{\it\zeta}^{4}\left({\it\zeta}^{2}-\frac{5}{2}\right)\sin ^{3}{\it\theta}_{{\it\zeta}}{\it\Phi}_{b}^{(1)}E\,\text{d}{\it\theta}_{{\it\zeta}}\,\text{d}{\it\zeta}. & \displaystyle\end{eqnarray}$$

Here, ${\it\zeta}=|{\bf\zeta}|$, and ${\it\theta}_{{\it\zeta}}$ ($0\leqslant {\it\theta}_{{\it\zeta}}\leqslant {\rm\pi}$) is the angle between ${\bf\zeta}$ and the radial direction, i.e. ${\it\zeta}_{r}=|{\bf\zeta}|\cos {\it\theta}_{{\it\zeta}}$ (see figure 1).

5.3. Asymptotic behaviour of ${\it\phi}_{S(1)}$ at large $r$

Based on the above decomposition, we shall now derive an asymptotic expression for ${\it\phi}_{S(1)}$ at large $r\gg 1$ that will be needed in the subsequent analysis. In this problem, the approach of ${\it\phi}_{S(1)}$ (or (${\it\Phi}_{a}^{(1)},{\it\Phi}_{b}^{(1)}$)) to the (linearised) uniform flow is very slow. If we assume it to be of $r^{-{\it\alpha}}$ (${\it\alpha}>0$), the length scale of variation of ${\it\phi}_{S(1)}$ will be of the order of $r^{\ast }$ for $r>r^{\ast }$, where $r^{\ast }$ is a sufficiently large number. In other words, the local Knudsen number, $k_{loc}=k/r^{\ast }$, is sufficiently small in the region where $r>r_{\ast }$, and, therefore, we can make use of the asymptotic theory of the Boltzmann equation for a slightly rarefied gas (Sone 2002, 2007) to investigate the asymptotic behaviour of ${\it\phi}_{S(1)}$ at large $r$. Because ${\it\phi}_{S(1)}$ is described by the linearised Boltzmann equation, we can apply the linear theory (or the Grad–Hilbert expansion) (Sone 2002). As a result, the flow velocity and the pressure are described by the Stokes equations for an incompressible fluid and the temperature is described by the Laplace equation, i.e. 

(5.21a−d)$$\begin{eqnarray}\frac{\partial u_{iS(1)}}{\partial x_{i}}=0,\quad \frac{\partial P_{S(1)}}{\partial x_{i}}-{\it\gamma}_{1}k\frac{\partial ^{2}u_{iS(1)}}{\partial x_{j}^{2}}=0,\quad \frac{\partial ^{2}{\it\tau}_{S(1)}}{\partial x_{j}^{2}}=0,\quad P_{S(1)}={\it\omega}_{S(1)}+{\it\tau}_{S(1)}.\end{eqnarray}$$

In addition, the velocity distribution function is expressed in terms of $P_{S(1)}$, $u_{iS(1)}$ and ${\it\tau}_{S(1)}$ as follows (Sone 2002):

(5.22)$$\begin{eqnarray}\displaystyle {\it\phi}_{S(1)}\! & = & \displaystyle \!P_{S(1)}+2{\it\zeta}_{i}u_{iS(1)}+\left(|{\bf\zeta}|^{2}-\frac{5}{2}\right){\it\tau}_{S(1)}\nonumber\\ \displaystyle & & \displaystyle -\,k{\it\zeta}_{i}{\it\zeta}_{j}B(|{\bf\zeta}|)\frac{\partial u_{iS(1)}}{\partial x_{j}}-k{\it\zeta}_{i}A(|{\bf\zeta}|)\frac{\partial {\it\tau}_{S(1)}}{\partial x_{i}}+\frac{k}{{\it\gamma}_{1}}{\it\zeta}_{i}D_{1}(|{\bf\zeta}|)\frac{\partial P_{S(1)}}{\partial x_{i}}\nonumber\\ \displaystyle & & \displaystyle +\,k^{2}{\it\zeta}_{i}{\it\zeta}_{j}{\it\zeta}_{k}D_{2}(|{\bf\zeta}|)\frac{\partial ^{2}u_{iS(1)}}{\partial x_{j}\partial x_{k}}-k^{2}{\it\zeta}_{i}{\it\zeta}_{j}F(|{\bf\zeta}|)\frac{\partial ^{2}{\it\tau}_{S(1)}}{\partial x_{i}\partial x_{j}}+\cdots \,,\end{eqnarray}$$

where the functions $A(|{\bf\zeta}|)$, $B(|{\bf\zeta}|),\dots ,F(|{\bf\zeta}|)$ as well as the constant ${\it\gamma}_{1}$ are defined in appendix B.

Since the angular dependence of the macroscopic quantities is explicit in (5.19), its substitution into (5.21) yields a set of ordinary differential equations for $\widetilde{u}_{rS(1)}(r)$, $\widetilde{u}_{{\it\theta}S(1)}(r)$, $\widetilde{P}_{S(1)}(r)$ and $\widetilde{{\it\tau}}_{S(1)}(r)$, which is easily integrated under the condition $\widetilde{u}_{rS(1)}\rightarrow 1$, $\widetilde{u}_{{\it\theta}S(1)}\rightarrow -1$, $\widetilde{P}_{S(1)}\rightarrow 0$, $\widetilde{{\it\tau}}_{S(1)}\rightarrow 0$ at $r\rightarrow \infty$, corresponding to the condition ${\it\phi}_{S(1)}\rightarrow 2{\it\zeta}_{1}$. As the result, we obtain

(5.23a)$$\begin{eqnarray}\displaystyle & \displaystyle \widetilde{u}_{rS(1)}=1+\frac{c_{1}}{r}+\frac{c_{2}}{r^{3}},\quad \widetilde{u}_{{\it\theta}S(1)}=-\left(1+\frac{c_{1}}{2r}-\frac{c_{2}}{2r^{3}}\right), & \displaystyle\end{eqnarray}$$

(5.23b)$$\begin{eqnarray}\displaystyle & \displaystyle \widetilde{{\it\omega}}_{S(1)}=\frac{{\it\gamma}_{1}kc_{1}-c_{3}}{r^{2}},\quad \widetilde{{\it\tau}}_{S(1)}=\frac{c_{3}}{r^{2}},\quad \widetilde{P}_{S(1)}=\frac{{\it\gamma}_{1}kc_{1}}{r^{2}}, & \displaystyle\end{eqnarray}$$

where $c_{1}$, $c_{2}$ and $c_{3}$ are constants introduced because the boundary condition at $r=1$ is not specified. Then, further substitution of (5.19) with (5.23) into (5.22) (written in the spherical coordinates) yields, after taking into account the decomposition (5.14), the following asymptotic expressions for ${\it\Phi}_{a}^{(1)}$ and ${\it\Phi}_{b}^{(1)}$ at large $r$:

(5.24a)$$\begin{eqnarray}\displaystyle {\it\Phi}_{a}^{(1)} & = & \displaystyle 2\left(1+\frac{c_{1}}{r}+\frac{c_{2}}{r^{3}}\right){\it\zeta}_{r}+\frac{c_{3}}{r^{2}}\left(|{\bf\zeta}|^{2}-\frac{5}{2}\right)\nonumber\\ \displaystyle & & \displaystyle +\,\frac{k}{r^{2}}\left[{\it\gamma}_{1}c_{1}+\frac{2c_{3}}{r}{\it\zeta}_{r}A-\frac{1}{2}\left(c_{1}+\frac{3c_{2}}{r^{2}}\right)(|{\bf\zeta}|^{2}-3{\it\zeta}_{r}^{2})B\right]\nonumber\\ \displaystyle & & \displaystyle -\,\frac{k^{2}}{r^{3}}\left\{2c_{1}{\it\zeta}_{r}D_{1}-\frac{3c_{3}}{r}(|{\bf\zeta}|^{2}-3{\it\zeta}_{r}^{2})F\right.\nonumber\\ \displaystyle & & \displaystyle \left.+\,2\left[c_{1}(2|{\bf\zeta}|^{2}-3{\it\zeta}_{r}^{2})+\frac{3c_{2}}{r^{2}}(3|{\bf\zeta}|^{2}-5{\it\zeta}_{r}^{2})\right]{\it\zeta}_{r}D_{2}\right\},\end{eqnarray}$$

(5.24b)$$\begin{eqnarray}\displaystyle {\it\Phi}_{b}^{(1)} & = & \displaystyle -2-\frac{c_{1}}{r}+\frac{c_{2}}{r^{3}}+\frac{k}{r^{3}}\left(c_{3}A+\frac{3c_{2}}{r}{\it\zeta}_{r}B\right)\nonumber\\ \displaystyle & & \displaystyle -\,\frac{k^{2}}{r^{3}}\left\{c_{1}D_{1}+\frac{6c_{3}}{r}{\it\zeta}_{r}F+\frac{1}{2}\left[c_{1}(|{\bf\zeta}|^{2}-3{\it\zeta}_{r}^{2})+\frac{9c_{2}}{r^{2}}(|{\bf\zeta}|^{2}-5{\it\zeta}_{r}^{2})\right]D_{2}\right\},\qquad\end{eqnarray}$$

where the arguments of the functions $A(|{\bf\zeta}|)$, $B(|{\bf\zeta}|)$, etc. have been omitted. It should be noted that the expressions (5.23) and (5.24) contain the constants $c_{1}$, $c_{2}$ and $c_{3}$ that still need to be determined. A way to determine them in a numerical computation will be given in § 7 below. For the moment, we shall assume that the solutions ${\it\Phi}_{a}^{(1)}$ and ${\it\Phi}_{b}^{(1)}$ as well as the constants $c_{1}$, $c_{2}$ and $c_{3}$ have been obtained and proceed to the next-order approximation.

5.4. Second-order outer problem for ${\it\phi}_{O(2)}$

Now we proceed to the second-order approximation in ${\it\epsilon}$. We first consider the outer problem for ${\it\phi}_{O(2)}$ here, and subsequently ${\it\phi}_{S(2)}$ in the following subsection. To begin with, we repeat the fluid-dynamic-type equations describing the overall behaviour of the macroscopic quantities associated with the second-order outer solution. That is, owing to (5.11) and (5.12), the pressure $P_{O(2)}$ is constant and the fluid-dynamic-type equations (4.10a–d) are simplified to

(5.25a)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial u_{iO(2)}}{\partial y_{i}}=0, & \displaystyle\end{eqnarray}$$

(5.25b)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial u_{iO(2)}}{\partial y_{1}}=-\frac{1}{2}\frac{\partial P_{O(3)}}{\partial y_{i}}+\frac{{\it\gamma}_{1}k}{2}\frac{\partial ^{2}u_{iO(2)}}{\partial y_{j}^{2}}, & \displaystyle\end{eqnarray}$$

(5.25c)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\tau}_{O(2)}}{\partial y_{1}}=\frac{{\it\gamma}_{2}k}{2}\frac{\partial ^{2}{\it\tau}_{O(2)}}{\partial y_{j}^{2}}, & \displaystyle\end{eqnarray}$$

(5.25d)$$\begin{eqnarray}\displaystyle & P_{O(2)}={\it\omega}_{O(2)}+{\it\tau}_{O(2)}. & \displaystyle\end{eqnarray}$$

The equations are supplemented by the boundary condition at infinity (4.18), i.e. 

(5.26a−c)$$\begin{eqnarray}P_{O(2)}\rightarrow 0,\quad u_{iO(2)}\rightarrow 0,\quad {\it\tau}_{O(2)}\rightarrow 0\quad \text{as}~{\it\eta}\rightarrow \infty ,\end{eqnarray}$$

as well as by an appropriate matching condition at ${\it\eta}\rightarrow 0$ derived below.

As we have seen in § 5.1, ${\it\phi}_{O(1)}$ coincides with ${\it\phi}_{S(1)}^{\ast }$ if we allow an error of the order of ${\it\epsilon}$ (cf. (5.8) and (5.13)). The leading-order term in the difference can be made precise if we refine ${\it\phi}_{S(1)}^{\ast }$ with the aid of (5.24) and (5.14), that is,

(5.27)$$\begin{eqnarray}{\it\phi}_{S(1)}^{\ast }-{\it\phi}_{O(1)}={\it\epsilon}\frac{c_{1}}{{\it\eta}}(2{\it\zeta}_{r}\cos {\it\theta}-{\it\zeta}_{{\it\theta}}\sin {\it\theta})+O({\it\epsilon}^{2}).\end{eqnarray}$$

Thus, it follows that, in order for ${\it\phi}_{O(1)}+{\it\phi}_{O(2)}{\it\epsilon}$ to coincide with ${\it\phi}_{S(1)}^{\ast }$ within an error of $O({\it\epsilon}^{2})$, the following condition should be imposed:

(5.28)$$\begin{eqnarray}{\it\phi}_{O(2)}\rightarrow \frac{c_{1}}{{\it\eta}}(2{\it\zeta}_{r}\cos {\it\theta}-{\it\zeta}_{{\it\theta}}\sin {\it\theta})\quad \text{as}~{\it\eta}\rightarrow 0.\end{eqnarray}$$

To further derive the corresponding condition for the macroscopic quantities, we recall that the functional form of ${\it\phi}_{O(2)}$ (with respect to ${\bf\zeta}$) is specified and is given by (4.11b), which, with the aid of (5.11) and (5.12), reduces to

(5.29)$$\begin{eqnarray}{\it\phi}_{O(2)}=P_{O(2)}+2{\it\zeta}_{i}u_{iO(2)}+\left(|{\bf\zeta}|^{2}-{\textstyle \frac{5}{2}}\right){\it\tau}_{O(2)}+2{\it\zeta}_{1}^{2}-1.\end{eqnarray}$$

Therefore, in order for (5.28) to be satisfied, the macroscopic quantities contained in (5.29) should take the following values as ${\it\eta}\rightarrow 0$:

(5.30a−d)$$\begin{eqnarray}{\it\eta}u_{rO(2)}\rightarrow c_{1}\cos {\it\theta},\quad {\it\eta}u_{{\it\theta}O(2)}\rightarrow -\frac{c_{1}}{2}\sin {\it\theta},\quad {\it\eta}u_{{\it\varphi}O(2)}\rightarrow 0,\quad {\it\eta}{\it\tau}_{O(2)}\rightarrow 0\quad \text{as}~{\it\eta}\rightarrow 0.\end{eqnarray}$$

Here, we have used the fact that $P_{O(2)}$ and $2{\it\zeta}_{1}^{2}-1$ are independent of ${\it\eta}$. In summary, the second and third conditions of (5.26) as well as (5.30) provide the boundary conditions for the fluid-dynamic-type equations (5.25a–c).

From the first condition of (5.26), $P_{O(2)}(=\text{const}.)$ must vanish identically, i.e.

(5.31)$$\begin{eqnarray}P_{O(2)}=0.\end{eqnarray}$$

On the other hand, (5.25a) and (5.25b) are, respectively, the continuity equation and the Oseen equation for an incompressible fluid, and (5.25c) is an energy equation with a linear convection term. We first observe that ${\it\tau}_{O(2)}=0$ is a trivial solution of (5.25c) satisfying the boundary condition (5.26) (the third condition) as well as the matching condition (5.30) (the fourth condition). Consequently, the density ${\it\omega}_{O(2)}$ also vanishes owing to (5.31) and (5.25d). On the other hand, the axisymmetric solution of the Oseen equation, vanishing at infinity and having the proper singularity (see (5.30)) at the origin, is given by the Oseenlet oriented along the $x_{1}$ axis (Imai 1973), that is,

(5.32a)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{u_{rO(2)}}{{\it\alpha}}=-\frac{1}{{\it\eta}^{2}}+\frac{1}{{\it\eta}^{2}}\left[1+\frac{{\it\eta}}{{\it\gamma}_{1}k}(1+\cos {\it\theta})\right]\exp \left(-\frac{{\it\eta}}{{\it\gamma}_{1}k}(1-\cos {\it\theta})\right), & \displaystyle\end{eqnarray}$$

(5.32b)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{u_{{\it\theta}O(2)}}{{\it\alpha}}=-\frac{1}{{\it\gamma}_{1}k}\frac{\sin {\it\theta}}{{\it\eta}}\exp \left(-\frac{{\it\eta}}{{\it\gamma}_{1}k}(1-\cos {\it\theta})\right), & \displaystyle\end{eqnarray}$$

(5.32c)$$\begin{eqnarray}\displaystyle & u_{{\it\varphi}O(2)}=0, & \displaystyle\end{eqnarray}$$

(5.32d)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{P_{O(3)}}{{\it\alpha}}=\frac{2\cos {\it\theta}}{{\it\eta}^{2}}, & \displaystyle\end{eqnarray}$$

where ${\it\alpha}$ is a constant whose magnitude is of the order of unity, introduced by virtue of the linearity of the equations. Note that $P_{O(3)}$ is determined only up to an additive constant at this point. However, it is not difficult to see that this constant is zero by using the boundary condition for ${\it\phi}_{O(3)}$ at infinity. This fact has already been used in (5.32d).

To determine ${\it\alpha}$, we shall make use of the matching condition (5.30). For this purpose, we expand (5.32a) and (5.32b) in terms of (small) ${\it\eta}$ to obtain

(5.33a)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{u_{rO(2)}}{{\it\alpha}}=\frac{2}{{\it\gamma}_{1}k}\frac{\cos {\it\theta}}{{\it\eta}}-\frac{(1-\cos {\it\theta})(1+3\cos {\it\theta})}{2({\it\gamma}_{1}k)^{2}}+O({\it\eta}), & \displaystyle\end{eqnarray}$$

(5.33b)$$\begin{eqnarray}\displaystyle & \displaystyle \frac{u_{{\it\theta}O(2)}}{{\it\alpha}}=-\frac{1}{{\it\gamma}_{1}k}\frac{\sin {\it\theta}}{{\it\eta}}+\frac{\sin {\it\theta}(1-\cos {\it\theta})}{({\it\gamma}_{1}k)^{2}}+O({\it\eta}). & \displaystyle\end{eqnarray}$$

Multiplying by ${\it\eta}$ and taking the limit ${\it\eta}\rightarrow 0$, we have

(5.34a,b)$$\begin{eqnarray}{\it\eta}u_{rO(2)}\rightarrow \frac{2{\it\alpha}}{{\it\gamma}_{1}k}\cos {\it\theta},\quad {\it\eta}u_{{\it\theta}O(2)}\rightarrow -\frac{{\it\alpha}}{{\it\gamma}_{1}k}\sin {\it\theta}\quad ({\it\eta}\rightarrow 0).\end{eqnarray}$$

Thus, a comparison with the first two conditions in (5.30) immediately yields

(5.35)$$\begin{eqnarray}{\it\alpha}=\frac{{\it\gamma}_{1}kc_{1}}{2}.\end{eqnarray}$$

To summarise, the second-order outer solution has been obtained as

(5.36a)$$\begin{eqnarray}\displaystyle & \displaystyle u_{rO(2)}=\frac{{\it\gamma}_{1}kc_{1}}{2{\it\eta}^{2}}\left[-1+\left(1+\frac{{\it\eta}}{{\it\gamma}_{1}k}(1+\cos {\it\theta})\right)\exp \left(-\frac{{\it\eta}}{{\it\gamma}_{1}k}(1-\cos {\it\theta})\right)\right], & \displaystyle\end{eqnarray}$$

(5.36b)$$\begin{eqnarray}\displaystyle & \displaystyle u_{{\it\theta}O(2)}=-c_{1}\frac{\sin {\it\theta}}{2{\it\eta}}\exp \left(-\frac{{\it\eta}}{{\it\gamma}_{1}k}(1-\cos {\it\theta})\right), & \displaystyle\end{eqnarray}$$

(5.36c)$$\begin{eqnarray}\displaystyle & u_{{\it\varphi}O(2)}={\it\omega}_{O(2)}={\it\tau}_{O(2)}=P_{O(2)}=0, & \displaystyle\end{eqnarray}$$

with the corresponding velocity distribution function (cf. (5.29)) given by

(5.37)$$\begin{eqnarray}{\it\phi}_{O(2)}=2{\it\zeta}_{r}u_{rO(2)}+2{\it\zeta}_{{\it\theta}}u_{{\it\theta}O(2)}+2{\it\zeta}_{1}^{2}-1.\end{eqnarray}$$

Note that ${\it\phi}_{O(2)}$ depends on the constant $c_{1}$ through $u_{rO(2)}$ and $u_{{\it\theta}O(2)}$. It may also be noted that, owing to (5.13) and (5.37), the outer solution is a Maxwellian up to the order ${\it\epsilon}^{2}$, i.e. 

(5.38)$$\begin{eqnarray}(1+{\it\phi}_{O})E=\frac{1}{{\rm\pi}^{3/2}}\text{exp}(-({\it\zeta}_{i}-{\it\epsilon}{\it\delta}_{i1}-{\it\epsilon}^{2}u_{iO(2)})^{2})+O({\it\epsilon}^{3}).\end{eqnarray}$$

5.5. Second-order inner problem for ${\it\phi}_{S(2)}$: connection condition

Finally, we consider the second-order inner problem. The equation and the boundary condition on the sphere are the same as (3.12) and (3.13), i.e.

(5.39)$$\begin{eqnarray}\displaystyle & \displaystyle {\it\zeta}_{i}\frac{\partial {\it\phi}_{S(2)}}{\partial x_{i}}=\frac{1}{k}\mathscr{L}[{\it\phi}_{S(2)}]+\frac{1}{k}\mathscr{J}[{\it\phi}_{S(1)},{\it\phi}_{S(1)}], & \displaystyle\end{eqnarray}$$

(5.40)$$\begin{eqnarray}\displaystyle & \displaystyle {\it\phi}_{S(2)}=\mathscr{K}[{\it\phi}_{S(2)}],\quad {\it\zeta}_{i}n_{i}>0\quad (r=1), & \displaystyle\end{eqnarray}$$

but the boundary condition at infinity (3.14) should be replaced by an appropriate connection condition. Below, we shall derive the connection condition and complete the boundary-value problem.

To derive the connection condition, we consider again a point where $r$ is large but ${\it\epsilon}r$ is small. Owing to (5.33) (with (5.35)) and (5.37), ${\it\phi}_{O(2)}$ has the following asymptotic representation at the point:

(5.41)$$\begin{eqnarray}\displaystyle {\it\phi}_{O(2)} & = & \displaystyle \frac{2c_{1}}{{\it\eta}}\left({\it\zeta}_{r}\cos {\it\theta}-\frac{{\it\zeta}_{{\it\theta}}}{2}\sin {\it\theta}\right)+\frac{c_{1}}{2{\it\gamma}_{1}k}\left[-(1-\cos {\it\theta})(1+3\cos {\it\theta}){\it\zeta}_{r}\right.\nonumber\\ \displaystyle & & \displaystyle \left.+\,2\sin {\it\theta}(1-\cos {\it\theta}){\it\zeta}_{{\it\theta}}\right]+2{\it\zeta}_{1}^{2}-1+O({\it\eta}).\end{eqnarray}$$

Thus, noting the relation ${\it\eta}={\it\epsilon}r$ (see also (5.13)), the outer solution ${\it\phi}_{O}\,(={\it\epsilon}{\it\phi}_{O(1)}+{\it\epsilon}^{2}{\it\phi}_{O(2)}+O({\it\epsilon}^{3}))$ is written in terms of the inner variable $r$ as

(5.42)$$\begin{eqnarray}\displaystyle {\it\phi}_{O}^{\ast } & = & \displaystyle {\it\epsilon}\left[2{\it\zeta}_{1}+\frac{2c_{1}}{r}\left({\it\zeta}_{r}\cos {\it\theta}-\frac{{\it\zeta}_{{\it\theta}}}{2}\sin {\it\theta}\right)\right]+{\it\epsilon}^{2}\left\{\frac{c_{1}}{2{\it\gamma}_{1}k}\left[-(1-\cos {\it\theta})(1+3\cos {\it\theta}){\it\zeta}_{r}\right.\right.\nonumber\\ \displaystyle & & \displaystyle \left.\left.+\,2\sin {\it\theta}(1-\cos {\it\theta}){\it\zeta}_{{\it\theta}}\right]+2{\it\zeta}_{1}^{2}-1\right\}+O({\it\epsilon}^{3}).\end{eqnarray}$$

Here, the asterisk has been attached to indicate that the outer solution is expressed in terms of the inner variable, i.e. ${\it\phi}_{O}^{\ast }={\it\phi}_{O}|_{{\it\eta}={\it\epsilon}r}$. It confirms that the terms of order ${\it\epsilon}$ coincide with ${\it\phi}_{S(1)}$ (at large $r\gg 1$) up to the terms of $r^{-1}$. On the other hand, the terms of ${\it\epsilon}^{2}$ order cannot be represented by ${\it\phi}_{S(1)}$. Therefore, in order for ${\it\phi}_{S}$ to match ${\it\phi}_{O}^{\ast }$, ${\it\phi}_{S(2)}$ should satisfy the following connection condition as $r\rightarrow \infty$:

(5.43)$$\begin{eqnarray}{\it\phi}_{S(2)}\rightarrow \frac{c_{1}}{2{\it\gamma}_{1}k}\left[-(1-\cos {\it\theta})(1+3\cos {\it\theta}){\it\zeta}_{r}+2\sin {\it\theta}(1-\cos {\it\theta}){\it\zeta}_{{\it\theta}}\right]+2{\it\zeta}_{1}^{2}-1\quad \text{as}~r\rightarrow \infty .\end{eqnarray}$$

In summary, the boundary-value problem for ${\it\phi}_{S(2)}$ consists of (5.39) and (5.40), supplemented by the connection condition (5.43). It may be noted that the term $2{\it\zeta}_{1}^{2}-1$ in (5.43) corresponds to the uniform flow (cf. (3.7) or (4.16)) and the rest is the contribution from the slowly varying outer solution.

5.6. Further transformation of the problem for ${\it\phi}_{S(2)}$

With the aid of the relation ${\it\zeta}_{1}={\it\zeta}_{r}\cos {\it\theta}-{\it\zeta}_{{\it\theta}}\sin {\it\theta}$, the matching condition (5.43) is transformed into

(5.44)$$\begin{eqnarray}\displaystyle {\it\phi}_{S(2)} & \rightarrow & \displaystyle -\frac{c_{1}}{{\it\gamma}_{1}k}{\it\zeta}_{1}-\frac{c_{1}}{2{\it\gamma}_{1}k}\left[(1-3\cos ^{2}{\it\theta}){\it\zeta}_{r}+2{\it\zeta}_{{\it\theta}}\cos {\it\theta}\sin {\it\theta}\right]\nonumber\\ \displaystyle & & \displaystyle +\,2({\it\zeta}_{r}\cos {\it\theta}-{\it\zeta}_{{\it\theta}}\sin {\it\theta})^{2}-1\quad \text{as}~r\rightarrow \infty .\end{eqnarray}$$

Then, by the linearity of the problem, we can set the solution in the form

(5.45)$$\begin{eqnarray}{\it\phi}_{S(2)}=-\frac{c_{1}}{2{\it\gamma}_{1}k}{\it\phi}_{S(1)}+{\it\Phi}^{(2)},\end{eqnarray}$$

where the first term satisfies the homogeneous equation and the (homogeneous) boundary condition on the sphere as well as the first term of the matching condition (see also (5.6)), while the second term ${\it\Phi}^{(2)}$ solves the following problem:

(5.46)$$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle {\it\zeta}_{i}\frac{\partial {\it\Phi}^{(2)}}{\partial x_{i}}=\frac{1}{k}\mathscr{L}[{\it\Phi}^{(2)}]+\frac{1}{k}\mathscr{J}[{\it\phi}_{S(1)},{\it\phi}_{S(1)}],\end{eqnarray}$$

(5.47)$$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle {\it\Phi}^{(2)}=\mathscr{K}[{\it\Phi}^{(2)}]\quad \text{for}~{\it\zeta}_{i}n_{i}>0~(r=1),\end{eqnarray}$$

(5.48)$$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle {\it\Phi}^{(2)}\rightarrow -\frac{c_{1}}{2{\it\gamma}_{1}k}\left[(1-3\cos ^{2}{\it\theta}){\it\zeta}_{r}+2{\it\zeta}_{{\it\theta}}\cos {\it\theta}\sin {\it\theta}\right]\nonumber\\ \displaystyle & & \displaystyle \quad +\,2({\it\zeta}_{r}\cos {\it\theta}-{\it\zeta}_{{\it\theta}}\sin {\it\theta})^{2}-1\quad (r\rightarrow \infty ).\end{eqnarray}$$

Further, thanks to the axial symmetry of $\mathscr{L}$, $\mathscr{J}$ and $\mathscr{K}$ (appendix C), we can seek ${\it\Phi}^{(2)}$ in the following form (similarity solution):

(5.49)$$\begin{eqnarray}{\it\Phi}^{(2)}={\it\Phi}_{a}^{(2)}(r,{\it\zeta}_{r},|{\bf\zeta}|)\cos ^{2}{\it\theta}+{\it\zeta}_{{\it\theta}}{\it\Phi}_{b}^{(2)}\cos {\it\theta}\sin {\it\theta}+\frac{{\it\zeta}_{{\it\theta}}^{2}-{\it\zeta}_{{\it\varphi}}^{2}}{2}{\it\Phi}_{c}^{(2)}\sin ^{2}{\it\theta}+{\it\Phi}_{d}^{(2)},\end{eqnarray}$$

where ${\it\Phi}_{J}^{(2)}$ ($J=a,b,c$ or $d$) are functions of $r$, ${\it\zeta}_{r}$ and $|{\bf\zeta}|$. Then, the problem is reduced to a (spatially one-dimensional) boundary-value problem for (${\it\Phi}_{a}^{(2)},{\it\Phi}_{b}^{(2)},{\it\Phi}_{c}^{(2)}$,${\it\Phi}_{d}^{(2)}$) as follows:

(5.50a)$$\begin{eqnarray}\displaystyle & \displaystyle {\it\zeta}_{r}\frac{\partial {\it\Phi}_{a}^{(2)}}{\partial r}+\frac{|{\bf\zeta}|^{2}-{\it\zeta}_{r}^{2}}{r}\frac{\partial {\it\Phi}_{a}^{(2)}}{\partial {\it\zeta}_{r}}+\frac{3}{2}\frac{|{\bf\zeta}|\sin {\it\theta}_{{\it\zeta}}}{r}{\it\Phi}_{b}^{(2)}=\frac{1}{k}\mathscr{L}[{\it\Phi}_{a}^{(2)}]+\frac{1}{k}I_{a}, & \displaystyle\end{eqnarray}$$

(5.50b)$$\begin{eqnarray}\displaystyle & \displaystyle {\it\zeta}_{r}\frac{\partial {\it\Phi}_{b}^{(2)}}{\partial r}+\frac{|{\bf\zeta}|^{2}-{\it\zeta}_{r}^{2}}{r}\frac{\partial {\it\Phi}_{b}^{(2)}}{\partial {\it\zeta}_{r}}-\frac{{\it\zeta}_{r}}{r}{\it\Phi}_{b}^{(2)}-\frac{2}{r}{\it\Phi}_{a}^{(2)}+\frac{|{\bf\zeta}|^{2}-{\it\zeta}_{r}^{2}}{r}{\it\Phi}_{c}^{(2)}=\frac{1}{k}\mathscr{L}_{1}[{\it\Phi}_{b}^{(2)}]+\frac{1}{k}I_{b}, & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(5.50c)$$\begin{eqnarray}\displaystyle & \displaystyle {\it\zeta}_{r}\frac{\partial {\it\Phi}_{c}^{(2)}}{\partial r}+\frac{|{\bf\zeta}|^{2}-{\it\zeta}_{r}^{2}}{r}\frac{\partial {\it\Phi}_{c}^{(2)}}{\partial {\it\zeta}_{r}}-2\frac{{\it\zeta}_{r}}{r}{\it\Phi}_{c}^{(2)}-\frac{{\it\Phi}_{b}^{(2)}}{r}=\frac{1}{k}\mathscr{L}_{2}[{\it\Phi}_{c}^{(2)}]+\frac{1}{k}I_{c}, & \displaystyle\end{eqnarray}$$

(5.50d)$$\begin{eqnarray}\displaystyle & \displaystyle {\it\zeta}_{r}\frac{\partial {\it\Phi}_{d}^{(2)}}{\partial r}+\frac{|{\bf\zeta}|^{2}-{\it\zeta}_{r}^{2}}{r}\frac{\partial {\it\Phi}_{d}^{(2)}}{\partial {\it\zeta}_{r}}-\frac{|{\bf\zeta}|^{2}-{\it\zeta}_{r}^{2}}{2r}{\it\Phi}_{b}^{(2)}=\frac{1}{k}\mathscr{L}[{\it\Phi}_{d}^{(2)}]+\frac{1}{k}I_{d}, & \displaystyle\end{eqnarray}$$

(5.51)$$\begin{eqnarray}\displaystyle \left[\begin{array}{@{}c@{}}{\it\Phi}_{a}^{(2)}\\ {\it\Phi}_{b}^{(2)}\\ {\it\Phi}_{c}^{(2)}\\ {\it\Phi}_{d}^{(2)}\end{array}\right]=\left[\begin{array}{@{}c@{}}\mathscr{K}[{\it\Phi}_{a}^{(2)}]\\ \mathscr{K}_{1}[{\it\Phi}_{b}^{(2)}]\\ \mathscr{K}_{2}[{\it\Phi}_{c}^{(2)}]\\ \mathscr{ K}[{\it\Phi}_{d}^{(2)}]\end{array}\right]\quad \text{for}~{\it\zeta}_{r}>0~(r=1), & & \displaystyle\end{eqnarray}$$

(5.52)$$\begin{eqnarray}\displaystyle \left[\begin{array}{@{}c@{}}{\it\Phi}_{a}^{(2)}\\ {\it\Phi}_{b}^{(2)}\\ {\it\Phi}_{c}^{(2)}\\ {\it\Phi}_{d}^{(2)}\end{array}\right]\rightarrow \left[\begin{array}{@{}c@{}}\displaystyle \frac{3}{2}\frac{c_{1}}{{\it\gamma}_{1}k}{\it\zeta}_{r}+3{\it\zeta}_{r}^{2}-|{\bf\zeta}|^{2}\\ \displaystyle -\frac{c_{1}}{{\it\gamma}_{1}k}-4{\it\zeta}_{r}\\ 2\\ \displaystyle -\frac{c_{1}}{2{\it\gamma}_{1}k}{\it\zeta}_{r}+|{\bf\zeta}|^{2}-{\it\zeta}_{r}^{2}-1\end{array}\right]\quad (r\rightarrow \infty ), & & \displaystyle\end{eqnarray}$$

where the operators $\mathscr{L}_{i}$ and $\mathscr{K}_{i}$ ($i=1$ or 2), as well as $\mathscr{J}_{i}$ ($i=1,2$ or 3) in (5.53) below, are defined in appendix C. The inhomogeneous terms $(I_{a},I_{b},I_{c},I_{d})$ are given by

(5.53a)$$\begin{eqnarray}\displaystyle \displaystyle I_{a} & = & \displaystyle \mathscr{J}[{\it\Phi}_{a}^{(1)},{\it\Phi}_{a}^{(1)}]-\frac{|{\bf\zeta}|^{2}-{\it\zeta}_{r}^{2}}{2}\mathscr{J}_{2}[{\it\Phi}_{b}^{(1)},{\it\Phi}_{b}^{(1)}]-\mathscr{J}_{3}[{\it\Phi}_{b}^{(1)},{\it\Phi}_{b}^{(1)}],\end{eqnarray}$$

(5.53b)$$\begin{eqnarray}\displaystyle I_{b} & = & \displaystyle 2\mathscr{J}_{1}[{\it\Phi}_{a}^{(1)},{\it\Phi}_{b}^{(1)}],\end{eqnarray}$$

(5.53c)$$\begin{eqnarray}\displaystyle I_{c} & = & \displaystyle \mathscr{J}_{2}[{\it\Phi}_{b}^{(1)},{\it\Phi}_{b}^{(1)}],\end{eqnarray}$$

(5.53d)$$\begin{eqnarray}\displaystyle \displaystyle I_{d} & = & \displaystyle \frac{|{\bf\zeta}|^{2}-{\it\zeta}_{r}^{2}}{2}\mathscr{J}_{2}[{\it\Phi}_{b}^{(1)},{\it\Phi}_{b}^{(1)}]+\mathscr{J}_{3}[{\it\Phi}_{b}^{(1)},{\it\Phi}_{b}^{(1)}].\end{eqnarray}$$

Further, if we substitute (5.45) with (5.49) into (3.15), we obtain the following expressions for the (second-order) macroscopic quantities:

(5.54a)$$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle {\it\omega}_{S(2)}=-\frac{c_{1}}{2k}\widetilde{{\it\omega}}_{S(1)}\cos {\it\theta}+\widetilde{{\it\omega}}_{S(2)}^{a}\cos ^{2}{\it\theta}+\widetilde{{\it\omega}}_{S(2)}^{d},\end{eqnarray}$$

(5.54b)$$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle u_{rS(2)}=-\frac{c_{1}}{2k}\widetilde{u}_{rS(1)}\cos {\it\theta}+\left(\widetilde{u}_{rS(2)}^{a}-\widetilde{{\it\omega}}_{S(1)}\widetilde{u}_{rS(1)}\right)\cos ^{2}{\it\theta}+\widetilde{u}_{rS(2)}^{d},\end{eqnarray}$$

(5.54c)$$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle u_{{\it\theta}S(2)}=-\frac{c_{1}}{2k}\widetilde{u}_{{\it\theta}S(1)}\sin {\it\theta}+\left(\widetilde{u}_{{\it\theta}S(2)}^{b}-\widetilde{{\it\omega}}_{S(1)}\widetilde{u}_{{\it\theta}S(1)}\right)\cos {\it\theta}\sin {\it\theta},\end{eqnarray}$$

(5.54d)$$\begin{eqnarray}\displaystyle & & \displaystyle u_{{\it\varphi}S(2)}=0,\end{eqnarray}$$

(5.54e)$$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle {\it\tau}_{S(2)}=-\frac{c_{1}}{2k}\widetilde{{\it\tau}}_{S(1)}\cos {\it\theta}+\left(\widetilde{{\it\tau}}_{S(2)}^{a}-\frac{2}{3}(\widetilde{u}_{rS(1)})^{2}-\widetilde{{\it\omega}}_{S(1)}\widetilde{{\it\tau}}_{S(1)}\right)\cos ^{2}{\it\theta}+\widetilde{{\it\tau}}_{S(2)}^{d}\nonumber\\ \displaystyle & & \displaystyle \qquad \quad -\,\frac{2}{3}(\widetilde{u}_{{\it\theta}S(1)})^{2}\sin ^{2}{\it\theta},\end{eqnarray}$$

(5.54f)$$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle P_{S(2)}=-\frac{c_{1}}{2k}\widetilde{P}_{S(1)}\cos {\it\theta}+\!\left(\widetilde{P}_{S(2)}^{a}-\frac{2}{3}(\widetilde{u}_{rS(1)})^{2}\right)\cos ^{2}{\it\theta}\!+\!\widetilde{P}_{S(2)}^{d}\!-\!\frac{2}{3}(\widetilde{u}_{{\it\theta}S(1)})^{2}\sin ^{2}{\it\theta},\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(5.54g)$$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle \unicode[STIX]{x1D617}_{rrS(2)}=-\frac{c_{1}}{2k}\widetilde{P}_{rrS(1)}\cos {\it\theta}+\left(\widetilde{P}_{rrS(2)}^{a}-2(\widetilde{u}_{rS(1)})^{2}\right)\cos ^{2}{\it\theta}+\widetilde{P}_{rrS(2)}^{d},\end{eqnarray}$$

(5.54h)$$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle \unicode[STIX]{x1D617}_{r{\it\theta}S(2)}=-\frac{c_{1}}{2k}\widetilde{P}_{r{\it\theta}S(1)}\sin {\it\theta}+\left(\widetilde{P}_{r{\it\theta}S(2)}^{b}-2\widetilde{u}_{rS(1)}\widetilde{u}_{{\it\theta}S(1)}\right)\cos {\it\theta}\sin {\it\theta},\end{eqnarray}$$

(5.54i)$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D617}_{r{\it\varphi}S(2)}=0,\end{eqnarray}$$

(5.54j)$$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle \unicode[STIX]{x1D617}_{{\it\theta}{\it\theta}S(2)}=-\frac{c_{1}}{2k}\widetilde{P}_{{\it\theta}{\it\theta}(1)}\cos {\it\theta}+\widetilde{P}_{{\it\theta}{\it\theta}S(2)}^{a}\cos ^{2}{\it\theta}+\left(\frac{1}{4}\widetilde{P}_{{\it\theta}{\it\theta}S(2)}^{c}-2(\widetilde{u}_{{\it\theta}S(1)})^{2}\right)\sin ^{2}{\it\theta}\nonumber\\ \displaystyle & & \displaystyle \qquad \quad +\,\widetilde{P}_{{\it\theta}{\it\theta}S(2)}^{d},\end{eqnarray}$$

(5.54k)$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D617}_{{\it\theta}{\it\varphi}S(2)}=0,\end{eqnarray}$$

(5.54l)$$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle \unicode[STIX]{x1D617}_{{\it\varphi}{\it\varphi}S(2)}=-\frac{c_{1}}{2k}\widetilde{P}_{{\it\varphi}{\it\varphi}S(1)}\cos {\it\theta}+\widetilde{P}_{{\it\varphi}{\it\varphi}S(2)}^{a}\cos ^{2}{\it\theta}-\frac{1}{4}\widetilde{P}_{{\it\varphi}{\it\varphi}S(2)}^{c}\sin ^{2}{\it\theta}+\widetilde{P}_{{\it\varphi}{\it\varphi}S(2)}^{d},\end{eqnarray}$$

(5.54m)$$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle Q_{rS(2)}=-\frac{c_{1}}{2k}\widetilde{Q}_{rS(1)}\cos {\it\theta}+\left(\widetilde{Q}_{rS(2)}^{a}+\frac{5}{2}\widetilde{{\it\omega}}_{S(1)}\widetilde{u}_{rS(1)}-\widetilde{u}_{rS(1)}\widetilde{P}_{rrS(1)}\right.\nonumber\\ \displaystyle & & \displaystyle \qquad \quad \left.-\,\frac{3}{2}\widetilde{P}_{S(1)}\widetilde{u}_{rS(1)}\right)\cos ^{2}{\it\theta}+\widetilde{Q}_{rS(2)}^{d}-\widetilde{u}_{{\it\theta}S(1)}\widetilde{P}_{r{\it\theta}S(1)}\sin ^{2}{\it\theta},\end{eqnarray}$$

(5.54n)$$\begin{eqnarray}\displaystyle & & \displaystyle Q_{{\it\theta}S(2)}=-\frac{c_{1}}{2k}\widetilde{Q}_{{\it\theta}S(1)}\sin {\it\theta}+\left(\widetilde{Q}_{{\it\theta}S(2)}^{b}+\frac{5}{2}\widetilde{{\it\omega}}_{S(1)}\widetilde{u}_{{\it\theta}S(1)}-\widetilde{u}_{rS(1)}\widetilde{P}_{r{\it\theta}S(1)}\right.\nonumber\\ \displaystyle & & \displaystyle \qquad \quad \left.-\,\widetilde{u}_{{\it\theta}S(1)}\widetilde{P}_{{\it\theta}{\it\theta}(1)}^{a}-\frac{3}{2}\widetilde{P}_{S(1)}\widetilde{u}_{{\it\theta}S(1)}\right)\cos {\it\theta}\sin {\it\theta},\end{eqnarray}$$

(5.54o)$$\begin{eqnarray}\displaystyle & & \displaystyle Q_{{\it\varphi}S(2)}=0,\end{eqnarray}$$

where the quantities with tilde are functions of $r$. The explicit forms of these quantities are summarised as ($J=a$ or $d$)

(5.55a)$$\begin{eqnarray}\displaystyle & \displaystyle \widetilde{{\it\omega}}_{S(2)}^{J}=2{\rm\pi}\int _{0}^{\infty }\int _{0}^{{\rm\pi}}{\it\zeta}^{2}\sin {\it\theta}_{{\it\zeta}}{\it\Phi}_{J}^{(2)}E\,\text{d}{\it\theta}_{{\it\zeta}}\,\text{d}{\it\zeta}, & \displaystyle\end{eqnarray}$$

(5.55b)$$\begin{eqnarray}\displaystyle & \displaystyle \widetilde{u}_{rS(2)}^{J}={\rm\pi}\int _{0}^{\infty }\int _{0}^{{\rm\pi}}{\it\zeta}^{3}\sin 2{\it\theta}_{{\it\zeta}}{\it\Phi}_{J}^{(2)}E\,\text{d}{\it\theta}_{{\it\zeta}}\,\text{d}{\it\zeta}, & \displaystyle\end{eqnarray}$$

(5.55c)$$\begin{eqnarray}\displaystyle & \displaystyle \widetilde{u}_{{\it\theta}S(2)}^{b}={\rm\pi}\int _{0}^{\infty }\int _{0}^{{\rm\pi}}{\it\zeta}^{4}\sin ^{3}{\it\theta}_{{\it\zeta}}{\it\Phi}_{b}^{(2)}E\,\text{d}{\it\theta}_{{\it\zeta}}\,\text{d}{\it\zeta}, & \displaystyle\end{eqnarray}$$

(5.55d)$$\begin{eqnarray}\displaystyle & \displaystyle \widetilde{{\it\tau}}_{S(2)}^{J}=\frac{4}{3}{\rm\pi}\int _{0}^{\infty }\int _{0}^{{\rm\pi}}{\it\zeta}^{2}\left({\it\zeta}^{2}-\frac{3}{2}\right)\sin {\it\theta}_{{\it\zeta}}{\it\Phi}_{J}^{(2)}E\,\text{d}{\it\theta}_{{\it\zeta}}\,\text{d}{\it\zeta}, & \displaystyle\end{eqnarray}$$

(5.55e)$$\begin{eqnarray}\displaystyle & \displaystyle \widetilde{P}_{S(2)}^{J}=\frac{4}{3}{\rm\pi}\int _{0}^{\infty }\int _{0}^{{\rm\pi}}{\it\zeta}^{4}\sin {\it\theta}_{{\it\zeta}}{\it\Phi}_{J}^{(2)}E\,\text{d}{\it\theta}_{{\it\zeta}}\,\text{d}{\it\zeta}=\widetilde{{\it\omega}}_{S(2)}^{J}+\widetilde{{\it\tau}}_{S(2)}^{J}, & \displaystyle\end{eqnarray}$$

(5.55f)$$\begin{eqnarray}\displaystyle & \displaystyle \widetilde{P}_{rrS(2)}^{J}=2{\rm\pi}\int _{0}^{\infty }\int _{0}^{{\rm\pi}}{\it\zeta}^{4}\cos {\it\theta}_{{\it\zeta}}\sin 2{\it\theta}_{{\it\zeta}}{\it\Phi}_{J}^{(2)}E\,\text{d}{\it\theta}_{{\it\zeta}}\,\text{d}{\it\zeta}, & \displaystyle\end{eqnarray}$$

(5.55g)$$\begin{eqnarray}\displaystyle & \displaystyle \widetilde{P}_{r{\it\theta}S(2)}^{b}=2{\rm\pi}\int _{0}^{\infty }\int _{0}^{{\rm\pi}}{\it\zeta}^{5}\cos {\it\theta}_{{\it\zeta}}\sin ^{3}{\it\theta}_{{\it\zeta}}{\it\Phi}_{b}^{(2)}E\,\text{d}{\it\theta}_{{\it\zeta}}\,\text{d}{\it\zeta}, & \displaystyle\end{eqnarray}$$

(5.55h)$$\begin{eqnarray}\displaystyle & \left[\begin{array}{@{}c@{}}\widetilde{P}_{{\it\theta}{\it\theta}S(2)}^{J}\\ \widetilde{P}_{{\it\theta}{\it\theta}S(2)}^{c}\end{array}\right]=\left[\begin{array}{@{}c@{}}\widetilde{P}_{{\it\varphi}{\it\varphi}S(2)}^{J}\\ \widetilde{P}_{{\it\varphi}{\it\varphi}S(2)}^{c}\end{array}\right]=2{\rm\pi}\displaystyle \int _{0}^{\infty }\displaystyle \int _{0}^{{\rm\pi}}{\it\zeta}^{4}\sin ^{3}{\it\theta}_{{\it\zeta}}\left[\begin{array}{@{}c@{}}{\it\Phi}_{J}^{(2)}\\ ({\it\zeta}\sin {\it\theta}_{{\it\zeta}})^{2}{\it\Phi}_{c}^{(2)}\end{array}\right]E\,\text{d}{\it\theta}_{{\it\zeta}}\,\text{d}{\it\zeta}, & \displaystyle\end{eqnarray}$$

(5.55i)$$\begin{eqnarray}\displaystyle & \displaystyle \widetilde{Q}_{rS(2)}^{J}={\rm\pi}\int _{0}^{\infty }\int _{0}^{{\rm\pi}}{\it\zeta}^{3}\left({\it\zeta}^{2}-\frac{5}{2}\right)\sin 2{\it\theta}_{{\it\zeta}}{\it\Phi}_{J}^{(2)}E\,\text{d}{\it\theta}_{{\it\zeta}}\,\text{d}{\it\zeta}, & \displaystyle\end{eqnarray}$$

(5.55j)$$\begin{eqnarray}\displaystyle & \displaystyle \widetilde{Q}_{{\it\theta}S(2)}^{b}={\rm\pi}\int _{0}^{\infty }\int _{0}^{{\rm\pi}}{\it\zeta}^{4}\left({\it\zeta}^{2}-\frac{5}{2}\right)\sin ^{3}{\it\theta}_{{\it\zeta}}{\it\Phi}_{b}^{(2)}E\,\text{d}{\it\theta}_{{\it\zeta}}\,\text{d}{\it\zeta}, & \displaystyle\end{eqnarray}$$

where ${\it\zeta}=|{\bf\zeta}|$ and ${\it\theta}_{{\it\zeta}}=\cos ^{-1}({\it\zeta}_{r}/|{\bf\zeta}|)$ (see figure 1).

In summary, to obtain the flow field in the vicinity of the sphere to the second order in ${\it\epsilon}$, it suffices to solve one more spatially one-dimensional problem, in addition to the leading-order linearised problem. More precisely, we solve (5.50)–(5.52) to obtain ${\it\Phi}_{J}^{(2)}$ ($J=a,b,c,d$) and compute $\widetilde{{\it\omega}}_{S(2)}^{J}$, $\widetilde{u}_{rS(2)}^{J}$, $\widetilde{u}_{{\it\theta}S(2)}^{b}$, etc. from (5.55). Then, the second-order macroscopic quantities are readily obtained from (5.54). In this paper, however, we do not carry out this analysis, leaving it to a subsequent study. Instead, we derive the expression for the drag to the second order in ${\it\epsilon}$ in the remaining part, by assuming the existence of the solution to the problem.

6. Drag exerted on the sphere

Let us denote by $(F_{D},0,0)$ the total force acting on the sphere and by $\hat{F}_{D}=F_{D}/p_{\infty }L^{2}$ its dimensionless counterpart. Then, $\hat{F}_{D}$ is given by

(6.1)$$\begin{eqnarray}\hat{F}_{D}=-\int \unicode[STIX]{x1D617}_{1jS}n_{j}\,\text{d}S=-\int (\unicode[STIX]{x1D617}_{rrS}\cos {\it\theta}-\unicode[STIX]{x1D617}_{r{\it\theta}S}\sin {\it\theta})\,\text{d}S,\end{eqnarray}$$

where $\text{d}S=\sin {\it\theta}\,\text{d}{\it\theta}\,\text{d}{\it\varphi}$ is the (dimensionless) surface element on the sphere and the integration is carried out over the whole surface (i.e. $0\leqslant {\it\theta}\leqslant {\rm\pi}$ and $0\leqslant {\it\varphi}<2{\rm\pi}$). Also, the macroscopic quantities are evaluated on the sphere ($r=1$) in this section. Substituting the expansion $\unicode[STIX]{x1D617}_{ijS}=\unicode[STIX]{x1D617}_{ijS(1)}{\it\epsilon}+\cdots \,$ of the inner solution, we obtain the following expansion of the force in terms of ${\it\epsilon}$, i.e.

(6.2)$$\begin{eqnarray}\hat{F}_{D}=\hat{F}_{D(1)}{\it\epsilon}+\hat{F}_{D(2)}{\it\epsilon}^{2}+\cdots \,,\end{eqnarray}$$

where

(6.3)$$\begin{eqnarray}\hat{F}_{D(m)}=-\int (\unicode[STIX]{x1D617}_{rrS(m)}\cos {\it\theta}-\unicode[STIX]{x1D617}_{r{\it\theta}S(m)}\sin {\it\theta})\,\text{d}S\end{eqnarray}$$

($m=1,2,\!\dots$). With the aid of (5.19d), the leading-order drag $\hat{F}_{D(1)}$ is readily calculated to give

(6.4)$$\begin{eqnarray}\hat{F}_{D(1)}=-{\textstyle \frac{4}{3}}{\rm\pi}(\widetilde{P}_{rrS(1)}-2\widetilde{P}_{r{\it\theta}S(1)}),\end{eqnarray}$$

where $\widetilde{P}_{rrS(1)}$ and $\widetilde{P}_{r{\it\theta}S(1)}$ are given by (5.20f) and (5.20g), respectively. This is nothing but the drag investigated in the existing literature on the basis of the linearised Boltzmann equation (or linearised kinetic models). Since $\hat{F}_{D(1)}$ depends on $k$, we denote it by $\hat{F}_{D(1)}=h_{D}(k)$.

Next, with the aid of (5.54g) and (5.54h), the second-order drag is calculated as

(6.5)$$\begin{eqnarray}\displaystyle \hat{F}_{D(2)} & = & \displaystyle \frac{2{\rm\pi}}{3}\frac{c_{1}}{{\it\gamma}_{1}k}(\widetilde{P}_{rrS(1)}-2\widetilde{P}_{r{\it\theta}S(1)})\nonumber\\ \displaystyle & = & \displaystyle -\frac{c_{1}}{2{\it\gamma}_{1}k}\hat{F}_{D(1)}.\end{eqnarray}$$

In this calculation, the contributions of the terms of $(\cos ^{2}{\it\theta},1)$ in $\unicode[STIX]{x1D617}_{rrS(2)}$ and those of $\cos {\it\theta}\sin {\it\theta}$ in $\unicode[STIX]{x1D617}_{r{\it\theta}S(2)}$ vanish by the integration with respect to ${\it\theta}$, and (6.4) has been used in the last equality. Thus, the second-order drag is a constant multiple of the first-order drag. Summarising the above results, we obtain the following expression for the drag:

(6.6)$$\begin{eqnarray}\hat{F}_{D}={\it\epsilon}\left(1-\frac{c_{1}(k)}{2{\it\gamma}_{1}k}{\it\epsilon}\right)h_{D}(k),\end{eqnarray}$$

where the terms smaller than $O({\it\epsilon}^{2})$ have been omitted. It should be recalled that the constant $c_{1}$ depends on $k$. This dependence has been made explicit in the above expression. It should also be noted that $h_{D}(k)$ and $c_{1}(k)$ also depend on the molecular model as well as the gas–surface interaction rule through $\mathscr{L}$ and $\mathscr{K}$.

7. Summary of the drag and numerical results

If we replace ${\it\epsilon}$ by the original $\hat{v}_{\infty }$, the drag obtained in the preceding section is written as

(7.1)$$\begin{eqnarray}\hat{F}_{D}=\hat{v}_{\infty }\left(1-\frac{c_{1}(k)}{2{\it\gamma}_{1}k}\hat{v}_{\infty }\right)h_{D}(k),\end{eqnarray}$$

or, using dimensional quantities,

(7.2)$$\begin{eqnarray}F_{D}=6{\rm\pi}{\it\mu}v_{\infty }L\left(1-\frac{c_{1}(k)}{4}\frac{{\it\rho}_{\infty }v_{\infty }L}{{\it\mu}}\right)H_{D}(k),\end{eqnarray}$$

where ${\it\mu}=(\sqrt{{\rm\pi}}/2){\it\gamma}_{1}p_{\infty }(2RT_{\infty })^{-1/2}\ell _{\infty }$ is the viscosity of the gas at the reference equilibrium state (see appendix B), and the function $H_{D}(k)$ is defined by

(7.3)$$\begin{eqnarray}H_{D}(k)=\frac{h_{D}(k)}{6{\rm\pi}{\it\gamma}_{1}k}.\end{eqnarray}$$

It should be emphasised that the above second-order drag contains only two fundamental functions $h_{D}(k)$ (or $H_{D}(k)$) and $c_{1}(k)$ depending on $k$, which are calculated merely from the information on the leading-order problem for ${\it\phi}_{S(1)}$. This greatly reduces the amount of numerical computation required to obtain actual values of the drag.

In order to obtain the actual values of the drag, we need numerical values of $h_{D}(k)$ and $c_{1}(k)$ for given $\mathscr{L}$ and $\mathscr{K}$ (i.e. the inter-molecular and molecular–surface interactions). The values of $h_{D}$ for a hard-sphere gas under the diffuse reflection condition were obtained by Takata et al. (1993) by an accurate numerical analysis of the problem (5.15)–(5.18) (table V in the reference; $k$ is denoted by $k_{\infty }$). The drag for the linearised kinetic models has also been investigated by many authors (e.g. Cercignani et al. 1968; Lea & Loyalka 1982; Law & Loyalka 1986; Beresnev et al. 1990; Loyalka 1992; Takata et al. 1993). Therefore, we can use these data for $h_{D}$. On the other hand, the value of $c_{1}$ is not available at present. In this paper, we present some results for $c_{1}$ on the basis of the Bhatnagar–Gross–Krook (BGK) model of the Boltzmann equation (Bhatnagar, Gross & Krook 1954; Welander 1954) with the diffuse reflection condition for the gas–surface interaction.

Table 1. $h_{D}$ and $c_{1}$ as a function of $k$ for the BGK model under the diffuse reflection condition.

To obtain $c_{1}$, we need to solve the boundary-value problem (5.15)–(5.18). The numerical method used in the present study is essentially the same as that used in Takata et al. (1993). That is, we employ a finite difference method which is capable of capturing the discontinuity in the velocity distribution function (Sone & Takata 1992). A method to determine $c_{1}$, $c_{2}$ and $c_{3}$ has also been introduced there, where the asymptotic form (5.24) was first derived and used for the purpose of reducing the size of the computational domain (for $r$ variable). To be more specific, let us denote by $r_{D}$ the upper limit of $r$ used in the numerical computation, which is carefully chosen depending on $k$. Then, we impose (5.24) at $r=r_{D}$ as the boundary condition for the incoming molecules with the molecular velocity ${\it\zeta}_{r}<0$, in place of (5.18). In this process, we need the numerical values of $c_{1}$, $c_{2}$ and $c_{3}$ to evaluate (5.24). They are obtained by connecting the numerical values of $\widetilde{u}_{rS(1)}$, $\widetilde{u}_{{\it\theta}S(1)}$ and $\widetilde{{\it\tau}}_{S(1)}$ (cf. (5.20)) with those of (5.23) at $r=r_{D}$. Our numerical results show that $c_{1}$, $c_{2}$ and $c_{3}$ converge to certain values, together with ${\it\Phi}_{(1)}^{a}$ and ${\it\Phi}_{(1)}^{b}$. We take the value of $c_{1}$ thus obtained as its approximate value. The reader is referred to Takata et al. (1993) for further details of the numerical method.

The numerical values of $h_{D}$ and $c_{1}$ obtained in the present computations are shown in table 1 and in figure 3. The values of $h_{D}$ have been recomputed in this study with higher accuracy. $h_{D}$ increases monotonically in $k$. On the other hand, $c_{1}$ takes negative values and its magnitude decreases with $k$.  With the aid of the asymptotic solution for small $k$ (Sone & Aoki 1977a; Aoki & Sone 1987; Sone 2007), one can obtain the asymptotic expressions for $h_{D}$ and $c_{1}$ for small $k$ as follows:

(7.4)$$\begin{eqnarray}\displaystyle & h_{D}=6{\rm\pi}{\it\gamma}_{1}k(1+k_{0}k+(3k_{0}^{2}-4a_{1}+a_{2}+a_{3}+2b_{1})k^{2}+\cdots \,), & \displaystyle\end{eqnarray}$$

(7.5)$$\begin{eqnarray}\displaystyle & \displaystyle c_{1}=-{\textstyle \frac{3}{2}}(1+k_{0}k+(3k_{0}^{2}-4a_{1}+a_{2}+a_{3}+2b_{1})k^{2}+\cdots \,), & \displaystyle\end{eqnarray}$$

where $k_{0}$, $a_{i}$ ($i=1,2$ or 3) and $b_{1}$ are the slip or jump coefficients. The numerical values of these coefficients are tabulated in table 2 for the BGK model as well as for a hard-sphere gas. The $h_{D}$ and $c_{1}$ values based on the formulae, as well as the limiting value $h_{D}\rightarrow 2\sqrt{{\rm\pi}}({\rm\pi}+8)/3$ as $k\rightarrow \infty$, are also shown in figure 3.

Figure 3. (a) $h_{D}$ versus $k$ and (b) $c_{1}$ versus $k$ for the BGK model under the diffuse reflection condition. The symbols represent the numerical results. In (a), the solid (or dashed) line shows (7.4) up to the order $k$ (or $k^{3}$). The dash-dotted line indicates the limiting value $h_{D}\rightarrow 2\sqrt{{\rm\pi}}({\rm\pi}+8)/3$ as $k\rightarrow \infty$. In (b), the solid (or dashed) line shows (7.5) up to the order $k^{0}$ (or $k^{2}$).

Table 2. Slip and jump coefficients for the BGK model and for a hard-sphere gas (HS) under the diffuse reflection condition. The data are taken from Sone (2007) and Takata & Hattori (2012a) (see also Takata & Hattori 2012b).

With the information on $h_{D}$ and $c_{1}$, the force exerted on the sphere is readily calculated from (7.1). The result is shown in figure 4 as a function of $\hat{v}_{\infty }$ for $k=0.1$, 0.5, 1 and 10. Note that the curves are also shown for relatively large $\hat{v}_{\infty }$. In the figure, the solid line shows the drag based on (7.1), whereas the dashed line the linear drag (the first term of (7.1)). The nonlinear effect becomes noticeable with the increase of $\hat{v}_{\infty }$. For $k=10$, there is almost no difference between the present result (solid line) and the linear drag (dashed line). The difference becomes larger as $k$ becomes smaller ($k=1\rightarrow 0.5\rightarrow 0.1$). However, it should be kept in mind that the present analysis assumes that the Mach number is small while the Knudsen number is of the order of unity. In other words, for the formula to be applicable, the Mach number (or $\hat{v}_{\infty }$) should be much smaller than the Knudsen number (or $k$), so that the correction term in (7.1) remains smaller than the leading-order term. This is the limitation of the present formula, when applied to the case of small Knudsen numbers. Incidentally, it may be noted that the formula (7.2) formally coincides with the Navier–Stokes result (Van Dyke 1975) if we take the limit $k\rightarrow 0$ keeping $\hat{v}_{\infty }$ fixed (see (7.4) and (7.5)). However, as mentioned, such a limit cannot be taken without violating the underlying assumption. We need further investigations to clarify the (nonlinear) force in the case of small Knudsen numbers.

Figure 4. Dimensionless force $\hat{F}_{D}$ exerted on the sphere as a function of $\hat{v}_{\infty }$ for various $k$ (the BGK model under the diffuse reflection condition). (a) $k=10$, (b) $k=1$, (c) $k=0.5$, (d) $k=0.1$. The solid line shows (7.1) and the dashed line (7.1) without the second term.

Finally, we compare the present result with some existing numerical and experimental results available in the literature. The result of the comparison is summarised in figure 5 using the Mach number ($\mathit{Ma}$) and the Knudsen number ($\mathit{Kn}$) rather than $\hat{v}_{\infty }$ and $k$. Here, the (dimensionless) drag based on (7.1) and the linear drag based on the first term of (7.1) are shown (as a function of the Mach number) by the solid and dashed lines, respectively, for $\mathit{Kn}=0.1$, 0.2, 0.4, 0.6, 1, 2, 4 and 6. The symbol $\circ$ represents the numerical result obtained by Volkov (2011) by the direct simulation Monte Carlo (DSMC) method for a hard-sphere gas under the diffuse reflection condition. Strictly speaking, these DSMC results are not for a resting sphere, but for a rotating sphere, spinning around an axis of revolution perpendicular to the uniform stream with a constant (dimensionless) angular velocity ${\it\Omega}L/(2RT_{\infty })^{1/2}=0.1$. However, the effect of rotation on the drag is weak when the rotation speed is small (A. N. Volkov, private communication). Therefore, we may consider the DSMC result as approximately the drag for a non-rotating sphere and expect that the comparison provides reasonably accurate information for the present purpose of ‘rough’ comparison. It is also noted that the DSMC result assumes a hard-sphere gas, whereas the present result (more precisely, the numerical values of $h_{D}(k)$ and $c_{1}(k)$) has been obtained on the basis of the BGK model. In the present comparison, we take the viscosity ${\it\mu}=(\sqrt{{\rm\pi}}/2){\it\gamma}_{1}p_{\infty }(2RT_{\infty })^{-1/2}\ell _{\infty }$ (see appendix B) as a common basic quantity, instead of $\ell _{\infty }$, and convert $k$ (or $\mathit{Kn}$) for the BGK model to that for the hard-sphere model using the relation

(7.6)$$\begin{eqnarray}k(\text{BGK})=1.270042427k(\text{hard sphere}).\end{eqnarray}$$

In the same figure, the drag based on the semi-empirical formula proposed by Henderson (1976) (subsonic formula) is also included and shown by the dash-dotted line.

Figure 5. Comparison of the dimensionless force $\hat{F}_{D}$ exerted on the sphere versus $\mathit{Ma}$. (a) $\mathit{Kn}=6$, (b) $\mathit{Kn}=4$, (c) $\mathit{Kn}=2$, (d) $\mathit{Kn}=1$, (e) $\mathit{Kn}=0.6$, (f) $\mathit{Kn}=0.4$, (g) $\mathit{Kn}=0.2$, (h) $\mathit{Kn}=0.1$. The solid line indicates the result based on (7.1) and the dashed line that based on (7.1) without the second term. The symbols ($\circ$) indicate the DSMC result by Volkov (2011) for a rotating sphere, assuming a hard-sphere gas and the diffuse reflection condition. The dash-dotted line represents the drag based on the semi-empirical formula by Henderson (1976).

The agreement between the present and the DSMC results is favourable. For large and moderate Knudsen numbers (figure 5a–f), the present formula (7.1) gives a result close to the numerical result even at a relatively large Mach number (i.e. $\mathit{Ma}=0.6$). For small Knudsen numbers (figure 5g,h), the agreement at $\mathit{Ma}=0.6$ becomes poorer. However, as discussed previously, the present formula is theoretically applicable only to the case where $\mathit{Ma}$ is much smaller than $\mathit{Kn}$. Therefore, this discrepancy does not conflict with the present theory. As a matter of fact, the agreement becomes satisfactory as $\mathit{Ma}$ becomes small enough (say, $\mathit{Ma}\lesssim \mathit{Kn}$) even at these small Knudsen numbers. For Henderson’s semi-empirical formula, a good agreement is observed at moderate Knudsen numbers (see e.g. figure 5f).

8. Conclusions

In this paper, we considered the slow motion of a rarefied gas past a sphere with a uniform temperature and investigated the steady behaviour of the gas on the basis of kinetic theory (the Boltzmann equation). The problem has been the subject of many studies in the past under the assumption that the linearisation of the system is permissible under the condition of small Mach numbers. The present study aimed to extend these works to include a weak nonlinear effect, still keeping the Mach number to be a small quantity.

The primary length scale of the problem is the size of the sphere (i.e. $L$). However, it was shown that, when the Mach number $\mathit{Ma}$ is small, another length scale appears which is much larger than $L$ and is characterised by $L/\mathit{Ma}$, in the far region. The solution with this length scale of variation is described by the fluid-dynamic-type equations derived from the Boltzmann equation (§ 4). Then, the slowly varying solution described by the fluid-dynamic-type equations and the solution in the near region described by the linearised Boltzmann equation were looked for in the form of a power series of Mach number, in such a way that they are joined in the intermediate overlapping region (§ 5). As a result, a drag formula has been derived up to the second order of the Mach number (§ 6). It contains two fundamental functions, $h_{D}(k)$ and $c_{1}(k)$, that are obtained from the information on the solution of the linearised problem, a well-known problem in the literature. To be more specific, $h_{D}$ corresponds to the drag obtained from the linearised problem and $c_{1}$ is related to the asymptotic behaviour of the solution at a large distance from the sphere. Finally, the numerical values of these functions, as well as the drag, were obtained on the basis of the BGK model of the Boltzmann equation under the diffuse reflection boundary condition in § 7.