text:
# Abstract

## Abstract

We propose a systematic expansion method which is applied to freely evolving granular fluids contained in sufficiently small systems. Restricting ourselves to small systems, we show that there exists a small parameter which characterizes a typical size of density and temperature inhomogeneities. A solution of the hydrodynamic equations for a fluid of inelastic hard spheres is expanded in the small parameter. It is shown that our method quantitatively describes the asymptotic state of the system, such as flow profiles, density and temperature inhomogeneities, and the decay law of the global temperature and the energy per particle.

## PACS numbers

45.70.Mg, 45.70.Qj, 47.85.Dh, 02.30.Mv.

## Key words

Granular fluid; instabilities; pattern formation; hydrodynamic equations; asymptotic state; perturbation expansion.

## 1 Introduction

Granular matter shows fluid like behavior under certain conditions , such as Couette flow, convection flow under vibration, and flow down an inclined chute. A freely cooling rapid granular flow is an idealized limiting case of granular fluids in the absence of gravity and without any energy input. By analyzing their molecular dynamics (MD) simulations of a system of smooth inelastic hard discs, Goldhirsch and Zanetti  discovered that the freely cooling granular fluids exhibit interesting instabilities. When the system is prepared in a spatially homogeneous state, it slowly develops patterns in the flow field (vortices) and in the density field (clusters).

A linear stability analysis  of hydrodynamic equations for inelastic hard spheres (IHS) has revealed that the initial spatially homogeneous cooling state is unstable to the formation of vortices and clusters, and provided a good description of the initial growth of unstable hydrodynamic modes. At large times, however, these unstable hydrodynamic modes grow so large that nonlinear hydrodynamic effects, such as convection and viscous heating, may dominate the evolution of the system. Nonlinear effects in freely cooling granular fluids have been discussed by several authors . In Refs. , the importance of viscous heating effects has been pointed out and the threshold of the clustering instability has been given in terms of the coefficient of restitution and a typical length scale of a fluctuation. Nonlinear coupling between hydrodynamic modes through viscous heating has been studied  by means of numerical analysis of a hydrodynamic model and direct Monte Carlo simulation (DMCS) of the Boltzmann equation. Recently, strong numerical evidence has been shown  that a one-dimensional freely cooling granular gas asymptotically behaves as a sticky gas, and that the inviscid Burgers equation is an appropriate continuum theory. These three works are concerned with unstable growth regime in large systems, which is the physically most interesting problem of hydrodynamics for the freely cooling granular fluids. This regime is characterized by intrinsic patterns in the flow field and in the density field where the typical correlation length of these patterns are small compared to the system size. Unfortunately in this regime only analytical large time results have been obtained from a mode coupling theory .

Recently some progress in the analysis of nonlinear hydrodynamic equations has been made for freely cooling granular fluids contained in systems that are sufficiently small, such that the clustering instability is suppressed . A linear stability analysis  predicts that if the smallest wave number of a system, $`k_0=2\pi /L`$, is smaller than a critical value $`k_{}^{}`$, the homogeneous cooling state becomes unstable and the flow field develops a shear flow pattern. Several groups have observed by MD simulations  and by DMCS  of IHS contained in two-dimensional systems with periodic boundaries that there is a regime of $`k_0`$ for a given coefficient of restitution where the system develops shearing patterns, but no strong density inhomogeneities. In this regime the flow field eventually reaches a stationary shearing pattern that has the following properties: (i) a stationary shear flow profile with a period equal to the system size in the direction perpendicular to the direction of the flow , (ii) exponential decay of the energy per particle as a function of $`\tau `$, the average number of collisions suffered per particle in time $`t`$ , (iii) $`t^2`$-decay of the energy per particle , (iv) inhomogeneities in density  and temperature  with a period that is half the period of the shear flow profile; the density concentrates into regions of the peaks of the flow profile, and (v) saturation of the density profile . The properties (i)-(iii) have also been observed by MD simulations in five and six dimensions . In Ref. , the evolution of shear flow patterns in systems with $`k_0k_{}^{}`$ ($`k_0<k_{}^{}`$) has been studied. In such systems, unstable shear modes grow so slowly that the remaining hydrodynamic modes are enslaved by the shear modes, and the amplitudes of the unstable shear modes remain small. On the basis of this consideration, they have derived amplitude equations for the unstable shear modes and explained the above properties in a quantitative manner. For similar small systems, in Ref.  a Landau-Ginzburg-type equation for a nonconserved order parameter, the shear rate tensor, has been derived under an approximation in which inhomogeneities in density and temperature are neglected. It has been shown that the Landau-Ginzburg-type equation describes saturation of unstable shear modes due to viscous heating effects, and describes the properties of shear flow profile (i)-(iii) at large times in a quantitative manner in spite of drastic simplification achieved by the approximation.

This paper is also concerned with the small system regime; we propose a systematic expansion method to obtain solutions which correspond to the asymptotic state of sufficiently small systems. The method is applied to the stationary shear flow profile with the properties (i)-(v), which appears in sufficiently small systems with periodic boundaries. Although the asymptotic state of the small systems with periodic boudaries is not a physical state, we concentrate in this paper only on this case, which has been studied well so far by the computer simulations  and by the theories . The method is based on an expansion in the Prandtl number $`P_r`$ of a solution of the hydrodynamic equations for IHS. At large times in small systems, viscous heating, which is proportional to the kinematic viscosity $`\nu `$, due to growing flow field patterns becomes dominating effects that induce inhomogeneities in temperature . In sufficiently small systems, the inhomogeneities in temperature is largely reduced by thermal conduction, which is proportional to the thermal conductivity $`\kappa `$. In these small systems, a typical size of temperature inhomogeneities can be determined by these two competing effects; hence it is proportional to the Prandtl number $`P_r\nu /\kappa `$. We will show that by expanding the solutions in the Prandtl number, the inhomogeneities in density and temperature field in the small system regime can be described; corrections to shear flow profile, and the energy and temperature decay laws caused by these density and temperature inhomogeneities are obtained in a systematic manner. This method has the approximation used in Ref.  that density and temperature inhomogeneities are neglected as its zeroth approximation. Our method can describe systems with finite amplitudes of unstable shear modes, to which the method in Ref.  can not be applied. The idea of expanding solutions of hydrodynamic equations in a small Prandtl number was originally developed for classical fluids by Busse  in order to describe an oscillatory instability of convection rolls that appear in a fluid layer heated from below.

This paper is organized as follows. In Sec.2, we introduce the hydrodynamic equations for IHS and show that we can eliminate the global temperature, defined as the spatial average of the local temperature, from the hydrodynamic equations by introducing rescaled hydrodynamic variables. Then, it is shown that in the small system regime some of the properties (i)-(v) can be associated with the existence of a stationary solution of the hydrodynamic equations for the rescaled variables. In Sec.3, the method of systematic expansion in the Prandtl number is presented. We first discuss the basic idea of the expansion scheme, then show how the expansion scheme works up to the second approximation. A quantitative description of the properties (i)-(v) is given in the first approximation. It presents qualitatively new predictions of the behavior of $`d`$-dimensional IHS fluids. The range of validity of our method is discussed in the end of Sec.3. We end with some conclusions in Sec.4.

## 2 Hydrodynamics of IHS

### 2.1 Hydrodynamic equations for IHS

We consider a system of $`N`$ inelastic hard spheres, contained in a $`d`$-dimensional cubic system with sides $`L`$ and volume $`V=L^d`$. The boundary conditions are periodic boundary conditions. At an initial time, $`t=0`$, the system is prepared in a homogeneous state with temperature $`T_0`$ and vanishing flow field. We assume that weekly inelastic hard sphere systems can be described by the hydrodynamic equations supplemented by a term $`\mathrm{\Gamma }`$ that expresses the rate of collisional energy loss ,

$`_tn+𝐮n`$ $`=`$ $`n𝐮,`$ (1)
$`_t𝐮+𝐮𝐮`$ $`=`$ $`{\displaystyle \frac{1}{mn}}p+{\displaystyle \frac{2}{mn}}\left(\eta \text{D}\right)+{\displaystyle \frac{1}{mn}}\left(\zeta 𝐮\right),`$ (2)
$`_tT+𝐮T`$ $`=`$ $`{\displaystyle \frac{2p}{dn}}𝐮+{\displaystyle \frac{2}{dn}}(\kappa T)`$ (3)
$`+`$ $`{\displaystyle \frac{4\eta }{dn}}(\text{D}:\text{D})+{\displaystyle \frac{2\zeta }{dn}}(𝐮)^2\mathrm{\Gamma }.`$
In this paper the inelasticity is always assumed to be small. The shear rate D is the symmetrized dyadic, $`\text{D}_{\alpha \beta }=(_\alpha u_\beta +_\beta u_\alpha \frac{2}{d}\delta _{\alpha \beta }𝐮)/2`$, and $`A:B=_{\alpha =1}^d_{\beta =1}^dA_{\alpha \beta }B_{\beta \alpha }`$. The transport coefficients, the shear viscosity $`\eta `$, the bulk viscosity $`\zeta `$, and the heat conductivity $`\kappa `$, depend on the local density and temperature.

We observe that the total energy per particle consists of a convective and internal energy

$`E=E_{conv}+E_{int}={\displaystyle \frac{1}{N}}{\displaystyle 𝑑𝐫\left[\frac{1}{2}mnu^2+\frac{d}{2}nT\right]},`$ (4)
with rate of change, derived from Eqs. (1)-(3)

$`{\displaystyle \frac{dE}{dt}}`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle 𝑑𝐫\frac{d}{2}n\mathrm{\Gamma }},`$ (5)
$`{\displaystyle \frac{dE_{int}}{dt}}`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle }d𝐫[p𝐮+2\eta \text{D}:\text{D}+\zeta (𝐮)^2{\displaystyle \frac{d}{2}}n\mathrm{\Gamma }].`$ (6)
The terms on the right hand side of the last equation represent the work done by the flow, the gain of $`E_{int}`$ due to viscous heating and the loss of $`E_{int}`$ due to collisional dissipation.

On the basis of kinetic theory one can derive that the rate of collisional energy loss, $`\mathrm{\Gamma }=2\gamma _0\omega T`$, is proportional to the collision frequency $`\omega (n,T)`$ multiplied by the fraction of energy $`ϵT`$ lost per collision , where $`\gamma _0=ϵ/2d=(1\alpha ^2)/2d`$ is the inelasticity parameter and $`\alpha `$ is the coefficient of restitution.

### 2.2 Elimination of the global temperature

In this section, we will introduce a rescaling of the hydrodynamic variables by means of the global granular temperature defined as

$`\overline{T}(t)={\displaystyle \frac{2}{d}}E_{int}(t)={\displaystyle \frac{1}{N}}{\displaystyle 𝑑𝐫n(𝐫,t)T(𝐫,t)},`$ (7)
or equivalently the mean random (thermal) velocity $`\overline{v}_0=\sqrt{2\overline{T}/m}`$. Similarly we will use the global density $`\overline{n}=N/V`$.

The freely cooling hard sphere fluid has no intrinsic energy scale: there is no energy scale in parameters which characterize interactions between particles and interactions between a particle and a boundary. Hence, if we rescale time by $`tat`$ at fixed parameters that characterize interactions between particles and interactions between a particle and a boundary, then velocities of all particles, energy, and temperature are rescaled as $`va^1v`$, $`Ea^2E`$, and $`Ta^2T`$, respectively, but the trajectories of particles are not affected by this rescaling. Therefore, we can always keep one of energy scales in the system constant by rescaling time in an appropriate way. This fact has been utilized by Soto et al.  to enhance accuracy of their MD simulations by keeping the total energy constant. In this article, keeping the global temperature constant is suitable for our purpose because the collision frequency and the transport coefficients in the hydrodynamic equations for IHS are given as functions of temperature and density. The rescaling of the hydrodynamic variables given below has been used in a linear stability analysis .

We first introduce the dimensionless time and the scaled hydrodynamic variables

$`d\tau `$ $`=`$ $`\overline{\omega }dt,\stackrel{~}{n}={\displaystyle \frac{n}{\overline{n}}},\stackrel{~}{𝐮}={\displaystyle \frac{𝐮}{\overline{v}_0}},\stackrel{~}{T}={\displaystyle \frac{T}{\overline{T}}},`$ (8)
where $`\overline{\omega }=\omega (\overline{n},\overline{T})`$ is the collision frequency taken at the global temperature $`\overline{T}(t)`$ and the global density $`\overline{n}`$, and $`\tau `$ is the average number of collisions, suffered by a particle in time $`t`$. As a result of the rescaling, the scaled global temperature $`(1/N)𝑑𝐫\stackrel{~}{n}\stackrel{~}{T}`$ always remains constant. The collision frequency $`\omega `$ and the transport coefficients $`\{\eta ,\zeta ,\kappa \}`$ in the hydrodynamic equations will be calculated from the Enskog theory for hard spheres, in principle for inelastic ones . Moreover, $`\overline{\eta }=\eta (\overline{n},\overline{T})`$, $`\overline{\zeta }=\zeta (\overline{n},\overline{T})`$, $`\overline{\kappa }=\kappa (\overline{n},\overline{T})`$ denote the corresponding Enskog transport coefficients taken at $`\overline{T}(t)`$ and $`\overline{n}`$. Note that in the Enskog theory for hard spheres, $`\overline{\omega },\overline{\eta },\overline{\zeta },\overline{\kappa }\sqrt{\overline{T}}`$. Then, we get a closed set of equations for these scaled variables,

$`_\tau \stackrel{~}{n}+\stackrel{~}{𝐮}l_0\stackrel{~}{n}`$ $`=`$ $`\stackrel{~}{n}l_0\stackrel{~}{𝐮},`$ (9)
$`_\tau \stackrel{~}{𝐮}+\stackrel{~}{𝐮}l_0\stackrel{~}{𝐮}`$ $`=`$ $`{\displaystyle \frac{1}{2\stackrel{~}{n}}}l_0\stackrel{~}{p}+D_{}{\displaystyle \frac{2}{\stackrel{~}{n}}}\left(\eta ^{}(\stackrel{~}{n})\sqrt{\stackrel{~}{T}}\stackrel{~}{\text{D}}\right)+{\displaystyle \frac{\overline{\zeta }}{\overline{\rho }\overline{\omega }}}{\displaystyle \frac{1}{\stackrel{~}{n}}}\left(\zeta ^{}(\stackrel{~}{n})\sqrt{\stackrel{~}{T}}\stackrel{~}{𝐮}\right)`$ (10)
$``$ $`{\displaystyle \frac{1}{2}}\left(_\tau \mathrm{ln}\overline{T}\right)\stackrel{~}{𝐮},`$
$`_\tau \stackrel{~}{T}+\stackrel{~}{𝐮}l_0\stackrel{~}{T}`$ $`=`$ $`{\displaystyle \frac{2}{d}}{\displaystyle \frac{\stackrel{~}{p}}{\stackrel{~}{n}}}l_0\stackrel{~}{𝐮}+{\displaystyle \frac{\overline{c}_p}{c_v}}D_T{\displaystyle \frac{1}{\stackrel{~}{n}}}\left(\kappa ^{}(\stackrel{~}{n})\sqrt{\stackrel{~}{T}}\stackrel{~}{T}\right)`$ (11)
$`+`$ $`{\displaystyle \frac{8}{d}}D_{}{\displaystyle \frac{\eta ^{}(\stackrel{~}{n})}{\stackrel{~}{n}}}\sqrt{\stackrel{~}{T}}(\stackrel{~}{\text{D}}:\stackrel{~}{\text{D}})+{\displaystyle \frac{4}{d}}{\displaystyle \frac{\overline{\zeta }}{\overline{\rho }\overline{\omega }}}{\displaystyle \frac{\zeta ^{}(\stackrel{~}{n})}{\stackrel{~}{n}}}\sqrt{\stackrel{~}{T}}(\stackrel{~}{𝐮})^2`$
$``$ $`2\gamma _0\omega ^{}(\stackrel{~}{n})\stackrel{~}{T}^{3/2}\left(_\tau \mathrm{ln}\overline{T}\right)\stackrel{~}{T},`$
where we have introduced the following notation:

$`{\displaystyle \frac{\eta }{\overline{\eta }}}`$ $`=`$ $`\eta ^{}(\stackrel{~}{n})\sqrt{\stackrel{~}{T}},{\displaystyle \frac{\zeta }{\overline{\zeta }}}=\zeta ^{}(\stackrel{~}{n})\sqrt{\stackrel{~}{T}},{\displaystyle \frac{\kappa }{\overline{\kappa }}}=\kappa ^{}(\stackrel{~}{n})\sqrt{\stackrel{~}{T}},`$ (12)
$`{\displaystyle \frac{\omega }{\overline{\omega }}}`$ $`=`$ $`\omega ^{}(\stackrel{~}{n})\sqrt{\stackrel{~}{T}},{\displaystyle \frac{p}{\overline{n}\overline{T}}}=\stackrel{~}{n}\stackrel{~}{T}p^{}(\stackrel{~}{n})=\stackrel{~}{p}(\stackrel{~}{n},\stackrel{~}{T}).`$ (13)
The functions $`f^{}(\stackrel{~}{n})`$ depend on $`\stackrel{~}{n}`$ and the coefficient of restitution $`\alpha `$, and $`\eta ^{}(1)=\zeta ^{}(1)=\kappa ^{}(1)=\omega ^{}(1)=1`$. The constants $`D_{}=\overline{\nu }/\overline{\omega }`$ and $`D_T=\overline{\kappa }/\overline{c}_p\overline{\rho }\overline{\omega }`$ are the rescaled shear diffusivity and the heat diffusivity, respectively, and $`c_v=d/2m`$ and $`\overline{c}_p=c_p(\overline{n})`$ are the specific heat per unit mass at constant volume and pressure. The mean free path $`l_0`$ is defined by $`l_0=\overline{v}_0/\overline{\omega }_0`$. It is important to note that the constants $`D_{},D_T,\overline{\zeta }/\overline{\rho }\overline{\omega }`$ and $`l_0`$ are independent of $`\overline{T}`$. The macroscopic equation for the global temperature $`\overline{T}`$ can be obtained from $`dE_{int}/dt`$ in Eq. (6) and subsequent rescaling, and yields,

$`_\tau \mathrm{ln}\overline{T}`$ $`=`$ $`{\displaystyle \frac{1}{V}}{\displaystyle }d𝐫[{\displaystyle \frac{2}{d}}\stackrel{~}{p}l_0\stackrel{~}{𝐮}+{\displaystyle \frac{8}{d}}D_{}\eta ^{}(\stackrel{~}{n})\sqrt{\stackrel{~}{T}}(\stackrel{~}{\text{D}}:\stackrel{~}{\text{D}})`$ (14)
$`+{\displaystyle \frac{4}{d}}{\displaystyle \frac{\overline{\zeta }}{\overline{\rho }\overline{\omega }}}\zeta ^{}(\stackrel{~}{n})\sqrt{\stackrel{~}{T}}(\stackrel{~}{𝐮})^22\gamma _0\omega ^{}(\stackrel{~}{n})\stackrel{~}{n}\stackrel{~}{T}^{3/2}].`$
Hence, the time derivative of $`\mathrm{ln}\overline{T}`$ is a functional of the scaled hydrodynamic variables $`\stackrel{~}{n},\stackrel{~}{𝐮}`$ and $`\stackrel{~}{T}`$, and the global temperature can be completely eliminated from the equations for these scaled hydrodynamic variables.

### 2.3 Asymptotic state

Equations (9)-(11) for the freely evolving IHS-fluid are supposedly describing the growing patterns in the density and flow fields starting from a spatially homogeneous initial state, and the basic objective is to describe their asymptotic evolution for small inelasticity, where the time scale of collisional cooling is well separated for the free time between collisions. Comparison of scaled field Eqs. (9)-(11) with Eq. (14) for the global temperature suggests that there might be two totally different scenarios leading to an asymptotic state. One small system scenario, where the fields equations describe a fast process with a typical scale $`\tau _{}L^2/D_{}`$. For time $`\tau `$ large compared to $`\tau _{}`$ the scaled hydrodynamic variables reach a stationary solution and the rate of change of the global temperature, i.e. the right hand side of Eq. (14), approaches a nonvanishing constant. In the large system scenario, the field equations contain arbitrary slow hydrodynamic modes. The system, in which the global temperature initially decays like $`\mathrm{exp}(2\gamma _0\tau )`$ according to Haff’s law , gradually reaches an asymptotic state with inhomogeneities in the flow field, where the viscous heating almost entirely compensates the collisional cooling, and where the right hand side of Eq. (14) becomes very small. Some consequences of the second scenario have been analyzed in the mode coupling calculation of Ref. .

Here we concentrate on a sufficiently large time scale, $`\tau \tau _{}`$, in the small system scenario. There are some important observations that can be deduced from the assumption of the existence of a stable stationary solution of the scaled hydrodynamic equations. The results under the approximation in which inhomogeneities in density and temperature are neglected have already been reported in Ref. . The assumption of the existence of a stable stationary solution implies that the right hand side of Eq. (14) is constant, $`(_\tau \mathrm{ln}\overline{T})2\gamma _a<0`$. This leads the following important consequences: First, for large $`\tau `$ the global temperature $`\overline{T}`$ in Eq. (14) decays exponentially in terms of $`\tau `$, i.e.

$`\overline{T}e^{2\gamma _a\tau }.`$ (15)
Second, the exponential decay (15) of the global temperature implies its $`t^2`$-decay in the real time $`t`$. According to the definition of $`\tau `$ given in Eq. (8), the evolution of the global temperature $`\overline{T}`$ in terms of $`t`$ becomes,

$`_t\overline{T}`$ $`=`$ $`{\displaystyle \frac{d\tau }{dt}}_\tau \overline{T}=\overline{\omega }\overline{T}(_\tau \mathrm{ln}\overline{T})={\displaystyle \frac{\omega _0}{\sqrt{T_0}}}(_\tau \mathrm{ln}\overline{T})\overline{T}^{3/2},`$ (16)
where $`\omega _0(\overline{n},T_0)`$ is the collision frequency at the initial time and $`T_0`$ is the initial temperature. Here we have used $`\overline{\omega }\sqrt{\overline{T}}`$. At sufficiently large times $`(_\tau \mathrm{ln}\overline{T})`$ can be replaced by its stationary value $`2\gamma _a`$, i.e.

$`_t\overline{T}`$ $`=`$ $`{\displaystyle \frac{\omega _0}{\sqrt{T_0}}}\mathrm{\hspace{0.17em}2}\gamma _a\overline{T}^{3/2}.`$ (17)
This equation implies that $`\overline{T}`$ decays as $`t^2`$ at large times, which has been indeed observed in simulations. More precisely

$`\overline{T}(t)={\displaystyle \frac{T_0}{\gamma _a^2\omega _0^2(tt_e)^2}}{\displaystyle \frac{T_0}{\gamma _a^2\omega _0^2}}t^2,`$ (18)
where $`t_e`$ is an integration constant. It should be noted that the prefactor of $`t^2`$ is a constant which is independent of transient behavior of the system and of the initial temperature. We compare this decay with the initial decay in the homogeneous cooling state (Haff’s law)

$`\overline{T}(t)=T_0e^{2\gamma _0\tau }={\displaystyle \frac{T_0}{(1+\gamma _0\omega _0t)^2}}{\displaystyle \frac{T_0}{\gamma _0^2\omega _0^2}}t^2,`$ (19)
where the decay rate in $`\tau `$ is larger because $`\gamma _0>\gamma _a`$ . Third, because the collision frequency $`\overline{\omega }`$ is proportional to the square root of the global temperature $`\omega \sqrt{\overline{T}}`$, the relation (8) between $`\tau `$ and $`t`$ in the asymptotic state is given by

$`dt`$ $``$ $`{\displaystyle \frac{1}{\sqrt{\overline{T}}}}d\tau e^{\gamma _a\tau }d\tau .`$ (20)
Finally, it can be shown that a stationary solution of the scaled hydrodynamic variables implies that the energy per particle is proportional to the global temperature:

$`E`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle 𝑑𝐫\left[\frac{m}{2}n𝐮^2+\frac{d}{2}nT\right]}=\overline{T}\left[{\displaystyle \frac{1}{V}}{\displaystyle 𝑑𝐫\stackrel{~}{n}\stackrel{~}{𝐮}^2}+{\displaystyle \frac{d}{2}}\right].`$ (21)
These results seem to be consistent with the simulation results (i)-(v) in Sec.1. Hence, we shall assume that for sufficiently small systems Eqs. (9)-(11) have a stable stationary solution, and concentrate on how to obtain the stationary solution and the exponent $`\gamma _a`$ by means of an systematic approximation method.

It is important to note that the above argument is valid also for systems surrounded by inelastic walls as well as for systems with periodic boundaries; both systems have no intrinsic energy scale.

## 3 Systematic expansion method

### 3.1 Basic idea

In this section we present the basic idea of a systematic expansion method to solve Eqs. (9)-(11) for the stationary state, i.e. with the vanishing time derivatives. In this argument we shall assume that viscous heating plays an essential role in inducing inhomogeneities in the density and temperature fields on the basis of the previous studies  which have shown importance of viscous heating effects in the formation of inhomogeneities. It will be shown later that the result obtained by the systematic expansion method indeed supports this assumption.

The heat locally generated by viscous heating effects is redistributed by convection and thermal conduction, and partially compensated by collisional dissipation. A systematic expansion method to be presented in this paper works when thermal conduction plays a dominant role in the reduction of temperature inhomogeneities through the redistribution and compensation of heat. The conditions necessary for this situation to be realized will be examined below.

We shall first discuss a typical size of temperature and density inhomogeneities when effects of convection and collisional dissipation are small compared with thermal conduction effects, so that they may be neglected. Denoting a typical size and length scale of the scaled flow field are given by $`\stackrel{~}{U}`$ and $`l_u`$, respectively, the viscous heating term in Eq.(11) is given, in order of magnitude, by

$`{\displaystyle \frac{8}{d}}D_{}{\displaystyle \frac{\eta ^{}(\stackrel{~}{n})}{\stackrel{~}{n}}}\sqrt{\stackrel{~}{T}}(\stackrel{~}{\text{D}}:\stackrel{~}{\text{D}})`$ $``$ $`{\displaystyle \frac{4}{d}}D_{}L^2\stackrel{~}{U}^2\left(1+O(\mathrm{\Delta }\stackrel{~}{T},\mathrm{\Delta }\stackrel{~}{n})\right),`$ (22)
where $`\mathrm{\Delta }\stackrel{~}{T}`$ and $`\mathrm{\Delta }\stackrel{~}{n}`$ are a typical size of $`\stackrel{~}{T}1`$ and a typical size of $`\stackrel{~}{n}1`$, respectively, and assumed to be small. We shall show later that this assumption is valid for $`\stackrel{~}{U}^2\stackrel{<}{}1`$. The thermal conduction term in Eq.(11) is similarly given by

$`{\displaystyle \frac{\overline{c}_p}{c_v}}D_T{\displaystyle \frac{1}{\stackrel{~}{n}}}\left(\kappa ^{}(\stackrel{~}{n})\sqrt{\stackrel{~}{T}}\stackrel{~}{T}\right){\displaystyle \frac{\overline{c}_p}{c_v}}D_Tl_T^2\mathrm{\Delta }\stackrel{~}{T}\left(1+O(\mathrm{\Delta }\stackrel{~}{T},\mathrm{\Delta }\stackrel{~}{n})\right),`$ (23)
where $`l_T`$ is a typical length scale of temperature inhomogeneities. Because the temperature inhomogeneities are induced by viscous heating, which is proportional to the square of $`\stackrel{~}{𝐮}`$, $`l_T`$ can be approximately related to $`l_u`$ by $`l_Tl_u/2`$. The typical size of the inhomogeneous part of the scaled temperature $`\mathrm{\Delta }\stackrel{~}{T}`$ is then obtained from the balance between these two terms (22) and (23), i.e.,

$`\mathrm{\Delta }\stackrel{~}{T}{\displaystyle \frac{c_v}{\overline{c}_p}}{\displaystyle \frac{P_r}{d}}\stackrel{~}{U}^2\left(1+O(\mathrm{\Delta }\stackrel{~}{T},\mathrm{\Delta }\stackrel{~}{n})\right),`$ (24)
where $`P_r=D_{}/D_T`$ is the Prandtl number. A typical size of the inhomogeneous part of the scaled density can be estimated as follows: In the asymptotic state we may assume that the pressure is homogeneous as a consequence of the mechanical balance. Hence the typical size of the inhomogeneous part of the scaled density $`\mathrm{\Delta }\stackrel{~}{n}`$ is related to $`\mathrm{\Delta }\stackrel{~}{T}`$ by

$`\mathrm{\Delta }\stackrel{~}{n}`$ $`=`$ $`\overline{T}\overline{\beta }_B\mathrm{\Delta }\stackrel{~}{T}=\overline{\beta }_0\mathrm{\Delta }\stackrel{~}{T},`$ (25)
where the constant $`\overline{\beta }_B`$ is the bulk expansion coefficient; its definition is given in Appendix A. It is shown in Appendix A that for inelastic hard sphere fluids $`\overline{\beta }_0=\overline{T}\overline{\beta }_B`$ is a constant independent of $`\overline{T}`$, and in two and three dimensions $`\overline{\beta }_0`$ is less than $`1`$, i.e., $`\mathrm{\Delta }\stackrel{~}{n}<\mathrm{\Delta }\stackrel{~}{T}`$.

Values of $`(c_v/\overline{c}_p)(P_r/d)`$ for elastic hard spheres of various densities are listed in Table 1. As shown in Table 1, the factor $`(c_v/\overline{c}_p)(P_r/d)`$ is a small quantity and hence $`\mathrm{\Delta }\stackrel{~}{T}1`$ if $`\stackrel{~}{U}^2\stackrel{<}{}1`$. This result justifies the approximation used in Ref , in which density and temperature inhomogeneities are neglected, for a finite size of shear flow profile with $`\stackrel{~}{U}^2\stackrel{<}{}1`$ in the asymptotic state.

The fact that the inhomogeneous part of the scaled temperature and density is small, and is proportional to the Prandtl number motivates to expand the stationary solutions in the Prandtl number, expecting that it will eventually lead an expansion in the small parameter $`(c_v/\overline{c}_p)(P_r/d)`$, which characterizes the small magnitude of density and temperature inhomogeneities when $`\stackrel{~}{U}^2\stackrel{<}{}1`$. It will be shown that up to a part of the second approximation studied in this paper, this is indeed the case.

In Sec.4, we will discuss the range of validity of this method and give a system size above which our analysis breaks down. It is noteworthy that temperature profile produced by the combination of viscous hating effects and thermal conduction effects appears also in flows of classical fluids such as the compressible Couette flow for a fluid confined between isothermal or adiabatic walls .

Finally, let us consider the conditions under which thermal conduction may be regarded as dominant in comparison with convection and collisional dissipation. The condition that the convection term in Eq. (11), $`l_0\stackrel{~}{𝐮}\stackrel{~}{T}l_0\stackrel{~}{U}l_T^1\mathrm{\Delta }\stackrel{~}{T}`$, is negligible compared with the thermal conduction term yields,

$`{\displaystyle \frac{l_0\stackrel{~}{U}l_T^1\mathrm{\Delta }\stackrel{~}{T}}{\frac{\overline{c}_p}{c_v}D_Tl_T^2\mathrm{\Delta }\stackrel{~}{T}}}=\left({\displaystyle \frac{c_v}{\overline{c}_p}}{\displaystyle \frac{P_r}{d}}\right){\displaystyle \frac{d}{2}}R_e1.`$ (26)
Here $`R_e`$ is the Reynolds number defined by $`R_eUl_u/\nu `$, where $`U=\overline{v}_0\stackrel{~}{U}`$ is a typical size of the flow field (not scaled). Smallness of the parameter $`(c_v/\overline{c}_p)(P_r/d)`$ suggests that this condition is fulfilled if

$`{\displaystyle \frac{d}{2}}R_e\stackrel{<}{}1.`$ (27)
The inhomogeneous part of the collisional dissipation term, which contributes to reduction of inhomogeneities in the temperature field, is in order of magnitude given by $`\gamma _0\mathrm{\Delta }\stackrel{~}{T}`$. The condition that this term is small compared with the thermal conduction term yields,

$`{\displaystyle \frac{\gamma _0\mathrm{\Delta }\stackrel{~}{T}}{\frac{\overline{c}_p}{c_v}D_Tl_T^2\mathrm{\Delta }\stackrel{~}{T}}}=\left({\displaystyle \frac{c_v}{\overline{c}_p}}{\displaystyle \frac{P_r}{d}}\right){\displaystyle \frac{d}{4}}{\displaystyle \frac{\gamma _0}{D_{}l_u^2}}1.`$ (28)
Since $`(c_v/\overline{c}_p)(P_r/d)1`$, this condition is fulfilled if

$`{\displaystyle \frac{d}{4}}{\displaystyle \frac{\gamma _0}{D_{}l_u^2}}\stackrel{<}{}1.`$ (29)
For the stationary solutions with the properties (i)-(v), which are to be studied in this paper, however, the property (iv) and the assumption of homogeneous pressure suggest that temperature are inhomogeneous only in directions perpendicular to the direction of flow; hence the convection term identically vanishes. Therefore the condition (27) is not required for this type of solutions. Besides, it will be shown in Sec.3.4 that for solutions with the properties (i)-(v), the condition (29) is fulfilled if $`\stackrel{~}{U}^2\stackrel{<}{}1`$.

On the basis of the consideration denoted above, we expand the stationary solution in the Prandtl number as follows:

$`\stackrel{~}{n}`$ $`=`$ $`1+ϵ\stackrel{~}{n}^{(1)}+ϵ^2\stackrel{~}{n}^{(2)}+\mathrm{},`$
$`\stackrel{~}{𝐮}`$ $`=`$ $`\stackrel{~}{𝐮}^{(0)}+ϵ\stackrel{~}{𝐮}^{(1)}+ϵ^2\stackrel{~}{𝐮}^{(2)}+\mathrm{},`$
$`\stackrel{~}{T}`$ $`=`$ $`1+ϵ\stackrel{~}{T}^{(1)}+ϵ^2\stackrel{~}{T}^{(2)}+\mathrm{}.`$ (30)
Here we have introduced a parameter $`ϵ`$ that is of $`O(P_r)`$.

We require that the conservation of the total number and the total momentum, $`N=𝑑𝐫n`$ and $`\mathrm{𝟎}=𝑑𝐫n𝐮`$, is satisfied in each order in $`ϵ`$. This gives rise to the following conditions:

$`0`$ $`=`$ $`{\displaystyle 𝑑𝐫\stackrel{~}{n}^{(l)}}\text{for }l=1,2,\mathrm{},`$ (31)
$`\mathrm{𝟎}`$ $`=`$ $`{\displaystyle 𝑑𝐫\stackrel{~}{𝐮}^{(0)}},\mathrm{𝟎}={\displaystyle 𝑑𝐫\left(\stackrel{~}{𝐮}^{(1)}+\stackrel{~}{n}^{(1)}\stackrel{~}{𝐮}^{(0)}\right)},\mathrm{}.`$ (32)
In addition, rewriting the definition of the global temperature as $`1=(1/V)𝑑𝐫\stackrel{~}{n}\stackrel{~}{T}`$ and requiring that this relation is satisfied in each order in $`ϵ`$, we get the following conditions for the scaled temperature:

$`0`$ $`=`$ $`{\displaystyle 𝑑𝐫\stackrel{~}{T}^{(1)}},0={\displaystyle 𝑑𝐫\left(\stackrel{~}{n}^{(1)}\stackrel{~}{T}^{(1)}+\stackrel{~}{T}^{(2)}\right)},\mathrm{}.`$ (33)
We substitute (30) into the hydrodynamic equations (9)-(11) with the vanishing time derivatives

$`\stackrel{~}{𝐮}\stackrel{~}{n}`$ $`=`$ $`\stackrel{~}{n}\stackrel{~}{𝐮},`$ (34)
$`\stackrel{~}{n}\stackrel{~}{𝐮}l_0\stackrel{~}{𝐮}`$ $`=`$ $`{\displaystyle \frac{l_0}{2}}\stackrel{~}{p}+2D_{}\left(\eta ^{}(\stackrel{~}{n})\sqrt{\stackrel{~}{T}}\stackrel{~}{\text{D}}\right)+{\displaystyle \frac{\overline{\zeta }}{\overline{\rho }\overline{\omega }}}\left(\zeta ^{}(\stackrel{~}{n})\sqrt{\stackrel{~}{T}}\stackrel{~}{𝐮}\right)`$ (35)
$``$ $`{\displaystyle \frac{1}{2}}\left(_\tau \mathrm{ln}\overline{T}\right)\stackrel{~}{n}\stackrel{~}{𝐮},`$
$`\stackrel{~}{n}\stackrel{~}{𝐮}l_0\stackrel{~}{T}`$ $`=`$ $`{\displaystyle \frac{2}{d}}\stackrel{~}{p}l_0\stackrel{~}{𝐮}+ϵ^1{\displaystyle \frac{\overline{c}_p}{c_v}}D_T\left(\kappa ^{}(\stackrel{~}{n})\sqrt{\stackrel{~}{T}}\stackrel{~}{T}\right)`$ (36)
$`+`$ $`{\displaystyle \frac{8}{d}}D_{}\eta ^{}(\stackrel{~}{n})\sqrt{\stackrel{~}{T}}(\stackrel{~}{\text{D}}:\stackrel{~}{\text{D}})+{\displaystyle \frac{4}{d}}{\displaystyle \frac{\overline{\zeta }}{\overline{\rho }\overline{\omega }}}\zeta ^{}(\stackrel{~}{n})\sqrt{\stackrel{~}{T}}(\stackrel{~}{𝐮})^2`$
$``$ $`2\gamma _0\omega ^{}(\stackrel{~}{n})\stackrel{~}{n}\stackrel{~}{T}^{3/2}\left(_\tau \mathrm{ln}\overline{T}\right)\stackrel{~}{n}\stackrel{~}{T}.`$
Here both sides of Eqs. (10) and (11) are multiplied by $`\stackrel{~}{n}`$ for later convenience. We put $`ϵ^1`$ in front of the thermal conduction term, because $`P_r^1`$ appears there by rescaling collision time $`\tau `$ by $`L^2/D_{}`$ and by rescaling space $`r`$ by $`L`$. We omit this procedure since this rescaling of variables by constants does not change the final result.

Collecting terms of the same order, we obtain a set of equations that determine a stationary solution of the scaled hydrodynamic variables in each order in $`ϵ`$. The smallness parameter $`ϵ`$ is set equal to $`1`$ in the end of calculation.

### 3.2 Zeroth approximation

The terms of $`O(ϵ^1)`$ gives an identity $`0=0`$. The zeroth approximation, $`\stackrel{~}{𝐮}^{(0)}`$, is determined by the terms of $`O(ϵ^0)`$, i.e.

$`0`$ $`=`$ $`\stackrel{~}{𝐮}^{(0)},`$ (37)
$`\stackrel{~}{𝐮}^{(0)}l_0\stackrel{~}{𝐮}^{(0)}`$ $`=`$ $`D_{}^2\stackrel{~}{𝐮}^{(0)}{\displaystyle \frac{1}{2}}\left(_\tau \mathrm{ln}\overline{T}\right)^{(0)}\stackrel{~}{𝐮}^{(0)},`$ (38)
where

$`\left(_\tau \mathrm{ln}\overline{T}\right)^{(0)}`$ $`=`$ $`{\displaystyle \frac{1}{V}}{\displaystyle }d𝐫[{\displaystyle \frac{4}{d}}D_{}(\stackrel{~}{𝐮}^{(0)})^{}:\stackrel{~}{𝐮}^{(0)}2\gamma _0],`$ (39)
where $`(\stackrel{~}{𝐮}^{(0)})_{\alpha \beta }^{}=_\beta \stackrel{~}{u}_\alpha ^{(0)}`$. We show in Appendix B that if $`k_0>k_{}^{}`$, the nonlinear equation (38) has only the trivial solution $`\stackrel{~}{𝐮}^{(0)}=\mathrm{𝟎}`$. Here $`k_0=2\pi /L`$ is the smallest wave number of the system and $`k_{}^{}=\sqrt{\gamma _0/D_{}}`$ is the critical wave number of the instability in shear modes, predicted in a linear stability analysis . In this case all higher order corrections vanish, suggesting that the system is at large times in a homogeneous cooling state that follows Haff’s law , $`\left(_\tau \mathrm{ln}\overline{T}\right)=\left(_\tau \mathrm{ln}\overline{T}\right)^{(0)}=2\gamma _0`$, independently of the initial state of the system.

Hereafter, we shall consider only the non-trivial case $`k_0<k_{}^{}`$. If $`k_0<k_{}^{}`$, there is a large set of solutions  that satisfy Eqs. (37)-(39). In this paper we do not pursue to find all possible solutions that satisfy Eqs. (37)-(39). We show that there is indeed a solution that corresponds to the flow profile observed in simulations  with a period of the system size $`L`$. An intuitive explanation of the flow profile with the period $`L`$ is that this is the mode which shows the slowest hydrodynamic decay . Assuming the flow profile on the basis of the observations, we concentrate on a quantitative description of the amplitude of the flow profile, which eventually leads us to a theoretical prediction of the amplitude of inhomogeneous temperature and density profile.

One can show that the solutions of Eqs. (37)-(39) with a period of the system size $`L`$ have to be a flow that is inhomogeneous in directions which are perpendicular to the direction of the flow \[see Appendix C\], i.e.

$`\stackrel{~}{𝐮}^{(0)}(𝐫)={\displaystyle \underset{\alpha =1}{\overset{d^{}}{}}}{\displaystyle \underset{\beta =d^{}+1}{\overset{d}{}}}A_{\alpha \beta }\widehat{𝐞}_\beta \mathrm{cos}(k_0r_\alpha +\theta _{\alpha \beta }),`$ (40)
where we have chosen $`\alpha =1,2,\mathrm{},d^{}`$ ($`d^{}1`$) as to be directions in which $`\stackrel{~}{𝐮}^{(0)}`$ is inhomogeneous, and $`\beta =d^{}+1,\mathrm{},d`$ as the remaining directions. $`\widehat{𝐞}_\beta `$ is a unit vector in the direction $`\beta `$, and $`\theta _{\alpha \beta }`$ is an arbitrary phase factor. The amplitudes $`A_{\alpha \beta }`$ are real numbers that satisfy the following conditions. First, from Eq. (37), we get $`A_{\alpha \alpha }=0`$. Second, if $`k_0<k_{}^{}`$, Eq. (38) give rise to a relation :

$`A_0^2{\displaystyle \underset{\alpha =1}{\overset{d^{}}{}}}{\displaystyle \underset{\beta =d^{}+1}{\overset{d}{}}}A_{\alpha \beta }^2={\displaystyle \frac{d(\gamma _0D_{}k_0^2)}{D_{}k_0^2}}.`$ (41)
The proof of the results (40) and (41) are given in Appendix C. From Eqs. (40) and (41) we obtain a typical size of the scaled flow field $`\stackrel{~}{U}=[(1/V)𝑑𝐫|\stackrel{~}{𝐮}^{(0)}(𝐫)|^2]^{1/2}=A_0/\sqrt{2}`$.

The achievement of one of the stationary solutions can be viewed as a process of spontaneous symmetry breaking . Good quantitative agreement about the amplitude $`A_0`$ between the theoretical prediction (41) and the result of MD simulation has been shown in Ref. .

Substituting the solution (40) into Eq. (39), we get

$`\left(_\tau \mathrm{ln}\overline{T}\right)^{(0)}`$ $`=`$ $`2D_{}k_0^2.`$ (42)
### 3.3 First approximation

First order correction $`\stackrel{~}{T}^{(1)}`$ can be obtained from the terms of $`O(ϵ^0)`$,

$`^2\stackrel{~}{T}^{(1)}`$ $`=`$ $`{\displaystyle \frac{c_v}{\overline{c}_p}}{\displaystyle \frac{4}{d}}P_r\{(\stackrel{~}{𝐮}^{(0)})^{}:\stackrel{~}{𝐮}^{(0)}{\displaystyle \frac{1}{V}}{\displaystyle 𝑑𝐫(\stackrel{~}{𝐮}^{(0)})^{}}:\stackrel{~}{𝐮}^{(0)}\},`$ (43)
where we have used a relation $`\stackrel{~}{𝐮}^{(0)}:\stackrel{~}{𝐮}^{(0)}=0`$, which is obtained from Eqs. (37) and (38). The relation (43) manifests that the temperature inhomogeneities are excited by the inhomogeneous part of the viscous heating term. Substituting the expression (40) and solving the Laplace equation under the condition (33), we get

$`\stackrel{~}{T}^{(1)}(𝐫)`$ $`=`$ $`{\displaystyle \frac{c_v}{\overline{c}_p}}{\displaystyle \frac{P_r}{2d}}{\displaystyle \underset{\alpha =1}{\overset{d^{}}{}}}A_\alpha ^{}_{}{}^{}2\mathrm{cos}\left[2\left(k_0r_\alpha +\theta _\alpha ^{}\right)\right],`$ (44)
where $`A_\alpha ^{}`$ and $`\theta _\alpha ^{}`$ are defined by $`A_\alpha ^{}_{}{}^{}2\mathrm{cos}\left[2\left(k_0r_\alpha +\theta _\alpha ^{}\right)\right]=_{\beta =d^{}+1}^dA_{\alpha \beta }^2\mathrm{cos}\left[2\left(k_0r_\alpha +\theta _{\alpha \beta }\right)\right]`$. It can be shown that if $`d3`$, there exist solutions with the constants $`A_{\alpha \beta }`$ and $`\theta _{\alpha \beta }`$ with vanishing $`A_\alpha ^{}_{}{}^{}2`$. However, in the following analysis we will disregard these solutions because, to the best of our knowledge, the asymptotic states corresponding to these solutions have not been reported. Then, we will chose $`\theta _\alpha ^{}`$ so that $`A_\alpha ^{}_{}{}^{}2=_{\beta =d^{}+1}^dA_{\alpha \beta }^2\mathrm{cos}\left[2\left(\theta _{\alpha \beta }\theta _\alpha ^{}\right)\right]`$ becomes positive.

The first approximation of the scaled density and flow field can be obtained from the following relation given by the terms of $`O(ϵ^1)`$:

$`\stackrel{~}{𝐮}^{(0)}\stackrel{~}{n}^{(1)}`$ $`=`$ $`\stackrel{~}{𝐮}^{(1)},`$ (45)
$`\stackrel{~}{𝐮}^{(0)}l_0\stackrel{~}{𝐮}^{(1)}+\stackrel{~}{𝐮}^{(1)}l_0\stackrel{~}{𝐮}^{(0)}`$ (46)
$`=`$ $`{\displaystyle \frac{l_0}{2}}\left[\left({\displaystyle \frac{\stackrel{~}{p}}{\stackrel{~}{n}}}\right)_0\stackrel{~}{n}^{(1)}+\left({\displaystyle \frac{\stackrel{~}{p}}{\stackrel{~}{T}}}\right)_0\stackrel{~}{T}^{(1)}\right]`$
$`+`$ $`D_{}\left[\left(\left({\displaystyle \frac{\eta ^{}}{\stackrel{~}{n}}}\right)_0\stackrel{~}{n}^{(1)}+{\displaystyle \frac{1}{2}}\stackrel{~}{T}^{(1)}\right)\left(\stackrel{~}{𝐮}^{(0)}+(\stackrel{~}{𝐮}^{(0)})^{}\right)\right]`$
$`+`$ $`D_{}\left[\stackrel{~}{𝐮}^{(1)}+(\stackrel{~}{𝐮}^{(1)})^{}{\displaystyle \frac{2}{d}}U(\stackrel{~}{𝐮}^{(1)})\right]+{\displaystyle \frac{\overline{\zeta }}{\overline{\rho }\overline{\omega }}}\stackrel{~}{𝐮}^{(1)}`$
$``$ $`{\displaystyle \frac{1}{2}}\left(_\tau \mathrm{ln}\overline{T}\right)^{(0)}\left(\stackrel{~}{n}^{(1)}\stackrel{~}{𝐮}^{(0)}+\stackrel{~}{𝐮}^{(1)}\right){\displaystyle \frac{1}{2}}\left(_\tau \mathrm{ln}\overline{T}\right)^{(1)}\stackrel{~}{𝐮}^{(0)},`$
where $`(\mathrm{})_0`$ represents that $`(\mathrm{})`$ is evaluated in the homogeneous state, i.e. $`n=\overline{n}`$ and $`T=\overline{T}`$; $`\left(_\tau \mathrm{ln}\overline{T}\right)^{(0)}`$ is given by Eq. (42) and

$`\left(_\tau \mathrm{ln}\overline{T}\right)^{(1)}`$ $`=`$ $`{\displaystyle \frac{1}{V}}{\displaystyle }d𝐫[{\displaystyle \frac{4}{d}}D_{}(\left({\displaystyle \frac{\eta ^{}}{\stackrel{~}{n}}}\right)_0\stackrel{~}{n}^{(1)}+{\displaystyle \frac{1}{2}}\stackrel{~}{T}^{(1)})(\stackrel{~}{𝐮}^{(0)})^{}:\stackrel{~}{𝐮}^{(0)}`$ (47)
$`+`$ $`{\displaystyle \frac{8}{d}}D_{}(\stackrel{~}{𝐮}^{(0)})^{}:\stackrel{~}{𝐮}^{(1)}].`$
These equations can be solved analytically for $`\stackrel{~}{n}^{{}_{}{}^{}(1)}`$ and $`\stackrel{~}{𝐮}^{{}_{}{}^{}(1)}`$ that are spatially averaged over the directions $`\beta `$ ($`\beta =d^{}+1,\mathrm{},d`$), in which $`\stackrel{~}{𝐮}^{(0)}`$ is homogeneous:

$`\stackrel{~}{n}^{{}_{}{}^{}(1)}={\displaystyle \frac{1}{L^{dd^{}}}}{\displaystyle \underset{\beta =d^{}+1}{\overset{d}{}}dr_\beta \stackrel{~}{n}^{(1)}},\stackrel{~}{𝐮}^{{}_{}{}^{}(1)}={\displaystyle \frac{1}{L^{dd^{}}}}{\displaystyle \underset{\beta =d^{}+1}{\overset{d}{}}dr_\beta \stackrel{~}{𝐮}^{(1)}}.`$ (48)
By taking the spatial average over the directions $`\beta `$ ($`\beta =d^{}+1,\mathrm{},d`$), of both sides of Eqs. (45) and (46), we get a closed set of equations for $`\stackrel{~}{n}^{{}_{}{}^{}(1)}`$ and $`\stackrel{~}{𝐮}^{{}_{}{}^{}(1)}`$:

$`0`$ $`=`$ $`{\displaystyle \underset{\alpha =1}{\overset{d^{}}{}}}_\alpha \stackrel{~}{u}_\alpha ^{{}_{}{}^{}(1)}^{}\stackrel{~}{𝐮}^{{}_{}{}^{}(1)},`$ (49)
$`0`$ $`=`$ $`{\displaystyle \frac{l_0}{2}}\left[\left({\displaystyle \frac{\stackrel{~}{p}}{\stackrel{~}{n}}}\right)_0_\alpha \stackrel{~}{n}^{{}_{}{}^{}(1)}+\left({\displaystyle \frac{\stackrel{~}{p}}{\stackrel{~}{T}}}\right)_0_\alpha \stackrel{~}{T}^{(1)}\right]+D_{}^{}_{}{}^{}2\stackrel{~}{u}_\alpha ^{{}_{}{}^{}(1)}+D_{}k_0^2\stackrel{~}{u}_\alpha ^{{}_{}{}^{}(1)},`$ (50)
$`\stackrel{~}{𝐮}^{{}_{}{}^{}(1)}l_0^{}\stackrel{~}{u}_\beta ^{(0)}=`$ (51)
$`+`$ $`D_{}\left[\left(\left({\displaystyle \frac{\eta ^{}}{\stackrel{~}{n}}}\right)_0^{}\stackrel{~}{n}^{{}_{}{}^{}(1)}+{\displaystyle \frac{1}{2}}^{}\stackrel{~}{T}^{(1)}\right)^{}\stackrel{~}{u}_\beta ^{(0)}+\left(\left({\displaystyle \frac{\eta ^{}}{\stackrel{~}{n}}}\right)_0\stackrel{~}{n}^{{}_{}{}^{}(1)}+{\displaystyle \frac{1}{2}}\stackrel{~}{T}^{(1)}\right)^{}_{}{}^{}2\stackrel{~}{u}_\beta ^{(0)}\right]`$
$`+`$ $`D_{}^{}_{}{}^{}2\stackrel{~}{u}_\beta ^{{}_{}{}^{}(1)}+D_{}k_0^2\stackrel{~}{u}_\beta ^{{}_{}{}^{}(1)}+D_{}k_0^2\stackrel{~}{n}^{{}_{}{}^{}(1)}\stackrel{~}{u}_\beta ^{(0)}{\displaystyle \frac{1}{2}}\left(_\tau \mathrm{ln}\overline{T}\right)^{(1)}\stackrel{~}{u}_\beta ^{(0)},`$
where $`\alpha =1,\mathrm{},d^{}`$ and $`\beta =d^{}+1,\mathrm{},d`$, and $`^{}`$ has been defined by $`^{}=(_1,_2,\mathrm{},_d^{})`$. Here Eq. (42) has been used and

$`\left(_\tau \mathrm{ln}\overline{T}\right)^{(1)}`$ $`=`$ $`{\displaystyle \frac{1}{L^d^{}}}{\displaystyle }{\displaystyle \underset{\alpha =1}{\overset{d^{}}{}}}dr_\alpha {\displaystyle \frac{4}{d}}D_{}[(\left({\displaystyle \frac{\eta ^{}}{\stackrel{~}{n}}}\right)_0\stackrel{~}{n}^{{}_{}{}^{}(1)}+{\displaystyle \frac{1}{2}}\stackrel{~}{T}^{(1)})(^{}\stackrel{~}{𝐮}^{(0)})^{}:^{}\stackrel{~}{𝐮}^{(0)}`$
$`+`$ $`2(^{}\stackrel{~}{𝐮}^{(0)})^{}:^{}\stackrel{~}{𝐮}^{{}_{}{}^{}(1)}].`$
(52)
The equation (50) can be solved without any difficulty. Combination of Eqs. (49) and (50) yields

$`0=\left({\displaystyle \frac{\stackrel{~}{p}}{\stackrel{~}{n}}}\right)_0^{}_{}{}^{}2\stackrel{~}{n}^{{}_{}{}^{}(1)}+\left({\displaystyle \frac{\stackrel{~}{p}}{\stackrel{~}{T}}}\right)_0^{}_{}{}^{}2\stackrel{~}{T}^{(1)}.`$ (53)
Because of the periodic boundary conditions and the conditions (31) and (33), we obtain

$`\stackrel{~}{n}^{{}_{}{}^{}(1)}(𝐫)`$ $`=`$ $`{\displaystyle \frac{\left(\frac{\stackrel{~}{p}}{\stackrel{~}{T}}\right)_0}{\left(\frac{\stackrel{~}{p}}{\stackrel{~}{n}}\right)_0}}\stackrel{~}{T}^{(1)}(𝐫)={\displaystyle \frac{\frac{1}{\overline{n}}\left(\frac{p}{T}\right)_0}{\frac{1}{\overline{T}}\left(\frac{p}{n}\right)_0}}\stackrel{~}{T}^{(1)}(𝐫)=\overline{\beta }_0\stackrel{~}{T}^{(1)}(𝐫).`$ (54)
This relation is consistent with the pressure balance achieved in the asymptotic state. Substituting Eqs. (44) and (54) into Eq. (50), we get

$`0=D_{}^{}_{}{}^{}2\stackrel{~}{u}_\alpha ^{{}_{}{}^{}(1)}+D_{}k_0^2\stackrel{~}{u}_\alpha ^{{}_{}{}^{}(1)}.`$ (55)
This relation suggests in terms of the Fourier modes $`\stackrel{~}{u}_{\alpha 𝐤}^{{}_{}{}^{}(1)}=0`$ if $`|𝐤|k_0`$; hence $`\stackrel{~}{u}_\alpha ^{{}_{}{}^{}(1)}`$ ($`\alpha =1,\mathrm{},d^{}`$) is given by

$`\stackrel{~}{u}_\alpha ^{{}_{}{}^{}(1)}(𝐫)={\displaystyle \underset{\alpha ^{}=1}{\overset{d^{}}{}}}C_{\alpha ^{}\alpha }\mathrm{cos}(k_0r_\alpha ^{}+\zeta _{\alpha ^{}\alpha }),`$ (56)
where the amplitudes $`C_{\alpha ^{}\alpha }`$ are real numbers that satisfy $`C_{\alpha \alpha }=0`$ because of Eq. (49), and the phase factors $`\zeta _{\alpha ^{}\alpha }`$ are also real.

Since $`\stackrel{~}{T}^{(1)}`$, $`\stackrel{~}{n}^{{}_{}{}^{}(1)}`$ and $`\stackrel{~}{u}_\alpha ^{{}_{}{}^{}(1)}`$ ($`\alpha =1,\mathrm{},d^{}`$) are now known, Eq. (51) is an inhomogeneous linear equation for $`\stackrel{~}{u}_\beta ^{{}_{}{}^{}(1)}`$ ($`\beta =d^{}+1,\mathrm{},d`$). We will show in Appendix D that the solubility conditions of Eq. (51) is satisfied if and only if

$`\theta _{\alpha \beta }`$ $`=`$ $`\theta _\alpha ^{}\text{or}\theta _\alpha ^{}+\pi ,{\displaystyle \underset{\beta =d^{}+1}{\overset{d}{}}}A_{\alpha \beta }^2=A_\alpha ^{}_{}{}^{}2={\displaystyle \frac{1}{d^{}}}A_0^2,`$ (57)
for $`\alpha =1,\mathrm{},d^{}`$ and $`\beta =d^{}+1,\mathrm{},d`$. The freedom of the choice $`\theta _{\alpha \beta }=\theta _\alpha ^{}\text{ or }\theta _\alpha ^{}+\pi `$ can be absorbed into the definition of $`A_{\alpha \beta }`$; hence Eqs. (40) and (44) can be rewritten as

$`\stackrel{~}{𝐮}^{(0)}(𝐫)`$ $`=`$ $`{\displaystyle \underset{\alpha =1}{\overset{d^{}}{}}}{\displaystyle \underset{\beta =d^{}+1}{\overset{d}{}}}A_{\alpha \beta }\widehat{𝐞}_\beta \mathrm{cos}(k_0r_\alpha +\theta _\alpha ^{}),`$ (58)
$`\stackrel{~}{T}^{(1)}(𝐫)`$ $`=`$ $`{\displaystyle \frac{\widehat{T}}{d^{}}}{\displaystyle \underset{\alpha =1}{\overset{d^{}}{}}}\mathrm{cos}\left[2\left(k_0r_\alpha +\theta _\alpha ^{}\right)\right],`$ (59)
where $`\widehat{T}(c_v/\overline{c}_p)(P_r/2d)A_0^2`$. According to the relation (54), $`\stackrel{~}{n}^{{}_{}{}^{}(1)}`$ is given by

$`\stackrel{~}{n}^{{}_{}{}^{}(1)}(𝐫)`$ $`=`$ $`{\displaystyle \frac{\widehat{n}}{d^{}}}{\displaystyle \underset{\alpha =1}{\overset{d^{}}{}}}\mathrm{cos}\left[2\left(k_0r_\alpha +\theta _\alpha ^{}\right)\right],`$ (60)
where $`\widehat{n}=\overline{\beta }_0\widehat{T}`$. Once the solubility conditions are satisfied, calculation of the solution of Eq. (51) which satisfies the condition (32) is straightforward. The final result is presented in Appendix E.

Finally, substituting Eqs. (58), (59), (60) and (A.35) into (52), we obtain the first order correction to $`(_\tau \mathrm{ln}\overline{T})`$,

$`\left(_\tau \mathrm{ln}\overline{T}\right)^{(1)}`$ $`=`$ $`D_{}k_0^2{\displaystyle \frac{1}{d^{}}}\left[\left(1+\left({\displaystyle \frac{\eta ^{}}{\stackrel{~}{n}}}\right)_0\right)\widehat{n}+{\displaystyle \frac{1}{2}}\widehat{T}\right].`$ (61)
Hence, from Eqs. (42) and (61), the exponent $`2\gamma _a`$ in Eq. (15) is given by

$`2\gamma _a`$ $`=`$ $`\left[\left(_\tau \mathrm{ln}\overline{T}\right)^{(0)}+\left(_\tau \mathrm{ln}\overline{T}\right)^{(1)}\right]`$ (62)
$`=`$ $`2D_{}k_0^2{\displaystyle \frac{c_v}{\overline{c}_p}}{\displaystyle \frac{P_r}{2d^{}}}\left(\gamma _0D_{}k_0^2\right)\left[\left(1+\left({\displaystyle \frac{\eta ^{}}{\stackrel{~}{n}}}\right)_0\right)\overline{\beta }_0{\displaystyle \frac{1}{2}}\right].`$
Notice that the first order correction to $`2\gamma _a`$ depends explicitly on the coefficient of restitution $`\alpha `$ through $`\gamma _0`$. The ratio of the energy per particle to the global temperature is independent of time and can be calculated from Eq. (21) to order $`ϵ`$:

$`{\displaystyle \frac{E}{\frac{d}{2}\overline{T}}}`$ $`=`$ $`{\displaystyle \frac{2}{d}}{\displaystyle \frac{1}{V}}{\displaystyle 𝑑𝐫\left(\stackrel{~}{𝐮}^{(0)2}+2\stackrel{~}{𝐮}^{(0)}\stackrel{~}{𝐮}^{(1)}+\stackrel{~}{n}^{(1)}\stackrel{~}{𝐮}^{(0)2}\right)}+1`$ (63)
$`=`$ $`{\displaystyle \frac{\gamma _0}{D_{}k_0^2}}\left\{1+{\displaystyle \frac{c_v}{\overline{c}_p}}{\displaystyle \frac{P_r}{4d^{}}}{\displaystyle \frac{\left(\gamma _0D_{}k_0^2\right)}{D_{}k_0^2}}\left[\left(1+\left({\displaystyle \frac{\eta ^{}}{\stackrel{~}{n}}}\right)_0\right)\overline{\beta }_0{\displaystyle \frac{1}{2}}\right]\right\},`$
where Eqs. (58), (60) and (A.35) have been used.

Calculation in the higher order approximation is rather involved. In Appendix F we discuss the second order correction to the scaled temperature, $`\stackrel{~}{T}^{(2)}`$, in order to illustrate how excitation of modes with larger wave numbers can be described in the higher order approximations.

### 3.4 The range of validity of the expansion method

A restriction to the range of validity of the expansion method is given by the condition that the amplitude of the temperature inhomogeneities $`|\widehat{T}|=(c_v/\overline{c}_p)(P_r/2d)A_0^2=(c_v/\overline{c}_p)(P_r/d)\stackrel{~}{U}^2`$, which is greater than the amplitude of the density inhomogeneities $`|\widehat{n}|`$ in two and three dimensions, has to be much smaller than $`1`$. In Sec.3.1 we have seen that $`(c_v/\overline{c}_p)(P_r/d)`$ is a small parameter; therefore the above condition is fulfilled if

$`\stackrel{~}{U}^2={\displaystyle \frac{A_0^2}{2}}={\displaystyle \frac{d}{2}}{\displaystyle \frac{\gamma _0D_{}k_0^2}{D_{}k_0^2}}\stackrel{<}{}1.`$ (64)
This condition can be rewritten as

$`\gamma _0\stackrel{<}{}\left(1+{\displaystyle \frac{2}{d}}\right)D_{}k_0^2,`$ (65)
and hence tells how small the system has to be in order that our systematic expansion method describes well the asymptotic state.

We have also seen in Sec.3.1 that the inhomogeneous part of the collisional dissipation term, which contributes to reduction of inhomogeneities in the temperature field, is negligible if the condition (29) is fulfilled. From the condition (65), however, we see that if $`\stackrel{~}{U}^2\stackrel{<}{}1`$ then the condition (29) is fulfilled. This result is consistent with the fact that the inhomogeneous part of the collisional dissipation term appears in the second order correction to the scaled temperature, as shown in Appendix F. Therefore the range of validity of our systematic expansion method applied to a solution with the properties (i)-(v) is given only by the condition $`\stackrel{~}{U}^2\stackrel{<}{}1`$.

## 4 Conclusion

In this paper, a systematic expansion method in the Prandtl number has been developed in order to describe the asymptotic state of freely cooling granular fluids contained in sufficiently small systems. Although the Prandtl number of IHS fluids itself is not a very small parameter, it eventually leads an expansion in the small parameter $`(c_v/\overline{c}_p)(P_r/d)`$, which characterizes the size of temperature and density inhomogeneities if a typical size of the scaled flow field $`\stackrel{~}{U}=A_0/\sqrt{2}`$ satisfies $`\stackrel{~}{U}^2\stackrel{<}{}1`$. The expansion in the Prandtl number yields a set of inhomogeneous linear equations. The zeroth approximation is a homogeneous state where only the scaled flow field is allowed to have finite amplitude. The resulting equations in the zeroth approximation are in agreement with those equations for the stationary state which are given in Ref.  with an approximation in which density and temperature inhomogeneities are neglected. It has been shown that inhomogeneities in density and temperature are systematically described in the higher order approximation; corrections to macroscopic quantities such as the energy per particle and the global temperature caused by these inhomogeneities are estimated up to the first approximation.

Our theory describes sufficiently small systems where thermal conduction plays a more important role than convection and collisional dissipation. As a consequence, the temperature inhomogeneities is determined mainly by viscous heating induced by the flow profile, and thermal conduction. Density profile is then determined in such a way that pressure averaged over the directions of the flow in the zeroth approximation remains homogeneous. If the size of system is so large that thermal conduction as a diffusion process is less important than collisional dissipation for large scale fluctuations of temperature, a nonlinear theory by Goldhirsch et al.  is more appropriate to describe temperature profile determined mainly by viscous heating and collisional dissipation; their theory describes also pressure imbalance that eventually leads a clustering process. Apparently the clustering process can not be described by our method because the condition of small inhomogeneities in density and temperature is necessary for the systematic expansion method to work well.

The condition that the amplitude of temperature inhomogeneities $`|\widehat{T}|=(c_v/\overline{c}_p)(P_r/2d)A_0^2=(c_v/\overline{c}_p)(P_r/d)\stackrel{~}{U}^2`$ has to be much smaller than $`1`$ yields a restriction on the system size $`L`$ for a given packing fraction and the coefficient of restitution $`\alpha `$ of IHS. Because $`(c_v/\overline{c}_p)(P_r/d)1`$ in two and three dimensions \[see Table 1\], our theory is expected to work well for $`\stackrel{~}{U}^2\stackrel{<}{}1`$; this is the case that can not be described by those amplitude equations in Ref.  which are derived assuming $`\stackrel{~}{U}^21`$. One can show that our results expanded in the small parameter in Ref. , $`\stackrel{~}{U}^2=A_0^2/2`$, coincide to order $`A_0^2P_r^2`$ with the results in Ref.  expanded in the Prandtl number $`P_r`$, except for the indefinite constants in our theory, which can not be fixed in the first approximation.

In this paper, stability of the stationary solutions has not been discussed. We have assumed the flow profiles that are observed in computer simulations (the property (i) in Sec.1) and concentrated on the results that can be deduced from this assumption, such as the amplitude of the scaled flow field, the amplitude of the density and temperature inhomogeneities, and the decay exponent of the global temperature and the energy per particle. The stability analysis may eliminate solutions that are not observed in computer simulations and therefore not discussed in this paper.

In this paper the systematic expansion method has been applied to freely evolving granular fluids contained in small systems with periodic boundary conditions; these systems have been studied in the earlier works by computer simulations and by theories. The asymptotic state of these systems is unphysical because it is dominated by the finite size effect due to periodic boundary conditions. Our method of systematic expansion, however, is expected to be useful to study other practical problems which occur in finite size systems with physical boundary conditions, such as Couette flow and flow under vibration, if the dimensions of the system are sufficiently small such that viscous heating and thermal conduction play dominant roles in formation of density and temperature inhomogeneities.

## 5 Acknowledgements

I thank M. H. Ernst, R. Brito, R. Soto and M. Mareschal for helpful discussions. Thanks are due to M. H. Ernst for reading the manuscript and making useful suggestions. I acknowledges support of the foundation ”Fundamenteel Onderzoek der Materie (FOM)”, which is financially supported by the Dutch National Science Foundation (NWO).

## Appendix A Bulk expansion coefficient for inelastic hard spheres

In this appendix, we present the expression for the bulk expansion coefficient $`\beta _B`$ for $`d`$-dimensional IHS fluids and show that $`\beta _0T\beta _B<1`$ in two and three dimensions. The bulk expansion coefficient $`\beta _B`$ is defined by

$`\beta _B={\displaystyle \frac{1}{n}}{\displaystyle \frac{n}{T}}|_p={\displaystyle \frac{1}{n}}{\displaystyle \frac{\frac{p}{T}|_n}{\frac{p}{n}|_T}}.`$ (A.1)
The expression for the local pressure in Enskog’s kinetic theory has been given in Refs.  for three-dimensional IHS fluids. A generalization to $`d`$ dimensions has been given in Ref. :

$`p=nT\left(1+{\displaystyle \frac{1+\alpha }{2}}\chi n\sigma ^d{\displaystyle \frac{\mathrm{\Omega }_d}{2d}}\right),`$ (A.2)
where $`\mathrm{\Omega }_d=2\pi ^{d/2}/\mathrm{\Gamma }(d/2)`$ is the surface area of a $`d`$-dimensional unit sphere, $`\sigma `$ is the diameter of the hard sphere, and $`\chi `$ is the pair correlation function for two particles in contact. In the Enskog approximation $`\chi `$ has the same dependence on the local density at the point of contact as the equilibrium pair correlation function at contact has on the equilibrium density. Substituting the expression (A.2) into Eq. (A.1), we get

$`\beta _B={\displaystyle \frac{1}{T}}{\displaystyle \frac{1+\frac{1+\alpha }{2}\chi n\sigma ^d\frac{\mathrm{\Omega }_d}{2d}}{1+\frac{1+\alpha }{2}\chi n\left(1+\frac{1}{\chi }\frac{}{n}\left(\chi n\right)\right)\sigma ^d\frac{\mathrm{\Omega }_d}{2d}}}.`$ (A.3)
From this expression, we see that $`\beta _0=T\beta _B`$ is independent of T, and if

$`{\displaystyle \frac{}{n}}\left(\chi n\right)>0,`$ (A.4)
then $`\beta _0<1`$. The condition (A.4) is satisfied in three dimensions for $`\chi `$ given by the Carnahan-Starling approximation , $`\chi =(2\varphi )/2(1\varphi )^3`$, and in two dimensions for $`\chi `$ given by the Verlet-Levesque approximation , $`\chi =(17\varphi /16)/(1\varphi )^2`$, where $`\varphi =n\mathrm{\Omega }_d(\sigma /2)^d/d`$ is the packing fraction in $`d`$ dimensions.

## Appendix B Homogeneous cooling state

In this appendix, we show that if $`k>k_{}^{}`$, Eq. (38) has only the trivial solution $`\stackrel{~}{𝐮}^{(0)}=\mathrm{𝟎}`$. We first take the inner product of $`\stackrel{~}{𝐮}^{(0)}`$ with Eq. (38) and integrate both sides over the whole space:

$`{\displaystyle 𝑑𝐫\stackrel{~}{𝐮}^{(0)}\left(\stackrel{~}{𝐮}^{(0)}l_0\stackrel{~}{𝐮}^{(0)}\right)}`$ $`=`$ $`D_{}{\displaystyle 𝑑𝐫\stackrel{~}{𝐮}^{(0)}^2\stackrel{~}{𝐮}^{(0)}}{\displaystyle \frac{1}{2}}\left(_\tau \mathrm{ln}\overline{T}\right)^{(0)}{\displaystyle 𝑑𝐫(\stackrel{~}{𝐮}^{(0)})^2}.`$ (A.5)
Using the relation $`\stackrel{~}{𝐮}^{(0)}(\stackrel{~}{𝐮}^{(0)}\stackrel{~}{𝐮}^{(0)})=\stackrel{~}{𝐮}^{(0)}(\stackrel{~}{𝐮}^{(0)})^2/2`$ and performing an integration by parts, it can be shown that the left hand side of Eq. (A.5) vanishes as a consequence of Eq. (37) and the periodic boundary conditions. Rewriting the right hand side in terms of Fourier modes defined as $`a_𝐤=𝑑𝐫e^{i𝐤𝐫}a(𝐫)`$, we get

$`0={\displaystyle \underset{𝐤}{}}\left(\gamma _0D_{}k^2{\displaystyle \frac{2}{d}}D_{}{\displaystyle \frac{1}{V^2}}{\displaystyle \underset{𝐤_1}{}}k_1^2|\stackrel{~}{𝐮}_{𝐤_1}^{(0)}|^2\right)|\stackrel{~}{𝐮}_𝐤^{(0)}|^2.`$ (A.6)
We notice that if $`k_0>k_{}^{}\sqrt{\gamma _0/D_{}}`$, the inside of the parenthesis on the right hand side of Eq. (A.6) is negative for any $`𝐤`$. Therefore, the equality is satisfied if and only if $`\stackrel{~}{𝐮}^{(0)}=\mathrm{𝟎}`$.

## Appendix C Flow profile given in (40)

The flow field $`\stackrel{~}{𝐮}^{(0)}`$ that consists only of the Fourier modes with the smallest wave number $`k_0`$ is written in general as

$`\stackrel{~}{𝐮}^{(0)}(𝐫)={\displaystyle \underset{\alpha =1}{\overset{d}{}}}{\displaystyle \underset{\beta =1}{\overset{d}{}}}A_{\alpha \beta }\widehat{𝐞}_\beta \mathrm{cos}(k_0r_\alpha +\theta _{\alpha \beta }),`$ (A.7)
where the coefficients $`A_{\alpha \beta }`$ and $`\theta _{\alpha \beta }`$ are arbitrary real constants. From Eq. (37), we get $`A_{\alpha \alpha }=0`$.

For the flow field given by Eq. (A.7), both sides of Eq. (38) have to vanish independently of each other, because the right hand side of Eq. (38) consists only of the Fourier modes with the smallest wave number $`k_0`$, while the left hand side of Eq. (38) does not contain the smallest wave number modes. We first show that the condition of the vanishing left hand side of Eq. (38) leads to the flow field (40). The condition of the vanishing left hand side can be written in terms of the Fourier modes,

$`\mathrm{𝟎}`$ $`=`$ $`{\displaystyle \underset{𝐤_1}{}}\stackrel{~}{𝐮}_{𝐤𝐤_1}^{(0)}𝐤_1\stackrel{~}{𝐮}_{𝐤_1}^{(0)},`$ (A.8)
for any $`𝐤`$. Because $`\stackrel{~}{𝐮}_𝐤^{(0)}=\mathrm{𝟎}`$ if $`|𝐤|k_0`$, the non-trivial conditions are given only for the wave numbers $`𝐤`$ that satisfy $`|𝐤𝐤_{0\alpha }|=k_0`$ $`(\alpha =1,2,\mathrm{},d)`$, or, equivalently, only for $`𝐤=𝐤_{0\alpha }\pm 𝐤_{0\beta }`$ $`(\alpha =1,2,\mathrm{},d,\beta =1,2,\mathrm{},d)`$. Here $`𝐤_{0\alpha }=k_0\widehat{𝐞}_\alpha `$ and $`\widehat{𝐞}_\alpha `$ is a unit vector in the direction $`\alpha `$. For a given $`𝐤=𝐤_{0\alpha }+𝐤_{0\beta }`$, the condition (A.8) is written as

$`\mathrm{𝟎}`$ $`=`$ $`\stackrel{~}{𝐮}_{𝐤_{0\alpha }}^{(0)}𝐤_{0\beta }\stackrel{~}{𝐮}_{𝐤_{0\beta }}^{(0)}+\stackrel{~}{𝐮}_{𝐤_{0\beta }}^{(0)}𝐤_{0\alpha }\stackrel{~}{𝐮}_{𝐤_{0\alpha }}^{(0)},`$ (A.9)
$`\mathrm{𝟎}`$ $`=`$ $`\stackrel{~}{𝐮}_{𝐤_{0\alpha }}^{(0)}𝐤_{0\beta }\stackrel{~}{𝐮}_{𝐤_{0\beta }}^{(0)}+\stackrel{~}{𝐮}_{𝐤_{0\beta }}^{(0)}𝐤_{0\alpha }\stackrel{~}{𝐮}_{𝐤_{0\alpha }}^{(0)}.`$ (A.10)
Substituting the Fourier transform of Eq. (A.7), $`\stackrel{~}{𝐮}_{𝐤_{0\alpha }}=(V/2)_{\beta =1}^d\widehat{𝐞}_\beta B_{\alpha \beta }`$, where $`B_{\alpha \beta }=A_{\alpha \beta }\mathrm{exp}(i\theta _{\alpha \beta })`$, into the above condition (A.10), we obtain

$`0=B_{\alpha \beta }B_{\beta \gamma }+B_{\beta \alpha }B_{\alpha \gamma },`$ (A.11)
$`0=B_{\alpha \beta }B_{\beta \gamma }^{}B_{\beta \alpha }^{}B_{\alpha \gamma },`$ (A.12)
for any set of $`(\alpha ,\beta ,\gamma )`$.

This conditions imply the following properties of the complex matrix $`B_{\alpha \beta }`$. First, considering the case when $`\gamma =\alpha `$ in Eq. (A.11), we get

$`0=B_{\alpha \beta }B_{\beta \alpha },`$ (A.13)
for any $`\alpha `$ and $`\beta `$. Here we have used property that $`B_{\alpha \alpha }=0`$. Second, suppose $`B_{\alpha \beta }0`$. Then, $`B_{\beta \alpha }=0`$ from Eq. (A.13). Then, the condition (A.11) is reduced to $`0=B_{\alpha \beta }B_{\beta \gamma }`$. Because $`B_{\alpha \beta }0`$, we get $`B_{\beta \gamma }=0`$ for any $`\gamma `$. Then suppose $`A_{\alpha \beta }0`$ for a set of $`(\alpha ,\beta )=(\mu _1,\nu _1),(\mu _2,\nu _2),\mathrm{},(\mu _m,\nu _m)`$ and $`A_{\alpha \beta }=0`$ otherwise. Let us label the directions in the subspace spanned by $`(\mu _1,\mu _2,\mathrm{},\mu _m)`$ as $`(1,2,\mathrm{},d^{})`$, where $`d^{}`$ is the number of dimensions of the subspace. Note that $`d^{}m`$ because some of $`\mu _i`$ may be the same direction. Then one finds that any $`\nu _i`$ can be equal to none of $`(1,2,\mathrm{},d^{})`$, because $`B_{\nu _i\gamma }=A_{\nu _i\gamma }=0`$ for any $`\gamma `$. The condition (A.12) gives no additional conditions. Hence, we may conclude that the direction of the flow, which is in the subspace spanned by $`(\nu _1,\nu _2,\mathrm{},\nu _m)`$ has to be perpendicular to the directions $`(1,2,\mathrm{},d^{})`$ in which the flow shows inhomogeneities. Therefore the flow can be in general written as Eq. (40).

Then it is straightforward to get the condition (41) by substituting the expression (40) into the condition of the vanishing right hand side of Eq. (38).

## Appendix D Solubility conditions of Eq. (51)

Equation (51) is the inhomogeneous linear equation with respect to $`\stackrel{~}{u}_\beta ^{{}_{}{}^{}(1)}`$ and rewritten as

$`\stackrel{~}{u}_\beta ^{{}_{}{}^{}(1)}=f_\beta ,`$ (A.14)
where

$`\stackrel{~}{u}_\beta ^{{}_{}{}^{}(1)}`$ $`=`$ $`D_{}^{}_{}{}^{}2\stackrel{~}{u}_\beta ^{{}_{}{}^{}(1)}+D_{}k_0^2\stackrel{~}{u}_\beta ^{{}_{}{}^{}(1)}k_0^2\left[{\displaystyle \frac{1}{L^d^{}}}{\displaystyle \underset{\alpha =1}{\overset{d^{}}{}}dr_\alpha \stackrel{~}{𝐮}^{(0)}\stackrel{~}{𝐮}^{{}_{}{}^{}(1)}}\right]\stackrel{~}{u}_\beta ^{(0)},`$ (A.15)
and

$`f_\beta `$ $`=`$ $`\stackrel{~}{𝐮}^{{}_{}{}^{}(1)}l_0^{}\stackrel{~}{u}_\beta ^{(0)}D_{}k_0^2\stackrel{~}{n}^{{}_{}{}^{}(1)}\stackrel{~}{u}_\beta ^{(0)}`$ (A.16)
$``$ $`D_{}\left[\left(\left({\displaystyle \frac{\eta ^{}}{\stackrel{~}{n}}}\right)_0^{}\stackrel{~}{n}^{{}_{}{}^{}(1)}+{\displaystyle \frac{1}{2}}^{}\stackrel{~}{T}^{(1)}\right)^{}\stackrel{~}{u}_\beta ^{(0)}+\left(\left({\displaystyle \frac{\eta ^{}}{\stackrel{~}{n}}}\right)_0\stackrel{~}{n}^{{}_{}{}^{}(1)}+{\displaystyle \frac{1}{2}}\stackrel{~}{T}^{(1)}\right)^{}_{}{}^{}2\stackrel{~}{u}_\beta ^{(0)}\right]`$
$`+`$ $`{\displaystyle \frac{2}{d}}D_{}[{\displaystyle \frac{1}{L^d^{}}}{\displaystyle }{\displaystyle \underset{\alpha =1}{\overset{d^{}}{}}}dr_\alpha (\left({\displaystyle \frac{\eta ^{}}{\stackrel{~}{n}}}\right)_0\stackrel{~}{n}^{{}_{}{}^{}(1)}+{\displaystyle \frac{1}{2}}\stackrel{~}{T}^{(1)})(^{}\stackrel{~}{𝐮}^{(0)})^{}:^{}\stackrel{~}{𝐮}^{(0)}]\stackrel{~}{u}_\beta ^{(0)},`$
where $`\stackrel{~}{𝐮}^{(0)}`$, $`\stackrel{~}{T}^{(1)}`$, $`\stackrel{~}{n}^{(1)}`$, and $`\stackrel{~}{𝐮}^{(1)}`$ are given by Eqs. (40),(44),(54),(A.35) and (56). Let us define a Hilbert space $`|a`$ to represent any function $`a_\beta (𝐫)`$ ($`\beta =d^{}+1,\mathrm{},d,𝐫=(r_1,r_2,\mathrm{},r_d^{}))`$ by

$`a_\beta (𝐫)=𝐫,\beta |a.`$ (A.17)
We define the scalar product of two vectors $`|a`$ and $`|b`$ by

$`a|b={\displaystyle \frac{1}{L^d^{}}}{\displaystyle \underset{\alpha =1}{\overset{d^{}}{}}dr_\alpha \underset{\beta =d^{}+1}{\overset{d}{}}a_\beta ^{}(𝐫)b_\beta (𝐫)}.`$ (A.18)
With this definition, the inhomogeneous linear equation (A.14) can be written as

$`|\stackrel{~}{u}^{{}_{}{}^{}(1)}>=|f,`$ (A.19)
and the linear operator $``$ is Hermitian:

$`a|b=b|a.`$ (A.20)
The inhomogeneous linear equation (A.19) has a nontrivial solution if and only if $`|f`$ is orthogonal to the solutions $`|\varphi ^l>`$ of the homogeneous problem $`|\varphi ^l>=0`$ , i.e.

$`<\varphi ^l|f>=0.`$ (A.21)
One can show that the linearly independent solutions $`|\varphi ^l>`$ of the homogeneous problem are given by

$`\varphi _{c\beta }^l(𝐫)`$ $`=`$ $`<𝐫,\beta |\varphi _c^l>={\displaystyle \underset{\alpha =1}{\overset{d^{}}{}}}c_{\alpha \beta }^l\mathrm{cos}(k_0r_\alpha +\theta _{\alpha \beta }),`$ (A.22)
$`\varphi _{s\beta }^m(𝐫)`$ $`=`$ $`<𝐫,\beta |\varphi _s^m>={\displaystyle \underset{\alpha =1}{\overset{d^{}}{}}}s_{\alpha \beta }^m\mathrm{sin}(k_0r_\alpha +\theta _{\alpha \beta }),`$ (A.23)
where the coefficients $`c_{\alpha \beta }^l`$ $`(l=1,2,\mathrm{},d^{}(dd^{})1)`$ and $`s_{\alpha \beta }^m`$ $`(m=1,2,\mathrm{},d^{}(dd^{}))`$ are real numbers. The coefficients $`c_{\alpha \beta }^l`$ satisfy the condition

$`{\displaystyle \underset{\alpha =1}{\overset{d^{}}{}}}{\displaystyle \underset{\beta =d^{}+1}{\overset{d}{}}}A_{\alpha \beta }c_{\alpha \beta }^l=0.`$ (A.24)
We also require the orthogonality property

$`<\varphi _c^l|\varphi _c^l^{}>`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{\alpha =1}{\overset{d^{}}{}}}{\displaystyle \underset{\beta =d^{}+1}{\overset{d}{}}}c_{\alpha \beta }^lc_{\alpha \beta }^l^{}=\delta _{ll^{}},`$ (A.25)
$`<\varphi _s^m|\varphi _s^m^{}>`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{\alpha =1}{\overset{d^{}}{}}}{\displaystyle \underset{\beta =d^{}+1}{\overset{d}{}}}s_{\alpha \beta }^ms_{\alpha \beta }^m^{}=\delta _{mm^{}},`$ (A.26)
$`<\varphi _c^l|\varphi _s^m>`$ $`=`$ $`<\varphi _s^m|\varphi _c^l>=0.`$ (A.27)
The solubility conditions (A.21), $`<\varphi _c^l|f>=0`$ and $`<\varphi _s^m|f>=0`$, yield

$`{\displaystyle \underset{\alpha =1}{\overset{d^{}}{}}}{\displaystyle \underset{\beta =d^{}+1}{\overset{d}{}}}A_\alpha ^{}_{}{}^{}2A_{\alpha \beta }c_{\alpha \beta }^l\mathrm{cos}[2(\theta _\alpha ^{}\theta _{\alpha \beta })]`$ $`=`$ $`0,`$ (A.28)
and

$`{\displaystyle \underset{\alpha =1}{\overset{d^{}}{}}}{\displaystyle \underset{\beta =d^{}+1}{\overset{d}{}}}A_\alpha ^{}_{}{}^{}2A_{\alpha \beta }s_{\alpha \beta }^m\mathrm{sin}[2(\theta _\alpha ^{}\theta _{\alpha \beta })]`$ $`=`$ $`0,`$ (A.29)
respectively. Here we have assumed $`[(1+(\eta ^{}/\stackrel{~}{n})_0)\overline{\beta }_01/2]0`$, which is satisfied in general. In deriving Eq. (A.29) we have used the fact that $`\stackrel{~}{u}_{\alpha 𝐤}^{(1)}=0`$ $`(\alpha =1,\mathrm{}d^{})`$ if $`|𝐤|k_0`$ as shown in Eq. (56).

Because the coefficient $`s_{\alpha \beta }^1`$ is arbitrary, let us choose it as

$`s_{\alpha \beta }^1A_{\alpha \beta }\mathrm{sin}[2(\theta ^{}\theta _{\alpha \beta })].`$ (A.30)
Substituting this expression into Eq. (A.29), we get

$`{\displaystyle \underset{\alpha =1}{\overset{d^{}}{}}}{\displaystyle \underset{\beta =d^{}+1}{\overset{d}{}}}A_\alpha ^{}_{}{}^{}2A_{\alpha \beta }^2\mathrm{sin}^2[2(\theta _\alpha ^{}\theta _{\alpha \beta })]`$ $`=`$ $`0.`$ (A.31)
Because $`A_\alpha ^{}_{}{}^{}2>0`$ by definition \[see note below Eq. (44)\], the condition (A.31) suggests that if $`A_{\alpha \beta }0`$, then $`(\theta _{\alpha \beta }\theta _\alpha ^{})`$ is either $`0,\pi /2,\pi `$, or $`3\pi /2`$. We note that if $`A_{\alpha \beta }=0`$ then the phase factor $`\theta _{\alpha \beta }`$ does not appear in our theory \[see, Eq. (40)\]. If the condition (A.31) is satisfied, the conditions (A.29) for $`s_{\alpha \beta }^m`$ ($`m=2,3,\mathrm{},d^{}(dd^{})`$) are automatically satisfied.

The condition (A.28) has to be satisfied for a set of constants $`c_{\alpha \beta }^l`$, which satisfies Eqs. (A.24) and (A.25). This implies that $`A_\alpha ^{}_{}{}^{}2A_{\alpha \beta }\mathrm{cos}[2(\theta _\alpha ^{}\theta _{\alpha \beta })]=\text{const.}\times A_{\alpha \beta }`$, i.e.

$`A_\alpha ^{}_{}{}^{}2\mathrm{cos}[2(\theta _\alpha ^{}\theta _{\alpha \beta })]=\text{const.}=𝒵,`$ (A.32)
where the constant $`𝒵`$ is independent of $`\alpha `$ and $`\beta `$. The condition that the left hand side of Eq. (A.32) is independent of $`\beta `$ suggests that $`(\theta _{\alpha \beta }\theta _\alpha ^{})`$ is either $`0`$ or $`\pi `$, or $`(\theta _{\alpha \beta }\theta _\alpha ^{})`$ is either $`\pi /2`$ or $`3\pi /2`$. Because of the definition of $`\theta _\alpha ^{}`$ \[see note below Eq. (44)\], however, only the former choice

$`\theta _{\alpha \beta }\theta _\alpha ^{}=0\text{ or }\pi `$ (A.33)
is allowed. Then, from Eqs. (A.32) and (A.33) we get $`A_\alpha ^{}_{}{}^{}2=_{\beta =d^{}+1}^dA_{\alpha \beta }^2=𝒵`$. Finally, the constant $`𝒵`$ is determined by the condition (41),

$`A_0^2={\displaystyle \underset{\alpha =1}{\overset{d^{}}{}}}A_\alpha ^{}_{}{}^{}2=𝒵d^{}.`$ (A.34)
## Appendix E First order correction $`\stackrel{~}{u}_\beta ^{{}_{}{}^{}(1)}`$

First order correction $`\stackrel{~}{u}_\beta ^{{}_{}{}^{}(1)}`$ can be obtained by solving Eq. (51) under the solubility conditions (57) and the condition (32):

$`\stackrel{~}{u}_\beta ^{{}_{}{}^{}(1)}(𝐫)=𝒰_1{\displaystyle \underset{\alpha =1}{\overset{d^{}}{}}}A_{\alpha \beta }\mathrm{cos}[3(k_0r_\alpha +\theta _\alpha ^{})]`$ (A.35)
$`+`$ $`𝒰_2{\displaystyle \underset{\alpha =1}{\overset{d^{}}{}}}{\displaystyle \underset{\genfrac{}{}{0pt}{}{\alpha ^{}=1}{\alpha ^{}\alpha }}{\overset{d^{}}{}}}A_{\alpha ^{}\beta }\left\{\mathrm{cos}\left[k_0(2r_\alpha r_\alpha ^{})+2\theta _\alpha ^{}\theta _\alpha ^{}^{}\right]+\mathrm{cos}\left[k_0(2r_\alpha +r_\alpha ^{})+2\theta _\alpha ^{}+\theta _\alpha ^{}^{}\right]\right\}`$
$`+`$ $`{\displaystyle \frac{l_0}{D_{}k_0}}{\displaystyle \underset{\alpha =1}{\overset{d^{}}{}}}A_{\alpha \beta }\stackrel{~}{u}_\alpha ^{{}_{}{}^{}(1)}(𝐫)\mathrm{sin}(k_0r_\alpha +\theta _\alpha ^{})`$
$`+`$ $`{\displaystyle \underset{\alpha =1}{\overset{d^{}}{}}}E_{\alpha \beta }\mathrm{cos}(k_0r_\alpha +\theta _\alpha ^{})+{\displaystyle \underset{\alpha =1}{\overset{d^{}}{}}}F_{\alpha \beta }\mathrm{sin}(k_0r_\alpha +\theta _\alpha ^{}),`$
where the coefficients $`E_{\alpha \beta }`$ are real numbers that satisfy the condition

$`{\displaystyle \underset{\alpha =1}{\overset{d^{}}{}}}{\displaystyle \underset{\beta =d^{}+1}{\overset{d}{}}}A_{\alpha \beta }E_{\alpha \beta }=𝒰_3,`$ (A.36)
and the coefficients $`F_{\alpha \beta }`$ are arbitrary real numbers. Here the constants $`𝒰`$’s are defined by

$`𝒰_1`$ $`=`$ $`{\displaystyle \frac{1}{16d^{}}}\left[\left(13\left({\displaystyle \frac{\eta ^{}}{\stackrel{~}{n}}}\right)_0\right)\widehat{n}{\displaystyle \frac{3}{2}}\widehat{T}\right],`$ (A.37)
$`𝒰_2`$ $`=`$ $`{\displaystyle \frac{1}{8d^{}}}\left[\left(1\left({\displaystyle \frac{\eta ^{}}{\stackrel{~}{n}}}\right)_0\right)\widehat{n}{\displaystyle \frac{1}{2}}\widehat{T}\right],`$ (A.38)
$`𝒰_3`$ $`=`$ $`{\displaystyle \frac{d}{4d^{}}}\left[\left(1+\left({\displaystyle \frac{\eta ^{}}{\stackrel{~}{n}}}\right)_0\right)\widehat{n}+{\displaystyle \frac{1}{2}}\widehat{T}\right]+{\displaystyle \frac{1}{4d^{}}}\left(\left({\displaystyle \frac{\eta ^{}}{\stackrel{~}{n}}}\right)_0\widehat{n}+{\displaystyle \frac{1}{2}}\widehat{T}\right)A_0^2.`$ (A.39)
It should be noted that all of the $`𝒰`$’s are proportional to the small parameter $`(c_v/\overline{c}_p)(P_r/d)`$. The results (56) and (A.35) show that $`\stackrel{~}{𝐮}^{{}_{}{}^{}(1)}`$ has indefinite constants that can not be fixed in the first approximation. We expect that they are fixed by the solubility conditions of the inhomogeneous linear equations in the higher order approximations. As shown in Eqs. (62) and (63), however, the expressions for macroscopic quantities such as the global temperature and the energy per particle up to the first approximation do not contain these indefinite constants.

## Appendix F Second order correction $`\stackrel{~}{T}^{(2)}`$

The purpose of this appendix is to illustrate how excitation of modes with larger wave numbers can be described in the higher order approximations.

The spatial average over the directions $`\beta =d^{}+1,\mathrm{},d`$ of both sides of the $`O(ϵ^1)`$ equation obtained from Eq. (36) yields

$`{\displaystyle \frac{\overline{c}_p}{c_v}}D_T^{}_{}{}^{}2\stackrel{~}{T}^{{}_{}{}^{}(2)}={\displaystyle \frac{\overline{c}_p}{c_v}}D_T^{}\left[\left(\left({\displaystyle \frac{\kappa ^{}}{\stackrel{~}{n}}}\right)_0\stackrel{~}{n}^{{}_{}{}^{}(1)}+{\displaystyle \frac{1}{2}}\stackrel{~}{T}^{(1)}\right)^{}\stackrel{~}{T}^{(1)}\right]`$ (A.40)
$``$ $`{\displaystyle \frac{4}{d}}D_{}\left[\left({\displaystyle \frac{\eta ^{}}{\stackrel{~}{n}}}\right)_0\stackrel{~}{n}^{{}_{}{}^{}(1)}+{\displaystyle \frac{1}{2}}\stackrel{~}{T}^{(1)}\right](^{}\stackrel{~}{𝐮}^{(0)})^{}:^{}\stackrel{~}{𝐮}^{(0)}{\displaystyle \frac{8}{d}}D_{}(^{}\stackrel{~}{𝐮}^{(0)})^{}:^{}\stackrel{~}{𝐮}^{{}_{}{}^{}(1)}`$
$`+`$ $`2\gamma _0\left[\left({\displaystyle \frac{\omega ^{}\stackrel{~}{n}}{\stackrel{~}{n}}}\right)_0\stackrel{~}{n}^{{}_{}{}^{}(1)}+{\displaystyle \frac{3}{2}}\stackrel{~}{T}^{(1)}\right]+\left(_\tau \mathrm{ln}\overline{T}\right)^{(0)}\left(\stackrel{~}{n}^{{}_{}{}^{}(1)}+\stackrel{~}{T}^{(1)}\right)+\left(_\tau \mathrm{ln}\overline{T}\right)^{(1)},`$
where $`\stackrel{~}{T}^{{}_{}{}^{}(2)}`$ is that second order correction to the scaled temperature which is spatially averaged over the directions $`\beta =d^{}+1,\mathrm{},d`$. This is again the inhomogeneous linear equation (the Laplace equation) for $`\stackrel{~}{T}^{{}_{}{}^{}(2)}`$ and it can be shown by substituting the result of the first approximation that the solubility conditions of this equation are satisfied. Because the general expression for $`\stackrel{~}{T}^{{}_{}{}^{}(2)}`$ is rather lengthy and complicated, we concentrate here only on the case that the zeroth approximation $`\stackrel{~}{𝐮}^{(0)}`$ is inhomogeneous only in one direction, say $`\alpha =1`$, i.e. $`d^{}=1`$; this is always true in two dimensions if $`k_0<k_{}^{}`$. Then $`\stackrel{~}{u}_1^{{}_{}{}^{}(1)}=0`$ and for $`\beta =2,3,\mathrm{},d`$,

$`\stackrel{~}{u}_\beta ^{{}_{}{}^{}(1)}(r_1)=𝒰_1A_{1\beta }\mathrm{cos}[3(k_0r_1+\theta _1^{})]+E_{1\beta }\mathrm{cos}(k_0r_1+\theta _1^{})+F_{1\beta }\mathrm{sin}(k_0r_1+\theta _1^{}),`$ (A.41)
where the coefficients $`E_{1\beta }`$ satisfy $`_{\beta =2}^dA_{1\beta }E_{1\beta }=𝒰_3`$. Substituting the expressions (58), (59), (60) and (A.41) into Eq. (A.40), and using the condition (33), we eventually get

$`\stackrel{~}{T}^{{}_{}{}^{}(2)}(r_1)`$ $`=`$ $`𝒯_1\mathrm{cos}[4(k_0r_1+\theta _1^{})]+𝒯_2\mathrm{cos}[2(k_0r_1+\theta _1^{})]`$ (A.42)
$`+`$ $`𝒯_3\mathrm{sin}[2(k_0r_1+\theta _1^{})]+𝒯_4,`$
where

$`𝒯_1`$ $`=`$ $`{\displaystyle \frac{\widehat{T}}{32}}\left[\left(38\left({\displaystyle \frac{\kappa ^{}}{\stackrel{~}{n}}}\right)_05\left({\displaystyle \frac{\eta ^{}}{\stackrel{~}{n}}}\right)_0\right)\widehat{n}{\displaystyle \frac{13}{2}}\widehat{T}\right],`$ (A.43)
$`𝒯_2`$ $`=`$ $`{\displaystyle \frac{\widehat{T}}{8}}\{(5+5\left({\displaystyle \frac{\eta ^{}}{\stackrel{~}{n}}}\right)_0+8\left({\displaystyle \frac{\omega ^{}}{\stackrel{~}{n}}}\right)_0)\widehat{n}+{\displaystyle \frac{29}{2}}\widehat{T}`$ (A.44)
$`+`$ $`4dA_0^2[(1+\left({\displaystyle \frac{\eta ^{}}{\stackrel{~}{n}}}\right)_0+2\left({\displaystyle \frac{\omega ^{}}{\stackrel{~}{n}}}\right)_0)\widehat{n}+{\displaystyle \frac{3}{2}}\widehat{T}]\},`$
$`𝒯_3`$ $`=`$ $`2\widehat{T}A_0^2{\displaystyle \underset{\beta =2}{\overset{d}{}}}A_{1\beta }F_{1\beta },`$ (A.45)
$`𝒯_4`$ $`=`$ $`{\displaystyle \frac{\widehat{T}\widehat{n}}{2}}.`$ (A.46)
Here $`𝒯_4`$ has been obtained from Eq. (33). It should be noted that $`𝒯_1`$, $`𝒯_2`$ and $`𝒯_4`$ are proportional to the square of the small parameter $`(c_v/\overline{c}_p)(P_r/d)`$; $`𝒯_3`$ contains the indefinite constants $`F_{1\beta }`$, which are expected to be fixed in the higher order approximations. Generalization of this result to the case of an arbitrary $`\stackrel{~}{𝐮}^{(0)}`$ given by Eq. (58) is straightforward.

In order to make analytical calculation executable in the even higher order approximations, it would be necessary to introduce a physically sensible assumption that the pressure is homogeneous in the asymptotic state, which is in accordance with the result (54). Consequences of this assumption, however, shall not be pursued in this paper.