## 3 Numerical Framework
### Implicit Time Integration  
The numerical integration of Equation (5) by means of suitable Gaussian quadrature rules leads to a semi-discrete system of nonlinear differential-algebraic equations that can be compactly rewritten as follows  
\[\mathbf{M}_{\boldsymbol{P}}(\mathbf{W})\frac{\mathrm{d}\mathbf{W}}{\mathrm{d}t}+ \mathbf{R}(\mathbf{W})=\mathbf{0}, \tag{17}\]  
where \(\mathbf{M}_{\boldsymbol{P}}(\mathbf{W})\) is the global block diagonal mass matrix derived from the discretization of the first integral in Equation (5) and \(\mathbf{R}(\mathbf{W})\) is the global residual vector pertaining to all of the terms in Equation (5) but the first. Due to the choice of the set of working variables (10), \(\mathbf{M}_{\boldsymbol{P}}\) takes into account the change of variable matrix \(\boldsymbol{P}\). Thus, in two space dimensions (\(d=2\)), the diagonal block of the global mass matrix related to the generic mesh element \(\mathrm{K}\) reads  
\[\mathbf{M}_{\boldsymbol{P}}(\mathbf{W})|_{\mathrm{K}}=\int_{\mathrm{K}}\phi_{i }D_{ja}P_{ak}(\boldsymbol{w}_{h})\phi_{l}\,\mathrm{d}x= \begin{bmatrix}\mathbf{ 0}&\mathbf{0}&\mathbf{0}&\mathbf{0}\\ \mathbf{0}&\mathbf{B}_{\rho}&\mathbf{0}&\mathbf{B}_{\rho\bar{\rho}u_{1}}\\ \mathbf{0}&\mathbf{0}&\mathbf{B}_{\rho}&\mathbf{B}_{\rho\bar{\rho}u_{2}}\\ \mathbf{0}&\mathbf{0}&\mathbf{0}&\mathbf{B}_{\rho\bar{\rho}}\end{bmatrix} ,\]  
with \(i,l=1,\ldots,N_{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{{                                                         } }}}}}}}^{}}}}}}\), and all \(i,l=1,\ldots,N_{{{{{{{{{{{{{{{{{{{{{{\left. }\left {{{{{             }}}}}}}}}}}}}}}}}}}}}\), we obtain the sub-block \(B_{\chi}\) associated with a generic quantity \(\chi\), defined as  
\[(B_{\chi})_{i,l}=\bigg{\{}\int_{\mathrm{K}}\chi\phi_{i}\phi_{l}\,\mathrm{d}x \bigg{\}}.\]  
Notice that due to the lack of the pressure time derivative, the global mass matrix is singular. For the sake of presentation, the matrix \(\left.\mathbf{M}_{\boldsymbol{P}}(\mathbf{W})\right|_{\mathrm{K}}\) can be written in the equivalent form  
\[\left.\mathbf{M}_{\boldsymbol{P}}(\mathbf{W})\right|_{\mathrm{K}}=\lim_{ \boldsymbol{\gamma}\to 0} \begin{bmatrix}\mathbf{B}_{\gamma}&\mathbf{0}& \mathbf{0}&\mathbf{0}\\ \mathbf{0}&\mathbf{B}_{\rho}&\mathbf{0}&\mathbf{B}_{\rho\bar{\rho}u_{1}}\\ \mathbf{0}&\mathbf{0}&\mathbf{B}_{\rho}&\mathbf{B}_{\rho\bar{\rho}u_{2}}\\ \mathbf{0}&\mathbf{0}&\mathbf{0}&\mathbf{B}_{\rho\bar{\rho}}\end{bmatrix} ,\]  
where \(\lim_{\boldsymbol{\gamma}\to 0}\mathbf{B}_{\gamma}=\mathbf{0}\). In this work, we advance the DG discretized equations (17) by using the explicit singly-diagonally implicit Runge-Kutta (ESDIRK) schemes, a class of high-order one-step multi-stage implicit time integration schemes. Such schemes can be designed to have optimal stability properties (\(A\) - and \(L\) - stability) for very high orders, see [27; 28] and, thus, they allow avoiding Courant restrictions typical of explicit schemes. Given the initial condition \(\mathbf{W}^{0}=\mathbf{W}\big{(}t^{0}\big{)}\), the time integration of (17) by means of an \(s-\)stages ESDIRK scheme leads to  
\[\mathbf{W}^{n+1}=\mathbf{W}^{n}-\Delta t^{n}\sum_{\mathrm{i}=1}^{s}b_{ \mathrm{i}}\Big{(}\mathbf{M}_{\boldsymbol{P}}^{-1}\mathbf{R}\Big{)}^{\mathrm{i }}, \tag{18}\]  
where \(\Delta t^{n}=t^{n+1}-t^{n}\) is the \(n\)-th time step size, with \(n\geq 0\). Each stage \(\mathrm{i}=1,\ldots,s\) requires solving the nonlinear system  
\[\mathbf{W}^{\mathrm{i}}=\mathbf{W}^{n}-\Delta t^{n}\sum_{\mathrm{j}=1}^{ \mathrm{i}}a_{\mathrm{ij}}\Big{(}\mathbf{M}_{\boldsymbol{P}}^{-1}\mathbf{R} \Big{)}^{\mathrm{j}}, \tag{19}\]  
with \(\mathbf{W}^{\mathrm{i}}\approx\mathbf{W}(t^{n}+c_{\mathrm{i}}\Delta t^{n})\) the \(\mathrm{i}\)-th stage solution.