Preventing pressure oscillations does not fix local linear stability issues of entropy-based split-form high-order schemes

Recently, it was discovered that the entropy-conserving/dissipative high-order split-form discontinuous Galerkin discretizations have robustness issues when trying to solve the simple density wave propagation example for the compressible Euler equations. The issue is related to missing local linear stability, i.e. the stability of the discretization towards perturbations added to a stable base flow. This is strongly related to an anti-diffusion mechanism, that is inherent in entropy-conserving two-point fluxes, which are a key ingredient for the high-order discontinuous Galerkin extension. In this paper, we investigate if pressure equilibrium preservation is a remedy to these recently found local linear stability issues of entropy-conservative/dissipative high-order split-form discontinuous Galerkin methods for the compressible Euler equations. Pressure equilibrium preservation describes the property of a discretization to keep pressure and velocity constant for pure density wave propagation. We present the full theoretical derivation, analysis, and show corresponding numerical results to underline our findings. In addition, we characterize numerical fluxes for the Euler equations that are entropy-conservative, kinetic-energy-preserving, pressure-equilibrium-preserving, and have a density flux that does not depend on the pressure. The source code to reproduce all numerical experiments presented in this article is available online (DOI: 10.5281/zenodo.4054366).

Unsurprisingly, the choice of (symmetric) two-point flux function used in the novel volume term formulation is a key ingredient and determines the properties of the resulting high-order discretization.It is a somewhat surprising result that properties of the two-point fluxes used in simple low-order finite volume formulations directly translate to the high-order volume integral terms in this formulation.When using an entropy-conserving two-point finite volume flux, the corresponding two-point flux volume integral term of the DG scheme is entropy-conserving as well [10,27].The same is for instance true for kinetic energy preservation [14,34]; as we will show in this paper, it also holds for pressure equilibrium preservation.We note that the simple arithmetic mean two-point flux function recovers exactly the original nodal DG operator, while other choices of two-point flux functions may result in non-linear split-form DG operators, even for linear advection problems.
Unfortunately, it was recently discovered that the novel entropy-conserving/dissipative (and many other split-form) DG schemes can have stability issues [12].While the DG discretization is equipped with a provably discrete entropy inequality, it turns out that the schemes might struggle to retain local linear stability, i.e. the linear stability of the non-linear operator when linearized around a base-flow.Investigations of the spectrum of the linearized high-order operators revealed modes with spurious exponential growth, attributed to anti-diffusion of entropy-conserving twopoint fluxes.A particular striking example is given in [12] for the compressible Euler equations with a simple density wave where the density is variable, but the velocity and pressure are constant.Such a density wave (1) is a simple and smooth exact solution to the compressible Euler equations with perfect gas law, when equipped with appropriate (e.g.periodic) boundary conditions.Surprisingly, the entropy-conserving/dissipative DG schemes and other split-form variants fundamentally struggle for this simple problem.It turns out that the linearized spectrum shows spurious modes with exponential growth, that may cause fatal crashing of the simulation.
In another recent paper, Shima et al. [44] investigated the capability of their kinetic-energypreserving two-point flux to retain what they call pressure equilibrium.Consider the compressible Euler equations with an ideal gas law, where is the total energy, the internal energy, 2 /2 the kinetic energy, and Pressure equilibrium is precisely the case, when velocity and pressure are both constant, e.g. the density-wave (1).We get from the evolution of the compressible Euler equations the evolution equations of the velocity and of the pressure (5) It follows that for constant velocity and pressure, the time derivatives = 0 and = 0, hence the coined term pressure equilibrium.Shima et al. [44] found exponential spurious growth for a similar density-wave test case when using their kinetic-energy-preserving two-point flux [26].When they modified the two-point flux to discretely preserve pressure equilibrium, they could demonstrate numerically that the novel scheme robustly solves the density-wave, even for very long simulation times.
In summary, the starting point of this paper are the works [12,44] and we view the current work as a direct continuation of the analysis and discussion presented therein.This brings us directly to the research questions we are adressing in the current paper: (RQ1) Is it possible to construct two-point flux functions that are not only kinetic-energy-preserving and pressure-equilibrium-preserving as the one proposed by Shima et al. [44], but also entropy-conserving (EC) according to Tadmor's condition [50,51]?
(RQ2) Is pressure equilibrium preservation a remedy for the local linear stability issues of the entropy-conserving/dissipative DG framework reported in [12]?(RQ3) Are there entropies, such that the EC two-point fluxes and corresponding EC volume integral terms are locally linearly stable?
The remainder of the paper is organized as follows: in the next section, Section 2, we investigate research question (RQ1) and discuss the construction and existence of entropy-conserving (EC), kinetic-energy-preserving (KEP), and pressure-equilibrium-preserving (PEP) two-point flux functions.In Section 3, we investigate research questions (RQ2) & (RQ3) and discuss the impact of pressure equilibrium preservation on local linear stability.As a by-product, we show that the PEP property of the two-point flux function carries over to the high-order split-form DG scheme in the Appendix A. In the final Section 4, we summarize our results and collect the answers to the research questions.

Structure preservation properties
The first goal of this subsection is to collect and define the properties of the compressible Euler equations we want to preserve with our discretization.For the definition of two-point fluxes, it suffices to concentrate on semi-discrete finite volume methods of the form In what follows, we drop the subscript + for the numerical flux function for convenience and assume an interface at location and + 1 if not stated otherwise.
Definition 2.1 (Entropy-conservation [50,51]).A numerical flux num and the corresponding finite volume method is EC if where = are the entropy variables, is the flux potential, and [[ ]] := +1 − denotes the jump operator.
⊳ Unless stated otherwise, we will use the entropy of the compressible Euler equations ( 2), with associated entropy variables and flux potential = .
⊳ Definition 2.3 (Pressure equilibrium preservation).A numerical flux num = ( num , num , num ) and the corresponding finite volume method is PEP if whenever the velocity and the pressure are constant throughout the domain.⊳ We motivate our definition of PEP fluxes with the following Lemma 2.4.Pressure equilibrium, i.e. ≡ const, ≡ const, is preserved by (6) if and only if num is PEP.
Proof.The semidiscrete evolution equation for the velocity is Similarly, for = 0, the pressure evolves according to 1 Thus, = 0 and = 0 if and only if (11) is satisfied.
We are ready to formulate the central theorem of this work and to give the answer to the first research question (RQ1) in the following Theorem 2.5.The numerical flux of Ranocha [33,34], with logarithmic mean and product mean for the compressible Euler equations (2) is EC, KEP, PEP, and has a density flux num that does not depend on the pressure.Moreover, it is the only numerical flux with these properties for ≡ const.
Remark 2.6.The motivation for the last property, i.e. that the density flux does not depend on pressure such as e.g. in the EC flux by Ismail and Roe [21], is due to the discussion presented in [8,32], where positivity failure could be identified for certain setups with large pressure jumps and constant densities.⊳ Remark 2.7.The numerical flux (14) can also be derived by reversing the role of energy and entropy in the compressible Euler equations [32, Section 5].Indeed, the flux (66) of [32] is the same as (14) developed in [34,Theorem 7.8].This numerical flux is essentially uniquely defined by its properties, cf.Remark 2.13.⊳

Proof of Theorem 2.5
We first investigate the necessary conditions for EC and PEP and get the following Lemma 2.8.For ≡ const, ≡ const, an EC numerical flux that is also KEP or PEP satisfies Proof.For ≡ const, ≡ const, the left-hand side of (7) reduces to Inserting num = num + from the KEP (10) or PEP (11) property and using the discrete chain rule num .(20) This expression has to vanish for arbitrary values of ± for an EC flux, resulting in (17).Lemma 2.9.For ≡ const, ≡ const, an EC and PEP numerical flux must be of the form Proof.Comparing ( 11) and ( 17), − 1 {{ }} log num (22) must be independent of ± .Hence, num must be of the form num = {{ }} log for ≡ const, ≡ const.Inserting the PEP property (11) for num results in the final form (21).
Lemma 2.10.For ≡ const, an EC and PEP numerical flux for which the density flux does not depend on the pressure must be of the form where ( ± , ± ) is some kind of mean value depending on ± , ± such that ∀ ± , > 0 : Proof.Because of Lemma 2.9, the general form of dependencies on for ≡ const are already determined.The remaining degrees of freedom for non-constant pressure can be described by two functions 1,2 , resulting in the numerical fluxes where 1,2 depend on ± , ± such that Inserting this form of the numerical flux in the left-hand side of (7) for ≡ const results in Since this has to vanish for arbitrary ± , ± , , Having established the lemmata above, we are prepared to prove Theorem 2.5.
Proof of Theorem 2.5.The KEP (10) property is satisfied by construction.Moreover, the numerical flux for the total energy satisfies the PEP property (11), since it can be written as where 1 whenever is constant.Finally, the flux is EC as shown in [33,34].It is the only numerical flux with all these properties for ≡ const, since the KEP property (10) requires the pressure mean in (23) to be = {{ }} .
Remark 2.11.The pressure mean in the momentum flux is determined uniquely by the KEP property (10), resulting in a pressure mean depending on the density in the energy flux.As required by the PEP property (11), this dependency occurs only for non-constant pressure.However, such a mixed dependency on and of an approximation to the pressure is necessary for EC and PEP fluxes because of Lemma 2.10.⊳ We have obtained a complete characterization of numerical fluxes for the compressible Euler equations that are EC, KEP, PEP, and have a density flux num that does not depend on the pressure for ≡ const in Theorem 2.5.The analogous characterization for ≡ const is a bit more involved and leaves a degree of freedom.Lemma 2.12.For fixed ≡ const, a (symmetric) EC, KEP, and PEP numerical flux must be of the form where is a function depending on ± , ± (symmetrically with respect to ±) such that ∀ ± , : Proof.Because of consistency, every numerical flux can be written as the sum of a given numerical flux and a perturbation that is consistent with zero.Using ( 14) as baseline flux for fixed ≡ const, every numerical flux can be written as where ∀ , : ( , , , ) = ( , , , ) = ( , , , ) = 0. Kinetic energy preservation (10) requires = {{ }} .Since the chosen baseline numerical flux ( 14) is EC, requiring entropy conservation for the perturbed numerical flux yields Hence, Pressure equilibrium preservation (11) requires + const( , ) for ≡ const.The first term 1   2 (( • )) satisfies this requirement.However, the second term −1 {{ }} log fits if and only if ∀ ± , : ( + , − , , ) = 0. Remark 2.13.Extending the numerical fluxes (30) developed for fixed pressure ≡ const to general variable pressures results in a pressure-dependent density flux unless the perturbation vanishes, i.e. = 0. Thus, the numerical flux ( 14) is also unique for general velocities in the class of continuous numerical fluxes with the properties given in Theorem 2.5.⊳ Proof.Using the ansatz (31) for a general pressure yields num = {{ }} log {{ }} + ( ± , ± ), where ∀ , , : ( , , , ) = ( , , , , , ) = ( , , , , , ) = 0. Kinetic energy preservation (10) requires again = {{ }} .Requiring entropy conservation additionally yields For arbitrary ± , ± , choosing ± such that [[ / ]] = 0 requires = 0. Hence, the perturbation must vanish if the density flux does not depend on the pressure.

The KEP and PEP two-point flux of Shima et al.
Shima et al. [44] introduced a modification to their KEP flux [26] and constructed a KEP flux with the PEP property, We note that the density flux and the general structure of the momentum and energy fluxes is very closely related to Ranocha's two-point flux (14), except for the EC property, because Shima et al. use the arithmetic mean in the density flux instead of the logarithmic mean.Although the numerical flux (36) is not EC, it has four desirable properties, namely KEP, PEP, and a pressure-independent density flux.As we realize later in Section 3, the fourth desirable property is the arithmetic mean of the density in the density flux function, as it enhances robustness for the density wave propagation.
In their paper, Shima et al. demonstrate numerically very good robustness of their novel KEP and PEP discretization, even for highly non-linear problems such as underresolved turbulence.Hence, an interesting question is whether there is an entropy function for the compressible Euler equations such that the two-point flux function of Shima et al. with the arithmetic mean happens to be an EC flux.This would be a possible explanation of the improved numerical robustness of this flux for non-linear problems.To partially answer this question, we consider next the family of entropy functions introduced by Harten [15].
Harten [15] discovered the family of entropy functions for the Euler equations (2) where ℎ is a sufficiently smooth function satisfying to ensure convexity of the entropy function , equation (37).In particular, Harten discovered the one-parameter family Up to now, we considered the entropy (8) given by ℎ( ) ∝ above, since it is the only convex entropy (37) which symmetrizes the compressible Navier-Stokes equations with heat flux [19].Nevertheless, it is interesting to know whether there are other entropies (37) of the compressible Euler equations (2) that result in a corresponding EC numerical density flux num , where the mean value of the density is arithmetic.Following the approach used in Section 2.2, we will make use of the entropy variables and the flux potential associated with the entropy (37).Lemma 2.8 is a special case of Lemma 2.14.For ≡ const, ≡ const, an EC numerical flux for the entropy (37) that is also KEP or PEP satisfies Proof.For ≡ const, ≡ const, the left-hand side of (7) reduces to ]. (43)  This term vanishes if and only if ]. (44)   Inserting num = num + from the KEP (10) or PEP (11) property results in This expression has to vanish for arbitrary values of ± for an EC flux, resulting in (42).
Comparing (11) and ( 42), a PEP flux that is also EC for (37) must contain an average of the density proportional to In general, ( 46) is not the arithmetic mean of ± .It becomes the linear mean proportional to {{ }} for ≡ const for the choice of ℎ as in (39) with = −2 .However, in this case, the resulting is not convex anymore.Hence, entropy-conservative and pressure equilibrium preserving numerical fluxes for the compressible Euler equations have to use nonlinear means of the density and we have demonstrated that the two-point flux of Shima et al. is not related to one of Harten's entropies.

Corollary 2.15.
There is no Harten entropy pair for the compressible Euler equations such that a corresponding EC two-point flux with the KEP and PEP property uses the arithmetic mean of the density in the density flux.
So far, we have not found any evidence, that there is another strictly convex entropy pair for which the EC flux with KEP and PEP might have an arithmetic mean and thus have the conjecture, that there is none.
Finally, it is interesting to see whether the arithmetic mean can be used in the density flux of an EC flux if the additional constraints are relaxed by not requiring the KEP/PEP property anymore.Considering again the family of entropies (37), we consider the case ≡ const, ℎ ( ) / ≡ const.Inserting the entropy variables (40) into the EC condition (7) Hence, the density flux must again contain an average of the density proportional to (46), but for the case of ℎ ( ) ≡ const instead of the case ≡ const discussed above.For a given entropy such as (8) or (39), it is easy to solve ℎ ( + ) + / + = ℎ ( − ) − / − for + and pick values of ± , − such that (47) is not satisfied by num = {{ }} .Hence, we arrive at

Corollary 2.16. There is no Harten entropy pair for the compressible Euler equations such that a corresponding EC two-point flux uses the arithmetic mean of density in the density flux.
Remark 2.17.We refer to an entropy pair if the entropy function is strictly convex (resulting in an invertible transformation from the conserved variables to the entropy variables).Linear functionals of the conserved variables are of course non-strictly convex and can be combined with a density flux using the arithmetic mean.⊳

On local linear stability of EC schemes with the PEP property
In this section, we consider the (local) linear stability [12] of high-order discretizations based on two-point fluxes.The extension to high-order accuracy is relatively straight forward when assuming the SBP property.Several classes of numerical methods can be formulated via (periodic) SBP operators, including finite difference [25,48], finite volume [28,29], continuous Galerkin [17,18], discontinuous Galerkin [13], and flux reconstruction methods [39].A brief review how to formulate these methods in the SBP framework with application to structure-preserving numerical methods can be found in [38].Further details and background information about SBP methods can be found e.g. in the review articles [9,49].
Building upon earlier works such as [27,45], Fisher and Carpenter [10] created conservative high-order semi-discretizations of hyperbolic conservation laws using a special class of two-point numerical fluxes.The final extension to general symmetric numerical fluxes was obtained in [7,14,32] and will be recalled briefly below.

Numerical investigation of the robustness of the split-form DG scheme
In this part, we consider the split-form DG approximation with the numerical fluxes discussed above and apply them to solve a simple density wave problem.In particular, we compare the results when using a central flux with arithmetic means, Ranocha's two-point flux function (14) that is EC, KEP and PEP, and the two-point flux by Shima et al. (36) that is KEP and PEP.
Following [12], we consider the two-dimensional compressible Euler equations with the initial condition and fully periodic boundary conditions.We use the split-form DG methods with spectral collocation on Legendre-Gauss-Lobatto nodes, with a polynomial degree of = 5 on a grid with 4 × 4 elements implemented in the open source code Trixi.jl[42,43].The semi-discretizations are integrated in time using the fourth-order, five-stage, low-storage Runge-Kutta method of [23] with a relative CFL number cfl = 0.05, which then gets additionally scaled by the choice of polynomial degree , as is common for DG.This ensures a negligible impact of the time integrator on the numerical solution.
Remark 3.1.The numerical methods applied in this article are written in Julia [1].The plots are created using Matplotlib [20].The source code necessary to reproduce all results shown in this article is available online [36].We use the numerically stable evaluation of the logarithmic mean proposed in [21] in Trixi.jl.⊳ The simulation with the pure central approximation with the arithmetic mean flux is stable for all times for the density wave (48).We emphasize that we use the central flux for both, the splitform volume integral and the surface integral fluxes -so there is no added numerical dissipation.This shows that in principle, the problem is very well resolved by the chosen DG discretization.We can further numerically confirm, that the modification of Shima et al. [44] gives a high-order split-form DG discretization, that is able to robustly run this test case for very long integration times ( > 100, corresponding to more than 935 000 time steps).Hence, at first, it seems that the added PEP property indeed solves the robustness issue.However, in accordance with the findings of [12] for other EC fluxes, the high-order DG discretization with Ranocha's EC flux ( 14) crashes because of negative density already at ≈ 0.55.We note that Ranocha's flux is KEP and PEP, i.e., it preserves the pressure equilibrium by construction and there are no fluctuations in velocity and pressure!Following these numerical results, we can already answer our second research question and state, that, unfortunately, the answer to (RQ2) is no, the PEP property is not a remedy for the stability issues of the EC fluxes, as the simple density wave problem still crashes after very short simulation times.Note, that these findings are not sensitive to the choice of the cfl number or the time integration method.
However, it is interesting that the Shima

Stability for linear advection: Investigation of the impact of the choice of the mean value
In this subsection, we focus on how the choice of the mean values in a split-form approximation of the linear advection equation influences the stability.For this purpose, we consider periodic high-order SBP discretizations of the linear advection equation.We use ⊳ Given a SBP operator with symmetric mass matrix and a symmetric two-point numerical flux num for the hyperbolic conservation law is a conservative approximation of ( 50) with at least the same order of accuracy as the SBP operator [7,10,32].Moreover, the semi-discretization conserves the entropy of (50) if the numerical flux is entropy-conservative for that entropy [10].Following [12], we compute the spectrum of a semi-discretization (51) of the linear advection equation ( , ) + ( , ) = 0, (0, ) = 2 + 1.9 sin( ), (52) in a periodic domain ∈ [0, 2].The scheme ( 51) is implemented in Julia [1] and the Jacobian of the semi-discretization is computed via forward-mode automatic differentiation (AD) [40].Note that AD is not necessary if a linear numerical flux is used for this linear PDE.However, we are also interested in nonlinear numerical fluxes, involving e.g. the logarithmic mean.For nonlinear discretizations, it is not as trivial to compute the Jacobian and AD becomes a valuable tool.We note that if the numerical flux is chosen as the arithmetic mean, num = {{ }}, the semi- discretization is linear and skew-symmetric (with respect to the scalar product induced by the mass matrix ).This semi-discretization conserves the 2 entropy ( ) = 2 /2.Hence, as expected, all eigenvalues are purely imaginary in our numerical test.
In contrast, choosing the numerical flux as the logarithmic mean num = {{ }} log results in a nonlinear semi-discretization, which conserves the entropy ( ) = log − with entropy flux ( ) = log − .Indeed, the corresponding entropy variables are ( ) = ( ) = log( ) and the flux potential is ( ) = .Hence, the associated entropy-conservative numerical flux is The resulting spectra of the Jacobian of the semi-discretization with logarithmic mean values are shown in Figure 1 for different choices of SBP operators.Clearly, all of the spectra have eigenvalues with positive real part (of order unity) that do not converge to zero under grid refinement.In particular, these eigenvalues with positive real part occur for all choices of semi-discretizations.We checked that the occurrence of eigenvalues with positive real parts does not depend on the parity of the number of nodes/elements or the polynomial degree.For the investigation of spectra based on discretizations with other mean value choices, we refer to Appendix B. All other mean values tested give discretizations where the spectra have significant positive real parts.This nicely suggests that the reason for the stability issues of the high-order split-form DG scheme with the EC flux of Ranocha for the simple density wave example is due to the logarithmic mean of the density in the density flux.Even with the PEP property, which guarantees that pressure and velocity stay constant throughout the simulation, the discretization of the density evolution is unstable when using the logarithmic mean, while the discretization with the flux of Shima et al. is based on an arithmetic mean of the density and hence runs the example robustly.It remains to discuss however, if the Shima et al. flux is locally linearly stable as defined in [12], i.e., if the spectrum of the linearized operator is stable towards perturbations.

Investigation of local linear stability
We want to dig deeper and analyze the respective spectra of the Jacobians of the different DG semidiscretizations for the two-dimensional compressible Euler equations.We observed in Section 3.1 that the EC scheme immediately crashes whereas the central scheme and the scheme powered by the Shima et al. flux run for very long times ( > 100) without any problems.
For the linearization, we use the initial condition as the linearization state and compute the Jacobians approximately with a central finite difference approach in our simulation framework     As expected, the central flux results in a spectra that is almost purely imaginary, with only small deviations of the eigenvalues from the imaginary axis that are within machine accuracy when considering the approximation of the Jacobian via finite differences and the conditioning of the associated eigenvector matrix for this problem.Remark 3.3.We would like to stress that eigenvalues of nonlinear right-hand sides in an ODE ( ) = ( ( )) do not necessarily predict the global behavior of solutions.For example, an energyconserving ODE with purely positive and negative eigenvalues is discussed in [35,37].However, eigenvalues of the linearized Jacobian predict the local behavior of the solution, e.g. the temporal development of initial perturbations.In fluid dynamics, it is well known that there are many flow states that are physically unstable, i.e. flow states such as shear layers where added initial perturbations grow exponentially in time, until they start to behave non-linearly and transition to turbulence.In the considered case of density propagation discussed in this section, we do not expect to find significant signs of exponential growth, as also indicated by the almost imaginary spectrum obtained for the central flux.In particular, the compressible Euler equations are reduced to linear advection equations and hence this particular flow state is physically stable with respect to density perturbations.⊳ As anticipated following the numerical investigations in Section 3.1, the EC, KEP, and PEP flux (14) of Ranocha yields eigenvalues with clearly positive real parts of order unity, which do not vanish under grid refinement, but shift to higher imaginary values (see [12] for a more detailed discussion on this effect).This underlines our conclusion to research question (RQ2), that the PEP property does not fix the local linear stability issue of the EC split-form DG scheme.
Surprisingly, the spectrum with the KEP and PEP flux (36) of Shima et al. [44] shows similar issues -it clearly has eigenvalues with positive real parts.This discretization is not locally linearly stable neither.However, as observed in Section 3.1, the discretization could robustly handle the density wave example (we made sure to test very long times > 100 as well).
It is important to point out that the density wave example is a specific test case particular well suited to the Shima et al. powered discretization where it works perfectly fine, as it preserves the pressure and velocity as constants down to machine precision and hence reduces to the central scheme in this particular case.The spectrum clearly shows that adding just a small perturbation to this state may lead to spurious exponential growth -hence it is not locally linearly stable.The spectrum with the central flux is (almost) purely imaginary and has no growth of any modes.Constant pressure and velocity within machine precision accuracy cause very small perturbations in the range of 10 −15 .Hence, it would take a really (really) long simulation run time until these machine accuracy fluctuations grow.Furthermore, at this very small perturbation scale, the artificial dissipation of the time integration is effective and the cfl would have to be drastically reduced.
We thus need a feasible setup to further assess the robustness of the DG split-form with the Shima et al. flux: we make a simulation that investigates the growth of medium scale perturbations added to the initial conditions (see [12] for additional details).We start with the same setup as above and compute the eigenvector ˜ 0 associated with the biggest real eigenvalue of the semi-discretization using the numerical flux (36) of Shima et al. [44].The numerically computed eigenvector is purely real-valued.Note that this eigenvector can always be chosen to be real-valued since the Jacobian and its corresponding eigenvalue are both real-valued.We normalize the eigenvector such that ˜ 0 ∞ = 1 and use the perturbed initial condition 0 + 10 −3 ˜ 0 , where 0 is the original initial condition given by (48).Thus, the perturbation scale is now 10 −3 , instead of 10 −15 .Further, the shape of the perturbation corresponds to the eigenmode of the spectrum.Thus, we are able to compare the growth of the fluctuations from the simulation, to the one predicted by the spectra using the real part of the eigenvalue as the growth rate.
To get the evolution of the perturbation, we subtract in each Runge-Kutta stage the semidiscretization applied to the unperturbed initial condition from the resulting semi-discretization of the perturbed initial state.We perform these numerical experiments using the flux (36) of Shima et al. [44] as surface flux for the DG scheme, and in addition also using the dissipative HLL flux [16] as surface flux, while both discretizations use the Shima et al. flux for the split-form volume integral.
The resulting discrete ∞ error of the perturbations in the conserved variables is visualized in Figure 3. Clearly, the fluctuations grow exponentially with a rate perfectly matching that of the real part of the eigenvalue.The simulation terminates at ≈ 4.6 because of negative densities for the case without surface dissipation and at ≈ 9 if the HLL flux is used.Clearly, the dissipative HLL flux reduces the growth of the fluctuations but only quantitatively, not qualitatively.Surface dissipation can not guarantee to control errors stemming from badly discretized split-form volume integrals.Consequently, in this case, the error still spuriously grows exponentially in time and finally results in unphysical solutions, which underlines the local stability issues of the Shima et al. flux and hence confirms the statement in [12], that many split-form discretizations have these problems.
We also considered a random perturbation of the initial condition, where each conserved variable is perturbed randomly at each point with a uniform distribution that is symmetric around zero.Such an approach is also used to estimate the Lyapunov exponent of a dynamical system [52].The resulting errors of the fluctuation simulation are visualized in Figure 4. Without adding   dissipation, the error grows approximately exponentially with a rate governed by the maximal real eigenvalue.When surface dissipation in form of an HLL flux is added, the perturbation grows slowly at first but shows the same exponential growth governed by the largest real eigenvalue later.

EC and local linear stability
Combining the results from our numerical investigations with Corollary 2.15 or Corollary 2.16 gives at least a partial answer to our third research question (RQ3): There are no Harten entropies for the compressible Euler equations such that the associated EC two-point fluxes result in locally linearly stable schemes.

A. The PEP property for high-order schemes
Extending Lemma 2.4, the PEP property (11) of a numerical flux num extends directly to a highorder semi-discretization (51).
Lemma A.1.A pressure equilibrium ≡ const, ≡ const is preserved by any general linear method applied to the semi-discretization (51) if the numerical flux num is PEP.
Proof.It suffices to consider linear combinations of numerical solutions as well as the addition of the semidiscrete operator to a numerical solution.Linear combinations preserve a pressure equilibrium, since , , and (53) Here, we used , = 0, which is a necessary condition for a consistent derivative operator .
This list is ordered in descending order of the size of the mean values [6].
The resulting spectra are shown in Figure 5. Clearly, all mean values except the arithmetic mean result in eigenvalues with positive real parts.The size of the maximal real part of the spectrum increases for mean values that deviate more from the arithmetic mean value.
et al. flux can indeed robustly run the density wave test case.Because of the PEP property, both the flux (36) of Shima et al. and the flux (14) of Ranocha reduce the density wave (1) for the compressible Euler equations to four linear advection equations.The main difference between the Shima et al. flux and Ranocha's flux is the mean values of the density used in the density flux: The former uses the arithmetic mean to approximate the linear advection, the latter uses the logarithmic mean to discretize the linear advection.
Finite difference methods of order ∈ {2, 4} using different numbers of nodes.Continuous Galerkin methods using polynomials of degree ∈ {3, 4} on Lobatto Legendre bases and different numbers of elements.
Discontinuous Galerkin methods using polynomials of degree ∈ {3, 4} on Lobatto Legendre bases and different numbers of elements.

Figure 1 .
Figure 1.: Spectra of nonlinear semi-discretizations (51) using different variants of SBP operators and the logarithmic mean as flux.
Trixi.jl[42].The resulting spectra for the central flux with arithmetic means, the flux by Shima et al., and Ranocha's EC flux are shown in Figure2.As in Section 3.1, we use the same numerical flux for the volume terms and the surface terms without any further dissipation.

Figure 2 .
Figure 2.: Spectra of split-form DG semi-discretizations of the compressible Euler equations in two space dimensions with linearization state(48).The DG methods use polynomials of degree = 5 and a uniform grid of 4 × 4 elements in the domain [−1, 1] 2 with periodic boundary conditions.The same numerical flux is used for both the volume terms and at surfaces.

Figure 3 .
Figure 3.: Evolution of an eigenvector perturbation of the initial condition (48).The DG methods use polynomials of degree = 5 and a uniform grid of 4 × 4 elements in the domain [−1, 1] 2 with periodic boundary conditions and the numerical flux of Shima et al. [44].

Figure 4 .
Figure 4.: Evolution of a random perturbation of the initial condition (48).The DG methods use polynomials of degree = 5 and a uniform grid of 4×4 elements in the domain [−1, 1] 2 with periodic boundary conditions and the numerical flux of Shima et al. [44].

= 2 / 2 +
/( − 1) are linear in .Given ∈ R, the scaled addition of the semidiscrete operator to a solution preserves the pressure equilibrium, since :