Transient behaviour of a rarefied gas around a sphere caused by impulsive rotation

Abstract The unsteady behaviour of a rarefied gas caused by a sudden change of the angular velocity of a sphere, placed in an otherwise quiescent gas, is investigated based on the linearized Bhatnagar–Gross–Krook model of the Boltzmann equation and the diffuse reflection boundary condition. The initial and boundary value problem is solved numerically by the method of characteristics, which is capable of tracking the discontinuity of the velocity distribution function moving in the phase space. The transient behaviour of the macroscopic quantities, such as the circumferential flow velocity and shear stress as well as the heat flow around the sphere, is clarified for a wide range of the Knudsen number. Furthermore, the long-time behaviour of the macroscopic quantities is elucidated, that is, they approach terminal values with a rate $t^{-3/2}$ for $t\gg 1$, with $t$ being a time variable. The analytical expression for the free molecular gas as well as for the slip flow is obtained. It is found that the circumferential heat flow reverses its direction as time proceeds when the Knudsen number is finite. More precisely, the direction is the same as that of the circumferential velocity of the sphere in the initial stage and opposite in the final stage, being reversed at some point of time depending on the distance from the sphere. This makes a clear contrast with the case of a free molecular gas, for which the heat flow is always in the direction of the sphere rotation in finite time and vanishes in the long-time limit.


Introduction
The transient behaviour of fluids caused by impulsive change of boundary data is of great interest in fluid mechanics. In rarefied gas dynamics, the most classical problem of this type is the linearized Rayleigh problem (e.g. Yang & Lees 1956;Gross & Jackson 1958;Cercignani & Sernagiotto 1964;Sone 1964), which is described as follows. Suppose that a rarefied gas is initially at equilibrium with an infinite plane at rest which is maintained at a uniform and constant temperature. At time t * = 0, the plane suddenly starts to move in its surface with a uniform velocity. We investigate the unsteady behaviour of the gas caused by the impulsive motion of the plane for t * > 0. The linearized Rayleigh problem describes the propagation of the transverse momentum into the gas and contains both the free-molecular-like and the continuum-flow-like behaviours in a single problem. These characters have been revisited in recent mathematical studies (Kuo 2011(Kuo , 2017. Later on, flows caused by a sudden change of boundary data have been investigated in various situations (e.g. Sone & Sugimoto 1990;Aoki et al. 1991;Doi 2016). However, most of the studies are limited to spatially one-dimensional problems. In the present paper, we extend the linearized Rayleigh problem to a spatially three-dimensional axisymmetric flow. More specifically, we consider a rigid sphere placed in a rarefied gas which is initially at thermal equilibrium with the sphere. Suppose that the sphere is suddenly set into a rotational motion around one of its diameters at t * = 0. We investigate the linear response of the gas to the impulsive rotation of the sphere for t * > 0.
Unlike the one-dimensional Rayleigh problem for a planer boundary, the sphere radius appears as a quantity having a physical dimension of length in addition to the molecular mean free path in the present problem. Thus, their quotient, known as the Knudsen number, plays an important role in characterizing unsteady response. Another important difference from the one-dimensional Rayleigh problem lies in the fact that the present problem admits a steady solution (Loyalka 1992;Andreev & Popov 2010;Taguchi, Saito & Takata 2019). Thus, the speed of approach to the final steady state is also of interest. A non-trivial response occurring in the heat flow will be discussed as well.
It should be mentioned that the instantaneous change in the boundary data introduces a ( jump) discontinuity in the velocity distribution function (VDF) at time t * = 0 on the sphere, which propagates in the phase space as time proceeds. Indeed, the abrupt change in the macroscopic quantities taking place in the initial stage is closely related to the discontinuity of the VDF. In addition, we have discontinuities originating from the fact that the sphere is a convex body, that is, due to the molecules making grazing collisions with the sphere. These effects are fully accounted for in the present study, devising a method of characteristics based on the integral form of our basic equation. To make the problem tractable, we employ the Bhatnagar-Gross-Krook (BGK) model (Bhatnagar, Gross & Krook 1954;Welander 1954) of the Boltzmann equation and the diffuse reflection boundary condition on the sphere.
The transient behaviour of a rarefied gas around a spherical body is important in many applications in microscience or nanoscience as well as in vacuum engineering. The present problem is expected to provide fundamental information on the unsteady behaviour of rarefied gas if combined with recent results focusing on oscillatory shear-driven flows (e.g. Park, Bahukudumbi & Beskok 2004;Hadjiconstantinou 2005;Sharipov & Kalempa 2008;Doi 2010;Yap & Sader 2012, 2016.
The rest of the paper is organized as follows. In § 2, the precise statement of the problem is given along with the formulation using the similarity solution. Then, we discuss the discontinuity of the VDF in § 3, followed by the introduction of a numerical method in § 4. We summarize analytical results for the free molecular gas as well as for the continuum flow in § 5. Section 6 presents numerical results, where we show the behaviours of the macroscopic quantities as well as that of the VDF. Section 7 is devoted to further discussion and § 8 concludes the paper.

Problem
Consider a sphere with radius L placed in a monatomic rarefied gas at rest with density ρ 0 , temperature T 0 , and pressure p 0 = ρ 0 RT 0 , where R is the specific gas constant (i.e. the Boltzmann constant divided by the mass of a molecule). There is no external force acting on both the gas and the sphere. At time t * = 0, the sphere starts to rotate impulsively around a diameter with constant angular velocity Ω * > 0. The situation is schematically described in figure 1. We investigate the unsteady behaviour of the gas for t * > 0 based on the BGK model of the Boltzmann equation and the diffuse reflection boundary condition, under the assumption that Ω * L is so small compared with the thermal speed (2RT 0 ) 1/2 that the equation and the boundary condition can be linearized about the initial equilibrium state at rest for t * < 0.
The linearized BGK equation and its boundary and initial conditions for the present problem are summarized as follows: where dζ = dζ 1 dζ 2 dζ 3 and the bracket notation in (2.1c) and in (2.4a-c), below, indicates The Ω and k in (2.1d) and (2.1a) are the dimensionless angular velocity and the scaled mean free path, respectively, and are defined by where 0 is the mean free path of the gas molecules in the equilibrium state at rest with density ρ 0 and temperature T 0 , and Kn = 0 /L is the Knudsen number. For the BGK model, 0 = (2/ √ π)(2RT 0 ) 1/2 /A c ρ 0 with A c being a constant such that A c ρ 0 is the collision frequency at the reference state. By our assumption, Ω is a small positive constant, i.e. Ω 1. We refer to § 1.11 of Sone (2007) for the linearized system of the BGK equation (or of the Boltzmann equation), where the formulation is given in a more general situation.
The (perturbed) pressure, the stress tensor and the heat-flow vector of the gas are defined as the moments of φ as follows: Similarity solution It is easy to verify that the following similarity solution is compatible with the present initial and boundary value problem: where ζ = (ζ 2 i ) 1/2 = |ζ | and is the angle of the molecular velocity measured from the radial direction.
The initial and boundary value problem for φ S is summarized as follows: Under the similarity solution (2.5), the macroscopic variables are expressed as It should be noted that the temperature of the gas is uniform and is equal to the sphere temperature, that is, τ = 0 everywhere. Finally, we introduce the torque (the moment of force) acting on the sphere. That is, if we denote by ( p 0 L 3 Ωh M , 0, 0) the moment of force around the origin, h M is expressed in terms ofP rϕ (r, t) as follows: (2.10) The h M in the steady state (t → ∞) was investigated in Loyalka (1992) and Taguchi et al. (2019). The time evolution of h M for various k is one of our interests in the present study.

Conservation equation
If we multiply (2.1a) (or (2.7a)) by ζ ϕ E (or ζ 2 ϕ sin θ E) and integrate the result with respect to the molecular velocity, we obtain In a steady state (∂/∂t = 0), we have r 3P rϕ = const. and thereforeP rϕ is inversely proportional to r 3 .
The backward molecular trajectory for given (r, θ ζ , ζ, t). Heret is the backward time and is related to a variable of integrationt in (3.3) ast = t −t. In panel (a), where θ ζ > Arcsin(r −1 ), the molecular trajectory can be traced back until the initial time without encountering the sphere. In panel (b), where θ ζ < Arcsin(r −1 ), the molecular trajectory can be traced back until the initial time only when the condition ζ t < η w is met; otherwise it hits the sphere surface at time The location of the discontinuity (Γ 1 and Γ 2 ) of VDF projected on the θ ζ -ζ t plane for various r. The Γ 1 and Γ 2 for the same r meet on a point on the curve ζ t = cot θ ζ .

Propagation of the discontinuity of VDF
Equation (2.1a) describes the variation of φ along the molecular trajectory (i.e. the characteristics of the equation). At a given point of the time t, let us consider a point x i in the gas region and trace back the molecular trajectory along the characteristic curve starting from x i in the backward direction i = −ζ i /ζ (and in the backward direction in time). Since the characteristic curve is a straight line for a given ζ i , it is clear that one can follow the backward trajectory until the initial time without encountering the sphere in the case where θ ζ of (2.6) satisfies the condition θ ζ > Arcsin(r −1 ) (see figure 2a). On the other hand, when θ ζ < Arcsin(r −1 ) as shown in figure 2(b), the trajectory can be traced back until the initial time without encountering the sphere only if the molecular speed ζ is smaller than a certain value; otherwise, the trajectory encounters the diffusely reflecting boundary (i.e. the sphere) before reaching the initial time.
To be more precise, in the case of θ ζ < Arcsin(r −1 ), let us consider a half-line drawn from the point x i in the direction i and let x wi be the nearest point (from x i ) that this half-line intersects with the sphere surface. We further denote the linear distance from x i to x wi by η w = η w (r, θ ζ ), which is given by η w (r, θ ζ ) = r cos θ ζ − 1 − r 2 sin 2 θ ζ , 0 θ ζ Arcsin(r −1 ). (3.1) With this definition, we conclude two cases: (i) when ζ t < η w , the trajectory is traced back to the initial time without intersecting the boundary; (ii) when ζ t > η w , it encounters the boundary before the initial time is reached. In the latter case, the time of intersection of the trajectory with the sphere surface is given by for given r, θ ζ , ζ and t satisfying θ ζ < Arcsin(r −1 ).
The behaviour of the macroscopic quantities as t → 0 + near the boundary deserves special attention. As shown above, φ S is discontinuous at ζ = η w (r, θ ζ )/t for a given (r, θ ζ , t) (see (3.5)). If we denote this value of ζ as ζ * , it is a function of (r, θ ζ , t), i.e. (3.7) With this notation, we write a part of the integral in (2.7c) as (3.8) Now consider a simultaneous limit t 0 and r 1 such that (r − 1)/t = const. Since θ ζ * π/2 as r 1, the third term vanishes, and we are left with the integral I = π/2 0 On the other hand, ζ * is expressed as near the boundary. Therefore, the upper or lower limit ζ * (= (r − 1)/t cos θ ζ ) appearing in the inner integrals depends on the value of (r − 1)/t. Since φ S has a jump discontinuity at ζ = ζ * , this implies that the whole integral, and thus the circumferential flow velocity u ϕ , takes different values depending on the ratio (r − 1)/t, namely, the speed of approach to the point (r, t) = (1, 0). Clearly, the same is true for other macroscopic variables P rϕ and Q ϕ . If we do not take the limit but consider a point close to (r, t) = (1, 0), the macroscopic variables are uniquely determined but undergo abrupt changes in the r-t plane due to the variation of ζ * . In this way, the abrupt change of the macroscopic quantities near the sphere shortly after the onset of the rotation is closely related to the discontinuity of VDF propagating in the phase space.

Numerical analysis
As we have seen in the previous section, the discontinuity of VDF moves in the phase space as time goes on. The situation is similar to those in moving boundary problems if we regard the four-dimensional surface Γ 1 ∪ Γ 2 as a movable boundary, changing its location in the phase space. A moving boundary problem for a rarefied gas was treated in Tsuji & Aoki (2013), where a numerical scheme based on the method of characteristics has been developed. Therefore, we adopt a similar strategy to solve the integral equation numerically.

Preliminary
For convenience, we rewrite (3.3) as follows. We multiply exp(t/k) to (3.3), and integrate it with respect to the time variable from t 0 to t. Then, we divide it by exp(t/k) and differentiate it with respect to t to obtain (4.1) where the integrand G(r, θ ζ , ζ, t,t) is given by withr being the one defined in (3.4a). The upper limit t end of the integral is given by that is, t end is t end = η w /ζ when the molecule hits the sphere, while t end = t if the molecule goes back to initial time (see figures 2a and 2b). The H(x) is the Heaviside step function andũ ϕ in (4.2) depends on φ S through (2.7c). The distancer(ζt, ·, ·) in (4.2) is computed from (3.4a) withη replaced by ζt. For ζ = 0, consider the limit ζ 0. The integral on the right-hand side of the integral equation (4.1) amounts to the integration ofũ ϕ , which depends on two variables (r, t). This reduces dramatically the difficulty of the numerical analysis.

Some remarks on the numerical method
In the numerical analysis, the range of r and that of ζ are restricted to finite intervals, that is, r ∈ [1, r max ] and ζ ∈ [0, ζ max ], where r max and ζ max are sufficiently large positive numbers. The appropriateness of the values r max and ζ max is judged from the numerical results. The range of the time variable t is also restricted to a finite range as t ∈ [0, t max ]. We introduce lattice points for (r, θ ζ , ζ, t) as follows: ζ , t (n) ) (see (3.6) and (3.7)), the functions in (4.5) are chosen in such a way that they satisfy It should be noted that VDF is discontinuous at θ ζ = θ (i) ζ * and at ζ = ζ (i,j,n) * . This is the reason why the whole intervals for θ ζ and ζ are divided at θ ζ = θ (i) ζ * and ζ = ζ (i,j,n) The lattice system {θ ). These dependencies are not shown explicitly in (4.5c) and (4.5d) to shorten the notation.
We introduce discretized variables for φ S andũ ϕ as follows: For a given n( 1), let us suppose that the circumferential flow velocityũ (i,n ) ϕ is known for all n = 0, . . . , n − 1 and i = 0, . . . , N r . Then, we determine φ (i,j,m,n) S at time t = t (n) for all i, j and m by the following scheme: and η w , G and t end are given by (3.1), (4.2) and (4.3), respectively. Note that the two cases 0 θ ζ θ ζ * and θ ζ * < θ ζ π in (4.1) are unified in the above formula under the convention η (i,j) w = ∞ (and ζ * = ∞) for the latter case. For the integration with respect tõ t in (4.8), we further introduce discrete points fort as follows: where gt(x) is a monotonically increasing function that defines lattice points fort. We then used the four-point Gauss-Legendre quadrature to evaluate the integral in (4.8), applying a suitable interpolation forũ ϕ (r, t −t) in the integrand (essentially the second-order Lagrange interpolation, applied first to the r variable and then to the t variable). In the case using the data at a few preceding time steps by extrapolation (basically using the second-order Lagrange polynomials), and then apply the same interpolation in the r-t plane.
In the above procedure, we obtain φ . The last step is to carry out the numerical integration with respect to ζ and θ ζ in (2.7c) to obtainũ ϕ . Again, the four-point Gauss-Legendre quadrature is applied for both ζ and θ ζ variables. Other macroscopic quantitiesP rϕ andQ ϕ in (2.9) are also updated in this step using the same quadrature. Then, we proceed to the next time step.
We have chosen our lattice system carefully, inspecting the behaviour of the integrand. We give further information on the numerical analysis in appendix A.

Analytical results for k → ∞ and for k 1
Prior to presenting the numerical result, we show some analytical results obtained for the cases of k → ∞ and for k 1.
5.1. Free molecular gas First, we consider the case of the free molecular (or collisionless) gas, i.e. k → ∞. Letting k → ∞ in (4.1), we have where ζ * and θ ζ * are given by (3.6) and (3.7). Substituting (5.1) into (2.7c) and (2.9), the macroscopic quantities are obtained as follows: where erfc(x) is the complementary error function defined by On the sphere, the macroscopic variables take the following values: which are time independent. Thus, the torque on the sphere is given by In other words, the torque on the sphere is constant in time in the free molecular gas.
For a large t such that r/t 1, the following expression can be derived: from which we conclude that, as t → ∞, corresponding to the steady solution (Taguchi et al. 2019). Note that for large r the steady flow velocity is further simplified to Asymptotic behaviour for k 1 and t = O(k −1 ) Now we turn our attention to the case where k 1 and t = O(k −1 ). The discussion in this section is based on Sone (2007) and Takata & Hattori (2012), and is not restricted to the BGK model.
According to Sone (2007) and Takata & Hattori (2012), the macroscopic quantities of interest h (h = u ϕ , P rϕ , Q ϕ ) can be sought in the form where h H is the overall solution (i.e. the Hilbert solution) subject to the conditions ∂h H /∂ x i = O(h H ) and ∂h H /∂t = O(kh H ) on the spatial and temporal variations, and h K is the so-called Knudsen-layer correction which is appreciable only in a thin layer (the Knudsen layer) adjacent to the boundary, whose thickness is of the order of the mean free path. The Knudsen-layer correction is required to satisfy the condition In other words, we seek the solution whose temporal variation is slow (i.e. the diffusion time scale), disregarding the abrupt temporal variation which takes place in the early stage of the evolution (the initial and subsequent acoustic layers).

If we further expand h H and h K
the circumferential flow velocity u ϕHm is seen to be described by the following partial differential equation (the Stokes equation):  Sone (2002Sone ( , 2007. and γ 1 is a constant (transport coefficient) whose numerical value depends on the particular molecular model under consideration, see table 1. Physically, γ 1 is the dimensionless viscosity; the viscosity μ 0 at the reference equilibrium state is given by Note also that variable s is used to describe the slow temporal variation and u Hm is regarded as u Hm = u Hm (x i , s). In the derivation of (5.12), we have used the fact that the pressure is uniform (P = 0).
In this study, we obtain the asymptotic expression of the flow velocity up to the order k. Then the appropriate boundary conditions for the flow velocity are where b (1) 1 is the shear-slip coefficient, whose numerical value (under the diffuse reflection condition) is listed in table 1.
The (5.12) for m = 0 and 1 with the above boundary conditions should be solved under the following natural initial condition: The solution to (5.12) with m = 0 or 1 subject to the initial and boundary conditions, (5.15) and (5.14), is obtained by introducing a vector potential (Landau & Lifshitz 1987) and applying the Laplace transform. We thus obtain u ϕH0 and u ϕH1 and, consequently, u ϕ (= u ϕH0 + k(u ϕH1 + u ϕK1 )) for k 1 is written as where s has been changed back to kt and the function Y (1) 1 (x) is explained later in this section. The solution is expected to describe the slow evolution of the gas for k 1 over the time t O(1/k), which comes after more abrupt changes with time scales t ∼ O(k) and t ∼ O(1). The lower-order formula or the formula without the Knudsen-layer correction can be derived from (5.16) by discarding some of the terms as listed in table 2.
The corresponding formulae for the shear stress and the circumferential component of the heat-flow vector in the gas, as well as that for the torque on the sphere, are obtained as Leading order TABLE 2. The correspondence of the asymptotic formulae (5.16)-(5.18) for k 1 and the expansions of u ϕ , P rϕ and Q ϕ in k. Note that P rϕH0 , P rϕK1 , P rϕK2 , Q ϕH0 and Q ϕH1 are identically zero and they do not appear in the second column.
where γ 3 is a constant (transport coefficient, see table 1), H (1) is defined by  (2007) . On the other hand, the numerical values of H (1) i (x) (i = 4, 5, 6) can be found in Takata & Hattori (2015). The Y (1) 1 , H (1) 1 , etc. are universal functions in the sense that they are not problem dependent but are determined once the set of molecular model and kinetic boundary condition is specified. Therefore, we can exploit the numerical data already known in the literature for the present study. We show the Knudsen-layer functions in figure 3 for the convenience of readers.
The leading-order formula for u ϕ (=u ϕH0 ) for k 1, which corresponds to the solution of the Stokes equation with a no-slip boundary condition, is given by (5.16) with b (1) 1 = Y (1) 1 = 0 (see table 2). The corresponding leading-order formulae for P rϕ (= kP rϕH1 ) and Q ϕ (= kQ ϕK1 ) are derived from (5.17) 6 (x) and H (1) (x) under the diffuse reflection condition: (a) a hard-sphere gas; (b) BGK model. Data reconstructed from Takata & Hattori (2015). the BGK model), and hence, the leading-order heat flow is negative as a whole. In this way, the asymptotic theory predicts the occurrence of a heat flow, which is in the opposite direction to the overall flow velocity. We will give a further discussion on the behaviour of the heat flow for small k later in § 7, where we compare the formula (5.18) with the numerical solutions.
For a sufficiently large t such that (γ 1 kt) 1/2 /(r − 1) 1, equations (5.16)-(5.19) with b (1) 1 = Y (1) 1 = γ 3 = H (1) = 0 (the leading terms) can be expanded to give (5.21d) Thus, the approach to the steady state is proportional to t −3/2 in the continuum flow both in the flow velocity and the shear stress as well as in the heat-flow vector.
We now present the numerical results for the macroscopic quantities. Figures 5-7 show the behaviour of the macroscopic quantities for three different values of k, i.e. k = 10, 1 and 0.1; figure 5 shows the circumferential flow velocity, figure 6 the shear (tangential) stress, and figure 7 the circumferential component of the heat-flow vector. In these figures, panels (a,b) show the results for k = 10, panels (c,d) those for k = 1 and panels (e, f ) those for k = 0.1. Panels (a,c,e) show the spatial profiles at several values of t, while panels (b,d, f ) the temporal profiles at several values of r. Note that u ϕ /Ω sin θ(=ũ ϕ ), P rϕ /Ω sin θ(=P rϕ ) and Q ϕ /Ω sin θ(=Q ϕ ) depend on r and t when k is specified.
The disturbance caused by the sphere rotation propagates into the gas as time goes on. The flow velocity on the sphere (r = 1) is exactly 0.5 for the free molecular gas, 1 for the continuum flow (the Stokes equation with no-slip boundary condition) and is between 0.5 and 1 for 0 < k < ∞. The flow velocity gradually increases with time in the whole gas region. The disturbance propagates into the gas more slowly when k is smaller (see panels b,d, f ). However, the influence of the sphere rotation penetrates deeper into the gas when k is smaller and the gas attains a larger rotational speed overall for smaller k.
The value of the shear stressP rϕ is larger for large k. It increases abruptly after the onset of the self-rotation near the sphere and then decreases gradually with time. The disturbance propagates into the gas more slowly when k is smaller, as in the case ofũ ϕ (see panels b,d, f ).
The heat flux is localized in the vicinity of the sphere immediately after the onset of the self-rotation; it flows in the same direction as the sphere rotation and forms a peak-like profile near the boundary. The region of positive heat flow gradually propagates into the gas and spreads out. At the same time,Q ϕ decreases near the boundary and changes its sign to negative values, implying that the heat flow changes its flow direction in the opposite sense as the sphere rotation. The negative heat flow is more prominent for k = 1 than for k = 10 and gradually develops in the whole region as time goes on. It should be noted that no negative heat flow occurs for the free molecular gas, as can be verified from (5.2c). When k = 0.1, the negative heat flow is even greater, and is more confined near the boundary. It may be noted that the heat flux is negative in the whole gas region in the steady state for 0 < k < ∞ (Taguchi et al. 2019).
These features, particularly those in the initial stage near the boundary, are summarized in figure 8, where the isolines of (a)ũ ϕ , (b)P rϕ and (c)Q ϕ are plotted in the r-t plane for k = 0.1. We clearly see that the contour lines accumulate near the origin, indicating that the steep spatial and temporal variations take place in the vicinity of (r, t) = (1, 0). The reason of this abrupt change was discussed in the last paragraph of § 3.

Long-time behaviour
In this section, we discuss the long-time behaviour of the macroscopic quantities. In figure 9 the behaviours of the macroscopic quantities are shown until relatively large time (t ≈ 2500) for k = 1. More specifically, figure 9(a-c) show the circumferential component of the flow velocityũ ϕ , the shear stressP rϕ , and the circumferential component of the heat-flow vectorQ ϕ as functions of r for various times. The steady solution obtained by a hybrid method combining the finite-difference method and the method of characteristics (Taguchi et al. 2019) is also included in each panel. We can observe that the profile of each macroscopic variable approaches the corresponding steady solution as time is increased.
In figure 9(c), where the semilogarithmic plot of |Q ϕ | versus r is shown, the positive and negative parts ofQ ϕ are represented by the solid and broken curves, respectively. Note that, in the steady state,Q ϕ is everywhere negative. At t = 2500, the curve almost overlaps with the steady solution in the part r 20. The positive part ofQ ϕ keeps travelling in the r direction as time goes on, while decreasing its height and increasing the width (see also figure 7a,c,e).
To obtain further information on the long-time behaviour, we show in figure 10 the time derivative of each macroscopic variable as a function of t for various r: (a) ∂ũ ϕ /∂t; (b) ∂P rϕ /∂t; and (c) ∂Q ϕ /∂t. From panels (a,b), it is clearly seen that ∂ũ ϕ /∂t and ∂P rϕ /∂t tend to decay in proportion to t −5/2 for large t. The time derivative of the heat flow ∂Q ϕ /∂t is also likely to decay in proportion to t −5/2 for t 1, although it is more difficult to see the tendency. Thus, the rate of approach of the macroscopic quantities to the steady profiles is t −3/2 . Figure 11 contains further information on the approach of the macroscopic quantities to the steady state for different k. Here, we show the slope of the double-logarithmic plot of ∂h/∂t versus t at r = 1 for various k, where h =ũ ϕ orQ ϕ . We computed the slope α by applying a standard finite difference. If we do so, however, we observed a fluctuation developed after some moment, which is attributed to rounding errors. The fluctuated raw data are shown for k = 0.1 in each panel by grey lines. Note that the slope forQ ϕ is more fluctuated than that ofũ ϕ because the magnitude of ∂Q ϕ /∂t is much smaller. To smooth out the data, we took a moving average over 100 dimensionless times. The data thus obtained are shown by the symbols in the figures. The figure for h =P rϕ , which is similar to that of u ϕ , is omitted here. It is seen from the figure that the slope tends to approach α = −5/2 as the time proceeds, irrespective of k. This confirms the exponent of t −3/2 in the long-time behaviour.
6.3. Torque acting on the sphere In this section, we discuss the transient behaviour of the torque acting on the rotating sphere. Figure 12 α + (5/2) FIGURE 11. The temporal variation of the slope α of ln |∂h/∂t| ∼ αt + const., where h =ũ ϕ or Q ϕ at r = 1 for various k: (a) h =ũ ϕ ; (b) h =Q ϕ . The moving averages over 100 dimensionless times are shown by the symbols, which are connected by the solid lines. The fluctuated raw data are shown for k = 0.1 by grey lines. (1) 1 = 0, which corresponds to the continuum flow (i.e. the Stokes equation with no-slip boundary condition), is shown by the dashed curve. The asymptotic formula for the steady solution ((4.4b) in Taguchi et al. (2019)), which includes second-order slip corrections, is shown by the horizontal long-dashed line. Note that the asymptotic formula (5.19) is valid only for t k −1 . and is monotonically decreasing in t for k < ∞. It attains its supremum at t = 0 + , which coincides with the value for the free molecular gas. In case of k < ∞, after the onset of the rotation, −h M undergoes a rapid decrease in proportion to t followed by a subsequent slower decay in proportion to t −3/2 approaching the final steady state.
In figure 12(b), we compare the numerical solution and the asymptotic formula (5.19) for small k for a wider range of t (k = 0.1 and 0.01). Here, the solid curves represent our numerical results, while the dash-dotted curves the asymptotic formula (5.19). A favourable agreement is observed for k = 0.01 except for small values of t( 1). In the case of k = 0.1, the discrepancy is noticeable but the shape of the curve is well represented by the asymptotic expression. For comparison, we also show the result based on (5.19) with b (1) 1 = 0 by the dashed curve, which corresponds to the continuum flow (i.e. the Stokes equation with no-slip boundary condition). The deviation is more prominent, particularly for k = 0.1. Note that the asymptotic formula is valid for t k −1 . Nevertheless, a reasonable agreement between the numerical and asymptotic results is observed even for smaller values of t.
The approach to the steady torque is algebraic and is proportional to t −3/2 for t 1. Thus, the tangential stress on the sphereP rϕ for t 1 can be approximated by a functioñ where P 0 and P 1 are constants which depend on k. The P 0 and P 1 can be obtained by fitting the first two terms of the above expression with the numerical data ofP rϕ (1, t) for large t. It turned out that adding a third term, P 2 t −5/2 , gives a better fitting result. Then, the coefficients (P 0 , P 1 , P 2 ) were obtained by the least squares method, where the numerical data for t 10 2 were used for k 0.1 and those for t 8 × 10 2 for k < 0.1. Note that the presence of the t −5/2 -term can be verified in the case of the continuum flow from (5.17). The results thus obtained are summarized in figure 13 and in table 3 for P 0 and P 1 . It should be noted that P 0 corresponds to lim t→∞Prϕ | r=1 which can also be obtained from the steady solution (Taguchi et al. 2019). The latter, obtained by a numerical method different from the present approach, is also included in table 3 for comparison. As seen from the table, P 0 obtained by the fitting agrees well with the stationary result. The solid curve in figure 13(a) is the asymptotic formula for k 1 for the steady solution ((4.4b) in Taguchi et al. (2019)). The P 1 is decreasing in k and, thus, decays faster for larger k. The coefficient P 1 is likely to behave as P 1 ∝ k −5/2 for k 1 (see figure 13b). On the other hand, P 1 for small k is approximated by (see (5.21d)), where γ 1 = 1 for the BGK model. Equation (6.2) is shown in figure 13(b) by the dashed line.

Discussion
We have seen that the circumferential heat flow is induced by the rotation of the sphere. Let us first discuss its initial peak-like behaviour localized near the boundary. We know that the heat flow is given by (5.2c) in the case of a free molecular gas. If we consider the limit t → 0 such that (r − 1)/t = const., this reduces tõ (7.1) Note that the above expression describes the heat flow induced by an impulsive onset of a uniform tangential motion of a flat plate in a collisionless gas (i.e. one-dimensional 10 -2 10 -1 10 0 10 1 10 2 10 -3 10 -2 10 -1 10 0 10 1 10 2 Asymptotic formula for the steady solution (Taguchi et al. 2019) FIGURE 13. The coefficients P 0 and P 1 in (6.1): (a) P 0 versus k; (b) P 1 versus k. The symbol • shows the numerical result. In panel (a), the dashed line indicates the value at k → ∞ (the free molecular flow), while the solid curve the asymptotic formula for k 1 for the steady solution ((4.4b) in Taguchi et al. (2019)). In panel (b), the solid line indicates (6.2), while the broken line the slope corresponding to P 1 ∝ k −5/2 .   linearized Rayleigh problem for a free molecular gas), if r − 1 is interpreted as the distance from the plate. On the other hand, in the case of a finite k < ∞, it can be shown that the heat flow is also described by (7.1) as the same limit t → 0 + is approached. To see this, let us consider the following change of variables: where ε is a small positive number that has been introduced to make y andτ of the order of unity. The φ 1 satisfies the equation where ζ = (ζ 2 x + ζ 2 y + ζ 2 z ) 1/2 . Applying the same coordinate transformation to the boundary condition at r = 1 as well as to the initial condition, we see that the problem is reduced to a one-dimensional Rayleigh problem for a collisionless gas so far as the leading-order terms in ε are concerned. In this way, the initial behaviour of the heat flow in the vicinity of the sphere is reduced to (7.1), irrespective of the values of k < ∞.
We note thatQ ϕ of (7.1) is self-similar in (r − 1)/t (or y/τ ), non-negative and approaches zero as r → 1 or r → ∞ when t > 0 is fixed. Consequently, the maximum value does not change in time, which is 1/2 √ 2eπ, while the profile becomes wider as t is increased. In this manner, the initially localized peak near the sphere spreads into the gas. As t (and r) is increased, ε becomes no longer small, and the curvature effect as well as molecular collision begins to play a role. As a result,Q ϕ starts to behave differently from that of the one-dimensional free molecular solution. This behaviour near the origin (r, t) = (1 + , 0 + ) is also consistent with the discussion given at the end of § 3.
Next, we discuss the behaviour of the heat flow Q ϕ at relatively large t by restricting ourselves to the case of small k. Figure 14(a) shows a comparison of the numerical result with the asymptotic formula (5.18) for k = 0.05 and 0.01 in the region close to r = 1 at t = 100. The agreement between the numerical solution and the asymptotic solution is good, motivating us to discuss the behaviour of Q ϕ for small k based on the asymptotic theory. In the present problem, the most significant term in the asymptotic expansion Q ϕ = Q ϕH0 + k(Q ϕH1 + Q ϕK1 ) + k 2 (Q ϕH2 + Q ϕK2 ) + · · · is kQ ϕK1 , since Q ϕH0 and Q ϕH1 are identically zero (see the caption of table 2). The heat conduction due to Fourier's law is irrelevant since the temperature is uniform. The contribution of kQ ϕK1 , describing the Knudsen layer adjacent to the boundary, is shown by the dashed curves in the figure. This heat flow in the Knudsen layer does not vanish as t → ∞ and contributes to the negative heat flows in the long-time limit, as mentioned in § 5.2. Note also that the magnitude of the heat flow is proportional to k. Therefore, it is of the same order as the shear stress, although there is an essential difference in that the heat flow is confined in the Knudsen layer while the shear stress extends over the gas region.
The second-order Hilbert part k 2 Q ϕH2 in the heat flow merits some attention. The origin of this term is the Laplacian of the flow velocity (see (16b) in Takata & Hattori (2012)), which reduces to Q ϕH2 = 1 2 (γ 3 /γ 1 )∂u ϕH0 /∂s, s = kt, by the use of (5.12) in the present case. Hence, this second-order gas-rarefaction effect is related to the transient behaviour (and therefore vanishes in the steady state, i.e. t → ∞). The term is included in the asymptotic solution shown in figure 14(a) (dash-dotted curve). However, its modulus is much smaller than the typical magnitude of kQ ϕK1 at t = 100, and it cannot be Profiles ofQ ϕ (=Q ϕ /Ω sin θ) for k = 0.01 at various t ranging from t = 100 to 1000. Note the difference between the scales of panels (a,b) in the ordinate. In panels (a,b), the solid curve represents the numerical solution of the BGK equations, the dash-dotted curve the asymptotic formula (5.18) and the dashed curve the leading-order term kQ ϕK1 in the asymptotic solution (see table 2). The inset in panel (a) is a close-up of the case k = 0.01 near r = 1. In panel (b) the profiles of kQ ϕK1 + k 2 Q ϕH2 (with no second-order Knudsen-layer correction included) are also shown by the (green) long-dashed curves, which overlap with the dash-dotted curves representing (5.18). and decreasing tendency of the circumferential heat flow outside the Knudsen layer is well described by the time derivative of the overall circumferential flow velocity u ϕH0 . We close this section by presenting the leading-order heat flow kQ ϕK1 , (5.18) with γ 3 = H (1) = b (1) 1 = 0 (see table 2), in a dimensional form for a steady state (t → ∞). Let us denote the (dimensional) heat-flow vector by q. Noting the relation between k and the viscosity μ 0 (see the sentence after (5.12)), q in the steady state is expressed as where Ω * and x * are the angular velocity of the sphere and the position vector, respectively. From this formula, we clearly see that the magnitude of the (steady) heat flow is proportional to the viscosity.

Concluding remarks
In this paper, we investigated the transient behaviour of a rarefied gas caused by an impulsive onset of the rotating motion of a sphere based on the linearized BGK equation and the diffuse reflection boundary condition. This problem can be viewed as an extension of the one-dimensional linearized Rayleigh problem to a three-dimensional axisymmetric flow. Special attention was paid to the propagation of the discontinuity of VDF in the phase space. For this purpose, we carried out a numerical analysis based on the integral form of the BGK equation. Our findings are summarized as follows.
(i) The behaviours of the macroscopic quantities have been clarified for a wide range of the Knudsen number including the free molecular flow and the continuum flow. (ii) The heat flow, which is initially localized near the boundary and flows in the same direction as the sphere rotation, spreads into the gas as time goes on. As the final steady state is approached, the region with negative heat flow, flowing in the opposite direction as the sphere rotation, is established. (iii) (The long-time behaviour.) The rate of approach of the macroscopic quantities to the steady state is in proportion to t −3/2 for t 1. In other words, the gas-rarefaction effect does not change the decaying rate and the degree of gas rarefaction (the Knudsen number) affects merely the proportionality constant. The free molecular flow is an exception. (iv) The transient behaviour of the torque on the sphere is clarified. For large t 1, the magnitude of the term proportional to t −3/2 becomes small in proportion to k −5/2 for k 1, indicating the faster approach to the steady torque for larger k.
The grid points for θ ζ , ζ andt are regenerated at every time step in an adaptive manner. We omit details on the grids for these variables.
We have checked the accuracy of our numerical computations in various ways, for instance, by changing the number of grid points, grid parameters, etc. Here, we provide some numerical data to assess the present computations based on the conservation law. As shown in (2.11), the quantity on the left-hand side should be identically zero, while it is not in an actual numerical computation due to numerical errors. Thus, the deviation from zero is a measure of computational accuracy. Let us introduce