## 2 Pseudopotential MRT LB Model  
The pseudopotential LB model with MRT operation collision was employed in the present work [21; 22]. The evolution equation can be given as [23]:  
\[f_{\alpha}(\mathbf{x}+\mathbf{e}_{\alpha}\delta_{t},t+\delta_{t})=f_{\alpha}( \mathbf{x},t)-(\mathbf{M}^{-1}\mathbf{\Lambda}\mathbf{M})_{\alpha\beta}\Big{(} f_{\beta}-f_{\beta}^{eq}\Big{)}+\delta_{t}F_{\alpha}^{\prime} \tag{1}\]  
where \(f_{\alpha}\) denotes the density distribution function, \(f_{\alpha}^{eq}\) is its equilibrium distribution, \(t\) is the time, \(\mathbf{x}\) denotes the spatial position, \(\mathbf{e}_{\alpha}\) is the discrete velocity along the \(\alpha\)th direction, \(\delta_{t}\) is the time step, \(F_{\alpha}^{\prime}\) is the component of forcing term \(F^{\prime}\) in the velocity space, \(\mathbf{M}\) is an orthogonal transformation matrix, and \(\mathbf{\Lambda}\) is a diagonal matrix given by (for the D2Q9 lattice) [24]  
\[\mathbf{\Lambda}=diag\Big{(}\tau_{\rho}^{-1},\tau_{e}^{-1},\tau_{\xi}^{-1}, \tau_{j}^{-1},\tau_{q}^{-1},\tau_{j}^{-1},\tau_{q}^{-1},\tau_{v}^{-1}\Big{)} \tag{2}\]  
The \(\tau_{\rho}^{-1}\), \(\tau_{e}^{-1}\), \(\tau_{\xi}^{-1}\), \(\tau_{j}^{-1}\), \(\tau_{q}^{-1}\) and \(\tau_{v}^{-1}\) is related to density, energy, energy squared, components of momentum flux, energy flux components and stress tensor. These relaxation times were chosen as \(\tau_{v}=\tau_{\rho}=\tau_{j}=1.0\), \(\tau_{e}^{-1}=\tau_{\xi}^{-1}=0.8\) and \(\tau_{q}^{-1}=1.1\) in this paper. Through the transformation matrix  
**M**[25], the \(\mathbf{f}_{\mathbf{\alpha}}\) and \(\mathbf{f}_{\mathbf{\alpha}}^{\mathbf{q}}\) can be projected onto the moment space via \(\mathbf{m}=\mathbf{M}\mathbf{f}\) and \(\mathbf{m}^{\mathbf{\alpha}\mathbf{q}}=\mathbf{M}\mathbf{f}^{\mathbf{\alpha}\mathbf{q}}\). Thus, the collision step of Equation (1) can be rewritten as:  
\[\mathbf{m}^{*}=\mathbf{m}-\mathbf{\Lambda}(\mathbf{m}-\mathbf{m}^{\mathbf{\alpha}\mathbf{q} })+\delta_{\mathbf{l}}\left(\mathbf{I}-\frac{\mathbf{\Lambda}}{2}\right)\mathbf{S} \tag{3}\]  
where the \(\mathbf{I}\) is the unite tensor, and the \(\mathbf{S}\) is the forcing term in the moment space with \((\mathbf{I}-0.5\mathbf{\Lambda})\mathbf{S}=\mathbf{M}\mathbf{F}^{\prime}\). In D2Q9 lattice, \(\mathbf{m}^{\mathbf{\alpha}\mathbf{q}}\) can be given by [20]:  
\[\mathbf{m}^{\mathbf{\alpha}\mathbf{q}}=\rho\Big{(}1,-2+3|\mathbf{v}|^{2},1-3|\mathbf{ v}|^{2},v_{x},-v_{x},v_{y},-v_{y},v_{x}^{2}-v_{y}^{2},v_{x}v_{y}\Big{)}^{T} \tag{4}\]  
where the \(\rho=\sum\limits_{\alpha}f_{\alpha}\) is the macroscopic density, \(\mathbf{v}\) is the macroscopic velocity calculated by:  
\[\mathbf{v}=\left(\sum\limits_{\alpha}\mathbf{e}_{\alpha}f_{\alpha}+\frac{ \delta_{\mathbf{l}}}{2}\mathbf{F}\right)/\rho \tag{5}\]  
Here, \(\mathbf{F}=\left(F_{x},F_{y}\right)\) is the total force acting on the system for two-dimensional (2D) space and \(|\mathbf{v}|^{2}=v_{x}^{2}+v_{x}^{2}\). Then the steaming step of the MRT LB equation can be formulated as:  
\[f_{\alpha}(\mathbf{x}+\mathbf{e}_{\alpha}\delta_{\mathbf{l}},\;t+\delta_{\mathbf{l}}) =f_{\alpha}^{*}(\mathbf{x},t) \tag{6}\]  
where \(\mathbf{f}^{*}=\mathbf{M}^{-1}\mathbf{m}^{*}\). In the pseudopotential LB model, the intermolecular interaction force is defined via a pseudopotential \(\psi\), which depends on the local density. It is given as [26]:  
\[\mathbf{F_{m}}=-G\psi(x)\sum\limits_{\alpha}w_{\alpha}\psi(\mathbf{x}+ \mathbf{e}_{\alpha})\mathbf{e}_{\alpha} \tag{7}\]  
where \(G\) is the interaction strength, and \(w_{\alpha}\) is the weight. The weight \(w_{\alpha}=1/3\) for \(|\mathbf{e}_{\alpha}|^{2}=1\) and the weight \(w_{\alpha}=1/12\) for \(|\mathbf{e}_{\alpha}|^{2}=2\) on D2Q9 lattice for the nearest-neighbor interactions. According to [PERSON] and [PERSON] [27], the pseudopotential \(\psi\) is taken as:  
\[\psi=\sqrt{2(p_{EOS}-\rho c_{\mathrm{s}}^{2})/Gc^{2}} \tag{8}\]  
in which \(p_{EOS}\) represents a prescribed equation of state.