text:
# Combinatorial Identities from the Spectral Theory of Quantum Graphs

## I Introduction

In the present paper we show that some questions arising in the study of spectral correlations for quantum graphs can be cast as combinatorial problems. Solving these problems for a particular system, we discovered the following novel combinatorial identities:

(i) Let $`n,q`$ be arbitrary integers with $`1q<n`$ and

$$F_{\nu ,\nu ^{}}(n,q)=\frac{(n1)n}{2}\frac{(1)^{\nu +\nu ^{}}}{\nu \nu ^{}}\left(\genfrac{}{}{0pt}{}{n}{\nu +\nu ^{}}\right)^1\left(\genfrac{}{}{0pt}{}{q1}{\nu 1}\right)\left(\genfrac{}{}{0pt}{}{q1}{\nu ^{}1}\right)\left(\genfrac{}{}{0pt}{}{nq1}{\nu 1}\right)\left(\genfrac{}{}{0pt}{}{nq1}{\nu ^{}1}\right).$$

Then

$$S(n,q)=\underset{\nu ,\nu ^{}=1}{\overset{\mathrm{min}(q,nq)}{}}F_{\nu ,\nu ^{}}(n,q)=1.$$
(1)
(ii) Let $`s,t`$ be arbitrary positive integers and

$$𝒩(s,t)=\underset{\nu =1}{\overset{\mathrm{min}(s,t)}{}}(1)^{t\nu }\left(\genfrac{}{}{0pt}{}{t}{\nu }\right)\left(\genfrac{}{}{0pt}{}{s1}{\nu 1}\right)=(1)^{s+t}\left(\genfrac{}{}{0pt}{}{s+t1}{s}\right)^{1/2}P_{s+t1,s}(t),$$
(2)
where $`P_{N,k}(x)`$ are the Kravtchouk polynomials to be defined in Eq. (55). Further, let $`x,y`$ be complex with $`|x|,|y|<1/\sqrt{2}`$. Then we have the generating functions

$$G_1(x)=\underset{s,t=1}{\overset{\mathrm{}}{}}𝒩^2(s,t)x^{s+t}=\frac{x}{2x1}\left(\frac{1}{\sqrt{4x^2+1}}\frac{1}{1x}\right),$$
(3)
$$G_2(x)=\underset{s,t=1}{\overset{\mathrm{}}{}}𝒩(s,t)𝒩(t,s)x^{s+t}=\frac{1}{2}\frac{4x^2+2x+1}{(2x+1)\sqrt{4x^2+1}}\frac{1}{2}$$
(4)
and

$$g(x,y)=\underset{s,t=1}{\overset{\mathrm{}}{}}𝒩(s,t)x^sy^t\frac{xy}{(1+y)(1x+y2xy)}.$$
(5)
(iii) Let $`m`$ be any positive integer. Then

$$4m^2\underset{q=1}{\overset{2m1}{}}\left(\frac{𝒩(s,t)}{q}\right)^2=2^{2m+1}+(1)^m\left(\genfrac{}{}{0pt}{}{2m}{m}\right)2$$
(6)
and

$$(2m+1)^2\underset{q=1}{\overset{2m}{}}\left(\frac{𝒩(s,t)}{q}\right)^2=2^{2m+2}2(1)^m\left(\genfrac{}{}{0pt}{}{2m}{m}\right)2.$$
(7)
(iv) Let $`0qn`$, and define

$$A(n,q)=\frac{1}{\sqrt{2^n}}\{\begin{array}{cc}1\hfill & \mathrm{for}q=0,n\hfill \\ (1)^q(n/q)𝒩(nq,q)\hfill & \mathrm{for}\mathrm{\hspace{0.33em}0}<q<n.\hfill \end{array}$$
(8)
Then for any positive integers $`0\kappa \nu `$ and an arbitrary integer $`n_0`$,

$$\underset{ϵ0}{lim}ϵ\underset{n=n_0}{\overset{\mathrm{}}{}}\mathrm{e}^{nϵ}\underset{q=0}{\overset{n}{}}A(n+\nu ,q+\kappa )A(n,q)=A(\nu ,\kappa ).$$
(9)
These identities establish a new connection between combinatorics and the theory of random ensembles of $`2\times 2`$ matrices. The physical background and motivations are described in a few recent publications , to which the interested reader is referred. An immediate application of (1) is also given in . Here, we shall provide the minimum background necessary for the understanding of the combinatorial aspects of the problem, and for a self-contained exposition of our results. It will also enable us to propose a conjecture which generalizes the combinatorial approach to random matrix theory for matrices of large dimensions.

We start with a few definitions: Graphs consist of $`V`$ vertices connected by $`B`$ bonds (or edges). The valency $`v_i`$ of a vertex $`i`$ is the number of bonds meeting at that vertex. We shall assume that two vertices can be connected by a single bond at most. We denote the bonds connecting the vertices $`i`$ and $`j`$ by $`b=[i,j]`$. The notation $`[i,j]`$ will be used whenever we do not need to specify the direction on the bond. Hence $`[i,j]=[j,i]`$. Directed bonds will be denoted by $`d=(i,j)`$, and we shall always use the convention that the bond is directed from the first to the second index. If $`d=(i,j)`$ we use the notation $`\widehat{d}=(j,i)`$. Let $`g^{(i)}`$ be the set of directed bonds $`(i,j)`$ which emanate from the vertex $`i`$, and $`\widehat{g}^{(i)}`$ the set of directed bonds $`(j^{},i)`$ which converge at $`i`$. The vertices $`i`$ and $`j`$ are connected if $`g^{(i)}\widehat{g}^{(j)}\mathrm{}`$. The bond $`d^{}`$ is connected to the bond $`d`$ if there is some vertex $`i`$ with $`dg^{(i)}`$ and $`d^{}\widehat{g}^{(i)}`$. $`d`$ and $`\widehat{d}`$ are always connected.

The Schrödinger operator on the graph is defined after the natural metric is assigned to the bonds, and the solutions of the one-dimensional Schrödinger equation on each bond $`(\text{i}\mathrm{d}_xA)^2\psi (x)=k^2\psi (x)`$ is given as a linear combination of counter-propagating waves. $`A`$ stands here for a magnetic flux. The (complex) amplitudes of the counter propagating waves are denoted by $`a_d`$, where the subscript $`d`$ stands for the directed bond along which the wave propagates, $`d=1,..,2B`$. Appropriate boundary conditions at the vertices are imposed, and the spectrum of the Schrödinger operator on the graph is determined as the (infinite, discrete) set of energies $`k_n^2`$, for which there exists a non-trivial set of $`a_d`$ which is consistent with the boundary conditions. The condition of consistency can be expressed by the requirement that

$$det(IS(k))=0$$
(10)
where the bond-scattering matrix $`S(k)`$ is a unitary operator in the Hilbert space of $`2B`$ dimensional vectors of coefficients $`a_d`$. The unitarity of $`S(k)`$ ensures that the spectrum of the Schrödinger operator is real. The matrix $`S(k)`$, which is the object of our study is defined as

$$S_{(i,j),(l,m)}(k)=\delta _{j,l}\mathrm{e}^{\text{i}\varphi _{(i,j)}(k)}\sigma _{i,m}^{(j)}(k).$$
(11)
The matrix elements $`S_{d,d^{}}(k)`$ vanish if the bonds are not connected. As a consequence, the unitarity of $`S`$ implies also the unitarity of the $`v_j`$-dimensional vertex-scattering matrices $`\sigma _{i,m}^{(j)}`$ and vice versa. The phases $`\varphi _{(i,j)}`$, are given in terms of the bond length $`L_{[i,j]}`$, and the magnetic flux $`A_{(i,j)}=A_{(j,i)}`$ ,

$$\varphi _{(i,j)}(k)=(k+A_{(i,j)})L_{[i,j]}.$$
(12)
The two phases pertaining to the same bond $`\varphi _d`$ and $`\varphi _{\widehat{d}}`$ are equal when $`A_d=A_{\widehat{d}}=0`$. In this case $`S`$ is symmetric and the Schrödinger operator on the graph is invariant under time reversal. Time-reversal symmetry is violated when some magnetic fluxes do not vanish.

$`S`$ can also be interpreted as a quantum time evolution operator describing the scattering of waves with wave number $`k`$ between connected bonds. The wave gains the phase $`\varphi _{(i,j)}(k)`$ during the propagation along the bond $`(i,j)`$, while the $`\sigma _{i,m}^{(j)}`$ describe the scattering at the vertices. In this picture, the unitarity of $`S`$ guarantees the conservation of probability during the time evolution.

We will avoid unnecessary technical difficulties and consider the matrices $`\sigma _{i,m}^{(j)}`$ to be $`k`$ independent constants. One may find explicit expressions for $`\sigma _{i,m}^{(j)}`$ by requiring besides unitarity that the wave function is continuous at the vertices. The resulting expression is

$$\sigma _{j,j^{}}^{(i)}=\frac{2}{v_i}\delta _{j,j^{}}\text{(Neumann b. c.)}.$$
(13)
Note that back-scattering ($`j=j^{}`$) is singled out both in sign and in magnitude. In all nontrivial cases ($`v_i>2`$) the back-scattering amplitude is negative, and the back-scattering probability $`|\sigma _{j,j}^{(i)}|^2`$ approaches $`1`$ as the valency $`v_i`$ increases.

Finally, a “classical analogue” of the quantum dynamics can be defined as a random walk on the directed bonds, in which the transition probability between bonds $`(i,j)`$, $`(j,l)`$ connected at vertex $`j`$ is $`|\sigma _{i,l}^{(j)}|^2`$. The resulting classical evolution operator with matrix elements

$$U_{d,d^{}}=|S_{d,d^{}}|^2$$
(14)
is probability conserving, since unitarity implies $`_i|\sigma _{i,l}^{(j)}|^2=1`$.

## II Spectral Two-Point Correlations and Periodic-Orbit Sums

The spectrum of $`S`$ consists of $`2B`$ eigenvalues $`\mathrm{e}^{i\theta _l}`$ which are confined to the unit circle. Their distribution is given in terms of the spectral density

$$d(\theta )\underset{l=1}{\overset{2B}{}}\delta _{2\pi }(\theta \theta _l)=\frac{2B}{2\pi }+\frac{1}{2\pi }\underset{n=1}{\overset{\mathrm{}}{}}s_n\text{e}^{\text{i}\theta n}+\mathrm{c}.\mathrm{c}.,$$
(15)
where $`\delta _{2\pi }`$ denotes the $`2\pi `$ periodic delta function. The first term on the r.h.s. is the average density $`\overline{d}=\frac{2B}{2\pi }`$. The coefficients of the oscillatory part

$$s_n=\mathrm{tr}S^n$$

will play an important rôle in the following. $`s_n=\mathrm{tr}S^n`$ is a sum over products of $`n`$ matrix elements of $`S`$, and because of (11) the bond indices of each summand describe a connected $`n`$-cycle ($``$ $`n`$-periodic orbit) on the graph

$$s_n=\underset{p𝒫_n}{}𝒜_p\mathrm{e}^{i\mathrm{\Phi }_p}.$$
(16)
In (16) $`𝒫_n`$ denotes the set of all $`n`$-periodic orbits (PO’s) on the graph. Note that for the convenience of presentation we will consider cycles differing only by a cyclic permutation as different PO’s. The phases $`\mathrm{\Phi }_p=_{j=0}^{n1}\varphi _{d_j}`$ can be interpreted as the action along the PO $`p`$. The amplitudes $`𝒜_p`$ are given by

$$𝒜_p=\underset{j=0}{\overset{n1}{}}S_{d_{j+1},d_j},$$
(17)
where $`j`$ is understood mod $`n`$. Sometimes it is useful to split the amplitude in its absolute value and a phase $`\mathrm{e}^{\text{i}\mu _p\pi }`$. For example, in the case of Neumann b. c. (13) $`\mu _p`$ is an integer counting the number of back-scatterings along $`p`$.

In complete analogy to (16) we can represent also the traces of powers of the classical evolution operator

$$u_n=\mathrm{tr}U^n$$
(18)
as sums of periodic orbits of the graph.

The two-point correlations in the spectrum of $`S`$ (15) can be expressed in terms of the average excess probability density $`R_2(r;\beta )`$ of finding two phases at a distance $`r`$, where $`r`$ is measured in units of the mean spacing $`\frac{2\pi }{2B}`$,

$`R_2(r;\beta )`$ $`=`$ $`{\displaystyle \frac{1}{2B}}{\displaystyle _\pi ^{+\pi }}d\theta d\left(\theta \right)d\left(\theta [\pi /B]r\right)_\beta {\displaystyle \frac{B}{\pi }}`$ (19)
$`=`$ $`{\displaystyle \frac{2}{2\pi }}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}\mathrm{cos}\left({\displaystyle \frac{\pi }{B}}nr\right){\displaystyle \frac{1}{2B}}|s_n|^2_\beta .`$ (20)
The bond scattering matrix depends parametrically on the phases $`\varphi _d`$ (12). We shall define two statistical ensembles for $`S`$ in the following way. The ensemble for which time-reversal symmetry is broken consists of $`S`$ matrices for which the $`\varphi _d`$ are all different, and we consider them as independent variables distributed uniformly on the $`2B`$ torus. Invariance under time reversal implies $`\varphi _d=\varphi _{\widehat{d}}`$ and the corresponding ensemble is defined in terms of $`B`$ independent and uniformly distributed phases. We shall distinguish between these ensembles by the value of the parameter $`\beta =\{`$number of independent phases$`\}/B`$. Expectation values with respect to these measures are denoted in (20) by triangular brackets,

$$\mathrm{}_\beta \underset{d}{\overset{\beta B}{}}\left(\frac{1}{2\pi }_\pi ^{+\pi }d\varphi _d\right)\mathrm{}.$$
(21)
The Fourier transform of $`R_2(r;\beta )`$ is the form factor

$$K(n/2B;\beta )=\frac{1}{2B}|s_n|^2_\beta $$
(22)
on which our interest will be focussed. If the eigenvalues of the $`S`$ were statistically independent and uniformly distributed on the unit circle, $`K(n/2B)=1`$ for all $`n`$. Any deviation of the form factor from unity implies spectral correlations. Using (16) the form factor (22) is expressed as a double sum over PO’s

$`K(n/2B;\beta )`$ $`=`$ $`{\displaystyle \frac{1}{2B}}\left|{\displaystyle \underset{p𝒫_n}{}}𝒜_p\text{e}^{\text{i}\mathrm{\Phi }_p}\right|^2_\beta `$ (23)
$`=`$ $`{\displaystyle \frac{1}{2B}}{\displaystyle \underset{p,p^{}𝒫_n}{}}𝒜_p𝒜_p^{}\text{e}^{\text{i}(\mathrm{\Phi }_p\mathrm{\Phi }_p^{})}_\beta `$ (24)
In order to perform the average over all the phases $`\varphi _d`$ in (23) we write

$`\mathrm{\Phi }_p={\displaystyle \underset{d}{}}n_d^{(p)}\varphi _d,`$ (25)
where $`n_d^{(p)}`$ counts the number of traversals of each directed bond such that $`_dn_d^{(p)}=n`$. According to (21) we have

$`\text{e}^{\text{i}(n\varphi _d+n^{}\varphi _d^{})}_{\beta =1}`$ $`=`$ $`\delta _{n,0}\delta _{n^{},0},`$ (26)
$`\text{e}^{\text{i}(n\varphi _d+\widehat{n}\varphi _{\widehat{d}})}_{\beta =2}`$ $`=`$ $`\delta _{n+\widehat{n},0}.`$ (27)
Thus, the double sum in (23) can be restricted to families of orbits. For $`\beta =2`$, let $``$ be the family of isometric PO’s which have the same integers $`n_d^{()}`$. That is, the family consists of all the PO’s which traverse the same directed bonds the same number of times, but not necessarily in the same order. In the case $`\beta =1`$, $``$ consists of all PO’s sharing $`n_b^{()}n_d^{()}+n_{\widehat{d}}^{()}`$. That is, the family contains all PO’s which traverse the same set of undirected bonds the same number of times, irrespective of direction or order. We find

$`K(n/2B;\beta )`$ $`=`$ $`{\displaystyle \frac{1}{2B}}{\displaystyle \underset{(\beta )_n}{}}|{\displaystyle \underset{p}{}}𝒜_p|^2.`$ (28)
$`(\beta )_n`$ denotes the set of all vectors $`=[n_d]`$ for $`\beta =2`$ ($`=[n_b]`$ for $`\beta =1`$) of $`\beta B`$ non-negative integers summing to $`n`$, for which at least one PO exists. For Neumann b. c. (13), e. g., (28) amounts to counting the PO’s in a given set $`=[n_d]`$ taking into account the number of back-scatterings along the orbit. The problem of spectral statistics is now reduced to a counting (combinatorial) problem which is, however, very complicated in general. Even the determination of the number of families $``$ for a given $`n`$ is difficult. For $`\beta =2`$ an obvious necessary condition for the existence of a PO with a given set of bond traversals $`=[n_d]`$ is that at any vertex the number of incoming and outgoing bonds is the same, i. e.

$$\underset{dg^{(i)}}{}n_d=\underset{d\widehat{g}^{(i)}}{}n_d(i=1,\mathrm{}V).$$

For $`\beta =1`$, the analogous condition reads

$$\underset{dg^{(i)}}{}n_d\underset{d\widehat{g}^{(i)}}{}n_d\text{ mod }2=0(i=1,\mathrm{}V),$$

i. e. the total number of traversals of adjacent bonds should be even at each vertex. However, it is not so easy to formulate a sufficient condition for the existence of a PO given a set of numbers $`n_d`$. In particular one must take care to exclude cases, in which the set of traversed bonds is a union of two or more disconnected groups (“composite orbits”).

Extensive numerical work  revealed that for fully connected graphs ($`v_jV1`$), and for $`V1`$, the form-factor (22) is well reproduced by the predictions of random matrix theory  for the Circular Orthogonal Ensemble (COE) ($`\beta =1`$) or the Circular Unitary Ensemble (CUE) ($`\beta =2`$). This leads us to expect that (28) approaches the corresponding random matrix prediction in the limit $`V\mathrm{}`$. This conjecture is proposed as a challenge to asymptotic combinatorial theory.

## III The Ring Graph

In the following we will evaluate explicitly the quantities introduced in the previous section for one of the simplest quantum graphs. It consists of a single vertex on a loop (see fig. 1). There are two directed bonds $`d=1`$ and $`\widehat{d}=2`$ with $`\varphi _1\varphi _2`$, i. e. time-reversal symmetry is broken. Since this graph would be trivial for Neumann b.c. the vertex-scattering matrix at the only vertex is chosen as

$$\sigma (\eta )=\left(\begin{array}{cc}\mathrm{cos}\eta \hfill & \mathrm{i}\mathrm{sin}\eta \hfill \\ \mathrm{i}\mathrm{sin}\eta \hfill & \mathrm{cos}\eta \hfill \end{array}\right),$$
(29)
with $`0\eta \pi /2`$. The corresponding bond-scattering matrix is

$$S(\eta )=\left(\begin{array}{cc}\mathrm{e}^{\varphi _1}\hfill & 0\hfill \\ 0\hfill & \mathrm{e}^{\varphi _2}\hfill \end{array}\right)\left(\begin{array}{cc}\mathrm{cos}\eta \hfill & \mathrm{i}\mathrm{sin}\eta \hfill \\ \mathrm{i}\mathrm{sin}\eta \hfill & \mathrm{cos}\eta \hfill \end{array}\right).$$
(30)
We shall compute the form factor for two ensembles. The first is defined by a fixed value of $`\eta =\pi /4`$, and the only average is over the phases $`\varphi _d`$ according to (26). The second ensemble includes an additional averaging over the parameter $`\eta `$. We will show that the measure for the integration over $`\eta `$ can be chosen such that the model yields exactly the CUE form factor for 2$`\times `$ 2 random matrices .

### A Periodic Orbit Representation of $`u_n`$

We will first illustrate our method of deriving combinatorial results from the ring graph in a case where a known identity is obtained. Consider the classical evolution operator $`U`$ of the ring graph. According to (14) we have

$$U(\eta )=\left(\begin{array}{cc}\mathrm{cos}^2\eta \hfill & \mathrm{sin}^2\eta \hfill \\ \mathrm{sin}^2\eta \hfill & \mathrm{cos}^2\eta \hfill \end{array}\right).$$
(31)
The spectrum of $`U`$ consists of $`\{1,\mathrm{cos}2\eta \}`$, such that

$$u_n(\eta )=1+\mathrm{cos}^n2\eta .$$
(32)
We will now show how this result can be obtained from a sum over the periodic orbits of the system, grouped into families of orbits as in (28). In the classical calculation it is actually not necessary to take the families into account, but we would like to stress the analogy to the quantum case considered below. The periodic orbit expansion of the classical return probability can easily be obtained from (31) by expanding all matrix products in (18). We find

$`u_n`$ $`=`$ $`{\displaystyle \underset{i_1=1,2}{}}\mathrm{}{\displaystyle \underset{i_n=1,2}{}}{\displaystyle \underset{j=0}{\overset{n1}{}}}U_{i_j,i_{j+1}}(\eta ),`$ (33)
where $`j`$ is again taken mod $`n`$. In the following the binary sequence $`[i_j]`$ ($`i_j\{1,2\}`$; $`j=0,\mathrm{},n1`$) is referred to as the code of the orbit. We will now sort the terms in the multiple sum above into families of isometric orbits. In the present case a family is completely specified by the integer $`qq_1`$ which counts the traversals of the loop $`1`$, i.e., the number of letters $`1`$ in the code word. Each of these $`q`$ letters is followed by an uninterrupted sequence of $`t_j0`$ letters $`2`$ with the restriction that the total number of letters $`2`$ is given by

$$\underset{j=1}{\overset{q}{}}t_j=nq.$$
(34)
We conclude that each code word in a family $`0<q<n`$ which starts with $`i_1=1`$ corresponds to an ordered partition of the number $`nq`$ into $`q`$ non-negative integers, while the words starting with $`i_1=2`$ can be viewed as partition of $`q`$ into $`nq`$ summands.

To make this step very clear, consider the following example: All code words of length $`n=5`$ in the family $`q=2`$ are $`11222`$, $`12122`$, $`12212`$, $`12221`$ and $`22211`$, $`22121`$, $`21221`$, $`22112`$, $`21212`$, $`21122`$. The first four words correspond to the partitions $`0+3=1+2=2+1=3+0`$ of $`nq=3`$ into $`q=2`$ terms, while the remaining $`5`$ words correspond to $`2=0+0+2=0+1+1=1+0+1=0+2+0=1+1+0=2+0+0`$.

In the multiple products in (33) a backward scattering along the orbit is expressed by two different consecutive symbols $`i_ji_{j+1}`$ in the code and leads to a factor $`\mathrm{sin}^2\eta `$, while a forward scattering contributes a factor $`\mathrm{cos}^2\eta `$ . Since the sum is over periodic orbits, the number of back scatterings is always even and we denote it with $`2\nu `$. It is then easy to see that $`\nu `$ corresponds to the number of positive terms in the partitions introduced above, since each such term corresponds to an uninterrupted sequence of symbols $`2`$ enclosed between two symbols $`1`$ or vice versa and thus contributes two back scatterings. For the codes starting with a symbol $`1`$ there are $`\left(\genfrac{}{}{0pt}{}{q}{\nu }\right)`$ ways to choose the $`\nu `$ positive terms in the sum of $`q`$ terms, and there are $`\left(\genfrac{}{}{0pt}{}{nq1}{\nu 1}\right)`$ ways to decompose $`nq`$ into $`\nu `$ positive summands. After similar reasoning for the codes starting with the symbol $`2`$ we find for the periodic orbit expansion of the classical return probability

$`u_n(\eta )`$ $`=`$ $`2\mathrm{cos}^{2n}\eta +{\displaystyle \underset{q=1}{\overset{n1}{}}}{\displaystyle \underset{\nu }{}}\left[\left({\displaystyle \genfrac{}{}{0pt}{}{q}{\nu }}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{nq1}{\nu 1}}\right)+\left({\displaystyle \genfrac{}{}{0pt}{}{nq}{\nu }}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{q1}{\nu 1}}\right)\right]\mathrm{sin}^{4\nu }\eta \mathrm{cos}^{2n4\nu }\eta `$ (35)
$`=`$ $`2\mathrm{cos}^{2n}\eta +{\displaystyle \underset{q=1}{\overset{n1}{}}}{\displaystyle \underset{\nu }{}}{\displaystyle \frac{n}{\nu }}\left({\displaystyle \genfrac{}{}{0pt}{}{q1}{\nu 1}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{nq1}{\nu 1}}\right)\mathrm{sin}^{4\nu }\eta \mathrm{cos}^{2n4\nu }\eta `$ (36)
The summation limits for the variable $`\nu `$ are implicit since all terms outside vanish due to the properties of the binomial coefficients. From the equivalence between (32) and (36) the combinatorial identity

$$\underset{q=1}{\overset{n1}{}}\left(\genfrac{}{}{0pt}{}{q1}{\nu 1}\right)\left(\genfrac{}{}{0pt}{}{nq1}{\nu 1}\right)=\left(\genfrac{}{}{0pt}{}{n1}{2\nu 1}\right)=\frac{2\nu }{n}\left(\genfrac{}{}{0pt}{}{n}{2\nu }\right).$$
(37)
could be deduced which indeed reduces (36) to a form

$`u_n(\eta )`$ $`=`$ $`2{\displaystyle \underset{\nu }{}}\left({\displaystyle \genfrac{}{}{0pt}{}{n}{2\nu }}\right)\mathrm{sin}^{4\nu }\eta \mathrm{cos}^{2n4\nu }\eta `$ (38)
$`=`$ $`(\mathrm{cos}^2\eta +\mathrm{sin}^2\eta )^n+(\mathrm{cos}^2\eta \mathrm{sin}^2\eta )^n,`$ (39)
which is obviously equivalent to (32).

(37) can also be derived by some straight forward variable substitutions from the identity

$$\underset{k=l}{\overset{nm}{}}\left(\genfrac{}{}{0pt}{}{k}{l}\right)\left(\genfrac{}{}{0pt}{}{nk}{m}\right)=\left(\genfrac{}{}{0pt}{}{n+1}{l+m+1}\right)$$
(40)
which is found in the literature .

### B Quantum Mechanics: Spacing Distribution and Form Factor

In the following two subsections we derive novel combinatorial identities by applying the reasoning which led to (37) to the quantum evolution operator (30) of the ring graph. We can write the eigenvalues of $`S(\eta )`$ as $`\mathrm{e}^{\text{i}(\varphi _1+\varphi _2)/2}\mathrm{e}^{\pm \text{i}\lambda /2}`$ with

$$\lambda =2\mathrm{arcos}\left[\mathrm{cos}\eta \mathrm{cos}\left(\frac{\varphi _1\varphi _2}{2}\right)\right]$$
(41)
denoting the difference between the eigenphases. For the two-point correlator we find

$`R_2(r,\eta )`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle _\pi ^{+\pi }}d\theta d\left(\theta +{\displaystyle \frac{\pi r}{2}}\right)d\left(\theta {\displaystyle \frac{\pi r}{2}}\right)_{\varphi _{1,2}}{\displaystyle \frac{1}{\pi }}`$ (42)
$`=`$ $`{\displaystyle \frac{\delta _2(r)1}{\pi }}{\displaystyle \frac{\delta _{2\pi }\left(\pi r+\lambda \right)+\delta _{2\pi }\left(\pi r\lambda \right)}{2}}_{\varphi _{1,2}}`$ (43)
$`=`$ $`{\displaystyle \frac{\delta _2\left(r\right)1}{\pi }}+{\displaystyle \frac{\mathrm{sin}|\pi r/2|}{2\pi }}{\displaystyle \frac{\mathrm{\Theta }(\mathrm{cos}^2\eta \mathrm{cos}^2(\pi r/2))}{\sqrt{\mathrm{cos}^2\eta \mathrm{cos}^2(\pi r/2)}}})`$ (44)
Here, $`\delta _2\left(r\right)`$ is the $`2`$-periodic $`\delta `$ function. In particular for equal transmission and reflection probability ($`\eta =\pi /4`$) we have

$`R_2(r,\pi /4)`$ $`=`$ $`{\displaystyle \frac{\delta _2\left(r\right)1}{\pi }}+{\displaystyle \frac{1}{2\pi }}\sqrt{{\displaystyle \frac{\mathrm{cos}(\pi r)1}{\mathrm{cos}(\pi r)}}}\mathrm{\Theta }\left({\displaystyle \frac{1}{2}}|r1|\right)`$ (45)
and, by a Fourier transformation, we can compute the form factor

$`K(n,\pi /4)`$ $`=`$ $`\pi {\displaystyle _0^2}dr\mathrm{cos}\left(n\pi r\right)R_2(r,\pi /4)`$ (46)
$`=`$ $`1+{\displaystyle \frac{(1)^{m+n}}{2^{2m+1}}}\left({\displaystyle \genfrac{}{}{0pt}{}{2m}{m}}\right){\displaystyle \frac{3}{2}}\delta _{n,0}`$ (47)
$``$ $`1+{\displaystyle \frac{(1)^{m+n}}{2\sqrt{\pi n}}}(n1),`$ (48)
where $`m=[n/2]`$ and $`[]`$ stands for the integer part.

Next we consider the ensemble for which transmission and reflection probabilities are uniformly distributed between $`0`$ and $`1`$. For the parameter $`\eta `$ this corresponds to the measure $`\mathrm{d}\mu (\eta )=2|\mathrm{cos}\eta \mathrm{sin}\eta |\mathrm{d}\eta `$. The main reason for this choice is that upon integrating (44) one gets

$`R_2^{(\mathrm{av})}(r)`$ $`=`$ $`{\displaystyle \frac{\delta _2\left(r\right)1}{\pi }}+{\displaystyle \frac{\mathrm{sin}^2(\pi r/2)}{\pi }}`$ (49)
which coincides with the CUE result for $`2\times 2`$ matrices. A Fourier transformation results in

$$K_2(n)=\{\begin{array}{cc}\frac{1}{2}\hfill & \mathrm{for}n=1\hfill \\ 1\hfill & \mathrm{for}n2\hfill \end{array}.$$
(50)
The form factors (47), (48) and (50) are displayed in Fig. 1 below.

### C Periodic Orbit Expansion of the Form Factor

An explicit formulation of (28) for the ring graph is found by labelling and grouping orbits as explained in the derivation of (36). We obtain

$`K_2(n;\eta )`$ $`=`$ $`\mathrm{cos}^{2n}\eta +{\displaystyle \frac{n^2}{2}}{\displaystyle \underset{q=1}{\overset{n1}{}}}\left[{\displaystyle \underset{\nu }{}}{\displaystyle \frac{(1)^\nu }{\nu }}\left({\displaystyle \genfrac{}{}{0pt}{}{q1}{\nu 1}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{nq1}{\nu 1}}\right)\mathrm{sin}^{2\nu }\eta \mathrm{cos}^{n2\nu }\eta \right]^2,`$ (51)
where $`q`$ denotes the number of traversals of the ring in positive direction and $`2\nu `$ is the number of backward scatterings along the orbit. The inner sum over $`\nu `$ can be written in terms of Kravtchouk polynomials as

$`K_2(n;\eta )`$ $`=`$ $`\mathrm{cos}^{2n}\eta +{\displaystyle \frac{1}{2}}{\displaystyle \underset{q=1}{\overset{n1}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{n1}{nq}}\right)\mathrm{cos}^{2q}\eta \mathrm{sin}^{2(nq)}\eta \left[{\displaystyle \frac{n}{q}}P_{n1,nq}^{(\mathrm{cos}^2\eta ,\mathrm{sin}^2\eta )}(q)\right]^2,`$ (52)
and the Kravtchouk polynomials are defined as in  by

$`P_{N,k}^{(u,v)}(x)=\left[\left({\displaystyle \genfrac{}{}{0pt}{}{N}{k}}\right)(uv)^k\right]^{1/2}{\displaystyle \underset{\nu =0}{\overset{k}{}}}(1)^{k\nu }\left({\displaystyle \genfrac{}{}{0pt}{}{x}{\nu }}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{Nx}{k\nu }}\right)u^{k\nu }v^\nu \left(\begin{array}{c}0kN\hfill \\ u+v=1\hfill \end{array}\right).`$ (55)
These functions form a complete system of orthogonal polynomials of integers $`x`$ with $`0xN`$. They have quite diverse applications ranging from the theory of covering codes  to the statistical mechanics of polymers . The same functions appear also as a building block in our periodic-orbit theory of Anderson localization on graphs . Unfortunately, we were not able to reduce the above expression any further by using the known sum rules and asymptotic representations for Kravtchouk polynomials. The main obstacle stems from the fact that in our case the three numbers $`N,k,x`$ in the definition (55) are constrained by $`N=k+x1`$.

We will now consider the special case $`\eta =\pi /4`$ for which we obtained in the previous subsection the solution (47). The result can be expressed in terms of Kravtchouk polynomials with $`u=v=1/2`$ which is also the most important case for the applications mentioned above. We adopt the common practice to omit the superscript $`(u,v)`$ in this special case and find

$`K_2(n;\pi /4)`$ $`=`$ $`{\displaystyle \frac{1}{2^n}}+{\displaystyle \frac{1}{2^{n+1}}}{\displaystyle \underset{q=1}{\overset{n1}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{n1}{nq}}\right)\left[{\displaystyle \frac{n}{q}}P_{n1,nq}(q)\right]^2.`$ (56)
It is convenient to introduce at this point the symbol $`𝒩(s,t)`$ defined in Eq. (2). It allows to rewrite (56) with the help of some standard transformations of binomial coefficients as

$`K_2(n;\pi /4)`$ $`=`$ $`{\displaystyle \frac{1}{2^n}}+{\displaystyle \frac{1}{2^{n+1}}}{\displaystyle \underset{q=1}{\overset{n1}{}}}\left[{\displaystyle \frac{n}{q}}𝒩(q,nq1)\right]^2`$ (57)
$`=`$ $`{\displaystyle \frac{1}{2^n}}+{\displaystyle \frac{1}{2^{n+1}}}{\displaystyle \underset{q=1}{\overset{n1}{}}}\left[𝒩(q,nq)+(1)^n𝒩(nq,q)\right]^2`$ (58)
This expression is displayed in Fig. 1 together with (47) in order to illustrate the equivalence of the two results for the form factor.

An independent proof for the equivalence of (47), (56) can be given by comparing the generating functions of $`K_2(n;\pi /4)`$ in the two representations . We define

$`G(x)`$ $`=`$ $`{\displaystyle \underset{x=1}{\overset{\mathrm{}}{}}}K_2(n;\pi /4)(2x)^n(|x|<1/2)`$ (59)
and find from (47)

$`G(x)`$ $`=`$ $`{\displaystyle \frac{2x}{12x}}{\displaystyle \frac{1}{2}}+{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(1)^m}{2}}\left({\displaystyle \genfrac{}{}{0pt}{}{2m}{m}}\right)x^{2m}(12x)`$ (60)
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{12x}{\sqrt{1+4x^2}}}{\displaystyle \frac{1}{2}}{\displaystyle \frac{16x}{12x}}.`$ (61)
On the other hand we have from (57)

$`G(x)={\displaystyle \frac{x}{1x}}+G_1(x)+G_2(x)`$ (62)
with $`G_{1,2}(x)`$ defined in the introduction. A convenient starting point to obtain the r.h.s. of (3), (4) is the integral representation

$`𝒩(s,t)={\displaystyle \frac{(1)^t}{2\pi \text{i}}}{\displaystyle dz(1+z^1)^t(1z)^{s1}},`$ (63)
where the contour encircles the origin. With the help of (63) we find

$`g(x,y)`$ $`=`$ $`{\displaystyle \underset{s,t=1}{\overset{\mathrm{}}{}}}𝒩(s,t)x^sy^t`$ (64)
$`=`$ $`{\displaystyle \frac{1}{2\pi \text{i}}}{\displaystyle \underset{s,t=1}{\overset{\mathrm{}}{}}}{\displaystyle dz\underset{s,t=1}{\overset{\mathrm{}}{}}(1+z^1)^t(1z)^{s1}x^s(y)^t}`$ (65)
$`=`$ $`{\displaystyle \frac{xy}{2\pi \text{i}}}{\displaystyle \underset{s,t=0}{\overset{\mathrm{}}{}}}{\displaystyle dz\frac{1}{1x(1z)}\frac{1+z}{z+y(1+z)}}`$ (66)
$`=`$ $`{\displaystyle \frac{xy}{(1+y)(1x+y2xy)}}(|x|,|y|<1/\sqrt{2}),`$ (67)
which was already stated in (5). The contour $`|1+z^1|=|1z|=\sqrt{2}`$ has been chosen such that both geometric series converge everywhere on it. Now we have

$`G_1(x^2)`$ $`=`$ $`{\displaystyle \frac{1}{(2\pi \text{i})^2}}{\displaystyle \frac{\mathrm{d}z\mathrm{d}z^{}}{zz^{}}\underset{s,t=1}{\overset{\mathrm{}}{}}\underset{s^{},t^{}=1}{\overset{\mathrm{}}{}}𝒩(s,t)𝒩(s^{},t^{})(xz)^s(x/z)^s^{}(xz^{})^t(x/z^{})^t^{}}`$ (68)
$`=`$ $`{\displaystyle \frac{x^4}{(2\pi \text{i})^2}}{\displaystyle dzdz^{}\frac{1}{(1+xz^{})(1+x[z^{}z]2x^2zz^{})}\frac{z^{}}{(z^{}+x)(zz^{}+x[zz^{}]2x^2)}},`$ (69)
where $`|x|<1/\sqrt{2}`$ and the contour for $`z,z^{}`$ is the unit circle. We perform the double integral using the residua inside the contour and obtain (3) and in complete analogy also (4) such that finally

$$G(x)=\frac{x}{1x}+\frac{x}{2x1}\left(\frac{1}{\sqrt{4x^2+1}}\frac{1}{1x}\right)+\frac{1}{2}\frac{4x^22x+1}{(12x)\sqrt{4x^2+1}}\frac{1}{2}.$$
(70)
The proof is completed by a straight forward verification of the equivalence between the rational functions (60) and (70).

The identities (6), (7) follow now by separating even and odd powers of $`n`$ in (47) and (56). In terms of Kravtchouk polynomials these identities can be written as

$$\underset{q=1}{\overset{2m1}{}}\left(\genfrac{}{}{0pt}{}{2m1}{2mq}\right)\left[\frac{2m}{q}P_{2m1,2mq}(q)\right]^2=2^{2m+1}+(1)^m\left(\genfrac{}{}{0pt}{}{2m}{m}\right)2$$
(71)
and

$$\underset{q=1}{\overset{2m}{}}\left(\genfrac{}{}{0pt}{}{2m}{2m+1q}\right)\left[\frac{2m+1}{q}P_{2m,2m+1q}(q)\right]^2=2^{2m+2}2(1)^m\left(\genfrac{}{}{0pt}{}{2m}{m}\right)2.$$
(72)
Finally we will derive the CUE result (50) for the ensemble of graphs defined in the previous subsection starting from the periodic-orbit expansion (51). We find

$`K_2(n)`$ $`=`$ $`{\displaystyle _0^{\pi /2}}d\mu (\eta )K_2(n;\eta ).`$ (73)
Inserting (51), expanding into a double sum and using

$$_0^{\pi /2}d\eta \mathrm{sin}^{2(\nu +\nu ^{})+1}\eta \mathrm{cos}^{2(n\nu \nu ^{})+1}\eta =\frac{1}{2(n+1)}\left(\genfrac{}{}{0pt}{}{n}{\nu +\nu ^{}}\right)^1$$
(74)
we get

$`K_2(n)`$ $`=`$ $`{\displaystyle \frac{1}{n+1}}+`$ (76)
$`+{\displaystyle \frac{n^2}{4(n+1)}}{\displaystyle \underset{q=1}{\overset{n1}{}}}{\displaystyle \underset{\nu ,\nu ^{}}{}}{\displaystyle \frac{(1)^{\nu +\nu ^{}}}{\nu \nu ^{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{n}{\nu +\nu ^{}}}\right)^1\left({\displaystyle \genfrac{}{}{0pt}{}{q1}{\nu 1}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{nq1}{\nu 1}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{q1}{\nu ^{}1}}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{nq1}{\nu ^{}1}}\right).`$
Comparing this to the equivalent result (50) we were led to the identity (1) involving a multiple sum over binomial coefficients. In this case, an independent computer-generated proof was found , which is based on the recursion relation

$$q^2F_{\nu ,\nu ^{}}(n,q)(nq1)^2F_{\nu ,\nu ^{}}(n,q+1)+(n1)(n2q1)F_{\nu ,\nu ^{}}(n+1,q+1)=0.$$
(77)
This recursion relation was obtained with the help of a Mathematica routine , but it can be checked manually in a straight forward calculation. By summing (77) over the indices $`\nu ,\nu ^{}`$, the same recursion relation is shown to be valid for $`S(n,q)`$  and the proof is completed by demonstrating the validity of (1) for a few initial values. Having proven (1) we can use it to perform the summation over $`\nu ,\nu ^{}`$ in (76) and find

$`K_2(n)={\displaystyle \frac{1}{n+1}}+{\displaystyle \underset{q=1}{\overset{n1}{}}}{\displaystyle \frac{n}{n^21}}={\displaystyle \frac{1}{n+1}}+{\displaystyle \frac{n}{n+1}}(1\delta _{n,1}),`$ (78)
which is now obviously equivalent to the random-matrix form factor (50).

### D Trace identities

Let $`S`$ be an arbitrary unitary matrix with a non degenerate spectrum and $`s_n=\mathrm{tr}S^n`$. Then

$`\underset{ϵ0}{lim}ϵ{\displaystyle \underset{n=n_0}{\overset{\mathrm{}}{}}}s_n^{}s_{n+\nu }\mathrm{e}^{nϵ}=s_\nu ,`$ (79)
for arbitrary integers $`n_0`$ and $`\nu `$ .

We shall now apply this identity to the ring graph with $`\eta =\pi /4`$ in order to prove (9). Again, the traces of $`S^n`$ can be expanded in periodic orbits which are grouped into $`n+1`$ families with equal phases $`\mathrm{\Phi }(n,q)=q\varphi _1+(nq)\varphi _2`$ ($`0qn`$). Thus, one can write

$$s_n(k)=\underset{q=0}{\overset{n}{}}A(n,q)\mathrm{e}^{\text{i}\mathrm{\Phi }(n,q)},$$
(80)
where $`A(n,q)`$ is the coherent sum of all amplitudes of PO’s in the corresponding set. $`A(n,q)`$ can be expressed in terms of Kravtchouk Polynomials as

$$A(n,q)=\frac{1}{\sqrt{2^n}}\{\begin{array}{cc}1\hfill & \mathrm{for}q=0\mathrm{or}n\hfill \\ (1)^{n+q}(n/q)\left(\genfrac{}{}{0pt}{}{n1}{nq}\right)^{1/2}P_{n1,nq}(q)\hfill & \mathrm{for}0<q<n\hfill \end{array}$$
(81)
(compare Eq. (52)). This is equivalent to (8). Substituting (80) into the trace identity (79), we get for arbitrary integers $`\nu `$ and $`n_0`$

$`\underset{ϵ0}{lim}ϵ`$ $`{\displaystyle \underset{n=n_0}{\overset{\mathrm{}}{}}}`$ $`\mathrm{e}^{nϵ}{\displaystyle \underset{q=0}{\overset{n}{}}}{\displaystyle \underset{p=0}{\overset{n+\nu }{}}}A(n+\nu ,p)A(n,q)\mathrm{e}^{\text{i}\left[(pq)\varphi _1+(\nu (pq))\varphi _2\right]}`$ (82)
$`=`$ $`{\displaystyle \underset{\kappa =0}{\overset{\nu }{}}}A(\nu ,\kappa )\mathrm{e}^{\text{i}\left[\kappa \varphi _1+(\nu \kappa )\varphi _2\right]}.`$ (83)
This is valid for arbitrary phases $`\varphi _1`$ and $`\varphi _2`$ and therefore the coefficients of the phase factors $`\mathrm{e}^{\text{i}\mathrm{\Phi }(\nu ,\kappa )}`$ on both sides are equal. Eq. (9) follows.

## IV Conclusions

We have shown how within periodic-orbit theory the problem of finding the form factor (the spectral two-point correlation function) for a quantum graph can be reduced exactly to a well-defined combinatorial problem. Even for the very simple graph model that we considered in the last section the combinatorial problems involved were highly non-trivial. In fact we encountered previously unknown identities which we could not have obtained if it were not for the independent method of computing the form factor directly from the spectrum. However, since the pioneering work documented in  the investigation of sums of the type we encountered in this paper is a rapidly developing subject, and it can be expected that finding identities like (1), (6) and (7) will shortly be a matter of computer power.

Numerical simulations in which the form factor was computed for fully connected graphs ($`v_i=V1i`$) and, e. g., Neumann boundary conditions indicate that the spectral correlations of $`S`$ match very well Dyson’s predictions for the circular ensembles COE or CUE, respectively . The agreement is improved as $`V`$ increases. We conjecture that this might be rigorously substantiated by asymptotic combinatorial theory. A first step towards this goal was taken in the present paper.

## V Acknowledgements

This research was supported by the Minerva Center for Physics of Nonlinear Systems, and by a grant from the Israel Science Foundation. We thank Uri Gavish for introducing us to the combinatorial literature, and Brendan McKay and Herbert Wilf for their interest and support. We are indebted to Gregory Berkolaiko for his idea concerning the proof of (6) and (7), and to Akalu Tefera for his kind help in obtaining a computer-aided proof of (1).


format: text

text:
# Monte Carlo simulation of nonlinear Couette flow in a dilute gas

## I Introduction

One of the most interesting states for analyzing transport processes far from equilibrium is the steady planar Couette flow. The physical situation corresponds to a system enclosed between two infinite, parallel plates in relative motion and, in general, kept at different temperatures. These boundary conditions lead to combined heat and momentum transport. If $`x`$ and $`y`$ denote the coordinates parallel to the flow and orthogonal to the plates, respectively, then the corresponding steady hydrodynamic balance equations are

$$\frac{P_{xy}}{y}=\frac{P_{yy}}{y}=0,$$
(1)
$$P_{xy}\frac{u_x}{y}+\frac{q_y}{y}=0,$$
(2)
where $`𝐮=u_x\widehat{𝐱}`$ is the flow velocity, $`𝖯`$ is the pressure tensor, and $`𝐪=q_x\widehat{𝐱}+q_y\widehat{𝐲}`$ is the heat flux. The presence of $`q_y`$ in Eq. (2) indicates that a thermal gradient $`T/y`$ is induced by the velocity gradient, even if both plates are kept at the same temperature. The balance equations (1) and (2) do not constitute a closed set unless the dependence of the pressure tensor and the heat flux on the hydrodynamic fields is known. If the gradients are small, the fluxes $`𝖯`$ and $`𝐪`$ are described by the Navier-Stokes (NS) constitutive relations, which in this problem yield

$$P_{xx}=P_{yy}=P_{zz},P_{xy}=\eta _0\frac{u_x}{y},$$
(3)
$$q_x=0,q_y=\kappa _0\frac{T}{y},$$
(4)
where $`\eta _0`$ and $`\kappa _0`$ are the NS shear viscosity and thermal conductivity coefficients, respectively. As a consequence of the absence of normal stress differences in the NS description, the hydrostatic pressure $`p=(P_{xx}+P_{yy}+P_{zz})/3`$ is a constant, on account of the balance equation (1).

Even in the linear regime described by the NS equations, one still needs to know the spatial dependence of the transport coefficients to obtain the exact solution of the hydrodynamic equations. The problem becomes tractable in the case of a low density gas, where the state of the system is completely specified by the velocity distribution function $`f(𝐫,𝐯;t)`$, which obeys the Boltzmann equation.  A relevant dimensionless quantity in a dilute gas is the Knudsen number $`\text{Kn}=\lambda /\mathrm{}_h`$, defined as the ratio of the mean free path $`\lambda `$ to the scale length of the hydrodynamic gradients $`\mathrm{}_h`$. In many laboratory conditions, $`\text{Kn}1`$ and so the Boltzmann equation can be solved by means of the Chapman-Enskog method as an expansion of the distribution function in powers of the Knudsen number. The zeroth order approximation leads to the Euler hydrodynamic equations, while the first order approximation yields the NS equations with explicit expressions for the transport coefficients $`\eta _0`$ and $`\kappa _0`$. The results show that the ratio $`\eta _0/\kappa _0`$ is a constant. Consequently, it then follows from the NS hydrodynamic equations (1)–(4) that the flow velocity profile is quasi-linear,

$$\eta _0\frac{u_x}{y}=\text{const},$$
(5)
and the temperature is quasi-parabolic,

$$\left(\kappa _0\frac{}{y}\right)^2T=\frac{\kappa _0}{\eta _0}\left(\eta _0\frac{u_x}{y}\right)^2=\text{const}.$$
(6)
Note that the profile of $`u_x`$ is not strictly linear, due to the space dependence of $`\eta _0`$ through the temperature. Analogously, the temperature profile is not strictly quadratic. In fact, the specific form of both profiles depends on the interaction potential under consideration. On the other hand, from Eqs. (5) and (6) it is easy to derive a nice result, namely that if the temperature $`T`$ is seen as a function of $`u_x`$ rather than as a function of the coordinate space $`y`$, then one has

$$\frac{^2T}{u_x^2}=\frac{\eta _0}{\kappa _0}.$$
(7)
This is a sort of nonequilibrium “equation of state,” according to which the temperature is a quadratic function of the flow velocity. Moreover, the “curvature” of the profile is practically universal, given the weak influence of the interaction potential on the Prandtl number $`\text{Pr}=5k_B\eta _0/2m\kappa _0\frac{2}{3}`$, where $`k_B`$ is the Boltzmann constant and $`m`$ is the mass of a particle.

Since the mean free path is inversely proportional to the density, the Knudsen number, at a given value of the scale length $`\mathrm{}_h`$, increases as the gas becomes more rarefied. In general, when the Knudsen number is not small, the NS relations are not expected to hold and the transport must be described by nonlinear constitutive equations. In the special case of Maxwell molecules (particles interacting via an $`r^4`$ potential), it has been shown  that the Boltzmann equation admits a consistent solution in the nonlinear Couette flow characterized by a constant pressure $`p`$ and profiles similar to those obtained in the NS regime, Eqs. (5)–(7), except that $`\eta _0`$ and $`\kappa _0`$ are replaced by a generalized shear viscosity coefficient $`\eta (a)=\eta _0F_\eta (a)`$ and a generalized thermal conductivity coefficient $`\kappa _{yy}(a)=\kappa _0F_\kappa (a)`$, respectively. Here, $`a=(\eta _0/p)u_x/y`$ is a constant (dimensionless) shear rate and $`F_\eta `$ and $`F_\kappa `$ are nonlinear functions of $`a`$. In addition, $`P_{xx}P_{yy}P_{zz}`$ and $`q_x0`$. In this problem, the hydrodynamic scale length can be identified as $`\mathrm{}_h\sqrt{k_BT/m}(u_x/y)^1`$, while the mean free path is $`\lambda \sqrt{k_BT/m}(\eta _0/p)`$. Thus, the reduced shear rate $`a`$ represents the Knudsen number in this problem, i.e. $`a\text{Kn}`$. Henceforth, we will use the reduced shear rate $`a`$ to refer to the Knudsen number Kn. The solution considered in Refs.  describes heat and momentum transport for arbitrary velocity and thermal gradients in the bulk domain, where boundary effects are negligible. On the other hand, the full nonlinear dependence of $`F_\eta (a)`$ and $`F_\kappa (a)`$ is not explicitly known, since it involves the infinite hierarchy of moment equations. Their knowledge is limited to super-Burnett order and the result is  $`F_\eta (a)=13.111a^2`$ and $`F_\kappa (a)=17.259a^2`$.

Consequently, if one wants to get the transport properties for arbitrary values of the shear rate and the thermal gradient, one must resort to approximate schemes or to computer simulations. In the first alternative, explicit expressions for the nonlinear transport coefficients in the Couette flow have been obtained from exact solutions of the Bhatnagar-Gross-Krook (BGK) model  and related models  for general interactions, as well as from the Grad method.  In the simulation side, Risso and Cordero  have recently studied the shear-rate dependence of $`F_\eta `$ and $`F_\kappa `$ by means of molecular dynamics simulations of a hard disk gas. Comparison between the different analytical results with those obtained from the simulation shows that the predictions given by kinetic models are in better agreement than those given by the Grad method, especially in the case of the thermal conductivity.  Nevertheless, given the difficulties associated with molecular dynamics simulations to achieve large shear rates in a dilute gas, the above comparison is restricted to a range of shear rates for which non-Newtonian effects are hardly significant. For instance, the shear viscosity has only decreased around $`10\%`$ with respect to its Navier-Stokes value for the largest value of the shear rate considered by Risso and Cordero.  In order to overcome such difficulties and extend the range of values of $`a`$, one may use the so-called Direct Simulation Monte Carlo (DSMC) method,  which is known to qualify as an efficient and accurate method to numerically solve the Boltzmann equation.

The aim of this paper is to solve the Boltzmann equation by means of the DSMC method for a gas subjected to the planar Couette flow. The motivation of this work is twofold. On the one hand, we want to test the reliability of the far from equilibrium results obtained from kinetic models and the Grad method by making a comparison with the Boltzmann solution in the case of Maxwell molecules, for which the form of the hydrodynamic profiles in the bulk region (far from the boundaries) is known. We will determine not only the hydrodynamic profiles but also the nonlinear transport coefficients and the velocity distribution function. On the other hand, we want to investigate whether the above results for a system of Maxwell molecules extend to other interaction potentials. This extension holds when the Boltzmann equation is replaced by kinetic models where, in terms of a conveniently scaled space variable, all the results are independent of the interaction law. Thus, we will also solve numerically the Boltzmann equation by the DSMC method for a hard-sphere gas.

Since we are interested in describing transport properties in the bulk region, i.e., free from finite-size effects, we need to use appropriate boundary conditions in the simulations. In the conventional boundary conditions,  the gas is assumed to be enclosed between two baths at equilibrium in relative motion and, in general, at different temperatures. Under these conditions, a particle leaving the system is formally replaced by a particle coming from the bath, so the incoming velocity is sampled from a local equilibrium distribution. As a consequence, there exists a mismatch between the velocity distribution function of the reemitted particles and that of those particles of the gas located near the wall and moving along the same direction. In this case, in order to inhibit the influence of boundary effects one needs to take very large systems (normal distance between the plates much larger than the mean free path), what is not practical from a computational point of view. To overcome this difficulty, a possibility is to assume that both baths are out of equilibrium in a state close to that of the actual gas. Since such a state is not known “a priori”, in this paper we assume that the state of the baths is described by the BGK solution of the planar Couette flow.  Although the boundary effects are still unavoidable, we expect that the above mismatch between reemitted and gas particles will be much smaller. As a matter of fact, the use of these conditions has been shown to be appropriate to analyze bulk transport properties in the special case of planar Fourier flow (both walls at rest).

The plan of the paper is as follows. In Sec. II we give a brief description of the planar Couette flow and a summary of the main results obtained from the Boltzmann equation and kinetic models. The boundary conditions used in the simulations and the DSMC method are described in Sec. III. Section IV presents the main results of the paper, where special attention is paid to the nonlinear transport coefficients. A comparison with the analytical results derived from the BGK and the ellipsoidal statistical (ES) models and from the Grad method is also carried out. The comparison shows in general a good agreement of the kinetic models with computer simulations, even for large shear rates. In addition, the velocity distribution function obtained from the simulation in the bulk domain is compared with the ones given by the kinetic models. It is shown again that the agreement is qualitatively good. We close the paper in Sec. V with some concluding remarks.

## II Description of the problem and summary of theoretical results

Let us consider a dilute gas. In this case, a kinetic description is sufficient to characterize the state of the system by means of the velocity distribution function $`f(𝐫,𝐯;t)`$. This distribution function obeys the nonlinear Boltzmann equation, which in the absence of external forces is given by

$$\frac{f}{t}+𝐯f=J[f,f],$$
(8)
where

$$J[f,f]=𝑑𝐯_1𝑑\widehat{𝐤}gI(g,\widehat{𝐤})\left[f(𝐯^{})f(𝐯_1^{})f(𝐯)f(𝐯_1)\right]$$
(9)
is the collision operator. In this equation, $`I(g,\widehat{𝐤})`$ is the differential cross section, $`g|𝐯𝐯_1|`$ being the relative velocity, and $`(𝐯^{},𝐯_1^{})`$ are precollisional velocities yielding postcollisional velocities $`(𝐯,𝐯_1)`$. From the distribution function, one may define the hydrodynamic quantities

$$n=𝑑𝐯f,$$
(10)
$$𝐮=\frac{1}{n}𝑑𝐯𝐯f,$$
(11)
$$\frac{3}{2}nk_BT=\frac{m}{2}𝑑𝐯𝐕^2f,$$
(12)
and the momentum and heat fluxes

$$𝖯=m𝑑𝐯\mathrm{𝐕𝐕}f,$$
(13)
$$𝐪=\frac{m}{2}𝑑𝐯𝐕^2𝐕f.$$
(14)
Here, $`n`$ is the number density, $`𝐮`$ is the flow velocity, $`T`$ is the temperature, $`𝖯`$ is the pressure tensor, $`𝐪`$ is the heat flux, and $`𝐕=𝐯𝐮`$ is the peculiar velocity. In addition, the equation of state is that of an ideal gas, i.e., $`p=nk_BT`$.

Most of the known solutions to Eq. (8) for spatially inhomogeneous states  correspond to the special case of Maxwell molecules, namely, a repulsive potential of the form $`V(r)r^4`$. For this potential, the collision rate $`gI(g,\widehat{𝐤})`$ is independent of the relative velocity and this allows the infinite hierarchy of velocity moments to be recursively solved in some specific situations. Furthermore, the NS transport coefficients $`\eta _0`$ and $`\kappa _0`$ can be exactly obtained from the Chapman-Enskog method.  They are given by

$$\eta _0=\frac{p}{\nu },\kappa _0=\frac{15}{4}\frac{k_B}{m}\eta _0,$$
(15)
where $`\nu =\theta n`$, $`\theta `$ being an eigenvalue of the linearized Boltzmann collision operator.

In this paper we are interested in studying the planar Couette flow for a dilute gas. We consider a gas enclosed between two parallel plates in relative motion and maintained at different temperatures. Under these conditions, the system is driven out of equilibrium by the combined action of the velocity and thermal gradients along the direction normal to the plates. After a transient period, the gas is expected to reach a steady state and the corresponding Boltzmann equation reads

$$v_y\frac{}{y}f=J[f,f],$$
(16)
where we have chosen the axis $`y`$ as the one normal to the plates. In general, this equation must be solved subjected to specific boundary conditions. Nevertheless, in the same spirit as in the Chapman-Enskog method,  one may look for “normal” solutions in which all the space dependence of the distribution function occurs through a functional dependence on the hydrodynamic fields. The normal solution describes the state of the gas in the hydrodynamic regime, namely, for times much longer than the mean free time and for distances from the walls much larger than the mean free path. As mentioned in the Introduction, it has been proved  that, in the case of Maxwell molecules, Eq. (16) admits an exact solution corresponding to the planar Couette flow problem. This solution belongs to the normal class, so that no explicit boundary conditions appear. In other words, all the space dependence of $`f`$ is given through the local density, the local velocity, the local temperature, and their gradients. The solution is characterized by hydrodynamic profiles that are a simple generalization of those predicted by the NS approximation, Eqs. (5)–(7), namely

$$p=\text{const},$$
(17)
$$\frac{1}{\nu (y)}\frac{u_x}{y}a=\text{const},$$
(18)
$$\left[\frac{1}{\nu (y)}\frac{}{y}\right]^2T=\text{Pr}\frac{2m}{k_B}\gamma (a).$$
(19)
The dimensionless parameter $`\gamma (a)`$ is a nonlinear function of the reduced shear rate $`a`$ that, by construction, behaves as $`\gamma a^2/5`$ in the limit $`a0`$. Again, the temperature can be seen as a quadratic function of the flow velocity, but now the coefficient $`\eta _0/\kappa _0`$ appearing in Eq. (7) is replaced by a shear-rate dependent coefficient $`(\eta _0/\kappa _0)5\gamma (a)/a^2`$. Furthermore, in this solution the pressure tensor is independent of the thermal gradient and the heat flux vector is linear in the thermal gradient, but all these fluxes are nonlinear functions of the shear rate. This nonlinear dependence can be characterized through five generalized transport coefficients. First, the shear stress $`P_{xy}`$ defines a generalized shear viscosity $`\eta (a)`$ as

$$P_{xy}=\eta (a)\frac{u_x}{y}\eta _0F_\eta (a)\frac{u_x}{y}.$$
(20)
Analogously, the component of the heat flux parallel to the thermal gradient defines a generalized thermal conductivity coefficient $`\kappa _{yy}(a)`$:

$$q_y=\kappa _{yy}(a)\frac{T}{y}\kappa _0F_\kappa (a)\frac{T}{y}.$$
(21)
The dimensionless functions $`F_\eta (a)`$ and $`F_\kappa (a)`$ are the most relevant quantities of the problem. They are related to the function $`\gamma (a)`$ by $`\gamma (a)=a^2F_\eta (a)/5F_\kappa (a)`$. Normal stress differences are different from zero and are measured by two viscometric functions $`\mathrm{\Psi }_{1,2}(a)`$ defined by

$$\frac{P_{yy}P_{xx}}{p}=\mathrm{\Psi }_1(a)a^2,$$
(22)
$$\frac{P_{zz}P_{yy}}{p}=\mathrm{\Psi }_2(a)a^2.$$
(23)
Another interesting quantity related to the pressure tensor is the friction coefficient $`\mu (a)=P_{xy}/P_{yy}`$. To NS order, we simply have $`\mu (a)=a`$. In the non-Newtonian regime, we generalize this coefficient as $`\mu (a)=aF_\mu (a)`$, where the friction function $`F_\mu (a)`$ is

$$F_\mu (a)=\frac{F_\eta (a)}{1[\mathrm{\Psi }_2(a)\mathrm{\Psi }_1(a)]a^2/3}.$$
(24)
Finally, there exists a nonzero component of the heat flux orthogonal to the thermal gradient given by

$$q_x=\kappa _{xy}(a)\frac{}{y}T\kappa _0\mathrm{\Phi }(a)a\frac{}{y}T.$$
(25)
The three functions $`\mathrm{\Psi }_{1,2}(a)`$ and $`\mathrm{\Phi }(a)`$ are generalizations of Burnett coefficients. In fact, $`\mathrm{\Psi }_1(0)=14/5`$, $`\mathrm{\Psi }_2(0)=4/5`$, and $`\mathrm{\Phi }(0)=7/2`$ for Maxwell molecules and for hard spheres in the first Sonine approximation. The determination of the nonlinear transport coefficients $`F_\eta `$, $`F_\kappa `$, $`\mathrm{\Psi }_{1,2}`$, and $`\mathrm{\Phi }`$ would imply the solution of an infinite hierarchy that cannot be solved in a recursive way.  This hierarchy can only be solved step by step when one performs a perturbation expansion in powers of the shear rate. In fact, Tij and Santos  determined the solution up to super-Burnett order:

$$F_\eta (a)=13.111a^2+𝒪(a^4),$$
(26)
$$F_\kappa (a)=17.259a^2+𝒪(a^4).$$
(27)
From Eqs. (26) and (27) it follows that $`\gamma (a)=(a^2/5)[1+4.148a^2+𝒪(a^4)]`$. In addition, $`F_\mu (a)=11.911a^2+𝒪(a^3)`$.

Although the above analyses are valuable, they have two main limitations. On the one hand, they are restricted to the special case of Maxwell molecules. For other interaction potentials (e.g., hard spheres), the collisional moments involve all the moments of the distribution function and, as a consequence, the hydrodynamic profiles (17)–(19) are not strictly true. On the other hand, even for Maxwell molecules, the above perturbative solution is not useful for finite shear rates. These two limitations can be overcome, in the context of analytical methods, by introducing additional approximations, such as the Grad method,  or by describing the system by means of kinetic models.  In these approaches, one looks for a solution having the same hydrodynamic profiles as in the case of the Boltzmann equation, cf. Eqs. (17)–(19). As before, this solution describes the properties of the system in the bulk region, which is insensitive to the details of the boundary conditions. From the Grad method and from the BGK and ES kinetic models, one explicitly gets the full nonlinear shear-rate dependence of the transport coefficients in a non-perturbative way. Moreover, the results are universal, namely the functions $`F_\eta (a)`$, $`F_\kappa (a)`$, …are independent of the interaction potential, provided that the reduced shear rate is defined as in Eq. (18) with $`\nu =p/\eta _0`$.

The thirteen-moment Grad method consists of replacing the actual distribution by

$$ff_L\left\{1+\frac{m}{n(k_BT)^2}\left[\left(\frac{mV^2}{5k_BT}1\right)𝐕𝐪+\frac{1}{2}\left(P_{ij}p\delta _{ij}\right)V_iV_j\right]\right\},$$
(28)
where

$$f_L=n\left(\frac{m}{2\pi k_BT}\right)^{3/2}\mathrm{exp}\left(\frac{mV^2}{2k_BT}\right)$$
(29)
is the local equilibrium distribution function. When the approximation (28) is inserted into the Boltzmann equation (8) and velocity moments are taken, one gets a closed set of equations for $`n`$, $`𝐮`$, $`𝖯`$, and $`𝐪`$. According to the geometry of the planar Couette flow, there are eight independent moments instead of the original thirteen moments appearing in Eq. (28). Risso and Cordero  found that the set of independent moment equations, neglecting nonlinear terms in the fluxes, admits a solution consistent with the profiles (17)–(19). In addition, they obtained explicit expressions for the transport coefficients introduced above as nonlinear functions of the reduced shear rate $`a`$. Those expressions are displayed in the Appendix.

Now we consider the kinetic model approach. The basic idea is to replace the detailed Boltzmann collision operator by a much simpler term that otherwise retains the main physical properties of the original operator $`J[f,f]`$. The most familiar choice is the BGK model,

$$J[f,f]\nu (ff_L),$$
(30)
where $`\nu `$ is an effective collision frequency that depends on the temperature, according to the interaction potential. The NS transport coefficients of the BGK model are $`\eta _0=p/\nu `$ and $`\kappa _0=5pk_B/2m\nu `$. Thus the BGK model has the drawback that it predicts an incorrect value for the Prandtl number, $`\text{Pr}=1`$. This is a consequence of the fact that $`\nu `$ is the only adjustable parameter in the model. This deficiency is corrected by the so-called ellipsoidal statistical (ES) model,  in which case

$$J[f,f]\zeta (ff_R),$$
(31)
where $`\zeta `$ is again an effective collision frequency and the reference function $`f_R`$ is

$$f_R=n\pi ^{3/2}(det𝖠)^{1/2}\mathrm{exp}\left(𝖠^1:\mathrm{𝐕𝐕}\right),$$
(32)
where $`A_{ij}=(2k_BT/m)\text{Pr}^1\delta _{ij}2(\text{Pr}^11)P_{ij}/mn`$. The NS coefficients are now $`\eta _0=p/(\zeta \text{Pr}^1)`$ and $`\kappa _0=5pk_B/2m\zeta `$. If, as in Eq. (15), we define a collision frequency $`\nu =p/\eta _0`$, then $`\nu =\zeta \text{Pr}^1`$ in the ES model. Note that the ES model reduces to the BGK model if we formally make $`\text{Pr}=1`$. Therefore, the ES model can be seen as an extension of the BGK model to account for the correct Prandtl number, which can be particularly relevant in the Couette flow where heat transport and momentum transport are coupled. The results derived from the BGK and the ES models for the Couette flow problem are also given in the Appendix. Apart from obtaining the transport properties, the velocity distribution function can be explicitly written. In particular, the BGK distribution is given by

$`f(𝐫,𝐯)`$ $`=`$ $`n\left({\displaystyle \frac{m}{2\pi k_BT}}\right)^{3/2}{\displaystyle \frac{2\alpha (1+\alpha )^{3/2}}{ϵ|\xi _y|}}{\displaystyle _{t_0}^{t_1}}𝑑t\left[2t(1\alpha )t^2\right]^{5/2}`$ (34)
$`\times \mathrm{exp}\left\{{\displaystyle \frac{2\alpha }{1+\alpha }}{\displaystyle \frac{1t}{ϵ\xi _y}}{\displaystyle \frac{1+\alpha }{2t(1\alpha )t^2}}\left[\left(\xi _x+{\displaystyle \frac{2a\alpha }{1+\alpha }}{\displaystyle \frac{1t}{ϵ}}\right)^2+\xi _y^2+\xi _z^2\right]\right\}.`$
Here, $`(t_0,t_1)=(0,1)`$ if $`\xi _y>0`$ and $`(t_0,t_1)=[1,2/(1\alpha )]`$ if $`\xi _y<0`$. Besides, $`𝝃(m/2k_BT)^{1/2}𝐕`$,

$$\alpha =\frac{ϵ}{\left(ϵ^2+8\gamma \right)^{1/2}},$$
(35)
and

$$ϵ=\left(\frac{2k_BT}{m}\right)^{1/2}\frac{1}{\nu }\frac{}{y}\mathrm{ln}T$$
(36)
is a (local) reduced thermal gradient. Equation (34) shows that the distribution function is a highly nonlinear function of the reduced gradients $`a`$ and $`ϵ`$.

In Ref.  a comparison between the analytical results obtained from the Grad method and the BGK and ES models with those obtained from molecular dynamics simulations  for hard disks was carried out. It was found that the three theories reproduced quite well the simulation data for $`F_\eta `$, but the Grad method failed for $`F_\kappa `$ and $`\mathrm{\Phi }`$. The latter quantity was reproduced better by the ES model than by the BGK model. Notwithstanding this, more definite conclusions could not be reached because the simulation data were restricted to rather small shear rates, namely $`a0.2`$. In this range of shear rates, non-Newtonian effects are not especially significant. In addition, although the molecular dynamics simulations correspond to a very dilute gas (area fraction $`\varphi 1\%`$), the collisional contributions to the transport coefficients (absent in a Boltzmann description) are not strictly zero. Finally, as a minor point, conclusions drawn in the context of two-dimensional systems should not be extrapolated without caution to the more realistic case of three dimensions. As said in the Introduction, the aim of this paper is to solve numerically the Boltzmann equation by means of the DSMC method for the planar Couette flow and compare the results with the theoretical predictions. We will consider three-dimensional systems of Maxwell molecules and hard spheres subjected to shear rates as large as $`a1.2`$. In addition to the nonlinear transport coefficients, the comparison will be extended to the level of the velocity distribution function itself.

## III Boundary conditions and Monte Carlo simulation

### A Boundary conditions

The goal now is to solve numerically the Boltzmann equation corresponding to the planar Couette flow by using the successful DSMC method.  The gas is enclosed between two parallel plates located at $`y=0`$ and $`y=L`$, which are moving along the $`x`$-direction with velocities $`𝐔_0=U_0\widehat{𝐱}`$ and $`𝐔_L=U_L\widehat{𝐱}`$, respectively. In addition, they are kept at temperatures $`T_0`$ and $`T_L`$, respectively. In order to solve Eq. (16), we need to impose the corresponding boundary conditions. They can be expressed in terms of the kernels $`K_{0,L}(𝐯,𝐯^{})`$ defined as follows. When a particle with velocity $`𝐯^{}`$ hits the wall at $`y=L`$, the probability of being reemitted with a velocity $`𝐯`$ within the range $`d𝐯`$ is $`K_L(𝐯,𝐯^{})d𝐯`$; the kernel $`K_0(𝐯,𝐯^{})`$ represents the same but at $`y=0`$. The boundary conditions are then

$`\mathrm{\Theta }(\pm v_y)|v_y|f(y=\{0,L\},𝐯)`$ $`=`$ $`\mathrm{\Theta }(\pm v_y){\displaystyle 𝑑𝐯^{}|v_y^{}|K_{0,L}(𝐯,𝐯^{})}`$ (38)
$`\times \mathrm{\Theta }(v_y^{})f(y=\{0,L\},𝐯^{},t).`$
In this paper we consider boundary conditions of complete accommodation with the walls, so that $`K_{0,L}(𝐯,𝐯^{})=K_{0,L}(𝐯)`$ does not depend on the incoming velocity $`𝐯^{}`$ and can be written as

$$K_{0,L}(𝐯)=A_{0,L}^1\mathrm{\Theta }(\pm v_y)|v_y|\varphi _{0,L}(𝐯),A_{0,L}=𝑑𝐯\mathrm{\Theta }(\pm v_y)|v_y|\varphi _{0,L}(𝐯),$$
(39)
where $`\varphi _{0,L}(𝐯)`$ represents the probability distribution of a fictitious gas in contact with the system at $`y=\{0,L\}`$. Equation (39) can then be interpreted as meaning that when a particle hits a wall, it is absorbed and then replaced by a particle leaving the fictitious bath. Of course, any choice of $`\varphi _{0,L}(𝐯)`$ must be consistent with the imposed wall velocities and temperatures, i.e.,

$$U_{0,L}=𝑑𝐯v_x\varphi _{0,L}(𝐯),$$
(40)
$$k_BT_{0,L}=\frac{1}{3}m𝑑𝐯(𝐯𝐔_{0,L})^2\varphi _{0,L}(𝐯).$$
(41)
The simplest and most common choice is that of a Maxwell-Boltzmann (MB) distribution:

$$\varphi _{0,L}^{\text{MB}}(𝐯)=\left(\frac{m}{2\pi k_BT_{0,L}}\right)^{3/2}\mathrm{exp}\left[\frac{m(𝐯𝐔_{0,L})^2}{2k_BT_{0,L}}\right].$$
(42)
Under these conditions, the system is understood to be enclosed between two independent baths at equilibrium. While the conditions (42) are adequate for analyzing boundary effects,  they are not very efficient when one is interested in obtaining the transport properties in the bulk region. In order to inhibit the influence of boundary effects, it is more convenient to imagine that the two fictitious baths are in nonequilibrium states resembling the state of the actual gas near the walls. More specifically, we can assume that the fictitious gases are described by the BGK equation, whose exact solution for the steady planar Couette flow is given by Eq. (34). In this case, the probability distributions $`\varphi _{0,L}(𝐯)`$ are

$`\varphi _{0,L}^{\text{BGK}}(𝐯)`$ $`=`$ $`\pi ^{3/2}{\displaystyle \frac{m}{k_BT_{0,L}}}{\displaystyle \frac{\alpha _{0,L}(1+\alpha _{0,L})^{3/2}}{ϵ_{0,L}|v_y|}}{\displaystyle _{t_0}^{t_1}}𝑑t\left[2t(1\alpha _{0,L})t^2\right]^{5/2}`$ (45)
$`\times \mathrm{exp}\{\left({\displaystyle \frac{2k_BT_{0,L}}{m}}\right)^{1/2}{\displaystyle \frac{2\alpha _{0,L}}{1+\alpha _{0,L}}}{\displaystyle \frac{1t}{ϵ_{0,L}v_y}}{\displaystyle \frac{m}{2k_BT_{0,L}}}{\displaystyle \frac{1+\alpha _{0,L}}{2t(1\alpha _{0,L})t^2}}`$
$`\times [(v_xU_{0,L}+{\displaystyle \frac{2a^{}\alpha _{0,L}}{1+\alpha _{0,L}}}{\displaystyle \frac{1t}{ϵ_{0,L}}})^2+v_y^2+v_z^2]\},`$
where $`(t_0,t_1)=(0,1)`$ if $`v_y>0`$ and $`(t_0,t_1)=[1,2/(1\alpha _{0,L})]`$ if $`v_y<0`$, and $`\alpha _{0,L}=ϵ_{0,L}/[ϵ_{0,L}^2+8\gamma _{\text{BGK}}(a^{})]^{1/2}`$. Here, $`a^{}`$ is the estimated value of the reduced shear rate, as predicted by the BGK model for specific values of the boundary parameters $`U_{0,L}`$ and $`T_{0,L}`$, and $`\gamma _{\text{BGK}}(a^{})`$ is obtained from Eq. (68). Given the values of the four independent boundary parameters $`U_{0,L}`$ and $`T_{0,L}`$ (as well as the distance $`L`$), the shear rate $`a^{}`$ and the local thermal gradients $`ϵ_{0,L}`$ are fixed by the conditions (17)–(19). Therefore,

$$a^{}=\frac{U_LU_0}{},$$
(46)
$$ϵ_{0,L}=\frac{1}{}\left(\frac{2k_BT_{0,L}}{m}\right)^{1/2}\left[\frac{T_LT_0}{T_{0,L}}\pm \frac{m\gamma _{\text{BGK}}(a^{})\text{Pr}}{k_BT_{0,L}}^2\right].$$
(47)
In these equations, $``$ is related to the actual separation $`L`$ between the plates through the nonlinear equation

$$L=_0^{}\frac{ds}{\nu (s)},$$
(48)
where $`s`$ is a variable in terms of which the temperature is a quadratic function, namely $`T(s)=T_0[1+ϵ_0(m/2k_BT_0)^{1/2}sm\gamma _{\text{BGK}}(a^{})\text{Pr}s^2/k_BT_0]`$, and the $`s`$-dependence of the collision frequency $`\nu (s)`$ appears only through the temperature (taking into account that $`p=\text{const}`$). The solution of the nonlinear set of equations (46)–(48) gives $`a^{}`$, $`ϵ_{0,L}`$, and $``$ for any choice of $`U_{0,L}`$, $`T_{0,L}`$, and $`L`$. Nevertheless, from a practical point of view, it is more convenient to fix $`U_0`$, $`T_{0,L}`$, $`a^{}`$, and $`ϵ_0`$ as independent parameters. Without loss of generality we take $`U_0=0`$ and $`T_L=1`$. In addition, we will choose $`ϵ_0=0`$. This implies that, if boundary effects were absent, the temperature would have a maximum at the lower plate $`y=0`$. Thus,

$$U_L=a^{},=\left[\frac{\mathrm{\Delta }}{2\text{Pr}\gamma _{\text{BGK}}(a^{})}\right]^{1/2},$$
(49)
$$ϵ_L=2\frac{\mathrm{\Delta }}{}.$$
(50)
In the above equations $`\mathrm{\Delta }T_01`$ and we have taken $`m=1`$ and $`k_B=\frac{1}{2}`$. Finally, the actual distance $`L`$ is given by Eq. (48). However, from the simulation point of view, it is more useful to express $`L`$ in terms of the collision frequency $`\overline{\nu }`$ corresponding to the temperature $`T_L`$ and the average density $`\overline{n}`$. In other words, instead of Eq. (48) we use

$$\overline{n}=\frac{1}{L}_0^L𝑑yn(y)=\frac{1}{L}_0^{}𝑑s\frac{n(s)}{\nu (s)}$$
(51)
to determine $`L`$. For the sake of concreteness, let us consider repulsive potentials, for which $`\nu =\overline{\nu }(n/\overline{n})(T/T_L)^\omega `$, where $`\omega `$ ranges from 0 (Maxwell molecules) to $`\frac{1}{2}`$ (hard spheres). In that case,

$$L=\frac{1}{\overline{\nu }}_0^{}𝑑s[T(s)]^\omega .$$
(52)
For Maxwell molecules, this simply reduces to $`L=/\overline{\nu }`$, while for hard spheres one has $`L=(/\overline{\nu })\mathrm{tan}^1(\sqrt{\mathrm{\Delta }})/\sqrt{\mathrm{\Delta }}`$. In summary, given the values of $`a^{}`$ and $`\mathrm{\Delta }`$, the separation $`L`$ between the plates and the velocity of the upper plate $`U_L`$ are uniquely determined. In addition, the thermal gradient at the upper plate $`ϵ_L`$ is also determined, while the value at the lower plate is fixed as $`ϵ_0=0`$. The knowledge of these boundary parameters allows one to obtain the distributions $`\varphi _{0,L}^{\text{BGK}}(𝐯)`$, according to Eq. (45).

### B The DSMC method

Now, we briefly describe the numerical algorithm we have employed to solve the Boltzmann equation by means of the so-called Direct Simulation Monte Carlo (DSMC) method.  In this method, the velocity distribution function is represented by the velocities $`\{𝐯_i\}`$ and positions $`\{y_i\}`$ of a sufficiently large number of particles $`N`$. Given the geometry of the problem, the physical system is split into layers of width $`\delta y`$, sufficiently smaller than the mean free path. The velocities and coordinates are updated from time $`t`$ to time $`t+\delta t`$, where the time step $`\delta t`$ is much smaller than the mean free time, by applying a streaming step followed by a collision step. In the streaming step, the particles are moved ballistically, $`y_iy_i+v_{iy}\delta t`$. In addition, those particles crossing the boundaries are reentered with velocities sampled from the corresponding probability distribution $`K_{0,L}(𝐯)`$. Suppose that a particle $`i`$ crosses the lower plate between times $`t`$ and $`t+\delta t`$, i.e., $`y_i+v_{iy}\delta t<0`$. Then, regardless of the incoming velocity $`𝐯_i`$, a new velocity $`\stackrel{~}{𝐯}_i`$ is assigned according to the following rules. First, a velocity $`\stackrel{~}{𝐯}_i`$ (with $`\stackrel{~}{v}_{iy}>0`$) is sampled with a probability proportional to $`|v_y|\varphi _0^{\text{MB}}(𝐯)`$. If one is considering the MB boundary conditions, this velocity is accepted directly. Otherwise, the above acts as a “filter” to optimize the acceptance-rejection procedure and the velocity $`\stackrel{~}{𝐯}_i`$ is accepted with a probability proportional to the ratio $`\varphi _0^{\text{BGK}}(𝐯)/\varphi _0^{\text{MB}}(𝐯)`$. If the velocity is rejected, a new velocity $`\stackrel{~}{𝐯}_i`$ is sampled and the process is repeated until acceptance. The new position is assigned as $`\stackrel{~}{v}_{iy}(\delta t+y_i/v_{iy})`$. The process is analogous in the case of the upper plate.

The collision step proceeds as follows For each layer $`\alpha `$, a pair of potential collision partners $`i`$ and $`j`$ are chosen at random with equiprobability. The collision between those particles is then accepted with a probability equal to the corresponding collision rate times $`\delta t`$. For hard spheres, the collision rate is proportional to the relative velocity $`|𝐯_i𝐯_j|`$, while it is independent of the relative velocity for Maxwell molecules (an angle cut-off is needed in the latter case). If the collision is accepted, the scattering direction is randomly chosen according to the interaction law and post-collisional velocities are assigned to both particles, according to the conservation of momentum and energy. After the collision is processed or if the pair is rejected, the routine moves again to the choice of a new pair until the required number of candidate pairs has been taken.

In the course of the simulations, the following “coarse-grained” local quantities are computed. The number density in layer $`\alpha `$ is

$$n_\alpha =\overline{n}\frac{N_\alpha }{(N/L)\delta y}=\frac{\overline{n}L}{N\delta y}\underset{i=1}{\overset{N}{}}\mathrm{\Theta }_\alpha (y_i),$$
(53)
where $`\mathrm{\Theta }_\alpha (y)`$ is the characteristic function of layer $`\alpha `$, i.e., $`\mathrm{\Theta }_\alpha (y)=1`$ if $`y`$ belongs to layer $`\alpha `$ and is zero otherwise. Similarly, the flow velocity, the temperature, the pressure tensor, and the heat flux of layer $`\alpha `$ are

$$𝐮_\alpha =\frac{1}{N_\alpha }\underset{i=1}{\overset{N}{}}\mathrm{\Theta }_\alpha (y_i)𝐯_i,$$
(54)
$$k_BT_\alpha =\frac{p_\alpha }{n_\alpha }=\frac{m}{3N_\alpha }\underset{i=1}{\overset{N}{}}\mathrm{\Theta }_\alpha (y_i)(𝐯_i𝐮_\alpha )^2,$$
(55)
$$𝖯_\alpha =\frac{L}{N\delta y}m\underset{i=1}{\overset{N}{}}\mathrm{\Theta }_\alpha (y_i)(𝐯_i𝐮_\alpha )(𝐯_i𝐮_\alpha ),$$
(56)
$$𝐪_\alpha =\frac{L}{N\delta y}\frac{m}{2}\underset{i=1}{\overset{N}{}}\mathrm{\Theta }_\alpha (y_i)(𝐯_i𝐮_\alpha )^2(𝐯_i𝐮_\alpha ).$$
(57)
From these quantities one can get local values of the gradients and of the transport coefficients. For instance, the reduced shear rate is

$$a_\alpha =\frac{\overline{n}}{\overline{\nu }n_\alpha }\left(\frac{T_L}{T_\alpha }\right)^\omega \frac{u_{\alpha +1,x}u_{\alpha ,x}}{\delta y}$$
(58)
and the viscosity function is

$$F_{\eta ,\alpha }=\frac{P_{\alpha ,xy}}{a_\alpha p_\alpha }.$$
(59)
As said before, in the simulations we take units such that $`m=1`$, $`T_L=1`$, and $`k_B=\frac{1}{2}`$. It remains to fix the time unit or, equivalently, the length unit. The standard definition of mean free path in the case of hard spheres is

$$\lambda =\frac{1}{\sqrt{2}n\pi \sigma ^2},$$
(60)
where $`\sigma `$ is the diameter of the spheres. The Navier-Stokes shear viscosity is (in the first Sonine approximation) $`\eta _0=5(mk_BT/\pi )^{1/2}/16\sigma ^2`$. Consequently, the effective collision frequency $`\nu =p/\eta _0`$ and the mean free path $`\lambda `$ are related as $`\nu =(8/5\sqrt{\pi })(2k_BT/m)^{1/2}/\lambda `$. As usual, we choose the mean free path corresponding to the average density $`\overline{n}`$ as the length unit. This in turn implies that $`\overline{\nu }=8/5\sqrt{\pi }0.903`$. For convenience, we take the latter value for Maxwell molecules as well. The typical values of the simulation parameters are $`N=2\times 10^5`$ particles, a layer width $`\delta y=0.02`$, and a time step $`\delta t=0.003`$.

The procedure to measure the relevant quantities of the problem is as follows. First the values of the imposed shear rate $`a^{}`$ and the temperature difference $`\mathrm{\Delta }`$ are chosen. This choice fixes the system size $`L`$, as well as the upper velocity $`U_L`$ and the upper thermal gradient $`ϵ_L`$, according to Eqs. (49), (50), and (52). Starting from an equilibrium initial state with $`T=T_0`$, the system evolves driven by the boundary conditions described before. After a transient period (typically up to $`t=25`$), the system reaches a steady state in which the quantities fluctuate around constant values. In this state, the balance equations predict that the quantities $`u_y`$, $`P_{xy}`$, $`P_{yy}`$, and $`u_xP_{xy}+q_y`$ are spatially uniform and this is used in the simulations as a test to make sure that the steady state has been achieved. Once the steady state is reached, the local quantities (53)–(59) are averaged over typically 100 snapshots equally spaced between $`t=25`$ and $`t=55`$, which corresponds to about 60 collisions per particle in the case of hard spheres.

### C Test of the numerical algorithm

Before closing this Section, it is worthwhile carrying out a test of the reliability of the numerical method. To that end, we have simulated the BGK equation by a DSMC-like method similar to the one described in Ref. . If the boundary conditions are implemented correctly and the simulation parameters are well chosen, then the simulation results should agree with the theoretical BGK predictions. We have checked that this indeed the case. As an example, consider the hard-sphere situation with $`a^{}=1`$ and $`\mathrm{\Delta }=5`$. This corresponds to $`\gamma _{\text{BGK}}=0.248`$, $`L=1.81`$, $`U_L=3.17`$, and $`ϵ_L=3.15`$, where in this case $`\text{Pr}=1`$. Figure 1 shows the marginal velocity distributions of particles reemitted by the walls,

$$𝒦_{0,L}(\xi _y)=\left(\frac{2k_BT_{0,L}}{m}\right)^{1/2}_{\mathrm{}}^{\mathrm{}}𝑑v_x_{\mathrm{}}^{\mathrm{}}𝑑v_zK_{0,L}(𝐯),$$
(61)
as functions of $`\xi _y=(m/2k_BT_{0,L})^{1/2}v_y`$. The case $`\xi _y>0`$ ($`\xi _y<0`$) corresponds to particles that are reemitted from the lower (upper) plate. The agreement with the imposed distribution is excellent. Note the strong asymmetry between the distributions corresponding to both plates, in contrast to the symmetry of the MB distributions obtained from (42). The temperature and velocity profiles are shown in Fig. 2. The simulation values overlap, within statistical fluctuations, with the theoretical predictions. Apart from the profiles, we have verified that the generalized transport coefficients obtained from the fluxes also agree with the theory. A more stringent test is provided in Fig. 3, where the marginal distribution function

$$\phi (\xi _y)=\frac{1}{n}\left(\frac{2k_BT}{m}\right)^{1/2}_{\mathrm{}}^{\mathrm{}}𝑑v_x_{\mathrm{}}^{\mathrm{}}𝑑v_zf(𝐯),\xi _y=\left(\frac{m}{2k_BT}\right)^{1/2}v_y$$
(62)
is plotted at the point $`y/L=0.5`$, which corresponds to a local thermal gradient $`ϵ=0.60`$. It is apparent again that an excellent agreement exists between simulation and theory.

## IV Results

This Section is devoted to a comparison between the simulation results for the Boltzmann equation obtained by the simulation method described in the previous Section and the theoretical predictions provided by the Grad method and the BGK and ES kinetic models. The comparison will be carried out at the levels of the transport coefficients and the velocity distribution, both for Maxwell molecules and hard spheres. Before that, the hydrodynamic profiles obtained from the simulations by using the two types of boundary conditions considered in Sec. III are presented.

### A Hydrodynamic profiles

As said in Secs. I and II, the Boltzmann equation for Maxwell molecules admits an exact solution characterized by Eqs. (17)–(19). This solution applies to the bulk region, i.e., the region where boundary effects are negligible. Obviously, in a simulation with a finite size of the system, it is not possible to avoid boundary effects completely. On the other hand, one can expect that the “nonequilibrium” boundary conditions based on the BGK distribution, Eq. (45), inhibit the boundary effects, as compared with the conventional “equilibrium” boundary conditions (42). We have confirmed that this is indeed the case. As an illustrative example, let us consider $`\mathrm{\Delta }=5`$ and $`a^{}=0.92`$ for Maxwell molecules. This corresponds to $`\gamma _{\text{BGK}}=0.21`$, $`L=4.68`$, $`U_L=3.90`$, and $`ϵ_L=2.37`$. Figure 4 shows the temperature and velocity profiles as obtained by using the MB and BGK boundary conditions. It is apparent that the velocity slips and the temperature jumps at the walls are much larger in the former case than in the latter. Note that the maximum temperature is not exactly located in the layer adjacent to the lower plate but it is slightly shifted. It is interesting to remark that when plotting the temperature $`T`$ as a function of the flow velocity $`u_x`$, a parabolic curve is observed with both boundary conditions, as expected from Eqs. (18) and (19). A more evident proof of the advantage of the BGK conditions is shown in Fig. 5. Since the pressure is a constant in the exact solution valid in the bulk, any deviation from $`p=\text{const}`$ can be attributed to boundary effects. The pressure obtained with the BGK boundary conditions is practically constant, except near the upper plate, while the one obtained with the MB conditions is only nearly constant in a small region around $`y/L0.75`$. Finally, Fig. 6 shows the ratio $`a/a^{}`$ between the actual (local) shear rate $`a`$ measured in the simulations, cf. Eq. (58), and the imposed shear rate $`a^{}`$. The ratio is closer to 1 in the case of the BGK boundary conditions than in that of the MB conditions. In addition, in the former case the region where $`a`$ is practically constant extends to layers closer to the lower plate.

In summary, the above example illustrates that the BGK boundary conditions are much more efficient than the MB ones to measure transport properties in the bulk. Therefore, in what follows we will only consider the BGK conditions. In each case, we identify a bulk domain comprised between the layers $`y=y_0`$ and $`y=y_1`$ where $`a\text{const}`$, $`p\text{const}`$, and $`\gamma \text{const}`$, and take averages of $`a`$, $`F_\eta `$, $`F_\kappa `$, $`\mathrm{\Psi }_{1,2}`$, and $`\mathrm{\Phi }`$ over those layers. Typical values are $`y_0/L0.2`$ and $`y_1/L0.8`$.

### B Nonlinear transport coefficients

In this subsection we compare the simulation results for Maxwell molecules and hard spheres obtained from the Monte Carlo simulations with the (universal) predictions of the Grad method and the BGK and ES kinetic models. As said in the Introduction, we are interested in situations where the Knudsen number is not small, so that nonlinear effects are important. An interesting quantity in the nonlinear Couette flow is the parameter $`\gamma `$, which is related to the curvature of the temperature profile. Its shear-rate dependence is shown in Fig. 7. A remarkable feature is that the simulation data for both interactions seem to lie on a common curve. This indicates that, as predicted by the models, the transport properties are hardly sensitive to the interaction potential, provided that the quantities are conveniently scaled. While the kinetic models exhibit a good quantitative (in the case of the ES model) or qualitative (BGK model) agreement, the Grad method fails, except for small shear rates. Since $`\gamma (a)=a^2F_\eta (a)/5F_\kappa (a)`$, one can interpret $`5\gamma (a)/a^2`$ as an effective, shear-rate dependent Prandtl number (relative to the usual Pr). This quantity is bounded between 1 for small shear rates and $`5/3`$ (in the BGK model) or $`5/2`$ (in the ES model) in the limit of large shear rates. This is consistent with the fact that the BGK model underestimates the value of $`\gamma `$.

The most important transport coefficient is the nonlinear viscosity represented by the function $`F_\eta (a)`$. According to Fig. 8, the three theories retain the qualitative trends of the simulation data, namely the decrease of $`F_\eta `$ with increasing $`a`$ (shear thinning effect). In general, however, the kinetic models (especially the BGK model) have a better agreement than the Grad method. It is also interesting to remark that a slight influence of the interaction potential seems to exist, the shear thinning effect being a little bit more significant for hard spheres than for Maxwell molecules. The nonlinear thermal conductivity $`F_\kappa (a)`$ is plotted in Fig. 9. It is quite apparent that Grad’s solution does not capture even the qualitative shape of $`F_\kappa `$, as was already noted in the case of hard disks. Again, the kinetic models present a good agreement, especially in the case of the ES model.

Normal stress differences are characterized by the viscometric functions $`\mathrm{\Psi }_{1,2}(a)`$. These quantities are well described by the kinetic models, as shown in Figs. 10 and 11. At a quantitative level, the agreement is better in the case of Maxwell molecules. In fact, the viscometric functions, especially the second one, exhibit a certain influence of the potential. Although the functions $`\mathrm{\Psi }_{1,2}(a)`$ were not evaluated from the Grad method in Refs.  and , the friction function $`F_\mu (a)`$ was calculated. This quantity is plotted in Fig. 12, showing an agreement between the theories and the simulation data similar to the one found in Fig. 8 for the viscosity function. The last transport coefficient is the cross-coefficient $`\mathrm{\Phi }`$ defined by Eq. (25). The comparison with simulation results for this quantity is a stringent test of the theories, since it is a generalization of a Burnett coefficient that measures the component of the heat flux orthogonal to the thermal gradient. Figure 13 shows that, as happened with the thermal conductivity function $`F_\kappa (a)`$, the Grad method gives a wrong prediction for the shear-rate dependence of $`\mathrm{\Phi }(a)`$. On the other hand, the kinetic models describe fairly well the nonlinear behavior of this function. In the case of the ES model, the agreement with the simulation data is practically perfect.

### C Velocity distribution function

Apart from the transport coefficients, the kinetic models provide the velocity distribution function. In the case of the BGK model, the solution is given by Eq. (34), while the ES distribution function can be found in Refs.  and . In order to assess their reliability, we have computed in the simulations the marginal distribution (62) at the layer $`y/L=0.5`$. The ratio $`R(\xi _y)\phi (\xi _y)/[\pi ^{1/2}\mathrm{exp}(\xi _y^2)]`$, where $`\phi (\xi _y)`$ is defined by Eq. (62), is a measure of the departure of the distribution function from the local equilibrium. This quantity is plotted in Fig. 14 for Maxwell molecules in the case of a reduced shear rate $`a=0.636`$ and a reduced (local) thermal gradient $`ϵ=0.272`$. Both theories capture the main features of the actual distribution. While the BGK distribution exhibits a better agreement near the maximum (around $`\xi _y0.5`$), the ES distribution seems to describe better the region $`\xi _y1`$. The case of hard spheres is illustrated in Fig. 15, which corresponds to $`a=0.419`$ and $`ϵ=0.195`$. In this case, the ES model shows a superiority over the BGK model both near the maximum and for large positive velocities.

## V Concluding remarks

This paper has dealt with the steady planar Couette flow in a dilute gas beyond the scope of the Navier-Stokes description. This nonlinear problem had been previously studied by means of kinetic theory tools, such as the Grad method, and the BGK and ES kinetic models. These theories predict momentum and heat fluxes characterized by five shear-rate dependent generalized transport coefficients: a viscosity function $`F_\eta (a)`$, Eq. (20), a thermal conductivity function $`F_\kappa (a)`$, Eq. (21), two viscometric functions $`\mathrm{\Psi }_{1,2}(a)`$, Eqs. (22) and (23), and a cross coefficient $`\mathrm{\Phi }(a)`$, Eq. (25). The main motivation of our study has been to perform DSMC simulations for Maxwell molecules and hard spheres in order to assess the reliability of the above theories in the non-Newtonian regime. Since we have been interested in the bulk properties, we have used “nonequilibrium” boundary conditions to inhibit the influence of finite-size effects.

An important outcome of the simulation results is that, as predicted by the kinetic theories considered here, the shear-rate dependence of the transport coefficients is practically insensitive to the interaction potential when the quantities are properly nondimensionalized. In particular, the actual shear rate has been reduced with respect to an effective collision frequency defined from the Navier-Stokes shear viscosity coefficient. The simulation results show, however, that the second viscometric function presents a non negligible influence of the interaction model, so that the normal stress difference $`P_{zz}P_{yy}`$ is smaller for Maxwell molecules than for hard spheres. The comparison with the theoretical predictions shows that the kinetic models give a fairly good description of the five transport coefficients. On the other hand, the Grad method yields a shear viscosity in qualitative agreement with the simulations but it dramatically fails for the coefficients measuring the heat flux. This is basically due to the truncation scheme of the Grad method at the level of the heat flux. The physical idea behind a kinetic model is quite different, since it consists of replacing the true Boltzmann collision operator by a simple relaxation term but otherwise all the velocity moments are taken into account. As a consequence, while in the Grad method one has to solve a closed set of coupled differential equations for the moments, in the case of the kinetic model one gets the velocity distribution function and determines the fluxes from it.

In the ES kinetic model the reference distribution function appearing in the collision operator is an anisotropic Gaussian parameterized by the pressure tensor. This allows the model to give the correct Prandtl number $`\text{Pr}=\frac{2}{3}`$ but at the expense of complicating its mathematical structure. In the case of the BGK model, however, the reference distribution is that of local equilibrium but the model leads to $`\text{Pr}=1`$. The agreement with simulation of the ES model is generally better than that of the BGK model, especially in the case of $`F_\kappa `$ and $`\mathrm{\Phi }`$. In spite of this, it is fair to say that the performance of the BGK model is quite good, given its simplicity relative to that of the ES model. Finally, the results reported in this paper clearly shows the usefulness of kinetic models to analyze nonlinear transport phenomena in the Couette flow problem. This complements previous conclusions drawn from other nonlinear problems, such as the uniform shear flow and the Fourier flow.

###### Acknowledgements.

The authors 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.

## Theoretical expressions for the transport coefficients

In this Appendix we list the explicit shear-rate dependence of the dimensionless transport coefficients defined in Sec. II, according to the Grad method, the BGK model, and the ES model.

### 1 The Grad method

From the Appendix of Ref. , corrected in Ref. , one has

$$F_\eta (a)=\frac{2}{1+\frac{72}{25}a^2+\mathrm{\Delta }(a)},$$
(63)
$$F_\kappa (a)=\frac{4}{1\frac{216}{25}a^2+3\mathrm{\Delta }(a)},$$
(64)
$$\mathrm{\Phi }(a)=7\frac{1\frac{36}{125}a^2}{1+\frac{6}{5}a^2+(1\frac{63}{25}a^2)\mathrm{\Delta }(a)},$$
(65)
where $`\mathrm{\Delta }(a)\sqrt{1+\frac{116}{25}a^2\frac{864}{625}a^4}`$. The viscometric functions are not evaluated in Refs.  and , although the friction function is provided. It is given by

$$F_\mu (a)=\frac{2}{1+\frac{12}{25}a^2+\mathrm{\Delta }(a)}.$$
(66)
Note that $`F_\eta (a)`$ and $`F_\mu (a)`$ become meaningless for $`a^225(\sqrt{1057}+29)/4323.56`$, while $`F_\kappa (a)`$ and $`\mathrm{\Phi }(a)`$ are unphysical for $`a^250/630.79`$. For small shear rates, the above transport coefficients behave as $`F_\eta 1\frac{13}{5}a^2`$, $`F_\kappa 1+\frac{21}{50}a^2`$, $`\mathrm{\Phi }\frac{7}{2}(1\frac{197}{250}a^2)`$, and $`F_\mu 1\frac{7}{5}a^2`$.

### 2 The BGK model

The derivation of the transport coefficients from the BGK model  implies the resummation of asymptotic series by means of the Borel method. As a consequence, the results are expressed in terms of the functions $`F_r(x)`$ defined by the recurrence relation $`F_r(x)=[(d/dx)x]^rF_0(x)`$, where

$$F_0(x)=\frac{2}{x}_0^{\mathrm{}}𝑑tte^{t^2/2}K_0(2t^{1/2}/x^{1/4}),$$
(67)
$`K_0`$ being the zeroth-order modified Bessel function. The curvature parameter $`\gamma (a)`$ is given by the solution to the following implicit equation

$$a^2=\gamma \frac{3F_1+2F_2}{F_1},$$
(68)
where the functions $`F_r`$ are evaluated at $`x=\gamma `$. Next, the transport coefficients are expressed in terms of $`a`$ and $`F_r(\gamma )`$ as

$$F_\eta (a)=F_0,$$
(69)
$$F_\kappa (a)=\frac{F_0}{5}\frac{3F_1+2F_2}{F_1},$$
(70)
$$\mathrm{\Psi }_1(a)=2F_1\frac{3F_1+4F_2}{3F_1+2F_2},$$
(71)
$$\mathrm{\Psi }_2(a)=4\frac{F_1F_2}{3F_1+2F_2},$$
(72)
$$\mathrm{\Phi }(a)=\left[5F_2+2F_3+a^2(F_2+5F_3+8F_4+4F_5)\right].$$
(73)
For small shear rates, one gets $`F_\eta 1\frac{18}{5}a^2`$, $`F_\kappa 1\frac{162}{25}a^2`$, $`\mathrm{\Psi }_1\frac{14}{5}(1\frac{1476}{175}a^2)`$, $`\mathrm{\Psi }_2\frac{4}{5}(1\frac{288}{25}a^2)`$, $`\mathrm{\Phi }\frac{14}{5}`$, and $`F_\mu 1\frac{12}{5}a^2`$.

### 3 The ES model

In Refs.  and  the solution of the ES model is worked out keeping the Prandtl number Pr as a free parameter. Here we particularize the results to the correct value $`\text{Pr}=\frac{2}{3}`$. It is convenient to express the transport coefficients in terms of an auxiliary parameter $`\beta `$, defined as the solution of the implicit equation

$$a^2=\frac{4\beta }{9}\frac{[2\beta (F_1+2F_2)3]^2[3F_1+2F_22\beta F_1(F_1+2F_2)]}{F_0^2[2\beta (F_1+F_2)3]+F_1[2\beta (F_1+2F_2)3]^2},$$
(74)
where now the functions $`F_r`$ are evaluated at $`x=\beta `$. The relationship between the curvature parameter $`\gamma `$ and $`\beta `$ is

$$\gamma (a)=\frac{2}{9}\beta [32\beta (F_1+2F_2)].$$
(75)
The transport coefficients are

$$F_\eta (a)=\frac{9F_0}{[2\beta (F_1+2F_2)3]^2},$$
(76)
$$F_\kappa (a)=\frac{a^2}{5\gamma (a)}F_\eta (a),$$
(77)
$$\mathrm{\Psi }_1(a)=\frac{12\beta }{a^2}\frac{3F_1+4F_22\beta F_1(F_1+2F_2)}{(32\beta F_1)[32\beta F_1(F_1+2F_2]},$$
(78)
$$\mathrm{\Psi }_2(a)=\frac{24\beta }{a^2}\frac{F_2}{(32\beta F_1)[32\beta F_1(F_1+2F_2]},$$
(79)
$`\mathrm{\Phi }(a)`$ $`=`$ $`{\displaystyle \frac{1}{40}}C_1\{a^2[C_1^4F_0^3(F_1+2F_2)9C_1^3F_0^2(F_2+2F_3)+54C_1^2F_0(F_2+4(F_3+F_4))`$ (82)
$`108C_1(F_2+5F_3+4(2F_4+F_5))]C_3(C_1F_0F_16F_2)+4C_1^2F_0(F_1+2F_2)`$
$`+4C_1[C_2F_0F_16(F_2+2F_3)]24C_2F_2\}.`$
In Eq. (82),

$$C_1\frac{3}{32\beta (F_1+2F_2)},$$
(83)
$$C_2\frac{3}{32\beta F_1},$$
(84)
$$C_33a^2C_1^3F_0^2+4C_2\{2\beta [C_1(F_1+4F_2)+3F_1]3C_1\}.$$
(85)
For small shear rates, one has $`F_\eta 1\frac{21}{5}a^2`$, $`F_\kappa 1\frac{197}{25}a^2`$, $`\mathrm{\Psi }_1\frac{14}{5}(1\frac{2126}{175}a^2)`$, $`\mathrm{\Psi }_2\frac{4}{5}(1\frac{413}{25}a^2)`$, $`\mathrm{\Phi }\frac{7}{2}`$, and $`F_\mu 13a^2`$.