Solutions of the moment hierarchy in the kinetic theory of Maxwell models

# Solutions of the moment hierarchy in the kinetic theory of Maxwell models

## Abstract

In the Maxwell interaction model the collision rate is independent of the relative velocity of the colliding pair and, as a consequence, the collisional moments are bilinear combinations of velocity moments of the same or lower order. In general, however, the drift term of the Boltzmann equation couples moments of a given order to moments of a higher order, thus preventing the solvability of the moment hierarchy, unless approximate closures are introduced. On the other hand, there exist a number of states where the moment hierarchy can be recursively solved, the solution generally exposing non-Newtonian properties. The aim of this paper is to present an overview of results pertaining to some of those states, namely the planar Fourier flow (without and with a constant gravity field), the planar Couette flow, the force-driven Poiseuille flow, and the uniform shear flow.

###### Keywords:
Boltzmann equation Maxwell molecules Moment hierarchy Non-Newtonian properties
###### pacs:
05.20.Dd 05.60.-k 51.10.+y 51.20.+d 47.50.-d
2

## 1 Introduction

As is well known, the classical kinetic theory of low-density gases began to be established as a mathematically sound statistical-physical theory with the work of James Clerk Maxwell (1831–1879) (15); (16); (17). Apart from obtaining the velocity distribution function at equilibrium (first in 1860 and then, in a more rigorous way, in 1867), Maxwell derived in 1866 and 1867 the transfer equations characterizing the rate of change of any quantity (such as mass, momentum, or energy) which can be defined in terms of molecular properties. This paved the way to Ludwig Boltzmann (1844–1906) in the derivation of his celebrated equation (1872) for the rate of change of the velocity distribution itself. As a matter of fact, the Boltzmann equation is formally equivalent to Maxwell’s infinite set of transfer equations.

In his work entitled “On the Dynamical Theory of Gases,” Maxwell departs from the hard-sphere model and writes (15); (17)

“In the present paper I propose to consider the molecules of a gas, not as elastic spheres of definite radius, but as small bodies or groups of smaller molecules repelling one another with a force whose direction always passes very nearly through the centres of gravity of the molecules, and whose magnitude is represented very nearly by some function of the distance of the centres of gravity. I have made this modification of the theory in consequence of the results of my experiments on the viscosity of air at different temperatures, and I have deduced from these experiments that the repulsion is inversely as the fifth power of the distance.”

Maxwell realized that the hypothesis of a force inversely proportional to the fifth power of the distance, or, equivalently, of an interaction potential , makes the nonequilibrium distribution function enter the transfer equations in such a way that the transport coefficients can be evaluated without the need of knowing the detailed form of . His results with this interaction model showed that the shear viscosity was proportional to the absolute temperature , whereas in the case of hard spheres it is proportional to the square root of the temperature. As said in the preceding quotation, Maxwell himself carried out a series of experiments to measure the viscosity of air as a function of temperature. He writes (15); (17)

“I have found by experiment that the coefficient of viscosity in a given gas is independent of the density, and proportional to the absolute temperature, so that if be the viscosity, .”

Modern measurements of the viscosity of air show that it is approximately proportional to over a wide range of temperatures (5); (6), so that the actual power is intermediate between the values predicted by the hard-sphere and Maxwell models.

In any case, the adoption of Maxwell’s interaction model allows one to extract some useful information from the nonlinear Boltzmann equation, or its associated set of moment equations, for states far from equilibrium. Apart from its intrinsic interest, exact results derived for Maxwell molecules are important as benchmarks to assess approximate moment methods, model kinetic equations, or numerical algorithms (deterministic or stochastic) to solve the Boltzmann equation. Moreover, it turns out that many consequences of the Boltzmann equation for Maxwell molecules, when properly rescaled with respect to the collision frequency, can be (approximately) extrapolated to more general interaction potentials (32).

The aim of this paper is to review a few examples of nonequilibrium states whose hierarchy of moment equations can be recursively solved for the Maxwell model. The structure and main properties of the moment equations for Maxwell molecules are recalled in Section 2. This is followed by a description of the solution for the planar Fourier flow (Section 3), the planar Couette flow (Section 4), the force-driven Poiseuille flow (Section 5), and the uniform shear flow (Section 6). The paper is closed with some concluding remarks in Section 7.

## 2 Moment equations for Maxwell molecules

### 2.1 The Boltzmann equation

Let us consider a dilute gas made of particles of mass interacting via a short-range pair potential . The relevant statistical-mechanical description of the gas is conveyed by the one-particle velocity distribution function . The number density , the flow velocity , and the temperature are related to through

 n(r,t)=∫dvf(r,v,t), (1)
 u(r,t)=1n(r,t)∫dvvf(r,v,t), (2)
 n(r,t)kBT(r,t)=p(r,t)=m3∫dvV2(v,r,t)f(r,v,t), (3)

where is the Boltzmann constant and is the so-called peculiar velocity. The fluxes of momentum and energy are measured by the pressure (or stress) tensor and the heat flux vector , respectively. Their expressions are

 P(r,t)=m∫dvV(v,r,t)V(v,r,t)f(r,v,t), (4)
 Extra open brace or missing close brace (5)

The time evolution of the velocity distribution function is governed by the Boltzmann equation (21); (23)

 ∂∂tf=−v⋅∇f−∂∂v⋅(Fmf)+J[v|f,f]. (6)

The first and second terms on the right-hand side represent the rate of change of due to the free motion and to the possible action of an external force , respectively. In general, such a force can be non-conservative and velocity-dependent, and in that case it must appear to the right of the operator . The last term on the right-hand side is the fundamental one. It gives the rate of change of due to the interactions among the particles treated as successions of local and instantaneous binary collisions. The explicit form of the (bilinear) collision operator is

 J[v|f,f]=∫dv1∫dΩwB(w,χ)[f(v′)f(v′1)−f(v)f(v1)], (7)

where is the pre-collisional relative velocity,

 v′=v−(w⋅ˆσ)ˆσ,v′1=v1+(w⋅ˆσ)ˆσ (8)

are the post-collisional velocities, is a unit vector pointing on the apse line (i.e., the line joining the centers of the two particles at their closest approach), is the scattering angle (i.e., the angle between the post-collisional relative velocity and ), and is an element of solid angle about the direction of . Finally, is the differential cross section (33). This is the only quantity that depends on the choice of the potential :

 B(w,χ)=b(w,χ)sinχ∣∣∣∂b(w,χ)∂χ∣∣∣, (9)

where the impact parameter is obtained from inversion of

 χ(w,b)=π−2∫∞r0(w,b)drb/r2[1−(b/r)2−4φ(r)/mw2]1/2, (10)

being the distance at closest approach, which is given as the root of .

In the special case of inverse power-law (IPL) repulsive potentials of the form , the scattering angle depends on both and through the scaled dimensionless parameter (, , and being constants), namely

 χ(β)=π−2∫β0(β)0dβ′[1−β′2−2ζ(β′β)ζ]−1/2, (11)

where is the root of the quantity enclosed by brackets. For this class of potentials the differential cross section has the scaling form

 B(w,χ)=σ2(2ζφ0mw2)2/ζB(χ),B(χ)=β(χ)sinχ∣∣∣dβ(χ)dχ∣∣∣. (12)

The hard-sphere potential is recovered in the limit , in which case , , , and . On the other hand, in the case of the Maxwell potential (), the collision rate is independent of the relative speed , namely

 wB(w,χ)=QB(χ),Q≡σ2√8φ0/m. (13)

Maxwell himself realized that if the integral in Eq. (11) can be expressed in terms of the complete elliptic integral of the first kind , namely

 χ(β)=π−2β(2+β4)1/4K(12−β22(2+β4)1/2). (14)

Obviously, the associated function for the IPL Maxwell potential becomes rather cumbersome and needs to be evaluated numerically. Nevertheless, it is sometimes convenient to depart from a strict adherence to the interaction potential by directly modeling the scattering angle dependence of the differential cross section (25). In particular, in the variable soft-sphere (VSS) model one has (47); (48)

 β(χ)=cosϑχ2,B(χ)=ϑ4cos2(ϑ−1)χ2. (15)

If one chooses , then the scattering is assumed to be isotropic (38) and one is dealing with the variable hard-sphere (VHS) model (12). However, this leads to a viscosity/self-diffusion ratio different from that of the true IPL Maxwell interaction. To remedy this, in the VSS model (47) one takes . In the context of inelastic Maxwell models (10), it is usual to take .

### 2.2 Moment equations

Let be a test velocity function. Multiplying both sides of the Boltzmann equation (6) by and integrating over velocity one gets the balance equation

 ∂Ψ∂t+∇⋅Φψ=σ(F)ψ+Jψ, (16)

where

 Ψ=∫dvψ(v)f(v)≡n⟨ψ⟩ (17)

is the local density of the quantity represented by ,

 Φψ=∫dvvψ(v)f(v) (18)

is the associated flux,

 σ(F)ψ=∫dv∂ψ∂v⋅[Fmf(v)] (19)

is a source term due to the external force, and

 Jψ = ∫dvψ(v)J[v|f,f] (20) = 14∫dv∫dv1∫dΩwB(w,χ)[ψ(v)+ψ(v1)−ψ(v′)−ψ(v′1)]f(v)f(v1)

is the source term due to collisions. The general balance equation (16), which is usually referred to as the weak form of the Boltzmann equation, is close to the approach followed by Maxwell in 1867 (17); (95). We say that is a moment of order if is a polynomial of degree . In that case, if the external force is independent of velocity, the source term is a moment of order . The hierarchical structure of the moment equations is in general due to the flux and the collisional moment . The former is a moment of order , while the latter is a bilinear combination of moments of any order because of the velocity dependence of the collision rate . In order to get a closed set of equations some kind of approximate closure needs to be applied. In the Hilbert and Chapman–Enskog (CE) methods (13); (21); (77) one focuses on the balance equations for the five conserved quantities, namely , so that . Next, an expansion of the velocity distribution function in powers of the Knudsen number (defined as the ratio between the mean free path and the characteristic distance associated with the hydrodynamic gradients) provides the Navier–Stokes (NS) hydrodynamic equations and their sequels (Burnett, super-Burnett, …). On a different vein, Grad proposed in 1949 (35); (36) to expand the distribution function in a complete set of orthogonal polynomials (essentially Hermite polynomials), the coefficients being the corresponding velocity moments. Next, this expansion is truncated by retaining terms up to a given order , so the (orthogonal) moments of order higher than are neglected and one finally gets a closed set of moment equations. In the usual 13-moment approximation, the expansion includes the density , the three components of the flow velocity , the six elements of the pressure tensor , and the three components of the heat flux . The method can be augmented to twenty moments by including the seven third-order moments apart from the heat flux. Variants of Grad’s moment method have been developed in the last few years (79), this special issue reflecting the current state of the art.

As said above, the collisional moment involves in general moments of every order. An important exception takes place in the case of Maxwell models, where is independent of the relative speed [cf. Eq. (12) with ]. This implies that, if is a polynomial of degree , then becomes a bilinear combination of moments of order equal to or less than (38); (52); (95). To be more precise, let us define the reduced orthogonal moments (19); (35); (52); (62); (95)

 ψkℓμ(c)=NkℓcℓL(ℓ+12)k(c2)Yμℓ(ˆc), (21)

where

 c=√m2kBTV (22)

is the peculiar velocity normalized with respect to the thermal velocity, are generalized Laguerre polynomials (1); (37) (also known as Sonine polynomials), are spherical harmonics, being the unit vector along the direction of , and

 Nkℓ=⎡⎣2π3/2k!Γ(k+ℓ+32)⎤⎦1/2 (23)

are normalization constants, being the gamma function. The polynomials form a complete set of orthonormal functions with respect to the inner product . Let us denote as

 Mα=1n∫dvψα(c)f(v) (24)

the (reduced) moment of order associated with the polynomial . In the case of Maxwell models one has

 Jα = 1n∫dvψα(c)J[v|f,f] (25) = −λkℓMα+†∑α′,α′′Cαα′α′′Mα′Mα′′.

The dagger in the summation means the constraints and . Moreover, it is necessary that form a triangle (i.e., ) and (52). The explicit expression of the coefficients is rather involved and can be found in Ref. (52). A computer-aided algorithm to generate the collisional moments for Maxwell models is described in Ref. (66).

Equation (25) shows that the polynomials are eigenfunctions of the linearized collision operator for Maxwell particles, the eigenvalues being given by

 λkℓ=2πnQ∫π0dχsinχB(χ)[1+δr0δℓ0−cos2k+ℓχ2Pℓ(cosχ2)−sin2k+ℓχ2Pℓ(sinχ2)], (26)

where are Legendre polynomials. Irrespective of the precise angular dependence of , the following relationships hold

 λk,1=λk+1,0,(2ℓ+1)λkℓ=(ℓ+1)λk−1,ℓ+1+ℓλk,ℓ−1. (27)

The eigenvalues can be expressed as linear combinations of the numerical coefficients

 A2a = 2π∫π0dχsinχB(χ)cosaχ2sinaχ2 (28) = 2π∫∞0dββcos2aχ(β)2sin2aχ(β)2.

In particular, . The reduced eigenvalues of order are listed in Table 1. For the IPL model the coefficients must be evaluated numerically. On the other hand, by assuming Eq. (15) one simply has

 A2a=πϑB(a+ϑ,a+1), (29)

where is the Euler beta function (1); (37). Table 2 gives the first few values of for the IPL, VSS, and VHS models. It can be observed that the latter two models provide values quite close to the correct IPL ones, especially as increases. A rather extensive table of the eigenvalues for the IPL Maxwell model can be found in Ref. (4).

The most important eigenvalues are and , which provide the NS shear viscosity and thermal conductivity, namely

 PNSij=pδij−ηNS(∇iuj+∇jui−23∇⋅uδij),ηNS% =pλ02, (30)
 qNS=−κNS∇T,κNS=52kBpmλ11. (31)

### 2.3 Solvable states

Even in the case of Maxwell models the moment hierarchy (16) couples moments of order to moments of order through the flux term . Therefore, the moment equations cannot be solved in general unless an approximate closure is introduced. There exist, however, a few steady states where the hierarchy (16) for Maxwell molecules can be recursively solved (32); (74). In those solutions one focuses on the bulk of the system (i.e., away from the boundary layers) and assumes that the velocity distribution function adopts a normal form, i.e., it depends on space only through the hydrodynamic fields (, , and ) and their gradients. On the other hand, it is not necessary to invoke that the Knudsen number (where, as said before, is the mean free path and is the characteristic distance associated with the hydrodynamic gradients) is small.

The aim of the remainder of the paper is to review some of those solutions. All of them have the common features of a planar or channel geometry (gas enclosed between infinite parallel plates) and a one-dimensional spatial dependence along the direction orthogonal to the plates. Figure 1 sketches the four steady states to be considered. In the planar Fourier flow (without gravity) the recursive solvability of the moment equations is tied to the fact that the (reduced) moments of order are just polynomials in the Knudsen number (here associated with the thermal gradient) of degree and parity . This simple polynomial dependence is broken down when a gravity field is added but the problem is still solvable by a perturbation expansion in powers of the gravity strength. When the gas is sheared by moving plates (planar Couette flow) the dependence of the reduced moments on the Knudsen number associated with the thermal gradient is still polynomial, but with coefficients that depend on the (reduced) shear rate. In the case of the force-driven Poiseuille flow the nonequilibrium hydrodynamic profiles are induced by the presence of a longitudinal body force only. Again a perturbation expansion allows one to get those profiles, which exhibit interesting non-Newtonian features.

## 3 Planar Fourier flow

### 3.1 Without gravity

In the planar Fourier flow the gas is enclosed between two infinite parallel plates kept at different temperatures [see Fig. 1(a)]. We assume that no external force is acting on the particles () and consider the steady state of the gas with gradients along the direction normal to the plates () and no flow velocity,

 u=0. (32)

Under these conditions, the Boltzmann equation (6) becomes

 vz∂∂zf(z,v)=J[v|f,f] (33)

and the conservation laws of momentum and energy yield

 Pzz=const, (34)
 qz=const, (35)

respectively.

In 1979 Asmolov, Makashev, and Nosik (9) proved that an exact normal solution of Eq. (33) for Maxwell molecules exists with

 p=nkBT=const, (36a) T(z)∂T(z)∂z=const. (36b)

The original paper is rather condensed and difficult to follow, so an independent derivation (70) is expounded here.

Since in this problem the only hydrodynamic gradient is the thermal one, i.e., , the obvious choice of hydrodynamic length is . As for the mean free path one can take . Both quantities are local and their ratio defines the relevant local Knudsen number of the problem, namely

 ϵ(z)=√2kBT(z)/mλ11(z)∂lnT(z)∂z. (37)

Note that , so that , where use has been made of Eqs. (36). Here it is assumed that the separation between the plates is large enough as compared with the mean free path to identify a bulk region where the normal solution applies. Such a solution can be nondimensionalized in the form

 ϕ(c;ϵ)=1n(z)[2kBT(z)m]3/2f(z,v), (38)

where is defined by Eq. (22). All the spatial dependence of is contained in its dependence on and . Thus

 ∂f∂z=∂T∂z∂f∂T, (39)

where

 ∂f∂T=n(2kBTm)−3/2(−52T−1ϕ+∂c∂T⋅∂ϕ∂c+∂ϵ∂T∂ϕ∂ϵ). (40)

Taking into account that and , one finally gets

 ∂∂zf(z,v)=−n(m2kBT)2λ11ϵ2(2+∂∂c⋅c+ϵ∂∂ϵ)ϕ(c;ϵ). (41)

Consequently, the Boltzmann equation (33) becomes

 Extra open brace or missing close brace = Extra open brace or missing close brace (42) ≡ ¯¯¯¯J[c|ϕ(ϵ),ϕ(ϵ)].

The orthogonal moments of are

 Mkℓ(ϵ)=∫dcψkℓ(c)ϕ(c;ϵ), (43)

so that

 ϕ(c;ϵ)=π−3/2e−c2∞∑k=0∞∑ℓ=0Mkℓ(ϵ)ψkℓ(c). (44)

Here,

 ψkℓ(c)≡ψkℓ0(c)=Nkℓ√2ℓ+14πL(ℓ+12)k(c2)cℓPℓ(cz/c). (45)

The moments are a subclass of the moments defined by Eq. (24). From the definition of temperature and density, and the fact that the flow velocity vanishes, it follows that

 M00=1,M10=M01=0. (46)

Taking into account the recurrence relations of the Laguerre and Legendre polynomials (1); (37) it is easy to prove that

 c⋅∂∂cψkℓ(c)=(2k+ℓ)ψkℓ(c)−2√k(k+ℓ+12)ψk−1,ℓ(c), (47)
 czψkℓ(c)=ϖkℓψk,ℓ+1(c)−ωk+1,ℓ−1ψk+1,ℓ−1(c)+ϖk,ℓ−1ψk,ℓ−1(c)−ωkℓψk−1,ℓ+1(c), (48)

where

 ϖkℓ≡(ℓ+1)⎡⎣k+ℓ+32(2ℓ+1)(2ℓ+3)⎤⎦1/2,ωkℓ≡(ℓ+1)[k(2ℓ+1)(2ℓ+3)]1/2. (49)

Multiplying both sides of Eq. (42) by and integrating over one gets

 ϵ2(2k+ℓ−1−ϵ∂∂ϵ)(ϖkℓMk,ℓ+1−ωk+1,ℓ−1Mk+1,ℓ−1+ϖk,ℓ−1Mk,ℓ−1−ωkℓMk−1,ℓ+1) (50) −ϵ√k(k+ℓ+12)(ϖk−1,ℓMk−1,ℓ+1−ωk,ℓ−1Mk,ℓ−1+ϖk−1,ℓ−1Mk−1,ℓ−1−ωk−1,ℓMk−2,ℓ+1) = −\lambdabarkℓMkℓ+1λ11†∑k′ℓ′k′′ℓ′′Ckℓk′ℓ′k′′ℓ′′Mk′ℓ′Mk′′ℓ′′.

The choice of the orthogonal moments simplifies the structure of the collisional terms but, on the other hand, complicates the convective terms. Alternatively, one can use the non-orthogonal moments

 ¯¯¯¯¯¯Mkℓ(ϵ)=∫dcc2kcℓzϕ(c;ϵ). (51)

In that case, one gets

 ϵ2(2k+ℓ−1−ϵ∂∂ϵ)¯¯¯¯¯¯Mk,ℓ+1(ϵ) = ∫dcc2kcℓz¯¯¯¯J[c|ϕ(ϵ),ϕ(ϵ)] (52) ≡ ¯¯¯¯Jkℓ(ϵ)

from Eq. (42). Obviously, is a bilinear combination of moments of order equal to or smaller than .

Equations (50) and (52) reveal the hierarchical character of the moment equations: moments of order are coupled to moments of lower order but also to moments of order . However, the hierarchy can be solved by following a recursive scheme. Of course, the moments of zeroth and first order, as well as one of the two moments of second order are known by definition [cf. Eq. (46)]. The key point is that the moments of order are polynomials in of degree and parity :

 Mkℓ(ϵ)=2k+ℓ−2∑j=0μ(kℓ)jϵj,μ(kℓ)j=0 if j+ℓ=odd, (53)

where are pure numbers to be determined. Let us see that Eq. (53) is consistent with Eq. (50). First note that all the moments on the left-hand side of Eq. (50) have a parity different from that of but the parity is restored when multiplying by or applying the operator . Similarly, the condition assures that the product is a polynomial of the same degree and parity as those of . Next, although the moments and are polynomials of degree , the action of the operator transforms those polynomials into polynomials of degree .

In order to complete the proof that (53) provides a solution to the moment equations (50), a recursive scheme must be devised to get the numerical coefficients . First, notice that the left-hand side of Eq. (50) contains at least the first power of . So, if one easily gets , provided that is known for . Next, if , the coefficient is determined from the previous knowledge of for and of for . Again, if and we know and for , as well as for , we can get , and so on. As a starting point (moments of zeroth degree) we have and , the latter being a consequence of . The recursive scheme is sketched in Fig. 2, where the arrows indicate the sequence followed in the determination of . The open (closed) circles represent the coefficients associated with even (odd) and . Notice that the number of coefficients needed to determine a given moment is finite. In fact, the coefficients represented in Fig. 2 are the ones involved in the evaluation of for . In practice, the number of coefficients actually needed is smaller since

 μ(kℓ)j=0 if j

To check the above property, note that Eq. (50) implies that

 μ(kℓ)j = L.C.{μ(k,ℓ+1)j−1,μ(k+1,ℓ−1)j−1,μ(k−1,ℓ+1)j−1,μ(k,ℓ−1)j−1,μ(k−2,ℓ+1)j−1,μ(k−1,ℓ−1)j−1} (55) +L.C.⎧⎨⎩j∑j′=0μ(k′ℓ′)j′μ(k′′ℓ′′)j−j′⎫⎬⎭2k′+ℓ′+2k′′+ℓ′′=2k+ℓ,|ℓ′−ℓ′′|≤ℓ≤ℓ′+ℓ′′,

where L.C. stands for “linear combination of”. If for then the first term on the right-hand side of Eq. (55) vanishes because ; the second term also vanishes because either or . Similarly, if for then both terms on the right-hand side vanish because and either or , respectively. The property for yields

 M0ℓ=0,ℓ≥1. (56)

Since the non-orthogonal moments are linear combinations of the orthogonal moments with and , the polynomial form (53) also holds for :

 ¯¯¯¯¯¯Mkℓ(ϵ)=2k+ℓ−2∑j=0¯¯¯¯¯m(kℓ)jϵj,¯¯¯¯¯m(kℓ