## 1 Introduction

Thermophoresis refers to the force and, potentially, motion experienced by solid particles or surfaces exposed to gas under rarefied conditions and in the presence of a temperature gradient. This phenomenon seems to have been first noted by Tyndall (Reference Tyndall1870) while observing the spatial redistribution of ambient dust in the proximity of a heated surface. Commonly, the thermophoretic force has been assumed to point from the hot to the cold region, that is, opposite to the temperature gradient. This condition has been labelled as ‘positive’ thermophoresis. On the other hand, a reversal of the force direction, known as ‘negative’ thermophoresis, has also been predicted and – although rarely – observed, as discussed below. Thermophoresis belongs to a class of phenomena promoted in a gas away from thermodynamic equilibrium, which occurs when the collisions between gas molecules are insufficient. Among the various effects in this class are the velocity slip and temperature jump at solid walls and liquid–gas interfaces, transpiration flow, thermal stress, Knudsen layers and heat flux without temperature gradients (Sone Reference Sone2007; Struchtrup & Torrilhon Reference Struchtrup and Torrilhon2008). Review papers (Talbot *et al.*
Reference Talbot, Cheng, Schefer and Willis1980; Bakanov Reference Bakanov1991; Zheng Reference Zheng2002) and book sections (Sone Reference Sone2007) have been written on the subject of thermophoresis. In particular, Young (Reference Young2011) presented a comprehensive examination of the various theories on particle thermophoresis at arbitrary Knudsen numbers under the light of experimental data, and concluded that the accuracy of the measurements and the interval of Knudsen numbers explored are such that confirmation of the validity of the theories could not be achieved.

The main measure of departure from local thermodynamic equilibrium in a gas is the Knudsen number, defined as the ratio of the mean free path between collisions and a characteristic macroscopic length scale. Typically, it is accepted that for Knudsen numbers below 0.01 the Navier–Stokes equations and the Fourier law for conductive heat flux are reliable. Beyond this threshold, because departure from equilibrium may be significant, predictions from these models become doubtful and models from kinetic theory, based on the Boltzmann equation for the molecular velocity distribution function, take precedence. For Knudsen numbers of order unity or higher, solution of the Boltzmann equation either directly or by means of stochastic techniques such as the direct simulation Monte Carlo method (DSMC) of Bird (1994) are computationally tractable, in general. On the other hand, in the interval of Knudsen numbers between, roughly, 0.01 and 1 such computations become increasingly expensive (Torrilhon Reference Torrilhon2016).

Since for gas flows confined in micro- or nano-devices the Knudsen number may lie in the transition regime or beyond, it is important to resort to models capable of describing, at least qualitatively, rarefaction effects when designing or analysing such systems. Other examples include the transport of particles or droplets by a gaseous stream, or of bubbles by a liquid, when the size of these solid or fluid objects lies in the micrometre or nanometre range. Although the Navier–Stokes–Fourier equations of classical hydrodynamics can be extended into the transition regime for Knudsen numbers beyond 0.01 by including slip and jump effects in the boundary conditions, not all rarefaction effects occurring in the bulk can be captured by this approximation (Mohammadzadeh *et al.*
Reference Mohammadzadeh, Rana and Struchtrup2015; Torrilhon Reference Torrilhon2016). Another aspect relevant to modelling micro- or nano-flows is that, for such small scales, the flow Reynolds number is typically smaller than one, so that neglecting inertia is an admissible assumption.

Efforts to model thermophoresis or related phenomena can be traced back for more than a century. The first theory on gas motion induced by temperature gradients was due to Maxwell (Reference Maxwell1879). His analysis stems from kinetic theory. Later, Epstein (Reference Epstein1929) presented a theory for the thermophoretic force on a spherical particle for small Knudsen numbers based on the continuum approach that takes into account thermal slip as well as the particle’s thermal conductivity. Waldmann (Reference Waldmann1959) derived a model valid in the free-molecule regime (large Knudsen numbers) that has been widely applied. Brock (Reference Brock1962) improved Epstein’s theory to develop an expression for the thermophoretic force for Knudsen numbers
${\lesssim}$
0.1. Talbot *et al.* (Reference Talbot, Cheng, Schefer and Willis1980) proposed a correlation for the transition regime by modifying the coefficients in Brock’s expression such that Waldmann’s free-molecule regime expression was approached at large Knudsen numbers.

A commonly used route to study gas flow in the transition regime from first principles, since numerical solutions of the Boltzmann equation can be computationally very costly, has been the analysis of linearized versions of the Boltzmann equation with simpler models for the collision term. Some of the simpler models applied to the problem of thermophoresis of a spherical particle include the Bhatnager, Gross, and Krook model (BGK) (Sone & Aoki Reference Sone and Aoki1983; Yamamoto & Ishihara Reference Yamamoto and Ishihara1988; Takata *et al.*
Reference Takata, Sone and Aoki1993; Sone Reference Sone2007) the S model (Beresnev & Chernyak Reference Beresnev and Chernyak1995) and the hard-sphere model (Sone Reference Sone2007).

An alternative approach to model rarefied gas flow is the use of macroscopic transport equations derived using moment methods. In the original moment method, introduced by Grad (Reference Grad1949), the distribution function in the Boltzmann equation was expanded in Hermite polynomials and the macroscopic variables describing the flow were represented as moments (integrals) of this distribution. For the first approximation beyond the Navier–Stokes–Fourier equations, in total, thirteen moments are needed for the same number of fields, namely, mass density, macroscopic velocity vector, temperature, heat-flux vector and deviatoric stress tensor (symmetric and trace free), yielding Grad’s 13-moment equations (G13). The pressure appears through an equation of state, typically, the ideal gas law. Dwyer (Reference Dwyer1967) presented a theory for the thermophoretic force on a spherical particle based on Grad’s moment method predicting, for the first time, reversed thermophoresis. Recently, Young (Reference Young2011) noted that Dwyer did not account for the totality of the stress and heat-flux coupling terms in the temperature jump boundary condition. Young corrected this and rederived the expression for thermophoretic force from the G13 equations. He then modified this result to present an interpolation formula that fits Waldmann’s expression for large Knudsen numbers and, by introducing values for the thermal creep, velocity slip, and temperature jump coefficients cited by Sharipov (Reference Sharipov2004) based on solutions of model Boltzmann equations, also matches results from kinetic theory for small Knudsen numbers.

A notorious deficiency of the G13 equations is their inability to describe Knudsen layers, that is, regions adjacent to solid surfaces where rarefaction effects are conspicuous. Struchtrup & Torrilhon (Reference Struchtrup and Torrilhon2003) regularized Grad’s 13-moment equations (R13) by adding second-order derivatives to the closures. Starting with the Boltzmann equation, R13 equations are best derived using the order of magnitude method (Struchtrup Reference Struchtrup2005*a*
; Struchtrup *et al.*
Reference Struchtrup, Beckmann, Rana and Frezzotti2017). The R13 equations are equipped with a set of boundary conditions which are valid at walls or gas–liquid interfaces without mass transfer (Torrilhon & Struchtrup Reference Torrilhon and Struchtrup2008) or at evaporating and condensing interfaces (Struchtrup & Frezzotti Reference Struchtrup and Frezzotti2016; Struchtrup *et al.*
Reference Struchtrup, Beckmann, Rana and Frezzotti2017). In contrast to G13, these equations can partially capture the structure and effects of the Knudsen layer. Because they involve the typical variables describing fluid flow and heat transfer at the macroscopic level, interpretation of results may be facilitated by inspecting specific terms in the differential equations, contributing to a good physical understanding. In addition, the extension of tested, well-established numerical techniques in computational fluid dynamics to this system of equations may be achieved straightforwardly. Also, because the boundary conditions have been derived for sharp surfaces of discontinuity, R13 equations can be applied to two-phase systems involving a rarefied gas and a liquid sharing an interface whose instantaneous position is another unknown in the problem. With R13, the limit of continuum models to give meaningful results when rarefaction effects are important has been pushed to a Knudsen number of approximately 0.5 in the transition regime (Torrilhon Reference Torrilhon2016; Struchtrup *et al.*
Reference Struchtrup, Beckmann, Rana and Frezzotti2017). Early advances in the development of moment methods, with emphasis on R13, can be found in the textbook by Struchtrup (Reference Struchtrup2005*b*
); more recent developments and applications of R13 have been compiled in the reviews by Struchtrup & Taheri (Reference Struchtrup and Taheri2011) and Torrilhon (Reference Torrilhon2016).

To the best of our knowledge, application of R13 to investigate transport phenomena involving spherical or near spherical particles in rarefied environments is limited to the analytical work of Torrilhon (Reference Torrilhon2010) on the slow flow of gas past a sphere, and to the numerical treatment of the same problem by Claydon *et al.* (Reference Claydon, Shrestha, Rana, Sprittles and Lockerby2017) using a mesh-free method. The application of R13 to model thermophoresis on a spherical particle has not yet been pursued. After recognizing the complexity of the R13 equations in comparison to G13, Young (Reference Young2011) carried out the modelling of this problem with the latter. In the conclusions of his article, he recommends ‘solving the R13-moment equations in order to study reversed thermophoresis in greater detail’.

The aim of this work is to obtain the thermophoretic force acting on a sphere surrounded by a rarefied gas exposed to a uniform temperature gradient far from the sphere using the R13 equations and taking into account heat conductivity inside the sphere. In addition, we compute the thermophoretic velocity of the sphere when this is free to move under the thermophoretic force. This velocity corresponds to the balance between the thermophoretic force and the drag caused on the sphere in its motion by the surrounding gas. To model the drag, we expand the scope in Torrilhon’s (Reference Torrilhon2010) work with R13 to include the thermal conductivity of the solid (i.e. a non-isothermal sphere), even though the drag has been shown to be fairly insensitive to changes in this parameter (Sone Reference Sone2007). Instead of using the form of the solution for classical Stokes flow past a sphere to write the ansatz for the equations (cf. Torrilhon Reference Torrilhon2010), we apply a somewhat more general approach, namely, the method of multipole potentials. We present closed-form expressions for the thermophoretic force and the drag resulting from R13. We compare predictions for the thermophoretic force from this model with results from simplified models of the Boltzmann equation used in kinetic theory and from other systems of moment equations. We include in the comparison the new experimental data by Bosworth *et al.* (Reference Bosworth, Ventura, Ketsdever and Gimelshein2016) for both positive and negative thermophoresis.

The content of this paper is organized as follows. In the next section, we separately formulate the problems of thermophoresis of a spherical particle and of slow flow of a rarefied gas past a sphere and introduce the main tool of analysis, the regularized 13-moment equations, or R13, and their boundary conditions for the gas–solid interface in linearized form. We then rewrite the system of equations in a different form, more convenient for the solution method of the following section. In § 3, we describe the solution of the system of equations; the procedure involves the method of multipole potentials. Next, § 4 begins by discussing results for the problem of thermophoresis on a sphere, including spatial profiles for the macroscopic field variables, contour plots and streamline patterns and the thermophoretic force, exploring the effect of the solid-to-gas thermal conductivity ratio. For the force, we compare results from R13 with recent experimental data showing reversed thermophoresis and with other models from the literature. Then, we present results for the drag force arising from the gas flow past a sphere from R13 considering the sphere’s heat conductivity and compare with theoretical predictions from the literature, including Torrilhon’s (Reference Torrilhon2010) R13 results in the case of an isothermal particle, which serves as validation of the solution method implemented here. The section closes presenting results for the thermophoretic velocity from R13 and other models. Finally, § 5 contains some concluding remarks.

## 2 Problem formulation

In this section we formulate mathematically two problems involving a sphere in a rarefied gas. First, we consider the problem of thermophoresis of a spherical particle with the gas far from the sphere at rest with a uniform temperature gradient. Second, we detail the problem of a uniform flow past a sphere with no temperature gradient imposed in the far field. The results of these two problems will be combined later to obtain the thermophoretic velocity of a sphere.

### 2.1 Thermophoresis on a sphere by a uniform temperature gradient in the far field

Consider a gas at rest with uniform, constant pressure $p_{0}^{\ast }$ and temperature $T_{0}^{\ast }$ surrounding a sphere of radius $a^{\ast }$ at the same temperature and motionless with respect to the laboratory frame. Under these conditions, the gas is in equilibrium with vanishing heat flux and deviatoric stress. In energy units, the temperature in this state is given by $\unicode[STIX]{x1D703}_{0}^{\ast }=R^{\ast }T_{0}^{\ast }$ , where $R^{\ast }$ is the gas specific constant. Suppose that this state of equilibrium is disturbed by imposing, far from the sphere, a temperature field, $T_{0}^{\ast }+z^{\ast }(\unicode[STIX]{x2202}T^{\ast }/\unicode[STIX]{x2202}z^{\ast })_{\infty }$ , with $(\unicode[STIX]{x2202}T^{\ast }/\unicode[STIX]{x2202}z^{\ast })_{\infty }$ constant, and where plane $z^{\ast }=0$ passes through the sphere’s centre (figure 1). The viscosity and thermal conductivity coefficient of the gas evaluated at the equilibrium state are denoted by $\unicode[STIX]{x1D707}_{0}^{\ast }$ and $k_{0}^{\ast }$ , respectively. The thermal conductivity coefficient for the sphere’s material is denoted by $k_{s}^{\ast }$ (constant). The gas is assumed to be ideal, so that in the equilibrium state the gas density is $\unicode[STIX]{x1D70C}_{0}^{\ast }=p_{0}^{\ast }/\unicode[STIX]{x1D703}_{0}^{\ast }$ . Assume also that the ratio of the gas molecules’ mean free path to the sphere radius is such that rarefaction effects cannot be ignored. On the basis of this ratio, we introduce a Knudsen number

Note that, from kinetic theory, a commonly used definition for the gas mean free path in the undisturbed state is $(\unicode[STIX]{x03C0}/2)^{1/2}\unicode[STIX]{x1D707}_{0}^{\ast }\,\unicode[STIX]{x1D703}_{0}^{\ast \,1/2}/p_{0}^{\ast }$ . In definition of (2.1), we have dropped the factor $(\unicode[STIX]{x03C0}/2)^{1/2}$ for simplicity.

In what follows, we consider the governing equations for a monatomic gas composed of Maxwell molecules. In this case, the Prandtl number
$Pr=2/3$
. From the well-known definition of
$Pr$
, a useful relationship between the gas thermal conductivity and dynamic viscosity can be obtained (Struchtrup Reference Struchtrup2005*b*
)

For a monatomic gas, the ratio of specific heats is $\unicode[STIX]{x1D6FE}=5/3$ (e.g. see Young Reference Young2011).

With the aim of modelling flow and heat transfer phenomena in a rarefied gas, we consider the conservation laws for mass, momentum and energy, supplemented by the constitutive equations for the deviatoric stress and heat flux from the R13 theory, and the associated augmented set of boundary conditions for a gas–solid interface. We are interested here in the steady-state gas flow and temperature fields, in the gas and solid, resulting from the far-field temperature gradient and gas rarefaction. Assuming that the dimensionless group $a^{\ast }(\unicode[STIX]{x2202}T^{\ast }/\unicode[STIX]{x2202}z^{\ast })_{\infty }/T_{0}^{\ast }\ll 1$ , we can model the transport phenomena using the linearized version of the governing equations written in terms of deviations from the equilibrium state. The non-dimensional form of these deviations can be written as

for the pressure, temperature, density, velocity, heat flux and deviatoric stress, respectively, in the gas. The deviatoric stress is symmetric and trace free. Length is non-dimensionalized with the sphere’s radius $a^{\ast }$ . The temperature deviation in the sphere $\unicode[STIX]{x1D703}_{s}^{\ast }$ is non-dimensionalized as that for the gas. The dimensionless form of the temperature gradient defines a new dimensionless group, the Epstein number – coined by Young (Reference Young2011) after P. S. Epstein, a pioneer in the study of thermophoresis of spherical particles, who presented the first theory on the subject (Epstein Reference Epstein1929). It is given by

Note that the non-dimensional pressure and temperature fields in the gas are given by $1+p$ and $1+\unicode[STIX]{x1D703}$ , respectively, whereas the temperature in the sphere is $1+\unicode[STIX]{x1D703}_{s}$ . On the other hand, $\text{ u}$ , $\text{q}$ and $\unicode[STIX]{x1D748}$ in (2.3) determine the actual velocity, heat flux and deviatoric stress in the gas.

The linearized, steady conservation equations of gas, momentum and energy are

*a*) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\text{u}=0, & \displaystyle\end{eqnarray}$$

*b*) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}p+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}=0, & \displaystyle\end{eqnarray}$$

*c*) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\text{q}=0, & \displaystyle\end{eqnarray}$$

*b*; Lockerby & Collyer Reference Lockerby and Collyer2016; Torrilhon Reference Torrilhon2016)

*a*) $$\begin{eqnarray}\displaystyle & \displaystyle \left(1-{\textstyle \frac{2}{3}}Kn^{2}\unicode[STIX]{x1D6E5}\right)\unicode[STIX]{x1D748}-{\textstyle \frac{4}{5}}Kn^{2}\overline{\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}}=-2Kn\overline{\unicode[STIX]{x1D735}\text{u}}-{\textstyle \frac{4}{5}}Kn\overline{\unicode[STIX]{x1D735}\text{q}}, & \displaystyle\end{eqnarray}$$

*b*) $$\begin{eqnarray}\displaystyle & \displaystyle \left(1-{\textstyle \frac{9}{5}}Kn^{2}\unicode[STIX]{x1D6E5}\right)\text{q}=-{\textstyle \frac{15}{4}}Kn\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}-{\textstyle \frac{3}{2}}Kn\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}. & \displaystyle\end{eqnarray}$$

Far from the sphere’s surface, the imposed disturbance is represented by the uniform temperature gradient $Ep\,\text{k}$ , where $\text{k}$ designates the unit vector pointing in the direction of the positive $z$ semi-axis. It is convenient to write the problem in such a way that deviations from the base (equilibrium) state vanish in the far field. For this purpose we introduce the transformations

*a*-

*c*) $$\begin{eqnarray}\unicode[STIX]{x1D703}=\check{\unicode[STIX]{x1D703}}+Ep\,z,\quad \unicode[STIX]{x1D703}_{s}=\check{\unicode[STIX]{x1D703}}_{s}+Ep\,z,\quad \text{q}=\check{\text{q}}-{\textstyle \frac{15}{4}}Kn\,Ep\,\text{k},\end{eqnarray}$$

whereas the rest of the variables are left unchanged. With (2.8), expressions (2.5)–(2.7) remain invariant in form. On the other hand, the boundary conditions do change after using (2.8). Note that subscript ‘ $s$ ’ will be used to denote quantities in the solid sphere.

Equations (2.5)–(2.7) will be solved subjected to the following boundary conditions taken from Struchtrup *et al.* (Reference Struchtrup, Beckmann, Rana and Frezzotti2017) in the absence of phase change (see also Struchtrup & Frezzotti Reference Struchtrup and Frezzotti2016) at the solid–gas interface
$r=1$
, with
$r=|\text{x}|$
, where
$\text{x}$
is the position vector with origin at the sphere’s centre. The quantities at the interface corresponding to the liquid in Struchtrup *et al.* (Reference Struchtrup, Beckmann, Rana and Frezzotti2017) are, in the present work, associated with the solid (sphere). Note that in this case, boundary conditions (37)–(41) in Struchtrup *et al.* (Reference Struchtrup, Beckmann, Rana and Frezzotti2017) reduce to the wall boundary conditions (33)–(37) in Torrilhon & Struchtrup (Reference Torrilhon and Struchtrup2008) if the gas flow is two-dimensional. In all these equations for the boundary the accommodation coefficient has been set equal to one and
$\unicode[STIX]{x1D703}$
,
$\unicode[STIX]{x1D703}_{s}$
and
$\text{q}$
have been substituted according to (2.8).

The generalized slip condition is

The generalized temperature jump condition is

The generalized interface conditions for higher moments are

and

Furthermore, the solution must also satisfy the non-penetration condition

and the interfacial linearized energy balance (Young Reference Young2011)

where $\unicode[STIX]{x1D6EC}$ is the solid-to-gas thermal conductivity ratio $k_{s}^{\ast }/k_{0}^{\ast }$ . Here, indices $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D6FD}$ can take a value of 1 or 2; $\text{n}$ is the unit vector normal to the interface pointing into the gas; $\text{t}_{1}$ and $\text{t}_{2}$ represent two mutually orthogonal unit vectors tangential to the interface, and subscripts $n$ , $t_{1}$ and $t_{2}$ denote components in the corresponding normal or any of the two tangential directions, respectively. In these interfacial conditions one must set $u_{s,t_{1}}=u_{s,t_{2}}=0$ . Regarding the far field, all deviations from the basic equilibrium state must vanish as $r\rightarrow \infty$ .

Boundary conditions (2.14) and (2.15) hold for the interface between a fluid and an impenetrable solid. Their counterparts for a fluid–fluid interface are given by the conditions for mass and energy conservation presented, for instance, in appendix A of Struchtrup *et al.* (Reference Struchtrup, Beckmann, Rana and Frezzotti2017) – note that in their expression (A5), the internal energy should be written instead of the enthalpy.

The interfacial conditions contain the components of the higher-order moments $\text{R}$ and $\text{m}$ , a rank-two and rank-three tensor, respectively. These are defined as

*a*) $$\begin{eqnarray}\displaystyle & \text{R}=-{\textstyle \frac{24}{5}}Kn\,\overline{\unicode[STIX]{x1D735}\check{\text{q}}}, & \displaystyle\end{eqnarray}$$

*b*) $$\begin{eqnarray}\displaystyle & \text{m}=-2Kn\,\overline{\unicode[STIX]{x1D735}\unicode[STIX]{x1D748}}. & \displaystyle\end{eqnarray}$$

*b*, see chapters 7 and 9) – there, a third (scalar) higher-order moment is also present that contributes nothing to the linear equations when $\text{q}$ is divergence free. Once these equations are linearized, relations (2.16) are then used to eliminate $\text{R}$ and $\text{m}$ resulting in (2.6).

It is important to note that because $\unicode[STIX]{x1D748}$ , $\text{R}$ and $\text{m}$ are trace-free tensors, we have that $\unicode[STIX]{x1D70E}_{nn}=-\unicode[STIX]{x1D70E}_{t_{1}t_{1}}-\unicode[STIX]{x1D70E}_{t_{2}t_{2}}$ , $\text{R}_{nn}=-\text{R}_{t_{1}t_{1}}-\text{R}_{t_{2}t_{2}}$ and $\text{m}_{nnn}=-\text{m}_{t_{1}t_{1}n}-\text{m}_{t_{2}t_{2}n}$ . Using these constraints, we can easily show that boundary condition (2.11) can be obtained by writing (2.12) twice, first for $\text{m}_{t_{1}t_{1}n}$ and then for $\text{m}_{t_{2}t_{2}n}$ , and adding the resulting expressions. Therefore, it suffices for the solution to satisfy only one of these equations provided (2.11) is also satisfied.

Since these boundary conditions will be applied to a spherical interface, it is fitting to introduce a system of spherical-polar coordinates $(r,\unicode[STIX]{x1D717},\unicode[STIX]{x1D711})$ with its origin, $r=0$ , located at the centre of the solid sphere (figure 1). Here, $0\leqslant r<\infty$ , $0\leqslant \unicode[STIX]{x1D717}\leqslant \unicode[STIX]{x03C0}$ , and $0\leqslant \unicode[STIX]{x1D711}<2\unicode[STIX]{x03C0}$ . Semi-axis $\unicode[STIX]{x1D717}=0$ coincides with the positive $z$ semi-axis. To these coordinate directions correspond unit vectors $\hat{\text{r}}$ , $\hat{\unicode[STIX]{x1D751}}$ and $\hat{\unicode[STIX]{x1D753}}$ , respectively; these triplet forms an orthogonal set. Thus $\{\hat{\text{r}},\hat{\unicode[STIX]{x1D751}},\hat{\unicode[STIX]{x1D753}}\}$ take the place of the set $\{\text{n},\text{t}_{1},\text{t}_{2}\}$ when writing the boundary conditions. The problems considered here are axisymmetric, so quantities do not vary in the $\unicode[STIX]{x1D711}$ direction.

Finally, if needed, the gas density deviation from its value at equilibrium can be computed from the linearized form of the ideal gas equation of state, $p=\unicode[STIX]{x1D70C}+\unicode[STIX]{x1D703}$ .

### 2.2 Uniform flow past a sphere

Suppose that instead of prescribing a far-field temperature gradient, the state of equilibrium of the gas described in the previous section is disturbed by imposing, far from the sphere, a uniform flow with constant velocity $U_{0}^{\ast }$ in the $z^{\ast }$ -direction. We can model this flow, including rarefaction effects, by using the linearized conservation laws, R13 constitutive relations and boundary conditions written for the disturbances, provided the dimensionless (Reynolds) number $\unicode[STIX]{x1D70C}_{0}^{\ast }U_{0}^{\ast }a^{\ast }/\unicode[STIX]{x1D707}_{0}^{\ast }\ll 1$ . Again, steady state will be assumed. Our work will extend Torrilhon’s (Reference Torrilhon2010) efforts on this problem by including heat conduction throughout the sphere. Furthermore, we will present an expression for the drag force over the sphere from our analytical solution – a closed-form expression for the drag was not given in Torrilhon’s work. Combining this expression for the drag with that for the thermophoretic force will result in the thermophoretic velocity when these two forces balanced each other.

For convenience, we transform the original problem to the equivalent one of the disturbance flow resulting from a sphere translating with dimensionless velocity $-Ma\,\text{k}$ in a fluid at rest far away from the sphere. Here, $Ma$ is a pseudo Mach number

The actual Mach number can be obtained by multiplying $Ma$ by $\unicode[STIX]{x1D6FE}^{-1/2}$ , with $\unicode[STIX]{x1D6FE}$ the specific heat ratio.

We thus set $\text{u}=\check{\text{u}}+Ma\text{k}$ , so that $\check{\text{u}}\rightarrow 0$ as $r\rightarrow \infty$ . The governing equations and boundary conditions for this problem are obtained from those in § 2.1 by writing $\check{\text{u}}$ instead of $\text{u}$ , dropping the ‘ $\check{~}\,$ ’ from $\check{\unicode[STIX]{x1D703}}$ , $\check{\unicode[STIX]{x1D703}}_{s}$ and $\check{\text{q}}$ and setting $Ep=0$ , ${\check{u}}_{s,t_{1}}={\check{u}}_{s,\unicode[STIX]{x1D717}}=-Ma\,\text{k}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D751}}$ and ${\check{u}}_{s,t_{2}}={\check{u}}_{s,\unicode[STIX]{x1D711}}=0$ . In addition, instead of (2.14) we must have

Hereinafter, we drop symbol ‘ $\check{~}\,$ ’ unless otherwise noted.

### 2.3 Equivalent Stokes–Fourier system of equations

We shall proceed to rewrite the linearized R13-moment equations introduced previously as an equivalent Stokes–Fourier set of equations supplemented by non-homogeneous elliptic equations for the pressure, heat flux and deviatoric stress. To this alternative system of equations, we can apply, somewhat straightforwardly, analytical tools used successfully in scalar and, more notably, vector equations appearing in slow flow hydrodynamics. By introducing the auxiliary variables

it can be readily shown that expressions (2.5)–(2.6) can be written as the equivalent Stokes–Fourier system

*a*) $$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D755}=0, & \displaystyle\end{eqnarray}$$

*b*) $$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D735}\unicode[STIX]{x1D6F1}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D72E}=0, & \displaystyle\end{eqnarray}$$

*c*) $$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D734}=0, & \displaystyle\end{eqnarray}$$

*a*) $$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D72E}=-2Kn\overline{\unicode[STIX]{x1D735}\unicode[STIX]{x1D755}}, & \displaystyle\end{eqnarray}$$

*b*) $$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D734}=-{\textstyle \frac{15}{4}}Kn\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}, & \displaystyle\end{eqnarray}$$

Furthermore, because of (2.21*a*
) and (2.21*c*
), (2.21*b*
) and (2.22*b*
) lead to

and, with (2.22*a*
), expression (2.21*b*
) becomes

Equations (2.21*a*
) and (2.21*c*
) result from (2.5*a*
), (2.5*c*
) and from the divergence of (2.6*b*
). Expression (2.21*b*
) is obtained by taking the divergence of (2.6*a*
), using (2.5*b*
) to eliminate
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}$
in favour of
$\unicode[STIX]{x1D735}p$
, collecting the terms containing this vector and introducing
$\unicode[STIX]{x1D6F1}$
defined in (2.23). Finally, expressions (2.24) and (2.25) are simply (2.6*a*
) and (2.6*b*
) written in terms of auxiliary variables
$\unicode[STIX]{x1D734}$
and
$\unicode[STIX]{x1D72E}$
, respectively. Rewriting the R13-moment equations in a way that includes the more familiar form of (2.21) or, more recognizably, (2.26)–(2.28), will facilitate the application to this set of equations of analytical methods used in the solution of the equations for Stokes flow in vector form. We shall extend these methods to (2.23)–(2.25), which include a rank-two tensor, second-order differential equation.

## 3 Solution applying the method of multipole potentials

In this section, we seek solutions for the flow disturbances in the gas that vanish in the far field alongside temperature profiles inside the sphere which are coupled through the set of conditions prescribed at the sphere’s surface.

Consider partial differential equations (2.26)–(2.28) and (2.23)–(2.25) for the gas, and (2.7) for the temperature in the solid particle. Based on this, scalar fields $\unicode[STIX]{x1D6F1}$ , $\unicode[STIX]{x1D701}$ and $\unicode[STIX]{x1D703}_{s}$ are harmonic functions. For fields $\unicode[STIX]{x1D755}$ , $p$ , $\text{q}$ , $\unicode[STIX]{x1D748}$ , their solutions will be written as the sum of the solution of the homogeneous equation plus a particular solution of (2.28) and (2.23)–(2.25), respectively (e.g. see Lamb Reference Lamb1932; Leal Reference Leal2007). The solutions of the homogeneous equations will be found with the method of multipole potentials. After these fields have been determined, the velocity and temperature fields $\text{u}$ and $\unicode[STIX]{x1D703}$ can be obtained from expressions (2.19) and (2.20), respectively.

### 3.1 Particular solutions

Starting with (2.28), to find a particular solution we let

where $\unicode[STIX]{x1D6FC}_{i}$ is a vector field to be determined and the subscript $i$ denotes Cartesian index notation. Computing

after using the fact that $\unicode[STIX]{x1D6F1}$ is harmonic, and comparing this result with the left-hand side of (2.28), namely, $Kn^{-1}\unicode[STIX]{x2202}\unicode[STIX]{x1D6F1}/\unicode[STIX]{x2202}x_{i}$ , we find $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FC}_{i}=0$ and $\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}_{i}/\unicode[STIX]{x2202}x_{j}=\unicode[STIX]{x1D6FF}_{ij}/2$ , where $\unicode[STIX]{x1D6FF}_{ij}$ is the Kronecker delta. Therefore, $\unicode[STIX]{x1D6FC}_{i}=x_{i}/2$ and (Leal Reference Leal2007)

with $\text{x}$ the position vector from the centre of the sphere. Next, for (2.23) and (2.24), solutions $p^{(p)}$ and $\text{q}^{(p)}$ can be readily found by inspection by using the fact that $\unicode[STIX]{x1D6F1}$ and $\unicode[STIX]{x1D701}$ are harmonic functions. This leads to

and

respectively.

Finally, solution
$\unicode[STIX]{x1D748}^{(p)}$
for (2.25), after eliminating
$\unicode[STIX]{x1D72E}$
with (2.22*a*
), can be found by letting

such that vector field $\unicode[STIX]{x1D737}$ and scalar field $\unicode[STIX]{x1D702}$ satisfy, respectively,

*a*) $$\begin{eqnarray}\displaystyle & \displaystyle (\unicode[STIX]{x1D6E5}-\unicode[STIX]{x1D706}_{3}^{2})\unicode[STIX]{x1D737}=\unicode[STIX]{x1D755}, & \displaystyle\end{eqnarray}$$

*b*) $$\begin{eqnarray}\displaystyle & \displaystyle (\unicode[STIX]{x1D6E5}-\unicode[STIX]{x1D706}_{3}^{2})\unicode[STIX]{x1D702}=p. & \displaystyle\end{eqnarray}$$

*b*), we seek solutions of (3.7) of the form

*a*) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D737}=-\unicode[STIX]{x1D755}/\unicode[STIX]{x1D706}_{3}^{2}+m\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F1}, & \displaystyle\end{eqnarray}$$

*b*) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D702}=p/(\unicode[STIX]{x1D706}_{1}^{2}-\unicode[STIX]{x1D706}_{3}^{2})+n\unicode[STIX]{x1D6F1}. & \displaystyle\end{eqnarray}$$

*a*) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D737}=-{\textstyle \frac{1}{\unicode[STIX]{x1D706}_{3}^{4}}}(\unicode[STIX]{x1D706}_{3}^{2}\unicode[STIX]{x1D755}+Kn^{-1}\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F1}), & \displaystyle\end{eqnarray}$$

*b*) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D702}=\frac{1}{\unicode[STIX]{x1D706}_{1}^{2}-\unicode[STIX]{x1D706}_{3}^{2}}\left(p-\frac{\unicode[STIX]{x1D706}_{1}^{2}}{\unicode[STIX]{x1D706}_{3}^{2}}\unicode[STIX]{x1D6F1}\right), & \displaystyle\end{eqnarray}$$

### 3.2 Solution of the homogeneous equations

The form of the solutions of the homogeneous equations associated with (2.28) and (2.23)–(2.25) is constructed by means of the method of multipole potentials. This method is presented in the context of classical hydrodynamics and applied to problems of low Reynolds number flows in Leal (Reference Leal2007, chap. 8) as well as thoroughly discussed and used in a wider range of problems of mathematical physics in Hess (Reference Hess2015, chap. 10). In the former, it is named as the method of superposition of vector harmonic functions. The most illustrative example discussed in those sources is perhaps the problem of Stokes flow past a sphere. As we shall see, a significant advantage of the method of multipole potentials concerning our problem is that we can obtain the solution of the deviatoric stress $\unicode[STIX]{x1D748}$ without going through the complexity of writing and solving differential equations for its scalar components in spherical coordinates. The general definitions of multipole potential tensors, the form of the first few of them and properties relevant to this work are presented in appendix A.

Because the governing equations and boundary conditions are linear in the perturbation fields, and the non-homogeneous terms in the boundary conditions must be linear in
$G\text{k}$
, with
$G=Ep$
or
$G=Ma$
for the problems in §§ 2.1 or 2.2, respectively, the homogeneous solutions are constructed by adding products of multipole potentials with vector
$G\text{k}$
. From the set of all multipole potential tensors, the method of multipole potential guides the selection of those potentials that should be used to construct the solution. Only those multipole potentials that conform to the rank of the unknown field (i.e. whether it is a scalar, a vector, a rank-two tensor, etc.), its symmetry, and its parity (i.e. whether it is a true scalar, true vector or true tensor or a pseudo-scalar, pseudo-vector or pseudo-tensor), after their product with
$G\text{k}$
can be part of the solution. Regarding the parity attribute, in general terms, the components of a true scalar, vector or tensor change sign when subjected to an improper rotation, such as a reflection or inversion – e.g. changing the coordinate system from right-handed to left-handed – in order to continue properly describing an unchanged physical situation (Arfken *et al.*
Reference Arfken, Weber and Harris2012, chap. 3). For a so-called pseudo-scalar, vector or tensor, their components do not change sign accordingly. For instance, the angular velocity is a pseudo-vector. A deeper look at this aspect and, in general, at the method of multipole potentials for constructing solutions of physical quantities is beyond the scope of the present article. For detailed discussions and a list of enlightening examples, the reader may referred to the cited books by Leal (Reference Leal2007) and Hess (Reference Hess2015).

Following these considerations, using Cartesian index notation, and letting $G_{i}$ be the components of vector $G\text{k}$ , we can write for the solutions

where terms with superscript ‘ $(p)$ ’ are given by (3.3)–(3.6); $A_{1}$ , $B_{1}$ , $C_{1}$ and $C_{2}$ , are constants to be determined by the boundary conditions, and fields $X_{0}$ , $X_{i}$ , $X_{ij},\ldots \,$ , denote the descending multipole potentials listed in appendix A. Note that in (3.13)–(3.15) we use functions of the radial coordinate $a(r)$ , $b(r)$ , $c(r)$ , $f(r)$ and $g(r)$ , instead of constants as for (3.10)–(3.12), because the solutions of the homogeneous equations for $p$ , $\text{q}$ and $\unicode[STIX]{x1D748}$ , must satisfy modified Helmholtz equations rather than Laplace equations. The descending multipole potentials satisfy Laplace equation for $r\neq 0$ . In the present problem, $r=0$ is outside the gas domain.

For the interior of the sphere, the temperature profile is given by the ascending multipole potential

which introduces another constant, $D_{1}$ . Both the equivalence in tensor rank and parity between the left- and right-hand sides of (3.10)–(3.16) determined the type of products employed between the external perturbation vector $G_{i}$ and the multipole potential tensors in these expressions.

Because
$\unicode[STIX]{x1D701}$
is harmonic, equation (3.5) yields
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\text{q}^{(p)}=0$
. With the divergence-free vector
$\text{q}=\text{q}^{(h)}+\text{q}^{(p)}$
, where superscript ‘(
$h$
)’ denotes the homogeneous solution, we are left with
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\text{q}^{(h)}=0$
. Similarly, one can show that taking the divergence of
$\unicode[STIX]{x1D748}^{(p)}$
using (3.6) results in
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}^{(p)}=-\unicode[STIX]{x1D735}p$
. Therefore, since
$\unicode[STIX]{x1D748}=\unicode[STIX]{x1D748}^{(h)}+\unicode[STIX]{x1D748}^{(p)}$
, and because of expression (2.5*b*
), we must also have
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}^{(h)}=0$
. Note in (3.15) that, by construction,
$\unicode[STIX]{x1D748}^{(h)}$
is also symmetric and trace free.

Taking the divergence of $\unicode[STIX]{x1D755}$ using (3.12) and (3.3), setting it equal to zero and using the definitions of the multipole potentials from appendix A results in

With $\text{q}^{(h)}$ from (3.14), setting $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\text{ q}^{(h)}=0$ yields this constraint for functions $b(r)$ and $c(r)$

Furthermore, the vector equation $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}^{(h)}=0$ , with $\unicode[STIX]{x1D748}^{(h)}$ from (3.15), leads to

*a*) $$\begin{eqnarray}\displaystyle & \displaystyle 2r^{2}f^{\prime }-rf+18g^{\prime }=0, & \displaystyle\end{eqnarray}$$

*b*) $$\begin{eqnarray}\displaystyle & 32r^{2}f^{\prime }+rf-18g^{\prime }=0, & \displaystyle\end{eqnarray}$$

*a*) $$\begin{eqnarray}\displaystyle & \displaystyle a^{\prime \prime }-2r^{-1}a^{\prime }-\unicode[STIX]{x1D706}_{1}^{2}a=0, & \displaystyle\end{eqnarray}$$

*b*) $$\begin{eqnarray}\displaystyle & \displaystyle b^{\prime \prime }-\unicode[STIX]{x1D706}_{2}^{2}b=0, & \displaystyle\end{eqnarray}$$

*c*) $$\begin{eqnarray}\displaystyle & \displaystyle c^{\prime \prime }-4r^{-1}c^{\prime }-\unicode[STIX]{x1D706}_{2}^{2}c=0, & \displaystyle\end{eqnarray}$$

*d*) $$\begin{eqnarray}\displaystyle & \displaystyle f^{\prime \prime }-2r^{-1}f^{\prime }-\unicode[STIX]{x1D706}_{3}^{2}f=0, & \displaystyle\end{eqnarray}$$

*e*) $$\begin{eqnarray}\displaystyle & \displaystyle g^{\prime \prime }-6r^{-1}g^{\prime }-\unicode[STIX]{x1D706}_{3}^{2}g=0. & \displaystyle\end{eqnarray}$$

*d*) yields $f(r)\equiv 0$ ; hence, from above, we have that $g^{\prime }(r)\equiv 0$ and (3.20

*e*) results in $g(r)\equiv 0$ . We thus find that $\unicode[STIX]{x1D748}^{(h)}(\text{x})\equiv 0$ . Therefore, $\unicode[STIX]{x1D748}=\unicode[STIX]{x1D748}^{(p)}$ , given in (3.6) which, in turn, makes use of expressions (3.10), (3.12) and (3.13). The solution of (3.20

*b*) is well known. Obtaining exact solutions for (3.20

*a*) and (3.20

*c*) is discussed in appendix B. These solutions can be written as

*a*) $$\begin{eqnarray}\displaystyle & \displaystyle a(r)=E_{1}\exp [-\unicode[STIX]{x1D706}_{1}(r-1)](1+\unicode[STIX]{x1D706}_{1}r), & \displaystyle\end{eqnarray}$$

*b*) $$\begin{eqnarray}\displaystyle & \displaystyle b(r)=F_{1}\exp [-\unicode[STIX]{x1D706}_{2}(r-1)], & \displaystyle\end{eqnarray}$$

*c*) $$\begin{eqnarray}\displaystyle & \displaystyle c(r)=F_{2}\exp [-\unicode[STIX]{x1D706}_{2}(r-1)](1+\unicode[STIX]{x1D706}_{2}r+\unicode[STIX]{x1D706}_{2}^{2}r^{2}/3), & \displaystyle\end{eqnarray}$$

*b*) and (3.21

*c*) into the divergence-free condition for the heat flux (3.18) results in

Note that one could have also written explicit expressions for
$f(r)$
and
$g(r)$
from solving (3.20*d*
) and (3.20*e*
) in a manner similar to that followed to obtain (3.21), only to find out that, after enforcing
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D748}^{(h)}=0$
, the associated integration constants would be equal to zero, thus yielding
$f=g\equiv 0$
.

According to the R13 theory, the structure of the Knudsen layer near a boundary is partially determined by the factors $\exp (-\unicode[STIX]{x1D706}_{i}^{2}r)$ in the solution ( $i=1,2,3$ for the present case). The contribution with $\unicode[STIX]{x1D706}_{3}$ vanishes in this case resulting from the vanishing of $\unicode[STIX]{x1D748}^{(h)}$ . The factor containing the exponential with $\unicode[STIX]{x1D706}_{1}$ is present in the pressure, temperature and, through the pressure, in the deviatoric stress. On the other hand, the factor containing the exponential with $\unicode[STIX]{x1D706}_{2}$ occurs in the heat flux and velocity. All these features were already noted by Torrilhon (Reference Torrilhon2010) in the solution for the special case of slow flow past an isothermal sphere.

### 3.3 Final field expressions – system of equations for the integration constants

We can now write the final form of the exact solutions for the field variables. The pressure deviation in the gas $p$ can be computed from (3.13), using (3.4) and (3.10). The temperature deviation $\unicode[STIX]{x1D703}$ is found from (2.20), with (3.11) and (3.13). The expression for the gas velocity $\text{u}$ is obtained from (2.19), using (3.12) and (3.14), together with (3.3), (3.10), (3.5) and (3.11). The heat-flux vector $\text{q}$ is computed using (3.14), (3.5) and (3.12), whereas the stress tensor $\unicode[STIX]{x1D748}$ is determined from expressions (3.15) and (3.6), with (3.9), (3.10), (3.12) and (3.13). The temperature deviation in the solid $\unicode[STIX]{x1D703}_{s}$ is given by (3.16). After substitution of (3.17), (3.21) and (3.22), the final expressions for the deviations are

*a*) $$\begin{eqnarray}\displaystyle p & = & \displaystyle A_{1}Gr^{-2}\cos \unicode[STIX]{x1D717}+E_{1}Gr^{-2}(\unicode[STIX]{x1D706}_{1}r+1)\exp [-\unicode[STIX]{x1D706}_{1}(r-1)]\cos \unicode[STIX]{x1D717},\end{eqnarray}$$

*b*) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D703} & = & \displaystyle {\textstyle \frac{1}{5}}(2A_{1}+5B_{1})Gr^{-2}\cos \unicode[STIX]{x1D717}+{\textstyle \frac{2}{5}}E_{1}Gr^{-2}(\unicode[STIX]{x1D706}_{1}r+1)\exp [-\unicode[STIX]{x1D706}_{1}(r-1)]\cos \unicode[STIX]{x1D717},\qquad\end{eqnarray}$$

*c*) $$\begin{eqnarray}\displaystyle u_{r} & = & \displaystyle Kn^{-1}Gr^{-3}(A_{1}r^{2}-3Kn^{2}B_{1}+2KnC_{2})\cos \unicode[STIX]{x1D717}\nonumber\\ \displaystyle & & \displaystyle -\,{\textstyle \frac{4}{5}}F_{2}Gr^{-3}(\unicode[STIX]{x1D706}_{2}r+1)\exp [-\unicode[STIX]{x1D706}_{2}(r-1)]\cos \unicode[STIX]{x1D717},\end{eqnarray}$$

*d*) $$\begin{eqnarray}\displaystyle u_{\unicode[STIX]{x1D717}} & = & \displaystyle -{\textstyle \frac{1}{2}}Kn^{-1}Gr^{-3}(A_{1}r^{2}+3Kn^{2}B_{1}-2KnC_{2})\sin \unicode[STIX]{x1D717}\nonumber\\ \displaystyle & & \displaystyle -\,{\textstyle \frac{2}{5}}F_{2}Gr^{-3}(\unicode[STIX]{x1D706}_{2}^{2}r^{2}+\unicode[STIX]{x1D706}_{2}r+1)\exp [-\unicode[STIX]{x1D706}_{2}(r-1)]\sin \unicode[STIX]{x1D717},\end{eqnarray}$$

*e*) $$\begin{eqnarray}\displaystyle q_{r} & = & \displaystyle {\textstyle \frac{15}{2}}KnB_{1}Gr^{-3}\cos \unicode[STIX]{x1D717}+2F_{2}Gr^{-3}(\unicode[STIX]{x1D706}_{2}r+1)\exp [-\unicode[STIX]{x1D706}_{2}(r-1)]\cos \unicode[STIX]{x1D717},\end{eqnarray}$$

*f*) $$\begin{eqnarray}\displaystyle q_{\unicode[STIX]{x1D717}} & = & \displaystyle {\textstyle \frac{15}{4}}KnB_{1}Gr^{-3}\sin \unicode[STIX]{x1D717}+F_{2}Gr^{-3}(\unicode[STIX]{x1D706}_{2}^{2}r^{2}+\unicode[STIX]{x1D706}_{2}r+1)\exp [-\unicode[STIX]{x1D706}_{2}(r-1)]\sin \unicode[STIX]{x1D717},\end{eqnarray}$$

*g*) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70E}_{rr} & = & \displaystyle \displaystyle \frac{2}{5}Gr^{-4}(5A_{1}r^{2}-12Kn^{2}A_{1}-30\unicode[STIX]{x1D706}_{3}^{-2}A_{1}+30KnC_{2})\cos \unicode[STIX]{x1D717}+\displaystyle \frac{8}{15}\frac{\unicode[STIX]{x1D706}_{3}^{2}}{\unicode[STIX]{x1D706}_{1}^{2}-\unicode[STIX]{x1D706}_{3}^{2}}\nonumber\\ \displaystyle & & \displaystyle \times \,Kn^{2}E_{1}Gr^{-4}(\unicode[STIX]{x1D706}_{1}^{3}r^{3}+4\unicode[STIX]{x1D706}_{1}^{2}r^{2}+9\unicode[STIX]{x1D706}_{1}r+9)\exp [-\unicode[STIX]{x1D706}_{1}(r-1)]\cos \unicode[STIX]{x1D717},\end{eqnarray}$$

*h*) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D717}r} & = & \displaystyle -{\textstyle \frac{6}{5}}Gr^{-4}(2Kn^{2}A_{1}+5\unicode[STIX]{x1D706}_{3}^{-2}A_{1}-5KnC_{2})\sin \unicode[STIX]{x1D717}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{4}{5}\frac{\unicode[STIX]{x1D706}_{3}^{2}}{\unicode[STIX]{x1D706}_{1}^{2}-\unicode[STIX]{x1D706}_{3}^{2}}Kn^{2}E_{1}Gr^{-4}(\unicode[STIX]{x1D706}_{1}^{2}r^{2}+3\unicode[STIX]{x1D706}_{1}r+3)\exp [-\unicode[STIX]{x1D706}_{1}(r-1)]\sin \unicode[STIX]{x1D717},\end{eqnarray}$$

*i*) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D703}_{s} & = & \displaystyle D_{1}Gr\cos \unicode[STIX]{x1D717}.\end{eqnarray}$$

To compute the higher-order moment tensors
$\text{R}$
and
$\text{m}$
defined in (2.16*a*
) and (2.16*b*
), respectively, we first calculate rank-two and rank-three tensors
$\unicode[STIX]{x1D735}\text{q}$
and
$\unicode[STIX]{x1D735}\unicode[STIX]{x1D748}$
, respectively. We compute these gradients by means of a popular symbolic algebra and calculus computer package (Wolfram Mathematica), leading to expressions that we do not reproduce here for the sake of space. An alternative path to
$\unicode[STIX]{x1D735}\unicode[STIX]{x1D748}$
is to extract components
$\unicode[STIX]{x1D70E}_{rr}$
,
$\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D717}r}$
and
$\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D717}\unicode[STIX]{x1D717}}$
from
$\unicode[STIX]{x1D748}$
, compute their gradients and combine them with the dyads and the gradients of the dyads formed by the spherical coordinate unit vectors. Another option is to apply the formulae for the components in spherical coordinates of the gradient of a rank-two tensor field given in Torrilhon (Reference Torrilhon2010). Once
$\unicode[STIX]{x1D735}\text{q}$
and
$\unicode[STIX]{x1D735}\unicode[STIX]{x1D748}$
has been determined, their corresponding symmetric and trace-free tensors are computed from the formulae in appendix C. The expressions for the components of
$\text{R}$
and
$\text{m}$
needed in the boundary conditions are

*a*) $$\begin{eqnarray}\displaystyle \text{R}_{rr} & = & \displaystyle 108Kn^{2}B_{1}Gr^{-4}\cos \unicode[STIX]{x1D717}+{\textstyle \frac{48}{5}}KnF_{2}Gr^{-4}\nonumber\\ \displaystyle & & \displaystyle \times \exp [-\unicode[STIX]{x1D706}_{2}(r-1)](\unicode[STIX]{x1D706}_{2}^{2}r^{2}+3\unicode[STIX]{x1D706}_{2}r+3)\cos \unicode[STIX]{x1D717},\end{eqnarray}$$

*b*) $$\begin{eqnarray}\displaystyle \text{R}_{\unicode[STIX]{x1D717}r} & = & \displaystyle 54Kn^{2}B_{1}Gr^{-4}\sin \unicode[STIX]{x1D717}+\frac{12}{5}KnF_{2}Gr^{-4}\nonumber\\ \displaystyle & & \displaystyle \times \exp [-\unicode[STIX]{x1D706}_{2}(r-1)](\unicode[STIX]{x1D706}_{2}^{3}r^{3}+3\unicode[STIX]{x1D706}_{2}^{2}r^{2}+6\unicode[STIX]{x1D706}_{2}r+6)\sin \unicode[STIX]{x1D717},\end{eqnarray}$$

*c*) $$\begin{eqnarray}\displaystyle \text{m}_{rrr} & = & \displaystyle -{\textstyle \frac{48}{5}}KnGr^{-5}(10\unicode[STIX]{x1D706}_{3}^{-2}A_{1}+4Kn^{2}A_{1}-r^{2}A_{1}-10KnC_{2})\cos \unicode[STIX]{x1D717}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{16}{25}\frac{\unicode[STIX]{x1D706}_{3}^{2}}{\unicode[STIX]{x1D706}_{1}^{2}-\unicode[STIX]{x1D706}_{3}^{2}}Kn^{3}E_{1}Gr^{-5}\exp [-\unicode[STIX]{x1D706}_{1}(r-1)]\nonumber\\ \displaystyle & & \displaystyle \times \,(\unicode[STIX]{x1D706}_{1}^{4}r^{4}+7\unicode[STIX]{x1D706}_{1}^{3}r^{3}+27\unicode[STIX]{x1D706}_{1}^{2}r^{2}+60\unicode[STIX]{x1D706}_{1}r+60)\cos \unicode[STIX]{x1D717},\end{eqnarray}$$

*d*) $$\begin{eqnarray}\displaystyle \text{m}_{\unicode[STIX]{x1D717}rr} & = & \displaystyle -{\textstyle \frac{8}{5}}KnGr^{-5}(30\unicode[STIX]{x1D706}_{3}^{-2}A_{1}+12Kn^{2}A_{1}-r^{2}A_{1}-30KnC_{2})\nonumber\\ \displaystyle & & \displaystyle \times \sin \unicode[STIX]{x1D717}+\frac{32}{25}\frac{\unicode[STIX]{x1D706}_{3}^{2}}{\unicode[STIX]{x1D706}_{1}^{2}-\unicode[STIX]{x1D706}_{3}^{2}}Kn^{3}E_{1}Gr^{-5}\nonumber\\ \displaystyle & & \displaystyle \times \exp [-\unicode[STIX]{x1D706}_{1}(r-1)](\unicode[STIX]{x1D706}_{1}^{3}r^{3}+6\unicode[STIX]{x1D706}_{1}^{2}r^{2}+15\unicode[STIX]{x1D706}_{1}r+15)\sin \unicode[STIX]{x1D717}.\end{eqnarray}$$

Regarding the boundary conditions, from the solutions we also find that $\text{R}_{\unicode[STIX]{x1D717}\unicode[STIX]{x1D717}}=\text{R}_{\unicode[STIX]{x1D719}\unicode[STIX]{x1D719}}$ , and $\text{m}_{\unicode[STIX]{x1D717}\unicode[STIX]{x1D717}r}=\text{m}_{\unicode[STIX]{x1D719}\unicode[STIX]{x1D719}r}$ . Because of this, and referring to the discussion at the end of § 2.1, we have that by enforcing (2.11), constraint (2.12) can be dropped entirely. Moreover, equations (2.9) and (2.13) provide constraints for $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D717}r}$ and $\text{R}_{\unicode[STIX]{x1D717}r}$ whereas they are trivially satisfied for $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D719}r}$ and $\text{R}_{\unicode[STIX]{x1D719}r}$ . Therefore, including (2.10), (2.14) and (2.15), in total, we are left with six scalar boundary conditions.

Substitution of the solutions for the various fields into these boundary conditions leads to a system of six linear algebraic equations for the six unknowns $A_{1}$ , $B_{1}$ , $C_{2}$ , $D_{1}$ , $E_{1}$ and $F_{2}$ . For the problem of a sphere exposed to a uniform temperature gradient in the far field, this system may be written as

*a*) $$\begin{eqnarray}\displaystyle & & \displaystyle \hspace{-10.0pt}-18Kn(256\sqrt{2}Kn^{4}+64\sqrt{\unicode[STIX]{x03C0}}Kn^{3}-8\sqrt{2}Kn^{2}+5\sqrt{2})A_{1}-135\sqrt{2}Kn^{3}B_{1}\nonumber\\ \displaystyle & & \displaystyle \hspace{-10.0pt}\quad +\,180Kn^{2}(24\sqrt{2}Kn^{2}+6\sqrt{\unicode[STIX]{x03C0}}Kn+\sqrt{2})C_{2}\nonumber\\ \displaystyle & & \displaystyle \hspace{-10.0pt}\quad -\,18Kn^{2} [\!216\sqrt{2}Kn^{3}+18(4\sqrt{15}+3\sqrt{\unicode[STIX]{x03C0}})Kn^{2}\nonumber\\ \displaystyle & & \displaystyle \hspace{-10.0pt}\quad +\,9\sqrt{2}(\sqrt{15\unicode[STIX]{x03C0}}+8)Kn+4\sqrt{15}+15\sqrt{\unicode[STIX]{x03C0}}\!]E_{1}\nonumber\\ \displaystyle & & \displaystyle \hspace{-10.0pt}\quad -\,4\sqrt{2}(9Kn^{2}+3\sqrt{5}Kn+5)F_{2}=-135\sqrt{2}Kn^{3},\end{eqnarray}$$

*b*) $$\begin{eqnarray}\displaystyle & & \displaystyle \hspace{-10.0pt}-84\sqrt{2}Kn(32Kn^{2}-9)A_{1}+30Kn(270\sqrt{2}Kn^{2}+105\sqrt{\unicode[STIX]{x03C0}}Kn+28\sqrt{2})B_{1}\nonumber\\ \displaystyle & & \displaystyle \hspace{-10.0pt}\quad +\,2520\sqrt{2}Kn^{2}C_{2}-840\sqrt{2}KnD_{1}-42(54\sqrt{2}Kn^{3}+18\sqrt{15}Kn^{2}+12\sqrt{2}Kn-\sqrt{15})E_{1}\nonumber\\ \displaystyle & & \displaystyle \hspace{-10.0pt}\quad +\,40[54\sqrt{2}Kn^{2}+3(6\sqrt{10}+7\sqrt{\unicode[STIX]{x03C0}})Kn+10\sqrt{2}+7\sqrt{5\unicode[STIX]{x03C0}}]F_{2}=1575\sqrt{\unicode[STIX]{x03C0}}Kn^{2},\end{eqnarray}$$

*c*) $$\begin{eqnarray}\displaystyle & & \displaystyle \hspace{-10.0pt}-42Kn(1280\sqrt{\unicode[STIX]{x03C0}}Kn^{3}+224\sqrt{2}Kn^{2}-120\sqrt{\unicode[STIX]{x03C0}}Kn-33\sqrt{2})A_{1}+30\sqrt{2}Kn(135Kn^{2}-7)B_{1}\nonumber\\ \displaystyle & & \displaystyle \hspace{-10.0pt}\quad +1260Kn^{2}(40\sqrt{\unicode[STIX]{x03C0}}Kn+7\sqrt{2})C_{2}+210\sqrt{2}KnD_{1}\nonumber\\ \displaystyle & & \displaystyle \hspace{-10.0pt}\quad -\,21[\!2160\sqrt{\unicode[STIX]{x03C0}}Kn^{4}+18\sqrt{2}(20\sqrt{15\unicode[STIX]{x03C0}}+21)Kn^{3}\nonumber\\ \displaystyle & & \displaystyle \hspace{-10.0pt}\quad +\,18(7\sqrt{15}+45\sqrt{\unicode[STIX]{x03C0}})Kn^{2}+\sqrt{2}(35\sqrt{15\unicode[STIX]{x03C0}}+144)Kn+13\sqrt{15}+25\sqrt{\unicode[STIX]{x03C0}}\!]E_{1}\nonumber\\ \displaystyle & & \displaystyle \hspace{-10.0pt}\quad +\,40\sqrt{2}(27Kn^{2}+9\sqrt{5}Kn+5)F_{2}=0,\end{eqnarray}$$

*d*) $$\begin{eqnarray}\displaystyle & & \displaystyle -18\sqrt{2}Kn(256Kn^{4}-8Kn^{2}-5)A_{1}+135Kn^{3}(72\sqrt{\unicode[STIX]{x03C0}}Kn+13\sqrt{2})B_{1}\nonumber\\ \displaystyle & & \displaystyle \hspace{-10.0pt}\quad +\,180\sqrt{2}Kn^{2}(24Kn^{2}-1)C_{2}\nonumber\\ \displaystyle & & \displaystyle \hspace{-10.0pt}\quad -\,72Kn^{2}(54\sqrt{2}Kn^{3}+18\sqrt{15}Kn^{2}+18\sqrt{2}Kn+\sqrt{15})E_{1}\nonumber\\ \displaystyle & & \displaystyle \hspace{-10.0pt}\quad +\,4[\!648\sqrt{\unicode[STIX]{x03C0}}Kn^{3}+9(13\sqrt{2}+24\sqrt{5\unicode[STIX]{x03C0}})Kn^{2}\nonumber\\ \displaystyle & & \displaystyle \hspace{-10.0pt}\quad +\,3(13\sqrt{10}+60\sqrt{\unicode[STIX]{x03C0}})Kn+65\sqrt{2}+20\sqrt{5\unicode[STIX]{x03C0}}\!]F_{2}=-1485\sqrt{2}Kn^{3},\end{eqnarray}$$

*e*) $$\begin{eqnarray}\displaystyle & & \displaystyle \hspace{-10.0pt}A_{1}-3Kn^{2}B_{1}+2KnC_{2}-{\textstyle \frac{4}{15}}(3Kn+\sqrt{5})F_{2}=0,\end{eqnarray}$$

*f*) $$\begin{eqnarray}\displaystyle & & \displaystyle \hspace{-10.0pt}{\textstyle \frac{15}{2}}Kn^{2}B_{1}+{\textstyle \frac{15}{4}}\unicode[STIX]{x1D6EC}Kn^{2}D_{1}+{\textstyle \frac{2}{3}}(3Kn+\sqrt{5})F_{2}=-{\textstyle \frac{15}{4}}(\unicode[STIX]{x1D6EC}-1)Kn^{2}.\end{eqnarray}$$

For the velocity problem, the left-hand sides of expressions (3.25) remain the same, whereas the right-hand sides for (3.25*a*
), (3.25*b*
), (3.25*d*
), (3.25*e*
) and (3.25*f*
) become, respectively,
$180\sqrt{2}Kn^{2}$
, 0,
$-180\sqrt{2}Kn^{2}$
,
$-Kn$
and 0. The solution of this linear system of equations is found by means of a computer symbolic algebra package. The expressions for the coefficients are exceedingly large and are not reproduced here. They depend on the parameters
$Kn$
and
$\unicode[STIX]{x1D6EC}$
. Once these coefficients are determined, the field variables can be computed at any position
$(r,\unicode[STIX]{x1D717})$
, for given values of parameters
$Kn$
,
$\unicode[STIX]{x1D6EC}$
and
$Ep$
(or
$Ma$
), using (3.23).

In general terms, the main attributes of the solution method proposed here are twofold. First, it replaces an *a priori* knowledge of the exact dependency of the various scalar fields, including spherical-coordinate components of vectors and tensors, on the polar angle – e.g. from the classical solution of Stokes flow past a sphere – with the rather more general constraint of a dependency of the various tensor fields (scalar, vectors and rank-two tensors) on the external perturbation vector
$\text{G}$
multiplied by a few multipole potentials appropriately chosen. Second, instead of dealing with the components of the Laplacian of the stress tensor in spherical coordinates and with the scalar differential equations carrying them, we solve a second-order differential equation for the stress in compact, tensor form.

## 4 Results and discussion

This section contains three parts. First, we discuss predictions from our solution of the R13 moment equations for the problem of thermophoresis of a spherical particle presented in § 2.1. In particular, we show flow and temperature profiles for the sphere’s surroundings contrasted with results from kinetic theory and then compare predictions for the thermophoretic force from our solution to those from other theories as well as recently published experimental measurements. Second, we present results from R13 and other theories for the problem described in § 2.2, that is, the drag force due to a uniform flow past a sphere taking into account the particle-to-gas heat conductivity ratio. In both sub-sections, we present profiles of velocity components, density and temperature in a neighbourhood of the spherical surface. Third, the solutions from these two problems are coupled in the balance of the thermophoretic force by the drag acting on the sphere when this is free to move giving rise to the sphere’s thermophoretic terminal velocity. We show predictions from R13 and other theories for this quantity.

Note that the results presented here for the temperature deviations in the gas and solid, $\unicode[STIX]{x1D703}$ and $\unicode[STIX]{x1D703}_{s}$ and the gas heat flux, $\text{q}$ , correspond to the variables on the left-hand sides of expressions in (2.8). These are obtained after transforming back from the fields given in (3.23), which represent the variables with ‘ $\check{\,}$ ’.

### 4.1 Thermophoresis on a heat-conducting sphere

#### 4.1.1 Flow and temperature fields

Examining the response of some of the variables modelled with the R13 equations to changes in the Knudsen number
$Kn$
and solid-to-gas conductivity ratio
$\unicode[STIX]{x1D6EC}$
near the surface of the sphere is of interest. In figure 2, we show the profiles of the velocity components, density and temperature deviations as functions of the radial coordinate starting at the surface of the sphere,
$r=1$
, obtained with the R13 exact solution derived here, the numerical solutions of the linearized Boltzmann equation by Sone (Reference Sone2007) for a hard-sphere gas, and his asymptotic expression for small Knudsen number. Curves are presented for an isothermal sphere, i.e.
$\unicode[STIX]{x1D6EC}\rightarrow \infty$
, and
$Kn=0$
, 0.090, 0.180, 0.269 and 0.539, corresponding to Sone’s
$k=0$
, 0.1, 0.2, 0.3 and 0.6, respectively – the relationship between
$Kn$
and
$k$
is given in the figure’s caption. Predictions from the theories depict typical rarefaction effects, namely, velocity slip (figure 2
*b*), and temperature jump (figure 2
*d*) at the sphere’s surface. When
$k\rightarrow 0$
, and
$k=0.1$
, R13 agrees well with the results from kinetic theory; however, as
$k$
increases the differences between the two models becomes noticeable, suggesting that R13’s quantitative description of the Knudsen layer near the surface of the sphere begins to become incomplete for
$Kn\gtrsim 0.1$
.

Contour plots of gas speed (normalized by $Ep$ ) and velocity streamlines are presented in figure 3 for combinations of $Kn=0.02$ and 0.2 and $\unicode[STIX]{x1D6EC}=4$ and $\unicode[STIX]{x1D6EC}\rightarrow \infty$ . The temperature gradient points from left to right. Note that except for the case of $Kn=0.02$ and $\unicode[STIX]{x1D6EC}\rightarrow \infty$ , the gas flow near the sphere is in the direction of the temperature gradient. Following the reasoning by Sone (Reference Sone2007), this flow is driven by the force exerted by the solid surface on the gas in this direction. This force is the reaction to the momentum transferred onto the sphere’s surface by the gas in the opposite direction, that is, from the hot to the cold region. This corresponds to a scenario of normal or ‘positive’ thermophoresis. For the combination $Kn=0.02$ and $\unicode[STIX]{x1D6EC}\rightarrow \infty$ , on the other hand, the situation is inverted, and the gas flow is from hot to cold. This case is thus known as reversed or ‘negative’ thermophoresis. We shall discuss this phenomenon below with reference to the macroscopic transport equations of the previous section and, in particular, the boundary condition for slip. This is followed by a quantitative investigation of the net thermophoretic force acting on the sphere with special attention to its change in direction. Another feature shown in figure 3 is that for $Kn=0.02$ , the maximum speed in the case $\unicode[STIX]{x1D6EC}\rightarrow \infty$ is approximately one order of magnitude smaller than in the case $\unicode[STIX]{x1D6EC}=4$ , and nearly two orders of magnitude smaller than the other two cases corresponding to $Kn=0.2$ . One may thus expect that, among the four cases considered, the smallest net force on the sphere should be attained in the case $Kn=0.02$ and $\unicode[STIX]{x1D6EC}\rightarrow \infty$ .

Contour plots of gas temperature deviation (normalized by $Ep$ ) with streamlines for the heat-flux vector in the gas and within the solid sphere are shown in figure 4 for four cases with $Kn=0.02$ and 0.2 and $\unicode[STIX]{x1D6EC}=4$ and $\unicode[STIX]{x1D6EC}\rightarrow \infty$ . The heat-flux vectors point to the left, opposite to the temperature gradient. In the cases where $\unicode[STIX]{x1D6EC}\rightarrow \infty$ , figures show no temperature deviation inside the sphere. In contrast, for $\unicode[STIX]{x1D6EC}=4$ , a non-zero temperature gradient is observed in the spherical particle, with the temperature varying linearly with the axial coordinate. In the two cases with $\unicode[STIX]{x1D6EC}=4$ and also in the case where $\unicode[STIX]{x1D6EC}\rightarrow \infty$ and $Kn=0.2$ , temperature jumps across the spherical surface are evident. Enhancing the gas rarefaction by increasing $Kn$ for fixed $\unicode[STIX]{x1D6EC}$ increases the temperature jump across the spherical surface. For instance, with $\unicode[STIX]{x1D6EC}=4$ , there are more contour lines inside the sphere for $Kn=0.02$ than for 0.2, pointing toward a smoother temperature gradient in the latter, whereas the opposite takes place in the gas, as the isothermal lines are more bent for $Kn=0.02$ than for 0.2. This yields a greater temperature jump at a given point of the sphere’s surface ( $z\neq 0$ ) in the case with higher $Kn$ .

Another notable feature in figure 4 is that lines of constant temperature intercept the surface of the sphere at various points in the two cases where
$\unicode[STIX]{x1D6EC}=4$
as well as for
$Kn=0.2$
and
$\unicode[STIX]{x1D6EC}\rightarrow \infty$
, indicating the presence of a temperature gradient in the tangential direction. Under rarefied conditions, it is known that this temperature gradient induces gas motion along its direction (i.e. from cold to hot) near the solid surface, an effect known as thermal creep or thermal transpiration (Maxwell Reference Maxwell1879; Kennard Reference Kennard1938; Sone Reference Sone2007; Mohammadzadeh *et al.*
Reference Mohammadzadeh, Rana and Struchtrup2015). This is the type of flow depicted by the streamlines in the corresponding plots of figure 3. On the other hand, for
$Kn=0.02$
and
$\unicode[STIX]{x1D6EC}\rightarrow \infty$
, the isothermal lines in the gas tend to wrap the sphere’s surface resulting in temperature gradients that essentially vanish in the tangential direction on the gas side of the sphere’s surface. This points toward a hindering of the thermal creep in comparison with the other cases, in accord with the small velocity magnitudes and, more importantly, the reversal in the direction of the streamlines shown for the same case in figure 3. In fact, this reversed flow, now from the hot to the cold region, occurs because another type of flow, the so-called thermal-stress slip flow (Sone Reference Sone2007; Young Reference Young2011), becomes dominant over the thermal creep. The thermal-stress slip flow represents a Knudsen number higher-order effect and, even though it is also induced by changes in the gas temperature distribution, is of a different nature from the thermal creep.

From his asymptotic analysis of the linearized Boltzmann equation for a gas that deviates slightly from a state of uniform equilibrium at rest, Sone (Reference Sone2007) identified that the slip flow in the slip boundary condition was determined by the term proportional to $-\text{t}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}\boldsymbol{\cdot }\text{n}$ multiplied by a positive constant when the boundary surface has a uniform temperature or in the absence of thermal creep (here $\text{t}$ and $\text{n}$ are unit vectors tangential and normal to the bounding surface, respectively). Because $\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}$ multiplied by a constant is one of the terms in his expression for the stress tensor, he designated this flow as thermal-stress slip flow. As discussed by Sone (Reference Sone2007), if the boundary temperature is constant then $-\text{t}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}\boldsymbol{\cdot }\text{n}=-\text{t}\boldsymbol{\cdot }\unicode[STIX]{x1D735}(\text{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D703})=-\unicode[STIX]{x2202}(\unicode[STIX]{x2202}\unicode[STIX]{x1D703}/\unicode[STIX]{x2202}n)/\unicode[STIX]{x2202}s$ , where $s$ is the arc length. Then, if the isothermal surfaces in the gas are not parallel to the solid boundary (also of constant temperature), i.e. the component of the temperature gradient normal to the boundary ( $\unicode[STIX]{x2202}\unicode[STIX]{x1D703}/\unicode[STIX]{x2202}n$ ) changes along it, when the boundary temperature is higher (lower) than that in the gas, a flow is promoted in the direction in which the isothermal surfaces converge (diverge). This is exemplified by the contour plots in figures 3 and 4 for $Kn=0.02$ and $\unicode[STIX]{x1D6EC}\rightarrow \infty$ , and figure 5 sketches this behaviour for clarity. The thermal-stress slip flow is the primary cause of a negative thermophoretic force.

Recalling the generalized slip boundary condition (2.9), rewritten here for convenience in a slightly different form

setting
$u_{s,\unicode[STIX]{x1D717}}=0$
and using the expressions from the analytical solution obtained with the R13-moment equations, we investigate the contribution of each of the terms on the right-hand side of (4.1) to the gas slip velocity
$u_{\unicode[STIX]{x1D717}}$
on the surface of the sphere,
$r=1$
, and selected
$\unicode[STIX]{x1D717}$
values. Considering for this exercise the conditions of figures 3 and 4, we have that for the cases of
$\unicode[STIX]{x1D6EC}=4$
and also for
$\unicode[STIX]{x1D6EC}\rightarrow \infty$
and
$Kn=0.2$
, where thermal creep is predominant, the greater contribution in absolute value is from the first term on the right-hand side, giving a negative value (flow from cold to hot), whereas for the remaining case of
$\unicode[STIX]{x1D6EC}\rightarrow \infty$
and
$Kn=0.02$
, exhibiting thermal-stress slip flow, the larger value is positive and results from the shear-stress term
$-(\unicode[STIX]{x03C0}/2)^{1/2}\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D717}r}$
. In all cases, the last term in (4.1), a Knudsen number higher-order term, results in a significantly smaller absolute value in comparison. Using now the balance equation for the heat flux (2.6*b*
) for
$q_{\unicode[STIX]{x1D717}}$
when the flow has the direction of the temperature gradient or the balance equation for the deviatoric stress (2.6*a*
) for
$\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D717}r}$
when the flow has the opposite direction – after moving the terms with derivatives of
$\text{q}$
and
$\unicode[STIX]{x1D748}$
, of higher order in
$Kn$
, to the right-hand side – we found that in the former, the term
$(3/4)Kn\,\unicode[STIX]{x2202}\unicode[STIX]{x1D703}/\unicode[STIX]{x2202}\unicode[STIX]{x1D717}$
carries the largest weight whereas in the latter, the term proportional to
$Kn(\overline{\unicode[STIX]{x1D735}\text{ q}})_{\unicode[STIX]{x1D717}r}$
is predominant. Note that this term and not the shear stress from classical hydrodynamics,
$-2Kn(\overline{\unicode[STIX]{x1D735}\text{ u}})_{\unicode[STIX]{x1D717}r}$
, is the leading contributor in this case. Examination of the term proportional to
$Kn(\overline{\unicode[STIX]{x1D735}\text{ q}})_{\unicode[STIX]{x1D717}r}$
by substitution of (2.6*b*
) for
$\text{q}$
when the slip flow is from hot to cold revealed that, as expected, the major role in this situation is played by the term
$-3(\unicode[STIX]{x03C0}/2)^{1/2}Kn^{2}(\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}\unicode[STIX]{x1D703})_{\unicode[STIX]{x1D717}r}$
(
${\approx}-3(\unicode[STIX]{x03C0}/2)^{1/2}Kn^{2}\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D703}/\unicode[STIX]{x2202}r\unicode[STIX]{x2202}\unicode[STIX]{x1D717}$
when
$\unicode[STIX]{x2202}\unicode[STIX]{x1D703}/\unicode[STIX]{x2202}\unicode[STIX]{x1D717}\approx 0$
) in agreement with the argument of Sone (Reference Sone2007) for the thermal-stress slip flow and the illustration in figure 5. From a microscopic perspective, Sone (Reference Sone2007) emphasizes that for either thermal creep or thermal-stress slip flow, the gas motion is induced by the difference in the velocity distribution functions of the molecules colliding with the solid boundary and those leaving it. The different nature of the flow and force between thermal creep and thermal-stress slip flow is affected by the fact that, before impinging onto the solid boundary, particles starting within a mean free path essentially keep the attributes of their origins, and the temperature distribution in the gas surrounding the boundary is notably different in one case or the other.

#### 4.1.2 Thermophoretic force

The thermophoretic force acting on the sphere is obtained by integrating over the surface of the sphere the projection of the total stress vector at $r^{\ast }=a^{\ast }$ onto direction $\text{k}$ , i.e. $(-p^{\ast }\hat{\text{r}}-\unicode[STIX]{x1D748}^{\ast }\boldsymbol{\cdot }\hat{\text{r}})\boldsymbol{\cdot }\text{k}$ . Young (Reference Young2011) introduced the dimensionless thermophoretic force

where $F_{T}^{\ast }$ denotes the dimensional thermophoretic force. From our solution of the R13 equations for the pressure and the deviatoric stress, we obtain the exact expression

and coefficients $\unicode[STIX]{x1D6FC}_{m}^{(0)}$ , $\unicode[STIX]{x1D6FC}_{m}^{(1)}$ , $\unicode[STIX]{x1D6FD}_{m}^{(0)}$ and $\unicode[STIX]{x1D6FD}_{m}^{(1)}$ are given in table 1. For the particular case of an isothermal sphere, we simply compute the limit $\unicode[STIX]{x1D6EC}\rightarrow \infty$ in this expression, resulting in $\unicode[STIX]{x1D6F7}=-12\unicode[STIX]{x03C0}\sum \unicode[STIX]{x1D6FC}_{m}^{(1)}Kn^{m}/\sum \unicode[STIX]{x1D6FD}_{m}^{(1)}Kn^{m}$ , with index $m$ spanning the same ranges as in (4.3).

Results from (4.3) are plotted as $-\unicode[STIX]{x1D6F7}/(2\unicode[STIX]{x03C0})$ versus $(\unicode[STIX]{x03C0}/2)^{1/2}Kn$ for $\unicode[STIX]{x1D6EC}=4$ , 10, and $22.4\times 10^{3}$ in figure 6. The highest $\unicode[STIX]{x1D6EC}$ value is chosen motivated by one of the experimental data sets depicted in the figure (see below). We also include predictions from various models based on the linearized Boltzmann equation, namely, by Sone & Aoki (Reference Sone and Aoki1983) for an isothermal sphere using the BGK equation – denoted in their work as the Boltzmann–Krook–Welander (BKW) equation – by Beresnev & Chernyak (Reference Beresnev and Chernyak1995) using the S model, and by Sone (Reference Sone2007) for a hard-sphere gas. Results from G13 and from its modification represented by Young’s (Reference Young2011) interpolation formula are also added. In addition, values from Waldmann (Reference Waldmann1959) formula, valid for large $Kn$ and insensitive to changes in $\unicode[STIX]{x1D6EC}$ , are presented. This formula as well as the expressions for the dimensionless thermophoretic force from G13 and Young (Reference Young2011) are presented, for completeness, in appendix D.

In figure 6 we included the experimental results on thermophoresis recently presented by Bosworth & Ketsdever (Reference Bosworth and Ketsdever2016) for spheres made of the polymer acrylonitrile butadiene styrene (ABS) and by Bosworth *et al.* (Reference Bosworth, Ventura, Ketsdever and Gimelshein2016) for spheres of copper. In the latter, data include observations of reversed or negative thermophoresis. In the experiments, they set the sphere in the mid-plane between two copper plates and placed the entire assembly in a vacuum chamber filled with argon. Rarefied conditions were attained by reducing the pressure to values significantly below the atmospheric pressure. The force was generated by keeping one plate at the ambient temperature while the other plate was heated uniformly to a higher temperature. To plot the data in figure 6 – originally presented as dimensional force versus pressure – we used an ambient temperature of 24.5
$^{\circ }$
C (Bosworth, R. & Ketsdever, A. 2017 Private communication.) and a temperature difference between plates of
$35~~\text{K}~\text{m}^{-1}$
, a distance between plates of 0.40 m, a particle radius of 0.0254 m, a gas specific constant of
$208~\text{J}~\text{kg}^{-1}\,\text{K}^{-1}$
and a gas viscosity of
$2.295\times 10^{-5}~\text{Pa}~\text{s}$
at the mid-plane temperature (NIST 2017). The solid-to-gas heat conductivity ratios were estimated to be
$\unicode[STIX]{x1D6EC}=10$
and
$22.4\times 10^{3}$
for the ABS and copper spheres, respectively. In practice, in the copper case, such a large value would correspond to an isothermal sphere (
$\unicode[STIX]{x1D6EC}\rightarrow \infty$
). To account for Knudsen number effects in the temperature gradient by the mid-plane between the plates, we adopted the model in chapter 4 of Sone (Reference Sone2007) based on the linearized Boltzmann equation for a hard-sphere gas. These gas rarefaction effects include a lower temperature gradient in the mid-plane with respect to the nominal value, gas–wall temperature jumps at the plates’ surfaces and noticeable deviations from linearity in the temperature profile in regions near the walls (Knudsen layers). These effects are evident, for instance, in figure 5 in Bosworth *et al.* (Reference Bosworth, Ventura, Ketsdever and Gimelshein2016), obtained using DSMC. This corrected value of the temperature gradient was used to compute the Epstein number
$Ep$
for each experimental point. We obtain values of
$Ep$
in the order of
$3.0\times 10^{-3}$
for the experimental data sets.

For $\sqrt{\unicode[STIX]{x03C0}/2}\,Kn\lesssim 0.1$ , figure 6 shows that the thermophoretic force from R13 tends very closely to the results from models of the linearized Boltzmann equation. In contrast, G13 significantly under-predicts these results. As demonstrated in the figure, this discrepancy is corrected in Young’s (Reference Young2011) interpolation formula, obtained from G13 after fitting the thermal creep, velocity slip and temperature jump coefficients from the so-called Maxwell–Smoluchowski values (e.g. see Kennard Reference Kennard1938; Nguyen & Wereley Reference Nguyen and Wereley2002; Gu & Emerson Reference Gu and Emerson2007; Young Reference Young2011) with those resulting from solutions of model Boltzmann equations compiled by Sharipov (Reference Sharipov2004). Young argues that changing these coefficients in the boundary conditions of G13 is needed to compensate for the inability of this model to reproduce the Knudsen layer. On the other hand, our analytical result from R13 gives the correct prediction without changing the original numerical factors appearing in the various terms of both the boundary conditions and bulk equations, i.e. with no fitting. Whereas the new set of coefficients adopted by Young (Reference Young2011) for G13 works well for the thermophoretic force on a sphere for very small $Kn$ , they might not be applicable to situations other than for which they were fitted (e.g. for different geometries).

When $Kn\rightarrow 0$ we have from (4.3) for R13 that

whereas in the special case $\unicode[STIX]{x1D6EC}\rightarrow \infty$ (isothermal sphere), we find the following asymptotic result when $Kn\rightarrow 0$ ,

Note from these expressions that for very small $Kn$ and finite $\unicode[STIX]{x1D6EC}$ , $-\unicode[STIX]{x1D6F7}$ tends to a positive constant, indicating positive or normal thermophoresis. On the other hand, after passing the limit $\unicode[STIX]{x1D6EC}\rightarrow \infty$ , the asymptotic result shows that $-\unicode[STIX]{x1D6F7}$ tends to zero linearly with $Kn$ through negative values corresponding to slightly reversed thermophoresis. A similar trend can be observed in the asymptotic formulae by Sone (Reference Sone2007).

Within the interval
$0.02\lesssim \sqrt{\unicode[STIX]{x03C0}/2}\,Kn\lesssim 0.2$
, figure 6(*c*) shows that the experimental data for the copper sphere exhibits negative or reversed thermophoresis,
$-\unicode[STIX]{x1D6F7}<0$
. The reversed net force direction, such that it is now from the cold to the hot region for large values of
$\unicode[STIX]{x1D6EC}$
, is mainly caused by the thermal-stress flow near the sphere’s surface. In this figure, R13, G13, Young’s formula and the models by Sone & Aoki (Reference Sone and Aoki1983) and Beresnev & Chernyak (Reference Beresnev and Chernyak1995) predict negative thermophoresis. In particular, the contour plots for
$Kn=0.02$
and
$\unicode[STIX]{x1D6EC}\rightarrow \infty$
in figures 3 and 4 show the various features that characterized a scenario of negative thermophoresis, such as the reversed flow direction and the absence of temperature changes in the gas side along the solid surface, in comparison to the other cases considered. Nonetheless, whilst R13 agrees well with the benchmark Boltzmann solutions (for hard spheres and the S-model), both predict force magnitudes significantly smaller than the experimental values; a discrepancy which at present we cannot explain. Referring again to the values of
$Kn$
and
$\unicode[STIX]{x1D6EC}$
considered in figure 3, the force magnitude in the case
$Kn=0.02$
and the largest
$\unicode[STIX]{x1D6EC}$
, according to figure 6, is notably smaller than for the other cases, which seems to correspond with the fact that the velocity magnitude in this case is also one or more orders of magnitude smaller than in the others.

Expressions of thermophoretic force for R13 and G13, as well as Young’s interpolation formula exhibit a critical $\unicode[STIX]{x1D6EC}$ above which reversed thermophoresis occurs for some bounded $Kn$ interval. For R13, this threshold is $\unicode[STIX]{x1D6EC}=87.17$ , whereas for G13 it is a much lower ratio, perhaps less likely to actually occur, of 10.93. In the case of Young’s interpolation formula the threshold is $\unicode[STIX]{x1D6EC}=26.05$ . For larger, but finite values of $\unicode[STIX]{x1D6EC}$ , the corresponding curve of $-\unicode[STIX]{x1D6F7}/(2\unicode[STIX]{x03C0})$ will intercept the $\unicode[STIX]{x1D6F7}=0$ axis at two points, and in both $Kn>0$ . The maximum negative value of $-\unicode[STIX]{x1D6F7}/(2\unicode[STIX]{x03C0})$ occurs when $\unicode[STIX]{x1D6EC}\rightarrow \infty$ . From the R13 model, this is $-\unicode[STIX]{x1D6F7}/(2\unicode[STIX]{x03C0})=-0.068$ at $Kn=0.024$ ; for G13 it is $-\unicode[STIX]{x1D6F7}/(2\unicode[STIX]{x03C0})=-0.321$ at $Kn=0.087$ , and for Young’s formula it is $-\unicode[STIX]{x1D6F7}/(2\unicode[STIX]{x03C0})=-0.180$ at $Kn=0.049$ . Negative thermophoresis cannot occur for $Kn\geqslant 0.054$ for R13, $Kn\geqslant 0.251$ for G13 and $Kn\geqslant 0.130$ for Young’s formula; these critical values are attained with $\unicode[STIX]{x1D6EC}\rightarrow \infty$ . This upper limit of $Kn$ will decrease for $\unicode[STIX]{x1D6EC}<\infty$ . The critical magnitudes of Knudsen numbers listed in this paragraph were instrumental in choosing $Kn=0.02$ and 0.2 for the contour plots in figures 3 and 4.

For $0.2\lesssim \sqrt{\unicode[STIX]{x03C0}/2}\,Kn\lesssim 1$ , R13 shows qualitatively the best performance among the models derived with the moments method for both the copper and polymer (ABS) data. In particular, for $\sqrt{\unicode[STIX]{x03C0}/2}\,Kn~\lesssim ~0.4$ , its predictions agree with the experimental points for the copper’s sphere. It is likely coincidental that R13 performs better with respect to the experimental data than the benchmark models of Beresnev & Chernyak (Reference Beresnev and Chernyak1995) and Sone (Reference Sone2007) from kinetic theory. G13 depicts a rather unsatisfactory performance, as its maximum is significantly shifted to higher $Kn$ values. When $Kn>1$ , Waldmann’s formula matches the experiments, as expected. For larger $Kn$ values, the various model Boltzmann equations approximate well the experiments. G13 significantly over-predicts these values, although it tends to zero in the same asymptotic manner as Waldmann’s formula, that is, as $Kn^{-1}$ when $Kn\rightarrow \infty$ (see appendix D), but evidently with an incorrect coefficient. By construction, Young’s formula matches Waldmann’s for large $Kn$ . On the other hand, for R13, as $Kn\rightarrow \infty$ , we obtain from (4.3) that

and thus goes to zero faster than Waldmann’s.

For large $Kn$ ( $Kn>1$ , say), heat conductivity ratio $\unicode[STIX]{x1D6EC}$ plays a fairly inconsequential role. In this regime, according to figure 6, R13 under-predicts both the experimental force as well as Waldmann’s results. It is known that for such large values of $Kn$ , models such as G13 and R13 are not expected to give the correct quantitative result (Torrilhon Reference Torrilhon2010).

### 4.2 Uniform flow past a heat-conducting sphere

Profiles of the gas velocity components, density and temperature deviations as functions of the radial coordinate computed with the R13 exact solution of the previous section are presented in figure 7 for the case of a uniform flow past a stationary sphere with a thermal conductivity ratio
$\unicode[STIX]{x1D6EC}\rightarrow \infty$
(isothermal sphere). For
$Kn=0.05$
and 0.3, our solution is validated as its predictions match very well the exact solution of Torrilhon (Reference Torrilhon2010) obtained from his most recent Wolfram Mathematica code (see § 5.3 in Claydon *et al.*
Reference Claydon, Shrestha, Rana, Sprittles and Lockerby2017). Good agreement is also observed, in general, between our solution and results from the asymptotic expressions of Sone (Reference Sone2007) when
$Kn\rightarrow 0$
(
$k\rightarrow 0$
) and his numerical solution of a kinetic theory model for
$Kn=0.090$
(
$k=0.1$
).

The drag force acting on the sphere is computed from the same surface integral mentioned in the previous sub-section but using the pressure and deviatoric stress fields obtained from the solution of the problem of slow flow past a sphere. In dimensionless form, this drag force from R13 may be written as

where $F_{\text{Stokes}}^{\ast }$ denotes Stokes formula for drag on a sphere, $F_{\text{Stokes}}^{\ast }=6\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}_{0}^{\ast }U_{0}^{\ast }a^{\ast }$ and coefficients $\unicode[STIX]{x1D6FC}_{m}^{(0)}$ , $\unicode[STIX]{x1D6FC}_{m}^{(1)}$ , $\unicode[STIX]{x1D6FD}_{m}^{(0)}$ and $\unicode[STIX]{x1D6FD}_{m}^{(1)}$ are given in Table 2. For the special case of the drag on an isothermal sphere due to a streaming flow, expression (4.7), in the limit of $\unicode[STIX]{x1D6EC}\rightarrow \infty$ , reduces to $F_{D}^{\ast }/F_{\text{Stokes}}^{\ast }=(1+\sum \unicode[STIX]{x1D6FC}_{m}^{(1)}Kn^{\,m})/(1+\sum \unicode[STIX]{x1D6FD}_{m}^{(1)}Kn^{\,m})$ , where $m$ takes the same values as in (4.7). It should be noted that Torrilhon (Reference Torrilhon2010) showed the curve of drag versus Knudsen number resulting from his analysis, but did not present a closed-form expression for this force.

As $Kn\rightarrow 0$ , the drag can be computed from

Furthermore, when $Kn\rightarrow \infty$ , we have from (4.7) the expression

Figure 8 shows the dimensionless drag force as a function of
$Kn$
for
$\unicode[STIX]{x1D6EC}=4$
and
$\unicode[STIX]{x1D6EC}\rightarrow \infty$
. We added results extracted from the work by Torrilhon (Reference Torrilhon2010) for an isothermal sphere. The curves match each other very well. This clearly indicates that the drag force is insensitive to changes in the solid-to-gas thermal conductivity ratio, a tendency already noted by Sone (Reference Sone2007) from his analysis based on the linearized Boltzmann equation in chapter 4 of his monograph. For small and large
$Kn$
, this trend becomes evident in expressions (4.8) and (4.9), respectively. The matching with Torrilhon’s results serves as validation of the solution method adopted in our analysis. In addition, we included predictions from Young’s (Reference Young2011) solution for an isothermal sphere of the G13 equations using Maxwell–Smoluchowski values for the thermal creep, velocity slip and temperature jump coefficients, as well as with the coefficients extracted by Young from the work of Sharipov (Reference Sharipov2004) (see appendix D). For the interval shown, no significant difference is noted between these two sets. We included also predictions from kinetic theory by Sone (Reference Sone2007) for a hard-sphere gas and an isothermal sphere. Whereas G13 over-predicts the results by Sone, the R13 tendency is to under-predict them, although with a smaller difference. For
$Kn\lesssim 0.4$
, R13 approximates Sone’s curve very well. For each of the models, as in the case of R13, changing
$\unicode[STIX]{x1D6EC}$
from infinity to 4 produced no significant changes. Figure 8 also includes experimental data by Goldberg (Reference Goldberg1954) (extracted from Torrilhon, Reference Torrilhon2010) and by Allen & Raabe (Reference Allen and Raabe1982) (extracted from Claydon *et al.*, Reference Claydon, Shrestha, Rana, Sprittles and Lockerby2017). Predictions from both R13 and Sone (Reference Sone2007) under-predict the drag, with the former resulting in the largest difference.

### 4.3 Thermophoretic velocity

The thermophoretic velocity is the velocity of a moving particle driven by the thermophoretic force, when this force is balanced by the drag; that is, when the net force on the particle vanishes. Setting $F_{T}^{\ast }+F_{D}^{\ast }=0$ , and using expression (4.2) for $F_{T}^{\ast }$ , we readily have for the dimensionless velocity

expressed as the pseudo Mach number. Previously, $Ma$ was a prescribed parameter and now it is the unknown. Here, $\unicode[STIX]{x1D6F9}$ , a dimensionless quantity introduced by Young (Reference Young2011), is equal to $6\unicode[STIX]{x03C0}F_{D}^{\ast }/F_{\text{Stokes}}^{\ast }$ . In general, it is a function of $Kn$ and $\unicode[STIX]{x1D6EC}$ . The expression for $\unicode[STIX]{x1D6F9}$ from the G13 theory is given in appendix D.

With (4.3) and (4.7) for the non-dimensional thermophoretic force and drag, respectively, relation (4.10) provides an expression for the thermophoretic velocity according to R13. Figure 9 shows results from this expression normalized with factor $\sqrt{\unicode[STIX]{x03C0}/2}\,Ep\,Kn$ as a function of $\sqrt{\unicode[STIX]{x03C0}/2}\,Kn$ for different values of $\unicode[STIX]{x1D6EC}$ . This figure includes predictions from G13 and from solutions of model Boltzmann equations (Beresnev & Chernyak Reference Beresnev and Chernyak1995; Sone Reference Sone2007), and from the combination in (4.10) of Waldmann’s formula for $\unicode[STIX]{x1D6F7}$ and Epstein’s formula for $\unicode[STIX]{x1D6F9}$ intended for the free-molecule regime of large $Kn$ (see appendix D). The velocity from this model is independent of $\unicode[STIX]{x1D6EC}$ .

From figure 9 we note that R13 predicts the correct values given by the model Boltzmann equations for $\sqrt{\unicode[STIX]{x03C0}/2}\,Kn\lesssim 0.2$ . For greater values of $\sqrt{\unicode[STIX]{x03C0}/2}\,Kn$ , the agreement deteriorates. For $Kn\rightarrow \infty$ , R13 predicts a constant value, as is the case with the Waldmann–Epstein formula and with the results from kinetic theory, but with a smaller magnitude. G13, on the other hand, although also tending to a constant value, shows significant discrepancies with the other theories considered. Even though R13 is outside its limits of applicability for $Kn>1$ , it provides a reasonable engineering tool in that it has the correct asymptotic behaviour and, quantitatively, its differences with the benchmark models from kinetic theory and Waldmann–Epstein’s formula are not significantly large.

## 5 Concluding remarks

In this work, we investigated theoretically an instance of the phenomenon of thermophoresis, a term that refers to the forces on and motions of objects caused by temperature gradients when these objects are exposed to rarefied gases. In particular, we considered the problem of thermophoresis of a spherical particle, obtaining an analytical solution by solving the R13 moment equations, a model that provides a macroscopic description of rarefied gas flows up to the transition regime for Knudsen numbers smaller than one. Besides writing expressions for the field variables describing the fluid flow and heat transfer problem, by integration of the total stress on the surface of the sphere, we obtained a closed-form expression for the thermophoretic force as a function of the Knudsen number, dimensionless temperature gradient (Epstein number), and solid-to-gas heat conductivity ratio.

Employing the closed-form expressions for the field variables in the gas obtained here, we plotted their profiles in a neighbourhood of the sphere alongside predictions from a model Boltzmann equation from the literature. This comparison revealed that up to a Knudsen number of about $0.1$ , based on the sphere’s radius, the agreement between these solutions is very good thereby showing that the solution based on R13 is capable of fully describing the Knudsen layer for this particular physical situation in this range of Knudsen numbers. We also showed contour plots of speed and temperature including velocity and heat flux streamlines for various conditions. These figures exhibit the interplay between the various mechanisms involved in the complex process of thermophoresis, namely, the more common gas thermal creep toward the hot region caused by a temperature gradient in the gas parallel to the solid surface and, under very specific conditions in the slip regime that include a highly thermally conductive solid, the reversal of the flow direction resulting from the now dominant mechanism of slip flow driven by thermal stresses.

We extended the same modelling approach with the R13 equations to obtain an analytical solution for the drag by considering the problem of a heat-conducting spherical particle in a uniform, slow gas stream, thereby extending the analysis by Torrilhon (Reference Torrilhon2010) for an isothermal sphere. Predictions from the new solution are insensitive to changes in the solid-to-gas thermal conductivity ratio and agree very well with the results of Torrilhon, who used a different method to obtain his solution. For Knudsen numbers smaller than $0.4$ , R13 approximates reasonably well kinetic theory results for hard spheres. We then computed the thermophoretic velocity of the spherical particle when this can move driven by the thermophoretic force. The expression for the velocity results from balancing the thermophoretic force with the drag resistance exerted by the surrounding gas on the sphere. For Knudsen numbers smaller than $0.2$ , approximately, values for the thermophoretic velocity from R13 show good agreement with those from two models of the Boltzmann equation.

Results from the new thermophoretic force model derived with the R13 equations were compared with results from G13 and from various models based on the linearized Boltzmann equation, such as BGK and the S model, and with data from very recent experiments, for a wide range of Knudsen numbers. For Knudsen numbers below approximately $0.1$ , R13 results match predictions from the model Boltzmann equations. In this interval, the various theories considered, including R13, under-predict significantly the experimental measurements of reversed thermophoretic force for a metallic sphere (highly thermally conductive). For Knudsen numbers lying between $0.1$ and $1$ , approximately, the graphs from R13 follow qualitatively the experimental curves for both low and high solid-to-gas thermal conductivity ratios, and in general its predictions show, although important, the least differences with the experimental data. Surprisingly, solutions of models of the Boltzmann equation taken as benchmark demonstrate larger discrepancies with the experiments.

There are more involved macroscopic models based on the moments method, such as the R26 equations of Gu & Emerson (Reference Gu and Emerson2009), that exhibit higher accuracy in the transition regime. Nevertheless, modelling phenomena in rarefied gas dynamics with the R13 equations may provide perhaps the best compromise between a sufficiently complex mathematical description capable of capturing the most significant features not only in the bulk but also near the boundaries (Knudsen layer) and an analytically tractable set of equations. In any case, the exact solution presented here can be useful in the step of validation of numerical tools developed for the R13 equations simulating flow and heat transfer phenomena occurring either in the interior of or external to complicated geometries.

## Acknowledgements

We acknowledge the support of the Engineering and Physical Sciences Research Council (grant nos EP/N016602/1, EP/P020887/1, and EP/P031684/1) and the Leverhulme Trust (Research Project Grant). We also acknowledge and are thankful for informative communications from R. Bosworth and A. Ketsdever from the University of Colorado, Colorado Springs, regarding their experiments. We are grateful to Dr L. Gibelli at the University of Warwick for fruitful discussions.

## Appendix A. Multipole potentials

We define in this appendix some concepts from the theory of multipole potentials and briefly review some of its properties relevant to the work presented in the main body of this paper. For a more detailed discussion of the fundamentals of the method and its application to a variety of problems, the reader is referred to the textbooks by Leal (Reference Leal2007, chap. 8) – who names it as the method of superposition of vector harmonic functions – and Hess (Reference Hess2015, chap. 10). The material included in this review is taken from the latter.

The multipole potentials can be defined as tensorial solutions of the Laplace equation. There are two classes of multipole potentials, namely, descending and ascending potentials. The descending multipole potentials tend to zero when $r\rightarrow \infty$ and diverge when $r\rightarrow 0$ , where $r\equiv |\text{x}|$ and $r^{2}=x_{m}x_{m}$ . The descending multipole potentials are defined by

The first two descending potentials are given by

and

with the second potential being known as the dipole potential. The quadrupole and octupole potential tensors are

and

respectively. The rank-four multipole potential tensor is

Multipole potentials satisfy

This property is helpful in solving differential equations with elliptic operators. Multipole potential tensors are symmetric in any pair of indexes and, because they solve Laplace equation, also vanish after contracting any pair of indexes.

Finally, for problems in interior domains, the so-called ascending multipoles are needed. They arise because in (A 7), the factor $g^{\prime \prime }-2\ell r^{-1}g^{\prime }=r^{2\ell }(r^{-2\ell }g^{\prime })^{\prime }=0$ not only for $g=1$ – leading to the descending multipoles – but also for $g=r^{(2\ell +1)}$ . Therefore, Laplace equation is also solved by

which are known as ascending multipole potentials. They are zero at $r=0$ and, for $\ell >0$ , diverge when $r\rightarrow \infty$ .

## Appendix B. Solution of the ordinary differential equations

Consider the real-valued function
$\unicode[STIX]{x1D711}(x)$
. The ordinary differential equations in (3.20*a*
) and (3.20*c*
) can be represented in the generic form

where $n=0,1,2,\ldots$ and $\unicode[STIX]{x1D706}$ is a given parameter. The exact solution of this differential equation can be extracted from the handbook of solutions of ordinary differential equations by Zaitsev & Polyanin (Reference Zaitsev and Polyanin2002) (page 219). It can be written as

where
$J_{n+1/2}$
and
$Y_{n+1/2}$
are the Bessel functions of half-integer order of the first and second kind, respectively,
$i$
is the imaginary number, and
${\mathcal{A}}_{0}$
and
${\mathcal{B}}_{0}$
are arbitrary constants. With the relations (e.g. Arfken *et al.*, Reference Arfken, Weber and Harris2012, chap. 14; or Abramowitz & Stegun, Reference Abramowitz and Stegun1972, chap. 9)

and

where $I_{\unicode[STIX]{x1D708}}$ and $K_{\unicode[STIX]{x1D708}}$ are the modified Bessel functions of the first and second kind, respectively, and $\unicode[STIX]{x1D708}$ may be a complex number, and introducing the modified spherical Bessel functions (notice the different scaling factors in these definitions)

we can write solution (B 2) as

This result indicates that the substitution
$\unicode[STIX]{x1D711}(x)=x^{\unicode[STIX]{x1D700}}\tilde{\unicode[STIX]{x1D711}}(\unicode[STIX]{x1D706}x)$
in (B 1) with the choice
$\unicode[STIX]{x1D700}=n+1$
leads to the modified spherical Bessel differential equation for
$\tilde{\unicode[STIX]{x1D711}}$
, an expression that, unlike (B 1), is rather well-known. Because
$i_{n}(x)$
grows unbounded whereas
$k_{n}(x)$
tends to zero when
$x\rightarrow \infty$
, we set
${\mathcal{A}}_{0}=0$
in order to use (B 7) to represent the solutions of (3.20*a*
) and (3.20*c*
) in the main body of the document. Finally, we can write
$k_{n}$
in terms of elementary functions with the relations

which can be extended with the recurrence formula
$k_{n-1}(x)-k_{n+1}(x)=-(2n+1)k_{n}(x)/x$
(Arfken *et al.*
Reference Arfken, Weber and Harris2012). Using (B 9) and (B 10), we can obtain the expressions in (3.21*a*
) and (3.21*c*
), respectively.

## Appendix C. Expressions for trace-free symmetric tensors

For a rank-two tensor $\text{Q}$ , the associated symmetric trace-free tensor is given by

where $\mathbf{1}$ is the (rank-two) identity tensor and $\text{e}_{\unicode[STIX]{x1D6FC}}$ , with $\unicode[STIX]{x1D6FC}=$ 1, 2 and 3, denotes a set of orthonormal basis vectors. A repeated subscript in a term implies summation. Note that $\text{e}_{\unicode[STIX]{x1D6FC}}\boldsymbol{\cdot }\text{Q}\text{e}_{\unicode[STIX]{x1D6FC}}$ is just the transpose of $\text{Q}$ and $\text{Q}:\mathbf{1}$ its trace.

For a rank-three tensor $\text{A}$ , the corresponding symmetric tensor $\widetilde{\text{A}}$ is

and the corresponding symmetric trace-free tensor $\overline{\text{A}}$ is

We have written these expressions in vector notation starting with the expressions in Cartesian index notation for symmetric and trace-free tensors given in appendix A.2 of Struchtrup (Reference Struchtrup2005*b*
). For completeness, we show them here. The components of
$\widetilde{\text{A}}$
are

and the components of $\overline{\text{A}}$

## Appendix D. Thermophoretic force and drag from Grad’s 13-moment method and from other approaches

Following Young’s (Reference Young2011) work, the non-dimensional thermophoretic force and velocity drag resulting from Grad’s 13-moment method G13 are given by

and

respectively, where $Kn^{\prime }=\sqrt{\unicode[STIX]{x03C0}/2}\,Kn$ . For the results presented in this paper, the thermal creep, velocity slip and temperature jump coefficients take, respectively, the Maxwell–Smoluchowski values, namely, $K_{tc}~=~3/4$ , $C_{m}~=~1$ and $C_{e}~=~15/8$ (Young Reference Young2011). It should be said that the first closing parenthesis in the numerator of (D 1) and the factor $Kn^{\prime }$ after coefficient $C_{m}$ in the numerator of (D 2) as well as the second closing parenthesis in the denominators of both (D 1) and (D 2) are missing in Young’s (Reference Young2011) article (see his formulae (32a) and (32b)). When passing the limit $\unicode[STIX]{x1D6EC}\rightarrow \infty$ in (D 2), we recover the factor obtained by Lockerby & Collyer (Reference Lockerby and Collyer2016) in their formula (5.5), who have already noted the typographical errors in Young’s paper for $\unicode[STIX]{x1D6F9}$ .

The interpolation formula presented by Young (Reference Young2011) for the thermophoretic force is

with $K_{tc}~=~1.10$ , $C_{m}~=~1.13$ , $C_{e}~=~2.17$ and $C_{int}=0.5$ .

For completeness, we add the expressions used in this paper for the free-molecule regime ( $Kn^{\prime }\gg 1$ ). For the non-dimensional thermophoretic force, Waldmann (Reference Waldmann1959) obtained

whereas, for the drag caused by a free stream past a sphere, Epstein (Reference Epstein1924) obtained

These expressions are extracted from table 1 of Young’s article.