## 1 Introduction  
Due to the rapid development of micro-electro-mechanical systems (MEMS) [1; 2] and the shale gas revolution in North America [3], extensive works have been devoted to constructing accurate physical models and efficient numerical methods for non-equilibrium gas flows [4; 5]. In general, the behavior of such non-equilibrium gas flows can be described from either the microscopic or macroscopic points of view. The Boltzmann equation is the fundamental equation for the dynamics of rarefied gases at the microscopic level, which uses a molecular velocity distribution function to describe the system state. Historically, there are two mainly categories of approaches to solve the Boltzmann equation. One is based on the probabilistic approach, well-known examples are the direct simulation Monte Carlo (DSMC) methods developed by [PERSON] and [PERSON] [6; 7]. DSMC method is an excellent approach for solving high-speed rarefied gas flows, however, it becomes time-consuming if the flow is in the slip- and transition-flow regimes. Another one is based on the deterministic approach, well-known examples are the discrete velocity methods (DVMs) [8], which use a finite set of discrete velocity points to approximate the continuous molecular velocity space. In the past two decades, most of the deterministic numerical methods are developed based on the DVM. Many full Boltzmann solvers, especially the Fourierspectral methods [9; 10] and fast-spectral methods [11; 12; 13; 14], provide accurate numerical results, which can serve as reference solutions. However, the structure of the collision operator in the full Boltzmann solver is rather complicated: It is a fivefold integral with three dimensions in velocity space and two dimensions in a unit sphere. Therefore, the Boltzmann equation is usually replaced by simplified kinetic model equations, such as the Bathnagar-Gross-Krook (BGK) [15], ES-BGK [16], and Shakhov [17] models. And most of DVMs are developed for these Boltzmann model equations. Among those DVM solvers, the gas-kinetic unified algorithm was developed and applied to study the hypersonic flows around spacecraft re-entry problems [18; 19]. The asymptotic preserving schemes, especially the unified gas-kinetic scheme (UGKS) [20; 21] and the discrete unified gas-kinetic scheme [22; 23], are proposed for problems involving different flow regimes. The performance of the discrete unified gas kinetic scheme and traditional Godunov DVM, in terms of accuracy and efficiency, has been evaluated with a wide range of Knudsen numbers [24]. Most recently, the unified gas-kinetic wave-particle method was developed under the UGKS framework, in which the evolution of simulation particles and distribution function was guided by evolution of macroscopic variables [25; 26]. Even though great effort has been devoted to reducing the computational cost in solving the Boltzmann equation, the multi-dimensional nature of the VDF and the complicated structure of the nonlinear collision term still make the Boltzmann equation impractical for many applications. Therefore, alternative macroscopic modeling and simulation strategies for rarefied gas flows have been developed [27]. The most classical macroscopic equations for rarefied flows are the Navier-Stokes-Fourier (NSF) equations in association with the velocity-slip and temperature-jump wall boundary conditions. This set of equations has found a wide range of applications in most of today's industrial problems where the Knudsen number \(Kn\), i.e., the ratio of the mean free path \(\lambda_{0}\) of gas molecules to the characteristic flow length \(L_{0}\)[27], is less than \(10^{-1}\). However, additional efforts are required as the Knudsen number increases, and many macroscopic equations are proposed to describe the rarefied gas dynamics beyond the Navier-Stokes level. [PERSON] and [PERSON] proposed a technique via a formal asymptotic expansion of the molecular VDF in powers of the Knudsen number [28; 29].