## 2 Methodology Details  
### Governing Equations and Numerical Methods  
LES was conducted to solve the flow process by filtering the turbulent motion and decomposing the turbulent flow into large scale pulsations and small-scale pulsations. The large-scale pulsations are solved by direct numerical simulation (DNS) and the mass, momentum, and energy transport of the small-scale pulsations on the large-scale motion are modeled to close the equations, where the model is called a subgrid model. The Navier-Stokes equations, filtered by Deardroff box filtering [30], are as follows:  
\[\frac{\partial\overline{p}}{\partial t}+\frac{\partial}{\partial x_{i}}( \overline{p}\widetilde{u}_{i})=0 \tag{1}\]  
\[\frac{\partial}{\partial t}(\overline{p}\widetilde{u}_{i})+\frac{\partial}{ \partial x_{i}}\Big{(}\overline{p}\widetilde{u}_{i}\widetilde{u}_{j}+ \overline{p}\delta_{ij}-\widetilde{u}_{ij}+\tau_{ij}^{ggs}\Big{)}=0 \tag{2}\]  
\[\frac{\partial\overline{p}\widetilde{E}}{\partial t}+\frac{\partial}{\partial x _{i}}\Big{[}\Big{(}\overline{p}\widetilde{E}+\overline{p}\Big{)}\widetilde{u }_{i}+\overline{q}_{i}-\widetilde{u}_{j}\widetilde{u}_{ij}+H_{i}^{ggs}+\sigma _{i}^{ggs}\Big{]}=0 \tag{3}\]  
In the above equations, the mark \({}^{\prime}\)\({}^{\prime}\)\({}^{\prime}\)\({}^{\prime}\) denotes spatial filtering and \({}^{\prime}\)\({}^{\prime}\)\({}^{\prime}\) denotes Favre filtering. Superscript'sgs' is variances associated with subgrid model. Here, \(\tau_{ij}\) is the viscous stress tensor, \(q_{i}\) is heat flux vector, and \(\delta_{ij}\) is unit tensor. \(E\) is energy, \(H_{i}\) is the enthalpy flux vector, and \(\sigma_{i}\) is the energy flux vector.  
The kinetic-energy transport sub-grid model is adopted to measure the influence of small-scale transport on large-scale pulsation. The sub-grid viscosity coefficient is defined as  
\[\mu_{t}=C_{k}\rho k_{sgs}^{1/2}\Delta f \tag{4}\]  
The subgrid viscous stress tensor meets the equation as following:  
\[\tau_{ij}-\frac{2}{3}\rho k_{sgs}\delta_{ij}=-2C_{k}\rho k_{sgs}^{1/2}\Delta f \overline{S}_{ij} \tag{5}\]  
In this equation,  
\[k_{sgs}=\frac{1}{2}\Big{(}\overline{u_{K}^{2}}-\overline{u}_{k}^{2}\Big{)} \tag{6}\]  
\[\Delta f\equiv V^{1/3} \tag{7}\]  
\(V\) is the volume of the grid. Moreover, subgrid kinetic energy can be obtained by the following equation:  
\[\rho\bigg{(}\frac{\partial\overline{k}_{\text{sgs}}}{\partial t}+\frac{ \partial\overline{u}_{j}\overline{k_{\text{sgs}}}}{\partial x_{j}}\bigg{)}= -\tau_{ij}\frac{\partial\overline{u}_{i}}{\partial x_{j}}-C_{\text{c}}\rho \frac{k_{sgs}^{3/2}}{\Delta f}+\frac{\partial}{\partial x_{j}}\bigg{(}\frac{ \mu_{t}}{\sigma_{k}}\frac{\partial k_{\text{sgs}}}{\partial x_{j}}\bigg{)} \tag{8}\]It should be noted that \(C_{k}\), \(C_{\epsilon}\) and \(\sigma_{k}\) are constants in above formula. In present work, the implicit method based on density is applied to numerical simulation. The filter equation is discretized by finite volume method and a second-order upwind scheme is utilized to calculate the convection terms. Besides, a second-order central difference scheme is applied to calculating the viscous term, and a second-order implicit scheme is employed to discretize the time. Such numerical schemes [17] are widely used in previous studies on mixing layers.