## 2 Basic Theory
### Common Approach of TDG Methods  
Consider a one-dimensional conservation law, defined on a domain \(\Omega\), with solution \(\mathbf{u}(x,t)\) and flux \(\mathbf{f}(\mathbf{u})\):  
\[\frac{\partial\mathbf{u}}{\partial t}+\frac{\partial\mathbf{f}}{\partial x}=\mathbf{0}. \tag{1}\]  
Nodal discontinuous Galerkin (nodal DG) schemes typically approximate the solution \(\mathbf{u}(x,t)\) by using its nodal values through a discretization process. However, DG methods offer a more flexible alternative, enabling the approximation of the solution \(\mathbf{u}(x,t)\) without discretization strictly based on nodal values. Specifically, in the \(K\)-th element, where \(x\) lies in the interval \(x\in\left[x^{(k-1)},x^{(k)}\right]\), the \(j\)-th component of the solution \(u_{j}(x,t)\) can be approximated by a dot product of the expansion coefficient vector \(\mathbf{u}_{j}^{(k)}(t)\) and a basis function vector \(\mathbf{\varphi}(x)\). In this study, the basis functions \(\mathbf{\varphi}(x)\) are chosen as Legendre polynomials [16], which provide orthogonality and optimal convergence properties for polynomial-based approximations. Given that the basis function vector \(\varphi\) has N components, the approximation achieves a polynomial order of \((N-1)\). Thus, the approximate solution \(\overline{u}(t)\) within the \(K\)-th element is expressed as:  
\[\overline{u}(x,t)=\left\{ \begin{matrix}\vdots\\ \overline{u}_{j}(x,t)\\ \vdots\end{matrix} \right\}=\left\{ \begin{matrix}\vdots\\ \sum\limits_{i=1}^{N}u_{ii}^{(k)}(t)\varphi_{i}^{(k)}(x)\\ \vdots\end{matrix} \right\}=\left\{ \begin{matrix}\vdots\\ \boldsymbol{u}_{j}^{(k)}\cdot\boldsymbol{\varphi}\\ \vdots\end{matrix} \right\}=\boldsymbol{\Phi}^{\text{T}}\boldsymbol{U}, \tag{2}\]  
where \(\boldsymbol{\Phi}\) is the basis function matrix:  
\[\boldsymbol{\Phi}=\left[ \begin{matrix}\ddots&&O\\ &\boldsymbol{\varphi}\\ O&&\ddots\end{matrix} \right]=\left[ \begin{matrix}\vdots\\ \boldsymbol{\varphi}\\ \vdots\end{matrix} \right]. \tag{3}\]  
The numerical solution and basis functions may have different numbers of components for each element and should therefore be distinguished element-wise. However, for better readability of the equations, the superscript \((k)\) is omitted. Specifically, the \(j\)-th component of the \(\overline{u}\) vector, \(\overline{u}_{j}\), should be expressed as \(\overline{u}_{j}^{(k)}(x,t)=\boldsymbol{u}_{j}^{(k)}\cdot\boldsymbol{\varphi} ^{(k)}\) for \(K\)-th element (interval \(x^{(k-1)}\sim x^{(k)}\)). For simplicity, it is instead represented as \(\overline{u}_{j}(x,t)=\boldsymbol{u}_{j}^{(k)}\cdot\boldsymbol{\varphi}\). To ensure that the approximation \(\overline{u}(x,t)\) satisfies the governing equation in a weak form, DG methods impose the following integral condition in each element [1]:  
\[\int_{x^{(k-1)}}^{x^{(k)}}\boldsymbol{\Phi}\left(\frac{\partial\overline{u}}{ \partial t}+\frac{\partial f}{\partial x}\right)\text{d}x=\boldsymbol{0}, \tag{4}\]  
Within this framework, introducing Taylor series expansions for the temporal derivative \(\partial\overline{u}/\partial t\) enables the construction of TDG methods. By applying high-order Taylor expansions in time, TDG schemes can achieve superior temporal accuracy compared to standard time-stepping methods. This section focuses on a conventional third-order TDG approach [12], which relies on Lax-Wendroff-type time integration to represent the inertia term \(\partial\overline{u}/\partial t\) at \(t=t_{n}\), in \(K\)-th element as:  
\[ \begin{array}{c}\overline{u}_{t}^{n}=\frac{\partial\overline{u}}{\partial t }\Big{|}_{t=t_{n}}=\frac{\overline{u}^{n+1}-\overline{u}^{n}}{2}-\frac{ \Delta t}{2\overline{u}^{2}}\Big{|}_{t=t_{n}}-\frac{\Delta t^{2}}{6}\frac{ \partial^{3}\overline{u}}{\partial s^{3}}\Big{|}_{t=t_{n}}\\ =\frac{\overline{u}^{n+1}-\overline{u}^{n}}{\Delta t}-\frac{\Delta t}{2 \overline{u}^{n}}\overline{u}^{n}_{tt}-\frac{\Delta t^{2}}{6}\overline{u}^{n} _{tt},\end{array}  \tag{5}\]  
Therefore, the weak formulation of the Euler equations is:  
\[