2.1. Linearisation and non-dimensionalisation of the Boltzmann equation

In the original Boltzmann equation, the distribution function $f(\boldsymbol {x}, \boldsymbol {\xi }, t)$ satisfies the following integro-differential equation:

(2.1) \begin{equation} \frac {\partial f}{\partial t} + \xi _k \frac {\partial f}{\partial x_k} = Q[f,f] := \int _{\mathbb {R}^3} \int _{\mathbb {S}^2} B(\boldsymbol {\xi } - \boldsymbol {\xi }_1, \boldsymbol {\omega }) [f(\boldsymbol {\xi }_1^{\prime}) f(\boldsymbol {\xi }') - f(\boldsymbol {\xi }_1) f(\boldsymbol {\xi })] \,\mathrm {d}\boldsymbol {\omega } \,\mathrm {d}\boldsymbol {\xi }_1. \end{equation}

Here $B(\cdot , \cdot )$ stands for the collision kernel, and

(2.2) \begin{equation} \boldsymbol {\xi }' = \frac {\boldsymbol {\xi } + \boldsymbol {\xi }_1}{2} + \frac {|\boldsymbol {\xi } - \boldsymbol {\xi }_1|}{2} \boldsymbol {\omega }, \qquad \boldsymbol {\xi }_1' = \frac {\boldsymbol {\xi } + \boldsymbol {\xi }_1}{2} - \frac {|\boldsymbol {\xi } - \boldsymbol {\xi }_1|}{2} \boldsymbol {\omega }. \end{equation}

In this work we are interested in the linear regime where the distribution function is close to a global Maxwellian with density $\rho _0$ , velocity $\boldsymbol {v}_0$ , temperature $\theta _0$ and mass of a single molecule $m$ :

(2.3) \begin{equation} \mathcal {M}_0(\boldsymbol {\xi }) = \frac {\rho _0}{m(2\pi \theta _0)^{3/2}} \exp \left ( -\frac {|\boldsymbol {\xi } - \boldsymbol {v}_0|^2}{2\theta _0} \right ). \end{equation}

We can then introduce a small quantity $\epsilon$ and write the distribution function $f$ as

(2.4) \begin{equation} f(\boldsymbol {x},\boldsymbol {\xi },t) = \mathcal {M}_0(\boldsymbol {\xi }) [1 + \epsilon \hat {f}(\hat {\boldsymbol {x}},\hat {\boldsymbol {\xi }},\hat {t})], \end{equation}

where $\hat {\boldsymbol {x}}$ , $\hat {\boldsymbol {\xi }}$ and $\hat {t}$ are dimensionless variables defined by

(2.5) \begin{equation} \hat {\boldsymbol {x}} = (\boldsymbol {x} - \boldsymbol {v}_0 t) / L, \qquad \hat {\boldsymbol {\xi }} = (\boldsymbol {\xi } - \boldsymbol {v}_0) / \sqrt {\theta _0}, \qquad \hat {t} = t / (L / \sqrt {\theta _0}). \end{equation}

Plugging (2.4) into the Boltzmann equation yields

(2.6) \begin{equation} \frac {\partial \hat {f}}{\partial \hat {t}} + \hat {\xi }_k \frac {\partial \hat {f}}{\partial \hat {x}_k} = \frac {\rho _0 L}{m \sqrt {\theta _0}} \frac {\mathcal {Q}[\hat {\mathcal {M}}_0,\hat {f} \hat {\mathcal {M}}_0] + \mathcal {Q}[\hat {f}\hat {\mathcal {M}}_0, \hat {\mathcal {M}}_0] + \epsilon \mathcal {Q}[\hat {f}\hat {\mathcal {M}}_0, \hat {f} \hat {\mathcal {M}}_0]}{\hat {\mathcal {M}}_0}, \end{equation}

where $\hat {\mathcal {M}}_0$ is the dimensionless Maxwellian distribution function

(2.7) \begin{equation} \hat {\mathcal {M}}_0(\hat {\boldsymbol {\xi }}) = \frac {1}{(2\pi )^{3/2}} \exp \left ( -\frac {|\hat {\boldsymbol {\xi }}|^2}{2} \right ). \end{equation}

The linearised equation can then be obtained by discarding the $O(\epsilon )$ term in (2.6) and defining

(2.8) \begin{equation} \hat {\mathcal {L}}[\hat {f}] := \frac {\mathcal {Q}[\hat {\mathcal {M}}_0,\hat {f} \hat {\mathcal {M}}_0] + \mathcal {Q}[\hat {f}\hat {\mathcal {M}}_0, \hat {\mathcal {M}}_0]}{B_0 \hat {\mathcal {M}}_0}, \qquad \textit{Kn} = \frac {m B_0 \sqrt {\theta _0}}{\rho _0 L}. \end{equation}

Here $B_0$ is a constant with the same dimension as the collision kernel $B(\cdot , \cdot )$ to guarantee that $\hat {\mathcal {L}}$ is dimensionless, and here we choose $B_0$ such that

(2.9) \begin{equation} \langle \hat {\xi }_1 \hat {\xi }_2, \hat {\mathcal {L}}(\hat {\xi }_1 \hat {\xi }_2) \rangle = -1, \end{equation}

where the inner product is defined by

(2.10) \begin{equation} \langle \hat {f}, \hat {g} \rangle = \int _{\mathbb {R}^3} \hat {f}(\hat {\boldsymbol {\xi }}) \hat {g}(\hat {\boldsymbol {\xi }}) \hat {\mathcal {M}}_0(\hat {\boldsymbol {\xi }}) \,\mathrm {d} \hat {\boldsymbol {\xi }}. \end{equation}

Note that the linearised Boltzmann equation satisfies the linearised H-theorem. In the original Boltzmann equation (2.1), the entropy density and the entropy flux are defined by

(2.11) \begin{equation} \eta = -k_B \int _{\mathbb {R}^3} f \log f \,\mathrm {d}\boldsymbol {\xi }, \qquad \psi _k = -k_B \int _{\mathbb {R}^3} \xi _k f \log f \,\mathrm {d}\boldsymbol {\xi }, \end{equation}

which satisfy the inequality

(2.12) \begin{equation} \frac {\partial \eta }{\partial t} + \frac {\partial \psi _k}{\partial x_k} \geqslant 0, \end{equation}

indicating the second law of thermodynamics. For the linearised Boltzmann equation, the corresponding definitions are

(2.13) \begin{equation} \hat {\eta } = -\int _{\mathbb {R}^3} |\hat {f}|^2 \,\mathrm {d}\hat {\boldsymbol {\xi }}, \qquad \hat {\psi }_k = -\int _{\mathbb {R}^3} \xi _k |\hat {f}|^2 \,\mathrm {d}\hat {\boldsymbol {\xi }}. \end{equation}

Since the linearised collision operator $\hat {\mathcal {L}}$ satisfies $\langle \hat {f}, \hat {\mathcal {L}}[\hat {f}] \rangle \leqslant 0$ for any $\hat {f}$ , the H-theorem turns out to be the $L^2$ stability of the linearised Boltzmann equation

(2.14) \begin{equation} \frac {\partial \hat {\eta }}{\partial \hat {t}} + \frac {\partial \hat {\psi }_k}{\partial \hat {x}_k} \geqslant 0. \end{equation}

Later in our derivation of moment equations, this property will be preserved.

2.2. Linearisation and non-dimensionalisation of the conservation laws

The laws of mass, momentum and energy conservation can be derived by taking moments of the Boltzmann equation, which yields

(2.15) \begin{align} \frac {\partial \rho }{\partial t} + \frac {\partial (\rho v_j)}{\partial x_j} &= 0, \nonumber \\ \frac {\partial (\rho v_i)}{\partial t} + \frac {\partial }{\partial x_j} (\rho v_i v_j + \rho \theta \delta _{ij} + \sigma _{ij}) &= 0, \nonumber \\ \frac {\partial }{\partial t} \left (\frac {1}{2} \rho v_i v_i + \frac {3}{2} \rho \theta \right ) + \frac {\partial }{\partial x_j} \left ( \frac {5}{2} \rho \theta v_j + \frac {1}{2} \rho v_i v_i v_j + \sigma _{ij} v_i + q_j \right ) &= 0. \end{align}

Here $\sigma _{ij}$ denotes the stress tensor and $q_j$ denotes the heat flux. To linearise these equations, we again assume that the fluid state is close to the equilibrium specified by density $\rho _0$ , velocity $v_{0,i}$ and temperature $\theta _0$ , and introduce dimensionless variables $\hat {\rho }$ , $\hat {\theta }$ , $\hat {v}_i$ , $\hat {\sigma }_{ij}$ and $\hat {q}_{i}$ by

(2.16) \begin{equation} \rho = \rho _0(1+\epsilon \hat {\rho }) , \!\!\quad\!\! \theta = \theta _0(1+\epsilon \hat {\theta }), \!\!\quad\!\! v_i = v_{0,i} + \epsilon \sqrt {\theta _0} \hat {v}_i, \!\!\quad\!\! \sigma _{ij} = \epsilon \rho _0 \theta _0 \hat {\sigma }_{ij}, \!\!\quad\!\! q_i = \epsilon \rho _0 \theta _0^{3/2}\hat {q}_i. \end{equation}

By substituting the above expressions together with (2.5) into the conservation laws, and then dropping all the higher-order terms, we reach the linearised and dimensionless conservation laws:

(2.17) \begin{equation} \begin{aligned} \frac {\partial \hat {\rho }}{\partial \hat {t}}+ \frac {\partial \hat {v}_j}{\partial \hat {x}_j} &= 0, \\ \frac {\partial \hat {v}_i}{\partial \hat {t}}+ \frac {\partial \hat {\rho }}{\partial \hat {x}_i} +\frac {\partial \hat {\theta }}{\partial \hat {x}_i}+ \frac {\partial \hat {\sigma }_{ij}}{\partial \hat {x}_j} &= 0, \\ \frac {\partial \hat {\theta }}{\partial \hat {t}}+ \frac {2}{3} \frac {\partial \hat {v}_j}{\partial \hat {x}_j} +\frac {2}{3}\frac {\partial \hat {q}_j}{\partial \hat {x}_j} &= 0. \end{aligned} \end{equation}

Another approach to derive the same equations is to directly take moments of the linearised Boltzmann equation. The relationship between the dimensionless moments ( $\hat {\rho }$ , $\hat {\theta }$ , etc.) and the dimensionless distribution function $\hat {f}$ is

(2.18) \begin{align} \hat {\rho } &= \langle 1, \hat {f} \rangle , \quad \hat {v}_i = \langle \hat {\xi }_i, \hat {f} \rangle , \quad \hat {\theta } = \left \langle \frac {1}{3} |\hat {\boldsymbol {\xi }}|^2 - 1, \hat {f} \right \rangle , \nonumber \\ \hat {\sigma }_{ij} &= \left \langle \hat {\xi }_i \hat {\xi }_j - \frac {1}{3} |\hat {\boldsymbol {\xi }}|^2 \delta _{ij}, \hat {f} \right \rangle , \qquad \hat {q}_i = \left \langle \frac {1}{2} |\hat {\boldsymbol {\xi }}|^2 \hat {\xi }_i, \hat {f} \right \rangle . \end{align}

Thus, using the fact that $\langle 1, \hat {\mathcal {L}} [\hat {f}] \rangle = \langle \hat {\xi }_i, \hat {\mathcal {L}} [\hat {f}] \rangle = \langle |\hat {\boldsymbol {\xi }}|^2, \hat {\mathcal {L}} [\hat {f}] \rangle = 0$ , one can take the inner product of the linearised equation (1.1) and the polynomials $1$ , $\hat {\xi }_i$ and $|\hat {\boldsymbol {\xi }}|^2$ to obtain the moment equations (2.17).

The two methods to formulate the conservation laws are essentially equivalent, since taking moments and linearising equations are two commutable operations. In our derivation of moment equations to be presented below, we follow the second approach that takes moments of the linearised equations (1.1). Due to the presence of $\hat {\sigma }_{ij}$ and $\hat {q}_i$ in (2.17), our purpose is to find a suitable closure that provides decent approximation to the linearised Boltzmann equation when the Knudsen number $\textit{Kn}$ is small. Hereafter, the accent ‘ $\hat {\cdot }$ ’ will no longer be added to dimensionless variables.

3.1. Choice of variables

The 13-moment equations were first proposed by Grad (Reference Grad1949), where the distribution function was expanded into an infinite series using Hermite polynomials. The coefficients of the polynomials were regarded as moments, and Grad chose to include the first 13 moments, namely the density ( $\rho$ ), momentum ( $\rho v_i$ ), pressure ( $\rho \theta$ ), stress tensor ( $\sigma _{ij}$ ) and heat flux ( $q_i$ ), into his equations. In later developments of 13-moment systems, such a choice of the 13 moments becomes a standard and is followed by most works (Myong Reference Myong1999; Struchtrup & Torrilhon Reference Struchtrup and Torrilhon2003; Zhu et al. Reference Zhu, Hong, Yang and Yong2015). In particular, for the regularised 13-moment equations for Maxwell molecules, such a choice of moments matches exactly with the analysis using the order of magnitude method (Struchtrup(Reference Struchtrup2005b , Chapter 8). It is seen that for linear equations, these moments, denoted by $\langle \mathbf {\phi }_{13}, f\rangle$ with $\mathbf {\phi }_{13}$ being a 13-dimensional vector of polynomials of $\boldsymbol {\xi }$ , fully describes the distribution function $f$ up to order $O(\textit{Kn})$ , which means for any polynomial $r(\boldsymbol {\xi })$ satisfying $\langle \mathbf {\phi }_{13}, r \rangle = 0$ , the moment $\langle r, f\rangle$ has at least order $O(\textit{Kn}^2)$ when performing asymptotic analysis. This property leads to the super-Burnett order of R13 equations for Maxwell molecules, meaning that the super-Burnett equations can be derived from the R13 equations by Chapman–Enskog expansion. However, it is also known that this no longer holds for other types of molecules (Struchtrup Reference Struchtrup2005a ). Inspired by this fact, in our study of 13-moment equations for general gas molecules, we adopt a different set of variables $\langle \tilde {\mathbf {\phi }}_{13}, f\rangle$ , such that we can again obtain $\langle r, f\rangle \sim o(\textit{Kn})$ once $\langle \tilde {\mathbf {\phi }}_{13}, r\rangle = 0$ . This will help achieve the super-Burnett order for the moment equations.

The details on the construction of $\tilde {\mathbf {\phi }}_{13}$ will be discussed in § 4. Here we only provide some brief information necessary for the formulation of the moment equations. Our 13 moments also include conservative quantities such as the density ( $\rho$ ), velocity ( $v_i$ ) and temperature ( $\theta$ ). Note that the velocity and the temperature are conservative only under linear settings. For the remaining 8 moments, they also include a trace-free 2-tensor and a vector like in the classical case. For consistency, we denote the 2-tensor by $\bar {\sigma }_{ij}$ and the vector by $\bar {q}_i$ . The relationship between these variables and the stress tensor ( $\sigma _{ij}$ ) and the heat flux ( $q_i$ ) will be given later in this section. The conservative variables $\rho$ , $v_i$ and $\theta$ are $O(\textit{Kn}^0)$ moments, while the higher-order moments $\bar {\sigma }_{ij}$ and $\bar {q}_i$ are $O(\textit{Kn})$ moments. Note that the conservative moments are always defined using the distribution function $f$ by

(3.1) \begin{equation} (\rho , v_i, \theta ) = \left \langle \left (1, \xi _i, \frac {1}{3} \xi _k \xi _k \right ), f \right \rangle , \end{equation}

where the definitions of $\bar {\sigma }_{ij}$ and $\bar {q}_i$ depend on the collision model. When Maxwell molecules are considered, the definitions of $\bar {\sigma }_{ij}$ and $\bar {q}_i$ reduce to those of $\sigma _{ij}$ and $q_i$ .

3.2. Moment equations and boundary conditions

In this section wel list the equations and boundary conditions derived for the aforementioned variables. Like the generalised 13-moment equations derived in Struchtrup (Reference Struchtrup2005a ), our equations will contain coefficients that depend on the type of gas molecules. In the equations below, these constants will be denoted by $k_i$ and $l_i$ . The equations are

(3.2) \begin{align} &\frac {\partial \rho }{\partial t}+ \nabla \cdot \boldsymbol {v} = 0, \end{align}

(3.3) \begin{align} &\frac {\partial \theta }{\partial t}+ \frac {2}{3} \nabla \cdot \boldsymbol {v} +\frac {2}{3}k_0\nabla \cdot \bar {\boldsymbol{q}}-k_1 \textit{Kn} \Delta \theta +k_2 \textit{Kn} \nabla \cdot (\nabla \cdot \bar {\boldsymbol {\sigma }} ) = 0 , \end{align}

(3.4) \begin{align} &\frac {\partial \boldsymbol {v}}{\partial t}+ \nabla \rho +\nabla \theta -k_3 \textit{Kn} \nabla \cdot (\nabla \boldsymbol {v})_{\textit{stf}} -k_4 \textit{Kn} \nabla \cdot (\nabla \bar{\boldsymbol{q}})_{\textit{stf}}+ k_5 \nabla \cdot \bar {\boldsymbol {\sigma }} = 0, \end{align}

(3.5) \begin{align} &\frac {\partial \bar{\boldsymbol{q}}}{\partial t}+ \frac {5}{2} k_0 \nabla \theta - \frac {5}{2}k_4 \textit{Kn}\nabla \cdot (\nabla \boldsymbol {v})_{\textit{stf}} -2 k_6 \textit{Kn} \nabla (\nabla \cdot \bar{\boldsymbol{q}}) \notag \\ & \quad {} -\frac {12}{5}k_7 \textit{Kn}\nabla \cdot (\nabla \bar{\boldsymbol{q}})_{\textit{stf}}+k_8 \nabla \cdot \bar {\boldsymbol {\sigma }} = - \frac {1}{\textit{Kn}} \frac {2}{3} l_1 \bar{\boldsymbol{q}}, \end{align}

(3.6) \begin{align} &\frac {\partial \bar {\boldsymbol {\sigma }}}{\partial t} + 3k_2 \textit{Kn} \left (\nabla ^2 \theta \right )_{\textit{stf}}+ 2k_5 (\nabla \boldsymbol {v})_{\textit{stf}} +\frac {4}{5}k_8 (\nabla \bar{\boldsymbol{q}})_{\textit{stf}} \notag \\ & \quad {}- 2k_9 \textit{Kn} \nabla \cdot (\nabla \bar {\boldsymbol {\sigma }})_{\textit{stf}} -k_{10} \textit{Kn} \left (\nabla (\nabla \cdot \bar {\boldsymbol {\sigma }}) \right )_{\textit{stf}}= - \frac {1}{\textit{Kn}} l_2 \bar {\boldsymbol {\sigma }}. \end{align}

Here $(\cdot )_{\textit{stf}}$ refers to the symmetric and trace-free part of a given tensor, which is defined entry wisely as $(\cdot )_{ij} \mapsto ((\cdot )_{\textit{stf}})_{ij} = (\cdot )_{\langle ij\rangle }$ for a matrix and $(\cdot )_{ijk} \mapsto ((\cdot )_{\textit{stf}})_{ijk} = (\cdot )_{\langle ijk\rangle }$ for a 3-tensor. The detailed definition can be found in Struchtrup (Reference Struchtrup2005b ). The coefficients $k_i$ and $l_i$ have been computed for gas molecules whose intermolecular potential satisfies inverse power laws. Values for these coefficients for some power indices are given in Appendix B. Although these coefficients look arbitrary, our derivation in § 4 will reveal that they are fully determined by the linearised and non-dimensionalised Boltzmann collision operator, which depends only on the potential energy between two gas molecules. In inverse-power-law models, only two parameters exist: the intensity coefficient of the potential and the power $\eta$ . After non-dimensionalisation, the intensity coefficient is integrated into the Knudsen number and, thus, $\eta$ is the only parameter that dictates all the coefficients.

The first three equations (3.2), (3.3) and (3.4) are conservation laws of mass, energy and momentum. By comparing (3.3) and (3.4) with the conservation laws expressed by the stress tensor and the heat flux (see (2.17)), i.e.

(3.7) \begin{equation} \frac {\partial \theta }{\partial t}+ \frac {2}{3} \nabla \cdot \boldsymbol {v} +\frac {2}{3}\nabla \cdot \boldsymbol{q}= 0, \qquad \frac {\partial \boldsymbol {v}}{\partial t}+\nabla \rho +\nabla \theta + \nabla \cdot \boldsymbol{\sigma } = 0, \end{equation}

one can observe the relationship between $ \bar {\boldsymbol {\sigma }}$ , $ \bar{\boldsymbol{q}}$ and $\boldsymbol{\sigma }$ , $\boldsymbol{q}$ :

(3.8) \begin{equation} \boldsymbol{\sigma } = k_5 \bar {\boldsymbol {\sigma }}-k_4 \textit{Kn} (\nabla \bar{\boldsymbol{q}})_{\textit{stf}} -k_3 \textit{Kn}(\nabla \boldsymbol {v})_{\textit{stf}}, \quad \boldsymbol{q} = k_0 \bar{\boldsymbol{q}}-\frac {3}{2}k_1 \textit{Kn} \nabla \theta +\frac {3}{2}k_2 \textit{Kn}\nabla \cdot \bar {\boldsymbol {\sigma }}. \end{equation}

Such a relationship allows us to calculate the stress tensor and heat flux once (3.2)–(3.6) are solved. We remark that the moment equations (3.2)–(3.6) we derived are Galilean invariant. The detailed proof can be found in Appendix A.

Another important ingredient of the moment equations is the wall boundary conditions, since the linear moment equations are often applied in channel flows, where the boundaries play crucial roles in rarefaction effects. The boundary conditions we present here are derived from the Maxwell boundary conditions for the Boltzmann equation, where the reflected gas flow is a linear combination of specular and diffusive reflections. The accommodation coefficient, denoted by $\chi$ , is the parameter describing the proportion of diffusive reflection. For simplicity, we are going to use $\tilde {\chi } = 2\chi / (2-\chi )$ in our formulation of boundary conditions.

We now focus on a specific point on the wall and let $\boldsymbol{n}$ be the outer unit normal vector at this point. Meanwhile, we let $\mathbf {\tau }_1$ and $\mathbf {\tau }_2$ be the two tangential unit vectors that are perpendicular to each other. When describing boundary conditions, we temporarily rotate our coordinate system such that the three axes become $\boldsymbol{n}$ , $\mathbf {\tau }_1$ and $\mathbf {\tau }_2$ . Correspondingly, the three components of the velocity will be denoted by $v_n$ , $v_{\tau _1}$ and $v_{\tau _2}$ (similar for $\bar {q}$ ), and the components of $\bar {\sigma }$ are written by $\bar {\sigma }_{nn}, \bar {\sigma }_{n \tau _1}, \bar {\sigma }_{n \tau _2}, \bar {\sigma }_{\tau _1 \tau _1}, \bar {\sigma }_{\tau _1 \tau _2}, \bar {\sigma }_{\tau _2 \tau _2}$ . Moreover, we assume that the temperature and velocity of the wall are $\theta ^W$ and $v_i^W$ at this point. Under these assumptions, the wall boundary conditions have the following form:

(3.9) \begin{align} & v_n = 0, \end{align}

(3.10) \begin{align} & \bar {q}_{n} = \tilde {\chi } \left [ m_{11} (\theta -\theta ^{W}) +m_{12} \bar {\sigma }_{nn} - m_{13} \textit{Kn} \frac {\partial \bar {q}_j}{\partial x_j} - m_{14} \textit{Kn} \frac {\partial \bar {q}_{\langle n}}{\partial x_{n \rangle }} - m_{15} \textit{Kn} \frac {\partial v_{\langle n}}{\partial x_{n \rangle }} \right ], \end{align}

(3.11) \begin{align} m_{26}\bar {q}_{n} + m_{27} \textit{Kn} \frac {\partial \theta }{\partial x_n} - m_{28} \textit{Kn} \frac {\partial \bar {\sigma }_{nj}}{\partial x_j} &=\tilde {\chi } \left [-m_{21} (\theta -\theta ^{W}) + m_{22} \bar {\sigma }_{nn} + m_{23} \textit{Kn} \frac {\partial \bar {q}_j}{\partial x_j}\notag\right. \\ & \qquad \left. + m_{24} \textit{Kn} \frac {\partial \bar {q}_{\langle n}}{\partial x_{n \rangle }} + m_{25} \textit{Kn} \frac {\partial v_{\langle n}}{\partial x_{n \rangle }} \right ], \end{align}

(3.12) \begin{align} \bar {\sigma }_{\tau _i n} &= \tilde {\chi } \left [m_{31} (v_{\tau _i} -v^{W}_{\tau _i}) +m_{32}\bar {q}_{\tau _i} - m_{33} \textit{Kn} \frac {\partial \bar {\sigma }_{\tau _i j}}{\partial x_j} - m_{34} \textit{Kn} \frac {\partial \bar {\sigma }_{\langle \tau _i n}}{\partial x_{n \rangle }} + m_{35} \textit{Kn} \frac {\partial \theta }{\partial x_{\tau _i}} \right ], \notag \\ i&=1,2, \end{align}

(3.13) \begin{align} m_{46}\bar {\sigma }_{\tau _in} &+ m_{47} \textit{Kn}\frac {\partial v_{\langle \tau _i}}{\partial x_{n\rangle }} + m_{48} \textit{Kn}\frac {\partial \bar {q}_{\langle \tau _i}}{\partial x_{n\rangle }} = - \tilde {\chi } \left [-m_{41} (v_{\tau _i} -v^{W}_{\tau _i}) +m_{42}\bar {q}_{\tau _i}\notag \right.\\ & \left. +\, m_{43} \textit{Kn} \frac {\partial \bar {\sigma }_{\tau _i j}}{\partial x_j} + m_{44} \textit{Kn} \frac {\partial \bar {\sigma }_{\langle \tau _i n}}{\partial x_{n \rangle }} - m_{45} \textit{Kn} \frac {\partial \theta }{\partial x_{\tau _i}} \right ],\quad i=1,2, \end{align}

(3.14) \begin{align} m_{56}\bar {\sigma }_{\tau _in} & + m_{57} \textit{Kn}\frac {\partial v_{\langle \tau _i}}{\partial x_{n\rangle }} + m_{58} \textit{Kn}\frac {\partial \bar {q}_{\langle \tau _i}}{\partial x_{n\rangle }} = - \tilde {\chi } \left [-m_{51} (v_{\tau _i} -v^{W}_{\tau _i}) +m_{52}\bar {q}_{\tau _i} \notag \right.\\ & \left.+\, m_{53} \textit{Kn} \frac {\partial \bar {\sigma }_{\tau _i j}}{\partial x_j} + m_{54} \textit{Kn} \frac {\partial \bar {\sigma }_{\langle \tau _i n}}{\partial x_{n \rangle }} - m_{55} \textit{Kn} \frac {\partial \theta }{\partial x_{\tau _i}} \right ],\quad i=1,2, \end{align}

(3.15) \begin{align} m_{66} \bar {q}_n & + \ \textit{Kn}\left (m_{67} \frac {\partial \bar {\sigma }_{\langle nn }}{\partial x_{n \rangle }} +m_{68}\frac {\partial \bar {\sigma }_{nj }}{\partial x_{j}} + m_{69}\frac {\partial \theta }{\partial x_n} \right ) = - \tilde {\chi } \left [-m_{61} (\theta -\theta ^{W}) +m_{62} \bar {\sigma }_{nn}\notag \right.\\ & \left. + \, m_{63} \textit{Kn} \frac {\partial \bar {q}_j}{\partial x_j} + m_{64} \textit{Kn} \frac {\partial \bar {q}_{\langle n}}{\partial x_{n \rangle }} + m_{65} \textit{Kn} \frac {\partial v_{\langle n}}{\partial x_{n \rangle }}\right ], \end{align}

(3.16) \begin{align} & \textit{Kn}\left (\frac {\partial \bar {\sigma }_{\langle \tau _i \tau _i}}{\partial x_{n \rangle }} + \frac {1}{2} \frac {\partial \bar {\sigma }_{\langle n n}}{\partial x_{n \rangle }} \right ) = -\tilde {\chi } m_{71} \left ( \bar {\sigma }_{\tau _i \tau _i} + \frac {1}{2} \bar {\sigma }_{nn} \right ), \quad i=1,2, \end{align}

(3.17) \begin{align} & \textit{Kn} \frac {\partial \bar {\sigma }_{\langle \tau _1 \tau _2}}{\partial x_{n\rangle }} = -\tilde {\chi } m_{81} \bar {\sigma }_{\tau _1 \tau _2}. \end{align}

In these equations, $m_{\!jk}$ are constants depending on the molecular interactions. As examples, we consider inverse-power-law models in which the repulsive force between two gas molecules is proportional to the $\eta$ th power of the distance between them. For such models, the values of the constants $k_i$ and $m_{ij}$ are tabulated in Appendix B for some power indices.

It will become clear by its derivation that this set of boundary conditions satisfies the general form of $L^2$ -stable boundary conditions (1.6), which implies the second law of thermodynamics. The external contribution $\boldsymbol {g}_{\textit{ext}}$ is given by the terms with wall velocity $v_i^W$ and wall temperature $\theta ^W$ .

3.3. Second law of thermodynamics

The derivation of the new R13 equations (3.2)—(3.6), to be detailed in the next section, maintains the structure required by (1.2). Therefore, according to (1.3), the second law of thermodynamics is automatically satisfied. In fact, the entropy inequality can also be observed from (3.2) to (3.6) directly by straightforward calculations. Assume that the spatial domain is $\mathbb {R}^3$ and all quantities decay to zero when $\boldsymbol {x}$ tends to infinity. Then, by calculating

(3.18) \begin{equation} \int _{\mathbb {R}^3} \left [\rho \cdot (3.1) + \frac {3}{2} \theta \cdot (3.2) + \boldsymbol {v} \cdot (3.3) + \frac {2}{5} \bar{\boldsymbol{q}} \cdot (3.4) + \frac {1}{2} \bar {\boldsymbol {\sigma }} : (3.5) \right ] \mathrm {d}\boldsymbol {x}, \end{equation}

one obtains

(3.19) \begin{align} & \frac {\mathrm {d}}{\mathrm {d}t} \int _{\mathbb {R}^3} \frac {1}{2} \left (\rho ^2 + \frac {3}{2} \theta ^2 + |\boldsymbol {v}|^2 + \frac {2}{5} | \bar{\boldsymbol{q}}|^2 + \frac {1}{2} \| \bar {\boldsymbol {\sigma }}\|^2 \right ) \mathrm {d}\boldsymbol {x} = \notag\\& -\ \textit{Kn}\int _{\mathbb {R}^3} \left ( \frac {3}{2} k_1 |\nabla \theta |^2 + k_3 \|(\nabla \boldsymbol {v})_{\textit{stf}}\|^2 + \frac {4}{5} k_6 \|\nabla \bar{\boldsymbol{q}}\|^2 + \frac {24}{25} k_7 \|(\nabla \bar{\boldsymbol{q}})_{\textit{stf}}\|^2 \right ) \mathrm {d}\boldsymbol {x} \notag \\ & -\ \textit{Kn}\int _{\mathbb {R}^3} \left ( k_9 {\vert \kern -0.25ex \vert \kern -0.25ex \vert (\nabla \bar {\boldsymbol {\sigma }})_{\textit{stf}} \vert \kern -0.25ex \vert \kern -0.25ex \vert }^2 + \frac {1}{2} k_{10} |\nabla \cdot \bar {\boldsymbol {\sigma }}|^2\right ) \mathrm {d}\boldsymbol {x} - \frac {1}{ \textit{Kn}}\int _{\mathbb {R}^3} \left ( \frac {4}{15} l_1 | \bar{\boldsymbol{q}}|^2 + \frac {1}{2} l_2 \| \bar {\boldsymbol {\sigma }}\|^2 \right ) \mathrm {d}\boldsymbol {x} \notag \\ & + 3\ \textit{Kn} k_2\int _{\mathbb {R}^3} \nabla \theta \cdot \left ( \nabla \cdot \bar {\boldsymbol {\sigma }} \right ) \mathrm {d}\boldsymbol {x} - 2 \ \textit{Kn} k_4 \int _{\mathbb {R}^3} \left ( \nabla \boldsymbol {v} \right )_{\textit{stf}} : \left ( \nabla \bar{\boldsymbol{q}} \right )_{\textit{stf}} \mathrm {d}\boldsymbol {x}, \end{align}

where we have used $|\cdot |$ , $\|\cdot \|$ and $\vert \kern -0.25ex \vert \kern -0.25ex \vert \cdot \vert \kern -0.25ex\vert \kern -0.25ex\vert$ to denote the Frobenius norms of vectors, 2-tensors and 3-tensors, respectively. We show later (at the end of § 4.3) that

(3.20) \begin{equation} \frac {(3k_2)^2}{4} \leqslant \frac {3}{2}k_1 \cdot \frac {1}{2} k_{10}, \qquad k_4^2 \leqslant k_3 \cdot \frac {24}{25} k_7, \end{equation}

which yields

(3.21) \begin{align} \left |3k_2\int _{\mathbb {R}^3}\nabla \theta \cdot (\nabla \cdot \bar {\boldsymbol {\sigma }}) \,\mathrm {d}\boldsymbol {x}\right | \leqslant \int _{\mathbb {R}^3} \left ( \frac {3}{2} k_1 |\nabla \theta |^2 + \frac {1}{2} k_{10} |\nabla \cdot \bar {\boldsymbol {\sigma }}|^2\right )\mathrm {d}\boldsymbol {x}, \nonumber \\ \left |2k_4\int _{\mathbb {R}^3}(\nabla \boldsymbol {v})_{\textit{stf}} : (\nabla \bar{\boldsymbol{q}})_{\textit{stf}} \,\mathrm {d}\boldsymbol {x}\right | \leqslant \int _{\mathbb {R}^3} \left ( k_3 \|(\nabla \boldsymbol {v})_{\textit{stf}}\|^2 + \frac {24}{25} k_7 \|(\nabla \bar{\boldsymbol{q}})_{\textit{stf}} \|^2\right )\mathrm {d}\boldsymbol {x}. \end{align}

Thus,

(3.22) \begin{align} & \frac {\mathrm {d}}{\mathrm {d}t} \int _{\mathbb {R}^3} \frac {1}{2} \left (\rho ^2 + \frac {3}{2} \theta ^2 + |\boldsymbol {v}|^2 + \frac {2}{5} | \bar{\boldsymbol{q}}|^2 + \frac {1}{2} \| \bar {\boldsymbol {\sigma }}\|^2 \right ) \mathrm {d}\boldsymbol {x}\leqslant \nonumber \\ & {-\,\ \textit{Kn}} \int _{\mathbb {R}^3} \left ( \frac {4}{5} k_6 \|\nabla \bar{\boldsymbol{q}}\|^2 + k_9 { \vert \kern -0.25ex \vert \kern -0.25ex \vert (\nabla \bar {\boldsymbol {\sigma }})_{\textit{stf}} \vert \kern -0.25ex \vert \kern -0.25ex \vert }^2 \right ) \mathrm {d}\boldsymbol {x} - \frac {1}{ \textit{Kn}}\int _{\mathbb {R}^3} \left ( \frac {4}{15} l_1 | \bar{\boldsymbol{q}}|^2 + \frac {1}{2} l_2 \| \bar {\boldsymbol {\sigma }}\|^2 \right ) \mathrm {d}\boldsymbol {x} , \end{align}

Since the coefficients $k_6$ , $k_9$ , $l_1$ and $l_2$ are all non-negative, the right-hand side of the equation above is non-positive. The entropy density can thus be defined by

(3.23) \begin{equation} H = H_0 - \frac {1}{2} \left ( \rho ^2 + \frac {3}{2} \theta ^2 + |\boldsymbol {v}|^2 + \frac {2}{5} | \bar{\boldsymbol{q}}|^2 + \frac {1}{2} \| \bar {\boldsymbol {\sigma }}\|^2 \right ), \end{equation}

where $H_0$ is an arbitrary constant. In the case of Maxwell molecules, this result echos the conclusion in Struchtrup & Torrilhon (Reference Struchtrup and Torrilhon2007).

For bounded domains, with boundary conditions (3.9)–(3.17), we can also show that the entropy may increase only due to the incoming fluxes from the boundary. The derivation is similar to the unbounded case, while boundary terms appear when performing integration by parts, for which the boundary conditions need to be applied to make further simplifications. For any bounded domain $\Omega$ , the conclusion has the form

(3.24) \begin{align} -\frac {\mathrm {d}}{\mathrm {d}t} \int _{\Omega } H \,\mathrm {d}\boldsymbol {x} \leqslant {} & \int _{\partial \Omega } \theta ^W (C_1\theta + C_2\bar {q}_n + C_3\bar {\sigma }_{nn}) \,\mathrm {d}s \notag \\ + \int _{\partial \Omega } [v_{\tau _1}^W & (C_4 v_{\tau _1} + C_5 \bar {q}_{\tau _1} + C_6 \bar {\sigma }_{n\tau _1}) + v_{\tau _2}^W (C_4 v_{\tau _2} + C_5 \bar {q}_{\tau _2} + C_6 \bar {\sigma }_{n\tau _2})] \,\mathrm {d}s, \end{align}

where the constants $C_1$ — $C_6$ depend only on the coefficients in the equations and the boundary conditions. In (3.24) the right-hand side only contains terms related to $\boldsymbol {v}^W$ and $\theta ^W$ , meaning that the entropy can increase only due to the velocities and temperatures of the walls, which does not violate the second law of thermodynamics.

4.1. Choice of variables

As described in the previous section, our choice of the 13 variables is different from Grad’s 13-moment equations. To explain how the new variables are defined, we start from the following definition of trace-free moments, which has been utilised in the derivation of steady-state R13 equations (Cai et al. Reference Cai, Torrilhon and Yang2024b ):

(4.1) \begin{equation} w_{i_1 \cdots i_l}^n = \langle \psi _{i_1 \cdots i_l}^n, f \rangle . \end{equation}

Here $\psi _{i_1 \cdots i_l}^n$ is a polynomial defined by

(4.2) \begin{equation} \psi _{i_1 \cdots i_l}^n(\boldsymbol {\xi }) = \bar {L}_n^{(l+1/2)} \left ( \frac {|\boldsymbol {\xi }|^2}{2} \right ) \xi _{\langle i_1} \cdots \xi _{i_l \rangle }, \end{equation}

and $\bar {L}_n^{(l+1/2)}$ is the normalised Laguerre polynomial

(4.3) \begin{equation} \bar {L}_n^{(l+1/2)}(x) = \sqrt {\frac {\sqrt {\pi }}{2^{l+1} n! \Gamma (n+l+3/2)}} x^{-(l+1/2)} \left ( \frac {\mathrm {d}}{\mathrm {d}x} - 1 \right )^n x^{n+l+1/2}. \end{equation}

These moments can be viewed as generalisations of Grad’s 13 moments, and they are related to Grad’s moments by

(4.4) \begin{equation} \rho = w^0, \quad v_i = \sqrt {3} w_i^0, \quad \theta = -\sqrt {\frac {2}{3}} w^1, \quad \sigma _{ij} = \sqrt {15} w_{ij}^0, \quad q_i = -\sqrt {\frac {15}{2}} w_i^1. \end{equation}

By taking moments of the Boltzmann equation (1.1), one can derive evolution equations for the moments $w_{i_1 \cdots i_l}^n$ :

(4.5) \begin{align} & \frac {\partial w_{i_1\cdots i_l}^n }{\partial t} + \left ( \sqrt {2(n+l)+3} \frac {\partial w_{i_1\cdots i_l j}^n}{\partial x_j} - \sqrt {2n} \frac {\partial w_{i_1\cdots i_l j}^{n-1}}{\partial x_j} \right )\notag \\ & \quad + \frac {l}{2l+1} \left ( \sqrt {2(n+l)+1} \frac {\partial w_{\langle i_1\cdots i_{l-1}}^n}{\partial x_{i_l \rangle }} - \sqrt {2(n+1)} \frac {\partial w_{\langle i_1\cdots i_{l-1}}^{n+1}}{\partial x_{i_l\rangle }} \right )= \frac {1}{\textit{Kn}} \sum _{n'=0}^{+\infty } a_{lnn'} w_{i_1\cdots i_l}^{n'}. \end{align}

Here the coefficients $a_{lnn'}$ on the right-hand side are given by

(4.6) \begin{equation} a_{lnn'} = \frac {(2l+1)!!}{l!} \langle \psi _{i_1 \cdots i_l}^n, \mathcal {L} \psi _{i_1 \cdots i_l}^{n'} \rangle , \end{equation}

where we do not take summation over the indices $i_1, \ldots , i_l$ , and the choice of these indices does not affect the value of $a_{lnn'}$ due to the rotational invariance of the linearised collision operator. By the self-adjointness of $\mathcal {L}$ , one can observe that $a_{lnn'} = a_{ln'n}$ . Due to the conservation of mass, momentum and energy, it holds that

(4.7) \begin{equation} a_{00n} = a_{0n0} = a_{01n} = a_{0n1} = a_{10n} = a_{1n0} = 0. \end{equation}

For inverse-power-law models, the computation of these coefficients are provided in Cai & Torrilhon (Reference Cai and Torrilhon2015).

To select appropriate moments in our equations, we apply the Chapman—Enskog expansion to the moment equations, which requires assuming that $\textit{Kn}$ is a small parameter and applying the following asymptotic expansion:

(4.8) \begin{equation} w_{i_1\cdots i_l}^n = w_{i_1\cdots i_l}^{n|0} + \ \textit{Kn} w_{i_1\cdots i_l}^{n|1} + \ \textit{Kn}^2 w_{i_1\cdots i_l}^{n|2} + \ \textit{Kn}^3 w_{i_1\cdots i_l}^{n|3} + \cdots . \end{equation}

The conservative variables $w^0$ , $w_i^0$ and $w^1$ are regarded as $O(1)$ variables. Straightforward analysis leads to the following results:

(4.9) \begin{align} w^{n|0} = w^{n|1} = 0, \quad n \geqslant 2; \qquad w_i^{n|0} = 0, \quad n \geqslant 1; \end{align}

(4.10) \begin{align} w_{i_1\cdots i_l}^{n|0} = w_{i_1\cdots i_l}^{n|1} = 0, \quad l = 3; \qquad w_{i_1\cdots i_l}^{n|0} = w_{i_1\cdots i_l}^{n|1} = w_{i_1\cdots i_l}^{n|2} = 0, \quad l \geqslant 4; \end{align}

(4.11) \begin{align} w_i^{n|1} = \beta _1^{(1),n} \frac {\partial w^1}{\partial x_i}, \quad n \geqslant 1; \qquad w_{ij}^{n|1} = \beta _2^{(1),n} \frac {\partial w_{\langle i}^0}{\partial x_{j\rangle }}; \end{align}

(4.12) \begin{align} w^{n|2} = \gamma _0^{(2),n} \frac {\partial w_j^{1|1}}{\partial x_j}, \quad n \geqslant 2; \qquad w_i^{n|2} = \gamma _{t1}^{(1),n} \frac {\partial w^{1|1}_i}{\partial t} + \gamma _{s1}^{(1),n} \frac {\partial w^{0|1}_{ij}}{\partial x_j}, \quad n \geqslant 2; \end{align}

(4.13) \begin{align} w_{ij}^{n|2} = \gamma _{t2}^{(1),n} \frac {\partial w^{0|1}_{ij}}{\partial t} + \gamma _{s2}^{(1),n} \frac {\partial w^{1|1}_{\langle i}}{\partial x_{j \rangle }}, \quad n \geqslant 1; \qquad w_{ijk}^{n|2} = \gamma _3^{(2),n} \frac {\partial w_{\langle ij}^{0|1}}{\partial x_{k \rangle }}. \end{align}

The coefficients $\beta$ and $\gamma$ are all constants depending on the collision term. Note that (4.11) implies the Navier–Stokes law and the Fourier law, and these equations lead to the following linear relationship between the moments:

(4.14) \begin{equation} w_i^{n|1} = \frac {\beta _1^{(1),n}}{\beta _1^{(1),1}} w_i^{1|1}, \quad n \geqslant 1; \qquad w_{ij}^{n|1} = \frac {\beta _2^{(1),n}}{\beta _2^{(1),0}} w_{ij}^{n|0}. \end{equation}

By the definition of $w_{i_1\cdots i_l}^n$ (see (4.1)), the above equations indicate that

(4.15) \begin{equation} \left \langle \psi _i^n - \frac {\beta _1^{(1),n}}{\beta _1^{(1),1}} \psi _i^1, f \right \rangle \sim O(\textit{Kn}^2), \quad n \geqslant 1; \qquad \left \langle \psi _{ij}^n - \frac {\beta _2^{(1),n}}{\beta _2^{(1),0}} \psi _{ij}^0, f \right \rangle \sim O(\textit{Kn}^2). \end{equation}

Meanwhile, (4.9) yields

(4.16) \begin{equation} \langle \psi ^n, f \rangle \sim O(\textit{Kn}^2), \quad n \geqslant 2; \qquad \langle \psi _{i_1\cdots i_l}^n, f \rangle \sim O(\textit{Kn}^2), \quad l \geqslant 3. \end{equation}

Recall that we want to choose the 13 variables of the our equations to be $\langle \tilde {\mathbf {\phi }}_{13}, f \rangle$ with $\tilde {\mathbf {\phi }}_{13}$ satisfying $\langle r, f\rangle \sim o(\textit{Kn})$ whenever $\langle \tilde {\mathbf {\phi }}_{13}, r\rangle = 0$ . A straightforward approach is to select $\tilde {\mathbf {\phi }}_{13}$ to be a basis of $(\mathbb {S}^{(2)})^{\perp }$ , with $\mathbb {S}^{(2)}$ being the linear span of all polynomials appearing in (4.15)and (4.16):

(4.17) \begin{align} \mathbb {S}^{(2)} = \operatorname {span} \left \{ \psi _i^n - \frac {\beta _1^{(1),n}}{\beta _1^{(1),1}} \psi _i^1 \,\Bigg \vert \, n \geqslant 1\right \} \oplus \operatorname {span}\left \{ \psi _{ij}^n - \frac {\beta _2^{(1),n}}{\beta _2^{(1),0}} \psi _{ij}^0 \,\Bigg \vert \, n \geqslant 0\right \} \notag \\ \quad \oplus \operatorname {span}\{ \psi ^n \mid n \geqslant 2\} \oplus \operatorname {span}\{ \psi _{i_1\cdots i_l}^n \mid l \geqslant 3\}. \end{align}

This leads to the following choice of $\tilde {\mathbf {\phi }}_{13}$ :

(4.18) \begin{equation} \psi ^0, \quad \psi ^1, \quad \psi _i^0, \quad \phi _i^1, \quad \phi _{ij}^0. \end{equation}

Here $\psi ^0$ , $\psi ^1$ and $\psi _i^0$ correspond to conservative moments, and $\phi _i^1$ and $\phi _{ij}^0$ , corresponding to non-equilibrium variables, are defined by

(4.19) \begin{equation} \phi _i^{1} = \sum _{n=1}^{+\infty } c_1^{1,n} \psi _i^n, \qquad \phi _{ij}^{0} = \sum _{n=0}^{+\infty } c_2^{0,n} \psi _{ij}^n, \end{equation}

with the coefficients $c_1^{1,n}$ and $c_2^{0,n}$ satisfying

(4.20) \begin{equation} c_1^{1,n} = \frac {\beta _1^{(1),n}}{\beta _1^{(1),1}} c_1^{1,1}, \qquad c_2^{0,n} = \frac {\beta _2^{(1),n}}{\beta _2^{(1),0}} c_2^{0,0}. \end{equation}

Thus, the 13 moments in our equations include

(4.21) \begin{equation} \rho , \quad \theta , \quad v_i, \quad \bar {q}_i = \langle \phi _i^1, f \rangle , \quad \bar {\sigma }_{ij} = \langle \phi _{ij}^0, f \rangle . \end{equation}

In (4.19) the constants $c_1^{1,1}$ and $c_2^{0,0}$ can be chosen as arbitrary non-zero numbers, and here we choose them such that when the collision model reduces Maxwell molecules ( $\eta = 5$ ), $\bar {q}_i$ and $\bar {\sigma }_{ij}$ reduce to $q_i$ and $\sigma _{ij}$ . More precisely, we choose negative $c_1^{1,1}$ and positive $c_2^{0,0}$ such that

(4.22) \begin{equation} \sum _{n=1}^{+\infty } \left ( c_1^{1,n} \right )^2 = \frac {15}{2}, \qquad \sum _{n=0}^{+\infty } \left ( c_2^{0,n} \right )^2 = 15. \end{equation}

Note that the definitions of $\bar {q}_i$ and $\bar {\sigma }_{ij}$ depend on the collision model since the coefficients $c_1^{1,n}$ and $c_2^{0,n}$ are determined by $\beta _1^{(1),n}$ and $\beta _2^{(1),0}$ , and these constants will vary when the collision model changes.

4.2. Second-order variables

To derive a model with the super-Burnett order, second-order moments are needed as auxiliary variables. The choice of these variables also follows the requirement of orthogonality, meaning that we want to find a vector function $\tilde {\psi }(\boldsymbol {\xi })$ such that:

1. (i) $\langle \tilde {\phi }_{13}, \tilde {\psi }\rangle$ = 0;

2. (ii) for any function $r(\boldsymbol {\xi })$ such that $\langle r,\tilde {\psi }\rangle = 0$ , it holds that $\langle r, f \rangle \sim o(\textit{Kn}^2)$ ;

3. (iii) the number of components of $\tilde {\psi }$ should be as low as possible.

To find $\tilde {\psi }$ , we need the following linear relationship between the second-order terms of moments:

(4.23) \begin{align} w^{n|2} = d^n_{02} w^{2|2}, \quad n \geqslant 2; \qquad w^{n|2}_i = d^n_{12} w^{2|2}_i + d^n_{13} w^{3|2}_i, \quad n \geqslant 2; \notag\\ w^{n|2}_{ij} = d^n_{21} w^{1|2}_{ij} + d^n_{22} w^{2|2}_{ij}, \quad n \geqslant 1; \qquad w^{n|2}_{ijk} = d^n_{30} w^{0|2}_{ijk}, \quad n \geqslant 0, \end{align}

where

(4.24) \begin{align} d^n_{02} &= \frac {\gamma _0^{(2),n}}{\gamma _0^{(2),2}}, \qquad\qquad\qquad\qquad\qquad d^n_{30} = \frac {\gamma _3^{(2),n}}{\gamma _3^{(2),0}}, \nonumber \\ d^n_{12} &= \frac {\gamma ^{(1),n}_{t1} \gamma ^{(1),3}_{s1}-\gamma ^{(1),n}_{s1} \gamma ^{(1),3}_{t1}}{\gamma ^{(1),3}_{s1} \gamma ^{(1),2}_{t1} - \gamma ^{(1),2}_{s1} \gamma ^{(1),3}_{t1}}, \qquad d^n_{13} = \frac {\gamma ^{(1),n}_{s1} \gamma ^{(1),2}_{t1}-\gamma ^{(1),n}_{t1} \gamma ^{(1),2}_{s1}}{\gamma ^{(1),3}_{s1} \gamma ^{(1),2}_{t1} - \gamma ^{(1),2}_{s1} \gamma ^{(1),3}_{t1}}, \nonumber \\ d^n_{21} &= \frac {\gamma ^{(1),n}_{t2} \gamma ^{(1),2}_{s2}-\gamma ^{(1),n}_{s2} \gamma ^{(1),2}_{t2}}{\gamma ^{(1),2}_{s2} \gamma ^{(1),1}_{t2} - \gamma ^{(1),1}_{s2} \gamma ^{(1),2}_{t2}},\qquad d^n_{22} = \frac {\gamma ^{(1),n}_{s2} \gamma ^{(1),1}_{t2}-\gamma ^{(1),n}_{t2} \gamma ^{(1),1}_{s2}}{\gamma ^{(1),2}_{s2} \gamma ^{(1),1}_{t2} - \gamma ^{(1),1}_{s2} \gamma ^{(1),2}_{t2}}. \end{align}

These equations can be derived from (4.12), (4.13) by cancelling the differential terms. The equations in (4.23) reveal that low-order parts of the moments can be cancelled by linear combinations:

(4.25) \begin{align} w^n - d_{02}^n w^2 \sim o\left(\textit{Kn}^2\right), \quad n \geqslant 2; \notag\\ \left(w_i^n - \ \textit{Kn} w_i^{n|1}\right) - d_{12}^n \left(w_i^2 - {\textit{Kn}} w_i^{2|1}\right) - d_{13}^n \left(w_i^3 - \ \textit{Kn} w_i^{3|1}\right) \sim o\left({\textit{Kn}}^2\right), \quad n \geqslant 2; \notag \\ \left(w_{ij}^n - {\textit{Kn}} w_{ij}^{n|1}\right) - d_{21}^n \left(w_{ij}^1 - {\textit{Kn}} w_{ij}^{1|1}\right) - d_{22}^n \left(w_{ij}^2 - {\textit{Kn}} w_{ij}^{2|1}\right) \sim o\left({\textit{Kn}}^2\right), \quad n \geqslant 1; \notag \\ w_{ijk}^n - d_{30} w_{ijk}^0 \sim o\left({\textit{Kn}}^2\right), \quad n \geqslant 0. \end{align}

Note that the correctness of (4.12) and (4.13) (and thus, (4.23) and (4.25)) replies only on the relation (4.14), and does not require the specific choice of $w_i^{n|1}$ and $w_{ij}^{n|1}$ as described in (4.11). Therefore, in (4.25) we can follow (4.15) to select

(4.26) \begin{equation} \textit{Kn} w_i^{n|1} = \left \langle \frac {\beta _1^{(1),n}}{\beta _1^{(1),1}} \psi _i^1, f \right \rangle , \quad n \geqslant 1; \qquad \textit{Kn} w_{ij}^{n|1} = \left \langle \frac {\beta _2^{(1),n}}{\beta _2^{(1),0}} \psi _{ij}^0, f \right \rangle , \quad n \geqslant 0, \end{equation}

so that we can find a linear space $\mathbb {S}^{(3)}$ such that $\langle r, f \rangle \sim o(\textit{Kn}^2)$ for all $r \in \mathbb {S}^{(3)}$ . The definition of $\mathbb {S}^{(3)}$ is

(4.27) \begin{align} \mathbb {S}^{(3)} &= \operatorname {span} \left \{ \psi ^n - d_{02}^n \psi ^2 \mid n \geqslant 3\right \} \notag\\ &\oplus \operatorname {span}\left \{ \left (\!\!\psi _i^n - \frac {\beta _1^{(1),n}}{\beta _1^{(1),1}} \psi _i^1 \!\right ) - d_{12}^n\left (\!\! \psi _i^2 - \frac {\beta _1^{(1),2}}{\beta _1^{(1),1}} \psi _i^1 \!\right ) - d_{13}^n \left (\!\! \psi _i^3 - \frac {\beta _1^{(1),3}}{\beta _1^{(1),1}} \psi _i^1\!\right )\Bigg \vert \, n \geqslant 4\right \} \notag\\ &\oplus \operatorname {span}\left \{ \left (\!\!\psi _{ij}^n - \frac {\beta _2^{(1),n}}{\beta _2^{(1),0}} \psi _{ij}^0 \!\right ) - d_{21}^n\left (\!\! \psi _{ij}^1 - \frac {\beta _2^{(1),1}}{\beta _2^{(1),0}} \psi _{ij}^0 \!\right ) - d_{22}^n \left (\! \psi _{ij}^2 - \frac {\beta _2^{(1),2}}{\beta _2^{(1),0}} \psi _{ij}^0\!\right ) \Bigg \vert \, n \geqslant 3\right \} \notag\\ & \oplus \operatorname {span}\left \{ \psi _{ijk}^n - d_{30} \psi _{ijk}^0 \mid n \geqslant 1\right \} \oplus \operatorname {span}\left \{ \psi _{i_1 \cdots i_l}^n \mid l \geqslant 4\right \}. \end{align}

The second-order variables are chosen to be a basis of $\mathbb {S}^{(2)} \cap (\mathbb {S}^{(3)})^{\perp }$ .

The function space $\mathbb {S}^{(2)} \cap (\mathbb {S}^{(3)})^{\perp }$ has the form

(4.28) \begin{equation} \mathbb {S}^{(2)} \cap (\mathbb {S}^{(3)})^{\perp } = \operatorname {span} \{ \phi ^{2}, \phi _{i}^{2}, \phi _{i}^{3}, \phi _{ij}^{1}, \phi _{ij}^{2},\phi _{ijk}^{0} \}, \end{equation}

where

(4.29) \begin{align} \phi ^{2} &= \sum _{n=2}^{+\infty } c_0^{2,n} \psi ^n, \quad \phi _{i}^{2} = \sum _{n=2}^{+\infty } c_{1}^{2,n} \left (\!\psi _{i}^{n} - \frac {\beta _1^{(1),n}}{\beta _1^{(1),1}} \psi ^1_i\!\right ), \!\!\quad \phi _{i}^{3} = \sum _{n=2}^{+\infty } c_{1}^{3,n} \left (\!\psi _{i}^{n} - \frac {\beta _1^{(1),n}}{\beta _1^{(1),1}} \psi ^1_i\!\right ), \notag\\ \phi _{ij}^{1} &= \sum _{n=1}^{+\infty }\! c_{2}^{1,n} \left (\!\psi _{ij}^{n} - \frac {\beta _2^{(1),n}}{\beta _2^{(1),0}} \psi ^0_{ij}\!\right ), \!\!\quad\!\! \phi _{ij}^{2} = \sum _{n=1}^{+\infty }\! c_{2}^{2,n} \left (\!\psi _{ij}^{n} - \frac {\beta _2^{(1),n}}{\beta _2^{(1),0}} \psi ^0_{ij}\!\right ), \!\!\quad\!\! \phi _{ijk}^{0} = \sum _{n=0}^{+\infty } c_3^{0,n} \psi _{ijk}^n. \end{align}

The coefficients $c_k^{\ell ,n}$ can be represented as functions of $\beta _k^{(\ell ), n}$ , and the details can be found in Appendix D. Here we only emphasise that the choice has been made such that $\langle \phi _i^2, \phi _i^3\rangle = \langle \phi _{ij}^1, \phi _{ij}^2 \rangle = 0$ , and in the case of Maxwell molecules, $\phi _i^{2,3} = \psi _i^{2,3}$ and $\phi _{ij}^{1,2} = \psi _{ij}^{1,2}$ . Equation (4.28) indicates that the minimum number of second-order variables is 24 ( $= 1+3+3+5+5+7$ ). These variables will be denoted by

(4.30) \begin{align} u^2 &:= \langle \phi ^2, f \rangle , \qquad u_i^2 := \langle \phi _i^2, f \rangle , \qquad u_i^3 := \langle \phi _i^3, f \rangle , \notag\\ u_{ij}^1 &:= \langle \phi _{ij}^1, f \rangle , \qquad u_{ij}^2 := \langle \phi _{ij}^2, f \rangle , \qquad u_{ijk}^0 := \langle \phi _{ijk}^0, f \rangle . \end{align}

4.3. Derivation of R13 equations

By the analysis above, the distribution function $f$ can be separated into four parts according to the order of magnitude

(4.31) \begin{equation} f = f^{(0)} + f^{(1)} + f^{(2)} + f^{(\mathrm {r})}, \end{equation}

where $f^{(k)}$ is the projection of $f$ onto the function space $\mathbb {V}^{(k)}$ , defined by

(4.32) \begin{align} \mathbb {V}^{(0)} &= \operatorname {span}\{\psi ^0, \psi ^1, \psi _i^0\}, \qquad \mathbb {V}^{(1)} = (\mathbb {V}^{(0)})^{\perp } \cap (\mathbb {S}^{(2)})^{\perp }, \nonumber \\ \mathbb {V}^{(2)} &= \mathbb {S}^{(2)} \cap (\mathbb {S}^{(3)})^{\perp }, \qquad \mathbb {V}^{(\mathrm {r})} = \mathbb {S}^{(3)}. \end{align}

It is clear that $f^{(0)}$ , $f^{(1)}$ and $f^{(2)}$ represent, respectively, the zeroth-, first- and second-order parts of $f$ , and $f^{(\mathrm {r})}$ represents the remaining part of $f$ that has order $o(\textit{Kn}^2)$ . For simplicity, we use $\mathcal {P}^{(k)}$ to denote the projection operator onto $\mathbb {V}^{(k)}$ . A key observation in our derivation is that the super-Burnett equations can be derived from the following equations without using $f^{(\mathrm {r})}$ :

(4.33) \begin{align} \frac {\partial }{\partial t} \begin{pmatrix} f^{(0)} \\ f^{(1)} \\ 0 \end{pmatrix} &+ \begin{pmatrix} \partial _{x_i} \mathcal {P}^{(0)} \xi _i f^{(0)} + \partial _{x_i} \mathcal {P}^{(0)} \xi _i f^{(1)} + \partial _{x_i} \mathcal {P}^{(0)} \xi _i f^{(2)} \\ \partial _{x_i} \mathcal {P}^{(1)} \xi _i f^{(0)} + \partial _{x_i} \mathcal {P}^{(1)} \xi _i f^{(1)} + \partial _{x_i} \mathcal {P}^{(1)} \xi _i f^{(2)} \\ \partial _{x_i} \mathcal {P}^{(2)} \xi _i f^{(0)} + \partial _{x_i} \mathcal {P}^{(2)} \xi _i f^{(1)} \end{pmatrix} \notag \\ &= \frac {1}{\textit{Kn}} \begin{pmatrix} 0 \\ \mathcal {P}^{(1)} \mathcal {L} f^{(1)} + \mathcal {P}^{(1)} \mathcal {L} f^{(2)} \\ \mathcal {P}^{(2)} \mathcal {L} f^{(1)} + \mathcal {P}^{(2)} \mathcal {L} f^{(2)} \end{pmatrix}. \end{align}

This form is similar to Grad’s moment equations with truncation up to the second order, but the temporal and spatial derivatives of $f^{(2)}$ in the last equation are removed, allowing us to express $f^{(2)}$ using $f^{(0)}$ and $f^{(1)}$ according to the last equation:

(4.34) \begin{equation} f^{(2)} = \left [\mathcal {P}^{(2)}\mathcal {L}|_{\mathbb {V}^{(2)}} \right ]^{-1} \left ( \partial _{x_i} \mathcal {P}^{(2)} \xi _i f^{(0)} + \partial _{x_i} \mathcal {P}^{(2)} \xi _i f^{(1)} - \mathcal {P}^{(2)} \mathcal {L} f^{(1)} \right ). \end{equation}

Thus, the second-order part $f^{(2)}$ can be eliminated from the equations, and the resulting equations contain only $f^{(0)}$ and $f^{(1)}$ , which can be represented as 13-moment equations.

The reason why (4.33) has the super-Burnett order will be detailed in the next subsection. Here we apply this result to derive the R13 equations. The derivation requires explicit formulations of $f^{(0)}$ , $f^{(1)}$ and $f^{(2)}$ :

(4.35) \begin{equation} \begin{aligned} f^{(0)} &= \rho \psi ^{0} - \sqrt {\frac {3}{2}} \theta \psi ^{1} + \sqrt {3} v_i \psi _i^{0}, \\ f^{(1)} &= \frac {2}{5} \bar {q}_i \phi _i^{1} + \frac {1}{2} \bar {\sigma }_{ij} \phi _{ij}^{0}, \\ f^{(2)} &= u^{2} \phi ^{2} + 3u_i^{2} \phi _i^{2} + 3u_i^{3} \phi _i^{3}+ \frac {15}{2}u_{ij}^{1} \phi _{ij}^{1} + \frac {15}{2}u_{ij}^{2} \phi _{ij}^{2} + \frac {35}{2}u_{ijk}^{0} \phi _{ijk}^{0}. \end{aligned} \end{equation}

Then, (4.33) can be converted to moment equations by the Galerkin method. For example, the equation of $\bar {\sigma }_{ij}$ can be obtained by

(4.36) \begin{equation} \langle \phi _{ij}^0, f^{(1)} \rangle + \frac {\partial }{\partial x_k} \langle \phi _{ij}^0, \xi _k (f^{(0)} + f^{(1)} + f^{(2)}) \rangle = \langle \phi _{ij}^0, \mathcal {L} (f^{(1)} + f^{(2)}) \rangle , \end{equation}

and the equation of $u_i^2$ is

(4.37) \begin{equation} \frac {\partial }{\partial x_k} \langle \phi _i^2, \xi _k(f^{(0)} + f^{(1)}) \rangle = \langle \phi _i^2, \mathcal {L}(f^{(1)} + f^{(2)}) \rangle . \end{equation}

All these inner products can be computed explicitly, and the complete 37-moment equations are

(4.38) \begin{align} &\frac {\partial \rho }{\partial t}+ \frac {\partial v_j}{\partial x_j} = 0, \end{align}

(4.39) \begin{align} &\frac {\partial \theta }{\partial t} + \frac {2}{3} \frac {\partial v_j}{\partial x_j} + \frac {2}{3} c_1^{1,1}\frac {\partial \bar {q}_j}{\partial x_j} - \sqrt {\frac {10}{3}}\left ( c^{2,1}_1 \frac {\partial u^{2}_j}{\partial x_j} + c^{3,1}_1 \frac {\partial u^{3}_j}{\partial x_j} \right ) = 0 , \end{align}

(4.40) \begin{align} &\frac {\partial v_i}{\partial t} + \frac {\partial \rho }{\partial x_i} + \frac {\partial \theta }{\partial x_i} + c_2^{0,0} \frac {\partial \bar {\sigma }_{ij}}{\partial x_j} + \sqrt {15} \left ( c^{1,0}_2 \frac {\partial u^{1}_{ij}}{\partial x_j} + c^{2,0}_2 \frac {\partial u^{2}_{ij}}{\partial x_j} \right ) = 0, \end{align}

(4.41) \begin{align} & \frac {\partial \bar {q}_i}{\partial t} + \frac {5}{2}c_1^{1,1} \frac {\partial \theta }{\partial x_i} - \frac {\sqrt {2}}{6} A_{45} \frac {\partial \bar {\sigma }_{ij}}{\partial x_j} - \sqrt {\frac {5}{6}} \left ( A_{49}\frac {\partial u^{1}_{ij}}{\partial x_j} + A_{4,10}\frac {\partial u^{2}_{ij}}{\partial x_j} + A_{46} \frac {\partial u^{2}}{\partial x_i} \right ) \notag \\ & \hspace {120pt}= \frac {1}{\textit{Kn}} \left( \frac {1}{3} \mathscr {L}^{\,(11)}_1 \bar {q}_{i} - \sqrt {\frac {5}{6}} \mathscr {L}^{\,(12)}_1 u^{2}_{i} - \sqrt {\frac {5}{6}} \mathscr {L}^{\,(13)}_1 u^{3}_{i} \right), \end{align}

(4.42) \begin{align} & \frac {\partial \bar {\sigma }_{ij}}{\partial t} + 2 c_2^{0,0} \frac {\partial v_{\langle i}}{\partial x_{j \rangle }} -\frac {2\sqrt {2}}{15} A_{45}\frac {\partial \bar {q}_{\langle i}}{\partial x_{j \rangle }} + \frac {2}{\sqrt {15}} \left ( A_{57}\frac {\partial u^{2}_{\langle i}}{\partial x_{j \rangle }} + A_{58}\frac {\partial u^{3}_{\langle i}}{\partial x_{j \rangle }} + A_{5,11} \frac {\partial u^{0}_{ijk}}{\partial x_k} \right ) \notag \\ & \hspace {100pt}= \frac {1}{\textit{Kn}}\left(\frac {2}{15} \mathscr {L}^{\,(00)}_2 \bar {\sigma }_{ij}+ \frac {2}{\sqrt {15}} \mathscr {L}^{\,(01)}_2 u^{1}_{ij}+ \frac {2}{\sqrt {15}} \mathscr {L}^{\,(02)}_2 u^{2}_{ij}\right), \end{align}

(4.43) \begin{align} & -\sqrt {\frac {2}{15}} A_{46} \frac {\partial \bar {q}_j}{\partial x_j} = \frac {1}{\textit{Kn}}\mathscr {L}^{\,(22)}_0 u^{2}, \end{align}

(4.44) \begin{align} & -\sqrt {\frac {15}{2}} c^{2,1}_1 \frac {\partial \theta }{\partial x_i} + \sqrt {\frac {1}{15}} A_{57} \frac {\partial \bar {\sigma }_{ij}}{\partial x_j} = \frac {1}{\textit{Kn}}\left(-\sqrt {\frac {2}{15}} \mathscr {L}^{\,(21)}_1 \bar {q}_{i}+\mathscr {L}^{\,(22)}_1 u^{2}_{i}+\mathscr {L}^{\,(23)}_1 u^{3}_{i}\right), \end{align}

(4.45) \begin{align} & -\sqrt {\frac {15}{2}} c^{3,1}_1 \frac {\partial \theta }{\partial x_i} + \sqrt {\frac {1}{15}} A_{58} \frac {\partial \bar {\sigma }_{ij}}{\partial x_j} = \frac {1}{\textit{Kn}}\left(-\sqrt {\frac {2}{15}}\mathscr {L}^{\,(31)}_1 \bar {q}_{i}+\mathscr {L}^{\,(32)}_1 u^{2}_{i}+\mathscr {L}^{\,(33)}_1 u^{3}_{i}\right), \end{align}

(4.46) \begin{align} & \sqrt {15} c^{1,0}_2 \frac {\partial v_{\langle i}}{\partial x_{j \rangle }} -\sqrt {\frac {2}{15}} A_{49}\frac {\partial \bar {q}_{\langle i}}{\partial x_{j \rangle }} = \frac {1}{\textit{Kn}}\left(\sqrt {\frac {1}{15}} \mathscr {L}^{\,(10)}_2 \bar {\sigma }_{ij}+\mathscr {L}^{\,(11)}_2 u^{1}_{ij}+\mathscr {L}^{\,(12)}_2 u^{2}_{ij}\right), \end{align}

(4.47) \begin{align} & \sqrt {15} c^{2,0}_2 \frac {\partial v_{\langle i}}{\partial x_{j \rangle }} -\sqrt {\frac {2}{15}} A_{4,10}\frac {\partial \bar {q}_{\langle i}}{\partial x_{j \rangle }} = \frac {1}{\textit{Kn}}\left(\sqrt {\frac {1}{15}} \mathscr {L}^{\,(20)}_2 \bar {\sigma }_{ij}+\mathscr {L}^{\,(21)}_2 u^{1}_{ij}+\mathscr {L}^{\,(22)}_2 u^{2}_{ij}\right), \end{align}

(4.48) \begin{align} & \sqrt {\frac {1}{15}} A_{5,11} \frac {\partial \bar {\sigma }_{\langle ij}}{\partial x_{k \rangle }} = \frac {1}{\textit{Kn}} \mathscr {L}^{\,(00)}_3 u^{0}_{ijk}. \end{align}

The coefficients $A_{ij}$ are a series of $c_k^{l,n}$ , with their expressions provided in Appendix E.

It is now clear that the second-order variables defined in (4.30) can be solved according to (4.43)–(4.48). Plugging the results into (4.39)–(4.42) yields our final equations (3.2)–(3.6). Here we would like to write down the explicit expressions for a few coefficients in (3.2)–(3.6):

(4.49) \begin{align} k_1 &= -5 \left(\begin{array}{l@{\quad}l} c_1^{2,1} & c_1^{3,1} \end{array} \right) \left(\begin{array}{l@{\quad}l} \mathscr {L}_1^{\,(22)} & \mathscr {L}_1^{\,(23)} \\[3pt] \mathscr {L}_1^{\,(32)} & \mathscr {L}_1^{\,(33)} \end{array} \right)^{\!-1}\begin{pmatrix} c_1^{2,1} \\[3pt] c_1^{3,1} \end{pmatrix}, \end{align}

(4.50) \begin{align} k_2 &= -\frac {\sqrt {2}}{3} \left(\begin{array}{l@{\quad}l} c_1^{2,1} & c_1^{3,1} \end{array}\right) \left(\begin{array}{l@{\quad}l}\mathscr {L}_1^{\,(22)} & \mathscr {L}_1^{\,(23)} \\[3pt] \mathscr {L}_1^{\,(32)} & \mathscr {L}_1^{\,(33)} \end{array}\right)^{\!-1}\begin{pmatrix} A_{57} \\ A_{58} \end{pmatrix}, \end{align}

(4.51) \begin{align} k_3 &= -15 \left(\begin{array}{l@{\quad}l} c_2^{1,0} & c_2^{2,0} \end{array}\right) \left(\begin{array}{l@{\quad}l} \mathscr {L}_2^{\,(11)} & \mathscr {L}_2^{\,(12)} \\[3pt] \mathscr {L}_2^{\,(21)} & \mathscr {L}_2^{\,(22)} \end{array} \right)^{\!-1}\begin{pmatrix} c_2^{1,0} \\[3pt] c_2^{2,0} \end{pmatrix}, \end{align}

(4.52) \begin{align} k_4 &= \sqrt {2} \left(\begin{array}{l@{\quad}l} A_{49} & A_{4,10} \end{array}\right) \left(\begin{array}{l@{\quad}l} \mathscr {L}_2^{\,(11)} & \mathscr {L}_2^{\,(12)} \\[3pt] \mathscr {L}_2^{\,(21)} & \mathscr {L}_2^{\,(22)} \end{array} \right)^{\!-1}\begin{pmatrix} c_2^{1,0} \\[3pt] c_2^{2,0} \end{pmatrix}, \end{align}

(4.53) \begin{align} k_7 &= -\frac {5}{36} \left(\begin{array}{l@{\quad}l} A_{49} & A_{4,10} \end{array}\right) \left(\begin{array}{l@{\quad}l} \mathscr {L}_2^{\,(11)} & \mathscr {L}_2^{\,(12)} \\[3pt] \mathscr {L}_2^{\,(21)} & \mathscr {L}_2^{\,(22)} \end{array}\right)^{\!-1}\begin{pmatrix} A_{49} \\ A_{4,10} \end{pmatrix}, \end{align}

(4.54) \begin{align} &= -\frac {2}{15} \left(\begin{array}{l@{\quad}l} A_{57} & A_{58} \end{array}\right) \left(\begin{array} {l@{\quad}l} \mathscr {L}_1^{\,(22)} & \mathscr {L}_1^{\,(23)} \\[3pt] \mathscr {L}_1^{\,(32)} & \mathscr {L}_1^{\,(33)} \end{array}\right)^{\!-1}\begin{pmatrix} A_{57} \\ A_{58} \end{pmatrix}. \end{align}

Since the matrices

(4.55) \begin{equation} \left(\begin{array}{l@{\quad}l} \mathscr {L}_1^{\,(22)} & \mathscr {L}_1^{\,(23)} \\[3pt] \mathscr {L}_1^{\,(32)} & \mathscr {L}_1^{\,(33)} \end{array} \right) \quad \text {and} \quad \left(\begin{array}{l@{\quad}l} \mathscr {L}_2^{\,(11)} & \mathscr {L}_2^{\,(12)} \\[3pt] \mathscr {L}_2^{\,(21)} & \mathscr {L}_2^{\,(22)} \end{array} \right) \end{equation}

are symmetric and negative definite, we can obtain (3.20) by the Cauchy–Schwarz inequality, which completes the proof of the H-theorem (3.22).

4.4. Verification of the super-Burnett order

To show that (4.33) have super-Burnett order, we rewrite the Boltzmann equation in the form

(4.56) \begin{align} & \frac {\partial }{\partial t} \begin{pmatrix} f^{(0)} \\ f^{(1)} \\ f^{(2)} \\ f^{(\mathrm {r})} \end{pmatrix} + \left(\begin{array}{l@{\quad}l@{\quad}l@{\quad}l} \mathcal {A}^{(00)} & \mathcal {A}^{(01)} & \mathcal {A}^{(02)} &\mathcal {A}^{(0\mathrm {r})} \\ \mathcal {A}^{(10)} & \mathcal {A}^{(11)} & \mathcal {A}^{(12)} & \mathcal {A}^{(1\mathrm {r})} \\ \mathcal {A}^{(20)} & \mathcal {A}^{(21)} & \mathcal {A}^{(22)} & \mathcal {A}^{(2\mathrm {r})} \\ \mathcal {A}^{(\mathrm {r}0)} & \mathcal {A}^{(\mathrm {r}1)} & \mathcal {A}^{(\mathrm {r}2)} & \mathcal {A}^{(\mathrm {rr})} \end{array} \right) \begin{pmatrix} f^{(0)} \\ f^{(1)} \\ f^{(2)} \\ f^{(\mathrm {r})} \end{pmatrix} \notag\\ & = \frac {1}{\textit{Kn}} \left(\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} 0 &\quad 0 &\quad 0 &\quad 0 \\ 0 & \mathcal {L}^{(11)} & \mathcal {L}^{(12)} & \mathcal {L}^{(1\mathrm {r})} \\ 0 & \mathcal {L}^{(21)} & \mathcal {L}^{(22)} & \mathcal {L}^{(2\mathrm {r})} \\ 0 & \mathcal {L}^{(\mathrm {r}1)} & \mathcal {L}^{(\mathrm {r}2)} & \mathcal {L}^{(\mathrm {rr})} \end{array}\right) \begin{pmatrix} f^{(0)} \\ f^{(1)} \\ f^{(2)} \\ f^{(\mathrm {r})} \end{pmatrix}, \end{align}

where

(4.57) \begin{equation} \mathcal {A}^{(kl)} f^{(l)} = \frac {\partial }{\partial x_i} \mathcal {A}_i^{(kl)} f^{(l)}, \qquad \mathcal {A}_i^{(kl)} = \mathcal {P}^{(k)} \xi _i \mathcal {P}^{(l)}, \qquad \mathcal {L}^{(kl)} f^{(l)} = \mathcal {P}^{(k)} \mathcal {L} f^{(l)}, \end{equation}

and they satisfy

(4.58) \begin{equation} \langle g^{(k)}, \mathcal {A}_i^{(kl)} g^{(l)} \rangle = \langle \mathcal {A}_i^{(lk)} g^{(k)}, g^{(l)} \rangle , \qquad \langle g^{(k)}, \mathcal {L}^{(kl)} g^{(l)} \rangle = \langle \mathcal {L}^{(lk)} g^{(k)}, g^{(l)} \rangle \end{equation}

for all $g^{(k)} \in \mathbb {V}^{(k)}$ and $g^{(l)} \in \mathbb {V}^{(l)}$ .

It is well known that the super-Burnett equations can be derived via three Maxwellian iterations applied to the Boltzmann equation, which begin with the initial value

(4.59) \begin{equation} f_0^{(1)} = f_0^{(2)} = f_0^{(r)} = 0, \end{equation}

and proceeds from the $j$ th iteration to the $(j+1)$ th iteration according to

(4.60) \begin{equation} \frac {\partial }{\partial t} \begin{pmatrix} f_j^{(1)} \\ f_j^{(2)} \\ f_j^{(\mathrm {r})} \end{pmatrix} + \begin{pmatrix} \mathcal {A}^{(10)} & \mathcal {A}^{(11)} & \mathcal {A}^{(12)} & \mathcal {A}^{(1\mathrm {r})} \\ \mathcal {A}^{(20)} & \mathcal {A}^{(21)} & \mathcal {A}^{(22)} & \mathcal {A}^{(2\mathrm {r})} \\ \mathcal {A}^{(\mathrm {r}0)} & \mathcal {A}^{(\mathrm {r}1)} & \mathcal {A}^{(\mathrm {r}2)} & \mathcal {A}^{(\mathrm {rr})} \end{pmatrix} \begin{pmatrix} f^{(0)} \\ f_j^{(1)} \\ f_j^{(2)} \\ f_j^{(\mathrm {r})} \end{pmatrix} = \frac {1}{\textit{Kn}} \begin{pmatrix} \mathcal {L}^{(11)} & \mathcal {L}^{(12)} & \mathcal {L}^{(1\mathrm {r})} \\ \mathcal {L}^{(21)} & \mathcal {L}^{(22)} & \mathcal {L}^{(2\mathrm {r})} \\ \mathcal {L}^{(\mathrm {r}1)} & \mathcal {L}^{(\mathrm {r}2)} & \mathcal {L}^{(\mathrm {rr})} \end{pmatrix} \begin{pmatrix} f_{j+1}^{(1)} \\[2pt] f_{j+1}^{(2)} \\[2pt] f_{j+1}^{(\mathrm {r})} \end{pmatrix}. \end{equation}

The super-Burnett equations can be written in the form

(4.61) \begin{equation} \frac {\partial f^{(0)}}{\partial t} + \mathcal {A}^{(00)} f^{(0)} + \mathcal {A}^{(01)} f_3^{(1)} + \mathcal {A}^{(02)} f_3^{(2)} + \mathcal {A}^{(02)} f_3^{(3)} = 0. \end{equation}

To show that (4.33) have the super-Burnett order, we just need to demonstrate that the derivation of super-Burnett equations does not involve any operators that do not appear in (4.33). Below, we carry out the derivation step by step in the following subsections.

4.4.1. First Maxwellian iteration

Setting $j = 0$ in (4.60) and using the initial data $f_0^{(1)} = f_0^{(2)} = f_0^{(\mathrm {r})} = 0$ , we obtain

(4.62) \begin{equation} \left(\begin{array}{c} \mathcal {A}^{(10)} f^{(0)} \\ \mathcal {A}^{(20)} f^{(0)} \\ \mathcal {A}^{(\mathrm {r}0)} f^{(0)} \end{array} \right) = \frac {1}{\textit{Kn}} \left(\begin{array}{c@{\quad}c@{\quad}c} \mathcal {L}^{(11)} & \mathcal {L}^{(12)} & \mathcal {L}^{(1\mathrm {r})} \\ \mathcal {L}^{(21)} & \mathcal {L}^{(22)} & \mathcal {L}^{(2\mathrm {r})} \\ \mathcal {L}^{(\mathrm {r}1)} & \mathcal {L}^{(\mathrm {r}2)} & \mathcal {L}^{(\mathrm {rr})} \end{array} \right) \begin{pmatrix} f_1^{(1)} \\[2pt] f_1^{(2)} \\[2pt] f_1^{(\mathrm {r})} \end{pmatrix}. \end{equation}

By Gaussian elimination,

(4.63) \begin{equation} \begin{pmatrix} \mathcal {A}^{(10)} f^{(0)} \\ \left [\mathcal {A}^{(20)} - \mathcal {L}^{(21)} (\mathcal {L}^{(11)})^{-1} \mathcal {A}^{(10)}\right ] f^{(0)} \\ \left [\mathcal {A}^{(\mathrm {r}0)} - \mathcal {L}^{(\mathrm {r}1)} (\mathcal {L}^{(11)})^{-1} \mathcal {A}^{(10)}\right ] f^{(0)} \\ \end{pmatrix} = \frac {1}{\textit{Kn}} \left(\begin{array}{l@{\quad}l@{\quad}l} \mathcal {L}^{(11)} & \mathcal {L}^{(12)} & \mathcal {L}^{(1\mathrm {r})} \\ & \mathcal {L}_*^{(22)} & \mathcal {L}_*^{(2\mathrm {r})} \\ & \mathcal {L}_*^{(\mathrm {r}2)} & \mathcal {L}_*^{(\mathrm {rr})} \end{array}\right) \begin{pmatrix} f_1^{(1)} \\[2pt] f_1^{(2)} \\[2pt] f_1^{(\mathrm {r})} \end{pmatrix}, \end{equation}

where

(4.64) \begin{align} \mathcal {L}_*^{(22)} &= \mathcal {L}^{(22)} - \mathcal {L}^{(21)} (\mathcal {L}^{(11)})^{-1} \mathcal {L}^{(12)}, \qquad \mathcal {L}_*^{(2\mathrm {r})} = \mathcal {L}^{(2\mathrm {r})} - \mathcal {L}^{(21)} (\mathcal {L}^{(11)})^{-1} \mathcal {L}^{(1\mathrm {r})}, \nonumber \\ \mathcal {L}_*^{(\mathrm {r}2)} &= \mathcal {L}^{(\mathrm {r}2)} - \mathcal {L}^{(\mathrm {r}1)} (\mathcal {L}^{(11)})^{-1} \mathcal {L}^{(12)}, \qquad \mathcal {L}_*^{(\mathrm {rr})} = \mathcal {L}^{(\mathrm {rr})} - \mathcal {L}^{(\mathrm {r}1)} (\mathcal {L}^{(11)})^{-1} \mathcal {L}^{(1\mathrm {r})}. \end{align}

Note that throughout the iterations, we always have $f_j^{(2)} \sim o(\textit{Kn})$ and $f_j^{(r)} \sim o(\textit{Kn}^2)$ , which requires

(4.65) \begin{equation} \mathcal {A}^{(20)} - \mathcal {L}^{(21)} (\mathcal {L}^{(11)})^{-1} \mathcal {A}^{(10)} = \mathcal {A}^{(\mathrm {r}0)} - \mathcal {L}^{(\mathrm {r}1)} (\mathcal {L}^{(11)})^{-1} \mathcal {A}^{(10)} = 0. \end{equation}

Thus, the result of the first Maxwellian iteration can be written as

(4.66) \begin{equation} f_1^{(1)} = {\textit{Kn}} (\mathcal {L}^{(11)})^{-1} \mathcal {A}^{(10)} f^{(0)}, \qquad f_1^{(2)} = 0, \qquad f_1^{(\mathrm {r})} = 0. \end{equation}

For simplicity, we define $\mathcal {S}_1^{(10)} = (\mathcal {L}^{(11)})^{-1} \mathcal {A}^{(10)}$ , so that $f_1^{(1)} = {\textit{Kn}} \mathcal {S}_1^{(10)} f^{(0)}$ .

4.4.2. Second Maxwellian iteration

Setting $j = 0$ in (4.60) and using the result of the first Maxwellian iteration, we get

(4.67) \begin{align} \frac {\partial }{\partial t} \begin{pmatrix} {\textit{Kn}} \mathcal {S}_1^{(10)} f^{(0)} \\ 0 \\ 0 \end{pmatrix} + \begin{pmatrix} \left ( \mathcal {A}^{(10)} + {\textit{Kn}} \mathcal {A}^{(11)} \mathcal {S}_1^{(10)} \right ) f^{(0)} \\[3pt] \left ( \mathcal {A}^{(20)} + {\textit{Kn}} \mathcal {A}^{(21)} \mathcal {S}_1^{(10)} \right ) f^{(0)} \\[3pt] \left ( \mathcal {A}^{(\mathrm {r}0)} + {\textit{Kn}} \mathcal {A}^{(\mathrm {r}1)} \mathcal {S}_1^{(10)} \right ) f^{(0)} \end{pmatrix}\notag \\ = \frac {1}{\textit{Kn}} \begin{pmatrix} \mathcal {L}^{(11)} & \mathcal {L}^{(12)} & \mathcal {L}^{(1\mathrm {r})} \\ \mathcal {L}^{(21)} & \mathcal {L}^{(22)} & \mathcal {L}^{(2\mathrm {r})} \\ \mathcal {L}^{(\mathrm {r}1)} & \mathcal {L}^{(\mathrm {r}2)} & \mathcal {L}^{(\mathrm {rr})} \end{pmatrix} \begin{pmatrix} f_2^{(1)} \\[2pt] f_2^{(2)} \\[2pt] f_2^{(\mathrm {r})} \end{pmatrix}. \end{align}

Next, we perform the following operations.

1. (i) Replace the time derivative with the spatial derivative using

   (4.68) \begin{equation} \frac {\partial f^{(0)}}{\partial t} + \mathcal {A}^{(00)} f^{(0)} = O(\textit{Kn}), \end{equation}
   where the $O(\textit{Kn})$ part can be discarded and it does not affect the order of accuracy.

2. (ii) Apply Gaussian elimination and (4.65) to obtain

   (4.69) \begin{align} \begin{pmatrix} \left [\mathcal {A}^{(10)} + {\textit{Kn}} \mathcal {A}^{(11)} \mathcal {S}_1^{(10)} - {\textit{Kn}} \mathcal {S}_1^{(10)} \mathcal {A}^{(00)}\right ] f^{(0)} \\[2pt] \textit{Kn} \left [ \mathcal {A}^{(21)} \mathcal {S}_1^{(10)} - \mathcal {L}^{(21)} (\mathcal {L}^{(11)})^{-1} \left ( \mathcal {A}^{(11)} \mathcal {S}_1^{(10)} - \mathcal {S}_1^{(10)} \mathcal {A}^{(00)} \right ) \right ] f^{(0)} \\[5pt] \textit{Kn} \left [ \mathcal {A}^{(\mathrm {r}1)} \mathcal {S}_1^{(10)} - \mathcal {L}^{(\mathrm {r}1)} (\mathcal {L}^{(11)})^{-1} \left ( \mathcal {A}^{(11)} \mathcal {S}_1^{(10)} - \mathcal {S}_1^{(10)} \mathcal {A}^{(00)} \right ) \right ] f^{(0)} \\ \end{pmatrix} \notag \\ = \frac {1}{\textit{Kn}} \begin{pmatrix} \mathcal {L}^{(11)} & \mathcal {L}^{(12)} & \mathcal {L}^{(1\mathrm {r})} \\ & \mathcal {L}^{(22)}_* & \mathcal {L}^{(2\mathrm {r})}_* \\ & \mathcal {L}^{(\mathrm {r}2)}_* & \mathcal {L}^{(\mathrm {rr})}_* \end{pmatrix} \begin{pmatrix} f_2^{(1)} \\[2pt] f_2^{(2)} \\[2pt] f_2^{(\mathrm {r})} \end{pmatrix}. \end{align}

3. (iii) Apply Gaussian elimination again to eliminate $\mathcal {L}_*^{(\mathrm {r}2)}$ . Then the last equation in the system becomes

   (4.70) \begin{equation} {\textit{Kn}} \mathcal {B}^{(\mathrm {r0})} f^{(0)} = \frac {1}{\textit{Kn}} \left ( \mathcal {L}_*^{(\mathrm {rr})} - \mathcal {L}_*^{(\mathrm {r2})} (\mathcal {L}_*^{(22)})^{-1} \mathcal {L}_*^{(\mathrm {2r})} \right ) f_2^{(\mathrm {r})}, \end{equation}
   with
   (4.71) \begin{align} \mathcal {B}^{(\mathrm {r}0)} &= \left ( \mathcal {A}^{(\mathrm {r1})} - \mathcal {L}_*^{(\mathrm {r2})} (\mathcal {L}_*^{(22)})^{-1} \mathcal {A}^{(\mathrm {21})} \right ) \mathcal {S}_1^{(10)} \notag \\ & \quad - \left ( \mathcal {L}^{(\mathrm {r1})} - \mathcal {L}_*^{(\mathrm {r2})} (\mathcal {L}_*^{(22)})^{-1} \mathcal {L}^{(\mathrm {21})} \right ) (\mathcal {L}^{(11)})^{-1} \left ( \mathcal {A}^{(11)} \mathcal {S}_1^{(10)} - \mathcal {S}_1^{(10)} \mathcal {A}^{(00)} \right ). \end{align}
   Due to the fact that $f_2^{(\mathrm {r})} \sim o(\textit{Kn}^2)$ , the operator $\mathcal {B}^{(\mathrm {r}0)}$ has to be zero, and thus $f_2^{(\mathrm {r})} = 0$ , and $f_2^{(1)}$ and $f_2^{(2)}$ can be solved from (4.69). The results are
   (4.72) \begin{align} f_2^{(1)} &= {\textit{Kn}} \mathcal {S}_1^{(10)} f^{(0)} - {\textit{Kn}}^2 (\mathcal {L}^{(11)})^{-1} \mathcal {L}^{(12)}(\mathcal {L}_*^{(22)})^{-1} \mathcal {A}^{(11)} \mathcal {S}_1^{(10)} f^{(0)} + {\textit{Kn}}^2 \left(\mathcal {L}^{(11)}\right)^{-1} \nonumber \\ & \quad\times \left ( \mathcal {I} - \mathcal {L}^{(12)} (\mathcal {L}_*^{(22)})^{-1} \mathcal {L}^{(21)} (\mathcal {L}^{(11)})^{-1} \right ) \left ( \mathcal {A}^{(11)} \mathcal {S}_1^{(10)} - \mathcal {S}_1^{(10)} \mathcal {A}^{(00)} \right ) f^{(0)}, \end{align}
   (4.73) \begin{align} f_2^{(2)} &= {\textit{Kn}}^2 (\mathcal {L}_*^{(22)})^{-1}\left [ \mathcal {A}^{(21)} \mathcal {S}_1^{(10)} - \mathcal {L}^{(21)} (\mathcal {L}^{(11)})^{-1} \left ( \mathcal {A}^{(11)} \mathcal {S}_1^{(10)} - \mathcal {S}_1^{(10)} \mathcal {A}^{(00)} \right ) \right ] f^{(0)}. \end{align}

For conciseness, below we write these equations in the following simpler form:

(4.74) \begin{equation} f_2^{(1)} = {\textit{Kn}} \mathcal {S}_1^{(10)} f^{(0)} + {\textit{Kn}}^2 \mathcal {S}_2^{(10)} f^{(0)}, \qquad f_2^{(2)} = {\textit{Kn}}^2 \mathcal {S}_2^{(20)} f^{(0)}. \end{equation}

4.4.3. Third Maxwellian iteration

Similarly, the third Maxwellian iteration requires to find $f_3^{(1)}$ , $f_3^{(2)}$ and $f_3^{(3)}$ from (4.60) with $j = 2$ . However, for the purpose of deriving super-Burnett equations, only the first equation is needed. Here we rewrite this equation below:

(4.75) \begin{align} &\frac {\partial }{\partial t} \left ({\textit{Kn}} \mathcal {S}_1^{(10)} f^{(0)} + {\textit{Kn}}^2 \mathcal {S}_2^{(20)} f^{(0)} \right ) \nonumber \\ &\qquad + \mathcal {A}^{(10)} f^{(0)} + {\textit{Kn}} \mathcal {A}^{(11)} \mathcal {S}_1^{(10)}f^{(0)} + {\textit{Kn}}^2 \left ( \mathcal {A}^{(11)} \mathcal {S}_2^{(10)} + \mathcal {A}^{(12)} \mathcal {S}_2^{(20)} \right ) f^{(0)} \nonumber \\ &\qquad \qquad \qquad = \frac {1}{\textit{Kn}} \left(\mathcal {L}^{(11)} f_3^{(1)} + \mathcal {L}^{(12)} f_3^{(2)} + \mathcal {L}^{(1\mathrm {r})} f_3^{(\mathrm {r})} \right). \end{align}

Further simplification again requires the time derivatives to be replaced without affecting the order of magnitude. This can be done by using

(4.76) \begin{align} & \frac {\partial }{\partial t} \left ({\textit{Kn}} \mathcal {S}_1^{(10)} f^{(0)} + {\textit{Kn}}^2 \mathcal {S}_2^{(20)} f^{(0)} \right ) = - {\textit{Kn}} \mathcal {S}_1^{(10)} \mathcal {A}^{(00)} f^{(0)} \notag\\ & \qquad- {\textit{Kn}}^2 (\mathcal {S}_1^{(10)} \mathcal {A}^{(01)} \mathcal {S}_1^{(10)} + \mathcal {S}_2^{(20)} \mathcal {A}^{(00)}) f^{(0)} + O({\textit{Kn}}^3), \end{align}

which leads to

(4.77) \begin{align} & \mathcal {L}^{(11)} f_3^{(1)} + \mathcal {L}^{(12)} f_3^{(2)} + \mathcal {L}^{(1\mathrm {r})} f_3^{(\mathrm {r})} = {\textit{Kn}} \mathcal {A}^{(10)} f^{(0)} + {\textit{Kn}}^2 \left (\mathcal {A}^{(11)} \mathcal {S}_1^{(10)} - \mathcal {S}_1^{(10)} \mathcal {A}^{(00)} \right ) f^{(0)} \notag \\ & \quad+ {\textit{Kn}}^3 \left ( \mathcal {A}^{(11)} \mathcal {S}_2^{(10)} + \mathcal {A}^{(12)} \mathcal {S}_2^{(20)} - \mathcal {S}_2^{(20)} \mathcal {A}^{(00)} - \mathcal {S}_1^{(10)} \mathcal {A}^{(01)} \mathcal {S}_1^{(10)} \right ) f^{(0)}. \end{align}

4.4.4. Super-Burnett equations

By the adjointness (4.58) and the equalities (4.65), we have

(4.78) \begin{equation} \mathcal {A}^{(02)} = \mathcal {A}^{(01)} (\mathcal {L}^{(11)})^{-1} \mathcal {L}^{(12)}, \qquad \mathcal {A}^{(0\mathrm {r})} = \mathcal {A}^{(01)} (\mathcal {L}^{(11)})^{-1} \mathcal {L}^{(1\mathrm {r})}. \end{equation}

Thus, the super-Burnett equations (4.61) can be reformulated as

(4.79) \begin{equation} \frac {\partial f^{(0)}}{\partial t} + \mathcal {A}^{(00)} f^{(0)} + \mathcal {A}^{(01)} (\mathcal {L}^{(11)})^{-1} \left (\mathcal {L}^{(11)}f_3^{(1)} + \mathcal {L}^{(12)} f_3^{(2)} + \mathcal {L}^{(1\mathrm {r})} f_3^{(\mathrm {r})} \right ) = 0. \end{equation}

The final super-Burnett equations can be found by plugging (4.77) into the above equation.

Such a derivation of super-Burnett equations shows that the final equations depend only on the operators

(4.80) \begin{equation} \mathcal {A}^{(00)}, \mathcal {A}^{(01)}, \mathcal {A}^{(02)}, \mathcal {A}^{(10)}, \mathcal {A}^{(11)}, \mathcal {A}^{(12)}, \mathcal {A}^{(20)}, \mathcal {A}^{(21)}, \mathcal {L}^{(11)}, \mathcal {L}^{(12)}, \mathcal {L}^{(21)}, \mathcal {L}^{(22)}, \end{equation}

and the operators

(4.81) \begin{equation} \mathcal {A}^{(0\mathrm {r})}, \mathcal {A}^{(1\mathrm {r})}, \mathcal {A}^{(22)}, \mathcal {A}^{(2\mathrm {r})}, \mathcal {A}^{(\mathrm {r}0)}, \mathcal {A}^{(\mathrm {r}1)}, \mathcal {A}^{(\mathrm {r}2)}, \mathcal {A}^{(\mathrm {rr})}, \mathcal {L}^{(1\mathrm {r})}, \mathcal {L}^{(2\mathrm {r})}, \mathcal {L}^{(\mathrm {r}1)}, \mathcal {L}^{(\mathrm {r}2)}, \mathcal {L}^{(\mathrm {rr})} \end{equation}

do not appear in the final equations. Therefore, we can set these operators to be zero, and the resulting equations still have the super-Burnett order. Meanwhile, the time derivative $\partial _t f_2$ does not participate in the derivation, which can also be set to zero without affecting the super-Burnett order.

5.1. Onsager boundary conditions

The abstract form of the R13 equations (4.33) shows that all operators in the Boltzmann equation are approximated by operators in $\mathbb {V} = \mathbb {V}^{(0)} \oplus \mathbb {V}^{(1)} \oplus \mathbb {V}^{(2)}$ . If we rewrite the Boltzmann equation (1.1) as

(5.9) \begin{equation} \frac {\partial }{\partial t} (\mathcal {I} f) + \frac {\partial }{\partial x_k} (\xi _k f) = \frac {1}{\textit{Kn}} \mathcal {L}f, \end{equation}

then the identity operator $\mathcal {I}$ is approximated by $\mathcal {P}^{(01)} := \mathcal {P}^{(0)} + \mathcal {P}^{(1)}$ , the velocity operator $\xi _k$ is approximated by

(5.10) \begin{equation} \mathcal {A}_k = \mathcal {P}^{(01)} \xi _k \mathcal {P}^{(01)} + \mathcal {P}^{(2)} \xi _k \mathcal {P}^{(01)} + \mathcal {P}^{(01)} \xi _k \mathcal {P}^{(2)} \end{equation}

and the collision operator $\mathcal {L}$ is approximated by $(\mathcal {P}^{(1)} + \mathcal {P}^{(2)}) \mathcal {L}(\mathcal {P}^{(1)} + \mathcal {P}^{(2)})$ . Similarly, when deriving boundary conditions for moment equations, operators in (5.8) are also approximated by using operators on $\mathbb {V}$ .

We first study the operator $\mathcal {P}_{\textit{even}}$ , which extracts the even part of a function with respect to $\xi _n$ . The discretisation of $\mathcal {P}_{\textit{even}}$ is simply its restriction on $\mathbb {V}$ . Note that all the basis functions of $\mathbb {V}$ are given in (4.35), from which it is not difficult to observe that

(5.11) \begin{align} & \mathbb {V}_{\textit{even}} := \mathcal {P}_{\textit{even}} \mathbb {V} = \operatorname {span} \{ \psi ^0, \psi ^1, \psi _{\tau _1}^0, \psi _{\tau _2}^0; \phi _{\tau _1}^1, \phi _{\tau _2}^1, \phi _{\tau _1 \tau _1}^0, \phi _{\tau _2 \tau _2}^0, \phi _{\tau _1 \tau _2}^0; \phi ^2, \phi _{\tau _1}^2, \phi _{\tau _2}^2,\notag\\ & \quad \phi _{\tau _1}^3, \phi _{\tau _2}^3, \phi _{\tau _1 \tau _1}^1, \phi _{\tau _2 \tau _2}^1, \phi _{\tau _1 \tau _2}^1, \phi _{\tau _1 \tau _1}^2, \phi _{\tau _2 \tau _2}^2, \phi _{\tau _1 \tau _2}^2, \phi _{\tau _1 \tau _1 \tau _1}^0, \phi _{\tau _1 \tau _1 \tau _2}^0,\phi _{\tau _1 \tau _2 \tau _2}^0,\phi _{\tau _2 \tau _2 \tau _2}^0\}. \end{align}

It is clear that $f_W \in \mathbb {V}_{\textit{even}}$ . Thus, the velocity operator $\xi _n$ should be discretised as the restriction of $\mathcal {A}_n$ on $\mathbb {V}_{\textit{even}}$ . We denote its range by $\mathbb {V}_{\textit{odd}} := \mathcal {A}_n\mathbb {V}_{\textit{even}}$ . By straightforward calculation,

(5.12) \begin{align} \mathbb {V}_{\textit{odd}} &= \operatorname {span}\Big \{ \psi _n^0; \phi ^1_n,\phi _{n\tau _1}^{0}, \phi _{n\tau _2}^{0}; \phi _{n\tau _1}^{1}, \phi _{n\tau _2}^{1},\phi _{n\tau _1}^{2}, \phi _{n\tau _2}^{2}, \phi _{n\tau _1 \tau _2}^{0}, \phi _{n\tau _1 \tau _1}^{0} + \frac {1}{2} \phi _{nnn}^{0}, \phi _n^{2} \notag\\ &\quad + \frac {c^{3,1}_1}{c^{2,1}_1} \phi _n^{3},\mu _1 \phi _{nnn}^{0} + \mu _2 \left ( \phi _n^{3} - \frac {c^{3,1}_1}{c^{2,1}_1} \phi _n^{2} \right )\Bigg \}, \end{align}

where

(5.13) \begin{equation} \mu _1 = 3 A_{5,11} \left [ 1 + \left ( \frac {c_1^{3,1}}{c_1^{2,1}} \right )^2 \right ], \qquad \mu _2 = 2 \left ( A_{58} - \frac {c_1^{3,1}}{c_1^{2,1}} A_{57} \right ). \end{equation}

Note that the basis functions listed in (5.12) are orthogonal basis functions, and the orthogonality will simplify our derivations of boundary conditions. The discretisation the operator $\mathcal {S}$ is then a straightforward application of Galerkin’s method with the subspace $\mathbb {V}_{\textit{odd}}$ .

We would like to comment that $\mathbb {V}_{\textit{odd}} \oplus \mathbb {V}_{\textit{even}}$ is only a proper subspace of $\mathbb {V}$ , or equivalently, $\mathbb {V}_{\textit{odd}} \neq \mathcal {P}_{\textit{odd}} \mathbb {V}$ . The reason is that the operator $\mathcal {A}_n$ has a null space being a subspace of $\mathcal {P}_{\textit{odd}} \mathbb {V}$ . This can be seen by focusing on the $x_n$ derivatives in (4.44), (4.45) and (4.48), and find that by the linear combination

(5.14) \begin{equation} A_{5,11} \left (c_1^{3,1} \times (4.29) - c_1^{2,1} \times (4.30) \right ) - \frac {5}{3} \left(c_1^{3,1} A_{57} - c_1^{2,1} A_{58} \right) \times (4.33) \end{equation}

the $x_n$ derivative will be cancelled out. This indicates that

(5.15) \begin{equation} \mathcal {A}_n \left ( 3A_{5,11} \left (c_1^{3,1} \phi _n^2 - c_1^{2,1} \phi _n^3 \right ) - \frac {175}{6} \left(c_1^{3,1} A_{57} - c_1^{2,1} A_{58}\right) \phi _{nnn}^0 \right ) = 0, \end{equation}

which reveals the null space of $\mathcal {A}_n$ . The function space $\mathbb {V}_{\textit{odd}}$ defined in (5.12) is orthogonal to the null space of $\mathcal {A}_n$ .

Practically, all the boundary conditions are written as

(5.16) \begin{equation} \langle \phi _{\textit{odd}}^i, f_{\textit{odd}} \rangle = \frac {2\chi }{2-\chi } \sum _j \frac {\langle \phi _{\textit{odd}}^i, \mathcal {C} \xi _n^{-1} \phi _{\textit{odd}}^j \rangle \langle \phi _{\textit{odd}}^j, \mathcal {A}_n(f_W - f_{\textit{even}}) \rangle }{\langle \phi _{\textit{odd}}^j, \phi _{\textit{odd}}^j \rangle }, \quad j = 1,\ldots , 12, \end{equation}

where $\phi _{\textit{odd}}^j$ refers to the $j$ th basis function of $\mathbb {V}_{\textit{odd}}$ (see (5.12)), and $f_{\textit{odd}}$ and $f_{\textit{even}}$ are, respectively, the projection of $f$ onto the function space $\mathbb {V}_{\textit{odd}}$ and $\mathbb {V}_{\textit{even}}$ . Note that this formula requires the orthogonality of the basis $\{\phi _{\textit{odd}}^j\}$ . Calculation of the integrals in (5.16) requires specifying the collision model. Here we only provide the general form

(5.17) \begin{equation} \begin{pmatrix} v_n \\ \bar {q}_n \\ u_n^2 + \frac {c_1^{3,1}}{c_1^{2,1}} u_n^3 \\ \bar {\sigma }_{\tau _1 n} \\ \bar {\sigma }_{\tau _2 n} \\ u_{\tau _1 n}^1 \\ u_{\tau _2 n}^1 \\ u_{\tau _1 n}^2 \\ u_{\tau _2 n}^2 \\ u_{nnn}^0 + \frac {\mu _2}{\mu _1} \left ( u_n^3 - \frac {c_1^{3,1}}{c_1^{2,1}} u_n^2 \right ) \\ u_{\tau _1 \tau _1 n}^0 + \frac {1}{2} u_{nnn}^0 \\ u_{\tau _1 \tau _2 n}^0 \end{pmatrix} = \frac {2\chi }{2-\chi } Z \begin{pmatrix} V_n \\ \bar {Q}_n \\ U_n^2 + \frac {c_1^{3,1}}{c_1^{2,1}} U_n^3 \\ \bar {\Sigma }_{\tau _1 n} \\ \bar {\Sigma }_{\tau _2 n} \\ U_{\tau _1 n}^1 \\ U_{\tau _2 n}^1 \\ U_{\tau _1 n}^2 \\ U_{\tau _2 n}^2 \\ U_{nnn}^0 + \frac {\mu _2}{\mu _1} \left ( U_n^3 - \frac {c_1^{3,1}}{c_1^{2,1}} U_n^2 \right ) \\ U_{\tau _1 \tau _1 n}^0 + \frac {1}{2} U_{nnn}^0 \\ U_{\tau _1 \tau _2 n}^0 \end{pmatrix}, \end{equation}

where the capitalised variables are given by

(5.18) \begin{align} V_n = &\ (\rho ^W - \rho ) + (\theta ^W - \theta ) - c_2^{0,0} \bar {\sigma }_{nn} - \sqrt {15} (c_2^{1,0} u_{nn}^1 + c_2^{2,0} u_{nn}^2), \notag\\ \bar {Q}_{n} = &\ \frac {5}{2} c^{1,1}_1(\theta ^W - \theta ) + \frac {\sqrt {2}}{6}A_{45}\bar {\sigma }_{nn}+\sqrt {\frac {5}{6}} (A_{46}u^{2}+A_{49}u^{1}_{nn} + A_{4,10}u^{2}_{nn}),\notag \\ U^{2}_{n} = &\ \frac {\sqrt {15}}{3}\left (\frac {\sqrt {2}}{2} c^{2,1}_1 (\theta - \theta ^W) - \frac {1}{15} A_{57}\bar {\sigma }_{nn} \right ),\notag\\ U^{3}_{n} = &\ \frac {\sqrt {15}}{3}\left (\frac {\sqrt {2}}{2} c^{3,1}_1 (\theta - \theta ^W) - \frac {1}{15} A_{58}\bar {\sigma }_{nn} \right ),\notag\\ \bar {\Sigma }_{\tau _in} = &\ c^{0,0}_2 \!\left(\!v^W_{\tau _i}-v_{\tau _i}\!\right) + \frac {\sqrt {2}}{15} A_{45} \bar {q}_{\tau _i} - \frac {1}{\sqrt {15}}\! \left(\!A_{57}u^{2}_{\tau _i} + A_{58}u^{3}_{\tau _i} + 2 A_{5,11} u^{0}_{\tau _inn}\!\right)\!, \!\!\!\quad i=1,2, \notag\\ U^{1}_{\tau _in} = &\ \frac {\sqrt {15}}{15} \left ( c^{1,0}_2 (v^W_{\tau _i}-v_{\tau _i}) + \frac {\sqrt {2}}{15} A_{49}\bar {q}_{\tau _i} \right ),\quad i=1,2,\notag \\ U^{2}_{\tau _in} = &\ \frac {\sqrt {15}}{15} \left ( c^{2,0}_2 (v^W_{\tau _i}-v_{\tau _i}) + \frac {\sqrt {2}}{15} A_{4,10}\bar {q}_{\tau _i} \right ),\quad i=1,2,\notag \\ U^{0}_{\tau _i\tau _jn} = &\ \frac {2\sqrt {15}}{1575} A_{5,11}\left ( \frac {2}{5} \delta _{ij} \bar {\sigma }_{nn} - \bar {\sigma }_{\tau _i\tau _j}\right ), \quad i,j=1,2, \notag\\ U^{0}_{nnn} = &\,{-\frac {2\sqrt {15}}{875}}A_{5,11}\bar {\sigma }_{nn}, \end{align}

and the structure of the matrix $Z$ is (“ $*$ ” means a non-zero entry)

(5.19) \begin{equation} Z = \left ( \begin{array}{@{}cccccccccccc@{}} * &\quad * &\quad * &\quad 0 &\quad 0 &\quad 0 &\quad 0 &\quad 0 &\quad 0 &\quad * &\quad * &\quad 0 \\ * & * & * & 0 & 0 & 0 & 0 & 0 & 0 & * & * & 0 \\ * & * & * & 0 & 0 & 0 & 0 & 0 & 0 & * & * & 0 \\ 0 & 0 & 0 & * & 0 & * & 0 & * & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & * & 0 & * & 0 & * & 0 & 0 & 0 \\ 0 & 0 & 0 & * & 0 & * & 0 & * & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & * & 0 & * & 0 & * & 0 & 0 & 0 \\ 0 & 0 & 0 & * & 0 & * & 0 & * & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & * & 0 & * & 0 & * & 0 & 0 & 0 \\ * & * & * & 0 & 0 & 0 & 0 & 0 & 0 & * & * & 0 \\ * & * & * & 0 & 0 & 0 & 0 & 0 & 0 & * & * & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & * \end{array} \right ) = \begin{pmatrix} z_0 &\quad \zeta _0^{\intercal } \\ \zeta _1 & Z_1 \end{pmatrix}. \end{equation}

The matrix $Z$ is the discrete form of the operator $\mathcal {S} = \mathcal {P}_{\textit{odd}} \mathcal {C} \xi _n^{-1}$ , and is therefore similar to a symmetric and positive definite matrix. To complete the derivation, we still need to take into account the condition $v_n = \langle \xi _n, f \rangle = 0$ , since $V_n$ depends on $\rho _W$ (see (5.3)), which has to be determined by this condition. Instead of solving $\rho _W$ , an equivalent way is to solve $V_n$ from the first equation of (5.17), and then plug the result into other equations. Such an equivalence is due to the fact that $\rho ^W$ appears only in $V_n$ in all the quantities defined in (5.18). This leads to boundary conditions that have the same form as (5.17), but the matrix $Z$ is replaced by

(5.20) \begin{equation} \begin{pmatrix} 0 & 0 \\ 0 & Z_1 - z_0^{-1} \zeta _1 \zeta _0^{\intercal } \end{pmatrix}. \end{equation}

Since the Schur complement of a symmetric and positive definite matrix is still symmetric and positive definite, we can guarantee that the above matrix can still be viewed as a positive semidefinite operator.

Such boundary conditions satisfy all the properties required by the $L^2$ stability. However, spurious boundary layers are observed in the solution, which will be detailed in the next subsection.

5.2. One-dimensional problems and spurious boundary layers

To show the boundary layers in the solution, we consider the steady-state one-dimensional flow in which

(5.21) \begin{equation} v_2 = v_3 = \bar {q}_2 = \bar {q}_3 = \bar {\sigma }_{13} = \bar {\sigma }_{23}= 0, \qquad \bar {\sigma }_{22} = \bar {\sigma }_{33}. \end{equation}

These conditions come from the assumption that the distribution function satisfies $f(\boldsymbol {x}, \xi _1, \xi _2, \xi _3, t) = \tilde {f}(\boldsymbol {x}, \xi _1, \sqrt {\xi _2^2 + \xi _3^2}, t)$ for every $\boldsymbol {x}$ and $t$ , meaning that $f$ is axisymmetric about the axis $\xi _2 = \xi _3 = 0$ . Under this assumption, only five variables remain and the equations have the form

(5.22) \begin{align} & \frac {\partial v}{\partial x} = 0, \end{align}

(5.23) \begin{align} & \frac {2}{3} \frac {\partial v}{\partial x} +\frac {2}{3}k_0\frac {\partial \bar q}{\partial x}-k_1 {\textit{Kn}} \frac {\partial ^2 \theta }{\partial x^2}+k_2 {\textit{Kn}}\frac {\partial ^2\bar {\sigma }}{\partial x^2} = 0, \end{align}

(5.24) \begin{align} & \frac {\partial \rho }{\partial x} +\frac {\partial \theta }{\partial x} - \frac {2}{3} k_3 {\textit{Kn}} \frac {\partial ^2 v}{\partial x^2} - \frac {2}{3} k_4 {\textit{Kn}} \frac {\partial ^2\bar {q}}{\partial x^2}+ k_5 \frac {\partial \bar {\sigma }}{\partial x} = 0, \end{align}

(5.25) \begin{align} & \frac {5}{2} k_0\frac {\partial \theta }{\partial x} -\frac {5}{3} k_4 {\textit{Kn}}\frac {\partial ^2v}{\partial x^2} -2 k_6 {\textit{Kn}}\frac {\partial ^2 \bar {q}}{\partial x^2}-\frac {8}{5}k_7 {\textit{Kn}}\frac {\partial ^2 \bar {q}}{\partial x^2}+k_8\frac {\partial \bar {\sigma }}{\partial x} = - \frac {2}{3} \frac {1}{\textit{Kn}}l_1\bar {q}, \end{align}

(5.26) \begin{align} & 2k_2 {\textit{Kn}}\frac {\partial ^2\theta }{\partial x^2}+ \frac {4}{3} k_5 \frac {\partial v}{\partial x} +\frac {8}{15}k_8\frac {\partial \bar {q}}{\partial x} - \frac {6}{5} k_9 {\textit{Kn}}\frac {\partial ^2 \bar {\sigma }}{\partial x^2} -\frac {2}{3} k_{10} {\textit{Kn}}\frac {\partial ^2 \bar {\sigma }}{\partial x^2}= - \frac {1}{\textit{Kn}}l_2\bar {\sigma }. \end{align}

Here time derivatives have been removed so that the solution of the equations correspond to the steady state of the fluid. In the one-dimensional case, the boundary conditions are reduced to

(5.27) \begin{align} & v = 0, \end{align}

(5.28) \begin{align} \bar {q} = &\, s_n \tilde {\chi } \left [ \hat {m}_{11} (\theta -\theta ^{W}) +\hat {m}_{12} \bar {\sigma } - \hat {m}_{13} {\textit{Kn}} \frac {\partial \bar {q}}{\partial x} - \frac {2}{3} \hat {m}_{14} {\textit{Kn}} \frac {\partial \bar {q}}{\partial x} - \frac {2}{3} \hat {m}_{15} {\textit{Kn}} \frac {\partial v}{\partial x} \right ], \\ & \quad\quad \hat {m}_{26}\bar {q} + \hat {m}_{27} {\textit{Kn}} \frac {\partial \theta }{\partial x} - \hat {m}_{28} {\textit{Kn}} \frac {\partial \bar {\sigma }}{\partial x} \nonumber \end{align}

(5.29) \begin{align} & \quad \!=s_n \tilde {\chi } \left [-\hat {m}_{21} (\theta -\theta ^{W}) + \hat {m}_{22} \bar {\sigma } + \hat {m}_{23} {\textit{Kn}} \frac {\partial \bar {q}}{\partial x} + \frac {2}{3} \hat {m}_{24} {\textit{Kn}} \frac {\partial \bar {q}}{\partial x} + \frac {2}{3} \hat {m}_{25} {\textit{Kn}} \frac {\partial v}{\partial x} \right ], \\ & \qquad \qquad {-\hat {m}_{66}} \bar {q} - {\textit{Kn}}\left (\frac {3}{5} \hat {m}_{67} \frac {\partial \bar {\sigma }}{\partial x} +\hat {m}_{68}\frac {\partial \bar {\sigma }}{\partial x} + \hat {m}_{69}\frac {\partial \theta }{\partial x} \right ) \nonumber \end{align}

(5.30) \begin{align} & \quad \! = s_n \tilde {\chi } \left [-\hat {m}_{61} (\theta -\theta ^{W}) +\hat {m}_{62} \bar {\sigma } + \hat {m}_{63} {\textit{Kn}} \frac {\partial \bar {q}}{\partial x} + \frac {2}{3} \hat {m}_{64} {\textit{Kn}} \frac {\partial \bar {q}}{\partial x} + \frac {2}{3} \hat {m}_{65} {\textit{Kn}} \frac {\partial v}{\partial x}\right ], \end{align}

where $s_n=1$ at the right boundary and $s_n=-1$ at the left boundary. The coefficient $\hat {m}_{ij}$ is computed from the Onsager boundary conditions (5.17), which are not equal to the $m_{ij}$ in our final boundary conditions, as shown in (3.9)–(3.17).

The one-dimensional equations (5.22)–(5.30) can be solved analytically. First, it is straightforward to get $v(x)=0$ from (5.22) and (5.27). Second, by noting that variable $\rho$ only appears in (5.24), we can solve other equations and then plug other variables into (5.24). Third, (5.23) reveals that $ {\partial \bar {q}}/{\partial x} = ({3}/{2 k_0}) {\textit{Kn}} ( k_1 ({\partial ^2 \theta }/{\partial x^2})) - (k_2 ({\partial ^2 \bar {\sigma }}/{\partial x^2} ))$ , which indicates that

(5.31) \begin{equation} \bar {q}(x) = \frac {3}{2 k_0} {\textit{Kn}} \left ( k_1 \frac {\partial \theta }{\partial x} - k_2 \frac {\partial \bar {\sigma }}{\partial x} \right ) + C_{\bar {q}}, \end{equation}

with the constant $C_{\bar {q}}$ to be determined. Plugging $\bar {q}(x)$ into (5.25) and integrating the equation, we find that

(5.32) \begin{equation} {\textit{Kn}}^2 \left (\! -k_6-\frac {4}{5} k_7 \!\right ) \frac {3}{k_0} \left (\! k_1 \frac {\partial ^2 \theta }{\partial x^2} - k_2 \frac {\partial ^2 \bar {\sigma }}{\partial x^2} \!\right ) + \frac {5}{2} k_0 \theta + k_8 \bar {\sigma } = - \frac {l_1}{k_0} \left ( k_1 \theta - k_2 \bar {\sigma } \!\right ) + C_{\bar {q}}^1 x + C_{\bar {q}}^2, \end{equation}

where $C_{\bar {q}}^1 = -( {2 l_1}/{3 {\textit{Kn}}}) C_{\bar {q}}$ and $C_{\bar {q}}^2$ are two constants to be determined. Similarly, plugging (5.31) into (5.26) yields

(5.33) \begin{equation} 2 k_2 {\textit{Kn}}\frac {\partial ^2\theta }{\partial x^2} + \frac {4 k_8 }{5 k_0} {\textit{Kn}} \left ( k_1 \frac {\partial ^2 \theta }{\partial x^2} - k_2 \frac {\partial ^2 \bar {\sigma }}{\partial x^2} \right ) - \left ( \frac {6}{5} k_9 + \frac {2}{3} k_{10} \right ) {\textit{Kn}}\frac {\partial ^2 \bar {\sigma }}{\partial x^2}= - \frac {1}{\textit{Kn}}l_2\bar {\sigma }. \end{equation}

We can then define

(5.34) \begin{equation} \xi = {\textit{Kn}} \frac {\partial \theta }{\partial x}, \qquad \xi ' = {\textit{Kn}} \frac {\partial \bar {\sigma }}{\partial x}, \end{equation}

and reformulate (5.32) and (5.33) as

(5.35) \begin{align} {\textit{Kn}} \left(\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ - \left ( k_6 + \frac {4}{5} k_7 \right ) \frac {3 k_1}{k_0} & \left ( k_6 + \frac {4}{5} k_7 \right ) \frac {3 k_2}{k_0} & 0 & 0 \\ 2 k_2 + \frac {4 k_8 k_1}{5 k_0} & - \left ( \frac {6}{5} k_9 + \frac {2}{3} k_{10} + \frac {4 k_8 k_2}{5 k_0} \right ) & 0 & 0 \\ \end{array} \right) \frac {\partial }{\partial x} \begin{pmatrix} \xi \\ \xi ' \\ \theta \\ \bar {\sigma } \end{pmatrix} \nonumber \\ + \left(\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} -1 & 0 & 0 & 0 \\ 0 & -1 & 0 & 0 \\ 0 & 0 & \frac {5}{2}k_0 + \frac {l_1 k_1}{k_0} & k_8 - \frac {l_1 k_2}{k_0} \\ 0 & 0 & 0 & l_2 \\ \end{array} \right) \begin{pmatrix} \xi \\ \xi ' \\ \theta \\ \bar {\sigma } \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ C_{\bar {q}}^1 x + C_{\bar {q}}^2 \\ 0 \end{pmatrix}. \end{align}

By multiplying the inverse of the matrix in front of the first-order derivatives, the above equation has the following structure:

(5.36) \begin{equation} {\textit{Kn}} \frac {\partial }{\partial x} \begin{pmatrix} \xi \\ \xi ' \\ \theta \\ \bar {\sigma } \end{pmatrix} + \left(\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} 0 & 0 & b_{11} & b_{12} \\ 0 & 0 & b_{21} & b_{22} \\ -1 & 0 & 0 & 0 \\ 0 & -1 & 0 & 0 \end{array} \right) \begin{pmatrix} \xi \\ \xi ' \\ \theta \\ \bar {\sigma } \end{pmatrix} = \begin{pmatrix} c_1 ( C_{\bar {q}}^1 x + C_{\bar {q}}^2 ) \\ c_2 ( C_{\bar {q}}^1 x + C_{\bar {q}}^2 ) \\ 0 \\ 0 \end{pmatrix}. \end{equation}

In general, given the following ordinary differential equation (ODE) system

(5.37) \begin{equation} {\textit{Kn}} \frac {\partial g}{\partial x} + A g = r(x), \end{equation}

its solution is

(5.38) \begin{equation} g(x) = e^{-A x / {\textit{Kn}}} g(0) + \frac {1}{\textit{Kn}} \int _0^x e^{A(s-x)/ {\textit{Kn}}} r(s) \mathrm {d}s. \end{equation}

In our case, the matrix in (5.36) has four distinct eigenvalues $\pm \lambda _1$ and $\pm \lambda _2$ with

(5.39) \begin{equation} \lambda _j = \frac {\sqrt {(-1)^j\sqrt {b_{11}^2-2 b_{22} b_{11}+b_{22}^2+4 b_{12} b_{21}}-b_{11}-b_{22}}}{\sqrt {2}}, \qquad j=1,2. \end{equation}

The solution of (5.36) therefore has the form

(5.40) \begin{equation} \kappa _1 \sinh \left ( \frac {\lambda _1 x}{\textit{Kn}} \right ) + \kappa _2 \cosh \left ( \frac {\lambda _1 x}{\textit{Kn}} \right ) + \kappa _3 \sinh \left ( \frac {\lambda _2 x}{\textit{Kn}} \right ) + \kappa _4 \cosh \left ( \frac {\lambda _2 x}{\textit{Kn}} \right ) + \kappa _5 x + \kappa _6, \end{equation}

where the hyperbolic sines and cosines denote the Knudsen layers. Since $\lambda _1 \neq \lambda _2$ , two boundary layers will be present in the general solution. In particular, when the magnitude of the eigenvalues is large or the Knudsen number is small, the solution would exhibit boundary layers.

For Maxwell molecules ( $\eta = 5$ ), the above procedure of solving one-dimensional equations can be further simplified from the third step (see (5.31). By noting that $k_1 = k_2 = 0$ , one can get $\ {\partial \bar {q}}/{\partial x} = 0$ and therefore, (5.31) becomes $\bar {q}(x) = C_{\bar {q}}$ . Equation (5.33) becomes a second-order ODE for $\bar {\sigma }$ , so that the general solution of $\bar {\sigma }(x)$ contains only one boundary layer. Furthermore, (5.32) is in absence of second-order derivatives, and thus,

(5.41) \begin{equation} \theta (x) = \frac {2}{5k_0} \left (C_{\bar {q}}^1 x + C_{\bar {q}}^2 - k_8 \bar {\sigma }(x) \right ), \end{equation}

showing that the boundary layer of $\theta$ has the same thickness as that of $\bar {\sigma }$ . In general, when the type of molecules gets closer to Maxwell molecules, the second boundary layer will get thinner, which explains why there is only one boundary layer left for Maxwell molecules.

The thinness of the second boundary layer also indicates the largeness of the corresponding eigenvalue $\lambda _2$ . To illustrate the magnitude of eigenvalues, we compute them according to (5.39) for various inverse-power-law models and summarise the results in table 1. For Maxwell molecules, only $\lambda _1$ is provided. The table clearly shows that all models share one similar eigenvalue $\lambda _1$ and each of the non-Maxwell cases has an additional larger eigenvalue $\lambda _2$ . This eigenvalue $\lambda _2$ approaches infinity when $\eta$ tends to $5$ .

Table 1. Eigenvalues in the general solution (5.40) for some inverse-power-law models with power index $\eta$ .

The one-dimensional solutions on the domain $(-0.5, 0.5)$ , with Knudsen number $\textit{Kn}=0.2$ and accommodation coefficients $\chi =1$ , are depicted in figure 1. The figure verifies that the larger eigenvalue $\lambda _2$ in non-Maxwell cases induces additional boundary layers compared with the Maxwell solutions. These boundary layers are non-physical according to previous results (see Hu, Yang & Cai Reference Hu, Yang and Cai2020; Cai et al. Reference Cai, Torrilhon and Yang2024b ) for the same problem. In the next subsection we propose a revision to the Onsager boundary conditions derived in § 5.1 to remove the boundary layers, leading to our final boundary conditions (3.9)–(3.17).

Figure 1. Plots of $\bar {q}$ and $\theta$ for the one-dimensional problem (5.22)–(5.30). Results are shown for (a) $\theta ^W=0$ and (b) $\theta ^{W}=0.2$ .

5.3. Removal of undesired boundary layers

Our discussion will be based on the following general form of one-dimensional symmetric hyperbolic equations:

(5.42) \begin{equation} \begin{pmatrix} 0 & \mathbf {A}_{\textit{oe}} \\ \mathbf {A}_{\textit{eo}} & 0 \end{pmatrix} \frac {\partial }{\partial x} \begin{pmatrix} \boldsymbol{u}_{\textit{o}} \\ \boldsymbol{u}_{\textit{e}} \end{pmatrix} = \begin{pmatrix} \mathbf {L}_{\textit{oo}} & 0 \\ 0 & \mathbf {L}_{\textit{ee}} \end{pmatrix} \begin{pmatrix} \boldsymbol{u}_{\textit{o}} \\ \boldsymbol{u}_{\textit{e}} \end{pmatrix}. \end{equation}

Here $\mathbf {A}_{\textit{eo}} = \mathbf {A}_{\textit{oe}}^{\intercal }$ , and $\mathbf {L}_{\textit{oo}}$ and $\mathbf {L}_{\textit{ee}}$ are symmetric and negative semidefinite matrices. As noted in § 5.1, the operator $\mathcal {A}_n$ has a non-trivial null space being a subspace of $\mathcal {P}_{\textit{odd}} \mathbb {V}$ . Accordingly, in order to mimic this property, we assume that the matrix $\mathbf {A}_{\textit{oe}}$ is rank deficient, and its singular value decomposition is

(5.43) \begin{equation} \mathbf {A}_{\textit{oe}} = \begin{pmatrix} \mathbf {U}_1 & \mathbf {U}_2 \end{pmatrix} \begin{pmatrix} \boldsymbol {\Lambda }_r & 0 \\ 0 & 0 \end{pmatrix} \begin{pmatrix} \mathbf {V}_1^{\intercal } \\ \mathbf {V}_2^{\intercal } \end{pmatrix}, \end{equation}

where $r$ denotes the rank of $\mathbf {A}_{\textit{oe}}$ and

(5.44) \begin{equation} \boldsymbol {\Lambda }_r = \begin{pmatrix} \lambda _1 & 0 & \cdots & 0 \\ 0 & \lambda _2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \lambda _r \end{pmatrix}, \qquad \lambda _1 \geqslant \lambda _2 \geqslant \cdots \geqslant \lambda _r \gt 0. \end{equation}

Inspired by the preliminary section in Jiang, Sun & Toh (Reference Jiang, Sun and Toh2013), we define

(5.45) \begin{equation} \boldsymbol {v}_{\textit{o}} = \mathbf {U}_1^{\intercal } \boldsymbol{u}_{\textit{o}}, \quad \boldsymbol {v}_{\textit{e}} = \mathbf {V}_1^{\intercal } \boldsymbol{u}_{\textit{e}}. \end{equation}

By straightforward calculation, it can be shown that $\boldsymbol {v}_{\textit{o}}$ and $\boldsymbol {v}_{\textit{e}}$ satisfy the differential equations

(5.46) \begin{equation} \frac {\partial }{\partial x} \begin{pmatrix} \boldsymbol {v}_{\textit{o}} \\ \boldsymbol {v}_{\textit{e}} \end{pmatrix} = \begin{pmatrix} 0 & \boldsymbol {\Lambda }_r^{-1} \tilde {\mathbf {L}}_{\textit{ee}} \\ \boldsymbol {\Lambda }_r^{-1} \tilde {\mathbf {L}}_{\textit{oo}} & 0 \end{pmatrix} \begin{pmatrix} \boldsymbol {v}_{\textit{o}} \\ \boldsymbol {v}_{\textit{e}} \end{pmatrix}, \end{equation}

where

(5.47) \begin{align} \tilde {\mathbf {L}}_{\textit{oo}} &= \mathbf {U}_1^{\intercal } \mathbf {L}_{\textit{oo}}\mathbf {U}_1 - \mathbf {U}_1^{\intercal } \mathbf {L}_{\textit{oo}}\mathbf {U}_2 (\mathbf {U}_2^{\intercal } \mathbf {L}_{\textit{oo}}\mathbf {U}_2)^{-1}\mathbf {U}_2^{\intercal } \mathbf {L}_{\textit{oo}}\mathbf {U}_1, \nonumber \\ \tilde {\mathbf {L}}_{\textit{ee}} &= \mathbf {V}_1^{\intercal } \mathbf {L}_{\textit{ee}}\mathbf {V}_1 - \mathbf {V}_1^{\intercal } \mathbf {L}_{\textit{ee}}\mathbf {V}_2 (\mathbf {V}_2^{\intercal } \mathbf {L}_{\textit{ee}}\mathbf {V}_2)^{-1}\mathbf {V}_2^{\intercal } \mathbf {L}_{\textit{ee}}\mathbf {V}_1. \end{align}

It is obvious that $\tilde {\mathbf {L}}_{\textit{oo}}$ and $\tilde {\mathbf {L}}_{\textit{ee}}$ are Schur complements of negative semidefinite matrices, and therefore, are also negative semidefinite matrices. Note that here we require the invertibility of $\mathbf {U}_2^{\intercal } \mathbf {L}_{\textit{oo}}\mathbf {U}_2$ and $\mathbf {V}_2^{\intercal } \mathbf {L}_{\textit{ee}}\mathbf {V}_2$ , which is satisfied in the R13 equations derived in previous sections.

The following lemma reveals the eigenpairs of the matrix in (5.46).

Proof. From

(5.48) \begin{equation} \begin{pmatrix} 0 & \boldsymbol {\Lambda }_r^{-1} \tilde {\mathbf {L}}_{\textit{ee}} \\ \boldsymbol {\Lambda }_r^{-1}\tilde {\mathbf {L}}_{\textit{oo}} & 0 \end{pmatrix} \begin{pmatrix} \boldsymbol {r}_{\textit{o}} \\ \boldsymbol {r}_{\textit{e}} \end{pmatrix} = \lambda \begin{pmatrix} \boldsymbol {r}_{\textit{o}} \\ \boldsymbol {r}_{\textit{e}} \end{pmatrix}, \end{equation}

one sees that $\boldsymbol {\Lambda }_r^{-1} \tilde {\mathbf {L}}_{\textit{ee}} \boldsymbol {r}_{\textit{e}} = \lambda \boldsymbol {r}_{\textit{o}}$ and $\boldsymbol {\Lambda }_r^{-1} \tilde {\mathbf {L}}_{\textit{oo}} \boldsymbol {r}_{\textit{o}} = \lambda \boldsymbol {r}_{\textit{e}}$ . Then

(5.49) \begin{equation} \begin{pmatrix} 0 & \boldsymbol {\Lambda }_1^{-1} \tilde {\mathbf {L}}_{\textit{ee}} \\ \boldsymbol {\Lambda }_1^{-1}\tilde {\mathbf {L}}_{\textit{oo}} & 0 \end{pmatrix} \begin{pmatrix} -\boldsymbol {r}_{\textit{o}} \\ \boldsymbol {r}_{\textit{e}} \end{pmatrix} = \begin{pmatrix} \boldsymbol {\Lambda }_r^{-1} \tilde {\mathbf {L}}_{\textit{ee}} \boldsymbol {r}_{\textit{e}} \\ -\boldsymbol {\Lambda }_r^{-1} \tilde {\mathbf {L}}_{\textit{oo}} \boldsymbol {r}_{\textit{o}} \end{pmatrix} = \begin{pmatrix} \lambda \boldsymbol {r}_{\textit{o}} \\ -\lambda \boldsymbol {r}_{\textit{e}} \end{pmatrix} = - \lambda \begin{pmatrix} -\boldsymbol {r}_{\textit{o}} \\ \boldsymbol {r}_{\textit{e}} \end{pmatrix}, \end{equation}

which completes the proof.

By the above lemma, the eigenvalue decomposition of the matrix in (5.46) is

(5.50) \begin{equation} \begin{pmatrix} 0 & \boldsymbol {\Lambda }_r^{-1} \tilde {\mathbf {L}}_{\textit{ee}} \\ \boldsymbol {\Lambda }_r^{-1} \tilde {\mathbf {L}}_{\textit{oo}} & 0 \end{pmatrix} = \mathbf{R}_r \begin{pmatrix} \tilde {\lambda }_1 \\ & -\tilde {\lambda }_1 \\ & & \ddots \\ & & & \tilde {\lambda }_r \\ & & & & -\tilde {\lambda }_r \end{pmatrix} \mathbf{R}_r^{-1}, \end{equation}

where $\tilde {\lambda }_1 \geqslant \tilde {\lambda }_2 \geqslant \cdots \geqslant \tilde {\lambda }_s \gt 0 = \tilde {\lambda }_{s+1} = \cdots = \tilde {\lambda }_r$ for some $s \leqslant r$ , and

(5.51) \begin{equation} \mathbf{R}_r = \begin{pmatrix} \mathbf{R}_s & \mathbf{R}_0 \end{pmatrix} \end{equation}

with

(5.52) \begin{equation} \mathbf{R}_s = \begin{pmatrix} \boldsymbol {r}_{\textit {o}1} & - \boldsymbol {r}_{\textit {o}1} & \ldots & \boldsymbol {r}_{\textit {os}} & -\boldsymbol {r}_{\textit {os}} \\ \boldsymbol {r}_{\textit {e}1} & \boldsymbol {r}_{\textit {e}1} & \ldots & \boldsymbol {r}_{\textit {es}} & \boldsymbol {r}_{\textit {es}} \end{pmatrix}. \end{equation}

For simplicity, we define $\mathbf {L}_r^{\intercal } = \mathbf{R}_r^{-1}$ so that $\mathbf {L}_r$ has the form $\begin{pmatrix} \mathbf {L}_s & \mathbf {L}_0 \end{pmatrix}$ with

(5.53) \begin{equation} \mathbf {L}_s = \begin{pmatrix} \boldsymbol {l}_{\textit {o}1} & - \boldsymbol {l}_{\textit {o}1} & \ldots & \boldsymbol {l}_{\textit {os}} & -\boldsymbol {l}_{\textit {os}} \\ \boldsymbol {l}_{\textit {e}1} & \boldsymbol {l}_{\textit {e}1} & \ldots & \boldsymbol {l}_{\textit {es}} & \boldsymbol {l}_{\textit {es}} \end{pmatrix} .\end{equation}

If we assume the simulation domain is $(0,1)$ without loss of generality, the solution of such a linear system has the form

(5.54) \begin{align} \boldsymbol {v}(x) = & \begin{pmatrix} \boldsymbol {v}_{\textit{o}}(x) \\ \boldsymbol {v}_{\textit{e}}(x) \end{pmatrix} = \mathbf{R}_r \begin{pmatrix} e^{\tilde {\lambda }_1 x} \\ & e^{-\tilde {\lambda }_1 x} \\ & & \ddots \\ & & & e^{\tilde {\lambda }_r x} \\ & & & & e^{-\tilde {\lambda }_r x} \end{pmatrix} \mathbf {L}_r^{\intercal } \boldsymbol {v}(0) \notag\\ & = \sum _{j=1}^{s} \left ( e^{\tilde {\lambda }_j x} \begin{pmatrix} \boldsymbol {r}_{\textit{oj}} \\ \boldsymbol {r}_{\textit{ej}} \end{pmatrix} (\boldsymbol {l}_{\textit{oj}}^{\intercal } , \boldsymbol {l}_{\textit{ej}}^{\intercal } ) + e^{-\tilde {\lambda }_j x} \begin{pmatrix} -\boldsymbol {r}_{\textit{oj}} \\ \boldsymbol {r}_{\textit{ej}} \end{pmatrix} (-\boldsymbol {l}_{\textit{oj}}^{\intercal } , \boldsymbol {l}_{\textit{ej}}^{\intercal } ) \right ) \boldsymbol {v}(0) + \mathbf{R}_0 \mathbf {L}_0^{\intercal } \boldsymbol {v}(0) \notag\\ & = \sum _{j=1}^{s} \left ( 2\cosh (\tilde {\lambda }_j x) \begin{pmatrix} \boldsymbol {r}_{\textit{oj}} \boldsymbol {l}_{\textit{oj}}^{\intercal } \boldsymbol {v}_{\textit{o}}(0) \\[3pt] \boldsymbol {r}_{\textit{ej}} \boldsymbol {l}_{\textit{ej}}^{\intercal } \boldsymbol {v}_{\textit{e}}(0) \end{pmatrix} + 2\sinh (\tilde {\lambda }_j x) \begin{pmatrix} \boldsymbol {r}_{\textit{oj}} \boldsymbol {l}_{\textit{ej}}^{\intercal } \boldsymbol {v}_{\textit{e}}(0) \\[3pt] \boldsymbol {r}_{\textit{ej}} \boldsymbol {l}_{\textit{oj}}^{\intercal } \boldsymbol {v}_{\textit{o}}(0) \end{pmatrix} \right ) + \mathbf{R}_0 \mathbf {L}_0^{\intercal } \boldsymbol {v}(0). \end{align}

The rank $r$ for Maxwell molecules is 3 less than the rank of the non-Maxwell cases. In fact, for non-Maxwell gases, the values of $\lambda _r$ , $\lambda _{r-1}$ and $\lambda _{r-2}$ are very close to $0$ , and thus $\tilde {\lambda }_{1,2,3} \gg 1$ , which exhibits additional boundary layers. Note that when the system is reduced to the one-dimensional case as in § 5.2, only one large eigenvalue is left.

To remove the boundary layers in the non-Maxwell cases, we revise the boundary condition to eliminate the corresponding terms in (5.54). For simplicity, we only introduce the removal of the boundary layer associated with $\tilde {\lambda }_1$ . The other two layers are removed similarly. The idea is to enforce the following equalities in our boundary conditions:

(5.55) \begin{equation} \boldsymbol {l}_{\textit{o}1}^{\intercal } \boldsymbol {v}_{\textit{o}}(0) = \boldsymbol {l}_{\textit{o}1}^{\intercal } \boldsymbol {v}_{\textit{o}}(1) =0. \end{equation}

This is sufficient to get rid of the terms involving $\cosh (\tilde {\lambda }_1 x)$ and $\sinh (\tilde {\lambda }_1 x)$ in (5.54). The reason will be explained as follows.

We can follow (5.54) to write the solution using boundary conditions at $x=1$ as

(5.56) \begin{align} \boldsymbol {v}(x) &= \sum _{j=1}^{r} \left ( 2\cosh (\tilde {\lambda }_j (x-1)) \begin{pmatrix} \boldsymbol {r}_{\textit{oj}} \boldsymbol {l}_{\textit{oj}}^{\intercal } \boldsymbol {v}_{\textit{o}}(1) \\[4pt] \boldsymbol {r}_{\textit{ej}} \boldsymbol {l}_{\textit{ej}}^{\intercal } \boldsymbol {v}_{\textit{e}}(1) \end{pmatrix} + 2\sinh (\tilde {\lambda }_j (x-1)) \begin{pmatrix} \boldsymbol {r}_{\textit{oj}} \boldsymbol {l}_{\textit{ej}}^{\intercal } \boldsymbol {v}_{\textit{e}}(1) \\[4pt] \boldsymbol {r}_{\textit{ej}} \boldsymbol {l}_{\textit{oj}}^{\intercal } \boldsymbol {v}_{\textit{o}}(1) \end{pmatrix} \right ) \notag\\ & \quad+ \mathbf{R}_0 \mathbf {L}_0^{\intercal } \boldsymbol {v}(1). \end{align}

Plugging the following equalities into (5.56) gives

(5.57) \begin{align} \cosh ( \tilde {\lambda }_1(x-1)) = \cosh (\tilde {\lambda }_1 x) \cosh (-\tilde {\lambda }_1) + \sinh (\tilde {\lambda }_1 x) \sinh (-\tilde {\lambda }_1), \notag\\ \sinh ( \tilde {\lambda }_1(x-1)) = \sinh (\tilde {\lambda }_1 x) \cosh (-\tilde {\lambda }_1) + \cosh (\tilde {\lambda }_1 x) \sinh (-\tilde {\lambda }_1), \end{align}

and matching the coefficients of $\cosh (\tilde {\lambda }_1 x)$ and $\sinh (\tilde {\lambda }_1 x)$ in (5.54), one can show that

(5.58) \begin{align} \begin{pmatrix} \boldsymbol {r}_{\textit{o}1} ( \boldsymbol {l}_{\textit{o}1}^{\intercal } \boldsymbol {v}_{\textit{o}}(1) \cosh (-\tilde {\lambda }_1) + \boldsymbol {l}_{\textit{e}1}^{\intercal } \boldsymbol {v}_{\textit{e}}(1) \sinh (-\tilde {\lambda }_1) ) \\[4pt] \boldsymbol {r}_{\textit{e}1} ( \boldsymbol {l}_{\textit{e}1}^{\intercal } \boldsymbol {v}_{\textit{e}}(1) \cosh (-\tilde {\lambda }_1) + \boldsymbol {l}_{\textit{o}1}^{\intercal } \boldsymbol {v}_{\textit{o}}(1) \sinh (-\tilde {\lambda }_1) ) \end{pmatrix} = \begin{pmatrix} \boldsymbol {r}_{\textit{o}1} \boldsymbol {l}_{\textit{o}1}^{\intercal } \boldsymbol {v}_{\textit{o}}(0) \\[4pt] \boldsymbol {r}_{\textit{e}1} \boldsymbol {l}_{\textit{e}1}^{\intercal } \boldsymbol {v}_{\textit{e}}(0) \end{pmatrix}, \notag \\ \begin{pmatrix} \boldsymbol {r}_{\textit{o}1} ( \boldsymbol {l}_{\textit{o}1}^{\intercal } \boldsymbol {v}_{\textit{o}}(1) \sinh (-\tilde {\lambda }_1) + \boldsymbol {l}_{\textit{e}1}^{\intercal } \boldsymbol {v}_{\textit{e}}(1) \cosh (-\tilde {\lambda }_1) ) \\[4pt] \boldsymbol {r}_{\textit{e}1} ( \boldsymbol {l}_{\textit{e}1}^{\intercal } \boldsymbol {v}_{\textit{e}}(1) \sinh (-\tilde {\lambda }_1) + \boldsymbol {l}_{\textit{o}1}^{\intercal } \boldsymbol {v}_{\textit{o}}(1) \cosh (-\tilde {\lambda }_1) ) \end{pmatrix} = \begin{pmatrix} \boldsymbol {r}_{\textit{o}1} \boldsymbol {l}_{\textit{e}1}^{\intercal } \boldsymbol {v}_{\textit{e}}(0) \\[4pt] \boldsymbol {r}_{\textit{e}1} \boldsymbol {l}_{\textit{o}1}^{\intercal } \boldsymbol {v}_{\textit{o}}(0) \end{pmatrix} . \end{align}

We can now apply $\boldsymbol {l}_{\textit{o}1}^{\intercal } \boldsymbol {v}_{\textit{o}}(0) = \boldsymbol {l}_{\textit{o}1}^{\intercal } \boldsymbol {v}_{\textit{o}}(1) =0$ to the above equation, and it is not difficult to derive $\boldsymbol {l}_{\textit{e}1}^{\intercal } \boldsymbol {v}_{\textit{e}}(0) = \boldsymbol {l}_{\textit{e}1}^{\intercal } \boldsymbol {v}_{\textit{e}}(1) = 0$ . Therefore, all the terms involving $\cosh (\tilde {\lambda }_1 x)$ and $\sinh (\tilde {\lambda }_1 x)$ in (5.54) or (5.56) disappear so that the corresponding boundary layers are not present in the solution.

The above derivation shows that (5.55) is a sufficient condition to suppress boundary layers. Moreover, it is straightforward to generalise such a condition to multidimensional cases by imposing $\boldsymbol {l}_{\textit{o}1}^{\intercal } \boldsymbol {v}_{\textit{o}}=0$ for $\boldsymbol {v}_{\textit{o}}$ at all boundary points. Recall that the Onsager boundary conditions take a form of (1.6), which can be reformulated as

(5.59) \begin{equation} \boldsymbol {v}_{\textit{o}} = \tilde {\mathbf {Q}} \boldsymbol {\Lambda }(\tilde {\boldsymbol {g}}_{\textit{ext}} - \boldsymbol {v}_{\textit{e}}). \end{equation}

To impose $\boldsymbol {l}_{\textit{o}1}^{\intercal } \boldsymbol {v}_{\textit{o}}=0$ , we choose to modify the symmetric matrix $\tilde {\mathbf {Q}}$ such that

(5.60) \begin{equation} \boldsymbol {l}_{\textit{o}1}^{\intercal } \tilde {\mathbf {Q}} = 0. \end{equation}

However, (5.60) is an underdetermined linear system. To make it uniquely solvable, we add the following constraints.

1. (i) The revised $\tilde {\mathbf {Q}}$ must be symmetric.

2. (ii) An element of $\tilde {\mathbf {Q}}$ should remain zero if it is zero for all non-Maxwell gas molecules.

3. (iii) An element $Q_{ij} = {2\chi }/{2-\chi } {\langle \phi _{\textit{odd}}^i, \mathcal {C} \xi _n^{-1} \phi _{\textit{odd}}^j \rangle }/{\langle \phi _{\textit{odd}}^j, \phi _{\textit{odd}}^j \rangle }$ (refer to (5.16)) of $\tilde {\mathbf {Q}}$ should remain unchanged if both $\phi _{\textit{odd}}^i$ and $\phi _{\textit{odd}}^j$ belong to $\mathbb {V}_{\textit{odd}} \cap (\mathbb {V}_{\textit{modified}})^{\perp }$ , where

   (5.61) \begin{equation} \mathbb {V}_{\textit{modified}} = \operatorname {span}\left \{ \phi _{n\tau _1}^{2}, \phi _{n\tau _2}^{2}, \phi _n^{2} + \frac {c^{3,1}_1}{c^{2,1}_1} \phi _n^{3}\right \} \end{equation}
   is a subspace of the $\mathbb {V}_{\textit{odd}}$ in (5.12).

The rationale for proposing constraint (iii) is that $\mathbb {V}_{\textit{modified}} = \emptyset$ in the case of Maxwell gas, which does not exhibit any non-physical boundary layers. The additional basis functions in $\mathbb {V}_{\textit{modified}}$ introduce these layers in the non-Maxwell cases. Therefore, we only modify those elements deduced from $\mathbb {V}_{\textit{modified}}$ .

After imposing these conditions, the matrix $\tilde {\mathbf {Q}}$ is uniquely determined. Currently, we do not have theoretical guarantee of the semidefiniteness of the revised $\tilde {\mathbf {Q}}$ , but this turns out to be true in all cases done in our tests. The symmetry and semidefiniteness of $\tilde {\mathbf {Q}}$ indicates that the structure (1.6) is well maintained after the modification, and thus, does not ruin the $L^2$ stability or the second law of thermodynamics.

To demonstrate the effect of this modification, we reconsider the example described in § 5.2, with the boundary conditions changed to the new ones introduced in this section by replacing $\hat {m}_{ij}$ with ${m}_{ij}$ in (5.27)–(5.30). The results are plotted in figure 2. A comparison with figure 1 clearly demonstrates that the non-physical boundary layers have been removed.

Figure 2. Plots of $\bar {q}$ and $\theta$ for the one-dimensional problem (5.22)–(5.30) with modified coefficient $\hat {m}_{ij} \to {m}_{ij}$ . Results are shown for (a) $\theta ^W=0$ and (b) $\theta ^{W}=0.2$ .

6.1. Steady-state examples

To compare with the example in §§ 5.2 and 5.3, we consider a steady-state problem here. In fact, (5.22)–(5.26) can be regarded as a portion of the one-dimensional steady-state channel flow considered here by using the following replacements: $x \to x_2$ , $v \to v_2$ , $\bar {q} \to \bar {q}_2$ and $\bar {\sigma } \to \bar {\sigma }_{22}$ . Boundary conditions (5.27)–(5.30) also share the same form of boundary conditions considered here by using the above replacements. We adopt the coefficients $\theta ^W_l = 0$ and $\theta ^W_r = 0.2$ , consistent with the example presented in § 5.2.

To compare our solutions with the DSMC results, we show the results by using the Knudsen number $\overline {\textit{Kn}}$ defined by the ratio of the mean free path to the characteristic length as in Bird (Reference Bird1994). The relationship between $\overline {\textit{Kn}}$ and the parameter $\textit{Kn}$ in this paper is

(6.2) \begin{equation} {\textit{Kn}} = \sqrt {\frac {\pi }{2}} \frac {15 \overline {\textit{Kn}}}{(5 - 2 \omega )(7-2 \omega )}, \qquad \omega = \frac {1}{2} + \frac {2}{\eta - 1}. \end{equation}

Moreover, we plot the actual heat flux $q_2$ instead of $\bar {q}_2$ (see (3.8)) in this example, since $q_2$ is obtained from the DSMC results. To observe the effect of the Knudsen number, we choose two different $\overline {\textit{Kn}}$ and plot the results for both original and fixed boundary conditions in figures 3 and 4. Our results are depicted as the blue, yellow and green solid lines corresponding to $\eta = 5, 10$ and $\infty$ . The relevant solutions obtained by DSMC are presented by dotted lines of the same colours.

Figure 3. Results of the steady-state example when $\overline {\textit{Kn}}=0.1$ . (a, b) Coefficients $\hat {m}_{ij}$ are used, (b, c) modified coefficients $m_{ij}$ are used. The DSMC solutions are given by dotted lines of the same colours.

Figure 4. Results of steady-state example when $\overline {\textit{Kn}}=0.05$ . (a, b) Coefficients $\hat {m}_{ij}$ are used, (b, c) modified coefficients $m_{ij}$ are used. The DSMC solutions are given by dotted lines of the same colours.

The solids lines of the R13 equations generally agree with the dotted lines of the DSMC results, and this agreement is more pronounced for smaller $\textit{Kn}$ , which validates our equations. Moreover, as the DSMC results display no boundary layer, the effectiveness of our modified Onsager boundary conditions is again validated by the successful removal of the undesired boundary layers. Furthermore, compared with the results in figure 1, the width of the undesired boundary layers decreases as the Knudsen number decreases, which is consistent with our analysis of the solution in (5.40). Also, it can be observed that the temperature jump gets more significant as the Knudsen number increases.

6.2. Time-dependent examples

After validating the steady-state equations, it would be more compelling to verify the derived model in the time-dependent case, as this is the primary focus of this work. However, it would be challenging to obtain the analytic time-dependent solution, therefore, we aim to solve the time-dependent equations using the finite element method.

The numerical set-up is described as follows. We use a uniform mesh with a mesh size of $1/1000$ and apply the finite element method with piecewise linear functions. The Crank–Nicolson method is used for temporal discretisation with a uniform time step of $1/4000$ . The initial value for $\theta$ is set to $ ({\theta ^W_l+\theta ^W_r})/{2}$ , while the initial values for all other moments are set to $0$ . Boundary coefficients are taken from Cai et al. (Reference Cai, Torrilhon and Yang2024b ), which will be specified in the following two cases.

6.2.1. Couette flow

In the planar Couette flow, the two plates move in opposite directions and maintain the same temperature. We choose $\theta ^W_l = \theta ^W_r = 0$ and $ v^W_l = -v^W_r = -0.2$ . Results are depicted in figure 5, where non-physical boundary layers can still be observed in the solution of time-dependent equations using the previous Onsager boundary conditions, although the width of these layers is very thin in this example. These thin layers are successfully removed at various time steps using our modified boundary conditions.

Figure 5. Results of the Couette flow. (a—d) Onsager boundary conditions are used. (e—h) Our modified boundary conditions are used.

6.2.2. Fourier flow

In the planar Fourier flow, the two plates are stationary, implying that $v^W_l = v^W_r = 0$ , and the gas dynamics between the plates are driven by the temperature difference. We set $\theta ^W_l = 0$ and $ \theta ^W_r = 0.2$ . The results in figure 6 demonstrate that the previous Onsager boundary conditions produce non-physical boundary layers, which can be eliminated using our modified boundary conditions. Additionally, these layers are more pronounced for longer computational times, suggesting our modified boundary conditions would be more important in long-time simulations of the proposed R13 equations.

Figure 6. Results of the Fourier flow. (a—d) Previous Onsager boundary conditions are used. (e—h) Our modified boundary conditions are used.