format: text

text:
# Nonlinear Couette flow in a low density granular gas*footnote **footnote *Suggested running title: Couette flow in a granular gas

## I Introduction

An idealized representation of rapid flow granular media is given by the Boltzmann equation for a gas of smooth, hard, inelastic spheres.  This kinetic equation provides a basis for the study of a wide range of transport properties, and significant progress has been made in recent years using numerical methods such as direct Monte Carlo simulation.  In contrast, analytic studies have been limited to homogeneous states or to those with small spatial gradients.  Even in the homogenous case there is no exact solution corresponding to the Maxwellian distribution for elastic collisions. Recently, the method of kinetic models has been proposed for practical access to transport in more complex states.  This method has proven successful for elastic collisions, where several exact solutions to the model kinetic equations have been obtained for spatially inhomogeneous states far from equilibrium.  Such methods are potentially more important for granular gases since the states of interest are typically driven by external boundary conditions, posing intractable difficulties for solution to the Boltzmann equation. One illustration is the exact solution of a kinetic model for uniform shear flow, where the rheological properties were calculated and shown to be in good agreement with those obtained from simulation of the Boltzmann equation at low density  and the Enskog equation at finite densities.  Uniform shear flow is perhaps the most extensively studied inhomogeneous state for inelastic particles.  On the other hand, a more realistic shearing state is the planar Couette flow,  where temperature and density gradients coexist with the velocity field. The objective here is to provide an exact solution to a model kinetic equation corresponding to Couette flow with arbitrarily large temperature and flow field gradients, and arbitrary inelasticity. The results extend a previous exact analysis of momentum and heat transport far from equilibrium for a gas with elastic collisions.  Comparison with simulation of the Boltzmann equation for inelastic collisions again shows good agreement, confirming that the kinetic model is not only instructive but practical as well.

As noted above, exact solutions in kinetic theory for spatially inhomogeneous states are exceedingly rare. When, furthermore, such a solution corresponds to a hydrodynamic state far from equilibrium a unique and important benchmark is obtained for both conceptual and computational issues. In the case of granular gases there are two interesting examples: uniform shear flow and Couette flow. Both allow controlled discussion of nonlinear rheological properties that are important for a wide class of real granular flows. Uniform shear flow is a useful idealization, but the combined heat and momentum transport of Couette flow considered here is more realistic and a stronger test of hydrodynamic transport. For uniform shear, the condition of stationarity imposes a relationship between the shear rate and the coefficient of restitution, at fixed temperature. However, for Couette flow the shear rate and the coefficient of restitution are independent variables since the temperature profile is allowed to change. The parameter space for testing hydrodynamics is now two dimensional, with uniform shear flow recovered as a limiting line corresponding to zero curvature for the temperature profile. This is discussed further below.

In light of some earlier speculation that a hydrodynamic description for granular flows might be limited to weak dissipation and/or weakly inhomogeneous states, the results of this paper provide an example to the contrary: a) a hydrodynamic description applies, since all space dependence of the heat and momentum fluxes occurs via explicit functionals of the hydrodynamic fields and the macroscopic balance equations become a closed set of hydrodynamic equations which determine the non-trivial space dependence of the hydrodynamic fields; b) the hydrodynamic description applies even at strong dissipation and strong inhomogeneity (i.e., beyond the Navier-Stokes limit). The exact heat and momentum fluxes are characterized by generalizations of Fourier’s law and Newton’s viscosity where the thermal conductivity and viscosity are functions of the shear rate. The viscometric functions (normal stresses) are spatially constant, but non-trivial functions of the shear rate as well.

These strong results are obtained in the context of a simplified model of the Boltzmann equation. The mathematical and physical basis for this model as a good representation of the Boltzmann equation is discussed in Ref. . However, application of the model far from equilibrium as done here raises the question of its limitations, both quantitative and qualitative. The Monte Carlo simulations of the Boltzmann equation provide a more secure basis for the study of this macroscopic state. The good agreement obtained in this paper for this complex hydrodynamic state extends that previously demonstrated for uniform shear flow, providing additional support for future applications of the kinetic model to address other realistic boundary driven problems (e.g., vibrated columns or vertical flow through a chute in a gravitational field). Additional concerns are associated with the stability of the stationary state obtained here, particularly for large spatial gradients and strong dissipation. No analysis of this problem is provided here (no indication of instability is seen in the Monte Carlo simulation, although the imposed planar symmetry may suppress some possible mechanisms).

This paper is organized as follows. In the next Section, the state of macroscopic Couette flow is first considered at the level of the balance equations for mass, energy, and momentum. Based on previous analysis of this problem for elastic collisions, it is postulated that the balance equations support a solution with constant pressure and constant (dimensionless) shear rate. For consistency, it is shown that the heat and momentum fluxes must be given by generalizations of Fourier’s heat law and Newton’s viscosity law, as noted above. The temperature and velocity profiles are then determined in terms of the coefficients in these fluxes. To confirm that this self-consistent ansatz for the fluxes is correct, and to determine the explicit forms for the transport coefficients, an exact solution to the kinetic equation for this macroscopic state is constructed in Section III. Only the main results are quoted in this Section, while the detailed analysis is provided in the Appendices. The primary results of this analysis are summarized as follows: The heat flux is given by a generalized version of Fourier’s law which is linear in the temperature gradient without any restriction on this gradient being small. The proportionality constant is a nonlinear function of the shear rate, including anisotropy effects inducing a heat flux in the direction of flow as well as in the direction of the temperature gradient. The momentum flux components are characterized by three scalar functions of the shear rate, a shear viscosity and two viscometric functions. All of these properties are applicable even far from the reference homogeneous state and for all values of the shear rate. The only restriction is to conditions such that viscous heating dominates collisional cooling so that the temperature in the bulk is higher than that at the walls. These exact results for the kinetic model are compared with direct Monte Carlo simulation results for the Boltzmann equation in Section IV. The agreement is found to be good. Finally, the results are summarized and discussed in Section V.

## II Hydrodynamic description

Consider a low density granular gas of smooth hard spheres ($`d=3`$) or disks ($`d=2`$) of diameter $`\sigma `$ and mass $`m`$. Collisions between particles are characterized through a constant coefficient of normal restitution $`\alpha `$ with values $`0<\alpha 1`$, the largest value corresponding to the elastic limit. The macroscopic balance equations for mass, energy, and momentum are

$$D_tn+n𝐮=0,$$
(1)
$$D_tT+\frac{2}{dn}\left(P_{ij}_ju_i+𝐪\right)=\zeta T,$$
(2)
$$D_tu_i+(mn)^1_jP_{ij}=0,$$
(3)
where $`n`$ is the density, $`T`$ is the granular temperature, $`𝐮`$ is the flow velocity, and $`D_t=_t+𝐮`$ is the material derivative. In addition, $`\zeta `$ is the cooling rate (related to the collisional energy dissipation), $`𝖯`$ is the pressure tensor (related to the transport of momentum), and $`𝐪`$ is the heat flux (related to the transport of energy). Whenever these fluxes can be expressed in terms of the hydrodynamic fields, Eqs. (1)–(3) become a closed set of hydrodynamic equations for these fields.

The specific problem considered in this paper is steady Couette flow \[cf. Fig. 1\]. The gas of inelastic hard spheres is enclosed between two parallel plates at $`y=\pm L`$ in relative motion along the $`x`$-axis and maintained, in general, at different temperatures. The resulting flow velocity is along the $`x`$ axis and, from symmetry, it is expected that the hydrodynamic fields vary only in the $`y`$ direction. The pressure tensor and the heat flux characterize a macroscopic state with combined heat and momentum transport. The objective here is to determine the hydrodynamic fields. Under the above conditions, Eq. (1) is identically satisfied, while the balance equations (2) and (3) yield

$$\frac{2}{dn}\left(P_{xy}\frac{u_x}{y}+\frac{q_y}{y}\right)=\zeta T,$$
(4)
$$\frac{P_{xy}}{y}=0,$$
(5)
$$\frac{P_{yy}}{y}=0.$$
(6)
Equations (4)–(6) are still exact. The temperature $`T`$ is defined in the usual way in terms of the average kinetic energy, and is related to the pressure $`p`$ through the ideal gas law

$$p=nT=\frac{1}{d}P_{ii}.$$
(7)
The second equality follows from the definition of the pressure tensor (see next Section). There are three independent hydrodynamic fields, which are taken here to be the pressure, temperature, and the $`x`$ component of the flow velocity (the other $`d1`$ components vanish). The boundary conditions impose a global shear rate given by $`\left[u_x\left(L\right)u_x\left(L\right)\right]/2L`$. It is useful to introduce a dimensionless local shear rate $`a`$

$$a=\frac{1}{\nu _0}\frac{u_x}{y}.$$
(8)
Here $`\nu _0`$ is a convenient collision frequency $`\nu _0=p/\eta _0`$,  where $`\eta _0`$ is the Navier-Stokes shear viscosity in the elastic limit. The variable $`a`$ will be one of the control parameters for the Couette flow state. It has the simple physical interpretation as the ratio of the mean free path $`\mathrm{}=v_0/\nu _0`$ (where $`v_0=\sqrt{2T/m}`$ is the thermal velocity) to the relevant hydrodynamic length $`h=v_0(u_x/y)^1`$. We now look for special solutions to the macroscopic balance equations characterized by constant pressure and constant $`a`$,

$$p=\text{constant},a=\text{constant}.$$
(9)
To provide a closed set of equations for the hydrodynamic fields, the heat and momentum fluxes must be expressed in terms of these fields. Forms for these fluxes are postulated in this Section, as a generalization of those obtained from a similar analysis of Couette flow for a fluid with elastic collisions. Their verification is given in the following Section and in Appendix C. Momentum transport is typically characterized by three scalar rheological functions, the shear viscosity and two viscometric functions. The $`xy`$ element for shear stresses is represented in the form of a generalized Newton’s viscosity law

$$P_{xy}=F_\eta (a,\alpha )\eta _0\frac{u_x}{y}.$$
(10)
The function $`F_\eta (a,\alpha )\eta _0`$ is the generalized shear viscosity, where the dissipation and all nonlinear rheological effects are included through the factor $`F_\eta (a,\alpha )`$. The normal stress components are specified in terms of the dimensionless viscometric functions

$$\mathrm{\Psi }_1(a,\alpha )=\frac{P_{yy}P_{xx}}{pa^2},\mathrm{\Psi }_2(a,\alpha )=\frac{P_{zz}P_{yy}}{pa^2}.$$
(11)
The heat flux is represented as

$$q_i=\lambda _{ij}_jT.$$
(12)
The tensor $`\lambda _{ij}`$ depends on two independent scalar “transport coefficients”

$$\lambda _{ij}=\lambda _0\left[F_\lambda (a,\alpha )\delta _{ij}+\mathrm{\Phi }(a,\alpha )a_{ij}\right],$$
(13)
where $`a_{ij}=a\delta _{ix}\delta _{jy}`$. The first term of Eq. (13) yields a generalization of Fourier’s law with the thermal conductivity $`\lambda _0F_\lambda (a,\alpha )`$ modified by the nonlinear rheological factor $`F_\lambda (a,\alpha )`$. The second term of (13) provides information on the anisotropy created by the Couette flow since it gives a heat flux along the $`x`$-axis due to a thermal gradient parallel to the $`y`$-axis. This effect is also nonlinear with no analogue at Navier-Stokes order, being first order in both the shear rate and the temperature gradient. In fact, Eqs. (12) and (13) are a generalization of the heat flux constitutive equation to Burnett order  for the geometry of the problem in the elastic case, according to which $`\mathrm{\Phi }(a=0,\alpha =1)=(\theta _5\theta _4\mathrm{ln}\eta _0/\mathrm{ln}T)\eta _0/2m\lambda _0=3.5\text{ (Boltzmann equation) or }2.8\text{ (BGK model)}`$.

Since $`p`$ and $`a`$ are constant, Eqs. (10) and (11) are consistent with Eqs. (5) and (6). Dimensional analysis shows that

$$\eta _0T^{1/2},\lambda _0T^{1/2},\zeta T^{1/2},$$
(14)
where it is understood that the independent variables are $`p`$ and $`T`$. Then Eq. (4) becomes

$$T^{1/2}\frac{}{y}T^{1/2}\frac{T}{y}=\frac{2mTp^2}{\eta _0^2}\text{Pr}\gamma ,$$
(15)
where $`\gamma (a,\alpha )`$ is the dimensionless constant defined by

$$\text{Pr}\gamma =\frac{\eta _0}{2m\lambda _0F_\lambda }\left(a^2F_\eta \frac{d}{2p}\eta _0\zeta \right)$$
(16)
and $`\text{Pr}=(d+2)\eta _0/2m\lambda _0`$ is the Prandtl number. Equation (15) is a closed equation that determines $`T=T(a,\alpha ,y)`$. With this known, the velocity field is determined from the condition $`a=`$ constant

$$\frac{u_x}{y}=a\frac{p}{\eta _0(T)}.$$
(17)
Before giving the solution, it is instructive to express the temperature in terms of the flow velocity: $`T(y)T(u_x)`$. Then Eqs. (15) and (17) give directly

$$\frac{^2T}{u_x^2}=\frac{2m}{a^2}\text{Pr}\gamma .$$
(18)
The sign of the constant $`\gamma `$ is a result of the competition between viscous heating and inelastic cooling (represented by the first and second terms, respectively, on the right-hand side of (16)). If the dissipation is sufficiently low and/or the shear rate is sufficiently large so that $`\zeta <(2pF_\eta /d\eta _0)a^2`$, then $`\gamma >0`$ and the temperature profile is convex. Equation (18) implies that the temperature is simply a quadratic function of the flow velocity,

$$T=T_0\left[1c^2(u_xu_0)^2\right],c^2\frac{m}{T_0a^2}\text{Pr}\gamma ,$$
(19)
where $`u_0`$ is an arbitrary constant and $`T_0`$ is the temperature at the point where $`T/u_x=0`$. Here and below attention is restricted to the case $`\gamma >0`$ which implies that the shear rate must be larger than a certain critical value, $`a>a_c(\alpha )`$. The relationship of $`T`$ and $`u_x`$ can be viewed as a nonequilibrium “equation of state.” The explicit $`y`$-dependence of $`u_x`$ and $`T`$ can be easily obtained from Eqs. (19) and (17), where the constants are determined by the boundary conditions. For instance, suppose that $`y=y_0`$ is the point where the temperature presents an extremum and that, without loss of generality, $`u_x(y_0)u_0=0`$. Each particular situation is then characterized by the uniform pressure $`p`$, the local temperature $`T(y_0)=T_0`$ and the shear rate $`a`$. Equation (17) becomes

$$\left(1c^2u_x^2\right)^{1/2}\frac{u_x}{y}=a\frac{p}{\eta _0(T_0)},$$
(20)
whose (implicit) solution is

$$y=y_0+\frac{\eta _0(T_0)}{2cpa}\left[\mathrm{sin}^1\left(cu_x\right)+cu_x\left(1c^2u_x^2\right)^{1/2}\right].$$
(21)
Note that Eqs. (19) and (21) are only meaningful for $`u_x^2c^2`$, which implies $`|yy_0|\pi \eta _0(T_0)/4cpa`$. The velocity profile (21) and the temperature profile (19) are shown in Fig. 1 for a representative case ($`\alpha =0.8`$, $`a=0.95`$) with $`y_0=0`$. The actual shear rate $`u_x/y`$ is practically constant across most of the system, but it rapidly increases as the temperature becomes much smaller than $`T_0`$ \[cf. Eq. (17)\].

To summarize, the conditions of constant pressure and shear rate, Eq. (9), together with the constitutive relations for the heat and momentum fluxes, Eqs. (10)–(13) allow an exact stationary solution to the macroscopic balance equations. In the next Section the assumed conditions and constitutive relations are confirmed by an exact solution to the kinetic equation. In addition, this solution provides explicit expressions for the transport coefficients $`F_\eta `$, $`\mathrm{\Psi }_{1,2}`$, $`F_\lambda `$, and $`\mathrm{\Phi }`$ as functions of $`a`$ and $`\alpha `$.

## III Kinetic equation and exact solution

Consider a low density granular gas of smooth hard spheres in $`d`$ dimensions. At sufficiently low density the phase space distribution function, $`f(𝐫,𝐯,t)`$, is determined from the Boltzmann kinetic equation modified appropriately for inelastic collisions. Due to the mathematical complexity of this equation, analysis has been limited to perturbative approximations for small spatial inhomogeneities.  To describe more general nonequilibrium states it is useful to replace the Boltzmann equation with a more tractable model kinetic equation constructed to preserve its most important qualitative features (e.g., macroscopic balance equations). The model kinetic equation chosen for analysis here is

$$_tf+𝐯f=\nu (ff_{\mathrm{}})+\frac{1}{2}\zeta _𝐯\left[\left(𝐯𝐮\right)f\right],$$
(22)
where $`f_{\mathrm{}}`$ is the local Maxwellian distribution

$$f_{\mathrm{}}(𝐫,𝐯,t)=n(𝐫,t)\left[\frac{m}{2\pi T(𝐫,t)}\right]^{d/2}\mathrm{exp}\left[\frac{m\left(𝐯𝐮(𝐫,t)\right)^2}{2T(𝐫,t)}\right].$$
(23)
This distribution function is parameterized by the nonequilibrium density $`n`$ , granular temperature $`T`$, and flow velocity $`𝐮`$, which are defined in terms of moments of $`f`$:

$$n=𝑑𝐯f,T=\frac{m}{dn}𝑑𝐯V^2f,𝐮=\frac{1}{n}𝑑𝐯𝐯f,$$
(24)
where $`𝐕=𝐯𝐮`$ is the velocity relative to the local flow. The right side of (22) is a model for the nonlinear Boltzmann collision operator. The first term describes collisional relaxation towards the local Maxwellian, with collision rate $`\nu `$. The second term describes the dominant collisional cooling effects, where $`\zeta `$ is the cooling rate. The necessity for this term to accurately represent the spectrum of the Boltzmann collision operator is discussed in Ref. . However, it can be viewed more simply as an effective “drag” force that produces the same energy loss rate as that produced by the inelastic collisions. The parameters $`\nu `$ and $`\zeta `$ are chosen for good quantitative agreement of the viscosity and cooling rate with those obtained from the Boltzmann equation,

$$\zeta =\frac{2\pi ^{(d1)/2}}{d\mathrm{\Gamma }(d/2)}\sigma ^{d1}n\left(\frac{T}{m}\right)^{1/2}(1\alpha ^2),\nu =\nu _0\zeta ,$$
(25)
$$\nu _0kn\sigma ^{d1}\left(\frac{\pi T}{m}\right)^{1/2},$$
(26)
where $`k=15/16`$ for $`d=3`$ (spheres) and $`k=2`$ for $`d=2`$ (disks). Both $`\nu `$ and $`\zeta `$ depend on the density and temperature, whose space and time dependence has been left implicit. Further details motivating these choices can be found in Ref. , where a somewhat more sophisticated model is described. Equations (22)–(25) define the kinetic equation to be applied in this work.

By taking moments with respect to $`1`$, $`𝐯`$, and $`v^2`$, this model kinetic equation yields the same form of the macroscopic balance equations for mass, energy, and momentum, Eqs. (1)–(3), as those given from the Boltzmann equation, thus confirming the interpretation of $`\zeta `$ as the cooling rate. The pressure tensor $`𝖯`$ and the heat flux $`𝐪`$ are given by

$$𝖯=m𝑑𝐯\mathrm{𝐕𝐕}f\text{ },𝐪=\frac{m}{2}𝑑𝐯V^2𝐕f.$$
(27)
We now consider the specific problem of steady Couette flow described in the previous Section. The main objective is to verify the constitutive equations (10)–(13) from the fundamental definitions (27), and to verify the assumed constancy for the pressure and shear rate. For the chosen geometry the kinetic equation becomes

$$\left(1\zeta ^{}+v_y\nu _0^1_y\right)f\frac{1}{2}\zeta ^{}_𝐯\left(𝐕f\right)=(1\zeta ^{})f_{\mathrm{}},$$
(28)
where

$$\zeta ^{}=\frac{\zeta }{\nu _0}=\frac{2\pi ^{(d2)/2}}{kd\mathrm{\Gamma }(d/2)}(1\alpha ^2)$$
(29)
is a constant that henceforth will characterize the inelasticity dependence. Equation (28) allows one to understand the structural simplification afforded by the kinetic model relative to the Boltzmann equation. Assuming the hydrodynamic fields of Section II, the local equilibrium distribution function $`f_{\mathrm{}}`$ is completely specified and Eq. (28) becomes a linear inhomogeneous partial first order differential equation that is readily solved with specified boundary conditions. This determines explicitly the entire velocity dependence. However, the solution is only formal since the assumed fields must satisfy (24). These are consistency conditions that are needed to justify the assumed forms for the hydrodynamic fields. Alternatively, since these fields have been shown to follow from the constitutive equations (10)–(13) for the heat and momentum fluxes in Section II, it is sufficient to show that these equations are verified by use of this formal solution in (27). Both the consistency conditions and the fluxes are verified explicitly in Appendix C.

The general solution to Eq. (28) is

$$f(y,𝐯)=f_B(y,𝐯)+(1\zeta ^{})_0^{\mathrm{}}𝑑te^{(1\zeta ^{}\frac{d}{2}\zeta ^{})t}e^{𝒟t}f_{\mathrm{}}(y,𝐯),$$
(30)
where $`𝒟`$ is the sum of a generator for translations of $`y`$ and the generator for scale transformation of $`𝐕`$,

$$𝒟v_y\nu _0^1_y\frac{\zeta ^{}}{2}𝐕_𝐯.$$
(31)
The explicit effects of this generator are described in Appendix A. The first term of (30), $`f_B(y,𝐯)`$, is a solution to the homogeneous kinetic equation obtained from (28) by setting $`f_{\mathrm{}}0`$. The detailed form of this contribution is determined by the chosen boundary conditions. The physical boundary conditions are specification of the half distributions for velocities directed away from the walls at $`y=\pm L`$, given explicitly or implicitly in terms of the distributions for velocities directed at the walls. The specific relationship characterizes the motion and temperature of the walls. One possibility for Couette flow is diffuse conditions, where the distributions for velocities away from the walls are given by a Maxwellian whose parameters specify the temperatures $`T_\pm `$ and velocities $`U_\pm `$ of the walls at $`y=\pm L`$.

In general $`f_B(y,𝐯)`$ has a detailed explicit dependence on the geometry of the system and is responsible for the “boundary layer” near the wall. In contrast, the second term of (30) is an example of a “normal” solution where all of its dependence on $`y`$ occurs only through the hydrodynamic fields. Here we consider idealized boundary conditions for which $`f_B(y,𝐯)0`$ so that the entire solution is given by the second term of (30). This idealized boundary condition eliminates boundary layers and admits the possibility of simple hydrodynamic profiles that are exact throughout the system. Furthermore, the solution to the kinetic equation is now normal so that the fluxes become functionals of the hydrodynamic fields and (1)–(3) become a closed set of hydrodynamic equations. The idealized boundary conditions correspond to the limit $`T_\pm 0`$ (as illustrated in Fig. 1). This is verified in Appendix B and will not be discussed further here.

Equation (30) with these idealized boundary conditions admits an exact solution for steady Couette flow \[see Eq. (C2) for its explicit form\] characterized by the solutions to Eqs. (8), (9), and (15) given in Section II, which are rewritten here for the sake of completeness:

$$p(y)=n(y)T(y)=\text{const},\nu _{0}^{}{}_{}{}^{1}_yu_x=a,(\nu _{0}^{}{}_{}{}^{1}_y)^2T=2m\text{Pr}\gamma .$$
(32)
The Prandtl number, defined following (16), is equal to 1 in the kinetic model while $`\text{Pr}=(d1)/d`$ for the Boltzmann equation. The fact that $`\text{Pr}=1`$ in the kinetic model is a consequence of the introduction of a single collision frequency, that is unable to reproduce simultaneously the exact shear viscosity and thermal conductivity Navier-Stokes coefficients. The components of the flow velocity in the $`y`$ and $`z`$ directions vanish. The parameter $`a`$ is the dimensionless shear rate and characterizes the relative velocities of the wall. Together with the coefficient of restitution $`\alpha `$ it is a given control parameter in terms of which all transport properties are expressed. The dimensionless parameter $`\gamma (a,\alpha )`$ characterizes the curvature of the temperature field and is a consequence of both the viscous heating due to the shear rate and the collisional dissipation characterized by $`\alpha `$. The consistency conditions and forms for the fluxes are verified in Appendix C, confirming that the distribution function (30) is an exact solution to the kinetic equation. This verification also fixes the functional dependence of $`\gamma (a,\alpha )`$ on the control parameters. The details of the analysis also are given in Appendix C, yielding the following implicit equation

$$1=\frac{4(1\zeta ^{})}{d\pi ^{1/2}}_0^{\mathrm{}}𝑑ue^{u^2}_0^{\mathrm{}}𝑑te^t\left[\frac{d1}{2}+\left(1+a^2t^2\right)u^2\right]\left(1+2\gamma w^2\right)^1.$$
(33)
where

$$w(u,t)=\frac{2u}{\zeta ^{}}\left(1e^{\frac{1}{2}\zeta ^{}t}\right).$$
(34)
Interestingly, the representation (33) exists only for $`\gamma 0`$ or, equivalently, for $`a`$ equal to or larger than some critical value $`a_c`$ of the shear rate. This critical shear rate corresponds to $`\gamma =0`$, or uniform temperature (see Eq. (19)). Setting $`\gamma =0`$ in Eq. (33), one gets

$$a_c^2=\frac{d}{2}\frac{\zeta ^{}}{1\zeta ^{}}.$$
(35)
For the limiting case $`a=a_c`$ the viscous heating is exactly balanced by collisional cooling and the gas is in a state of uniform shear flow. This is further discussed at the end of this Section. In the regime of low shear rates and low dissipation, Eq. (33) yields $`\gamma (a^2d\zeta ^{}/2)/(d+2)`$ and consequently $`a_c^2d\zeta ^{}/2`$.

The explicit expressions for the momentum and heat fluxes are also derived in Appendix C. From them it is possible to identify the generalized transport coefficients defined in Eqs. (10)–(13). The viscosity coefficient is

$$F_\eta (a,\alpha )=\frac{4(1\zeta ^{})}{\pi ^{1/2}}_0^{\mathrm{}}𝑑ue^{u^2}u^2_0^{\mathrm{}}𝑑te^tt\left(1+2\gamma w^2\right)^1.$$
(36)
The viscometric functions characterizing normal stresses are

$$\mathrm{\Psi }_1(a,\alpha )=\frac{2(1\zeta ^{})}{a^2\pi ^{1/2}}_0^{\mathrm{}}𝑑ue^{u^2}_0^{\mathrm{}}𝑑te^t\left[12\left(1a^2t^2\right)u^2\right]\left(1+2\gamma w^2\right)^1,$$
(37)
$$\mathrm{\Psi }_2(a,\alpha )=\frac{2(1\zeta ^{})}{a^2\pi ^{1/2}}_0^{\mathrm{}}𝑑ue^{u^2}\left(12u^2\right)_0^{\mathrm{}}𝑑te^t\left(1+2\gamma w^2\right)^1.$$
(38)
The thermal conductivity coefficient is given in terms of $`\gamma `$ and $`F_\eta `$ by Eq. (16)

$$F_\lambda =\frac{\eta _0}{2m\lambda _0\text{Pr}\gamma }\left(a^2F_\eta \frac{d}{2p}\eta _0\zeta \right).$$
(39)
Finally, the expression for the cross coefficient $`\mathrm{\Phi }`$ in the heat flux is

$`\mathrm{\Phi }(a,\alpha )`$ $`=`$ $`{\displaystyle \frac{8(1\zeta ^{})}{\left(d+2\right)\pi ^{1/2}}}{\displaystyle _0^{\mathrm{}}}𝑑ue^{u^2}u{\displaystyle _0^{\mathrm{}}}𝑑te^{\left(1+\frac{1}{2}\zeta ^{}\right)t}tw\left[{\displaystyle \frac{d+1}{2}}+\left(1+a^2t^2\right)u^2\right]`$ (41)
$`\times \left(1+2\gamma w^2\right)^2.`$
In summary, the distribution function, hydrodynamic fields, and transport coefficients for the heat and momentum fluxes have been determined exactly in terms of the imposed shear rate $`a`$ and the restitution coefficient $`\alpha `$. For small $`a`$ and $`\alpha `$ the results agree with predictions of the Navier-Stokes hydrodynamics. More generally, they extend the description of Couette flow to large spatial inhomogeneity and strong dissipation. The final results are still only implicit, but the entire problem has been brought to quadratures. A numerical evaluation of these expressions is provided in the next Section for comparison with numerical simulation of the Boltzmann kinetic equation.

It is easy to check that all results presented in this Section reduce to those previously derived in the elastic limit $`\alpha =1`$ by using the Bhatnagar-Gross-Krook (BGK) model.  In addition, the results of this Section include as a limiting case those corresponding to the uniform shear flow for an inelastic gas.  This happens when the shear rate and the inelasticity combine to yield a zero curvature for the temperature profile ($`\gamma 0`$). Equation (35) can be seen as the relationship between the reduced shear rate and the coefficient of restitution, which are not independent in the uniform shear flow problem. The transport coefficients in this limiting case are given in Appendix D. It is also interesting to note that this uniform shear flow state is mathematically equivalent to the corresponding case of elastic collisions with an external thermostat force $`\frac{1}{2}m\zeta 𝐕`$ adjusted to control the viscous heating.  In this context, Eq. (35) (with the adequate change of units) gives the shear rate dependence of the thermostat and this coincides with the previous result derived in the uniform shear flow problem for an elastic gas, as expected.

## IV Comparison with Monte Carlo simulations

In this Section we compare the predictions of the kinetic model with computer simulations of the Boltzmann equation by means of the Direct Simulation Monte Carlo (DSMC) method. Although this method was originally devised for elastic particles,  its extension to the inelastic case is straightforward.  This method has proven to be an efficient and reliable tool for solving numerically the Boltzmann equation.

We have used the adaptation of the DSMC method to the steady Couette flow for a granular gas following the same steps as those recently worked out in the elastic case,  so that the technical details of the simulations are omitted here. We have considered a system of hard spheres ($`d=3`$) with three different values for the coefficient of restitution: $`\alpha =1`$ (elastic case), $`\alpha =0.9`$, and $`\alpha =0.8`$. For each case, typically four or five different shearing states have been taken. Once a steady state is reached, the profiles of the hydrodynamic quantities and of the fluxes are measured in the bulk region. From these profiles the transport coefficients are identified.

The two parameters defining the kinetic model are the collision frequency $`\nu `$ and the cooling rate $`\zeta `$, given by (25). These choices were made to provide a quantitative representation of the shear viscosity at the Navier-Stokes level, while retaining the relative simplicity of the kinetic equation. To verify the extent to which these choices are valid for larger spatial gradients we compare the cooling rate of the kinetic model with that obtained from the Boltzmann equation. Both are proportional to $`(1\alpha ^2)`$. In Fig. 2 we plot the ratio $`\zeta ^{}/(1\alpha ^2)`$ obtained from the simulation versus the shear rate for different values of $`\alpha `$. It is apparent that $`\zeta ^{}`$ is well approximated by the kinetic model form, although the latter tends to underestimate the correct value, especially as the shear rate increases. This supports the expectation that the kinetic model retains quantitative as well as qualitative validity even for the extreme conditions admitted by the exact solution of the previous Sections.

The simulation results indicate that the pressure $`p`$ is practically constant, in agreement with the first assumption in (32). Also, the actual shear rate $`_yu_x`$ multiplied by $`T^{1/2}`$ is practically constant, allowing identification of a reduced shear rate $`a`$ in agreement with the second equality in (32). We have also checked that the temperature profile is well characterized by the third equation in (32) (with $`\text{Pr}=2/3`$). The parameter $`\gamma `$ characterizing the curvature of the temperature profile is plotted in Fig. 3 as a function of the shear rate for $`\alpha =0.8`$, $`0.9`$, and $`1`$. The agreement between theory and simulation is quite good. Surprisingly, this is true even for strong inelasticity. There is a remarkable influence of the degree of dissipation on $`\gamma `$ at a given value of the shear rate. As the coefficient of restitution $`\alpha `$ decreases, so does the temperature variation across the system, so that the state tends to resemble the uniform shear flow, in which case the kinetic model is known to be quite accurate.  Figure 3 also confirms that, for a given value of $`\alpha `$, there exists a critical shear rate $`a_c`$ such that $`\gamma 0^+`$ as $`a`$ tends to $`a_c`$ from above. For shear rates smaller than the critical value $`a_c`$ the collisional cooling effect dominates over the viscous heating and the temperature presents a minimum instead of a maximum, which could be characterized by a negative $`\gamma `$.  However, the solution obtained here does not extend to negative $`\gamma `$, as noted above. The dependence of $`a_c`$ on the coefficient of restitution $`\alpha `$ is shown in Fig. 4. The simple prediction (35) reproduces quite well the simulation data. Figure 4 can also be interpreted as the two dimensional parameter space of the inelastic Couette flow. The elastic Couette problem corresponds to the line $`\alpha =1`$ and arbitrary $`a`$. On the other hand, the uniform shear flow for a dilute granular gas is represented by the line $`a=a_c(\alpha )`$ for arbitrary $`\alpha `$. The case studied in this paper applies to the region enclosed between both lines, namely, for arbitrary $`\alpha `$ and $`a`$, provided that $`aa_c(\alpha )`$.

Now we compare the simulation results and the theoretical predictions for the five generalized transport coefficients in the domain $`aa_c`$. For the sake of clarity, we have also plotted the curves representing the transport coefficients for the uniform shear flow when eliminating the coefficient of restitution in favor of the shear rate. In Fig. 5 we plot the shear rate dependence of the viscosity function $`F_\eta `$ defined by Eq. (10). Needless to say, this is the most relevant quantity in a shearing state. Regardless of the value of $`\alpha `$, shear thinning effects are present, i.e., $`F_\eta `$ decreases as the shear rate increases. This rheological behavior is quantitatively well described by the kinetic model, again especially for large dissipation. In general, however, the model tends to overestimate the value of $`F_\eta `$. The two viscometric functions (11) are plotted in Figs. 6 and 7. While the agreement between the theoretical predictions and the simulation results in the case of $`\mathrm{\Psi }_1`$ is similar to the one observed with $`F_\eta `$, the agreement in the case of $`\mathrm{\Psi }_2`$ is only qualitative. In particular, the model succeeds in capturing the non-monotonic behavior of $`\mathrm{\Psi }_2`$ in the inelastic case.

The transport of energy is characterized by the coefficients $`F_\lambda `$ and $`\mathrm{\Phi }`$, defined by Eq. (13). They are plotted in Figs. 8 and 9, respectively. As happens with $`F_\eta `$, the generalized thermal conductivity $`F_\lambda `$ decreases as the shear rate increases. The accuracy of the kinetic model predictions for this quantity is less satisfactory than in the case of the viscosity function. This is a consequence of the fact that the model slightly underestimates the cooling rate $`\zeta ^{}`$ (Fig. 2) and the curvature parameter $`\gamma `$ (Fig. 3), while it slightly overestimates $`F_\eta `$ (Fig. 5) and so, according to Eq. (39), all these discrepancies contribute to magnify the inaccuracy of the value of $`F_\lambda `$ predicted by the model as well. In the case of the coefficient $`\mathrm{\Phi }`$, the agreement between theory and simulation is quite good. This is rather satisfactory if one takes into account that this is a coefficient measuring complex coupling effects between the velocity and temperature gradients, which are absent at the Navier-Stokes regime. Also note that the kinetic model yields $`\mathrm{\Phi }(a=0,\alpha =1)=2.8`$, while $`\mathrm{\Phi }(a=0,\alpha =1)=3.5`$ in the Boltzmann equation. This is again due to the fact that relaxation effects are described by a single frequency in the BGK model.

## V Discussion

The exact macroscopic balance equations for mass, energy, and momentum become closed hydrodynamic equations when the energy and momentum fluxes are specified in terms of the hydrodynamic fields. In Section II specific forms for these fluxes were postulated to hold for Couette flow, and the resulting hydrodynamic equations were solved exactly for the temperature, pressure, and flow velocity. In Section III an exact solution to a model kinetic equation was obtained for Couette flow and the assumed forms for the fluxes was verified. In addition, it is shown in Appendix C that the hydrodynamic fields obtained directly from this distribution function are the same as those obtained from solution to the hydrodynamic equations. Thus, both hydrodynamic and kinetic descriptions have been verified as exact and self-consistent. The kinetic theory provides the distribution function in phase space, Eq. (C2), as an explicit functional of the hydrodynamic fields and therefore allows in addition the calculation of non-hydrodynamic properties in terms of them. These fields are given by Eqs. (19) and (21) with the parameter $`\gamma (a,\alpha )`$ determined by the consistency condition (33).

Some comments to summarize these results and provide perspective are as follows:

* The model kinetic equation considered is a generalization of the familiar BGK equation used to model the Boltzmann equation for elastic collisions, and reduces to it for $`\alpha =1`$. The utility of the BGK equation to address complex states not accessible via the Boltzmann equation is well-established. The results here and those of Refs.  and  confirm its value for granular media as well. Although these kinetic model equations are structurally simple, they are highly nonlinear due to the implicit dependence of $`f_{\mathrm{}}`$ on the distribution function via the hydrodynamic fields.
* The kinetic model equation provides a coupled set of singular nonlinear integral equations for the hydrodynamic fields. Incorporation of specific boundary conditions is straightforward and an advantage of this description. However, explicit results generally require approximations near the homogeneous state or numerical solution. Exact analytic results valid for large spatial gradients and arbitrary inelasticity are rare even for the kinetic model equations. They require idealized boundary conditions such that boundary layer complications are avoided and “simple” profiles are possible in the appropriate variables.
* The idealized boundary conditions considered here are homogeneous, such that the solution is “normal” — all space dependence occurs entirely through the hydrodynamic fields. Since the fluxes in macroscopic balance equations (1)–(3) are determined from moments of this distribution, they can be calculated as functions of the hydrodynamic fields. Use of these in the macroscopic balance equations then gives a closed set of equations for the fields. The analysis here provides a non trivial example of the existence of hydrodynamics for a strongly inhomogeneous state together with the solution to those equations. The solution is valid even for strong dissipation subject to the condition $`\gamma 0`$ or $`aa_c=\sqrt{d\zeta ^{}\left(\alpha \right)/2[1\zeta ^{}(\alpha )]}`$.
* The results derived here include as a limiting case those previously obtained for the uniform shear flow problem.  This corresponds to an exact balance between the viscous heating and the inelastic cooling, resulting in a uniform temperature ($`\gamma 0^+`$). From that point of view, the critical shear rate $`a_c`$ gives the relationship between the shear rate and the coefficient of restitution in the uniform shear flow.
* The existence of a critical shear rate $`a_c`$ below which the analysis does not apply restricts the conditions such that viscous heating dominates collisional cooling. The results do not apply, therefore, to the simple case of zero shear rate and constant temperature walls. In this latter case, simulations and asymptotic analysis  suggest that the temperature profile is very close to that of (32) but with $`\gamma 0`$. The proof for this is still not available.
* Although the kinetic model is only a crude representation of the Boltzmann equation, it does preserve the most important features for transport, such as the reference homogeneous state and the macroscopic conservation laws. Otherwise there are only two adjustable parameters, the collision rate $`\nu `$ and the cooling rate $`\zeta `$, to improve correspondence with the Boltzmann equation. Here they have been chosen to fit the Navier-Stokes order viscosity. Surprisingly, this is sufficient to allow good agreement with the Boltzmann results for rheological properties far from equilibrium. This is true for elastic collisions  and now has been confirmed here for inelastic collisions.

###### Acknowledgements.

This work has been done under the auspices of the Agencia Española de Cooperación Internacional (Programa de Cooperación Interuniversitaria Hispano-Marroquí). J. M. M., V. G., and A.S. acknowledge partial support from the DGES (Spain) through grant No. PB97-1501 and from the Junta de Extremadura (Fondo Social Europeo) through grant No. IPR99C031. A.S. is also grateful to the DGES (Spain) for a sabbatical grant (No. PR2000-0117 0027235927). The research of J. D. W. was supported in part by a grant from the National Science Foundation, NSF PHY 9722133.

## A Action of the exponential operator

In this Appendix we evaluate the action of the operator $`e^{𝒟t}`$, where $`𝒟`$ is given by Eq. (31). To do this, it is convenient to refer the velocities of the particles to a local Lagrangian frame $`𝐕=𝐯𝐮`$ and to introduce a scaled space variable as

$$s=_{y_0}^y𝑑y^{}\nu _0(y^{}).$$
(A1)
In terms of this variable, the velocity and temperature profiles become simply

$$u_x\left(s\right)=as,T\left(s\right)=T_0m\gamma s^2,$$
(A2)
where we have taken into account that $`\lambda _0=(d+2)\eta _0/2m`$ in the BGK kinetic model. In the new Lagrangian frame the operator $`𝒟`$ becomes

$$𝒟=V_y_saV_y_{V_x}\frac{\zeta ^{}}{2}𝐕_𝐕,$$
(A3)
where the derivative $`_s`$ is taken at constant $`𝐕`$.

Let us consider the operators $`AV_y_s`$, $`B\frac{1}{2}\zeta ^{}𝐕_𝐕`$, and $`CaV_y_{V_x}`$. It is easy to see that the operators $`A+B`$ and $`C`$ commute, namely,

$$[A+B,C](A+B)CC(A+B)=0.$$
(A4)
Thus, the action of the operator $`e^{𝒟t}e^{(A+B+C)t}`$ on a given function $`F(s,𝐕)`$ is

$`e^{(A+B+C)t}F(s,𝐕)`$ $`=`$ $`e^{(A+B)t}e^{Ct}F=e^{Ct}e^{(A+B)t}F`$ (A5)
$`=`$ $`e^{(A+B)t}F(s,𝐕+t𝖺𝐕),`$ (A6)
where $`𝖺`$ is the matrix with elements $`a_{ij}=a\delta _{ix}\delta _{jy}`$. It remains to evaluate the action of the operator $`e^{(A+B)t}`$ on an arbitrary function of the velocity, $`h(s,𝐕)`$:

$$e^{(A+B)t}h(s,𝐕)e^{Bt}H(t,s,𝐕).$$
(A7)
An equation for $`H(t,s,𝐕)`$ follows by differentiating both sides to get

$$_tH(t,s,𝐕)e^{Bt}Ae^{Bt}H(t,s,𝐕)=0.$$
(A8)
Using the property for scale transformation

$$e^{\frac{1}{2}\zeta ^{}t𝐕_𝐕}\chi (𝐕)=\chi \left(e^{\frac{1}{2}\zeta ^{}t}𝐕\right),$$
(A9)
(A8) simplifies to

$$_tH(t,s,𝐕)+e^{\frac{1}{2}\zeta ^{}t}V_y_sH(t,s,𝐕)=0.$$
(A10)
Integrating this with the initial condition $`H(0,s,𝐕)=h(s,𝐕)`$ gives

$$H(t,s,𝐕)=\mathrm{exp}\left[\frac{2}{\zeta ^{}}\left(1e^{\frac{1}{2}\zeta ^{}t}\right)V_y_s\right]h(s,𝐕).$$
(A11)
Use of this in (A7) gives the desired result

$`e^{(A+B)t}h(s,𝐕)`$ $`=`$ $`\mathrm{exp}\left({\displaystyle \frac{1}{2}}\zeta ^{}t𝐕_𝐕\right)\mathrm{exp}\left[{\displaystyle \frac{2}{\zeta ^{}}}\left(1e^{\frac{1}{2}\zeta ^{}t}\right)V_y_s\right]h(s,𝐕)`$ (A12)
$`=`$ $`\mathrm{exp}\left[\tau (t)V_y_s\right]h(s,e^{\frac{1}{2}\zeta ^{}t}𝐕)=h(s\tau (t)V_y,e^{\frac{1}{2}\zeta ^{}t}𝐕),`$ (A13)
where

$$\tau (t)=\frac{2}{\zeta ^{}}\left(e^{\frac{1}{2}\zeta ^{}t}1\right).$$
(A14)
Finally, combining Eqs. (A6) and (A13) gives the result

$$e^{𝒟t}F(s,𝐕)=e^{(A+B+C)t}F(s,𝐕)=F(s\tau (t)V_y,e^{\frac{1}{2}\zeta ^{}t}\left(𝐕+t𝖺𝐕\right)).$$
(A15)
## B Idealized boundary conditions

The boundary conditions for the general solution (30) are imposed by the choice of $`f_B`$, which is a solution to

$$\left(1\zeta ^{}\frac{d}{2}\zeta ^{}+𝒟\right)f_B=0$$
(B1)
or, equivalently,

$$_sf_B=V_y^1\left(1\zeta ^{}\frac{d}{2}\zeta ^{}aV_y_{V_x}\frac{\zeta ^{}}{2}𝐕_𝐕\right)f_B,$$
(B2)
where we have made use of Eq. (A3). This can be integrated to get

$`f_B(s,𝐕)`$ $`=`$ $`\mathrm{\Theta }(V_y)\mathrm{exp}\left[{\displaystyle \frac{s+_{}}{V_y}}\left(1\zeta ^{}{\displaystyle \frac{d}{2}}\zeta ^{}aV_y_{V_x}{\displaystyle \frac{\zeta ^{}}{2}}𝐕_𝐕\right)\right]\mathrm{\Phi }_{}(𝐕)`$ (B4)
$`+\mathrm{\Theta }(V_y)\mathrm{exp}\left[{\displaystyle \frac{s_+}{V_y}}\left(1\zeta ^{}{\displaystyle \frac{d}{2}}\zeta ^{}aV_y_{V_x}{\displaystyle \frac{\zeta ^{}}{2}}𝐕_𝐕\right)\right]\mathrm{\Phi }_+(𝐕).`$
Here $`\mathrm{\Theta }(x)`$ is the Heaviside step function, $`s=\pm _\pm `$ correspond to the location of the walls at $`y=\pm L`$, and $`\mathrm{\Phi }_\pm (𝐕)`$ are arbitrary functions. These functions are fixed by the boundary conditions at $`s=\pm _\pm `$:

$$\mathrm{\Theta }(\pm V_y)f(_{},𝐕)=\mathrm{\Theta }(\pm V_y)\mathrm{\Phi }_{}(𝐕).$$
(B5)
Idealized boundary conditions are constructed as follows. First we write the boundary conditions as

$$\mathrm{\Theta }(\pm V_y)f(_{},𝐕)=\mathrm{\Theta }(\pm V_y)\frac{n_{}}{C_{}}(m/2T_{})^{d/2}\phi _{}\left(𝐕/\sqrt{2T_{}/m}\right),$$
(B6)
where

$$n_{}\left(2T_{}/m\right)^{1/2}𝑑𝐕\mathrm{\Theta }(V_y)|V_y|f(_{},𝐕),$$
(B7)
$$C_{}𝑑𝝃\mathrm{\Theta }(\pm \xi _y)|\xi _y|\phi _{}(𝝃),$$
(B8)
$`\phi _{}(𝝃)`$ being the distributions of velocities away from the walls normalized as

$$𝑑𝝃\phi _{}(𝝃)=\frac{2}{d}𝑑𝝃\xi ^2\phi _{}(𝝃)=1.$$
(B9)
The usual choice is the Maxwell-Boltzmann distribution $`\phi _+(𝝃)=\phi _{}(𝝃)=\pi ^{d/2}\mathrm{exp}(\xi ^2)`$. Now consider the limit of infinitely cold walls, $`T_\pm 0`$. In that limit $`\mathrm{\Phi }_\pm (𝐕)0`$ since $`\phi _\pm (𝝃)`$ must vanish if $`|𝝃|\mathrm{}`$. Consequently,

$$\underset{T_\pm 0}{lim}f_B(s,𝐕)=0.$$
(B10)
## C Consistency conditions and calculation of the fluxes

### 1 Generating function

The consistency conditions for the assumed forms of the hydrodynamic fields are the $`d+2`$ equations of (24). Evaluation of the right hand sides of these equations requires the explicit form for the distribution function in terms of these fields. Using the results of Appendix A, the solution (30) with the idealized boundary conditions discussed above becomes

$`f(s,𝐕)`$ $`=`$ $`(1\zeta ^{}){\displaystyle _0^{\mathrm{}}}𝑑te^{(1\zeta ^{}\frac{d}{2}\zeta ^{})t}e^{atV_y_{V_x}}e^{\tau (t)V_y_s}f_{\mathrm{}}(s,e^{\frac{1}{2}\zeta ^{}t}𝐕)`$ (C1)
$`=`$ $`(1\zeta ^{}){\displaystyle _0^{\mathrm{}}}𝑑te^{(1\zeta ^{}\frac{d}{2}\zeta ^{})t}f_{\mathrm{}}(sV_y\tau (t),e^{\frac{1}{2}\zeta ^{}t}\left(𝐕+t𝖺𝐕\right)),`$ (C2)
where the local equilibrium distribution is

$$f_{\mathrm{}}(s,𝐕)=\frac{2p}{m}\mathrm{\Theta }(T(s))v_0^{(d+2)}(s)\pi ^{d/2}\mathrm{exp}[\left(V/v_0(s)\right)^2]$$
(C3)
with $`v_0^2(s)=2T(s)/m`$.

The velocity integrals for the consistency conditions and for the fluxes are all low order moments of this distribution function. They can be obtained from appropriate derivatives of the generating function

$$G(s,𝐤)=𝑑𝐕e^{i𝐤𝐕}f(s,𝐕).$$
(C4)
Making use of the solution (C2), the generating function becomes

$`G(s,𝐤)`$ $`=`$ $`(1\zeta ^{}){\displaystyle _0^{\mathrm{}}}𝑑te^{(1\zeta ^{}\frac{d}{2}\zeta ^{})t}{\displaystyle 𝑑𝐕e^{i𝐤𝐕}e^{atV_y_{V_x}}e^{\tau (t)V_y_s}f_{\mathrm{}}(s,e^{\frac{1}{2}\zeta ^{}t}𝐕)}`$ (C5)
$`=`$ $`(1\zeta ^{}){\displaystyle _0^{\mathrm{}}}𝑑te^{(1\zeta ^{})t}{\displaystyle 𝑑𝐕e^{i𝐤^{}(t)𝐕}e^{\tau _1(t)V_y_s}f_{\mathrm{}}(s,𝐕)},`$ (C6)
where

$$𝐤^{}(t)=e^{\frac{1}{2}\zeta ^{}t}e^{atk_x_{k_y}}𝐤,\tau _1(t)=\frac{2}{\zeta ^{}}\left(1e^{\frac{1}{2}\zeta ^{}t}\right).$$
(C7)
In the last step of Eq. (C5) we have made the change $`𝐕e^{\frac{1}{2}\zeta ^{}t}𝐕`$ and have used the general property

$$𝑑𝐕F_1(𝐕)e^{atV_y_{V_x}}F_2(𝐕)=𝑑𝐕F_2(𝐕)e^{atV_y_{V_x}}F_1(𝐕).$$
(C8)
Integration over $`𝐕_{}𝐕V_y\widehat{𝐲}`$ yields

$$G(s,𝐤)=G_+(s,𝐤)+G_{}(s,𝐤),$$
(C9)
where

$$G_\pm (s,𝐤)=\frac{2p}{m}(1\zeta ^{})_0^{\mathrm{}}𝑑V_y_0^{\mathrm{}}𝑑te^{(1\zeta ^{})t}e^{\pm ik_y^{}(t)V_y}e^{\tau _1(t)V_y_s}e^{V_y^2/v_0^2(s)}F(s,𝐤_{}^{}(t)),$$
(C10)
with

$$F(s,𝐤_{}^{}(t))=\mathrm{\Theta }(v_0^2(s))v_0^3(s)\pi ^{1/2}\mathrm{exp}\left[\frac{1}{4}k_{}^{}{}_{}{}^{2}(t)v_0^2(s)\right].$$
(C11)
We will focus on $`G_+(s,𝐤)`$ since $`G_{}(s,𝐤)`$ can be evaluated as $`G_{}(s,𝐤)=G_+(s,𝐤)`$.

Next, define the change of variables in the $`t`$ integration

$$z=V_y\tau _1(t),dz=\left(V_y\frac{\zeta ^{}z}{2}\right)dt,t=\frac{2}{\zeta ^{}}\mathrm{ln}\left(1\frac{\zeta ^{}z}{2V_y}\right).$$
(C12)
The function $`G_+(s,𝐤)`$ becomes

$`G_+(s,𝐤)`$ $`=`$ $`{\displaystyle \frac{2p}{m}}(1\zeta ^{}){\displaystyle _0^{\mathrm{}}}𝑑V_y{\displaystyle _0^{\mathrm{}}}𝑑ze^{(1\zeta ^{})t}e^{ik_y^{}(t)V_y}\left(V_y{\displaystyle \frac{\zeta ^{}}{2}}z\right)^1\mathrm{\Theta }\left(V_y{\displaystyle \frac{\zeta ^{}}{2}}z\right)`$ (C14)
$`\times e^{V_y^2/v_0^2(sz)}F(sz,𝐤_{}^{}(t))`$
$`=`$ $`{\displaystyle \frac{2p}{m}}(1\zeta ^{}){\displaystyle _0^{\mathrm{}}}𝑑ue^{u^2}u^1{\displaystyle _0^{\mathrm{}}}𝑑ze^{(1\zeta ^{})t}e^{ik_y^{}(t)uv_0(sz)}\left[1{\displaystyle \frac{\zeta ^{}}{2}}{\displaystyle \frac{z}{uv_0(sz)}}\right]^1`$ (C16)
$`\times \mathrm{\Theta }\left(1{\displaystyle \frac{\zeta ^{}}{2}}{\displaystyle \frac{z}{uv_0(sz)}}\right)F(sz,𝐤_{}^{}(t)).`$
A change of variables in the $`V_y`$ integral has been made, $`V_yuv_0(sz)`$. Accordingly, the variable $`t`$ becomes

$$t=\frac{2}{\zeta ^{}}\mathrm{ln}\left[1\frac{\zeta ^{}z}{2uv_0(sz)}\right].$$
(C17)
Next, for the $`z`$ integral change variables to

$$w=\frac{z}{v_0(sz)}.$$
(C18)
As a consequence, $`z`$ as a function of $`w`$ is $`z=z_+(s,w)`$, where $`z_+(s,w)`$ is the positive root of the quadratic equation obtained by squaring (C18):

$$\left(1+2\gamma w^2\right)z^24\gamma w^2szw^2v_0^2(s)=0.$$
(C19)
It can be easily shown that

$$\frac{dz}{dw}=v_0^2(s)\left(\frac{z_+}{w}\right)^3\left[1\frac{z_+}{2}_s\mathrm{ln}T(s)\right]^1.$$
(C20)
The generating function is now

$`G_+(s,𝐤)`$ $`=`$ $`{\displaystyle \frac{2p}{mv_0^2(s)\pi ^{1/2}}}(1\zeta ^{}){\displaystyle _0^{\mathrm{}}}𝑑ue^{u^2}u^1{\displaystyle _0^{\mathrm{}}}𝑑we^{(1\zeta ^{})t}\left(1{\displaystyle \frac{\zeta ^{}}{2}}{\displaystyle \frac{w}{u}}\right)^1\mathrm{\Theta }\left(1{\displaystyle \frac{\zeta ^{}}{2}}{\displaystyle \frac{w}{u}}\right)`$ (C22)
$`\times e^{ik_y^{}(t)z_+u/w}\left[1{\displaystyle \frac{z_+}{2}}_s\mathrm{ln}T(s)\right]^1e^{k_{}^{}{}_{}{}^{2}(t)(z_+/2w)^2},`$
where $`t`$ is now a function of $`u`$ and $`w`$,

$$t=\frac{2}{\zeta ^{}}\mathrm{ln}\left(1\frac{\zeta ^{}w}{2u}\right).$$
(C23)
Finally, change variables in the $`w`$ integration to $`t`$, so that

$$w=\frac{2u}{\zeta ^{}}\left(1e^{\frac{1}{2}\zeta ^{}t}\right),dt=u^1\left(1\frac{\zeta ^{}w}{2u}\right)^1dw.$$
(C24)
This leads to the explicit result

$$G_+(s,𝐤)=\frac{p}{T(s)\pi ^{1/2}}(1\zeta ^{})_0^{\mathrm{}}𝑑ue^{u^2}_0^{\mathrm{}}𝑑te^{(1\zeta ^{})t}A_+(s,𝐤,u,t),$$
(C25)
where

$$A_+(s,𝐤,u,t)=\left[1\frac{z_+(s,w)}{2}_s\mathrm{ln}T(s)\right]^1\mathrm{exp}\left[i\left(k_yatk_x\right)\frac{z_+(s,w)}{\tau (t)}\frac{k_{}^2}{4}\frac{z_+^2(s,w)}{u^2\tau ^2(t)}\right].$$
(C26)
It must be recalled that $`w`$ is a function of $`u`$ and $`t`$ given by (C24).

The function $`G_{}(s,𝐤)=G_+(s,𝐤)`$ is then

$$G_{}(s,𝐤)=\frac{p}{T(s)\pi ^{1/2}}(1\zeta ^{})_0^{\mathrm{}}𝑑ue^{u^2}_0^{\mathrm{}}𝑑te^{(1\zeta ^{})t}A_{}(s,𝐤,u,t),$$
(C27)
with

$$A_{}(s,𝐤,u,t)=\left[1\frac{z_{}(s,w)}{2}_s\mathrm{ln}T(s)\right]^1\mathrm{exp}\left[i\left(k_yatk_x\right)\frac{z_{}(s,w)}{\tau (t)}\frac{k_{}^2}{4}\frac{z_{}^2(s,w)}{u^2\tau ^2(t)}\right],$$
(C28)
where $`z_{}(s,w)=z_+(s,w)`$ is the negative root of the quadratic equation (C19). Real roots require

$$\left(4\gamma sw^2\right)^2+4\left(1+2\gamma w^2\right)2w^2\frac{T(s)}{m}0,$$
(C29)
which is satisfied for positive $`T(s)`$ and positive $`\gamma `$.

### 2 Consistency conditions

The consistency condition for the density becomes

$`n(s)`$ $`=`$ $`{\displaystyle 𝑑𝐯f}=G(s,𝐤=\mathrm{𝟎})`$ (C30)
$`=`$ $`{\displaystyle \frac{p}{T(s)\pi ^{1/2}}}(1\zeta ^{}){\displaystyle _0^{\mathrm{}}}𝑑ue^{u^2}{\displaystyle _0^{\mathrm{}}}𝑑te^{(1\zeta ^{})t}\left\{\left[1{\displaystyle \frac{z_+}{2}}_s\mathrm{ln}T(s)\right]^1+\left[1{\displaystyle \frac{z_{}}{2}}_s\mathrm{ln}T(s)\right]^1\right\}`$ (C31)
$`=`$ $`{\displaystyle \frac{p}{T(s)}}.`$ (C32)
Use has been made of the identity

$$\left[1\frac{z_+}{2}_s\mathrm{ln}T(s)\right]^1+\left[1\frac{z_{}}{2}_s\mathrm{ln}T(s)\right]^1=2,$$
(C33)
which follows from the explicit forms for the roots of (C19). The result (C32) is consistent with the required equation of state defining $`p`$.

The consistency condition for the temperature is

$$dn(s)T(s)=m𝑑𝐯V^2f=m\left[_𝐤^2G(s,𝐤)\right]_{𝐤=\mathrm{𝟎}}.$$
(C34)
Direct evaluation gives

$`1`$ $`=`$ $`{\displaystyle \frac{m}{dT(s)\pi ^{1/2}}}(1\zeta ^{}){\displaystyle _0^{\mathrm{}}}𝑑ue^{u^2}{\displaystyle _0^{\mathrm{}}}𝑑te^t\left[{\displaystyle \frac{d1}{2}}+\left(1+a^2t^2\right)u^2\right]w^2`$ (C36)
$`\times \left\{z_+^2\left[1{\displaystyle \frac{z_+}{2}}_s\mathrm{ln}T(s)\right]^1+z_{}^2\left[1{\displaystyle \frac{z_{}}{2}}_s\mathrm{ln}T(s)\right]^1\right\}.`$
Using the identity

$$z_+^2\left[1\frac{z_+}{2}_s\mathrm{ln}T(s)\right]^1+z_{}^2\left[1\frac{z_{}}{2}_s\mathrm{ln}T(s)\right]^1=\frac{4T(s)}{m}\frac{w^2}{1+2\gamma w^2},$$
(C37)
one gets the final result

$$1=\frac{4}{d\pi ^{1/2}}(1\zeta ^{})_0^{\mathrm{}}𝑑ue^{u^2}_0^{\mathrm{}}𝑑te^t\left[\frac{d1}{2}+\left(1+a^2t^2\right)u^2\right]\left(1+2\gamma w^2\right)^1.$$
(C38)
This condition is enforced by using it as the definition of $`\gamma `$.

Finally, the consistency conditions for the flow velocity components are

$$0=m𝑑𝐯V_if=\left[_{k_i}G(s,𝐤)\right]_{𝐤=\mathrm{𝟎}}.$$
(C39)
This is clearly satisfied for the $`z`$ component. The $`x`$ and $`y`$ components are proportional to integrals including the term

$$z_+\left[1\frac{z_+}{2}_s\mathrm{ln}T(s)\right]^1+z_{}\left[1\frac{z_{}}{2}_s\mathrm{ln}T(s)\right]^1,$$
(C40)
which identically vanishes.

This completes confirmation of the consistency conditions for the hydrodynamic fields. In summary, the distribution function (C2) with the hydrodynamic fields (32) constitutes an exact solution to the kinetic equation.

### 3 Calculation of the fluxes

The macroscopic transport properties of the steady Couette flow are given by the momentum and heat fluxes, Eq. (27). In terms of the generating function these are

$$P_{ij}=m\left[_{k_i}_{k_j}G(s,𝐤)\right]_{𝐤=\mathrm{𝟎}},𝐪(s)=i\frac{1}{2}m\left[_𝐤_𝐤^2G(s,𝐤)\right]_{𝐤=\mathrm{𝟎}}.$$
(C41)
Following similar steps as in the consistency condition for the temperature, the nonzero elements of the pressure tensor can be evaluated. They are given by

$$P_{xx}=\frac{2p}{\pi ^{1/2}}(1\zeta ^{})_0^{\mathrm{}}𝑑ue^{u^2}_0^{\mathrm{}}𝑑te^t\left[1+2\left(atu\right)^2\right]\left(1+2\gamma w^2\right)^1,$$
(C42)
$$P_{yy}=\frac{4p}{\pi ^{1/2}}(1\zeta ^{})_0^{\mathrm{}}𝑑ue^{u^2}u^2_0^{\mathrm{}}𝑑te^t\left(1+2\gamma w^2\right)^1,$$
(C43)
$$P_{zz}=\frac{2p}{\pi ^{1/2}}(1\zeta ^{})_0^{\mathrm{}}𝑑ue^{u^2}_0^{\mathrm{}}𝑑te^t\left(1+2\gamma w^2\right)^1,$$
(C44)
$$P_{xy}=\frac{4ap}{\pi ^{1/2}}(1\zeta ^{})_0^{\mathrm{}}𝑑ue^{u^2}u^2_0^{\mathrm{}}𝑑te^tt\left(1+2\gamma w^2\right)^1.$$
(C45)
It is readily verified that these results satisfy the relationship $`p=\frac{1}{d}\left[P_{xx}+P_{yy}+(d2)P_{zz}\right]`$.

Let us consider now the heat flux vector. Note that, although the temperature gradient is only directed along the $`y`$ direction, the presence of the shear flow induces a nonzero $`x`$ component of the heat flux. Thus the nonzero components are $`q_x`$ and $`q_y`$. Using (C41), $`q_y`$ is found to be

$`q_y(s)`$ $`=`$ $`{\displaystyle \frac{mp}{2T(s)\pi ^{1/2}}}(1\zeta ^{}){\displaystyle _0^{\mathrm{}}}𝑑ue^{u^2}u{\displaystyle _0^{\mathrm{}}}𝑑te^{\left(1+\frac{1}{2}\zeta ^{}\right)t}w^3\left[{\displaystyle \frac{d1}{2}}+\left(1+a^2t^2\right)u^2\right]`$ (C47)
$`\times \left\{z_+^3\left[1{\displaystyle \frac{z_+}{2}}_s\mathrm{ln}T(s)\right]^1+z_{}^3\left[1{\displaystyle \frac{z_{}}{2}}_s\mathrm{ln}T(s)\right]^1\right\}`$
$`=`$ $`{\displaystyle \frac{8\gamma ps}{\pi ^{1/2}}}(1\zeta ^{}){\displaystyle _0^{\mathrm{}}}𝑑ue^{u^2}u{\displaystyle _0^{\mathrm{}}}𝑑te^{\left(1+\frac{1}{2}\zeta ^{}\right)t}\left[{\displaystyle \frac{d1}{2}}+\left(1+a^2t^2\right)u^2\right]w\left(1+2\gamma w^2\right)^2,`$ (C48)
where use has been made of the identity

$$z_+^3\left[1\frac{z_+}{2}_s\mathrm{ln}T(s)\right]^1+z_{}^3\left[1\frac{z_{}}{2}_s\mathrm{ln}T(s)\right]^1=16\gamma s\frac{T(s)}{m}\frac{w^4}{(1+2\gamma w^2)^2}.$$
(C49)
This can be written as

$`q_y`$ $`=`$ $`{\displaystyle \frac{2ps}{\pi ^{1/2}}}(1\zeta ^{}){\displaystyle _0^{\mathrm{}}}𝑑ue^{u^2}{\displaystyle _0^{\mathrm{}}}𝑑te^t\left[{\displaystyle \frac{d1}{2}}+\left(1+a^2t^2\right)u^2\right]{\displaystyle \frac{}{t}}\left(1+2\gamma w^2\right)^1`$ (C50)
$`=`$ $`{\displaystyle \frac{dp\zeta ^{}}{2}}s+{\displaystyle \frac{4a^2ps}{\pi ^{1/2}}}(1\zeta ^{}){\displaystyle _0^{\mathrm{}}}𝑑ue^{u^2}u^2{\displaystyle _0^{\mathrm{}}}𝑑te^tt\left(1+2\gamma w^2\right)^1.`$ (C51)
An integration by parts in the $`t`$ integral and the consistency condition (C38) have been used to obtain the last equality. Comparison of this result with that for $`P_{xy}`$ above shows

$$q_y=\frac{1}{2m\gamma }\left(\frac{dp\zeta ^{}}{2}+aP_{xy}\right)_sT.$$
(C52)
This result is in fact required by the hydrodynamic equations to support the forms for the hydrodynamic fields.

In a similar way $`q_x(s)`$ is calculated,

$`q_x(s)`$ $`=`$ $`{\displaystyle \frac{amp}{2T(s)\pi ^{1/2}}}(1\zeta ^{}){\displaystyle _0^{\mathrm{}}}𝑑ue^{u^2}u{\displaystyle _0^{\mathrm{}}}𝑑te^{\left(1+\frac{1}{2}\zeta ^{}\right)t}tw^3\left[{\displaystyle \frac{d+1}{2}}+\left(1+a^2t^2\right)u^2\right]`$ (C54)
$`\times \left\{z_+^3\left[1{\displaystyle \frac{z_+}{2}}_s\mathrm{ln}T(s)\right]^1+z_{}^3\left[1{\displaystyle \frac{z_{}}{2}}_s\mathrm{ln}T(s)\right]^1\right\}`$
$`=`$ $`{\displaystyle \frac{4p}{m\pi ^{1/2}}}a(1\zeta ^{}){\displaystyle _0^{\mathrm{}}}𝑑ue^{u^2}u{\displaystyle _0^{\mathrm{}}}𝑑te^{\left(1+\frac{1}{2}\zeta ^{}\right)t}tw\left[{\displaystyle \frac{d+1}{2}}+\left(1+a^2t^2\right)u^2\right]`$ (C56)
$`\times \left(1+2\gamma w^2\right)^2_sT.`$
## D Transport coefficients in the uniform shear flow limit

In this Appendix we derive the explicit expressions for the transport coefficients along the critical shear rate $`a_c(\alpha )`$. They are obtained by taking the limit $`\gamma 0^+`$ in the corresponding expressions of Sec. III. Equations (36)–(41) simply yield

$$F_\eta (a_c,\alpha )=1\zeta ^{}=\frac{d}{d+2a_c^2},$$
(D1)
$$\mathrm{\Psi }_1(a_c,\alpha )=2(1\zeta ^{})=\frac{2d}{d+2a_c^2},$$
(D2)
$$\mathrm{\Psi }_2(a_c,\alpha )=0,$$
(D3)
$$\mathrm{\Phi }(a_c,\alpha )=\frac{2}{d+2}\frac{(1\zeta ^{})(4+3\zeta ^{})}{(1+\zeta ^{})^4(2+\zeta ^{})^4}\left[7(2+3\zeta ^{}+\zeta _{}^{}{}_{}{}^{2})^2+18a_c^2(8+12\zeta ^{}+5\zeta _{}^{}{}_{}{}^{2})\right].$$
(D4)
In order to obtain $`F_\lambda `$ from Eq. (39), we need to evaluate $`a^2`$ and $`F_\eta `$ to first order in $`\gamma `$. From Eqs. (33) and (36) one gets

$$a^2=\frac{d}{2}\frac{\zeta ^{}}{1\zeta ^{}}+A(\zeta ^{})\gamma +\mathrm{},$$
(D5)
$$F_\eta =1\zeta ^{}B(\zeta ^{})\gamma +\mathrm{},$$
(D6)
where

$`A(\zeta ^{})`$ $`=`$ $`{\displaystyle \frac{2}{(1+\zeta ^{})^3(2+\zeta ^{})^3}}[(d+2)(2+3\zeta ^{}+\zeta _{}^{}{}_{}{}^{2})^2`$ (D8)
$`+6a_c^2(24+48\zeta ^{}+33\zeta _{}^{}{}_{}{}^{2}+9\zeta _{}^{}{}_{}{}^{3}+\zeta _{}^{}{}_{}{}^{4})],`$
$$B(\zeta ^{})=12\frac{1\zeta ^{}}{(1+\zeta ^{})^2(2+\zeta ^{})^2}(6+6\zeta ^{}+\zeta _{}^{}{}_{}{}^{2}).$$
(D9)
Thus, the thermal conductivity is

$$F_\lambda (a_c,\alpha )=\frac{(1\zeta ^{})A(\zeta ^{})a_c^2B(\zeta ^{})}{d+2}.$$
(D10)
Equations (D1)–(D3) coincide with those previously derived in Ref. . However, Eqs. (D4) and (D10) are new results. Despite the fact that there is no heat flux in the uniform shear flow, Eqs. (D4) and (D10) are intrinsic transport coefficients characterizing the state of the system.