## 2 High-Order Accurate Numerical Method  
### Governing Equations  
The three-dimensional compressible Euler equations governing inviscid flows can be expressed in a dimensionless conservative form [20] as follows:  
\[\frac{\partial Q}{\partial t}+\
abla\cdot\mathbf{F}=0, \tag{1}\]where the unknown vector \(\mathbf{Q}\) of conserved variables and the Cartesian components of the flux function \(\mathbf{F}=\left(\mathbf{f}_{1},\mathbf{f}_{2},\mathbf{f}_{3}\right)\) are given by  
\[\mathbf{Q}= \begin{pmatrix}\rho\\ pu\\ \rho v\\ \rho w\\ \rho E\end{pmatrix} ,\mathbf{f}_{1}= \begin{pmatrix}\rho u\\ \rho u^{2}+p\\ \rho uv\\ u(\rho E+p)\end{pmatrix} ,\mathbf{f}_{2}= \begin{pmatrix}\rho v\\ \rho uv\\ \rho v^{2}+p\\ \rho vw\\ v(\rho E+p)\end{pmatrix} ,\mathbf{f}_{3}= \begin{pmatrix}\rho w\\ \rho uv\\ \rho w\\ \rho w^{2}+p\\ w(\rho E+p)\end{pmatrix} , \tag{2}\]  
where \(\rho\), \(p\), and \(E\) are the density, pressure, and total energy per unit mass, respectively, and \(u\), \(v\), and \(w\) are the Cartesian components of the velocity vector \(\mathbf{V}\). This system can be closed up by adding the following equation of state:  
\[\rho E=\frac{p}{\Gamma-1}+\frac{1}{2}\rho\Big{(}u^{2}+v^{2}+w^{2}\Big{)}, \tag{3}\]  
which is valid for perfect gas, \(\Gamma=\frac{C_{p}}{C_{p}}\) is the ratio of specific heats and taken as \(\Gamma=1.4\) for air. The Equation (1) in conservative variables are reported to be ill-conditioned for the flows at low Mach numbers [21]. In order to single out the propagation of acoustic waves [22], the governing equations in primitive variables \(\mathbf{q}=\left(p,u,v,w,T\right)^{T}\) are considered here. By transforming the dependent variables in Equation (1) from conserved quantities to such primitive variables, Equation (1) can be written as  
\[\mathbf{P}\frac{\partial\mathbf{q}}{\partial t}+\
abla\cdot\mathbf{F}=0, \tag{4}\]  
where \(\mathbf{P}\) is the transformation matrix  
\[\mathbf{P}=\frac{\partial Q}{\partial q}= \begin{bmatrix}\rho_{p}&0&0&0&\rho_{T} \\ \rho_{p}u&\rho&0&0&\rho_{T}u\\ \rho_{p}v&0&\rho&0&\rho_{T}v\\ \rho_{p}w&0&0&\rho&\rho_{T}w\\ \rho_{p}H-1&\rho u&\rho v&\rho w&\rho_{T}H+\rho C_{p}\end{bmatrix}  \tag{5}\]  
here, \(H=E+p/\rho\) is total enthalpy per unit mass of the fluid. For inviscid flows of interest, the derivatives that appear in Equation (5) can be directly determined by the equation of state, \(p=\frac{1}{T}\rho T\), which reads:  
\[\rho_{p}=\frac{\partial\rho}{\partial p}\Big{|}T=const=\frac{\Gamma}{T}= \frac{\Gamma}{c^{2}}\,,\quad\rho_{T}=\frac{\partial\rho}{\partial T}\Big{|}p =const=-\frac{\rho}{T}\,,\]  
where \(c\) is the acoustic velocity. The derivative \(\rho_{p}\) that multiplies the pressure time derivative in the continuity equation controls the speed of propagation of acoustic waves [23] in the system (4). This is consistent with the notion of infinite pressure wave speeds in an incompressible flow due to \(\rho_{p}=0\) for constant density flows. The disparity of wave speeds of the system (4), mainly between acoustic waves and entropy waves, becomes very large for low Mach number flows, and can be ameliorated through adopting techniques such as preconditioning [10; 23]. Following the well-known preconditioning technique of [PERSON] and [PERSON] [23], the matrix \(\mathbf{P}\) in Equation (4) is modified through replacing the derivative \(\rho_{p}\), with one proportional to the inverse of the local velocity squared to control the disparity of the wave speeds of the system (4). The resultant matrix, termed as \"preconditioning matrix \(\mathbf{\Gamma}\)\", can be expressed as follows:  
\[\mathbf{\Gamma}= \begin{bmatrix}\Theta&0&0&0&\rho_{T}\\ \Theta u&\rho&0&0&\rho_{T}u\\ \Theta v&0&\rho&0&\rho_{T}v\\ \Theta w&0&0&\rho&\rho_{T}w\\ \Theta H-1&\rho u&\rho v&\rho w&\rho_{T}H+\rho C_{p}\end{bmatrix} , \tag{6}\]  
where \(\Theta\) is computed by  
\[\Theta=\frac{1}{U_{r}^{2}}-\frac{\rho_{T}}{\rho C_{p}} \tag{7}\]  
Here \(U_{r}\) is a reference velocity defined as follows:  
\[U_{r}=\left\{