1. Introduction
The improvement of manufacturing technology and computational capacity has led to increasing attention to the importance of rarefied gas flows in various practical applications such as high-speed hypersonic flows (Ivanov & Gimelshein Reference Ivanov and Gimelshein1998), miniaturized microchannel flows (Agrawal, Kushwaha & Jadhav Reference Agrawal, Kushwaha and Jadhav2020; Akhlaghi, Roohi & Stefanov Reference Akhlaghi, Roohi and Stefanov2023) and rarefied flow conditions (Akintunde & Petculescu Reference Akintunde and Petculescu2014). In such flow situations, the Knudsen number (
$Kn= \lambda /H$), which is the ratio of the mean free path (
$\lambda$) to the characteristic length scale (
$H$), plays a crucial role in governing the flow physics. For the case of a microchannel, the flow generally spans the slip (
$10^{-3}< Kn< 10^{-1}$) or transition (
$10^{-1}< Kn< 10$) regimes for which the well-known Navier–Stokes (N-S) equations based on the continuum formulation fail to predict the primary variables correctly (García-Colín, Velasco & Uribe Reference García-Colín, Velasco and Uribe2008; Agrawal et al. Reference Agrawal, Kushwaha and Jadhav2020). Further, several analytical and numerical studies have demonstrated that the N-S equations do not capture, even qualitatively, non-equilibrium phenomena occurring in the plane Couette and Poiseuille flows and lid-driven cavity problems in the slip flow regime (Tij & Santos Reference Tij and Santos1994; Mansour, Baras & Garcia Reference Mansour, Baras and Garcia1997; Zheng, Garcia & Alder Reference Zheng, Garcia and Alder2002; Benzi, Gu & Emerson Reference Benzi, Gu and Emerson2010, Reference Benzi, Gu and Emerson2013). However, the N-S equations can capture the flow physics in the early-slip flow regime by employing a 
$Kn$-dependent slip and temperature jump boundary conditions (Agrawal et al. Reference Agrawal, Kushwaha and Jadhav2020); such modifications still yield inaccurate predictions of the pressure drop, shear stress, heat flux and mass flow rate in the late-slip and transition flow regimes. Thus, improved flow modelling capabilities are needed to accurately describe flows, particularly in the transition flow regime.
Among the available analytical approaches to describe late-slip/transition flow regimes, the kinetic theory-based Boltzmann equation describes the non-equilibrium dynamics of dilute gaseous flows most effectively. However, due to the complicated nonlinear integro-differential nature of the collision integral, it is exceptionally challenging, if at all possible, to solve the Boltzmann equation for even simple flow configurations. The Boltzmann equation is, more often than not, made tractable by incorporating kinetic models such as the Bhatnagar–Gross–Krook (BGK) model (Bhatnagar, Gross & Krook Reference Bhatnagar, Gross and Krook1954). Indeed, the simplified Boltzmann–BGK equation has been used to numerically solve the simple isothermal pressure-driven Poiseuille flow problem and has demonstrated the ability to quantitatively capture the experimentally observed Knudsen minima phenomenon (Cercignani & Daneri Reference Cercignani and Daneri1963). Similarly, Xu (Reference Xu2003) solved the force-driven plane Poiseuille flow problem by employing the BGK–Burnett and the 
$\mathcal {O}$(
$Kn^3$) BGK–super-Burnett schemes in the Boltzmann equations and demonstrated better results than those obtained using the N-S equations. An alternate way of solving the Boltzmann–BGK equation directly is through the lattice Boltzmann method (LBM) route; recent attempts to employ LBM schemes to solve high-
$Kn$ flows have been detailed in Jonnalagadda, Sharma & Agrawal (Reference Jonnalagadda, Sharma and Agrawal2023). More recently, the ability of higher-order transport equations to address canonical wall-bounded flows in the slip and transition regimes has been investigated. It should be noted that several higher-order transport equations with different degrees of accuracy and complexity have been proposed based on the single-particle distribution function obtained from either of the following three approaches: the Chapman–Enskog expansion (Chapman & Cowling Reference Chapman and Cowling1970; Cercignani Reference Cercignani1975), moment-based methods (Grad Reference Grad1949) or equations incorporating principles of non-equilibrium thermodynamics (Singh & Agrawal Reference Singh and Agrawal2016; Singh, Jadhav & Agrawal Reference Singh, Jadhav and Agrawal2017; Jadhav, Yadav & Agrawal Reference Jadhav, Yadav and Agrawal2023).
The first category of higher-order transport equations, namely that of the Burnett-type equations, originated from the original Burnett equations (Burnett Reference Burnett1936). However, due to their complicated and cumbersome nature, the original Burnett equations did not gain widespread adoption; indeed, numerical studies using the original Burnett equations were only recently conducted for the Couette flow problem (Lockerby & Reese Reference Lockerby and Reese2003; Singh, Gavasane & Agrawal Reference Singh, Gavasane and Agrawal2014). In order to make the original Burnett equations more tractable, Chapman & Cowling (Reference Chapman and Cowling1970) presented an alternate form, namely the conventional Burnett equations, wherein the material derivatives appearing in the second-order distribution function were replaced by the spatial gradients described by the Euler equations. These modifications resulted in significantly simplified representations of the stress tensor and heat flux vector, and consequently, a wider use of the conventional Burnett equations followed. Note that both the original and conventional Burnett equations are 
$\mathcal {O}$(
$Kn^2$) approximations of the Boltzmann equation. Uribe & Garcia (Reference Uribe and Garcia1999) numerically solved the conventional Burnett equation for the Poiseuille flow problem in the transition regime at 
$Kn =0.1$ and obtained good agreement for most conserved and non-conserved quantities with direct simulation Monte Carlo (DSMC) results. Fang (Reference Fang2003) numerically investigated the plane Poiseuille flow with the conventional Burnett equations at 
$Kn = 0.02$ with two slip models. Later, Bao, Lin & Shi (Reference Bao, Lin and Shi2007) employed the conventional Burnett equations and reported that, in contrast to the N-S equations, convergent results could be obtained for the plane Couette flow problem at any Knudsen number. Similarly, analytical solutions for the pressure fields were presented for the pressure-driven Poiseuille flow problem using the conventional Burnett equations by employing a perturbation analysis (Rath, Singh & Agrawal Reference Rath, Singh and Agrawal2018) and exactly solving the third-order partial differential cross-stream momentum equation (Rath, Yadav & Agrawal Reference Rath, Yadav and Agrawal2021). Several wall-bounded flow studies were also conducted using 
$\mathcal {O}$(
$Kn^2$) BGK–Burnett equations. Aoki, Takata & Nakanishi (Reference Aoki, Takata and Nakanishi2002) employed an asymptotic analysis to obtain a numerically solvable form of the BGK–Burnett equations for the force-driven plane Poiseuille unidirectional flow problem in the slip regime. For microchannel flows, Xu & Li (Reference Xu and Li2004) demonstrated that the pressure, velocity and mass flow rate obtained using the BGK–Burnett equations quantitatively agree with experimental and DSMC results. Further, using 
$\mathcal {O}$(
$Kn^3$) augmented Burnett equations, which introduced certain super-Burnett terms into the original Burnett equations to improve numerical stability, Agarwal, Yun & Balakrishnan (Reference Agarwal, Yun and Balakrishnan2001) solved the plane Poiseuille flow problem numerically by employing Beskok's and Langmuir's boundary conditions for an inlet Knudsen number of 
$Kn=0.088$. Bao & Lin (Reference Bao and Lin2008) used the augmented Burnett equations to numerically solve the plane Poiseuille flow problem by employing a relaxation method on the boundary values, in line with that previously used by Lockerby & Reese (Reference Lockerby and Reese2003) to present comparisons of the velocity and pressure fields against experimental and DSMC data. More recently, the linearized form of the 
$\mathcal {O}$(
$Kn^3$) super-Burnett equations (Shavaliyev Reference Shavaliyev1993) was solved analytically for the Couette flow problem in the transition regime (Singh et al. Reference Singh, Gavasane and Agrawal2014). It should also be noted that, despite widespread use and success, the above-mentioned Burnett variants are associated with several limitations (Bobylev Reference Bobylev1982; Shavaliyev Reference Shavaliyev1993; Uribe & Garcia Reference Uribe and Garcia1999; García-Colín et al. Reference García-Colín, Velasco and Uribe2008; Dadzie Reference Dadzie2013) which have prompted the development of additional Burnett-like variants. These variants include the thermo-mechanically consistent Burnett equations (Dadzie Reference Dadzie2013) and the simplified Burnett equations (Zhao, Chen & Agarwal Reference Zhao, Chen and Agarwal2014).
More recently, a separate class of Burnett-like equations, namely the OBurnett equations, have been derived from the perspective of linear irreversible non-equilibrium thermodynamics (Singh et al. Reference Singh, Jadhav and Agrawal2017). It is evident that the OBurnett equations and the Burnett equations diverge in their foundational principles and origins. Consequently, the OBurnett equations are not plagued by the specific limitations inherent to the Burnett equations (Bobylev Reference Bobylev1982; Comeaux et al. Reference Comeaux, Chapman, Ma and Cormack1995; Uribe, Velasco & Garcia-Colin Reference Uribe, Velasco and Garcia-Colin2000; García-Colín et al. Reference García-Colín, Velasco and Uribe2008; De Groot & Mazur Reference De Groot and Mazur2013; Agrawal et al. Reference Agrawal, Kushwaha and Jadhav2020), as confirmed in Singh et al. (Reference Singh, Jadhav and Agrawal2017), Agrawal et al. (Reference Agrawal, Kushwaha and Jadhav2020), Jadhav & Agrawal (Reference Jadhav and Agrawal2020b, Reference Jadhav and Agrawal2021) and Jadhav et al. (Reference Jadhav, Yadav and Agrawal2023). Readers are directed to the above-mentioned studies for a more detailed exploration of the differences and advantages of the OBurnett equations relative to the Burnett equations. This set of equations has been validated for several canonical problems involving non-equilibrium flow conditions (Jadhav, Singh & Agrawal Reference Jadhav, Singh and Agrawal2017; Singh et al. Reference Singh, Jadhav and Agrawal2017; Jadhav & Agrawal Reference Jadhav and Agrawal2020b, Reference Jadhav and Agrawal2021; Yadav & Agrawal Reference Yadav and Agrawal2021); a consolidated account of the different validation cases can be found in Jadhav et al. (Reference Jadhav, Yadav and Agrawal2023). Specifically, the validation of the OBurnett equations for the force-driven Poiseuille flow problem (Jadhav et al. Reference Jadhav, Singh and Agrawal2017) is particularly relevant in the context of wall-bounded flows in the late-slip/transition regime. However, as is the case with a large portion of the existing literature employing Burnett-type higher-order transport equations for this problem, Jadhav et al. (Reference Jadhav, Singh and Agrawal2017) employ numerical methods to solve the resulting system of equations.
In this work, we employ an alternative approach to express the second-order approximation of the distribution function in terms of an iterative refinement expression of the Chapman–Enskog multiscale expansion technique that is further modified to ensure compliance with the additive invariance property. With this approach, we derive a more generalized form of the distribution function which encapsulates the formulations previously obtained by Singh et al. (Reference Singh, Jadhav and Agrawal2017) and Yadav, Jonnalagadda & Agrawal (Reference Yadav, Jonnalagadda and Agrawal2023). Using these representations of the distribution function, constitutive relationships for the stress tensor and heat flux vector are obtained at the Burnett and super-Burnett levels. It is noteworthy that this derivation procedure yields constitutive relationships, which include several additional linear terms, which are necessary to accurately account for near-wall non-equilibrium effects (Taheri, Torrilhon & Struchtrup Reference Taheri, Torrilhon and Struchtrup2009). Further, due to the additional linear terms, the linearized form of the proposed equations differs from the N-S equations and thus could provide better representations of wall-bounded flows in the transition regime. We demonstrate that the proposed equations are unconditionally linearly stable against any spatial disturbances for two-dimensional (2-D) plane wave flows. Apart from linear stability analysis, analytical solutions for rarefied long microchannel flow can be extremely instructive. Hence, one of the primary objectives of the present work is to present the complete analytical solution of the 2-D isothermal, pressure-driven plane Poiseuille flow problem using newly derived extended-OBurnett (EOBurnett) and a third-order accurate super-OBurnett (SOBurnett) equations. For this purpose, we obtain, for the first time, the mathematical expression for the pressure and velocity by analytically solving the axial and normal momentum equations using the perturbation method. This analytical solution is obtained under the assumption of low Mach and Reynolds numbers relevant to long microchannels and are shown to demonstrate better agreement with DSMC data as compared with those obtained from the N-S equations. In summary, the contributions of this work are threefold: (i) we derive new sets of higher-order continuum transport equations, namely the 
$\mathcal {O}(Kn^2)$ EOBurnett equations and 
$\mathcal {O}(Kn^3)$ SOBurnett equations, (ii) perform a 2-D linear stability analysis of the proposed equations and (iii) derive the analytical solution for the 2-D plane Poiseuille flow problem.
The paper is organized as follows: § 2 presents the derivation of the distribution function, which satisfies the additive invariance property and is consistent with Onsager's reciprocity principle. In § 3, the distribution function derived in § 2 has been used to evaluate the constitutive relationships. We further obtain a more generalized form of the second- and third-order accurate transport equations in this section. Section 4 presents a stability analysis, and we further show that the proposed second- and third-order equations are unconditionally stable to small disturbances in a 2-D flow. In § 5, we present the complete analytical solution for the 2-D pressure-driven plane Poiseuille flow problem. The solution has been derived using the perturbation method. Further, the first-, second- and third-order solutions of pressure and velocity are compared with DSMC results for flow in the transition regime under low Reynolds and Mach number conditions. Section 6 compares the proposed equation with the existing variants of the Burnett equations and highlights the important findings concerning the significance, novelty and benefit of the present work. Section 7 summarizes the important contributions made in this work.
2. Derivation of complete second- and third-order Burnett-level representations of the distribution function
The OBurnett equations (Singh & Agrawal Reference Singh and Agrawal2016), which are obtained as the 
$\mathcal {O}$(
$Kn^2$) approximation to the Boltzmann equation, require the 
$\mathcal {O}$(
$Kn^2$) representation of the single-particle distribution function, 
$f(\boldsymbol {x}, \boldsymbol {c}, t)$, where 
$\boldsymbol {x}$, 
$\boldsymbol {c}$, 
$t$ are the position vector, molecular velocity and time, respectively. This 
$\mathcal {O}$(
$Kn^2$) distribution function can be expressed as
\begin{equation} f_{2} = f_0 + \bar{f}_{1} + \bar{f}_{2}, \end{equation}
where 
$f_0$ is the Maxwell–Boltzmann equilibrium distribution function
\begin{equation} f_0 = \frac{\rho}{m}\left(\frac{\beta}{\rm \pi}\right)^{{3}/{2}} \exp{[-\beta (\lvert \boldsymbol{c}-\boldsymbol{u} \rvert)^2]}, \end{equation} with 
$\beta = 1/(2 R T)$ and where u is the macroscopic velocity, 
$\rho $ is mass density, R is specific gas constant and T is absolute temperature. The remaining symbols have their usual meaning. The first (N-S) order correction to 
$f_0$ is obtained from the following relationship (Mahendra & Singh Reference Mahendra and Singh2013):
\begin{equation} \bar{f}_{1} ={-}\sum_j\varUpsilon_j \odot X_j = t_{r(j)} \left[\frac{\partial f_0}{\partial t}+ \frac{\partial}{\partial x_j} (c_j f_0)\right]_{X_i=0,\quad \forall i \neq j},\enspace j\in\{\tau, q\}, \end{equation}
where the quantities 
$X_j$ and 
$\varUpsilon _j$ correspond to the destabilizing thermodynamic forces and fluxes, respectively, while the subscripts 
$\tau$ and 
$q$ correspond to irreversible viscous and thermal thermodynamic processes, respectively (McCourt et al. Reference McCourt, Beenakker, Köhler and Kušcer1991; Mahendra & Singh Reference Mahendra and Singh2013; Singh & Agrawal Reference Singh and Agrawal2016; Yadav et al. Reference Yadav, Jonnalagadda and Agrawal2023). The symbol 
$\odot$ represents a complete tensorial contraction of tensors of the same order.
For monatomic gasses, simplifying (2.3) yields,
\begin{equation} {\varUpsilon}_i ={-} f_0 t_{r(i)}\bar{\varUpsilon}_i, \end{equation}where
$$\begin{gather} {\bar{\varUpsilon}}_\tau ={-}\left[ \boldsymbol{C}\otimes \boldsymbol{C} -\frac{1}{3} \lvert \boldsymbol{C} \rvert ^2 \boldsymbol{\delta}\right], \end{gather}$$
$$\begin{gather}{\bar{\varUpsilon}}_q ={-} \left( \frac{5}{2 \beta} - \lvert \boldsymbol{C} \rvert^2 \right) \boldsymbol{C}, \end{gather}$$and
$$\begin{gather} {X}_\tau = \beta [\boldsymbol{\nabla} \otimes \boldsymbol{u} + (\boldsymbol{\nabla} \otimes \boldsymbol{u})^{\rm T}], \end{gather}$$
$$\begin{gather}{X}_q = \boldsymbol{\nabla} \beta. \end{gather}$$
The quantity 
$\boldsymbol {\delta }$ and symbol 
$\otimes$ represent the Kronecker delta and the outer product, respectively. The peculiar velocity is defined as 
$\boldsymbol {C} = (\boldsymbol {c}-\boldsymbol {u})$. Here, we bring to the reader's attention the use of different relaxation times for momentum (
$t_{r(\tau )} = \mu /p$) and energy transport (
$t_{r(q)} = \kappa (\gamma -1)/(R \gamma p)= t_{r(\tau )}/Pr$) which not only serves to separate the viscous and thermal time scales but also intrinsically ensures the correct value of Prandtl number 
$Pr$ for gases. The dynamic viscosity and thermal conductivity are treated as temperature-dependent functions and have the functional forms 
$\mu (T) = \mu _0(T/T_0)^\varphi$ and 
$\kappa (T) = \kappa _0(T/T_0)^\varphi$, respectively. Here, 
$\mu _0$ and 
$\kappa _0$ are the dynamic viscosity and thermal conductivity at a reference temperature 
$T_0$, and 
$\varphi$ is the interaction potential between two gaseous molecules. The quantities 
$\gamma$ and 
$p$ appearing in 
$t_{r(q)}$ are the adiabatic index and thermodynamic pressure, respectively.
For the second-order correction to 
$f_0$, 
$\bar {f}_{2}$ a detailed account of the derivation procedure was recently presented along with the observation that the additive invariance constraint imposed on 
$\bar {f}_{2}$ can yield several non-unique forms for 
$f_2$ (Yadav et al. Reference Yadav, Jonnalagadda and Agrawal2023). One realization of the 
$\mathcal {O}$(
$Kn^2$) representation was reported as (Yadav et al. Reference Yadav, Jonnalagadda and Agrawal2023)
\begin{equation} \bar{f}_{2} = [({\varUpsilon'_{\tau\tau}} \odot {X}_\tau) \odot X_\tau] + [({\varUpsilon'_{qq}} \odot {X}_q) \odot X_q], \end{equation}where
\begin{align} ({\varUpsilon'_{\tau\tau}} \odot
{X}_\tau) &= ({\varUpsilon}_{\tau\tau} \odot {X}_\tau) +
t_{r(\tau)}^2 f_0 \nonumber\\ \quad
&\quad \times\left\{\begin{gathered} C_l \left[C_i\frac{\partial
u_j}{\partial x_l} + \left(C_j\frac{\partial u_i}{\partial
x_l}\right)^{\rm T}\right] - \left[\frac{1}{3 \beta} ({C
\otimes C}) : {X}_\tau\right]\delta_{ij}\\
-\bar{\varUpsilon}_\tau
(\bar{\varUpsilon}_\tau: {X}_\tau) -\left(1
+\frac{1}{\varphi}\right)\bar{\varUpsilon}_\tau
\left(\frac{{\varphi}}{\beta} C_l \frac{\partial \beta }
{\partial x_l}\right)\end{gathered}\right\},
\end{align}
and
\begin{equation}
({\varUpsilon'_{qq}}\odot{X}_q) =
({\varUpsilon}_{qq}\odot{X}_q)+ t_{r(q)}^2 f_0
\left\{-\left(1+\frac{2}{\varphi} \right)
\bar{\varUpsilon}_q \frac{{\varphi}}{\beta} C_l
\frac{\partial \beta } {\partial x_l}\right\}.
\end{equation}
Here, the terms 
${\varUpsilon }_{\tau \tau } \odot {X}_\tau$ and 
${\varUpsilon }_{qq}\odot {X}_q$ are, respectively, defined as
\begin{align}
{\varUpsilon}_{\tau\tau}\odot{X}_{\tau} &=
t_{r({\tau})}^2 f_0 \nonumber\\
&\quad \times\left\{\begin{gathered} \underbrace{ - C_l \left[C_i
\frac{\partial u_j}{\partial x_l} + \left(C_j\frac{\partial
u_i}{\partial x_l}\right)^{\rm T}\right]}_{\omega_1} \quad
+\underbrace{\frac{1}{2\beta}\left[C_i \frac{\partial
g}{\partial x_j} + \left(C_j\frac{\partial g}{\partial
x_i}\right)^{\rm T}\right]}_{\omega_2}\\
-\left[\underbrace{\frac{1}{3 \beta} C_k \frac{\partial g
}{\partial x_k}}_{\omega_3} - \underbrace{\frac{1}{3 \beta}
(\boldsymbol{C \otimes C}) : {X}_\tau}_{\omega_4}\right]
\delta_{ij} + \bar{\varUpsilon}_\tau
[\underbrace{\bar{\varUpsilon}_\tau: {X}_\tau}_{\omega_{5}}
+ \underbrace{\bar{\varUpsilon}_q \cdot
{X}_q}_{\omega_{6}}]\\
+\bar{\varUpsilon}_\tau
\left[\underbrace{\frac{2\varphi-5}{3} \frac{\partial u_l
}{\partial x_l}}_{\omega_{7}} +
\underbrace{\frac{{\varphi}}{\beta} C_l \frac{\partial
\beta } {\partial x_l}}_{\omega_{8}} + \underbrace{ C_l
\frac{\partial g }{\partial x_l}}_{\omega_{9}}\right]
\end{gathered}\right\},\end{align}
and
\begin{align} {\varUpsilon}_{qq}\odot
{X}_{q} &= t_{r(q)}^2 f_o \nonumber\\ &\quad
\times\left\{\begin{gathered} \underbrace{
\bar{\varUpsilon}_q[\bar{\varUpsilon}_\tau: {X}_\tau +
\bar{\varUpsilon}_q \cdot {X}_q]}_{\xi_1} - C_i
\left[\underbrace{\frac{1}{\beta} C_l \frac{\partial
g}{\partial x_l}}_{\xi_2} - \underbrace{\frac{1}{\beta}
(\boldsymbol{C}\otimes\boldsymbol{C}) :{X}_\tau
}_{\xi_3}\right] \\
-C_i \left[\underbrace{
\frac{5}{3\beta} \frac{\partial u_k}{\partial x_k}
}_{\xi_4} + \underbrace{\frac{5}{2\beta^2} C_l
\frac{\partial \beta}{\partial x_l} }_{\xi_5}\right] +
\left(\frac{5}{2 \beta} -\lvert C \rvert^2 \right)
\left[\underbrace{\frac{1}{2\beta} \frac{\partial
g}{\partial x_i} }_{\xi_6} - \underbrace{C_l \frac{\partial
u_i}{\partial x_l} }_{\xi_7}\right]\\
+ \bar{\varUpsilon}_q
\left[\underbrace{\frac{2\varphi -5}{3} \frac{\partial u_l
}{\partial x_l} }_{\xi_{8}} +
\underbrace{\frac{{\varphi}}{\beta} C_l \frac{\partial
\beta } {\partial x_l}}_{\xi_{9}} + \underbrace{C_l
\frac{\partial g }{\partial x_l} }_{\xi_{10}}\right]
\end{gathered}\right\}.
\end{align}
The advantage of the form presented in (2.5)–(2.7) is that, in addition to satisfying the compatibility conditions, this realization also adheres to Onsager's symmetry principle. Furthermore, building upon (2.5)–(2.7), this work extends and presents linearly stable macroscopic higher-order continuum transport equations, a property that the conventional Burnett and super-Burnett equations do not admit.
Similarly, the 
$\mathcal {O}$(
$Kn^3$) additive invariant representation of the viscous and thermal corrections appearing in 
$\bar {f}_{2}$ was obtained as (Yadav et al. Reference Yadav, Jonnalagadda and Agrawal2023)
\begin{equation} \bar{f}^{'}_{2} = [({\varUpsilon''_{\tau\tau}} \odot {X}_\tau) \odot X_\tau] + [({\varUpsilon''_{qq}} \odot {X}_q) \odot X_q], \end{equation}where
\begin{align} & ({\varUpsilon''_{\tau\tau}}\odot {X}_{\tau}) = ({\varUpsilon'_{\tau\tau}} \odot {X}_\tau) + t_{r({\tau})}^2 f_0 \left\{- \bar{\varUpsilon}_\tau \left[\frac{4 \beta^2}{3 \rho } \varOmega\left(\frac{3}{2\beta} -C_l^{2}\right) -\frac{2 \beta C_l}{\rho} \frac{\partial \sigma^{NS}_{lk}}{\partial x_k}\right] \right.\nonumber\\ &\quad \left.+ \,\bar{\varUpsilon}_\tau\left[(\varphi -1) \frac{4 \beta}{3\rho} \varOmega\right] +\left[\frac{1}{\rho} \frac{\partial \sigma^{NS}_{ik}}{\partial x_k}C_j + \frac{1}{\rho} C_i \frac{\partial \sigma^{NS}_{jk}}{\partial x_k} - \frac{2}{3} C_k \frac{1}{\rho} \frac{\partial \sigma^{NS}_{kl}}{\partial x_l} \delta_{ij}\right]\right\}, \end{align}and
\begin{align} & ({\varUpsilon''_{qq}}\odot {X}_{q}) = ({\varUpsilon'_{qq}} \odot {X}_q) + t_{r({q})}^2 f_0 \left\{-\bar{\varUpsilon}_q \left[\frac{4 \beta^2}{3 \rho } \varOmega \left(\frac{3}{2\beta} -C_l^{2}\right) -\frac{2 \beta C_l}{\rho} \frac{\partial \sigma^{NS}_{lk}}{\partial x_k}\right] \right.\nonumber\\ &\quad \left.+\,\bar{\varUpsilon}_q \left[(\varphi -1) \frac{4 \beta}{3\rho} \varOmega\right] + \left[C_i \left( -\frac{10}{3 \rho} \varOmega - \frac{2}{\rho} C_k \frac{\partial \sigma^{NS}_{kj}}{\partial x_j} \right) + \left(\frac{5}{2 \beta} -\lvert C \rvert^2 \right) \left(\frac{1}{\rho} \frac{\partial \sigma^{NS}_{ij}}{\partial x_j}\right)\right]\right\}, \end{align}
where 
$\varOmega$, 
$\sigma _{ij}^{{NS}}$ and 
$q_i^{{NS}}$ are, respectively, defined as
\begin{equation} \varOmega = \left[\frac{\partial q^{NS}_l}{\partial x_l} + \sigma^{NS}_{lk} \frac{\partial u_l}{\partial x_k}\right] , \end{equation}
$$\begin{gather} \sigma_{ij}^{{NS}}= \sigma_{ij}^{\mathcal{O}(Kn^1)} =\langle C_{\langle i}C_{j\rangle}, {\bar{f}_1} \rangle ={-} 2 \mu \left(\frac{1}{2}\left[\frac {\partial u_i}{\partial x_j}+ \frac{\partial u_j}{\partial x_i}\right] -\frac{1}{3} \frac {\partial u_l}{\partial x_l} \delta_{ij}\right) = {-2 \mu \frac {\partial u_{\langle i}}{\partial x_{j\rangle}}}, \end{gather}$$
$$\begin{gather}q_i^{{NS}} = q_{i}^{\mathcal{O}(Kn^1)} = \frac{1}{2} \langle \lvert\boldsymbol{C}\rvert^2 C_i, {\bar{f}_1} \rangle ={-}\kappa \frac{\partial T}{\partial x_i}. \end{gather}$$ The terms 
$({\varUpsilon }_{kj}\odot {X}_{k})$ appearing in equations (2.7a) and (2.7b) are obtained using the relation (Mahendra & Singh Reference Mahendra and Singh2013)
\begin{equation} {\varUpsilon}_{kj}\odot {X}_{k} = t_{r(j)} \left[\frac{\partial \varUpsilon_j}{\partial t}+ \frac{\partial}{\partial x_m} (c_m \varUpsilon_j)\right]_{X_j=0,\ \forall j \neq k}, \end{equation}which was proposed as a modification to the iterative refinement procedure for obtaining the second-order correction to the distribution function through the Chapman–Enskog multi-scale expansion procedure given as
\begin{equation} \bar{f}_{2} = t_{r} \left[\frac{\partial}{\partial t}+ c_m\frac{\partial}{\partial x_m}\right]\bar{f}_{1}. \end{equation}
Note that (2.13) contains a single relaxation time for both the viscous and thermal processes. In this work, we revert back to the original iterative representation of the multi-scale Chapman–Enskog expansion procedure, i.e. (2.13), and incorporate the individual relaxation times that separate the time scales associated with each irreversible thermodynamic process. This modification amounts to replacing 
$t_r$ with 
$t_r(j)$ given as
\begin{align} \left.\begin{gathered} \bar{f}_{2} = t_{r(j)} \left[\frac{\partial}{\partial t}+ \frac{\partial}{\partial x_m} c_m\right] \bar{f}_{1} = t_{r(j)} \left[\frac{\partial}{\partial t}+ \frac{\partial}{\partial x_m} c_m\right] \left\{\sum_j \varUpsilon_j \odot X_j \right\} \\ = \left\lbrace \sum_{j} \left[ t_{r(j)} \left(\frac{\partial \varUpsilon_j }{\partial t} + c_k \frac{\partial \varUpsilon_j } {\partial x_k}\right) \odot X_j\right] \right.\\ \left.+\sum_{j} \left[t_{r(j)} \varUpsilon_j \odot \left(\frac{\partial X_j }{\partial t} + c_k \frac{\partial X_j }{\partial x_k}\right)\right]\right\rbrace_{X_i=0,\ \forall i \neq j} \\ \implies \bar{f}_{2} = \left\{\underbrace{\sum_{j} [{\varUpsilon}_{kj}\odot {X}_{k} \odot X_j]}_{\bar{f}_{2,1}} + \underbrace{\sum_{j} \left[t_{r(j)} \varUpsilon_j \odot \left(\frac{\partial X_j }{\partial t} + c_k \frac{\partial X_j}{\partial x_k}\right)\right]}_{\bar{f}_{2,2}}\right\}_{X_i=0,\ \forall i \neq j}. \end{gathered}\right\} \end{align}
Here, we remark that (2.14) contains (2.12) and includes additional analytically obtained 
$\mathcal {O}$(
$Kn^{2}$) terms. Form (2.8), we can further define 
$\bar {f}_{2,1}$ as
\begin{equation} \bar{f}_{2,1}= [({\varUpsilon}_{\tau\tau} \odot {X}_\tau) \odot X_\tau] + [({\varUpsilon}_{qq} \odot {X}_q) \odot X_q]. \end{equation}Equation (2.15) is the original form of (2.8), represented here as
\begin{equation} \bar{f}^{'}_{2,1}= [({\varUpsilon''_{\tau\tau}} \odot {X}_\tau) \odot X_\tau] + [({\varUpsilon''_{qq}} \odot {X}_q) \odot X_q]. \end{equation}
The term 
$\bar {f}_{2,2}$ is now expanded below as
$$\begin{gather} {\bar{f}_{2,2}} = \sum_{j} \left[t_{r(j)} \varUpsilon_j \odot \left(\frac{\partial X_j}{\partial t} + c_k \frac{\partial X_j } {\partial x_k}\right)\right]_{X_i=0,\ \forall i \neq j}, \end{gather}$$
$$\begin{gather}\implies {\bar{f}_{2,2}} = [{\varUpsilon}_{\tau}\odot({X}_{\tau\tau}\odot {X}_{\tau})]+ [{\varUpsilon}_{q}\odot({X}_{qq}\odot {X}_{q})]. \end{gather}$$
In order to obtain the explicit expression represented in (2.17b), we first represent (2.17a) in terms of the material derivative 
${\textrm {D}}/{\textrm {D} t}$ and the peculiar velocity as
\begin{equation} \bar{f}_{2,2} = t_{r(j)} \varUpsilon_j \odot \left[{\frac{{\rm D} X_j}{{\rm D} t}}+ {C_k \frac{\partial}{\partial x_k} ( X_j)}\right]_{X_i=0,\ \forall i \neq j}. \end{equation}To further simplify (2.18), the material derivative present in (2.18) can be expressed in terms of spatial derivatives using the Euler and N-S equations.
Upon using the Euler equations to represent the material derivatives, we obtain
\begin{align}
{X}_{\tau\tau}\odot {X}_{\tau} &= t_{r(\tau)}
\left(\underbrace{\left[\frac{1}{\beta} \frac{\partial
\beta }{\partial x_j}\frac{\partial g}{\partial x_i} -
\frac{\partial}{\partial x_j} \frac{\partial g}{\partial
x_i}\right]}_{\omega_{10}}
-2\beta\underbrace{\frac{\partial u_i}{\partial x_k}
\frac{\partial u_k}{\partial x_j}}_{\omega_{11}}+
\underbrace{\frac{4}{3} \beta\frac{\partial u_l}{\partial
x_l} \frac{\partial u_i}{\partial x_j} }_{\omega_{12}}
\right. \nonumber\\ &\quad + \left.C_l
\left[\underbrace{2\beta \frac{\partial } {\partial x_l}
\frac{\partial u_i}{\partial x_j} +2\frac{\partial \beta }
{\partial x_l} \frac{\partial u_i}{\partial
x_j}}_{\omega_{13}}\right]\right),
\end{align}and
\begin{equation} {X}_{qq}\odot {X}_{q} = t_{r(q)} \left(\underbrace{\frac{2}{3} \left[\frac{\partial \beta}{\partial x_i} \frac{\partial u_l}{\partial x_l} + \beta \frac{\partial }{\partial x_i} \frac{\partial u_l }{\partial x_l}\right]}_{\xi_{11}} - \underbrace{\frac{\partial \beta}{\partial x_k} \frac{\partial u_k}{\partial x_i}}_{\xi_{12}} + \underbrace{C_l \left[\frac{\partial } {\partial x_l} \frac{\partial \beta}{\partial x_i}\right]}_{\xi_{13}}\right). \end{equation}
With (2.19) and (2.20), we now have a complete 
$\mathcal {O}$(
$Kn^{2}$) extended representation of the distribution function. Upon checking if this representation satisfies the constraint of additive invariance, we find that only terms 
$\omega _{13}$ and 
$\xi _{13}$ do not satisfy the last two compatibility conditions, and thus require modification. Following previous works (Balakrishnan, Agarwal & Yun Reference Balakrishnan, Agarwal and Yun1999; Agarwal et al. Reference Agarwal, Yun and Balakrishnan2001; Yadav et al. Reference Yadav, Jonnalagadda and Agrawal2023), we represent (2.19) and (2.20) as
\begin{equation} \langle \boldsymbol{c}, \bar{f}_{2,2}\rangle + t_{r(\tau)} \left(\sum^{13}_{i=10} \alpha_i \langle c_i, \varUpsilon_{\tau} \odot \omega_i \rangle \right) +t_{r(q)} \left(\sum^{13}_{i=11} \beta_i \langle c_i, \varUpsilon_{q} \odot \xi_i \rangle\right)=0 ,\end{equation}and
\begin{align} \left\langle \frac{\lvert\boldsymbol{C}\rvert^2}{2}, \bar{f}_{2,2} \right\rangle + t_{r(\tau)} \left(\sum^{13}_{i=10} \alpha_i \left \langle \frac{\lvert \boldsymbol{C} \rvert^2}{2}, \varUpsilon_{\tau} \odot \omega_i \right \rangle \right) +t_{r(q)} \left(\sum^{13}_{i=11} \beta_i \left \langle \frac{\lvert \boldsymbol{C} \rvert^2}{2}, \varUpsilon_{q} \odot \xi_i \right \rangle\right)=0, \end{align}
where the running index 
$i$ represents each term appearing in (2.19) and (2.20), and 
$\alpha _i$ and 
$\beta _i$ are moment closure coefficients. Note that, similar to the equations presented in Yadav et al. (Reference Yadav, Jonnalagadda and Agrawal2023), (2.21)–(2.22) form an unconstrained system of algebraic equations and can therefore result in multiple acceptable solutions. One solution is evaluated as follows:
\begin{equation} \alpha_i = \begin{cases} \varPsi_1\lvert\boldsymbol{C}\rvert^2 \beta+\varPsi_2 & \text{if}\ i=13 \\ 0 & \text{otherwise}, \end{cases} \end{equation}and
\begin{equation} \beta_i = \begin{cases} \dfrac{-1}{\zeta_c}\left(\zeta_1\dfrac{7}{2\beta} - \zeta_2 \lvert \boldsymbol{C} \rvert^2 + \zeta_3\dfrac{\lvert \boldsymbol{C} \rvert}{\sqrt{\beta}} \right) & \text{if}\ i=13 \\ 0 & \text{otherwise}, \end{cases} \end{equation}
where 
$\varPsi _1 = -\frac {2}{7} (\varPsi _{2} + 1)$, 
$\zeta _c = ( \lvert \boldsymbol {C} \rvert ^2 -{5}/{2 \beta })$, 
$\zeta _1 = - {8 \zeta _{3}}/{21 \sqrt {{\rm \pi} }} - \frac {5}{7}$ and 
$\zeta _2 = {8 \beta _{3}}/ {15 \sqrt {{\rm \pi} }} - 1$.
Consequently, based on the additive invariance condition and stability criteria, we take 
$\varPsi _2 = -1$ and 
$\zeta _3 = {27\sqrt {{\rm \pi} }}/{2}$, resulting in 
$\alpha _{13} = -1$ and 
$\beta _{13} = ({1}/{\zeta _c})({31 C^{2}}/{5} - {27 \sqrt {{\rm \pi} } C}/{2 \sqrt {\beta }} + {41}/{2 \beta })$, and then obtain the modified forms of the 
$\mathcal {O}(Kn^2)$ representations of (2.19) and (2.20) as
$$\begin{gather} \boldsymbol{X}^{'}_{\tau\tau}\odot \boldsymbol{X}_{\tau} = \boldsymbol{X}_{\tau\tau}\odot \boldsymbol{X}_{\tau} + t_{r(\tau)} \left(C_l \underbrace{\left[2\beta \frac{\partial}{\partial x_l} \frac{\partial u_i}{\partial x_j} +2\frac{\partial \beta}{\partial x_l} \frac{\partial u_i}{\partial x_j}\right]}_{\omega_{13}}\right)\alpha_{13}, \end{gather}$$
$$\begin{gather}\boldsymbol{X}^{'}_{qq}\odot \boldsymbol{X}_{q}= \boldsymbol{X}_{qq}\odot \boldsymbol{X}_{q}+ t_{r(q)} \left(\underbrace{C_l \left[\frac{\partial } {\partial x_l} \frac{\partial \beta}{\partial x_i}\right]}_{\xi_{13}}\right) \beta_{13}. \end{gather}$$
We note that, by replacing the material derivatives appearing in (2.18) using the N-S equations, we obtain 
$\mathcal {O}(Kn^3$) representations of 
$\bar {f}_{2,2}$, which are given as
$$\begin{gather} \boldsymbol{X}^{''}_{\tau\tau}\odot \boldsymbol{X}_{\tau} = \boldsymbol{X}^{'}_{\tau\tau}\odot \boldsymbol{X}_{\tau}+ t_{r(\tau)} \left\{2\beta \frac{\partial }{\partial x_j}\left(-\frac{1}{\rho} \frac{\partial \sigma^{NS}_{ik}} {\partial x_k}\right) + \frac{8\beta^2}{3\rho} \varOmega \frac{\partial u_i}{\partial x_j}\right\}, \end{gather}$$
$$\begin{gather}\boldsymbol{X}^{''}_{qq}\odot \boldsymbol{X}_{q} = \boldsymbol{X}^{'}_{qq}\odot \boldsymbol{X}_{q}+t_{r(q)} \left\{\frac{\partial }{\partial x_i} \left(\frac{4\beta^2}{3\rho} \varOmega \right)\right\}. \end{gather}$$Summation of (2.26a) and (2.26b) results in the final form as
\begin{equation} \bar{f}^{'}_{2,2} = \boldsymbol{X}^{''}_{\tau\tau}\odot \boldsymbol{X}_{\tau}+ \boldsymbol{X}^{''}_{qq}\odot \boldsymbol{X}_{q}.\end{equation}This completes the derivation of the complete second- and third-order representations of the distribution function at the Burnett level, which final form of mathematical expression has been obtained by summation of (2.2), (2.3) along with (2.16) and (2.27) as
\begin{equation} {{f}_{2}} = f_0 + \bar{f}_1+ \bar{f}^{'}_{2,1} + \bar{f}^{'}_{2,2}. \end{equation}In summary, we presented the generalized form of a single-particle distribution function (2.28) by employing the iterative refinement technique. This distribution function core lies in the principles of non-equilibrium thermodynamics and satisfies the required physics (Singh & Agrawal Reference Singh and Agrawal2016; Singh et al. Reference Singh, Jadhav and Agrawal2017; Yadav et al. Reference Yadav, Jonnalagadda and Agrawal2023). Therefore, we combined both the Chapman–Enskog expansion (Chapman & Cowling Reference Chapman and Cowling1970; Cercignani Reference Cercignani1975) and Onsager's reciprocity principle (Onsager Reference Onsager1931a,Reference Onsagerb) for the first time. This function will be used to evaluate the super-Burnett-level stress tensor and heat flux vector to close the governing equations in the subsequent section.
3. Higher-order closure of governing equations
The distribution function derived in § 2 will now be used to obtain the higher-order constitutive relationships for the stress tensor and heat flux vector appearing in the mass, momentum and energy conservation equations given below
$$\begin{gather} \frac{\partial \rho}{\partial t}+\frac{\partial \rho u_k}{\partial x_k} =0, \end{gather}$$
$$\begin{gather}\rho \frac{\partial u_i}{\partial t} +\rho u_k \frac{\partial u_i}{\partial x_k} +\frac{\partial p}{\partial x_i} +\frac{\partial \sigma_{ik}} {\partial x_k}=\rho G_i, \end{gather}$$
$$\begin{gather}\rho \frac{\partial\epsilon}{\partial t} + \rho u_k \frac{\partial\epsilon}{\partial x_k}+\frac{\partial q_k}{\partial x_k} +p\frac{\partial u_k}{\partial x_k} + \sigma_{ij}\frac{\partial u_i}{\partial x_j}=0, \end{gather}$$
where 
$\varepsilon = (3/2) R T$ represents the internal energy for monatomic gases, 
$G_i$ is the external body force and 
$p$ is the pressure. Here, we evaluate the required constitutive relationships for the stress tensor and heat flux vector accurate up to 
$\mathcal {O}(Kn^3)$ using the final form of the single-particle distribution function presented in (2.28). The kinetic theory definition of the non-conserved stress tensor and heat flux vectors are given as
$$\begin{gather} \sigma_{ij} = m \int C_{\langle i}C_{j\rangle} f \,{\rm d}\boldsymbol{c}, \end{gather}$$
$$\begin{gather}q_i = \frac{m}{2} \int \lvert{\boldsymbol{C}}\rvert^2 C_i f \,{\rm d}\boldsymbol{c}. \end{gather}$$Upon solving the integrals, explicit constitutive relationships for the stress tensor and heat flux vector are obtained as
\begin{align} \sigma^{\mathcal{O}(Kn^3)}_{ij} &={-}2\mu\frac{\partial u_{\langle i}}{\partial x_{i \rangle}}+ t_{r_\tau}^2{2 p}\left[\left(\frac{5}{3}-\frac{2\varphi}{3}\right) + \frac{2}{3}\right]\frac{\partial u_{\langle i}}{\partial x_{j \rangle}} \frac{\partial u_k}{\partial x_k} - t_{r_\tau}^2 2R \frac{\partial T }{\partial x_{\langle i}} \frac{\partial p}{\partial x_{j \rangle}} \nonumber\\ &\quad - t_{r_\tau}^2{2pRT}\frac{\partial}{\partial x_{\langle i}} \frac{1}{p}\frac{\partial p} {\partial x_{j \rangle}} -t_{r_\tau}^2{2p} \frac{\partial u_{k}}{\partial x_{\langle i}} \frac{\partial u_{j \rangle}}{\partial x_{k}} + \eta \left( t_{r_q}^2 \frac{ 4 p R }{9} \frac{-1}{2} \frac{\partial}{\partial x_{\langle i}} \frac{\partial T}{\partial x_{j \rangle}}\right) \nonumber\\ &\quad + \eta \left( t_{r_q}^2 \frac{ 4 \rho R^2 }{9} \frac{\partial T }{\partial x_{\langle i}} \frac{\partial T}{\partial x_{j \rangle}} \right) + \underbrace{t_{r_\tau}^2 {2 p } \left\{\frac{\partial }{\partial x_{j \rangle}}\left(-\frac{1}{\rho} \frac{\partial \sigma^{NS}_{\langle ik}} {\partial x_k}\right) \right\}}_{\mathcal{O}(Kn^3)}, \end{align}
\begin{align} q^{\mathcal{O}(Kn^3)}_i &={-}\kappa \frac{\partial T}{\partial x_i} -t_{r_\tau}^2 {2 R T} \frac{\partial p}{\partial x_k}\frac{\partial u_{\langle k}}{\partial x_{i \rangle}} + t_{r_\tau}^2 {4 \rho R^2 T}\left(\frac{7}{4}-\frac{q_{c1}}{2}\right)\frac{\partial T}{\partial x_k} \frac{\partial u_{\langle k}}{\partial x_{i \rangle}} \nonumber\\ &\quad +t_{r_q}^2\frac{7\rho R^2 T}{2}\frac{\partial u_{k}}{\partial x_{i}}\frac{\partial T}{\partial x_k} + t_{r_q}^2 {\rho R^2 T}\frac{\partial u_{i}}{\partial x_{k}} \frac{\partial T}{\partial x_k}- {4 \rho R^2 T }t_{r_q}^2\left(\frac{10\varphi-1}{24} + \frac{5}{8}\right)\frac{\partial T}{\partial x_i}\frac{\partial u_j}{\partial x_j} \nonumber\\ &\quad +t_{r_\tau}^2 {q_{c2} 2 \rho R^2 T^2 }\frac{\partial }{\partial x_k} \frac{\partial u_{\langle i}}{\partial x_{k \rangle}} - t^2_{r_q}\frac{10 p R^2 T^2 }{3}\frac{\partial }{\partial x_i} \left(\beta\frac{\partial u_j}{\partial x_j}\right) \nonumber\\ &\quad - \underbrace{t^2_{r_q}{5 p R^2T^2} \left\{\frac{\partial }{\partial x_i} \left(\frac{1}{3 p R T} \varOmega \right)\right\}}_{\mathcal{O}(Kn^3)}, \end{align}where terms with angular brackets represent trace-free symmetric tensor quantities and
$$\begin{gather} \eta =\frac{9 (35 \sqrt{\rm \pi} \beta_{1} - 35 \sqrt{\rm \pi} \beta_{2} + 32 \beta_{3} - 10 \sqrt{\rm \pi})}{10 \sqrt{\rm \pi}}, \end{gather}$$
$$\begin{gather}q_{c1} = q_{c2} = \frac{7}{2} \left(\frac{9 \varPsi_{1}}{2} + \varPsi_{2} + 1 \right). \end{gather}$$
The coefficients 
$\eta$ and 
$q_{c1}$ = 
$q_{c2}$ depend upon the values of 
$\varPsi _1$ and 
$\varPsi _2$ as mentioned in § 2. The first terms in (3.6)–(3.7) are of 
$\mathcal {O}(Kn)$ while those appearing at the end are of 
$\mathcal {O}(Kn^3)$; the remaining terms are of 
$\mathcal {O}(Kn^2)$. Note that the stress tensor (3.6), similar to heat flux (3.7), includes both relaxation times corresponding to the momentum and energy transport, respectively. These constitutive relationships have several linear terms of the second- and third-order, which is not true in the case of the OBurnett equations (Singh et al. Reference Singh, Jadhav and Agrawal2017). In (3.7), nonlinear third-order terms are also present. Due to the extreme complexity, these nonlinear terms will not be considered in the stability and validation analysis presented in subsequent §§ 4–5; a similar consideration was employed earlier by Shavaliyev (Reference Shavaliyev1993), in which only linear terms of 
$\mathcal {O}(Kn^3)$ have been considered in the constitutive relationships.
Substituting the obtained expressions of the SOBurnett stress tensor (3.6) and heat flux vector (3.7) in (3.1)–(3.3) closes and completes the proposed governing set of SOBurnett equations. In the subsequent section, we analyse the linear stability properties of these equations.
4. Linear stability analysis of the derived equations
In this section, we perform the linear stability analysis of the derived equations, which is an important exercise because it was noted in Uribe et al. (Reference Uribe, Velasco and Garcia-Colin2000), Bobylev (Reference Bobylev2006) and Welder, Chapman & Maccormack (Reference Welder, Chapman and Maccormack1993) that the linearized forms of governing equations become unstable as the mesh size gradually reduces. This instability was noticed during numerical simulations of rarefied flow in the continuum–transition regimes, which stems from additional terms added to the constitutive relationships. The earlier equations, therefore, fail to capture the thermodynamic aspect of the flow accurately. As a result, the newly derived equations in the present work, having several additional terms in the constitutive relationships, must undergo a linearized stability analysis before being used to analytically solve the pressure-driven plane Poiseuille flow problem in §,5.
For this purpose, we first simplify the governing equations (3.1)–(3.3) by substituting constitutive relationships in the 2-D form from (A1)–(A2). In the second step, for the linearized and non-dimensionalized equations (3.1)–(3.3), we substitute these variables
\begin{equation} \left.\begin{gathered} \rho = \rho_o(1+\bar{\rho}),\quad T=T_o(1+\bar{T}),\quad p = \rho_o R T_o(1+\bar{p}),\quad u=\sqrt{RT_o}\bar{u},\\ v=\sqrt{RT_o}\bar{v},\quad x = H \bar{x},\quad y = H \bar{y},\quad t = \frac{H}{\sqrt{RT_o}}\bar{t}, \end{gathered}\right\} \end{equation}
where subscript ‘
$o$’ represents the variable at the equilibrium state, 
$H$ represents the characteristic length scale and superscript ‘
$\bar {}$’ denotes small perturbations around the equilibrium (rest) state. Finally, only linear terms are considered for the analysis because the perturbations have been assumed to be small, thereby simplifying the analysis to a great extent. This simplification results in the following reduced linearized and non-dimensionalized equations:
\begin{gather} \frac{\partial \bar{\rho}}{\partial \bar{t}} + \frac{\partial \bar{u}}{\partial \bar{x}} + \frac{\partial \bar{v} }{\partial \bar{y}} = 0, \end{gather}
\begin{gather} \begin{aligned}& 3 \frac{\partial \bar{T}}{\partial \bar{x}} - 4 \frac{\partial^{3}\bar{T}}{\partial \bar{x}^{3}} + 3 \frac{\partial \bar{\rho}}{\partial \bar{x}} - 4 \frac{\partial^{3} \bar{\rho}}{\partial \bar{x}^{3}} + 3 \frac{\partial \bar{u}}{\partial \bar{t}} - 4 \frac{\partial^{2} \bar{u}}{\partial \bar{x}^{2}} + \underline{\frac{16}{3} \frac{\partial^4 \bar{u}}{\partial \bar{x}^4}}- 3 \frac{\partial^{2} \bar{u} }{\partial \bar{y}^{2}}- 4 \frac{\partial^{3} \bar{T}}{\partial \bar{x}\partial \bar{y}^{2}} - 4 \frac{\partial^{3} \bar{\rho}}{\partial \bar{x}\partial \bar{y}^{2}} \nonumber\\ &\quad -\frac{\partial^{2} \bar{v}}{\partial \bar{x}\partial \bar{y}}- \underline{\frac{8}{3} \frac{\partial^4 \bar{v}}{\partial \bar{x} \partial \bar{y}^3}}= 0, \end{aligned}\end{gather}
\begin{gather} \begin{aligned}& 3 \frac{\partial \bar{T}}{\partial \bar{y}} - 4 \frac{\partial^{3} \bar{T}}{\partial \bar{y}^{3}} + 3 \frac{\partial \bar{\rho}}{\partial \bar{y}} - 4 \frac{\partial^{3} \bar{\rho}}{\partial \bar{y}^{3}} + 3 \frac{\partial \bar{v}}{\partial \bar{t}} - 3 \frac{\partial^{2} \bar{v}}{\partial \bar{x}^{2}} - 4 \frac{\partial^{2} \bar{v}}{\partial \bar{y}^{2}} +\underline{\frac{16}{3} \frac{\partial^4 \bar{v}}{\partial \bar{y}^4}} - 4 \frac{\partial^{3} \bar{T}}{\partial \bar{x}^{2}\partial \bar{y}} - 4 \frac{\partial^{3} \bar{\rho}}{\partial \bar{y}\partial \bar{x}^{2}} \nonumber\\ &\quad -\frac{\partial^{2} \bar{u}}{\partial \bar{y}\partial \bar{x}} -\underline{\frac{8}{3} \frac{\partial^4 \bar{u}}{\partial \bar{x}^3 \partial \bar{y}}} = 0, \end{aligned}\end{gather}
\begin{gather} \begin{aligned}& 18 \frac{\partial \bar{T}}{\partial \bar{t}} - 45 \frac{\partial^{2} \bar{T}}{\partial \bar{x}^{2}} + \underline{\frac{675}{4}\frac{\partial^{4} \bar{T}}{\partial \bar{x}^{4}}} - 45 \frac{\partial^{2} \bar{T}}{\partial \bar{y}^{2}} + \underline{\frac{675}{4} \frac{\partial^{4} \bar{T}}{\partial \bar{y}^{4}}} + 12 \frac{\partial \bar{u}}{\partial \bar{x}} - 45 \frac{\partial^{3}\bar{u}}{\partial \bar{x}^{3}}+ 12 \frac{\partial \bar{v} }{\partial \bar{y}}- 45 \frac{\partial^{3} \bar{v}}{\partial \bar{y}^{3}} \nonumber\\ &\quad - 45 \frac{\partial^{3} \bar{u}}{\partial \bar{y}^{2}\partial \bar{x}} - 45 \frac{\partial^{3} \bar{v}}{\partial \bar{y}\partial \bar{x}^{2}} = 0. \end{aligned}\end{gather}Note that super-Burnett-level terms, denoted by the underlined terms, have also been considered in (4.2)–(4.5). We now apply the method of normal modes to the perturbation
\begin{equation} \left.\begin{gathered} \bar{\rho}=\rho_A \exp{(\omega \bar{t}+{\rm i} k \bar{x}+{\rm i} k \bar{y})},\quad \bar{u}= u_A \exp{(\omega \bar{t}+{\rm i} k \bar{x}+{\rm i} k \bar{y})},\\ \bar{v}= v_A \exp{(\omega \bar{t}+{\rm i} k \bar{x}+{\rm i} k \bar{y})}, \quad \bar{T}= T_A \exp{(\omega \bar{t}+{\rm i} k \bar{x}+{\rm i} k \bar{y})}, \end{gathered}\right\} \end{equation}
where 
$k$, 
$\omega$ and subscript 
$A$ represent the wavenumber, wave frequency and complex amplitude of the plane wave, respectively. The substitution of these solutions in (4.2)–(4.5) results in a relation between 
$k$ and 
$\omega$, known as the dispersion relation, as
\begin{align} & 162 \omega^{4} + \left(\frac{7227 \kappa^{4}}{2} + 1566 \kappa^{2}\right) \omega^{3}+ (11\,184 \kappa^{8} + 29\,727 \kappa^{6} + 10\,296 \kappa^{4} + 540 \kappa^{2})\omega^{2} \nonumber\\ &\quad + (7200 \kappa^{12} + 58\,872 \kappa^{10} + 77\,136 \kappa^{8} + 27\,459 \kappa^{6} + 2700 \kappa^{4})\omega \nonumber\\ &\quad + (43\,200 \kappa^{12} + 60\,120 \kappa^{10} + 25\,110 \kappa^{8} + 3240 \kappa^{6})=0. \end{align}Additionally, following similar steps as mentioned above and neglecting super-Burnett-order terms gives the dispersion relation for the EOBurnett equations as
\begin{align} & 162 \omega^{4} + 1566 \omega^{3} \kappa^{2} + (11\,232 \kappa^{6} + 10\,296 \kappa^{4} + 540 \kappa^{2})\omega^{2} \nonumber\\ &\quad + (22\,464 \kappa^{8} + 19\,944\kappa^{6} + 2700 \kappa^{4}) \omega + (8640 \kappa^{8} + 3240 \kappa^{6})=0. \end{align}
Equations (4.7)–(4.8) are fourth-order polynomials in 
$\omega$ and give four roots having both real and complex parts. Using the solution of (4.7)–(4.8), we now test the stability of the proposed equations for a disturbance in space. Since we consider disturbance in space, the wavenumber is real, and the complex frequency is given by 
$\omega = \omega _r(k) + \textrm {i} \omega _i(k)$. The stability of the solution requires 
$\omega _r(k)\le 0$ so that the local amplitude of primary variables decreases with an increase in time.
Figure 1 demonstrates that the obtained results satisfy the stability condition, which requires that the solutions always have a negative real part. Similar to Balakrishnan et al. (Reference Balakrishnan, Agarwal and Yun1999), we further consider the definition of the Knudsen number as 
$Kn \approx k$ for a given wavelength. Consequently, we present the variation of the attenuation coefficient (
$\omega _r$) against the Knudsen number in figure 2, showing that all four roots of the solutions remain negative, unlike the conventional and super-Burnett equations (Bobylev Reference Bobylev1982).
From stability analysis, we demonstrate that the equations derived in this study are unconditionally stable for disturbances in space. This stability is a crucial characteristic that will be utilized in § 5 to address the pressure-driven plane Poiseuille flow problem to validate these equations in the transition flow regime.
5. Analytical solution of pressure-driven Poiseuille flow
In this section, we investigate the long micro-channel gas flow problem (see figure 3) using the derived sets of both EOBurnett and SOBurnett equations.
Figure 3. Microchannel schematic showing the 
$x$-
$y$ axis, streamwise velocity for steady 2-D and isothermal flow.
We assume that the flow is steady, two-dimensional, isothermal and free of any external body force. These assumptions simplify the present problem so that an analytical solution can be obtained. Consequently, we first provide the simplified form of the governing equations (3.1)–(3.2) by employing the aforementioned assumptions, which is followed by a discussion on the boundary conditions applied. Subsequently, the analytical solutions of the pressure and velocity fields have been obtained using the simplified governing equations, which are finally validated against DSMC data.
5.1. Simplified governing equations
Similar to Rath et al. (Reference Rath, Singh and Agrawal2018, Reference Rath, Yadav and Agrawal2021), we first substitute 2-D expressions of the stress tensor from (A1) upon applying the aforementioned assumptions in (3.1)–(3.2) and then non-dimensionalize the same equations using 
$\bar {x}={x}/{L}$, 
$\bar {y}={y}/{H}$, 
$\bar {u}={u}/{u_{out}}$, 
$\bar {p}={p}/{p_{out}}$ and 
$\bar {\rho }={\rho }/{\rho _{out}}$, where 
$L$ is length and 
$H$ is the height of the microchannel and ‘out’ represents the reference location considered at the microchannel's exit. Additionally, the reference velocity 
$u_{out}$ and reference pressure 
$p_{out}$ are utilized to normalize the pressure and velocity, respectively. Next, a smallness parameter 
$\epsilon ={H}/{L}$ has been defined. Further, the assumptions and parameters mentioned above have been used to express the governing equations in terms of the Mach number (
$Ma$) and Reynolds number (
$Re$) given in (A4)–(A5). While preserving the Knudsen number at order unity, we study gaseous flow under low Mach number (
$O(\epsilon )$) and low Reynolds number (
$O(\epsilon )$) conditions. These crucial assumptions have been further applied to the continuity and momentum equations (A4)–(A5) to eliminate terms of order larger than 
$\epsilon$, resulting in the following simplified stream and cross-stream momentum equations, respectively:
\begin{align} & \frac{Ma^{2}
\gamma}{Re^{2} \bar{p}^{2}}
\frac{\partial^{4}\bar{u}}{\partial \bar{y}^{4}} - \frac{4
\epsilon }{3 Re \bar{p}^{2}}\frac{\partial^{3}\bar{p}
}{\partial \bar{x}\partial \bar{y}^{2}} \nonumber\\ &\quad
+ \frac{2 \epsilon }{3 Re
\bar{p}^{3}}\frac{\partial\bar{p}}{\partial \bar{x}}
\frac{\partial^{2}\bar{p}}{\partial \bar{y}^{2}} + \frac{14
\epsilon }{3 Re \bar{p}^{3}}
\frac{\partial\bar{p}}{\partial \bar{y}}
\frac{\partial^{2}\bar{p}}{\partial \bar{x}\partial
\bar{y}} \nonumber\\ &\quad - \frac{4 \epsilon }{Re
\bar{p}^{4}}\frac{\partial\bar{p}}{\partial \bar{x}}
\left(\frac{\partial \bar{p}}{\partial \bar{y}}\right)^{2}
+ \frac{Re \epsilon }{Ma^{2} \gamma} \frac{\partial
\bar{p}}{\partial \bar{x}} =
\frac{\partial^{2}\bar{u}}{\partial \bar{y}^{2}},
\end{align}
\begin{equation} - \frac{4 }{3 Re \bar{p}^{2}}\frac{\partial^{3}\bar{p} }{\partial \bar{y}^{3}} + \frac{16 }{3 Re \bar{p}^{3}}\frac{\partial\bar{p}}{\partial \bar{y}} \frac{\partial^{2}\bar{p}}{\partial \bar{y}^{2}} - \frac{4 }{Re \bar{p}^{4}} \left(\frac{\partial \bar{p}}{\partial \bar{y}} \right)^{3}+ \frac{Re }{Ma^{2} \gamma} \frac{\partial\bar{p}}{\partial \bar{y}} =0. \end{equation}
The coupled equations presented in (5.1)–(5.2) involve two variables, namely 
$\bar {u}$ and 
$\bar {p}$, with the presence of a fourth-order derivative of velocity and third-order derivative of pressure. As a result, four and three integration constants are produced while solving for the velocity and pressure, respectively. These integration constants need to be evaluated using appropriate boundary conditions, as discussed in the subsequent subsection.
5.2. Boundary/initial conditions
It is crucial to specify accurate boundary and/or initial conditions to evaluate the integration constants. To describe these conditions, a model for the interaction between the fluid and the solid wall is needed, which is still unknown even for the Burnett equations (García-Colín et al. Reference García-Colín, Velasco and Uribe2008). In such a case, Uribe & Garcia (Reference Uribe and Garcia1999) has demonstrated a method to deal with the scatter in DSMC data and obtain these unknown conditions that are useful for both initial and boundary value problems. Therefore, instead of getting bogged down by the unavailability of accurate boundary conditions, treating the problem as an initial value problem has been suggested in the literature (Singh et al. Reference Singh, Gavasane and Agrawal2014, Reference Singh, Jadhav and Agrawal2017; Jadhav et al. Reference Jadhav, Singh and Agrawal2017; Rath et al. Reference Rath, Singh and Agrawal2018, Reference Rath, Yadav and Agrawal2021; Yadav et al. Reference Yadav, Jonnalagadda and Agrawal2023). This helps to bypass the unknown boundary conditions when initial conditions are known from the DSMC method. Therefore, we follow the same approach demonstrated in these references for evaluating the integration constants using DSMC data. This is both logical and essential in the present case since the derivation of boundary conditions is beyond the scope of the present work. However, we have also included a few results for the case when the present problem has been treated as a boundary value problem.
There are several advantages to considering the current problem as an initial value problem. First, the solution should automatically capture the DSMC data near the wall in the Knudsen layer, as we do not assume anything about the slip model. Second, measuring any quantity at the centre is easier than near the wall in an experiment. Therefore, this justifies the approach adopted in our present study and helps us to test the proposed equations for the problems considered in the present work, which otherwise would not be possible. Consequently, all integration constants coming from (5.2) are being obtained using initial conditions 
${\partial \hat {p}}/{\partial \bar {y}}|_{\bar {y}=0}=0$, 
$\hat {p}|_{\bar {y}=0} =\hat {p}_c(\bar {x})$ and 
${{\partial ^2 \hat {p}}/{\partial \bar {y}^2}|_{y=0}}=\hat {a}(\bar {x})$, where 
$\hat {p}_c$ is the centreline pressure and 
$\hat {a}$ is the curvature of the pressure at the centreline. Similarly, integration constants coming from (5.1) have been obtained using initial conditions 
${\partial \bar {u}}/{\partial \bar {y}}|_{\bar {y}=0}$, 
$\bar {u}|_{\bar {y}=0},\ {\partial ^2 \bar {u}}/{\partial \bar {y}^2}|_{\bar {y}=0}$ and 
${\partial ^3 \bar {u}}/{\partial \bar {y}^3}|_{\bar {y}=0}$ at the centre of a channel from DSMC data. These initial conditions will be used in the subsequent subsection presenting the pressure and velocity field profiles in the transition flow regime.
5.3. Derivation of analytical solution
To solve the coupled nonlinear partial differential equations (5.1)–(5.2) for the pressure and velocity fields, we use the perturbation method to simplify these two equations further, as demonstrated in Rath et al. (Reference Rath, Singh and Agrawal2018). As a result, we perturb pressure (
$\bar {p}$) around 
$\bar {p}_{NS}(\bar {x})$
\begin{equation} \bar{p}(\bar{x},\bar{y}) = \bar{p}_{NS}(\bar{x})+\hat{p}(\bar{x},\bar{y}), \end{equation}
where 
$\hat {p}(\bar {x},\bar {y})$ is the perturbed pressure term, which is small in magnitude. This suggests that 
$\hat {p}/\bar {p}_{NS} \ll 1$, which is true in the case of microchannel flow (Zheng et al. Reference Zheng, Garcia and Alder2002; Xu Reference Xu2003; Xu & Li Reference Xu and Li2004; Rath et al. Reference Rath, Singh and Agrawal2018, Reference Rath, Yadav and Agrawal2021). We further substitute the pressure in terms of a perturbed quantity using (5.3) in (5.1)–(5.2) and then consider only the first-order term of the perturbed pressure 
$\hat {p}$ and neglect any nonlinear term involving the derivative of both pressure and velocity with respect to the normal direction, as followed in Rath et al. (Reference Rath, Singh and Agrawal2018). These assumptions are reasonable in the case of microchannel gaseous flow when the flow is not far from equilibrium and the perturbed terms are rather small. Utilizing these assumptions finally results in a simplified version of the momentum equations
\begin{align} & \frac{Ma^{2} \gamma}{Re^{2} \bar{p}^{2}} \frac{\partial^{4}\bar{u}}{\partial \bar{y}^{4}} - \frac{\partial^{2} \bar{u}}{\partial \bar{y}^{2}} - \frac{4 \epsilon }{3 Re {\bar{p}_{NS}^{2}}} \frac{\partial^{3} \hat{p}}{\partial \bar{y}^{2}\partial \bar{x}} + \frac{2 \epsilon }{3 Re {\bar{p}_{NS}^{3}}} \frac{{\rm d} {\bar{p}_{NS}}}{{\rm d} \bar{x}} \frac{\partial^{2}\hat{p} }{\partial\bar{y}^{2}} \nonumber\\ &\quad + \frac{Re \epsilon }{Ma^{2} \gamma} \frac{{\rm d} {\bar{p}_{NS}}}{{\rm d} \bar{x}} + \frac{Re \epsilon }{Ma^{2} \gamma} \frac{\partial \hat{p}}{\partial \bar{x}} =0, \end{align}
\begin{equation} - \frac{2 \hat{\kappa}_{2} }{3 {\bar{p}^{2}_{NS}}}\frac{\partial^{3} \hat{p}}{\partial \bar{y}^{3}} + \frac{1}{\gamma}\frac{\partial\hat{p}}{\partial \bar{y}} = 0, \end{equation}
where 
$\hat {\kappa }_2$ is 
$2 Ma^2/Re^2$ (
$= 4Kn^2/({\rm \pi} \gamma )$). The solution of the cross-stream momentum equation (5.5) is first obtained as
\begin{equation} \hat{p} = {D_{1}} + {D_{2}} \exp({- A \bar{y} {\bar{p}_{NS}}}) + {D_{3}} \exp({A \bar{y} {\bar{p}_{NS}}}), \end{equation}
where 
$A = \sqrt {3/(2\hat {\kappa }_2\gamma )}$ and 
${D_{1}}({\bar {x}})-{D_{3}}({\bar {x}})$ are integration constants, which have been evaluated using initial conditions as discussed in § 5.2. The substitution of these integration constants in (5.6) then results in a closed-form expression of the pressure as
\begin{equation} \hat{p} = \hat{p}_{c} + \frac{\hat{a} \exp({A \bar{y} {\bar{p}_{NS}}})}{2 A^{2} {\bar{p}_{NS}^{2}}} - \frac{\hat{a}}{A^{2} {\bar{p}_{NS}^{2}}} + \frac{\hat{a} \exp({- A \bar{y} {\bar{p}_{NS}}})}{2 A^{2} {\bar{p}_{NS}^{2}}}. \end{equation}
Upon obtaining the solution of 
$\hat {p}$, substituting the result of (5.7) into (5.4) gives an analytical expression of the streamwise velocity, which read as
\begin{align} \bar{u} &= \bar{y} {C_2} + \frac{1}{2} \left(- \frac{4 \epsilon \exp({A \bar{y} {\bar{p}_{NS}}})}{3 A^{2} Re {\bar{p}_{NS}^{4}}} - \frac{4 \epsilon \exp({- A \bar{y} {\bar{p}_{NS}}})}{3 A^{2} Re {\bar{p}_{NS}^{4}}} - \frac{Re \epsilon \bar{y}^{2}}{A^{2} Ma^{2} \gamma {\bar{p}_{NS}^{2}}} + \frac{Re \epsilon \exp({A \bar{y} {\bar{p}_{NS}}})}{A^{4} Ma^{2} \gamma {\bar{p}_{NS}^{4}}} \right. \nonumber\\ &\quad \left.+ \frac{Re \epsilon \exp({- A \bar{y} {\bar{p}_{NS}}})}{A^{4} Ma^{2} \gamma {\bar{p}_{NS}^{4}}}\right) \frac{{\rm d} \hat{a} }{{\rm d} \bar{x}} + \frac{1}{2} \left( - \frac{4 \epsilon \bar{y} \hat{a} \exp({A \bar{y} {\bar{p}_{NS}}})}{3 A Re {\bar{p}_{NS}^{4}}} + \frac{4 \epsilon \bar{y} \hat{a} \exp({- A \bar{y} {\bar{p}_{NS}}})}{3 A Re {\bar{p}_{NS}^{4}}} \right. \nonumber\\ &\quad + \frac{10 \epsilon \hat{a} \exp({A \bar{y} {\bar{p}_{NS}}})}{3 A^{2} Re {\bar{p}_{NS}^{5}}} + \frac{10 \epsilon \hat{a} \exp({- A \bar{y} {\bar{p}_{NS}}})}{3 A^{2} Re {\bar{p}_{NS}^{5}}} + \frac{2 Re \epsilon \bar{y}^{2} \hat{a}}{A^{2} Ma^{2} \gamma {\bar{p}_{NS}^{3}}} + \frac{Re \epsilon \bar{y} \hat{a} \exp({A \bar{y} {\bar{p}_{NS}}})}{A^{3} Ma^{2} \gamma {\bar{p}_{NS}^{4}}} \nonumber\\ &\quad \left.-\frac{Re \epsilon \bar{y} \hat{a} \exp({- A \bar{y} {\bar{p}_{NS}}})}{A^{3} Ma^{2} \gamma {\bar{p}_{NS}^{4}}} -\frac{4 Re \epsilon \hat{a} \exp({A \bar{y} {\bar{p}_{NS}}})}{A^{4} Ma^{2} \gamma {\bar{p}_{NS}^{5}}} - \frac{4 Re \epsilon \hat{a} \exp({- A \bar{y} {\bar{p}_{NS}}})}{A^{4} Ma^{2} \gamma {\bar{p}_{NS}^{5}}}\right) \frac{{\rm d} {\bar{p}_{NS}}}{{\rm d} \bar{x}} \nonumber\\ &\quad +{C_{1}} +\frac{Re \epsilon \bar{y}^{2} }{2 Ma^{2} \gamma} \frac{{\rm d} \hat{p}_{c}}{{\rm d} \bar{x}}+ \frac{Re \epsilon \bar{y}^{2} }{2 Ma^{2} \gamma} \frac{{\rm d} {\bar{p}_{NS}}}{{\rm d} \bar{x}}, \end{align}
where SOBurnett-order terms have not been considered, and 
$C_1$–
$C_2$ are integration constants appearing due to the integration of the second-order differential equation. In contrast, considering SOBurnett-order terms present in (5.4) yields (B1), which is presented in Appendix B due to its long expression.
To find the integration constants (
$C_1$–
$C_2$) in (5.8) and additional integration constants (
$C_3$–
$C_4$) in (B1), we have employed the inherent symmetry 
${\partial \bar {u}}/{\partial \bar {y}}|_{\bar {y}=0}$ and other initial conditions such as 
$\bar {u}|_{\bar {y}=0},\ {\partial ^2 \bar {u}}/{\partial \bar {y}^2}|_{\bar {y}=0}$ and 
${\partial ^3 \bar {u}}/{\partial \bar {y}^3}|_{\bar {y}=0}$ at the centre of the channel, as discussed in § 5.2.
5.4. Validation of analytical solution using DSMC data
It is worth noting that the DSMC technique is already established as an essential tool for addressing rarefied gas flow problems. Akhlaghi, Roohi & Stefanov (Reference Akhlaghi, Roohi and Stefanov2012), Varade et al. (Reference Varade, Duryodhan, Agrawal, Pradeep, Ebrahimi and Roohi2015), Goshayeshi, Roohi & Stefanov (Reference Goshayeshi, Roohi and Stefanov2015) and Ebrahimi & Roohi (Reference Ebrahimi and Roohi2017) have conducted comprehensive studies on the DSMC technique applied to various micro- and nanochannel flow challenges, as well as the Knudsen pump. Baier et al. (Reference Baier, Hardt, Shahabi and Roohi2017) examined a Knudsen pump design influenced by the Crookes radiometer, utilizing vanes positioned within a temperature gradient. Additionally, a novel collision scheme, known as the simplified Bernoulli trial, has been introduced in Stefanov (Reference Stefanov2011). The DSMC technique has also been employed in Gavasane et al. (Reference Gavasane, Sachdev, Mittal, Bhandarkar and Agrawal2011, Reference Gavasane, Agrawal, Pradeep and Bhandarkar2017), Jadhav et al. (Reference Jadhav, Singh and Agrawal2017), Shah, Agrawal & Bhandarkar (Reference Shah, Agrawal and Bhandarkar2018a), Shah et al. (Reference Shah, Gavasane, Agrawal and Bhandarkar2018b) and Jadhav, Gavasane & Agrawal (Reference Jadhav, Gavasane and Agrawal2021) for different microchannel flow and shock wave problems. That is why we also resort to DSMC data to validate the analytical solution obtained above. Bird's DSMC code (Bird Reference Bird1994) is used to generate the DSMC data for the present problem, in which the sample size is 9 700 000. The smallness parameter 
$\epsilon = 1/10$, particle diameter 
$d = 4.17 \times 10^{-10}$ m, Boltzmann constant 
$k_b = 1.3806 \times 10^{-23}$ J K
$^{-1}$, dynamic viscosity 
$\mu = 2.274333 \times 10^{-5}$ N s m
$^{-2}$, gas constant for argon 
$R$ = 
$208$ J kg K
$^{-1}$. The exit velocity 
$u_{out}$, pressure 
$p_{out}$ and temperature 
$T_{out}$ are 
$122.5951$ m s
$^{-1}$, 
$1.3551\times 10^{4}$ Pa and 
$270.9344$ K, respectively. The inlet-to-outlet pressure 
$P_r$ and temperature 
$T_r$ ratios are 
$4.9466$ and 
$1.0335$, respectively. Moreover, the Knudsen number at the exit (
$Kn_{out}$) is 0.4994. Using these DSMC data, we evaluate the pressure and streamwise velocity for the present problem at two streamwise locations, 
$x/L = 0.1$ and 
$x/L = 0.8$, where 
$Kn|_{\bar {x}=0.1} = 0.1129$ and 
$Kn|_{\bar {x}=0.8} = 0.2578$. Moreover, as expected, the centreline pressure obtained from DSMC data 
$\bar {p}_c $ closely approximates 
$\bar {p}_{NS} $, with a difference of less than 1.25 observed at both locations, as indicated by 
$\frac{|\bar{p}_c - \bar{p}_{NS}|}{\bar{p}_c} \times 100$.
Upon evaluating all the required parameters, we first present the comparison of the perturbed pressure term (
$\hat {p}$) in figure 4 and total pressure (
$\bar {p}$) in figure 5 against DSMC data. In figures 4(a) and 5(a), 
$\hat {a}$ is 
$0.005898$, while in figures 4(b) and 5(b), 
$\hat {a}$ is 
$0.09176$. For figure 4, the value of 
$\hat {p}_{c}$ is 
$0.0$, while 
$\bar {p}_c$ is 
$4.6056$ and 
$2.0007$, respectively, for figures 5(a) and 5(b). These conditions have been obtained from DSMC data. Figures 4–5 reveal that pressure also varies along the normal direction of the microchannel and has a maximum value at the boundary and a minimum value at the centre, which is consistent with DSMC data. These results also emphasize that this pressure variation at a high Knudsen number is beyond the reach of the N-S equations (Arkilic, Schmidt & Breuer Reference Arkilic, Schmidt and Breuer1997; Zheng et al. Reference Zheng, Garcia and Alder2002; Rath et al. Reference Rath, Yadav and Agrawal2021).
Figure 4. Comparison of perturbed pressure 
$\hat {p}$ against the DSMC results and N-S equations solution across the microchannel at (a) 
$x/L = 0.1$ and (b) 
$x/L = 0.8$. The pressure has been normalized by the exit centreline pressure (
$p_{out}$). Therefore, the value of pressure is not close to unity at the centreline.
Figure 5. Comparison of total pressure 
$\bar {p}$ against the DSMC results and N-S equations solution across the microchannel at (a) 
$x/L = 0.1$ and (b) 
$x/L = 0.8$. The pressure has been normalized by the exit centreline pressure (
$p_{out}$). Therefore, the value of pressure is not close to unity at the centreline.
Figure 6 demonstrates the agreement between the present solution and DSMC data for the streamwise velocity field. Initial conditions employed for the calculations are as follows: at 
$x/L = 0.1$; 
$\bar {u}|_{\bar {y}=0}=0.3086,\ {\partial ^2 \bar {u}}/{\partial \bar {y}^2}|_{\bar {y}=0}=-1.626$, 
${\partial ^3 \bar {u}}/{\partial \bar {y}^3}|_{\bar {y}=0}=0.0005694$ and at 
$x/L = 0.8$; 
$\bar {u}|_{\bar {y}=0}=0.6549,\ {\partial ^2 \bar {u}}/{\partial \bar {y}^2}|_{\bar {y}=0}=-2.32$, 
${\partial ^3 \bar {u}}/{\partial \bar {y}^3}|_{\bar {y}=0}= -0.005458$. These conditions have been obtained from DSMC data. Notably, the solution incorporating third-order terms exhibits superior performance compared with the first-order and second-order solutions. This improvement is evident as both the first-order and second-order solutions tend to overestimate the velocity in the bulk region under the same condition, which becomes significant at 
$Kn|_{\bar {x}=0.8} = 0.2578$, as shown in figure 6(b).
Figure 6. Comparison of streamwise velocity (
$\bar {u}$) obtained using the EOBurnett and SOBurnett equations against the DSMC results and N-S equations solution across the microchannel at 
$x/L = 0.1$ and 
$x/L = 0.8$.
To elucidate how the present solution varies when considering the boundary value problem, in contrast to the aforementioned initial value problem, we have illustrated the variation of pressure and velocity in figure 7. The pressure variation remains the same as observed in figure 7(a). From figure 7(b), we note that the qualitative behaviour of the velocity field remains consistent, though there are quantitative differences compared with when the problem is addressed as an initial value problem. These quantitative differences are evident in the solutions derived from the N-S and EOBurnett equations since the integration constants have been obtained at different locations of the microchannel. However, as anticipated, the solutions from the SOBurnett equations remain indistinguishable.
Figure 7. Comparison of (a) pressure and (b) streamwise velocity (
$\bar {u}$) obtained using the EOBurnett and SOBurnett equations against the DSMC results and N-S equations solution across the microchannel at 
$x/L = 0.8$ when the present problem has been treated as a boundary value problem.
6. Discussion
The objectives of this section are: first, to present new insights from our study; second, to demonstrate the consistency of our results with prior analytical and simulation studies; third, to comment on the novelty of the work; and finally, to bring out the complexity of higher-order equations. We begin by pointing out the importance of our study regarding the issue with the conventional Burnett and super-Burnett equations. This includes comparing our proposed equations with these and the remaining existing equations and highlighting any differences or improvements. We discuss the stability analysis and compare it with the reported result of the Burnett and super-Burnett equations. Subsequently, we demonstrate that our present analytical solution can easily reproduce the solution of the N-S equations, providing additional analytical validation of our obtained solution. This is followed by a discussion of the results presented in § 5 and the novelty of the work. Finally, we examine the usefulness and complexity of the higher-order transport equations.
6.1. Comparison with existing equations
The linear forms of the Wood, conventional and super-Burnett equations are unstable when the wavenumber is larger than a critical value (Bobylev Reference Bobylev1982, Reference Bobylev2006; Welder et al. Reference Welder, Chapman and Maccormack1993; Agrawal et al. Reference Agrawal, Kushwaha and Jadhav2020), while the proposed equations are unconditionally stable at Burnett as well as a super-Burnett level, as shown in § 4. Additionally, in comparison with one-dimensional flow, the value of the critical Knudsen number for the onset of instability becomes smaller in the 2-D flow for the original and conventional Burnett equations (Bao & Lin Reference Bao and Lin2005), whereas there is no such limitation in the proposed equations. Due to the stability issue linked with these equations, the augmented Burnett equations were derived by Zhong, MacCormack & Chapman (Reference Zhong, MacCormack and Chapman1993) by adding some super-Burnett-order terms in an ad hoc manner to make the conventional Burnett equations stable. In contrast, the proposed equations are unconditionally stable with and without any third-order terms in (3.6) and (3.7) because these equations have been derived using first principles based on the newly proposed Onsager-principle-consistent approach.
It is commonly known that stability concerns arise when dealing with the conventional Burnett equations due to the additional linear terms. To address this issue, Zhao et al. (Reference Zhao, Chen and Agarwal2014) formulated simplified conventional Burnett equations to compute rarefied hypersonic flows and eliminated two linear terms from the stress constitutive relationship. The two terms were discarded based on an order of magnitude analysis. From these two terms, one term, a second-order derivative of pressure that is responsible for the analytical solution of pressure in the present solution, was also discarded in that analysis. However, we have observed that introducing this additional linear term does not disrupt the linear stability of the proposed second- and third-order equations, which is a significant achievement of the present study. From figure 1, it is clear that the proposed set of equations satisfies the stability condition and is unconditionally stable for any wavelength disturbances, even for 2-D flow.
Further, the constitutive relationships for the proposed equations, similar to the conventional Burnett equations, depend on the interaction coefficient, which varies according to the type of molecules involved. The results presented in § 5 demonstrate that the specific type of molecules present in the flow influences the solution of the microchannel flow problem. This observation holds true for both the conventional and proposed equations. This implies that the proposed equations, similar to conventional Burnett equations, apply to monatomic gases regardless of their molecular composition. Moreover, two different relaxation times have been used in the present work; hence the proposed set of equations provides the correct Prandtl number for monatomic gases (Singh & Agrawal Reference Singh and Agrawal2016), unlike the BGK–Burnett equations.
Unlike the OBurnett equations from Singh et al. (Reference Singh, Jadhav and Agrawal2017), which utilized only the first three terms of the distribution function, our equations integrate the entire distribution function (2.28). Furthermore, our methodology employs the N-S equations to replace the material derivatives within the distribution function, introducing third-order terms in the Knudsen number. This makes our proposed equations third-order accurate in Knudsen number, a step above the second-order accuracy of the OBurnett equations. The present equations form a super-set of the OBurnett equations, containing additional linear and nonlinear terms in the relationships of the stress tensor and heat flux vector (3.6)–(3.7). Moreover, neglecting nonlinear terms will result in a simplified set of equations (4.2)–(4.5) having additional higher-order linear terms, while this is not true in the case of the OBurnett equations (Singh et al. Reference Singh, Jadhav and Agrawal2017). Consequently, the linear stability of our proposed equations also differs from that of the OBurnett equations.
6.2. Comparison of analytical solution
The problem solved using the proposed set of equations in § 5 has been also solved by Arkilic et al. (Reference Arkilic, Schmidt and Breuer1997) and Rath et al. (Reference Rath, Singh and Agrawal2018, Reference Rath, Yadav and Agrawal2021) using the N-S equations and conventional Burnett equations, respectively. Arkilic et al. (Reference Arkilic, Schmidt and Breuer1997) reported that the pressure field only depends upon the axial direction, while Rath et al. (Reference Rath, Singh and Agrawal2018, Reference Rath, Yadav and Agrawal2021) presented a more generalized form of a solution of the cross-stream momentum equation analytically in which the pressure field depends upon both the normal and axial directions. In the same work, Rath et al. (Reference Rath, Singh and Agrawal2018, Reference Rath, Yadav and Agrawal2021) presented the numerical solution of the streamwise velocity field due to the complexity in analytically solving the streamwise momentum equation. However, the analytical solutions provide a clearer insight into rarefaction effects. Their generalized nature makes them broader in scope and more time efficient than DSMC simulations. Therefore, in § 5, we successfully presented the analytical solution of the streamwise velocity field for both EOBurnett and SOBurnett equations using the perturbation method. The analytical solution achieved in this study represents a notable advancement as no prior works exist, to the best of the authors’ knowledge, regarding an analytical solution for the streamwise velocity field with any variant of the Burnett equations for the specific problem under investigation. This represents the complexity of the considered problem with Burnett and super-Burnett-type equations and the significance of the present solution obtained in this work. In this sense, the present study is an important progress in the search for accurate higher-order transport equations.
To compare (5.8) and (B1) further with the solution of the N-S equations (Arkilic et al. Reference Arkilic, Schmidt and Breuer1997), we substitute pressure initial conditions to be zero in (5.8) and (B1). As a result, we obtain a simpler form of the solution for pressure and streamwise velocity as
$$\begin{gather} \bar{p} = \bar{p}_{NS} \end{gather}$$
$$\begin{gather}\bar{u} = \bar{y} {V_{2}} + {V_{1}} + \frac{Re \epsilon \bar{y}^{2} }{2 Ma^{2} \gamma} \frac{{\rm d} {\bar{p}_{NS}}}{{\rm d} \bar{x}}. \end{gather}$$It is exciting to note that (6.1) is the same as the solution presented in Arkilic et al. (Reference Arkilic, Schmidt and Breuer1997). This similarity strengthens our confidence in and validates the proposed equations at the level of the N-S equations.
Upon comparison of (5.1)–(5.2) in the present work with the corresponding equations in Rath et al. (Reference Rath, Singh and Agrawal2018), we have noted two major differences in the streamwise momentum equation; first, the second term in the streamwise momentum equation is not present in the corresponding (5.1) in the present work and, second, the presence of an additional fourth-order derivative term at the super-Burnett level in (5.1). Conversely, the cross-streamwise momentum equation remains unchanged. Therefore, the solution of pressure will be the same, while the solution of velocity will be different under the same assumptions at both Burnett and super-Burnett levels.
Following Rath et al. (Reference Rath, Singh and Agrawal2018), the perturbation method has been used to obtain the analytical solution for pressure from the cross-streamwise momentum in the present work, while the same equation has been solved in the exact form in Rath et al. (Reference Rath, Yadav and Agrawal2021). The comparison of the pressure profile presented in Rath et al. (Reference Rath, Yadav and Agrawal2021) reveals that the results obtained in these two studies (Rath et al. Reference Rath, Singh and Agrawal2018, Reference Rath, Yadav and Agrawal2021) are almost the same at 
$Kn =0.1$. This implies that the pressure solution obtained in the present work will also provide the same result. Conversely, the solution of the velocity is not available in the literature for Burnett- and super-Burnett-order equations. That is why any comparison with the solution of existing equations is not possible.
The obtained result has been compared with DSMC data to prove the above claim. In this context, figure 5 presents a very important outcome of the present work since this result shows the limitation of the solution of N-S equations, which fail to capture the variation of pressure in the normal direction of the microchannel. This improvement in the present solution near the wall can be attributed to the presence of these exponential terms, which capture the influence of the Knudsen layer and its effect on pressure and velocity. The present result of pressure matches well with the DSMC data even in the transition regime (
$Kn>0.1$), where most of the available equations fail to predict the pressure profile even in the qualitative sense. In figure 4, we compare the perturbed pressure term variation against the DSMC data to emphasize the explicit effect of perturbed pressure term variation over the solution while figure 5 presents the variation of total pressure against the DSMC results. These two results show that the curvature of the pressure is the same in both figures when the Knudsen number is the same. Moreover, the perturbed and total pressures also have the same qualitative variations across the microchannel cross-section. This pressure variation in the normal direction is critically important since it is the contribution of additional correction terms added to the constitutive relations of the N-S equations.
From figure 6, it is evident that the discrepancies between the results of DSMC and (B1) are less than those with (5.8) and (6.1b) in the Knudsen layer near the wall even though we have only provided the initial conditions at the centre of the channel. This discrepancy in result becomes significant at a higher Knudsen number, as shown in figure 6(b), which implies that the velocity results obtained from the N-S equations show the largest discrepancy. On the other hand, the SOBurnett equations exhibit the smallest deviation from the DSMC data, which proves the supremacy of this set of equations’ solution in the transition regime. Hence, the obtained solutions are more generalized and supposed to work at a wide range of Knudsen number even in the transition regime of flow.
As a result, the present result is more consistent with the DSMC result and more accurate than the N-S equations’ solution across the microchannel, as shown in figures 4–6. However, some discrepancies near the wall can be observed, which might be due to the constant viscosity or isothermal assumption. Perhaps for the first time, we are presenting agreement of any theory with DSMC data in the transition regime.
6.3. Significance, novelty and benefit of the present work
In the present work, we introduced a novel single-particle distribution function derived through an iterative refinement technique in which the core is grounded in non-equilibrium thermodynamics. Therefore, it uniquely combines the Chapman–Enskog expansion and Onsager's reciprocity principle for the first time, underscoring its novelty and significance in the field. We then introduced the EOBurnett and SOBurnett equations (3.6)–(3.7), derived from the novel distribution function, consistent with Onsager's reciprocity principle. These equations contain second- and third-order linear terms, distinguishing them from the OBurnett equations (Singh et al. Reference Singh, Jadhav and Agrawal2017), while maintaining linear stability (figures 1–2). The need for having a linear term of higher order in Knudsen number motivated Zhong et al. (Reference Zhong, MacCormack and Chapman1993) to propose the augmented Burnett equations. Towards this, they added a linear term of super-Burnett order in the conventional Burnett equations in an ad hoc manner, whereas we base our derivation on sound physical principles. This distinction underscores the significance of our work, as it avoids the need for introducing extraneous linear terms that cannot be derived from first principles.
We further present an analytical solution for pressure and velocity for the plane Poiseuille flow problem. The solution is applicable in the transition flow regime and is consistent with the DSMC data. Recently, only an analytical solution of pressure has been reported using the conventional Burnett equations by our group (Rath et al. Reference Rath, Singh and Agrawal2018, Reference Rath, Yadav and Agrawal2021). Owing to the well-known limitations and complexity of existing higher-order transport equations such as the conventional, augmented and super-Burnett equations, a complete solution including pressure and velocity under isothermal conditions could not have been derived in the past. In this sense, the present study is an important advancement in the case of higher-order transport equations. It is noteworthy that the analytical solution of the pressure adeptly captures the normal pressure variation within the microchannel, attributed to the presence of an additional term responsible for modelling the Knudsen layer. Moreover, the velocity solution also has several additional terms that are responsible for capturing the Knudsen layer, which becomes important in the transition flow regime. Description of the Knudsen layer effect, a manifestation of rarefaction, is, however, beyond the capabilities of the N-S equations (Arkilic et al. Reference Arkilic, Schmidt and Breuer1997; Zheng et al. Reference Zheng, Garcia and Alder2002; Rath et al. Reference Rath, Yadav and Agrawal2021). Therefore, deriving these analytical solutions is an important step in the field of higher-order transport equations. Furthermore, the analytical solutions provide a deeper understanding of the rarefaction effects. These solutions are not only more comprehensive but also more time efficient as compared with DSMC simulations. Hence, the current analytical solution readily enables the comprehension of the variation of conserved field quantities, non-conserved field quantities and the influence of the Knudsen layer on the solution.
6.4. Usefulness and complexity of higher-order transport equations
Whereas we here have been able to provide an analytical solution of the EOBurnett and SOBurnett equations, working with higher-order transport equations remains a challenge. The issue of boundary conditions remains a longstanding challenge even at the Burnett level, let alone at the super-Burnett level. Nevertheless, clues may be taken from the work done for moment-type higher-order transport equations, such as the R13 moment equations (Torrilhon & Struchtrup Reference Torrilhon and Struchtrup2008), to evaluate the boundary conditions for Burnett-type equations as well. Similarly, it may be possible to explore data-driven techniques to estimate the models for the boundaries under such conditions. However, the complexity of the equations at the Burnett and super-Burnett levels renders the particularly hard question of obtaining physically accurate descriptions of boundary conditions as an open research problem.
Notwithstanding the problem of boundary conditions, it is important to note that such equations find utility and prove indispensable in capturing non-equilibrium effects for problems such as shock waves (Jadhav & Agrawal Reference Jadhav and Agrawal2020b), standing shear waves (Lockerby & Reese Reference Lockerby and Reese2008) and Grad's second problem (Jadhav & Agrawal Reference Jadhav and Agrawal2020a), where explicit boundary conditions may not be required due to the absence of physical boundaries. With advancements in computational methods and increasing computational power, it has become feasible to tackle the complexities associated with these equations. For instance, adaptive mesh refinement techniques can be used to focus computational efforts in regions where the higher-order terms of the super-Burnett-level equations are significant. Nonetheless, the theoretical framework presented in this study will attain its full completeness upon the derivation of boundary conditions. This accomplishment will extend the applicability of the equations to intricate higher-dimensional flow scenarios. Consequently, based on the proposed equations, a comprehensive computational fluid dynamics solver can be developed to tackle rarefied gas flow problems effectively.
7. Conclusions
In this work, we modify the Onsager-consistent distribution function by adding several additive-invariant-property-consistent correction terms of Burnett and super-Burnett order by employing the iterative refinement technique. Further, we employ this modified density function to obtain the constitutive relationships for the stress tensor and heat flux vector having additional linear and nonlinear terms to close the proposed set of second- and third-order transport equations. These sets of second- (EOBurnett) and third-order (SOBurnett) equations have been demonstrated to be unconditionally stable for 2-D flow for any disturbances in space. We then solve a long microchannel pressure-driven Poiseuille flow problem by assuming small Mach and Reynolds numbers to validate the proposed equations. This resulted in the analytical expression of the velocity and pressure as functions of both the normal and axial directions for both second- and third-order equations, which has not been obtained using any Burnett- and super-Burnett-level equations in the past. We demonstrate that the present solutions quantitatively agree with existing DSMC results in the transition regime. The present solution for pressure successfully captures the normal variation of pressure with enough accuracy due to the added additional term in the constitutive relationships. Due to these terms, the proposed equations also accurately capture the velocity profile near the wall. Therefore, the proposed equations successfully capture the Knudsen layer phenomenon in the transition regime of rarefied flow.
The obtained results demonstrate that the proposed equations are more precise than the N-S equations and are expected to provide a better description of the flow physics at a large Knudsen number. Unlike the N-S equations, the proposed equations can capture the non-Newtonian stress. The proposed equation can also be used on a very fine grid while applying the numerical method. Although the results demonstrated in the present work suggest that the model would perform robustly within the transition flow regime, a comprehensive investigation to determine its precise limits of accuracy in Knudsen number needs to be undertaken.
In conclusion, the scope of the present work encompasses a wide range of Knudsen number. We discuss how this work can be employed in the context of issues related to the conventional, augmented and super-Burnett equations. Moreover, the approach followed in the present work can also be used to solve the existing higher-order set of equations and test their validity in the transition regime of gaseous rarefied flow, which has been hindered due to the complex nature of the equations. These findings have important implications for accurately predicting fluid flows in microchannels and other rarefied gas systems for their design.
Acknowledgements
We thank Dr A. Gavasane for providing the DSMC code and data for the pressure-driven plane Poiseuille flow and Dr R.S. Jadhav for the helpful and stimulating discussions.
Funding
This research received no specific grant from any funding agency, commercial or not-for-profit sectors.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The data that support the findings of this study are available from the corresponding author upon reasonable request.
Appendix A. Constitutive relations and governing equations
The explicit expressions of the 2-D stress tensor and heat flux vector of third-order in Knudsen number are required to perform stability analysis and derive the analytical solution in the present work. Therefore, these expressions obtained from (3.6)–(3.7) are expressed as follows:
\begin{align} \sigma_{11} &= \frac{\mu}{
27 R^{2} T^{2} \rho^{3} } (18 R^{2} \mu [- 2 T_{x x} + T_{y
y}] T \rho^{2} + 18 R^{2} \mu [2 \rho_{x}^{2} -
\rho_{y}^{2}] T^{2} \nonumber\\ &\quad + 18 R^{2} \mu [- 2
T_{x} \rho_{x} + T_{y} \rho_{y}] T \rho + 18 R^{2} [- 2
u_{x} + v_{y}] T^{2} \rho^{3} + 6 R \mu^{2} [8 u_{x x x} +
5 u_{y x y}] T \rho \nonumber\\ &\quad + 45 R \mu^{2} [- 2
T_{x x} u_{x} + T_{x x} v_{y} - 2 T_{y y} u_{x} + T_{y y}
v_{y}] \rho + 6 R \mu^{2} [- 8 \rho_{x} u_{x x} - 6
\rho_{x} u_{y y} \nonumber\\ &\quad -2\rho_{x} + 3 \rho_{y}
v_{x x} + 4 \rho_{y} v_{y y} + \rho_{y}] T + 6 R \mu \big[- 4
u_{x}^{2} \varphi + 8 u_{x}^{2} - 2 u_{x} v_{y} \varphi + 7
u_{x} v_{y} \nonumber\\ &\quad - 3 u_{y} v_{x} + 2v_{y}^{2}
\varphi - 4 v_{y}^{2}\big] T \rho^{2} + 4 \mu^{2} [- 8
u_{x}^{3} + 12 u_{x}^{2} v_{y}- 6 u_{x} u_{y}^{2} - 12
u_{x} u_{y} v_{x} - 6 u_{x} v_{x}^{2} \nonumber\\ &\quad -
12 u_{x} v_{y}^{2}+ 3 u_{y}^{2} v_{y} + 6 u_{y} v_{x} v_{y}
+ 3 v_{x}^{2} v_{y} + 4 v_{y}^{3}] \rho)
\end{align}
\begin{align} \sigma_{21} &=
\frac{\mu}{18 R^{2} T^{2} \rho^{3} } (36 R^{2} \mu \rho_{x}
\rho_{y} T^{2} - 18 R^{2} \mu [T_{x} \rho_{y} + T_{y}
\rho_{x}] T \rho - 18 R^{2}[u_{y} + v_{x}] T^{2} \rho^{3}
\nonumber\\ &\quad + 6 R \mu^{2} [3 u_{y y y} + 5 v_{y x
y}] T \rho - 45 R \mu^{2} [T_{x x} u_{y} + T_{x x} v_{x} +
T_{y y} u_{y} + T_{y y} v_{x}] \rho \nonumber\\ &\quad - 6
R \mu^{2} [3 \rho_{x} v_{x x} + 4 \rho_{x} v_{y y} +
\rho_{x} + 4 \rho_{y} u_{x x} + 3 \rho_{y} u_{y y} +
\rho_{y}] T \nonumber\\ &\quad + 12 R \mu [- u_{x} u_{y}
\varphi + 2 u_{x} u_{y} - u_{x} v_{x} \varphi + 2 u_{x}
v_{x} - u_{y} v_{y} \varphi \nonumber\\ &\quad + 2 u_{y}
v_{y} - v_{x} v_{y}\varphi + 2 v_{x} v_{y}] T \rho^{2} + 4
\mu^{2} [- 4 u_{x}^{2} u_{y} - 4 u_{x}^{2} v_{x}
\nonumber\\ &\quad + 4 u_{x} u_{y} v_{y} +4 u_{x} v_{x}
v_{y} - 3 u_{y}^{3} - 9 u_{y}^{2} v_{x} - 9 u_{y} v_{x}^{2}
- 4 u_{y} v_{y}^{2} - 3 v_{x}^{3} - 4 v_{x} v_{y}^{2}]
\rho)\end{align}
\begin{align} \sigma_{22} &= \frac{\mu}{27 R^{2} T^{2} \rho^{3}} (18 R^{2} \mu [T_{x x} - 2 T_{y y}] T \rho^{2} + 18 R^{2} \mu [- \rho_{x}^{2} + 2 \rho_{y}^{2}] T^{2} \nonumber\\ &\quad + 18 R^{2} \mu [T_{x} \rho_{x} - 2 T_{y} \rho_{y}] T \rho + 18 R^{2} [u_{x} - 2 v_{y}] T^{2} \rho^{3} + 6 R \mu^{2} [- 4 u_{x x x} - u_{y x y}] T \rho \nonumber\\ &\quad + 45 R \mu^{2} [T_{x x} u_{x} - 2 T_{x x} v_{y} + T_{y y} u_{x} - 2 T_{y y} v_{y}]\rho \nonumber\\ &\quad + 6 R \mu^{2} [4 \rho_{x} u_{x x} + 3 \rho_{x} u_{y y} + \rho_{x} - 6 \rho_{y} v_{x x} - 8 \rho_{y} v_{y y} - 2 \rho_{y}] T \nonumber\\ &\quad + 6 R \mu [2 u_{x}^{2}\varphi - 4 u_{x}^{2} - 2 u_{x} v_{y} \varphi + 7 u_{x} v_{y} - 3 u_{y} v_{x} - 4 v_{y}^{2} \varphi + 8 v_{y}^{2}] T \rho^{2} \nonumber\\ &\quad + 4 \mu^{2} [4 u_{x}^{3} - 12 u_{x}^{2} v_{y} + 3 u_{x} u_{y}^{2} + 6 u_{x} u_{y} v_{x} + 3 u_{x} v_{x}^{2} + 12 u_{x} v_{y}^{2} - 6 u_{y}^{2} v_{y} - 12 u_{y} v_{x} v_{y} \nonumber\\ &\quad - 6 v_{x}^{2} v_{y} - 8 v_{y}^{3}]\rho) \end{align}
\begin{align} q_x &= \frac{\mu}{48 R T^{2} \rho^{3}} (- 180 R^{2} T_{x} T^{2} \rho^{3} - 1350 R T_{x} \mu^{2} (T_{x x} + T_{y y}) \rho - 675 R \mu^{2} \rho_{x} (T_{x x} + T_{y y}) T \nonumber\\ &\quad - 180 R \mu (u_{x x} + 1) T^{2} \rho^{2} + 16 R \mu (- 4 \rho_{x} u_{x} + 2 \rho_{x} v_{y} - 3 \rho_{y} u_{y} - 3 \rho_{y} v_{x}) T^{2} \rho \nonumber\\ &\quad + 2 R \mu (- 90 T_{x} u_{x} \varphi + 287 T_{x} u_{x} - 90 T_{x} v_{y} \varphi + 59 T_{x} v_{y} + 249 T_{y} u_{y} - 21 T_{y} v_{x}) T \rho^{2} \nonumber\\ &\quad + 120 T_{x} \mu^{2} (- 4 u_{x}^{2} + 4 u_{x} v_{y} - 3 u_{y}^{2} - 6 u_{y} v_{x} - 3 v_{x}^{2} - 4 v_{y}^{2}) \rho \nonumber\\ &\quad + 60 \mu^{2} \rho_{x} (- 4 u_{x}^{2} + 4 u_{x} v_{y} - 3 u_{y}^{2} - 6 u_{y} v_{x} - 3 v_{x}^{2} - 4 v_{y}^{2}) T \nonumber\\ &\quad + 120 \mu^{2} (4 u_{x} u_{x x} - 2 u_{x} - 2 u_{x x} v_{y} + 3 u_{y} v_{x x} + 3 u_{y} + 3 v_{x} v_{x x} + 3 v_{x} + 4 v_{y}) T \rho) \end{align}
\begin{align} q_y &= \frac{\mu}{48 R
T^{2} \rho^{3}} (- 180 R^{2} T_{y} T^{2} \rho^{3} - 1350 R
T_{y} \mu^{2} (T_{x x} + T_{y y}) \rho - 675 R \mu^{2}
\rho_{y} (T_{x x} + T_{y y}) T \nonumber\\ &\quad - 180 R
\mu (v_{y y} + 1) T^{2} \rho^{2} + 16 R \mu (- 3 \rho_{x}
u_{y} - 3 \rho_{x} v_{x} + 2 \rho_{y} u_{x} - 4 \rho_{y}
v_{y}) T^{2} \rho \nonumber\\ &\quad + 2 R \mu (- 21 T_{x}
u_{y} + 249 T_{x} v_{x} - 90 T_{y} u_{x} \varphi + 59 T_{y}
u_{x} - 90 T_{y} v_{y} \varphi + 287 T_{y} v_{y}) T
\rho^{2} \nonumber\\ &\quad + 120 T_{y} \mu^{2} (- 4
u_{x}^{2} + 4 u_{x} v_{y} - 3 u_{y}^{2} - 6 u_{y} v_{x} - 3
v_{x}^{2} - 4 v_{y}^{2}) \rho \nonumber\\ &\quad + 60
\mu^{2} \rho_{y} (- 4 u_{x}^{2} + 4 u_{x} v_{y} - 3
u_{y}^{2} - 6 u_{y} v_{x} - 3 v_{x}^{2} - 4 v_{y}^{2}) T
\nonumber\\ &\quad + 120 \mu^{2} (- 2 u_{x} v_{y y} + 4
u_{x} + 3 u_{y} u_{y y} + 3 u_{y} + 3 u_{y y} v_{x} + 3
v_{x} + 4 v_{y} v_{y y} - 2 v_{y}) T \rho),
\end{align}
where 
$()_x = {\partial }/{\partial x}$ and 
$()_y = {\partial }/{\partial y}$ have been used to express the derivative terms in the more compact and simpler form to avoid the complexity due to the presence of a large number of terms in the constitutive relationships and equations. Note that the conventional form of the derivative has been used in the main part of this paper since fewer terms are present in the reduced equations presented in §§ 4 and 5.
Upon substituting the above stress tensor in the momentum, (3.2), the continuity and (stream and cross-stream) momentum balance can be expressed in non-dimensional form as
A.1. Continuity equations
\begin{equation} \epsilon (\bar{\rho}\bar{u})_{\bar{x}}+ (\bar{\rho}\bar{v})_{\bar{y}} = 0 .\end{equation}A.2. Stream momentum equation
A.3. Cross-stream momentum equation
where 
$\eta _x =\eta _y = {1}/({135 Ma^{2} Re^{2} \bar {p}^{4}})$.
Appendix B. Solution for velocity field
Solving the streamwise momentum equation yields the velocity field
\begin{align} \bar{u}&=\left\{\left[6
A^{4} Ma^{2} \gamma (A^{2} Ma^{2} \gamma - Re^{2}) (\bar{y}
{C_{2}} + {C_{1}} + {C_{4}} \exp({B \bar{y}
{\bar{p}_{NS}}})) (A^{4} Ma^{4} \gamma^{2}
\vphantom{\frac{{\rm d} {\bar{p}_{NS}}}{{\rm d}
\bar{x}}}\right.\right. \nonumber\\ &\quad - 2 A^{2} Ma^{2}
Re^{2} \gamma + Re^{4}) {\bar{p}_{NS}^{5}} \exp({\bar{y} (2
A + B) {\bar{p}_{NS}}}) \nonumber\\ &\quad +6 A^{4} Ma^{2}
\gamma (A^{2} Ma^{2} \gamma - Re^{2}) (A^{4} Ma^{4}
\gamma^{2}- 2 A^{2} Ma^{2} Re^{2} \gamma \nonumber\\ &\quad
+ Re^{4}) {C_{3}} {\bar{p}_{NS}^{5}} \exp({2 A \bar{y}
{\bar{p}_{NS}}}) + 3 A^{4} Re \epsilon \bar{y}^{2} (A^{2}
Ma^{2} \gamma - Re^{2}) \left(\frac{{\rm d}
{\bar{p}_{NS}}}{{\rm d} \bar{x}} + \frac{{\rm d}
\hat{p}_{c}}{{\rm d} \bar{x}} \right) \nonumber\\ &\quad
\times (A^{4} Ma^{4} \gamma^{2} - 2 A^{2} Ma^{2} Re^{2}
\gamma + Re^{4}) {\bar{p}_{NS}^{5}} \exp({\bar{y} (2 A + B)
{\bar{p}_{NS}}}) \nonumber\\ &\quad + Re \epsilon(A^{2}
Ma^{2} \gamma - Re^{2}) (4 A^{4} Ma^{4} \gamma^{2} - 7
A^{2} Ma^{2} Re^{2} \gamma \nonumber\\ &\quad - 3 A^{2}
\bar{y}^{2}(A^{4} Ma^{4} \gamma^{2} - 2 A^{2} Ma^{2} Re^{2}
\gamma + Re^{4}) {\bar{p}_{NS}^{2}} \exp({A \bar{y}
{\bar{p}_{NS}}}) + 3 Re^{4} \nonumber\\ &\quad + (4 A^{4}
Ma^{4} \gamma^{2} - 7 A^{2} Ma^{2} Re^{2} \gamma + 3
Re^{4}) \exp({2 A \bar{y} {\bar{p}_{NS}}})) {\bar{p}_{NS}}
\exp({\bar{y} (A + B) {\bar{p}_{NS}}}) \frac{{\rm d}
\hat{a}}{{\rm d} \bar{x}} \nonumber\\ &\quad - Re \epsilon
(18 A^{4} Ma^{4} \gamma^{2} (A^{2} Ma^{2} \gamma - Re^{2})
+ 4 A^{3} Ma^{2} \gamma \bar{y} (A^{4} Ma^{4} \gamma^{2} -
2 A^{2} Ma^{2} Re^{2} \gamma \nonumber\\ &\quad + Re^{4})
{\bar{p}_{NS}} - 28 A^{2} Ma^{2} Re^{2} \gamma (A^{2}
Ma^{2} \gamma - Re^{2}) - 6 A^{2} \bar{y}^{2} (A^{2} Ma^{2}
\gamma - Re^{2}) \nonumber\\ &\quad \times (A^{4} Ma^{4}
\gamma^{2}- 2 A^{2} Ma^{2} Re^{2} \gamma + Re^{4})
{\bar{p}_{NS}^{2}} \exp({A \bar{y} {\bar{p}_{NS}}}) - 3 A
Re^{2} \bar{y} (A^{4} Ma^{4} \gamma^{2} \nonumber\\ &\quad
- 2 A^{2} Ma^{2} Re^{2} \gamma + Re^{4}){\bar{p}_{NS}} + 12
Re^{4} (A^{2} Ma^{2} \gamma - Re^{2}) \nonumber\\ &\quad
+[18 A^{4} Ma^{4} \gamma^{2} (A^{2} Ma^{2} \gamma - Re^{2})
- 4 A^{3} Ma^{2} \gamma \bar{y} (A^{4} Ma^{4} \gamma^{2}
\nonumber\\ &\quad - 2 A^{2} Ma^{2} Re^{2} \gamma + Re^{4})
{\bar{p}_{NS}} - 28 A^{2} Ma^{2} Re^{2} \gamma (A^{2}
Ma^{2} \gamma - Re^{2}) \nonumber\\ &\quad +3 A Re^{2}
\bar{y} (A^{4} Ma^{4} \gamma^{2} - 2 A^{2} Ma^{2} Re^{2}
\gamma + Re^{4}) {\bar{p}_{NS}} \nonumber\\ &\quad \left.+
12 Re^{4}(A^{2} Ma^{2} \gamma - Re^{2})] \exp({2 A \bar{y}
{\bar{p}_{NS}}})) \hat{a} \exp({\bar{y} (A + B)
{\bar{p}_{NS}}})\frac{{\rm d} {\bar{p}_{NS}}}{{\rm d}
\bar{x}}\right] \nonumber\\ &\quad
\left.\vphantom{\frac{{\rm d} {\bar{p}_{NS}}}{{\rm d}
\bar{x}}} \times\exp({- \bar{y} (2 A + B)
{\bar{p}_{NS}}})\right\} \frac{1}{\zeta} ,
\end{align}
where
\begin{equation} \left.\begin{gathered} \zeta=[6 A^{4} Ma^{2} \gamma (A^{2} Ma^{2} \gamma - Re^{2}) (A^{4} Ma^{4} \gamma^{2} - 2 A^{2} Ma^{2} Re^{2} \gamma + Re^{4}) {\bar{p}_{NS}^{5}}] \\ B = Re \sqrt{3/5}/Ma, \end{gathered}\right\} \end{equation}
and 
$C_3$–
$C_4$ are additional integration constants appearing due to the presence of the fourth-order derivative term in (5.4).
Open access
