Appendix A: Shallow Water State Values in Time
If the bathymetry is assumed to be constant, the shallow water equations are
$$\begin{aligned} \frac{\partial }{\partial t}\begin{bmatrix} h\\ h\vec {v}\\ \end{bmatrix} + \vec {\nabla }_x\cdot \begin{bmatrix} h\vec {v}\\ h\vec {v}\vec {v}^T + \frac{g}{2}h^2\mathcal {\mathcal {I}}\\ \end{bmatrix} =\begin{bmatrix} 0\\ \vec {0}\\ \end{bmatrix}, \end{aligned}$$
(A.1)
where h is the water height, \(\vec {v} = (v_1,v_2)\) are the fluid velocities, g is the gravitational constant, and \(\mathcal {\mathcal {I}}\) is the \(2\times 2\) identity matrix. Here, we ignore the bottom topography, b, as it is constant in time. Thus, any jumps that arise in the current analysis will simply vanish. However, complete details on the development of the spatial entropy stability analysis for the shallow water equations with a bottom topography can be found in, e.g., [17, 24, 49].
Here, we focus on the temporal entropy analysis. The condition to design a discretely entropy conservative temporal state is (2.32). We collect the necessary shallow water quantities for the discrete temporal entropy analysis
$$\begin{aligned} \begin{aligned} \mathbf{U}&= (h,\,hv_1,\,hv_2)^T,\\ \mathbf{W}&= \left( gh-\frac{1}{2}\left( v_1^2+v_2^2\right) ,\,v_1,\,v_2\right) ^T,\\ s&= \frac{hv^2}{2} + \frac{gh^2}{2} + ghb,\\ \varPhi&= \frac{gh^2}{2},\\ \end{aligned} \end{aligned}$$
(A.2)
with the gravitational constant g. Now, from the properties of the jump operator (3.8), we compute the components of \(\mathbf{U}^{\#}\) from
$$\begin{aligned} \llbracket {\mathbf{W}}\rrbracket ^T\mathbf{U}^{\#} = g\llbracket {h}\rrbracket U_1^\# - \left\{ \!\left\{ v_1\right\} \!\right\} \llbracket {v_1}\rrbracket U_1^\#- \left\{ \!\left\{ v_2\right\} \!\right\} \llbracket {v_2}\rrbracket U_1^\# + \llbracket {v_1}\rrbracket U_2^\# + \llbracket {v_2}\rrbracket U_3^{\#} = g\left\{ \!\left\{ h\right\} \!\right\} \llbracket {h}\rrbracket . \end{aligned}$$
(A.3)
Solving, we determine the entropy conservative temporal state to be
$$\begin{aligned} \mathbf{U}^\# = \begin{bmatrix} \left\{ \!\left\{ h\right\} \!\right\} \\ \left\{ \!\left\{ h\right\} \!\right\} \left\{ \!\left\{ v_1\right\} \!\right\} \\ \left\{ \!\left\{ h\right\} \!\right\} \left\{ \!\left\{ v_2\right\} \!\right\} \\ \end{bmatrix}. \end{aligned}$$
(A.4)
It is easy to verify that the temporal state \(\mathbf{U}^{\#}\) is consistent and symmetric with respect to its arguments.
Just as in the Euler case, the upwind temporal state is used in practice to make the space–time DG numerical scheme computationally tractable. The upwind temporal state is given by
$$\begin{aligned} \mathbf{U}^* = \begin{bmatrix} h_{-}\\ h_{-}(v_1)_{-}\\ h_{-}(v_2)_{-}\\ \end{bmatrix}. \end{aligned}$$
(A.5)
Appendix B: Ideal Magnetohydrodynamics Time State Evaluation
Here, we consider a slightly modified version of the ideal magnetohydrodynamic (MHD) equations
$$\begin{aligned} \frac{\partial }{\partial t}\begin{bmatrix} \rho \\ \rho \vec {v}\\ E\\ \vec {B}\\ \end{bmatrix} + \vec {\nabla }_x\cdot \begin{bmatrix} \rho \vec {v}\\ \rho \vec {v}\vec {v}^T + \left( p + \frac{1}{2}|\vec {B}|^2\right) \mathcal {\mathcal {I}} - \vec {B}\vec {B}^T\\ \vec {v}\left( E+p+\frac{1}{2}|\vec {B}|^2\right) -\vec {B}\left( \vec {v}\cdot \vec {B}\right) \\ \vec {v}\vec {B}^T - \vec {B}\vec {v}^T\\ \end{bmatrix} = \left( \vec {\nabla }_x\cdot \vec {B}\right) \begin{bmatrix} 0\\ \vec {B}\\ \vec {v}\cdot \vec {B}\\ \vec {v}\\ \end{bmatrix} = (\vec {\nabla }_x\cdot \vec {B})\varvec{\varUpsilon }, \end{aligned}$$
(B.1)
where \(\rho \) is the density, \(\vec {v} = (v_1,\,v_2,\,v_3)^T\) are the fluid velocities, \(\vec {B}=(B_1,\,B_2,\,B_3)^T\) are the magnetic fields, E is the total energy, and \(\mathcal {\mathcal {I}}\) is the \(3\times 3\) identity matrix. With the assumption of an ideal gas, the pressure is
$$\begin{aligned} p = (\gamma -1)\left( E - \frac{\rho }{2}|\vec {v}|^2-\frac{1}{2}|\vec {B}|^2\right) . \end{aligned}$$
(B.2)
Also, the divergence-free constraint for magnetized fluids is incorporated into the model (B.1) as a nonconservative term [1, 42]. This nonconservative term is an important aspect of the spatial entropy analysis particularly when \(\vec {\nabla }_x\cdot \vec {B}\ne 0\) see, e.g., [2, 8, 12, 35]. We address how this nonconservative term, proportional to \(\vec {\nabla }_x\cdot \vec {B}\), affects the spatial part of the DG approximation in “Appendix B.1”.
For now, we concern ourselves with the temporal entropy analysis and collect the necessary quantities for the ideal MHD equations
$$\begin{aligned} \begin{aligned} \mathbf{U}&= (\rho ,\,\rho v_1,\,\rho v_2,\,\rho v_3,\,E,\,B_1,\,B_2,\,B_3)^T,\\ \mathbf{W}&= \left( \frac{\gamma -\varsigma }{\gamma -1}-\beta |\vec {v}|^2,\,2\beta v_1,\,2\beta v_2,\,2\beta v_3,\,-2\beta ,\,2\beta B_1,\,2\beta B_2,\,2\beta B_3\right) ^T,\\ s&= -\frac{\rho \varsigma }{\gamma -1},\\ \varPhi&= \rho + \beta \left| \mathbf{B}\right| ^2\\ \end{aligned} \end{aligned}$$
(B.3)
where the physical entropy \(\varsigma \) and \(\beta \) are defined as in the Euler case in Sect. 3.1.
To determine an entropy conservative temporal state, we must find a two-point function \(\mathbf{u}^\#\) such that
$$\begin{aligned} \llbracket {\mathbf{W}}\rrbracket ^T\mathbf{U}^\# = \llbracket {\varPhi }\rrbracket = \llbracket {\rho + \beta \left| \mathbf{B}\right| ^2}\rrbracket . \end{aligned}$$
(B.4)
To do so, we first compute the jump in the entropy variables
$$\begin{aligned} \llbracket {\mathbf{W}}\rrbracket = \begin{bmatrix} \frac{\llbracket {\rho }\rrbracket }{\rho ^{\ln }}+\frac{\llbracket {\beta }\rrbracket }{\beta ^{\ln }(\gamma -1)} - 2\left\{ \!\left\{ \beta \right\} \!\right\} \left( \left\{ \!\left\{ v_1\right\} \!\right\} \llbracket {v_1}\rrbracket +\left\{ \!\left\{ v_2\right\} \!\right\} \llbracket {v_2}\rrbracket +\left\{ \!\left\{ v_3\right\} \!\right\} \llbracket {v_3}\rrbracket \right) -\llbracket {\beta }\rrbracket \left( \left\{ \!\left\{ v_1^2\right\} \!\right\} +\left\{ \!\left\{ v_2^2\right\} \!\right\} +\left\{ \!\left\{ v_3^2\right\} \!\right\} \right) \\ 2\left\{ \!\left\{ v_1\right\} \!\right\} \llbracket {\beta }\rrbracket +2\left\{ \!\left\{ \beta \right\} \!\right\} \llbracket {v_1}\rrbracket \\ 2\left\{ \!\left\{ v_2\right\} \!\right\} \llbracket {\beta }\rrbracket +2\left\{ \!\left\{ \beta \right\} \!\right\} \llbracket {v_2}\rrbracket \\ 2\left\{ \!\left\{ v_3\right\} \!\right\} \llbracket {\beta }\rrbracket +2\left\{ \!\left\{ \beta \right\} \!\right\} \llbracket {v_3}\rrbracket \\ -2\llbracket {\beta }\rrbracket \\ 2\left\{ \!\left\{ B_1\right\} \!\right\} \llbracket {\beta }\rrbracket +2\left\{ \!\left\{ \beta \right\} \!\right\} \llbracket {B_1}\rrbracket \\ 2\left\{ \!\left\{ B_2\right\} \!\right\} \llbracket {\beta }\rrbracket +2\left\{ \!\left\{ \beta \right\} \!\right\} \llbracket {B_2}\rrbracket \\ 2\left\{ \!\left\{ B_3\right\} \!\right\} \llbracket {\beta }\rrbracket +2\left\{ \!\left\{ \beta \right\} \!\right\} \llbracket {B_3}\rrbracket \\ \end{bmatrix}.\nonumber \\ \end{aligned}$$
(B.5)
Then we compute the left hand side of (B.4) to be
$$\begin{aligned} \llbracket {\mathbf{W}}\rrbracket ^T\mathbf{U}^\#= & {} U_1^\#\Bigg (\frac{\llbracket {\rho }\rrbracket }{\rho ^{\ln }}+\frac{\llbracket {\beta }\rrbracket }{\beta ^{\ln }(\gamma -1)} - 2\left\{ \!\left\{ \beta \right\} \!\right\} \big (\left\{ \!\left\{ v_1\right\} \!\right\} \llbracket {v_1}\rrbracket +\left\{ \!\left\{ v_2\right\} \!\right\} \llbracket {v_2}\rrbracket +\left\{ \!\left\{ v_3\right\} \!\right\} \llbracket {v_3}\rrbracket \big )\nonumber \\&\quad -\llbracket {\beta }\rrbracket \left( \left\{ \!\left\{ v_1^2\right\} \!\right\} +\left\{ \!\left\{ v_2^2\right\} \!\right\} +\left\{ \!\left\{ v_3^2\right\} \!\right\} \right) \Bigg )\nonumber \\&+\,U_2^\#\left( 2\left\{ \!\left\{ v_1\right\} \!\right\} \llbracket {\beta }\rrbracket +2\left\{ \!\left\{ \beta \right\} \!\right\} \llbracket {v_1}\rrbracket \right) + U_3^\#\left( 2\left\{ \!\left\{ v_2\right\} \!\right\} \llbracket {\beta }\rrbracket +2\left\{ \!\left\{ \beta \right\} \!\right\} \llbracket {v_2}\rrbracket \right) \nonumber \\&\quad + U_4^\#\left( 2\left\{ \!\left\{ v_3\right\} \!\right\} \llbracket {\beta }\rrbracket +2\left\{ \!\left\{ \beta \right\} \!\right\} \llbracket {v_3}\rrbracket \right) \nonumber \\&-\,2U_5^\#\llbracket {\beta }\rrbracket + U_6^\#\left( 2\left\{ \!\left\{ B_1\right\} \!\right\} \llbracket {\beta }\rrbracket +2\left\{ \!\left\{ \beta \right\} \!\right\} \llbracket {B_1}\rrbracket \right) + U_7^\#\left( 2\left\{ \!\left\{ B_2\right\} \!\right\} \llbracket {\beta }\rrbracket +2\left\{ \!\left\{ \beta \right\} \!\right\} \llbracket {B_2}\rrbracket \right) \nonumber \\&+\, U_8^\#\left( 2\left\{ \!\left\{ B_3\right\} \!\right\} \llbracket {\beta }\rrbracket +2\left\{ \!\left\{ \beta \right\} \!\right\} \llbracket {B_3}\rrbracket \right) . \end{aligned}$$
(B.6)
Next, we expand the right hand side of (B.4)
$$\begin{aligned} \llbracket {\rho + \beta \left| \mathbf{B}\right| ^2}\rrbracket= & {} \llbracket {\rho }\rrbracket + \left( \left\{ \!\left\{ B_1^2\right\} \!\right\} +\left\{ \!\left\{ B_2^2\right\} \!\right\} +\left\{ \!\left\{ B_3^2\right\} \!\right\} \right) \llbracket {\beta }\rrbracket \nonumber \\&+\, 2\left\{ \!\left\{ \beta \right\} \!\right\} \left( \left\{ \!\left\{ B_1\right\} \!\right\} \llbracket {B_1}\rrbracket +\left\{ \!\left\{ B_2\right\} \!\right\} \llbracket {B_2}\rrbracket +\left\{ \!\left\{ B_3\right\} \!\right\} \llbracket {B_3}\rrbracket \right) . \end{aligned}$$
(B.7)
Next, we collect the individual jump terms to create conditions on the entropy conservative temporal state function \(\mathbf{U}^\#\):
$$\begin{aligned} \begin{aligned} \llbracket {\rho }\rrbracket :\quad&\frac{U_1^\#}{\rho ^{\ln }} = 1\\ \llbracket {v_1}\rrbracket :\quad&-2\left\{ \!\left\{ \beta \right\} \!\right\} \left\{ \!\left\{ v_1\right\} \!\right\} U_1^\# + 2\left\{ \!\left\{ \beta \right\} \!\right\} U_2^\# = 0\\ \llbracket {v_2}\rrbracket :\quad&-2\left\{ \!\left\{ \beta \right\} \!\right\} \left\{ \!\left\{ v_2\right\} \!\right\} U_1^\# + 2\left\{ \!\left\{ \beta \right\} \!\right\} U_3^\# = 0\\ \llbracket {v_3}\rrbracket :\quad&-2\left\{ \!\left\{ \beta \right\} \!\right\} \left\{ \!\left\{ v_3\right\} \!\right\} U_1^\# + 2\left\{ \!\left\{ \beta \right\} \!\right\} U_4^\# = 0\\ \llbracket {B_1}\rrbracket :\quad&2\left\{ \!\left\{ \beta \right\} \!\right\} U_6^\# = 2\left\{ \!\left\{ \beta \right\} \!\right\} \left\{ \!\left\{ B_1\right\} \!\right\} \\ \llbracket {B_2}\rrbracket :\quad&2\left\{ \!\left\{ \beta \right\} \!\right\} U_7^\# = 2\left\{ \!\left\{ \beta \right\} \!\right\} \left\{ \!\left\{ B_2\right\} \!\right\} \\ \llbracket {B_3}\rrbracket :\quad&2\left\{ \!\left\{ \beta \right\} \!\right\} U_8^\# = 2\left\{ \!\left\{ \beta \right\} \!\right\} \left\{ \!\left\{ B_3\right\} \!\right\} \\ \llbracket {\beta }\rrbracket :\quad&U_1^\#\left( \frac{1}{\beta ^{\ln }(\gamma -1)}-\left\{ \!\left\{ v_1^2\right\} \!\right\} -\left\{ \!\left\{ v_2^2\right\} \!\right\} -\left\{ \!\left\{ v_3^2\right\} \!\right\} \right) +2U_2^\#\left\{ \!\left\{ v_1\right\} \!\right\} +2U_3^\#\left\{ \!\left\{ v_2\right\} \!\right\} +2U_4^\#\left\{ \!\left\{ v_3\right\} \!\right\} \\&-2U_5^\# + 2U_6^\#\left\{ \!\left\{ B_1\right\} \!\right\} + 2U_7^\#\left\{ \!\left\{ B_2\right\} \!\right\} + 2U_8^\#\left\{ \!\left\{ B_3\right\} \!\right\} = \left\{ \!\left\{ B_1^2\right\} \!\right\} +\left\{ \!\left\{ B_2^2\right\} \!\right\} +\left\{ \!\left\{ B_3^2\right\} \!\right\} \end{aligned} \end{aligned}$$
(B.8)
It is then straightforward to determine
$$\begin{aligned} \mathbf{U}^\# = \begin{bmatrix} \rho ^{\ln }\\ \rho ^{\ln }\left\{ \!\left\{ v_1\right\} \!\right\} \\ \rho ^{\ln }\left\{ \!\left\{ v_2\right\} \!\right\} \\ \rho ^{\ln }\left\{ \!\left\{ v_3\right\} \!\right\} \\ U_5^\#\\ \left\{ \!\left\{ B_1\right\} \!\right\} \\ \left\{ \!\left\{ B_2\right\} \!\right\} \\ \left\{ \!\left\{ B_3\right\} \!\right\} \\ \end{bmatrix}, \end{aligned}$$
(B.9)
with the fifth term
$$\begin{aligned} \begin{aligned} U_5^\#&= \frac{\rho ^{\ln }}{2\beta ^{\ln }(\gamma -1)}+\rho ^{\ln }\left( \left\{ \!\left\{ v_1\right\} \!\right\} ^2+\left\{ \!\left\{ v_2\right\} \!\right\} ^2+\left\{ \!\left\{ v_3\right\} \!\right\} ^2-\frac{1}{2}\left( \left\{ \!\left\{ v_1^2\right\} \!\right\} +\left\{ \!\left\{ v_2^2\right\} \!\right\} +\left\{ \!\left\{ v_3^2\right\} \!\right\} \right) \right) \\&+\left\{ \!\left\{ B_1\right\} \!\right\} ^2+\left\{ \!\left\{ B_2\right\} \!\right\} ^2+\left\{ \!\left\{ B_3\right\} \!\right\} ^2-\frac{1}{2}\left( \left\{ \!\left\{ B_1^2\right\} \!\right\} +\left\{ \!\left\{ B_2^2\right\} \!\right\} +\left\{ \!\left\{ B_3^2\right\} \!\right\} \right) . \end{aligned} \end{aligned}$$
(B.10)
We see that the function \(\mathbf{U}^{\#}\) is symmetric with respect to its arguments and is consistent to the state vector if the left and right states are the same as
$$\begin{aligned} \mathbf{U}^{\#} = \begin{bmatrix} \rho \\ \rho v_1\\ \rho v_2\\ \rho v_3\\ \frac{\rho }{2\beta (\gamma -1)} + \frac{\rho }{2}\Vert \vec {v}\Vert ^2+\frac{1}{2}\Vert \vec {B}\Vert ^2\\ B_1\\ B_2\\ B_3\\ \end{bmatrix} = \begin{bmatrix} \rho \\ \rho v_1\\ \rho v_2\\ \rho v_3\\ \frac{p}{\gamma -1} + \frac{\rho }{2}\Vert \vec {v}\Vert ^2+\frac{1}{2}\Vert \vec {B}\Vert ^2\\ B_1\\ B_2\\ B_3\\ \end{bmatrix} = \begin{bmatrix} \rho \\ \rho v_1\\ \rho v_2\\ \rho v_3\\ E\\ B_1\\ B_2\\ B_3\\ \end{bmatrix} = \mathbf{U}. \end{aligned}$$
(B.11)
Again, the entropy conservative temporal state fully couples all time levels. Thus, it is not a computationally tractable option for approximating the solution of the ideal MHD equations. The entropy stable upwind temporal state given by
$$\begin{aligned} \mathbf{U}^*= & {} \begin{bmatrix} \rho _{-}\\ \rho _{-}(v_1)_{-}\\ \rho _{-}(v_2)_{-}\\ \rho _{-}(v_3)_{-}\\ E_{-}\\ (B_1)_{-}\\ (B_2)_{-}\\ (B_3)_{-}\\ \end{bmatrix},\nonumber \\&E_{-} = \frac{\rho _{-}}{2\beta _{-}(\gamma -1)} + \frac{\rho _{-}}{2}\left( (v_1)_{-}^2+(v_2)_{-}^2+(v_3)_{-}^2\right) + \frac{1}{2}\left( (B_1)_{-}^2+(B_2)_{-}^2+(B_3)_{-}^2\right) \nonumber \\ \end{aligned}$$
(B.12)
is a more feasible choice.
1.1 B.1 Spatial Part for Magnetohydrodynamics
For the ideal magnetohydrodynamic (MHD) equations, the entropy analysis is slightly more complicated due to the divergence-free constraint on the magnetic field
$$\begin{aligned} \vec {\nabla }_x\cdot \vec {B}=0. \end{aligned}$$
(B.13)
The divergence-free constraint must be incorporated into the ideal MHD equations as a nonconservative term in order for the discrete entropy analysis to mimic the continuous entropy analysis, e.g., [1, 2, 8, 12, 35], as it can be that \(\vec {\nabla }_x\cdot \vec {B}\ne 0\) for a numerical discretization.
The inclusion of this nonconservative term alters the derivation of the entropy conservative (or entropy stable) numerical flux in the three spatial directions. Moreover, the entropy flux potential \(\psi \) as shown in (2.38) contains an influence from this nonconservative term and has the form
$$\begin{aligned} \llbracket {\mathbf{W}}\rrbracket ^{T}{\mathbf{F}}_i^{\text {EC}} = \llbracket {\varPsi _i}\rrbracket \quad \text {with} \quad \varPsi _i =\mathbf{W}^{T} \mathbf{F}_i - {F}_i^{s} +\theta B_i, \end{aligned}$$
(B.14)
for \(i=1,2,3\) and
$$\begin{aligned} \theta = \mathbf{w}^T\varvec{\varUpsilon }= 2\beta (\vec {v}\cdot \vec {B}). \end{aligned}$$
(B.15)
See [8, 12] for complete details on the derivation of entropy stable spatial flux functions for the ideal MHD equations.
The DGSEM approximation on general curvilinear meshes discussed briefly in Sect. 2 and fully presented in Gassner et al. [22] remains largely the same for the ideal MHD equations. The only issue is how to compute, at high-order, the volume contributions and surface coupling of the nonconservative term proportional to the divergence-free condition present in (B.1). We will present and briefly discuss the modified spatial operator and how it fits into the present analysis in this work. We note that complete details on the DG approximation for the MHD equations can be found in Bohm et al. [2].
The spatial operator from (2.44) for the ideal MHD equations takes the form
$$\begin{aligned} \begin{aligned} A_{S}\!\left( {\mathbf {U}},\varvec{\varphi }\right) :=&\quad \frac{\varDelta t}{2}\left\langle \vec {{\mathbb {D}}}^{N}\!\cdot \overset{\leftrightarrow }{\tilde{{\mathbf {F}}}}{}^{\text {EC}},\varvec{\varphi }\right\rangle _{\!N\times M}+\frac{\varDelta t}{2}\int \limits _{\partial E^{3},N}\left\langle \varvec{\varphi }^{T}\left( \tilde{{\mathbf {F}}}^{*}-\left( \overset{\leftrightarrow }{\tilde{{\mathbf {F}}}}\cdot {\hat{n}}\right) \right) ,1\right\rangle _{\!M}\,\mathrm {d}S\\&+\frac{\varDelta t}{2}\left\langle \varvec{\varUpsilon }\vec {{\mathbb {D}}}^{N}\!\cdot \vec {{\tilde{B}}},\varvec{\varphi }\right\rangle _{\!N\times M}+\frac{\varDelta t}{2}\int \limits _{\partial E^{3},N}\left\langle \varvec{\varphi }^{T}\left( \left( \varvec{\varUpsilon }\left( \vec {B}\cdot {\hat{n}}\right) \right) ^{\!\!\Diamond }-\varvec{\varUpsilon }\left( \vec {{B}}\cdot {\hat{n}}\right) \right) ,1\right\rangle _{\!M}\,\mathrm {d}S, \end{aligned} \end{aligned}$$
(B.16)
where we have additional contributions in the volume and at the surface of the nonconservative terms. The nonconservative volume contribution is
$$\begin{aligned} \begin{aligned} \varvec{\varUpsilon }\vec {{\mathbb {D}}}^{N}\!\cdot \vec {{\tilde{B}}}_{\sigma ijk} :=2\sum _{m=0}^{N}\quad \,&{\mathcal {D}}_{im}\left( \varvec{\varUpsilon }_{ijk}\left( \vec {B}_{mjk}\cdot \left\{ \!\left\{ J\vec {a}^{1}\right\} \!\right\} _{\left( i,m\right) jk}\right) \right) \\ \quad +&\,{\mathcal {D}}_{jm}\left( \varvec{\varUpsilon }_{ijk}\left( \vec {B}_{imk}\cdot \left\{ \!\left\{ J\vec {a}^{2}\right\} \!\right\} _{i(j,m)k}\right) \right) \\ \quad +&\, {\mathcal {D}}_{km}\left( \varvec{\varUpsilon }_{ijk}\left( \vec {B}_{ijm}\cdot \left\{ \!\left\{ J\vec {a}^{3}\right\} \!\right\} _{ij(k,m)}\right) \right) , \end{aligned} \end{aligned}$$
(B.17)
where \(\sigma = 0,\ldots ,M\). If we select the test function \(\varvec{\phi }\) to be the interpolant of the entropy variables, \(\mathbf{W}\), then the volume contributions become the entropy flux at the boundary [2], i.e.,
$$\begin{aligned} \left\langle \vec {{\mathbb {D}}}^{N}\!\cdot \overset{\leftrightarrow }{\tilde{{\mathbf {F}}}}{}^{\text {EC}},{\mathbf {W}}\right\rangle _{\!N\times M} + \left\langle \varvec{\varUpsilon }\vec {{\mathbb {D}}}^{N}\!\cdot \left( \vec {{\tilde{B}}}\right) ,\mathbf{W}\right\rangle _{\!N\times M}=\int \limits _{\partial E^{3},N}\left\langle \vec {{\tilde{F}}}_{{\hat{n}}}^{s},1\right\rangle _{\!M}\,\mathrm {d}S. \end{aligned}$$
(B.18)
For the nonconservative surface contributions, we define the interface coupling to be
$$\begin{aligned} \left( \varvec{\varUpsilon }\left( \vec {B}\cdot {\hat{n}}\right) \right) ^{\!\!\Diamond } = \varvec{\varUpsilon }^{-}\left\{ \!\left\{ \vec {B}\right\} \!\right\} \cdot {\hat{n}}, \end{aligned}$$
(B.19)
where the “−” denotes that the value is coming from the primary element. We typically suppress this notation unless it is necessary to make the manipulations clear. Furthermore, it holds
$$\begin{aligned} \left\{ \!\left\{ \vec {B}\right\} \!\right\} -\vec {B}_{-}=\frac{1}{2}\llbracket {\vec {B}}\rrbracket . \end{aligned}$$
(B.20)
Then the last term in (B.16) becomes
$$\begin{aligned} \begin{aligned} \int \limits _{\partial E^{3},N}\left\langle \mathbf{W}^{T}\left( \left( \varvec{\varUpsilon }\left( \vec {B}\cdot {\hat{n}}\right) \right) ^{\!\!\Diamond }-\varvec{\varUpsilon }\left( \vec {B}\cdot {\hat{n}}\right) \right) ,1\right\rangle _{\!M}\,\mathrm {d}S&=\int \limits _{\partial E^{3},N}\left\langle \left( \mathbf{W}^{-}\right) ^{T}\varvec{\varUpsilon }^{-}\left( \left\{ \!\left\{ \vec {B}\right\} \!\right\} -\vec {B}\right) \cdot {\hat{n}},1\right\rangle _{\!M}\,\mathrm {d}S\\&=\frac{1}{2}\int \limits _{\partial E^3,N}\left\langle \theta \llbracket {\vec {B}}\rrbracket \cdot {\hat{n}},1\right\rangle _{\!M}\,\mathrm {d}S, \end{aligned} \end{aligned}$$
(B.21)
where we use the definition of \(\theta \) from (B.15). Thus, the property (2.41) becomes
$$\begin{aligned} A_{S}\!\left( \overset{\leftrightarrow }{{\mathbf {F}}},\mathbf{W}\right) :=\quad \frac{\varDelta t}{2}\int \limits _{\partial E^{3},N}\left\langle {\mathbf {W}}^{T}\tilde{{\mathbf {F}}}_{{\hat{n}}}^{*}-{\mathbf {W}}^{T}\tilde{{\mathbf {F}}}_{{\hat{n}}}+{\tilde{F}}_{{\hat{n}}}^{s}+\frac{\theta }{2}\llbracket {\vec {B}}\rrbracket \cdot {\hat{n}},1\right\rangle _{\!M}\,\mathrm {d}S, \end{aligned}$$
(B.22)
in the case of the ideal MHD equations.
If we again introduce the contravariant flux notation (2.53) then summing over all the interior faces gives us
$$\begin{aligned} \sum _{n=1}^{K_{T}}\underset{\text {faces}}{\sum _{\text {Interior}}}\frac{\varDelta t^{n}}{2}\int \limits _{\partial E^{3},N}\left\langle \llbracket {{\mathbf {W}}^{n,k}}\rrbracket ^{T}\tilde{{\mathbf {F}}}_{{\hat{n}}}^{n,k,*}-\llbracket {\left( {\mathbf {W}}^{n,k}\right) ^{T}\tilde{{\mathbf {F}}}_{{\hat{n}}}^{n,k}}\rrbracket +\llbracket {{\tilde{F}}_{{\hat{n}}}^{s}}\rrbracket +\left\{ \!\left\{ \theta \right\} \!\right\} \llbracket {\vec {B}}\rrbracket \cdot {\hat{n}},1\right\rangle _{\!M}\,\mathrm {d}S=0, \end{aligned}$$
(B.23)
due to the entropy stable numerical fluxes for the MHD equations in each Cartesian direction (B.14). We note the additional surface term for the magnetic field components \(\vec {B}\) comes from the orientation of the jump operator and that the physical normal vector \({\hat{n}}\) is defined uniquely at the interface [2]. By summing over all space–time elements and applying the inequality (2.54), we obtain
$$\begin{aligned} \begin{aligned}&\sum _{n=1}^{K_{T}}\sum _{k=1}^{K_{S}}A_{S}\!\left( \overset{\leftrightarrow }{{\mathbf {F}}}{}^{n,k},{\mathbf {W}}^{n,k}\right) \\&\quad = \sum _{n=1}^{K_{T}}\underset{\text {faces}}{\sum _{\text {Boundary}}}\frac{\varDelta t^{n}}{2}\int \limits _{\partial E^{3},N}\left\langle \left( {\mathbf {W}}^{n,k}\right) ^{T}\tilde{{\mathbf {F}}}_{{\hat{n}}}^{n,k,*}-\left( {\mathbf {W}}^{n,k}\right) ^{T}\tilde{{\mathbf {F}}}_{{\hat{n}}}+{\tilde{F}}_{{\hat{n}}}^{n,k,s}\right. \\&\qquad \left. +\frac{\theta ^{-}}{2}\left( \vec {B}^{\text {bndy}}-\vec {B}^{-}\right) \cdot {\hat{n}},1\right\rangle _{\!M}\,\mathrm {d}S\\&\qquad -\sum _{n=1}^{K_{T}}\underset{\text {faces}}{\sum _{\text {Interior}}}\frac{\varDelta t^{n}}{2}\int \limits _{\partial E^{3},N}\left\langle \llbracket {{\mathbf {W}}^{n,k}}\rrbracket ^{T}\tilde{{\mathbf {F}}}_{{\hat{n}}}^{n,k,*}-\llbracket {\left( {\mathbf {W}}^{n,k}\right) ^{T}\tilde{{\mathbf {F}}}_{{\hat{n}}}^{n,k}}\rrbracket +\llbracket {{\tilde{F}}_{{\hat{n}}}^{n,k,s}}\rrbracket \right. \\&\qquad \left. +\left\{ \!\left\{ \theta \right\} \!\right\} \llbracket {\vec {B}}\rrbracket \cdot {\hat{n}},1\right\rangle _{\!M}\,\mathrm {d}S\\&\quad = \sum _{n=1}^{K_{T}}\underset{\text {faces}}{\sum _{\text {Boundary}}}\frac{\varDelta t^{n}}{2}\int \limits _{\partial E^{3},N}\left\langle \left( {\mathbf {W}}^{n,k}\right) ^{T}\tilde{{\mathbf {F}}}_{{\hat{n}}}^{n,k,*}-\left( {\mathbf {W}}^{n,k}\right) ^{T}\tilde{{\mathbf {F}}}_{{\hat{n}}}+{\tilde{F}}_{{\hat{n}}}^{n,k,s}\right. \\&\qquad \left. +\frac{\theta ^{-}}{2}\left( \vec {B}^{\text {bndy}}-\vec {B}^{-}\right) \cdot {\hat{n}},1\right\rangle _{\!M}\,\mathrm {d}S, \end{aligned} \end{aligned}$$
(B.24)
where \(\vec {B}^{\text {bndy}}\) is the magnetic field evaluated from an appropriate boundary state. So, we obtain a discrete entropy equality of a similar form as (2.66) up to the prescription of proper boundary conditions for the ideal MHD equations.
We note that the entropy preserving MHD fluxes are not enough to ensure a stable discretization for discontinuous solutions. A dissipation operator, like the one found in [2], must be added to the entropy preserving flux functions. Then, the equations (B.23) and (B.24) become inequalities. However, even with these inequalities, a discrete entropy inequality of a similar form as (2.71) can be proven.
Appendix C: Proofs for the Temporal Entropy Analysis
In this section, we apply the following identities which result from the properties of the SBP operator \(\mathcal {\mathcal {Q}}\)
$$\begin{aligned} \sum _{i,j=0}^{K}\mathcal {\mathcal {Q}}_{ij}\llbracket {a}\rrbracket _{\left( i,j\right) }=&\qquad \sum _{i,j=0}^{n}\mathcal {\mathcal {B}}_{ij}a_{j}=-\left. a\right| _{-1}^{\,1}, \end{aligned}$$
(C.1)
$$\begin{aligned} \sum _{i,j=0}^{K}\mathcal {\mathcal {Q}}_{ij}\llbracket {a}\rrbracket _{\left( i,j\right) }\left\{ \!\left\{ b\right\} \!\right\} _{\left( i,j\right) }=&\qquad \sum _{i,j=0}^{K}\mathcal {\mathcal {Q}}_{ij}a_{i}b_{j}-\left. ab\right| _{-1}^{\,1}=-\sum _{i,j=0}^{K}\mathcal {\mathcal {Q}}_{ij}a_{j}b_{i}, \end{aligned}$$
(C.2)
$$\begin{aligned} \sum _{i,j=0}^{K}\mathcal {\mathcal {Q}}_{ij}\llbracket {a}\rrbracket _{\left( i,j\right) }\left\{ \!\left\{ b\right\} \!\right\} _{\left( i,j\right) }\left\{ \!\left\{ c\right\} \!\right\} _{\left( i,j\right) } =&-\frac{1}{2}\sum _{i,j=0}^{K}\mathcal {\mathcal {Q}}_{ij}a_{j}b_{i}c_{j}+\frac{1}{2}\sum _{i,j=0}^{n}\mathcal {\mathcal {Q}}_{ij}a_{i}b_{i}c_{j}-\frac{1}{2}\sum _{i,j=0}^{n}\mathcal {\mathcal {Q}}_{ij}a_{j}b_{i}c_{i}, \end{aligned}$$
(C.3)
where \(K=M;N\) and \(\left\{ a\right\} _{i=0}^{K}\), \(\left\{ b\right\} _{i=0}^{K}\) and \(\left\{ c\right\} _{i=0}^{K}\) are generic nodal values. These identities can be proven in a similar way as the discrete split forms in Lemma 1 in [23]. Thus, we skip a proof in this paper.
1.1 C.1 Proof of the Inequality (2.59)
The coefficients of the SBP operator \(\mathcal {\mathcal {Q}}\) are given by \(\mathcal {\mathcal {Q}}_{\sigma \theta }=\omega _{\sigma }\mathcal {\mathcal {D}}_{\sigma \theta }\) and the SBP property (2.30) supplies \(2\mathcal {\mathcal {Q}}_{\sigma \theta }=\mathcal {\mathcal {Q}}_{\sigma \theta }-\mathcal {\mathcal {Q}}_{\theta \sigma }+\mathcal {\mathcal {B}}_{\sigma \theta }\). Thus, it follows
$$\begin{aligned} \begin{aligned} \left\langle J\mathbf{{\mathbb {D}}}^{M}{\mathbf {U}}^{\text {EC}},{\mathbf {W}}\right\rangle _{\!N\times M} =&\quad \sum _{i,j,k=0}^{N}\omega _{ijk}J_{ijk}\sum _{\sigma ,\theta =0}^{M}2\mathcal {\mathcal {Q}}_{\sigma \theta }{\mathbf {W}}_{\sigma ijk}^{T}{\mathbf {U}}^{\#}\left( {\mathbf {U}}_{\sigma ijk},{\mathbf {U}}_{\theta ijk}\right) \\ =&\quad \sum _{i,j,k=0}^{N}\omega _{ijk}J_{ijk}\sum _{\sigma ,\theta =0}^{M}\mathcal {\mathcal {Q}}_{\sigma \theta }{\mathbf {W}}_{\sigma ijk}^{T}{\mathbf {U}}^{\#}\left( {\mathbf {U}}_{\sigma ijk},{\mathbf {U}}_{\theta ijk}\right) \\&-\sum _{i,j,k=0}^{N}\omega _{i,j,k}J_{ijk}\sum _{\sigma ,\theta =0}^{M}\mathcal {\mathcal {Q}}_{\theta \sigma }{\mathbf {W}}_{\sigma ijk}^{T}{\mathbf {U}}^{\#}\left( {\mathbf {U}}_{\sigma ijk},{\mathbf {U}}_{\theta ijk}\right) \\&+\sum _{i,j,k=0}^{N}\omega _{ijk}J_{ijk}\sum _{\sigma ,\theta =0}^{M}\mathcal {\mathcal {B}}_{\sigma \theta }{\mathbf {W}}_{\sigma ijk}^{T}{\mathbf {U}}^{\#}\left( {\mathbf {U}}_{\sigma ijk},{\mathbf {U}}_{\theta ijk}\right) . \end{aligned} \end{aligned}$$
(C.4)
The two-point state \({\mathbf {U}}^{\#}\) is symmetric, hence we obtain by the substitution \(\sigma \longleftrightarrow \theta \) in the penultimate sum in Eq. (C.4)
$$\begin{aligned} \begin{aligned} \left\langle J\mathbf{{\mathbb {D}}}^{M}{\mathbf {U}}^{\text {EC}},{\mathbf {W}}\right\rangle _{\!N\times M} =&\quad \sum _{ijk=0}^{N}\omega _{ijk}J_{ijk}\sum _{\sigma ,\theta =0}^{M}2\omega _{\sigma }\mathcal {\mathcal {D}}_{\sigma \theta }{\mathbf {W}}_{\sigma ijk}^{T}{\mathbf {U}}^{\#}\left( {\mathbf {U}}_{\sigma ijk},{\mathbf {U}}_{\theta ijk}\right) \\ =&\quad \sum _{ijk=0}^{N}\omega _{ijk}J_{ijk}\sum _{\sigma ,\theta =0}^{M}\mathcal {\mathcal {Q}}_{\sigma \theta }\llbracket {{\mathbf {W}}}\rrbracket _{\left( \sigma ,\theta \right) ijk}^{T}{\mathbf {U}}^{\#}\left( {\mathbf {U}}_{\sigma ijk},{\mathbf {U}}_{\theta ijk}\right) \\&+\sum _{ijk=0}^{N}\omega _{ijk}J_{ijk}\sum _{\sigma ,\theta =0}^{M}\mathcal {\mathcal {B}}_{\sigma \theta }{\mathbf {W}}_{\sigma ijk}^{T}{\mathbf {U}}^{\#}\left( {\mathbf {U}}_{\sigma ijk},{\mathbf {U}}_{\theta ijk}\right) . \end{aligned} \end{aligned}$$
(C.5)
Next, since \({\mathbf {U}}^{\#}\) is consistent, it follows
$$\begin{aligned} \sum _{ijk=0}^{N}\omega _{ijk}J_{ijk}\sum _{\sigma ,\theta =0}^{M}\mathcal {\mathcal {B}}_{\sigma \theta }{\mathbf {W}}_{\sigma ijk}^{T}{\mathbf {U}}^{\#}\left( {\mathbf {U}}_{\sigma ijk},{\mathbf {U}}_{\theta ijk}\right) =\left. \left\langle {\mathbf {W}}^{T}{\mathbf {U}},J\right\rangle _{\!N}\right| _{-1}^{\,1}. \end{aligned}$$
(C.6)
Furthermore, since \({\mathbf {U}}^{\#}\) satisfies the property (2.32), we obtain by the identity (C.1)
$$\begin{aligned} \begin{aligned}&\sum _{ijk=0}^{N}\omega _{ijk}J_{ijk}\sum _{\sigma ,\theta =0}^{M}\mathcal {\mathcal {Q}}_{\sigma \theta }\llbracket {{\mathbf {W}}}\rrbracket _{\left( \sigma ,\theta \right) ijk}^{T}{\mathbf {U}}^{\#}\left( {\mathbf {U}}_{\sigma ijk},{\mathbf {U}}_{\theta ijk}\right) \\&\quad = \sum _{i,j,k=0}^{N}\omega _{ijk}J_{ijk}\sum _{\sigma ,\theta =0}^{M}\mathcal {\mathcal {Q}}_{\sigma \theta }\llbracket {\varPhi \left( {\mathbf {U}}\right) }\rrbracket _{\left( \sigma ,\theta \right) ijk} \\&\qquad =-\sum _{i,j,k=0}^{N}\omega _{ijk}J_{ijk}\left. \varPhi \left( {\mathbf {U}}\right) \right| _{-1}^{\, 1} =-\left. \left\langle \varPhi \left( {\mathbf {U}}\right) ,J\right\rangle _{\!N}\right| _{-1}^{\,1}. \end{aligned} \end{aligned}$$
(C.7)
Next, by plugging the equations (C.6) and (C.7) in (C.5), we obtain the identity
$$\begin{aligned} \left\langle J\mathbf{{\mathbb {D}}}^{M}{\mathbf {U}}^{\text {EC}},{\mathbf {W}}\right\rangle _{\!N\times M} = \left. \left\langle {\mathbf {W}}^{T}{\mathbf {U}},J\right\rangle _{\!N}\right| _{-1}^{\,1}-\left. \left\langle \varPhi \left( {\mathbf {U}}\right) ,J\right\rangle _{\!N}\right| _{-1}^{\,1}=\left. \left\langle s\left( {\mathbf {U}}\right) ,J\right\rangle _{\!N}\right| _{-1}^{\,1}. \end{aligned}$$
(C.8)
Finally, a summation of the temporal part (2.43) over all space–time elements and the identity (C.8) provide
$$\begin{aligned} \begin{aligned} \sum _{n=1}^{K_{T}}\sum _{k=1}^{K_{S}}A_{T}\!\left( {\mathbf {U}}^{n,k},{\mathbf {W}}^{n,k}\right) =&\quad {\bar{S}}\left( {\mathbf {U}}\left( T\right) \right) + \sum _{k=1}^{K_{S}}\left. \left\langle \left( {\mathbf {W}}^{K_{T},k}\right) ^{T}\left( {\mathbf {U}}^{K_{T},k,*}-{\mathbf {U}}^{K_{T},k}\right) ,J^{k}\right\rangle _{\!N}\right| _{\tau =1} \\&\qquad -\sum _{n=2}^{K_{T}}\sum _{k=1}^{K_{S}}\left. \left\langle \llbracket {{\mathbf {W}}^{n,k}}\rrbracket ^{T}{\mathbf {U}}^{n,k,*}-\llbracket {\varPhi \left( {\mathbf {U}}^{n,k}\right) }\rrbracket ,J^{k}\right\rangle _{\!N}\right| _{\tau =-1} \\&\quad -{\bar{S}}\left( {\mathbf {U}}\left( 0\right) \right) -\sum _{k=1}^{K_{S}}\left. \left\langle \left( {\mathbf {W}}^{1,k}\right) ^{T}\left( {\mathbf {U}}^{1,k,*}-{\mathbf {U}}^{1,k}\right) ,J^{k}\right\rangle _{\!N}\right| _{\tau =-1}, \end{aligned} \end{aligned}$$
(C.9)
where the quantity \({\bar{S}}\left( {\mathbf {U}}\left( T\right) \right) \) is defined in (2.57) and the quantity \({\bar{S}}\left( {\mathbf {U}}\left( 0\right) \right) \) is defined in (2.60).
This completes the proof for the identity (2.59).
Appendix D: Proofs for the Temporal Kinetic Energy Analysis
1.1 D.1 Proof of Theorem 4 (Kinetic Energy Preservation)
We choose \(\varvec{\varphi }={\mathbf {V}}\) as test function in the Eq. (2.42) and sum over all space–time elements to obtain
$$\begin{aligned} \sum _{n=1}^{K_{T}}\sum _{k=1}^{K_{S}}\left( A_{T}\!\left( {\mathbf {U}}^{n,k},{\mathbf {V}}^{n,k}\right) +A_{S}\!\left( \overset{\leftrightarrow }{{\mathbf {F}}}{}^{n,k},{\mathbf {V}}^{n,k}\right) \right) =0. \end{aligned}$$
(D.1)
First, we investigate the temporal part of the space–time DGSEM. We proceed as in the proof of the identity (2.59) and obtain by the SBP property (2.30)
$$\begin{aligned} \begin{aligned} \left\langle J\mathbf{{\mathbb {D}}}^{M}{\mathbf {U}}^{\text {KEP}},{\mathbf {V}}\right\rangle _{N\times M} =&\quad \sum _{ijk=0}^{N}\omega _{ijk}J_{ijk}\sum _{\sigma ,\theta =0}^{M}2 \mathcal {\mathcal {Q}}_{\sigma \theta }{\mathbf {V}}_{\sigma ijk}^{T}{\mathbf {U}}^{\text {KEP}}\left( {\mathbf {U}}_{\sigma ijk},{\mathbf {U}}_{\theta ijk}\right) \\ =&\quad \sum _{ijk=0}^{N}\omega _{ijk}J_{ijk}\sum _{\sigma ,\theta =0}^{M}\mathcal {\mathcal {Q}}_{\sigma \theta }\llbracket {{\mathbf {V}}}\rrbracket _{\left( \sigma ,\theta \right) ijk}^{T}{\mathbf {U}}^{\text {KEP}}\left( {\mathbf {U}}_{\sigma ijk},{\mathbf {U}}_{\theta ijk}\right) \\&+\sum _{ijk=0}^{N}\omega _{ijk}J_{ijk}\sum _{\sigma ,\theta =0}^{M}\mathcal {\mathcal {B}}_{\sigma \theta }{\mathbf {V}}_{\sigma ijk}^{T}{\mathbf {U}}^{\text {KEP}}\left( {\mathbf {U}}_{\sigma ijk},{\mathbf {U}}_{\theta ijk}\right) . \end{aligned} \end{aligned}$$
(D.2)
The definition of the quantity \({\mathbf {V}}\) provides
$$\begin{aligned} \begin{aligned}&\sum _{ijk=0}^{N}\omega _{ijk}J_{ijk}\sum _{\sigma ,\theta =0}^{M}\mathcal {\mathcal {Q}}_{\sigma \theta }\llbracket {\mathbf{V}}\rrbracket _{\left( \sigma ,\theta \right) ijk}^{T}{\mathbf {U}}^{\text {KEP}}\left( {\mathbf {U}}_{\sigma ijk},{\mathbf {U}}_{\theta ijk}\right) \\&\quad = -\sum _{ijk=0}^{N}\omega _{ijk}J_{ijk}\sum _{\sigma ,\theta =0}^{M}\mathcal {\mathcal {Q}}_{\sigma \theta }\frac{1}{2}\llbracket {\left| \vec {v}\right| ^{2}}\rrbracket _{\left( \sigma ,\theta \right) ijk}{\mathbf {U}}_{1}^{\text {KEP}}\left( {\mathbf {U}}_{\sigma ijk},{\mathbf {U}}_{\theta ijk}\right) \\&\qquad +\sum _{ijk=0}^{N}\omega _{ijk}J_{ijk}\sum _{\sigma ,\theta =0}^{M}\mathcal {\mathcal {Q}}_{\sigma \theta }\llbracket {v_{1}}\rrbracket _{\left( \sigma ,\theta \right) ijk}{\mathbf {U}}_{2}^{\text {KEP}}\left( {\mathbf {U}}_{\sigma ijk},{\mathbf {U}}_{\theta ijk}\right) \\&\qquad +\sum _{ijk=0}^{N}\omega _{ijk}J_{ijk}\sum _{\sigma ,\theta =0}^{M}\mathcal {\mathcal {Q}}_{\sigma \theta }\llbracket {v_{2}}\rrbracket _{\left( \sigma ,\theta \right) ijk}{\mathbf {U}}_{3}^{\text {KEP}}\left( {\mathbf {U}}_{\sigma ijk},{\mathbf {U}}_{\theta ijk}\right) \\&\qquad +\sum _{ijk=0}^{N}\omega _{ijk}J_{ijk}\sum _{\sigma ,\theta =0}^{M}\mathcal {\mathcal {Q}}_{\sigma \theta }\llbracket {v_{3}}\rrbracket _{\left( \sigma ,\theta \right) ijk}{\mathbf {U}}_{4}^{\text {KEP}}\left( {\mathbf {U}}_{\sigma ijk},{\mathbf {U}}_{\theta ijk}\right) . \end{aligned} \end{aligned}$$
(D.3)
Next, it follows by the properties (3.24) of the state \({\mathbf {U}}^{\text {KEP}}\) and the identity (3.8)
$$\begin{aligned} \sum _{ijk=0}^{N}\omega _{ijk}J_{ijk}\sum _{\sigma ,\theta =0}^{M}\mathcal {\mathcal {Q}}_{\sigma \theta }\llbracket {{\mathbf {V}}}\rrbracket _{\left( \sigma ,\theta \right) ijk}^{T}{\mathbf {U}}^{\text {KEP}}\left( {\mathbf {U}}_{\sigma ijk},{\mathbf {U}}_{\theta ijk}\right) =0. \end{aligned}$$
(D.4)
Furthermore, since \({\mathbf {U}}^{\text {KEP}}\) is consistent, we obtain
$$\begin{aligned} \sum _{ijk=0}^{N}\omega _{ijk}J_{ijk}\sum _{\sigma ,\theta =0}^{M}\mathcal {\mathcal {B}}_{\sigma \theta }{\mathbf {V}}_{\sigma ijk}^{T}{\mathbf {U}}^{\#}\left( {\mathbf {U}}_{\sigma ijk},{\mathbf {U}}_{\theta ijk}\right) =\left. \left\langle J{\mathbf {U}},{\mathbf {V}}\right\rangle _{N}\right| _{-1}^{\ \,1}. \end{aligned}$$
(D.5)
Next, we obtain by (D.2), (D.4) and (D.5)
$$\begin{aligned} A_{T}\!\left( {\mathbf {U}},{\mathbf {V}}\right) =\left. \left\langle J{\mathbf {U}}^{*},{\mathbf {V}}\right\rangle _{N}\right| _{-1}^{\,1}. \end{aligned}$$
(D.6)
Since, the space–time DGSEM is applied with Dirichlet boundary conditions in time, a summation of the temporal part (2.43) over all space–time elements and the (D.6) provide
$$\begin{aligned} \begin{aligned} \sum _{k=1}^{K_{T}}\sum _{k=1}^{K_{S}}A_{T}\left( {\mathbf {U}}^{n,k},{\mathbf {V}}^{n,k}\right) =&\quad {\bar{K}}\left( {\mathbf {U}}\left( T\right) \right) -\sum _{k=2}^{K_{T}}\sum _{k=1}^{K_{S}}\left. \left\langle \llbracket {{\mathbf {V}}^{n,k}}\rrbracket ^{T}{\mathbf {U}}^{n,k,*},J^{k}\right\rangle _{N}\right| _{\tau =-1} \\&\qquad \qquad \qquad -\sum _{k=1}^{K_{S}}\left. \left\langle \left( {\mathbf {V}}_{+}^{1,k}\right) ^{T}{\mathbf {U}}^{1,k,*},J^{k}\right\rangle _{N}\right| _{\tau =-1}, \end{aligned} \end{aligned}$$
(D.7)
where \({\bar{K}}\left( {\mathbf {U}}\left( T\right) \right) \) is given by (3.21). The temporal numerical state functions \({\mathbf {U}}^{*}\) have the property (3.24) at the interior temporal boundary points. Thus, we obtain
$$\begin{aligned} \sum _{k=2}^{K_{T}}\sum _{k=1}^{K_{S}}\left. \left\langle \llbracket {{\mathbf {V}}^{n,k}}\rrbracket ^{T}{\mathbf {U}}^{n,k,*},J^{k}\right\rangle _{N}\right| _{\tau =-1}=0. \end{aligned}$$
(D.8)
Furthermore, at the exterior temporal boundary points, the temporal numerical state functions \({\mathbf {U}}^{*}\) satisfy the property (2.78). This provides the identity
$$\begin{aligned} \left. {\mathbf {U}}^{1,k,*}\right| _{\tau =-1}=\left. {\mathbf {U}}_{-}^{1,k}\right| _{\tau =-1}=\left. {\mathbf {u}}^{1,k}\right| _{\tau =-1}, \qquad k=1,\dots ,K_{S}, \end{aligned}$$
(D.9)
where \(\left. {\mathbf {u}}^{1,k}\right| _{\tau =-1}={\mathbf {u}}^{k}\left( 0\right) \) is the initial condition prescribed to the conservation law in the spatial cell \(e_{k}\), \(k=1,\dots ,K_{S}\). Thus, we obtain
$$\begin{aligned} -\sum _{k=1}^{K_{S}}\left. \left\langle \left( {\mathbf {V}}_{+}^{1,k}\right) ^{T}{\mathbf {U}}^{1,k,*},J^{k}\right\rangle _{N}\right| _{\tau =-1} = -{\bar{K}}\left( {\mathbf {u}}\left( 0\right) \right) -\sum _{k=1}^{K_{S}}\left. \left\langle \llbracket {{\mathbf {V}}^{1,k}}\rrbracket ^{T}{\mathbf {u}}^{1,k},J^{k}\right\rangle _{N}\right| _{\tau =-1}, \end{aligned}$$
(D.10)
where \({\bar{K}}\left( {\mathbf {u}}\left( 0\right) \right) \) is given by (3.22). Therefore, the Eq. (D.7) becomes
$$\begin{aligned} \begin{aligned} \sum _{k=1}^{K_{T}}\sum _{k=1}^{K_{S}}A_{T}\left( {\mathbf {U}}^{n,k},{\mathbf {V}}^{n,k}\right) =&\quad {\bar{K}}\left( {\mathbf {U}}\left( T\right) \right) -{\bar{K}}\left( {\mathbf {u}}\left( 0\right) \right) \\&-\sum _{k=1}^{K_{S}}\left. \left\langle \llbracket {{\mathbf {V}}^{1,k}}\rrbracket ^{T}{\mathbf {u}}^{1,k},J^{k}\right\rangle _{N}\right| _{\tau =-1}. \end{aligned} \end{aligned}$$
(D.11)
Next, we consider the spatial part of the space–time DGSEM. In “Appendix D.2”, we prove the identity
$$\begin{aligned} \left\langle \vec {{\mathbb {D}}}^{N}\!\cdot \overset{\leftrightarrow }{\tilde{{\mathbf {F}}}}{}^{\text {KEP}},{\mathbf {V}}\right\rangle _{N\times M}=-\left\langle \vec {\nabla }_{\xi }\cdot {\mathbb {I}}^{N}\!\left( \vec {{\tilde{v}}}\right) ,p\right\rangle _{N\times M}+\int \limits _{\partial E^{3},N}\left\langle {\tilde{F}}_{{\hat{n}}}^{\kappa }+p{\tilde{v}}_{{\hat{n}}},1\right\rangle _{M}\,\mathrm {d}S, \end{aligned}$$
(D.12)
where \({\tilde{F}}^{\kappa }_{{\hat{n}}}:=\vec {{\tilde{F}}}^{\kappa }\cdot {\hat{n}}\) and \(\vec {{\tilde{F}}}^{\kappa }\) is the polynomial approximation of the kinetic energy flux \(\vec {{\tilde{f}}}^{\kappa }\). It follows by (2.45)
$$\begin{aligned} \begin{aligned} {\mathbf {V}}^T\tilde{{\mathbf {F}}}_{{\hat{n}}} =&\quad \; {\hat{s}}n_{1}\left( -\frac{1}{2}\left| \vec {v}\right| ^{2} +\rho v_{1}^{3}+ \rho p v_{1} +\rho v_{1} v_{2}^{2}+\rho v_{1} v_{3}^{2} \right) \\&+{\hat{s}}n_{2}\left( -\frac{1}{2}\left| \vec {v}\right| ^{2} +\rho v_{1}^{2}v_{2} +\rho v_{2}^{3}+ \rho P v_{2}+\rho v_{2} v_{3}^{2} \right) \\&+{\hat{s}}n_{3}\left( -\frac{1}{2}\left| \vec {v}\right| ^{2} +\rho v_{1}^{2}v_{3} +\rho v_{2}^{2}v_{3}+\rho v_{3}^{3}+ \rho p v_{3} \right) \\ =&\frac{1}{2}\left| \vec {v}\right| ^{2}\left( {\hat{s}}\vec {n}\cdot \vec {v}\right) +p\left( {\hat{s}}\vec {n}\cdot \vec {v}\right) ={\tilde{F}}_{{\hat{n}}}^{\kappa }+p{\tilde{v}}_{{\hat{n}}}. \end{aligned} \end{aligned}$$
(D.13)
Thus, we obtain by (D.12) and (D.13)
$$\begin{aligned} \begin{aligned} A_{S}\!\left( \overset{\leftrightarrow }{{\mathbf {F}}},{\mathbf {V}}\right) =&-\frac{\varDelta t}{2}\left\langle \vec {\nabla }_{\xi }\cdot {\mathbb {I}}^{N}\!\left( \vec {{\tilde{v}}}\right) ,p\right\rangle _{N\times M} \\&+\frac{\varDelta t}{2}\int \limits _{\partial E^{3},N}\left\langle {\tilde{F}}_{{\hat{n}}}^{\kappa }+p{\tilde{v}}_{{\hat{n}}},1\right\rangle _{M}\,\mathrm {d}S+\frac{\varDelta t}{2}\int \limits _{\partial E^{3},N}\left\langle {\mathbf {V}}^{T}\tilde{{\mathbf {F}}}_{{\hat{n}}}^{*}-{\tilde{F}}_{{\hat{n}}}^{\kappa }-p{\tilde{v}}_{{\hat{n}}},1\right\rangle _{M} \\ =&-\frac{\varDelta t}{2}\left\langle \vec {\nabla }_{\xi }\cdot {\mathbb {I}}^{N}\!\left( \vec {{\tilde{v}}}\right) ,p\right\rangle _{N\times M} +\frac{\varDelta t}{2}\int \limits _{\partial E^{3},N}\left\langle {\mathbf {V}}^{T}\tilde{{\mathbf {F}}}_{{\hat{n}}}^{*},1\right\rangle _{M}. \end{aligned} \end{aligned}$$
(D.14)
Additionally, the conditions (3.31) and the properties (3.8) of the jump operator provide
$$\begin{aligned} \begin{aligned} \llbracket {{\mathbf {V}}}\rrbracket ^{T}\tilde{{\mathbf {F}}}_{{\hat{n}}}^{*} =&\quad \; {\hat{s}}n_{1}\left( -\frac{1}{2}\llbracket {\left| \vec {v}\right| ^{2}}\rrbracket {\mathbf {F}}_{1}^{1,\text {KEP}}+\llbracket {v_{1}}\rrbracket {\mathbf {F}}_{1}^{2,\text {KEP}}+\llbracket {v_{2}}\rrbracket {\mathbf {F}}_{1}^{3,\text {KEP}}+\llbracket {v_{3}}\rrbracket {\mathbf {F}}_{1}^{4,\text {KEP}}\right) \\&+{\hat{s}}n_{2}\left( -\frac{1}{2}\llbracket {\left| \vec {v}\right| ^{2}}\rrbracket {\mathbf {F}}_{2}^{1,\text {KEP}}+\llbracket {v_{1}}\rrbracket {\mathbf {F}}_{2}^{2,\text {KEP}}+\llbracket {v_{2}}\rrbracket {\mathbf {F}}_{2}^{3,\text {KEP}}+\llbracket {v_{3}}\rrbracket {\mathbf {F}}_{2}^{4,\text {KEP}}\right) \\&+{\hat{s}}n_{3}\left( -\frac{1}{2}\llbracket {\left| \vec {v}\right| ^{2}}\rrbracket {\mathbf {F}}_{3}^{1,\text {KEP}}+\llbracket {v_{1}}\rrbracket {\mathbf {F}}_{3}^{2,\text {KEP}}+\llbracket {v_{2}}\rrbracket {\mathbf {F}}_{3}^{3,\text {KEP}}+\llbracket {v_{1}}\rrbracket {\mathbf {F}}_{3}^{4,\text {KEP}}\right) \\ =&\quad \; {\hat{s}}n_{1}{\mathbf {F}}_{1}^{1,\text {KEP}}\left( -\frac{1}{2}\llbracket {\left| \vec {v}\right| ^{2}}\rrbracket +\llbracket {v_{1}}\rrbracket \left\{ \!\left\{ v_{1}\right\} \!\right\} +\llbracket {v_{1}}\rrbracket \left\{ \!\left\{ p\right\} \!\right\} +\llbracket {v_{2}}\rrbracket \left\{ \!\left\{ v_{2}\right\} \!\right\} +\llbracket {v_{3}}\rrbracket \left\{ \!\left\{ v_{3}\right\} \!\right\} \right) \\&+{\hat{s}}n_{2}{\mathbf {F}}_{1}^{1,\text {KEP}}\left( -\frac{1}{2}\llbracket {\left| \vec {v}\right| ^{2}}\rrbracket +\llbracket {v_{1}}\rrbracket \left\{ \!\left\{ v_{1}\right\} \!\right\} +\llbracket {v_{2}}\rrbracket \left\{ \!\left\{ v_{2}\right\} \!\right\} +\llbracket {v_{2}}\rrbracket \left\{ \!\left\{ p\right\} \!\right\} +\llbracket {v_{3}}\rrbracket \left\{ \!\left\{ v_{3}\right\} \!\right\} \right) \\&+{\hat{s}}n_{3}{\mathbf {F}}_{3}^{1,\text {KEP}}\left( -\frac{1}{2}\llbracket {\left| \vec {v}\right| ^{2}}\rrbracket +\llbracket {v_{1}}\rrbracket \left\{ \!\left\{ v_{1}\right\} \!\right\} +\llbracket {v_{2}}\rrbracket \left\{ \!\left\{ v_{2}\right\} \!\right\} +\llbracket {v_{3}}\rrbracket \left\{ \!\left\{ v_{3}\right\} \!\right\} +\llbracket {v_{3}}\rrbracket \left\{ \!\left\{ p\right\} \!\right\} \right) \\ =&\left\{ \!\left\{ p\right\} \!\right\} \left( {\hat{s}}\vec {n}\cdot \llbracket {\vec {v}}\rrbracket \right) =\left\{ \!\left\{ p\right\} \!\right\} \llbracket {{\tilde{v}}_{{\hat{n}}}}\rrbracket . \end{aligned} \end{aligned}$$
(D.15)
Next, we sum the Eq. (D.14) over all space–time elements and apply the identity (D.15). This results in
$$\begin{aligned} \sum _{n=1}^{K_{T}}\sum _{k=1}^{K_{S}}A_{S}\!\left( \overset{\leftrightarrow }{{\mathbf {F}}}{}^{n,k},{\mathbf {V}}^{n,k}\right)= & {} -\sum _{n=1}^{K_{T}}\sum _{k=1}^{K_{S}}\frac{\varDelta t^{n}}{2}\left\langle \vec {\nabla }_{\xi }\cdot {\mathbb {I}}^{N}\!\left( \vec {{\tilde{v}}}^{n,k}\right) ,p^{n,k}\right\rangle _{N\times M} \nonumber \\&+\sum _{n=1}^{K_{T}}\underset{\text {faces}}{\sum _{\text {Boundary}}}\frac{\varDelta t^{n}}{2}\int \limits _{\partial E^{3},N}\left\langle \left( {\mathbf {V}}^{n,k}\right) ^{T}\tilde{{\mathbf {F}}}_{{\hat{n}}}^{n,k,*},1\right\rangle _{M}\,\mathrm {d}S\nonumber \\&-\sum _{n=1}^{K_{T}}\underset{\text {faces}}{\sum _{\text {Interior}}}\frac{\varDelta t^{n}}{2}\int \limits _{\partial E^{3},N}\left\langle \llbracket {{\mathbf {V}}^{n,k}}\rrbracket ^{T}\tilde{{\mathbf {F}}}_{{\hat{n}}}^{n,k,*},1\right\rangle _{M}\,\mathrm {d}S\nonumber \\= & {} -\sum _{n=1}^{K_{T}}\sum _{k=1}^{K_{S}}\frac{\varDelta t^{n}}{2}\left\langle \vec {\nabla }_{\xi }\cdot {\mathbb {I}}^{N}\!\left( \vec {{\tilde{v}}}^{n,k}\right) ,p^{n,k}\right\rangle _{N\times M} \nonumber \\&-\sum _{n=1}^{K_{T}}\underset{\text {faces}}{\sum _{\text {Interior}}}\frac{\varDelta t^{n}}{2}\int \limits _{\partial E^{3},N}\left\langle \left\{ \!\left\{ p^{n,k}\right\} \!\right\} \llbracket {{\tilde{v}}_{{\hat{n}}}^{n,k}}\rrbracket ,1\right\rangle _{M}\,\mathrm {d}S+{\mathbf{B }}{\mathbf{C }},\nonumber \\ \end{aligned}$$
(D.16)
where
$$\begin{aligned} \mathbf BC :=\sum _{n=1}^{K_{T}}\underset{\text {faces}}{\sum _{\text {Boundary}}}\frac{\varDelta t^{n}}{2}\int \limits _{\partial E^{3},N}\left\langle \left( {\mathbf {V}}^{n,k}\right) ^{T}\tilde{{\mathbf {F}}}_{{\hat{n}}}^{n,k,*},1\right\rangle _{M}\,\mathrm {d}S. \end{aligned}$$
(D.17)
The term \(\mathbf BC \) in (D.17) vanished, since the space–time DGSEM is applied with periodic boundary conditions in space. Finally, we substitute the results (D.11) and (D.16) in (D.1), use that \(\mathbf BC =0\) and rearrange terms. This results in the final equation
$$\begin{aligned} \begin{aligned} {\bar{K}}\left( {\mathbf {U}}\left( T\right) \right)&= {\bar{K}}\left( {\mathbf {U}}\left( 0\right) \right) +\sum _{n=1}^{K_{T}}\sum _{k=1}^{K_{S}}\frac{\varDelta t^{n}}{2}\left\langle \vec {\nabla }_{\xi }\cdot {\mathbb {I}}^{N}\!\left( \vec {{\tilde{v}}}^{n,k}\right) ,p^{n,k}\right\rangle _{N\times M} \\&\quad \qquad \qquad \ \ +\sum _{n=1}^{K_{T}}\underset{\text {faces}}{\sum _{\text {Interior}}}\frac{\varDelta t^{n}}{2}\int \limits _{\partial E^{3},N}\left\langle \left\{ \!\left\{ p^{n,k}\right\} \!\right\} \llbracket {{\tilde{v}}_{{\hat{n}}}^{n,k}}\rrbracket ,1\right\rangle _{M}\,\mathrm {d}S\\&\quad \qquad \qquad \ \ +\sum _{k=1}^{K_{S}}\left. \left\langle \llbracket {{\mathbf {V}}^{1,k}}\rrbracket ^{T}{\mathbf {u}}^{1,k},J^{k}\right\rangle _{N}\right| _{\tau =-1}. \end{aligned} \end{aligned}$$
(D.18)
1.2 D.2 Proof of the Identity (D.12)
The definition of the spatial derivative projection operator (3.26) supplies
$$\begin{aligned} \begin{aligned}&\left\langle \vec {{\mathbb {D}}}^{N}\!\cdot \overset{\leftrightarrow }{\tilde{{\mathbf {F}}}}{}^{\text {KEP}},{\mathbf {V}}\right\rangle _{N\times M} \\&\quad = \sum _{\sigma =0}^{M}\omega _{\sigma } \sum _{i,j,k,m=0}^{N} \bigg \{\omega _{jk}2\mathcal {\mathcal {Q}}_{im}{\mathbf {V}}_{\sigma ijk}^{T}\left( \overset{\leftrightarrow }{{\mathbf {F}}}{}^{\text {KEP}}\left( {\mathbf {U}}_{\sigma ijk},{\mathbf {U}}_{\sigma mjk}\right) \cdot \left\{ \!\left\{ J\vec {a}^{1}\right\} \!\right\} _{\left( i,m\right) jk}\right) \\&\qquad +\,\omega _{ik}2\mathcal {\mathcal {Q}}_{jm}{\mathbf {V}}_{\sigma ijk}^{T}\left( \overset{\leftrightarrow }{{\mathbf {F}}}{}^{\text {KEP}}\left( {\mathbf {U}}_{\sigma ijk},{\mathbf {U}}_{\sigma imk}\right) \cdot \left\{ \!\left\{ J\vec {a}^{2}\right\} \!\right\} _{i\left( j,m\right) k}\right) \\&\qquad +\,\omega _{ij}2\mathcal {\mathcal {Q}}_{km}{\mathbf {V}}_{\sigma ijk}^{T}\left( \overset{\leftrightarrow }{{\mathbf {F}}}{}^{\text {KEP}}\left( {\mathbf {U}}_{\sigma ijk},{\mathbf {U}}_{\sigma ijm}\right) \cdot \left\{ \!\left\{ J\vec {a}^{3}\right\} \!\right\} _{ij\left( k,m\right) }\right) \bigg \}. \end{aligned} \end{aligned}$$
(D.19)
The numerical flux functions satisfy the symmetry condition (3.28). Thus, by the same arguments as in the appendix (C.1) and the SBP property, the first sum on the the right hand side in (D.19) can be written as
$$\begin{aligned} \begin{aligned}&\sum _{\sigma =0}^{M}\omega _{\sigma }\sum _{i,j,k,m=0}^{N} \omega _{jk}2\mathcal {\mathcal {Q}}_{im}{\mathbf {V}}_{\sigma ijk}^{T}\left( \overset{\leftrightarrow }{{\mathbf {F}}}{}^{\text {KEP}}\left( {\mathbf {U}}_{\sigma ijk},{\mathbf {U}}_{\sigma mjk}\right) \cdot \left\{ \!\left\{ J\vec {a}^{1}\right\} \!\right\} _{\left( i,m\right) jk}\right) \\&\quad = \quad \sum _{\sigma =0}^{M}\omega _{\sigma }\sum _{i,j,k,m=0}^{N}\omega _{jk}\mathcal {\mathcal {Q}}_{im}\llbracket {{\mathbf {V}}}\rrbracket _{\sigma ijk}^{T}\left( \overset{\leftrightarrow }{{\mathbf {F}}}{}^{\text {KEP}}\left( {\mathbf {U}}_{\sigma ijk},{\mathbf {U}}_{\sigma mjk}\right) \cdot \left\{ \!\left\{ J\vec {a}^{1}\right\} \!\right\} _{\left( i,m\right) jk}\right) \\&\qquad +\sum _{\sigma =0}^{M}\omega _{\sigma }\sum _{i,j,k,m=0}^{N}\omega _{jk}\mathcal {\mathcal {B}}_{im}{\mathbf {V}}_{\sigma ijk}^{T}\left( \overset{\leftrightarrow }{{\mathbf {F}}}{}^{\text {KEP}}\left( {\mathbf {U}}_{\sigma ijk},{\mathbf {U}}_{\sigma mjk}\right) \cdot \left\{ \!\left\{ J\vec {a}^{1}\right\} \!\right\} _{\left( i,m\right) jk}\right) . \end{aligned} \end{aligned}$$
(D.20)
For the last sum on the right hand side in (D.20), it follows
$$\begin{aligned}&\sum _{\sigma =0}^{M}\omega _{\sigma }\sum _{i,j,k,m=0}^{N}\omega _{jk}\mathcal {\mathcal {B}}_{im}{\mathbf {V}}_{\sigma ijk}^{T}\left( \overset{\leftrightarrow }{{\mathbf {F}}}{}^{\text {KEP}}\left( {\mathbf {U}}_{\sigma ijk},{\mathbf {U}}_{\sigma mjk}\right) \cdot \left\{ \!\left\{ J\vec {a}^{1}\right\} \!\right\} _{\left( i,m\right) jk}\right) \nonumber \\&= \sum _{\sigma =0}^{M}\omega _{\sigma }\sum _{j,k=0}^{N}\omega _{jk}\left\{ \left( \vec {F}_{\sigma Njk}^{\kappa }+p_{\sigma Njk}\vec {v}_{\sigma Njk}\right) ^{T}J\vec {a}_{Njk}^{1}-\left( \vec {F}_{\sigma 0jk}^{\kappa }+p_{\sigma 0jk}\vec {v}_{\sigma 0jk}\right) ^{T}J\vec {a}_{0jk}^{1}\right\} ,\nonumber \\ \end{aligned}$$
(D.21)
since the numerical flux functions are consistent with \(\overset{\leftrightarrow }{{\mathbf {F}}}\) and the equation
$$\begin{aligned} {\mathbf {V}}^{T}\left( \overset{\leftrightarrow }{{\mathbf {F}}}\cdot J\vec {a}^{1}\right) =\left( \vec {F}^{\kappa }+p\vec {v}\right) ^{T} J\vec {a}^{1}, \end{aligned}$$
(D.22)
is satisfied at the nodes. Moreover, we obtain by Jameson’s conditions (3.29) and the property (3.8) of the jump operator the identity
$$\begin{aligned} \begin{aligned}&\llbracket {{\mathbf {V}}}\rrbracket _{\sigma ijk}^{T}\left( \overset{\leftrightarrow }{{\mathbf {F}}}{}^{\text {KEP}}\left( {\mathbf {U}}_{\sigma ijk},{\mathbf {U}}_{\sigma mjk}\right) \cdot \left\{ \!\left\{ J\vec {a}^{1}\right\} \!\right\} _{\left( i,m\right) jk}\right) \\&\quad = \sum _{s=1}^{3}\bigg \{-\frac{1}{2}\llbracket {\left| \vec {v}\right| ^2}\rrbracket _{\sigma \left( i,m\right) jk}{\mathbf {F}}_{s}^{1,\text {KEP}}\left( {\mathbf {U}}_{\sigma ijk},{\mathbf {U}}_{\sigma mjk}\right) \left\{ \!\left\{ J{a}_{s}^{1}\right\} \!\right\} _{\left( i,m\right) jk} \\&\qquad + \left\{ \!\left\{ v_{1}\right\} \!\right\} _{\sigma \left( i,m\right) jk}\llbracket {v_{1}}\rrbracket _{\sigma \left( i,m\right) jk}{\mathbf {F}}_{s}^{1,\text {KEP}}\left( {\mathbf {U}}_{\sigma ijk},{\mathbf {U}}_{\sigma mjk}\right) \left\{ \!\left\{ J{a}_{s}^{1}\right\} \!\right\} _{\left( i,m\right) jk} \\&\qquad +\left\{ \!\left\{ v_{2}\right\} \!\right\} _{\sigma \left( i,m\right) jk}\llbracket {v_{2}}\rrbracket _{\sigma \left( i,m\right) jk}{\mathbf {F}}_{s}^{1,\text {KEP}}\left( {\mathbf {U}}_{\sigma ijk},{\mathbf {U}}_{\sigma mjk}\right) \left\{ \!\left\{ J{a}_{s}^{1}\right\} \!\right\} _{\left( i,m\right) jk} \\&\qquad +\left\{ \!\left\{ v_{3}\right\} \!\right\} _{\sigma \left( i,m\right) jk}\llbracket {v_{3}}\rrbracket _{\sigma \left( i,m\right) jk}{\mathbf {F}}_{s}^{1,\text {KEP}}\left( {\mathbf {U}}_{\sigma ijk},{\mathbf {U}}_{\sigma mjk}\right) \left\{ \!\left\{ J{a}_{s}^{1}\right\} \!\right\} _{\left( i,m\right) jk} \\&\qquad + \llbracket {v_{s}}\rrbracket _{\sigma \left( i,m\right) jk}\left( p_{s1}^{\star }\right) _{\sigma \left( i,m\right) jk}\left\{ \!\left\{ J{a}_{s}^{1}\right\} \!\right\} _{\left( i,m\right) jk} \bigg \} \\&\quad = \quad \llbracket {\vec {v}}\rrbracket _{\sigma \left( i,m\right) jk}^{T}\left( 2\left\{ \!\left\{ p\right\} \!\right\} _{\sigma \left( i,m\right) jk}\left\{ \!\left\{ J\vec {a}^{1}\right\} \!\right\} _{\left( i,m\right) jk}-\left\{ \!\left\{ pJ\vec {a}^{1}\right\} \!\right\} _{\sigma \left( i,m\right) jk}\right) \end{aligned} \end{aligned}$$
(D.23)
for all nodal values, since the Eq. (3.30) provides
$$\begin{aligned}&\left( p_{s1}^{\star }\right) _{\sigma \left( i,m\right) jk}\left\{ \!\left\{ J{a}_{s}^{1}\right\} \!\right\} _{\left( i,m\right) jk}\nonumber \\&\quad =2\left\{ \!\left\{ p\right\} \!\right\} _{\sigma \left( i,m\right) jk}\left\{ \!\left\{ J{a}_{s}^{1}\right\} \!\right\} _{\left( i,m\right) jk}-\left\{ \!\left\{ pJ{a}_{s}^{1}\right\} \!\right\} _{\sigma \left( i,m\right) jk},\quad s=1,2,3. \end{aligned}$$
(D.24)
By combining the identity (D.23) with (C.2) and (C.3), we obtain
$$\begin{aligned} \begin{aligned}&\sum _{\sigma =0}^{M}\omega _{\sigma }\sum _{i,j,k,m=0}^{N}\omega _{jk}\mathcal {\mathcal {Q}}_{im}\llbracket {{\mathbf {V}}}\rrbracket _{\sigma ijk}^{T}\left( \overset{\leftrightarrow }{{\mathbf {F}}}{}^{\text {KEP}}\left( {\mathbf {U}}_{\sigma ijk},{\mathbf {U}}_{\sigma mjk}\right) \cdot \left\{ \!\left\{ J\vec {a}^{1}\right\} \!\right\} _{\left( i,m\right) jk}\right) \\&\quad = \sum _{\sigma =0}^{M}\omega _{\sigma }\sum _{i,j,k,m=0}^{N}\omega _{jk}\mathcal {\mathcal {Q}}_{im} \llbracket {\vec {v}}\rrbracket _{\sigma \left( i,m\right) jk}^{T}\left( 2\left\{ \!\left\{ p\right\} \!\right\} _{\sigma \left( i,m\right) jk}\left\{ \!\left\{ J\vec {a}^{1}\right\} \!\right\} _{\left( i,m\right) jk}-\left\{ \!\left\{ pJ\vec {a}^{1}\right\} \!\right\} _{\sigma \left( i,m\right) jk}\right) \\&\quad = \left\langle p\vec {v},\frac{\partial {\mathbb {I}}^{N}\!\left( J\vec {a}^{1}\right) }{\partial \xi ^{1}}\right\rangle _{N\times M} -\left\langle \frac{\partial {\mathbb {I}}^{N}\!\left( \vec {v}^T J\vec {a}^{1}\right) }{\partial \xi ^{1}},p\right\rangle _{N\times M}. \end{aligned} \end{aligned}$$
(D.25)
Then, we substitute (D.21) and (D.25) in (D.20). This results in the equation
$$\begin{aligned} \begin{aligned}&\sum _{\sigma =0}^{M}\omega _{\sigma }\sum _{i,j,k,m=0}^{N} \omega _{jk}2\mathcal {\mathcal {Q}}_{jm}{\mathbf {V}}_{\sigma ijk}^{T}\left( \overset{\leftrightarrow }{{\mathbf {F}}}{}^{\text {KEP}}\left( {\mathbf {U}}_{\sigma ijk},{\mathbf {U}}_{\sigma mjk}\right) \cdot \left\{ \!\left\{ J\vec {a}^{1}\right\} \!\right\} _{\left( i,m\right) jk}\right) \\&\quad = \left\langle p\vec {v},\frac{\partial {\mathbb {I}}^{N}\!\left( J\vec {a}^{1}\right) }{\partial \xi ^{1}}\right\rangle _{N\times M} -\left\langle \frac{\partial {\mathbb {I}}^{N}\!\left( \vec {v}^T J\vec {a}^{1}\right) }{\partial \xi ^{1}},p\right\rangle _{N\times M} \\&\qquad + \sum _{\sigma =0}^{M}\omega _{\sigma }\sum _{j,k=0}^{N}\omega _{jk}\left\{ \left( \vec {F}_{\sigma Njk}^{\kappa }+p_{\sigma Njk}\vec {v}_{\sigma Njk}\right) ^{T}J\vec {a}_{\sigma Njk}^{1}-\left( \vec {F}_{\sigma 0jk}^{\kappa }+p_{\sigma 0jk}\vec {v}_{\sigma 0jk}\right) ^{T}J\vec {a}_{0jk}^{1}\right\} . \end{aligned} \end{aligned}$$
(D.26)
In the same way, the second and third sum on the right hand side in (D.19) are evaluated. Thus, we also have the identities
$$\begin{aligned} \begin{aligned}&\sum _{\sigma =0}^{M}\omega _{\sigma }\sum _{i,j,k,m=0}^{N} \omega _{ik}2\mathcal {\mathcal {Q}}_{jm}{\mathbf {V}}_{\sigma ijk}^{T}\left( \overset{\leftrightarrow }{{\mathbf {F}}}{}^{\text {KEP}}\left( {\mathbf {U}}_{\sigma ijk},{\mathbf {U}}_{\sigma mjk}\right) \cdot \left\{ \!\left\{ J\vec {a}^{2}\right\} \!\right\} _{i\left( j,m\right) k}\right) \\&\quad = \left\langle p\vec {v},\frac{\partial {\mathbb {I}}^{N}\!\left( J\vec {a}^{2}\right) }{\partial \xi ^{2}}\right\rangle _{N\times M} -\left\langle \frac{\partial {\mathbb {I}}^{N}\!\left( \vec {v}^T J\vec {a}^{2}\right) }{\partial \xi ^{2}},p\right\rangle _{N\times M} \\&\qquad + \sum _{\sigma =0}^{M}\omega _{\sigma }\sum _{i,k=0}^{N}\omega _{ik}\left\{ \left( \vec {F}_{\sigma iNk}^{\kappa }+p_{\sigma iNk}\vec {v}_{\sigma iNk}\right) ^{T}J\vec {a}_{iNk}^{2}-\left( \vec {F}_{\sigma i0k}^{\kappa }+p_{\sigma i0k}\vec {v}_{\sigma i0k}\right) ^{T}J\vec {a}_{i0k}^{2}\right\} \end{aligned} \end{aligned}$$
(D.27)
for the second sum on the right hand side in (D.19) and
$$\begin{aligned} \begin{aligned}&\sum _{\sigma =0}^{M}\omega _{\sigma }\sum _{i,j,k,m=0}^{N} \omega _{ij}2\mathcal {\mathcal {Q}}_{km}{\mathbf {V}}_{\sigma ijk}^{T}\left( \overset{\leftrightarrow }{{\mathbf {F}}}{}^{\text {KEP}}\left( {\mathbf {U}}_{\sigma ijk},{\mathbf {U}}_{\sigma mjk}\right) \cdot \left\{ \!\left\{ J\vec {a}^{3}\right\} \!\right\} _{ij\left( k,m\right) }\right) \\&\quad = \left\langle p\vec {v},\frac{\partial {\mathbb {I}}^{N}\!\left( J\vec {a}^{3}\right) }{\partial \xi ^{3}}\right\rangle _{N\times M} -\left\langle \frac{\partial {\mathbb {I}}^{N}\!\left( \vec {v}^T J\vec {a}^{3}\right) }{\partial \xi ^{3}},p\right\rangle _{N\times M} \\&\qquad + \sum _{\sigma =0}^{M}\omega _{\sigma }\sum _{i,j=0}^{N}\omega _{ij}\left\{ \left( \vec {F}_{\sigma ijN}^{\kappa }+p_{\sigma ikN}\vec {v}_{\sigma ijN}\right) ^{T}J\vec {a}_{ijN}^{3}-\left( \vec {F}_{\sigma ij0}^{\kappa }+p_{\sigma ik0}\vec {v}_{\sigma ij0}\right) ^{T}J\vec {a}_{ij0}^{3}\right\} \end{aligned} \end{aligned}$$
(D.28)
for the third sum on the right hand side in (D.19).
It follows that
$$\begin{aligned} \begin{aligned} \sum _{s=1}^{3}\left\langle \frac{\partial {\mathbb {I}}^{N}\!\left( \vec {v}^T J\vec {a}^{s}\right) }{\partial \xi ^{s}},p\right\rangle _{N\times M}= \left\langle \vec {\nabla }_{\xi }\cdot {\mathbb {I}}^{N}\!\left( \vec {{\tilde{v}}}\right) ,p\right\rangle _{N\times M} \end{aligned} \end{aligned}$$
(D.29)
and by the definition of the discrete surface integrals (2.28)
$$\begin{aligned} \begin{aligned}&\int \limits _{\partial E^{3},N}\left\langle {\tilde{F}}_{{\hat{n}}}^{\kappa }+p{\tilde{v}}_{{\hat{n}}},1\right\rangle _{M}\,\mathrm {d}S\\&\quad = \sum _{\sigma =0}^{M}\omega _{\sigma }\bigg [\sum _{j,k=0}^{N}\omega _{jk}\left\{ \left( \vec {F}_{\sigma Njk}^{\kappa }+p_{\sigma Njk}\vec {v}_{\sigma Njk}\right) ^{T}J\vec {a}_{Njk}^{1}-\left( \vec {F}_{\sigma 0jk}^{\kappa }+p_{\sigma 0jk}\vec {v}_{\sigma 0jk}\right) ^{T}J\vec {a}_{0jk}^{1}\right\} \\&\qquad + \sum _{i,k=0}^{N}\omega _{ik}\left\{ \left( \vec {F}_{\sigma iNk}^{\kappa }+p_{\sigma iNk}\vec {v}_{\sigma iNk}\right) ^{T}J\vec {a}_{\sigma iNk}^{2}-\left( \vec {F}_{\sigma i0k}^{\kappa }+p_{\sigma i0k}\vec {v}_{\sigma i0k}\right) ^{T}J\vec {a}_{i0k}^{2}\right\} \\&\qquad + \sum _{i,j=0}^{N}\omega _{ij}\left\{ \left( \vec {F}_{\sigma ijN}^{\kappa }+p_{\sigma ijN}\vec {v}_{\sigma ijN}\right) ^{T}J\vec {a}_{ijN}^{3}-\left( \vec {F}_{\sigma ij0}^{\kappa }+p_{\sigma ij0}\vec {v}_{\sigma ij0}\right) ^{T}J\vec {a}_{ij0}^{3}\right\} \bigg ]. \end{aligned} \end{aligned}$$
(D.30)
Moreover, we note that the contravariant coordinate vectors are discretized by (2.21) and thus, the discrete metric identities (2.22) are satisfied [30]. Hence, it follows
$$\begin{aligned} \begin{aligned}&\sum _{s=1}^{3}\left\langle p\vec {v},\frac{\partial {\mathbb {I}}^{N}\!\left( J\vec {a}^{s}\right) }{\partial \xi ^{s}}\right\rangle _{N\times M} \\&\quad =\sum _{\sigma =0}^{M}\omega _{\sigma }\bigg [\sum _{i,j,k=0}^{N}\omega _{ijk}p_{\sigma ijk}\left( v_{1}\right) _{\sigma ijk}\left\{ \sum _{m=0}^{N}\mathcal {\mathcal {D}}_{im}\left( J{a}_{1}^{1}\right) _{mjk}+\mathcal {\mathcal {D}}_{jm}\left( J{a}_{1}^{2}\right) _{imk}+\mathcal {\mathcal {D}}_{km}\left( J{a}_{1}^{3}\right) _{ijm}\right\} \\&\qquad +\omega _{ijk}p_{\sigma ijk}\left( v_{2}\right) _{\sigma ijk}\left\{ \sum _{m=0}^{N}\mathcal {\mathcal {D}}_{im}\left( J{a}_{2}^{1}\right) _{mjk}+\mathcal {\mathcal {D}}_{jm}\left( J{a}_{2}^{2}\right) _{imk}+\mathcal {\mathcal {D}}_{km}\left( J{a}_{2}^{3}\right) _{ijm}\right\} \\&\qquad +\omega _{ijk}p_{\sigma ijk}\left( v_{3}\right) _{\sigma ijk}\left\{ \sum _{m=0}^{N}\mathcal {\mathcal {D}}_{im}\left( J{a}_{3}^{1}\right) _{mjk}+\mathcal {\mathcal {D}}_{jm}\left( J{a}_{3}^{2}\right) _{imk}+\mathcal {\mathcal {D}}_{km}\left( J{a}_{3}^{3}\right) _{ijm}\right\} \bigg ]=0. \end{aligned} \end{aligned}$$
(D.31)
Therefore, by plugging (D.26), (D.27) and (D.28) in (D.19), we obtain the desired identity
$$\begin{aligned} \left\langle \vec {{\mathbb {D}}}^{N}\!\cdot \overset{\leftrightarrow }{\tilde{{\mathbf {F}}}}{}^{\text {KEP}},{\mathbf {V}}\right\rangle _{N\times M}=-\left\langle \vec {\nabla }_{\xi }\cdot {\mathbb {I}}^{N}\!\left( \vec {{\tilde{v}}}\right) ,p\right\rangle _{N\times M}+\int \limits _{\partial E^{3},N}\left\langle {\tilde{F}}_{{\hat{n}}}^{\kappa }+p{\tilde{v}}_{{\hat{n}}},1\right\rangle _{M}\,\mathrm {d}S. \end{aligned}$$
(D.32)
Appendix E: Entropy Conservative and Kinetic Energy Preserving Flux of Ranocha
Here, we state the entropy conservative and kinetic energy preserving (ECKEP) flux for the Euler equations recently developed by Ranocha [43]. The flux satisfies the entropy conservative spatial condition (2.38) as well as the KEP conditions of Jameson (3.31) on Cartesian meshes. The flux can be made entropy stable by adding a matrix dissipation of the form (2.55) where complete details are given in Gassner et al. [22]. We rewrite the form of the ECKEP originally presented by Ranocha [43] through the use of the identity
$$\begin{aligned} \frac{1}{4}\llbracket {a}\rrbracket \llbracket {b}\rrbracket = \left\{ \!\left\{ ab\right\} \!\right\} - \left\{ \!\left\{ a\right\} \!\right\} \left\{ \!\left\{ b\right\} \!\right\} . \end{aligned}$$
(E.1)
We do so such that the ECKEP flux only involves arithmetic and logarithmic means. The flux in the \(x-\)direction is given by
$$\begin{aligned} \mathbf{F}_1^{\text {ECKEP}} = \begin{bmatrix} \rho ^{\ln }\left\{ \!\left\{ v_1\right\} \!\right\} \\ \rho ^{\ln }\left\{ \!\left\{ v_1\right\} \!\right\} ^2 + \left\{ \!\left\{ p\right\} \!\right\} \\ \rho ^{\ln }\left\{ \!\left\{ v_1\right\} \!\right\} \left\{ \!\left\{ v_2\right\} \!\right\} \\ \rho ^{\ln }\left\{ \!\left\{ v_1\right\} \!\right\} \left\{ \!\left\{ v_3\right\} \!\right\} \\ \left( \rho ^{\ln }\left( \left\{ \!\left\{ v_1\right\} \!\right\} ^2{+}\left\{ \!\left\{ v_2\right\} \!\right\} ^2{+}\left\{ \!\left\{ v_3\right\} \!\right\} ^2{-}\frac{1}{2}\left( \left\{ \!\left\{ v_1^2\right\} \!\right\} {+}\left\{ \!\left\{ v_2^2\right\} \!\right\} {+}\left\{ \!\left\{ v_3^2\right\} \!\right\} \right) \right) {+}\frac{\rho ^{\ln }}{(\gamma -1)\left( \frac{\rho }{p}\right) ^{\ln }}\right) \left\{ \!\left\{ v_1\right\} \!\right\} {+} 2\left\{ \!\left\{ p\right\} \!\right\} \left\{ \!\left\{ v_1\right\} \!\right\} -\left\{ \!\left\{ p v_1\right\} \!\right\} \\ \end{bmatrix}. \end{aligned}$$
(E.2)
The ECKEP flux in the other Cartesian directions is then easily recovered by rotation.
