A Evaluation of \Phi_{n_{1}n_{2}n_{3}}^{(\text{I})} and \Phi_{n_{1}n_{2}n_{3}}^{(\text{II})}

# Non-Newtonian Couette-Poiseuille flow of a dilute gas

## Abstract

The steady state of a dilute gas enclosed between two infinite parallel plates in relative motion and under the action of a uniform body force parallel to the plates is considered. The Bhatnagar–Gross–Krook model kinetic equation is analytically solved for this Couette–Poiseuille flow to first order in the force and for arbitrary values of the Knudsen number associated with the shear rate. This allows us to investigate the influence of the external force on the non-Newtonian properties of the Couette flow. Moreover, the Couette–Poiseuille flow is analyzed when the shear-rate Knudsen number and the scaled force are of the same order and terms up to second order are retained. In this way, the transition from the bimodal temperature profile characteristic of the pure force-driven Poiseuille flow to the parabolic profile characteristic of the pure Couette flow through several intermediate stages in the Couette–Poiseuille flow are described. A critical comparison with the Navier–Stokes solution of the problem is carried out.

Bhatnagar–Gross–Krook kinetic model, Couette flow, Poiseuille flow, non-Newtonian properties.
\theoremstyle

definition Couette–Poiseuille flow] Non-Newtonian Couette–Poiseuille flow of
a dilute gas Mohamed Tij and Andrés Santos] \subjclassPrimary: 76P05, 82B40; Secondary: 82C40, 82C05.

Mohamed Tij

Département de Physique, Université Moulay Ismaïl

Meknès, Morocco

Andrés Santos

## 1 Introduction

Two paradigmatic stationary nonequilibrium flows are the plane Couette flow and the Poiseuille flow. In the plane Couette flow the fluid (henceforth assumed to be a dilute gas) is enclosed between two infinite parallel plates in relative motion, as sketched in Fig. 1(a). The walls can be kept at different or equal temperatures but, even if both wall temperatures are the same, viscous heating induces a temperature gradient in the steady state. If the Knudsen number associated with the shear rate is small enough the Navier–Stokes (NS) equations provide a satisfactory description of the Couette flow. On the other hand, as shearing increases, non-Newtonian effects (shear thinning and viscometric properties) and deviations of Fourier’s law (generalized thermal conductivity and streamwise heat flux component) become clearly apparent [20]. These nonlinear effects have been derived from the Boltzmann equation for Maxwell molecules [14, 28, 34, 40, 57], from the Bhatnagar–Gross–Krook (BGK) kinetic model [6, 19, 43], and also from generalized hydrodynamic theories [49, 51]. A good agreement with computer simulations [16, 17, 25, 30, 31, 37] has been found. The plane Couette flow has also been analyzed in the context of granular gases [60, 65]. In the case of plates at rest but kept at different temperatures, the Couette flow becomes the familiar plane Fourier flow, which also presents interesting properties by itself [4, 16, 17, 24, 29, 33, 40, 41, 42].

The Poiseuille flow, where a gas is enclosed in a channel or slab and fluid motion is induced by a longitudinal pressure gradient, is a classical problem in kinetic theory [10, 35]. Essentially the same type of flow field is generated when the pressure gradient is replaced by the action of a uniform longitudinal body force (e.g., gravity), as illustrated in Fig. 1(b). This force-driven Poiseuille flow has received a lot of attention both from theoretical [2, 3, 15, 18, 21, 32, 38, 39, 44, 49, 51, 55, 56, 58, 59, 64, 67] and computational [1, 12, 22, 23, 27, 38, 61, 62, 68] points of view. This interest has been mainly motivated by the fact that the force-driven Poiseuille flow provides a nice example illustrating the limitations of the NS description in the bulk domain (i.e., far away from the boundary layers). In particular, while the NS equations predict a temperature profile with a flat maximum at the center, computer simulations [27] and kinetic theory calculations [55, 56] show that it actually has a local minimum at that point.

Obviously, the Couette and Poiseuille flows can be combined to become the Couette–Poiseuille (or Poiseuille–Couette) flow [9, 36, 50, 53]. To the best of our knowledge, all the studies on the Couette–Poiseuille flow assume that the Poiseuille part is driven by a pressure gradient, not by an external force. This paper intends to fill this gap by considering the steady state of a dilute gas enclosed between two infinite parallel plates in relative motion, the particles of the gas being subject to the action of a uniform body force. This Couette–Poiseuille flow is sketched in Fig. 1(c). We will study the problem by the tools of kinetic theory by solving the BGK model for Maxwell molecules.

The aim of this work is two-fold. First, we want to investigate how the fully developed non-Newtonian Couette flow is distorted by the action of the external force. To that end we will assume a finite value of the Knudsen number related to the shear rate and perform a perturbation expansion to first order in the force. As a second objective, we will study how the non-Newtonian force-driven Poiseuille flow is modified by the shearing. This is done by assuming that the shear-rate Knudsen number and the scaled force are of the same order and neglecting terms of third and higher order. In both cases we are interested in the physical properties in the central bulk region of the slab, outside the influence of the boundary layers.

The organization of the paper is as follows. The Boltzmann equation for the Couette–Poiseuille flow is presented in Sec. 2. Section 3 deals with the NS description of the problem. The main part of the paper is contained in Sec. 4, where the kinetic theory approach is worked out. Some technical calculations are relegated to Appendix A. The results are graphically presented and discussed in Sec. 5. The paper ends with some concluding remarks in Sec. 6.

## 2 The Couette–Poiseuille flow. Symmetry properties

Let us consider a dilute monatomic gas enclosed between two infinite parallel plates located at . The plates are in relative motion with velocities along the axis and are kept at a common temperature . The imposed shear rate is therefore . Besides, an external body force , where is the mass of a particle and is a constant acceleration, is applied. The geometry of the problem is sketched in Fig. 1(c). In the absence of the external force () this problem reduces to the plane Couette flow [see Fig. 1(a)]. On the other hand, if the plates are at rest (), one is dealing with the force-driven Poiseuille flow [see Fig. 1(b)]. The general problem with and defines the Couette–Poiseuille flow analyzed in this paper.

In the steady state only gradients along the axis are present and thus the Boltzmann equation becomes

 (vy∂∂y+g∂∂vx)f(y,v|ω,g)=J[f,f], (1)

where is the one-particle velocity distribution function and is the Boltzmann collision operator [7, 8, 11, 13, 26, 48], whose explicit expression will not be written down here. The notation emphasizes the fact that, apart from its spatial and velocity dependencies, the distribution function depends on the independent external parameters and . As said above, and correspond to the Couette and Poiseuille flows, respectively.

In general, Eq. (1) must be solved subjected to specific boundary conditions, which can be expressed in terms of the kernels defined as follows. When a particle with velocity hits the wall at , the probability of being reemitted with a velocity within the range is ; the kernel represents the same but at . The boundary conditions are then [13]

 Θ(±vy)|vy|f(y=∓L/2,v)=Θ(±vy)∫dv′|v′y|K∓(v,v′)Θ(∓v′y)f(y=∓L/2,v′), (2)

where is Heaviside’s step function. In the case of boundary conditions of complete accommodation with the walls, so that does not depend on the incoming velocity , the kernels can be written as

 K∓(v)=A−1∓Θ(±vy)|vy|φ∓(v),A∓≡∫dvΘ(±vy)|vy|φ∓(v), (3)

where represents the probability distribution of a fictitious gas in contact with the system at . Equation (3) can then be interpreted as meaning that when a particle hits a wall, it is instantaneously absorbed and replaced by a particle leaving the fictitious bath. Of course, any choice of must be consistent with the imposed wall velocities and temperatures, i.e.

 U∓=∫dvvxφ∓(v), (4)
 kBTw=13m∫dv(v−U∓)2φ∓(v). (5)

Inserting Eq. (3) into Eq. (2), one has

 Θ(±vy)f(y=∓L/2,v)=Θ(±vy)n∓φ∓(v), (6)

where

 n∓≡∫dv′Θ(∓v′y)|v′y|f(y=∓L/2,v′)∫dv′Θ(±v′y)|v′y|φ∓(v′). (7)

The boundary conditions (6) are usually referred to as diffuse (or stochastic) boundary conditions. The simplest and most common choice is that of a Maxwell–Boltzmann (MB) distribution [52]:

 φMB∓(v)=(m2πkBTw)3/2exp[−m(v−U∓)22kBTw]. (8)

The first few moments of define the densities of conserved quantities (mass, momentum, and temperature) and the associated fluxes. More explicitly,

 n(y|ω,g)=∫dvf(y,v|ω,g), (9)
 n(y|ω,g)u(y|ω,g)=∫dvvf(y,v|ω,g), (10)
 n(y|ω,g)kBT(y|g)=p(y|ω,g)=m3∫dvV2(y,v|ω,g)f(y,v|ω,g), (11)
 Pij(y|ω,g)=m∫dvVi(y,v|ω,g)Vj(y,v|ω,g)f(y,v|ω,g), (12)
 q(y|ω,g)=m2∫dvV2(y,v|ω,g)V(y,v|ω,g)f(y,v|ω,g). (13)

In these equations is the number density, is the flow velocity,

 V(y,v|a,ω,)≡v−u(y|ω,g) (14)

is the peculiar velocity, is the temperature, is the Boltzmann constant, is the hydrostatic pressure, is the pressure tensor, and is the heat flux. Taking velocity moments in both sides of Eq. (1) one gets the following exact balance equations

 ∂yPyy=0, (15)
 ∂yPxy=mng, (16)
 ∂yqy+Pxy∂yux=0. (17)

Henceforth, without loss of generality, we will assume . In other words, we will adopt a reference frame solidary with the flow at the midpoint .

The symmetry properties of the Couette–Poiseuille flow imply the following invariance properties of the velocity distribution function:

 f(y,vx,vy,vz|ω,g) = f(−y,−vx,−vy,vz|ω,−g) (18) = f(−y,vx,−vy,vz|−ω,g) = f(y,vx,vy,−vz|ω,g),

As a consequence, if denotes a hydrodynamic variable or a flux, one has

 χ(y|ω,g) = Sgχ(−y|ω,−g) (19) = Sωχ(−y|−ω,g),

where and . The parity factors and for the non-zero hydrodynamic fields and fluxes are displayed in Table 1. In general, if is a moment of the form

 χ(y|ω,g)=∫dvVkxx(y|ω,g)vkyyv2kzzf(y,v|ω,g) (20)

then and .

The general solution to the stationary Boltzmann equation (1) with the boundary conditions (6) can be split into two parts [45, 46, 47]:

 f(y,v|ω,g)=fH(y,v|ω,g)+fB(y,v|ω,g). (21)

Here, represents the hydrodynamic, Hilbert-class, or normal contribution to the distribution function. This means that depends on space only through a functional dependence on the hydrodynamic fields, i.e.

 fH(y)=fH[n,T,ux]. (22)

The contribution represents the boundary-layer correction to , so that verifies the specified boundary conditions. The correction is appreciably different from zero only in a thin layer (the so-called boundary layer or Knudsen layer), adjacent to the plates, of thickness of the order of the mean free path. Consequently, if the separation between the plates is much larger than the characteristic mean free path, there exists a well defined bulk region where the boundary correction vanishes and the distribution function is fully given by its hydrodynamic part. In the boundary layers the hydrodynamic profiles are much less smooth than in the bulk domain. The values of the flow velocity near the walls are different from the velocity of the plates (velocity slip phenomenon), i.e. . Besides, the extrapolation of the velocity profile in the bulk to the boundaries, , is also different from both the actual values and the wall velocities . Of course, an analogous temperature jump effect takes place with the temperature profile. The boundary contribution for small Knudsen numbers has been analyzed elsewhere [8, 47].

In the remaining of this paper we will focus on the hydrodynamic part (and will drop the subscript H), with special emphasis on the corresponding hydrodynamic contributions to the momentum and heat fluxes. In order to nondimensionalize the problem, we choose quantities evaluated at the central plane as units:

 f∗(s,v∗|a,g∗)≡v3T(0)n(0)f(y,v|ω,g),v∗≡vvT(0),vT(0)≡√kBT(0)m, (23)
 n∗(s|a,g∗)≡n(y|ω,g)n(0),T∗(s|a,g∗)≡T(y|ω,g)T(0),p∗(s|a,g∗)≡p(y|ω,g)p(0), (24)
 P∗ij(s|a,g∗)≡Pij(y|ω,g)p(0),q∗(s|a,g∗)≡q(y|ω,g)p(0)vT(0), (25)
 a≡1ν(0)∂ux∂y∣∣∣y=0,g∗≡gvT(0)ν(0),ν∗≡νν(0). (26)

In the above equations we have found it convenient to introduce the dimensionless scaled spatial variable

 s(y)≡1vT(0)∫y0dy′ν(y′), (27)

where is an effective collision frequency. For the sake of concreteness, we choose it as

 ν(y)=p(y)η(y), (28)

where is the NS shear viscosity. The change from the boundary-imposed shear rate to the reduced local shear rate is motivated by our goal of focusing on the central bulk region of the system, outside the boundary layers. Note that represents the Knudsen number associated with the velocity gradient at . Likewise, measures the strength of the external field on a particle moving with the thermal velocity along a distance on the order of the mean free path.

The relationship (27) can be inverted to yield

 y∗(s)=∫s0ds′ν∗(s′),y∗≡yvT(0)/ν(0). (29)

The invariance properties (18) translate into

 f∗(s,v∗x,v∗y,v∗z|a,g∗) = f∗(−s,−v∗x,−v∗y,v∗z|a,−g∗) (30) = f∗(−s,v∗x,−v∗y,v∗z|−a,g∗) = f∗(−s,v∗x,v∗y,−v∗z|a,g∗).

Given the symmetry properties (30), we can restrict ourselves to and without loss of generality.

## 3 Navier–Stokes description

To gain some insight into the type of fields one can expect in the Couette–Poiseuille flow, it is instructive to analyze the solution provided by the NS level of description. In the geometry of the problem, the NS constitutive equations are

 Pxx=Pyy=Pzz=p, (31)
 Pxy=−η∂yux, (32)
 qx=0, (33)
 qy=−κ∂yT, (34)

where is the shear viscosity, as said above, and is the thermal conductivity. Inserting the NS approximate relations (31)–(34) into the exact conservation equations (15)–(17) one gets

 p=const, (35)
 (η∂y)2ux=−ηmng, (36)
 5kB2mPr(η∂y)2T=−(η∂yux)2, (37)

where is the Prandtl number. In dimensionless form, Eqs. (36) and (37) can be rewritten as

 ∂2su∗x(s)=−n∗(s)ν∗(s)g∗, (38)
 ∂2sT∗(s)=−2Pr5[∂su∗x(s)]2. (39)

For simplicity, let us assume that the particles are Maxwell molecules [7, 11, 63], so and . In that case, Eqs. (38) and (39) allow for an explicit solution:

 u∗x(s|a,g∗)=as−12g∗s2, (40)
 T∗(s|a,g∗)=1−Pr30s2(6a2−4ag∗s+g∗2s2). (41)

Here we have applied the Galilean choice and the symmetry property .

Equation (40) shows that, according to the NS approximation, the velocity field in the Couette–Poiseuille flow is simply the superposition of the (quasi) linear Couette profile and the (quasi) parabolic Poiseuille profile. In the case of the temperature field, however, apart from the (quasi) parabolic Couette profile and the (quasi) quartic Poiseuille profile, a (quasi) cubic coupling term is present. Here we use the term “quasi” because the simple polynomial forms in Eqs. (40) and (41) refer to the scaled variable . To go back to the real spatial coordinate one needs to make use of the relationship (27), taking into account that for Maxwell molecules . Instead of expressing as a function of it is more convenient to proceed in the opposite sense by using Eq. (29). Since one simply has

 y∗(s)=s[1−Pr30s2(2a2−ag∗s+15g∗2s2)]. (42)

For further use, note that, according to Eq. (41),

 ∂2T∗∂y∗2∣∣∣y∗=0=−Pr25a2. (43)

Thus, the NS temperature profile presents a maximum at the midpoint .

Before closing this section, let us write the pressure tensor and the heat flux profiles provided by the NS description:

 P∗xx(s|a,g∗)=P∗yy(s|a,g)=P∗zz(s|a,g)=1, (44)
 P∗xy(s|a,g∗)=−a+g∗s, (45)
 q∗x(s|a,g∗)=0, (46)
 q∗y(s|a,g∗)=s(a2−ag∗s+13g∗2s2). (47)

## 4 Kinetic theory description. Perturbation solution.

Now we want to get the hydrodynamic and flux profiles in the bulk domain of the system from a purely kinetic approach, i.e., without assuming a priori the applicability of the NS constitutive equations. To that end, instead of considering the detailed Boltzmann operator we will make use of the celebrated BGK kinetic model [5, 7, 66]. In the BGK model Eq. (1) is replaced by

 (vy∂∂y+g∂∂vx)f(y,v|ω,g)=−ν(y|ω,g)[f(y,v|ω,g)−M(y,v|ω,g)], (48)

where is the effective collision frequency defined by Eq. (28) and

 M(v)=n(m2πkBT)3/2exp(−mV22kBT) (49)

is the local equilibrium distribution function. In terms of the dimensionless variables introduced in Eqs. (23)–(27), Eq. (48) can be rewritten as

 (1+v∗y∂s)f∗(s,v∗|a,g∗)=M∗(s,v∗|a,g∗)−g∗ν∗(s|a,g∗)∂∂v∗xf∗(s,v∗|a,g∗). (50)

Its formal solution is

 f∗(v∗) = (1+v∗y∂s)−1[M∗(v∗)−g∗ν∗∂∂v∗xf∗(v∗)] (51) = ∞∑k=0(−v∗y∂s)k[M∗(v∗)−g∗ν∗∂∂v∗xf∗(v∗)].

The formal character of the solution (51) is due to the fact that appears on the right-hand side explicitly and also implicitly through and . The solvability (or consistency) conditions are

 Missing or unrecognized delimiter for \right (52)

Let us assume now that is a small parameter so the solution to Eq. (50) can be expanded as

 f∗(s,v∗|a,g∗)=f∗0(s,v∗|a)+f∗1(s,v∗|a)g∗+f∗2(s,v∗|a)g∗2+⋯. (53)

Likewise,

 χ∗(s|a,g∗)=χ∗0(s|a)+χ∗1(s|a)g∗+χ∗2(s|a)g∗2+⋯, (54)

where denotes a generic velocity moment of . The expansions of , , and induce the corresponding expansion of . The expansion in powers of allows the iterative solution of Eq. (51) by a scheme similar to that followed in Ref. [54] in the case of an external force normal to the plates.

### 4.1 Zeroth order in g. Pure Couette flow

#### Finite shear rates

To zeroth order in , Eqs. (50) and (51) become

 (1+v∗y∂s)f∗0(s,v∗|a)=M∗0(s,v∗|a), (55)
 f∗0(s,v∗|a)=∞∑k=0(−v∗y∂s)kM∗0(s,v∗|a), (56)

where

 M∗0(v∗)=p∗0(2π)3/2T∗05/2exp(−V∗022T∗0),V∗0≡v∗−u∗0. (57)

These are just the equations corresponding to the pure Couette flow. The complete solution has been obtained elsewhere [6, 20, 25] and so here we only quote the final results. The hydrodynamic profiles are

 p∗0(s|a)=1, (58)
 u∗x,0(s|a)=as, (59)
 T∗0(s|a)=1−γ(a)s2, (60)

where the dimensionless parameter is a nonlinear function of the reduced shear rate given implicitly through the equation [6, 20]

 a2=γ[3+2F2(γ)F1(γ)], (61)

where the mathematical functions are defined by

 F0(x)=2x∫∞0dtte−t2/2K0(2x−1/4t1/2),Fr(x)=(ddxx)rF0(x), (62)

being the zeroth-order modified Bessel function. Equation (62) clearly shows that has an essential singularity at and thus its expansion in powers of ,

 Fr(x)=∞∑k=0(k+1)r(2k+1)!(2k+1)!!(−x)k, (63)

is asymptotic and not convergent. However, the series representation (63) is Borel summable [6, 25], the corresponding integral representation being given by Eq. (62). The functions with can be easily expressed in terms of , , and as

 F3(x)=1−F0(x)8x−F2(x)−14F1(x), (64)
 Fr(x)=18xr−3∑m=0(r−3m)(−1)m+rFm(x)−Fr−1(x)−14Fr−2(x),r≥4. (65)

It is interesting to compare the hydrodynamic profiles with the results obtained from the Boltzmann equation at NS order (see Sec. 3). We observe that Eq. (58) agrees with Eq. (35) and Eq. (59) agrees with Eq. (40) for . On the other hand, Eq. (41) with differs from Eq. (60), except in the limit of small shear rates, in which case (Note that in the BGK model).

The relevant transport coefficients of the steady Couette flow are obtained from the pressure tensor and the heat flux. They are highly nonlinear functions of the reduced shear rate given by [6, 19, 20, 30]

 P∗xx,0(s|a)=1+4γ[F1(γ)+F2(γ)], (66)
 P∗yy,0(s|a)=1−2γ[F1(γ)+2F2(γ)], (67)
 P∗zz,0(s|a)=1−2γF1(γ), (68)
 P∗xy,0(s|a)=−aF0(γ), (69)
 q∗x,0(s|a)=a2{F0(γ)−1−10γF1(γ)−8γF2(γ)[1−F2(γ)F1(γ)]}s, (70)
 q∗y,0(s|a)=a2F0(γ)s. (71)

Notice that, although the temperature gradient is only directed along the axis (so that there is a response in this direction through ), the shear flow induces a nonzero component of the heat flux [19, 20, 30, 37]. Furthermore, normal stress differences (absent at NS order) are present. Equations (69) and (71) can be used to identify generalized nonlinear shear viscosity and thermal conductivity coefficients.

In general, the velocity moments of degree of are polynomial functions of the spatial variable of degree . An explicit expression for the velocity distribution function has also been derived [20, 25].

#### Limit of small shear rates

The coefficient characterizing the profile of the zeroth-order temperature is a complicated nonlinear function of the reduced shear rate , as clearly apparent from Eq. (61). Obviously, the zeroth-order pressure tensor and heat flux given by Eqs. (66)–(71) inherit this nonlinear character.

It is illustrative to assume that the reduced shear rate is small so one can express the quantities of interest as the first few terms in a (Chapman–Enskog) series expansion. From Eqs. (61)–(71) one obtains

 γ(a)=a25(1+7225a2+⋯), (72)
 P∗xx,0(s|a)=1+8a25(1−19825a2+⋯), (73)
 P∗yy,0(s|a)=1−6a25(1−22825a2+⋯), (74)
 P∗zz,0(s|a)=1−2a25(1−10825a2+⋯), (75)
 P∗xy,0(s|a)=−a(1−185a2+⋯), (76)
 q∗x,0(s|a)=−14a35(1−1836175a2+⋯)s, (77)
 q∗y,0(s|a)=a2(1−185a2+⋯)s. (78)

The terms of order , , and in Eqs. (72), (76), and (78), respectively, agree with the corresponding NS expressions, Eqs. (41), (45), and (47). On the other hand, as noted above, the normal stress differences ( and ) and the streamwise heat flux component reveal non-Newtonian effects of orders and , respectively.

### 4.2 First order in g. Couette–Poiseuille flow

#### Finite shear rates

To first order in Eq. (51) yields

 f∗1(v∗)−M∗1(v∗)=Λ(% I)(v∗)+Λ(II)(v∗), (79)

where

 Λ(I)(v∗)≡∞∑k=1(−v∗y∂s)kM∗1(v∗),Λ(II)(v∗)≡−∞∑k=0(−v∗y∂s)kT∗0∂∂v∗xf∗0(v∗), (80)
 M∗1(v∗)=M∗0(v∗)[p∗1+T∗12T∗0(V∗02T∗0−5)+u∗x,1T∗0V∗x,0], (81)

and we have already specialized to Maxwell molecules, so that . In order to apply the consistency conditions (52) in the derivation of the hydrodynamic fields , , and , it is convenient to define the moments

 Φn1n2n3=Φ(I)n1n2n3+Φ(II)n1n2n3, (82)
 Φ(I)n1n2n3=∫dv∗V∗x,0n1v∗yn2v∗zn3Λ(I)(v∗), (83)
 Φ(II)n1n2n3=∫dv∗V∗x,0n1v∗yn2v∗zn3Λ(II)(v∗). (84)

Therefore, the consistency conditions are

 Φ000=0, (85)
 Φ100=0, (86)
 Φ010=0, (87)
 Φ200+Φ020+Φ002=0. (88)

The evaluation of and is carried out in Appendix A. The first-order profiles are

 p∗1(s|a)=p(1)1(a)s, (89)
 u∗x,1(s|a)=u(2)x,1(a)s2, (90)
 T∗1(s|a)=T(3)1(a)s3, (91)

where

 p(1)1(a)=−1γ(1−F1F0)T(3)1(a), (92)
 u(2)x,1(a)=2γ(F1+2F2)−36F1−aγ(F1F0−F2F1)T(3)1(a), (93)
 T(3)1(a)=43aγ2F04γF2(F1+2F2)−F1(1−F0)−6F2D(a), (94)

with

 D(a) ≡ 2γF1[F20−2F0+F1−2γ(F0−2F1)(F1+2F2)]+a2F0F1(1−F0) (95) −2a2γ[F0(F21+6F1F2+8F22)−2F21(F1+2F2)].

In the above equations the functions are understood to be evaluated at .

As shown in Appendix A, the moment is a polynomial function of of degree . In particular, the non-zero elements of the first-order pressure tensor are

 P∗ij,1(s|a)=P(1)ij,1(a)s, (96)

with

 P(1)xx,1(a)=3p(1)1(a)−P(1)zz,1(a), (97)
 P(1)yy,1(a)=0, (98)
 P(1)zz,1(a)=−(F0−12γ+F1+2F2)T(3)1−[2γ(F1