Von Neumann Stability Analysis of DG-Like and PNPM-Like Schemes for PDEs with Globally Curl-Preserving Evolution of Vector Fields

This paper examines a class of involution-constrained PDEs where some part of the PDE system evolves a vector field whose curl remains zero or grows in proportion to specified source terms. Such PDEs are referred to as curl-free or curl-preserving, respectively. They arise very frequently in equations for hyperelasticity and compressible multiphase flow, in certain formulations of general relativity and in the numerical solution of Schrödinger’s equation. Experience has shown that if nothing special is done to account for the curl-preserving vector field, it can blow up in a finite amount of simulation time. In this paper, we catalogue a class of DG-like schemes for such PDEs. To retain the globally curl-free or curl-preserving constraints, the components of the vector field, as well as their higher moments, must be collocated at the edges of the mesh. They are updated using potentials collocated at the vertices of the mesh. The resulting schemes: (i) do not blow up even after very long integration times, (ii) do not need any special cleaning treatment, (iii) can operate with large explicit timesteps, (iv) do not require the solution of an elliptic system and (v) can be extended to higher orders using DG-like methods. The methods rely on a special curl-preserving reconstruction and they also rely on multidimensional upwinding. The Galerkin projection, highly crucial to the design of a DG method, is now conducted at the edges of the mesh and yields a weak form update that uses potentials obtained at the vertices of the mesh with the help of a multidimensional Riemann solver. A von Neumann stability analysis of the curl-preserving methods is conducted and the limiting CFL numbers of this entire family of methods are catalogued in this work. The stability analysis confirms that with the increasing order of accuracy, our novel curl-free methods have superlative phase accuracy while substantially reducing dissipation. We also show that PNPM-like methods, which only evolve the lower moments while reconstructing the higher moments, retain much of the excellent wave propagation characteristics of the DG-like methods while offering a much larger CFL number and lower computational complexity. The quadratic energy preservation of these methods is also shown to be excellent, especially at higher orders. The methods are also shown to be curl-preserving over long integration times.


Introduction
Novel PDEs, important to science and engineering problems are routinely discovered. Many of those PDE systems have involution constraints whose deeper study leads to mimetic numerical solution strategies that retain greater fidelity with the physics of those systems. While involution constraints can come in many forms, the last couple of years have seen the emergence of novel classes of PDEs that place constraints on the evolution of the curl of a vector field. The structure of the evolutionary equation for the vector field of interest is such that it either keeps the curl of the vector field zero, or evolves it in proportion to certain source terms. The former PDEs are referred to as curl-free whereas the latter class of PDEs are referred to as curl-preserving. Many of the hyperbolic systems resulting from the Godunov-Peshkov-Romenski (GPR) formulation for hyperelasticity and compressible multiphase flow with and without surface tension have such curl-preserving update equations (Godunov and Romenski [32], Romenski [40], Romenski et al. [41], Peshkov and Romenski [37,38], Dumbser et al. [27,30,31], Schmidmayer et al. [42]). The equations of general relativity, when cast in the FO-CCZ4 formulation, also have such a structure (Alic et al. [1,2], Brown et al. [18], Dumbser et al. [29], Dumbser et al. [28]). Similarly, it has recently become possible to recast Schrödinger's equation in first-order hyperbolic form, and the time-evolution of this very important equation also has curl-preserving constraints (Dhaouadi et al. [25], Busto et al. [19]). A motivating PDE system would help the reader a lot here. Let us take a simple example involving a fluid with thermal conduction in the GPR formulation. Let us denote the density by , the fluid velocity by , the fluid pressure by P , the fluid temperature by T , the internal thermal energy density by e , the total energy density by E ≡ e + 2 2 , the thermal impulse by a vector , the heat flux by a vector , and the thermal stress by the second rank tensor . The equations for a fluid with thermal conduction can be written as Here is the vector field of interest and is some specified velocity field. By " " we generically denote other variables that might be part of a larger set of equations to which Eq. (2) belongs. As a result, the potential ≡ ⋅ + ( ) can also depend on other variables and the source term ( , ) depends on " " as well as " ". By setting the source term " " to zero and taking the curl of the above equation, it is easy to see that if we initially have (∇ × ) = 0 , then the structure of the above equation retains (∇ × ) = 0 for all time; in other words, the evolution of the vector field is curl-free. If the above equation has a nonzero source term, the evolution of the vector field would be curl-preserving. In principle, any curl-preserving scheme should be able to reach the curl-free limit as the importance of the source term goes to zero. For the purpose of this study, it is adequate to study curl-free evolution of vector fields because the inclusion or exclusion of the source term does not change the discussion that follows. In other words, the inclusion of source terms requires a separate study of how stiff terms are included in the time-evolution of any hyperbolic PDE, and that is not the point of focus for this paper. Equation (2) shows up quite frequently as a part of many of the hyperbolic systems mentioned before. When Eq.
(2) appears in larger systems, it is usually the source of many numerical difficulties. Equation (2) and the larger systems that include it cannot be evolved with a traditional higher-order Godunov methodology, even when we have accounted for non-conservative products. Indeed, Dumbser et al. [28] have shown that if a classical higher-order Godunov scheme is applied to Eq. (2), the solution blows up rapidly. The blow-up manifests itself as an explosive increase in the curl of the vector field when the simulation is run over long periods of time; please see Fig. 5 of Dumbser et al. [28] which shows the explosive blow-up of the curl of the vector field. The blow-up occurs even when the source terms are zero, indicating that the source terms are not the cause of the difficulty. The problem is not specific to the PDE system considered in that paper because as shown in Boscheri et al. [17], a simple model system based on Eq. (2) will also blow up if treated with a classical higher-order Godunov scheme. A GLM-based cleaning procedure has been developed in Dumbser et al. [28] for suppressing the fictitious numerical build up of the circulation, but it requires greatly increasing the signal speed that is used in the cleaning equations and also adds many more vector fields than necessary. Indeed, Dumbser et al. [28] in their study of a general relativistic system had to use signal speeds in their GLM-style approach that were substantially larger than the speed of light. It is unusual to introduce superluminal signal speeds in a theory that is based on recognizing the speed of light as the maximal signal speed. As a result, GLM cleaning-based schemes force a very strong reduction in timestep, inconsistent with our notions of how the timestep should evolve in a higher-order Godunov scheme. Boscheri et al. [17] obtained curl constraint-preservation but only at the expense of solving an elliptic PDE system at every timestep. Furthermore, the method was restricted to second order of accuracy. This too is inconsistent with our notion that a higher-order Godunov scheme should be time-explicit, easy to solve, and extensible to all higher orders.
In Balsara et al. [12] we first realized that the constrained evolution mandated by Eq. (2) requires a special treatment. Two innovations were introduced in that paper. First, a curlconstraint preserving reconstruction was proposed which requires us to start with a vector field whose components are collocated at the edges of the mesh, and are indeed aligned with the edge directions. Second, it was realized that a multidimensional Riemann solver that is invoked at the vertices of the mesh can indeed give us stabilization through multidimensional upwinding. Coupled with strong stability preserving Runge-Kutta (SSP-RK) (2) t + ∇( ⋅ + ( )) − × (∇ × ) = ( , ).
timestepping (Shu and Osher [44,45], Shu [43], Spiteri and Ruuth [46,47], Gottlieb et al. [33]), we obtained a robust finite volume-based numerical scheme. As a result, both ingredients in the design of a higher-order Godunov scheme-i.e., the reconstruction as well as the Riemann solver-had to be fundamentally rethought when dealing with Eq. (2). These two ingredients proved highly beneficial because Balsara et al. [12] were able to obtain stable finite volume-like schemes for systems that used Eq. (2) which: (i) did not blow up even after very long integration times, (ii) did not need any GLM-style cleaning with its deleterious side-effect of needing very high signal speeds, (iii) could operate with large explicit timesteps, (iv) did not require the solution of an elliptic system and (v) could be extended to higher orders using WENO-like methods. The WENO-like methods draw on ideas from weighted essentially non-oscillatory schemes (Jiang and Shu [35], Balsara and Shu [15], Balsara et al. [9]). As a result, Balsara et al. [12] were able to integrate non-linear hybridization into their novel curl-preserving reconstruction, thus making the curl-preserving reconstruction suitable for use with higher-order Godunov schemes. In Balsara et al. [12] curl-preserving WENO-like methods were developed for evolving systems that used Eq. (2). We called such methods WENO-like because the reconstruction used many insights from WENO schemes, while being substantially different from traditional, finite volume-based WENO schemes. Recall that a conservation law has a flux form ensuring that the integrated conserved variables in any zone (or connected set of zones) evolve in response to the fluxes at the boundary of that volume. This telescoping property for the fluxes gives the discrete version of the conservation law a globally conservative property. In an entirely analogous fashion, the curl-free schemes presented in Balsara et al. [12] are such that the discrete circulation evaluated over the edges of any face (or collection of faces) depends only on the potentials at the vertices of that facial area. In that sense, the schemes developed in Balsara et al. [12] are globally curl-preserving because they have a telescoping property on the potentials. In that paper, we also presented a von Neumann stability analysis of WENO-like globally curl-preserving schemes, showing that with the increasing order of accuracy the schemes became progressively less dissipative and their dispersion error was also reduced. But it is well-known (Reed and Hill [39], Cockburn and Shu [22][23][24], Cockburn et al. [20,21], Liu et al. [36], Zhang and Shu [48]) that finite volume DG schemes have superior wave propagation properties relative to finite volume WENO schemes of comparable order. In their study of DG-like schemes for magnetohydrodynamics and Maxwell's equations that have one or more constrained, divergence-preserving vector fields, Balsara and Käppeli [10,11] found a similar trend, and that trend was numerically confirmed by Hazra et al. [34] and Balsara et al. [13]. We may, therefore, expect that curl-preserving DG-like schemes for Eq. (2) should have wave propagation characteristics that are substantially better than curl-preserving WENO-like schemes for Eq. (2) at comparable orders of accuracy. This paper, therefore, has the following three goals.
i) The first goal of this paper is to lay out the conceptual foundations for DG-like schemes that preserve the global curl constraint. Recall that a classical DG scheme starts with zone-centered mean values for the flow and endows it with higher moments. These higher moments are then evolved in time by a classical DG scheme. The time evolution of the higher moments is performed consistently with the governing equations. In an exactly analogous fashion, the globally curl-free DG-like schemes start with edgecentered components. These components are then endowed with higher moments whose time-evolution is carried out consistently with the governing equations.

3
ii) The second goal of this paper is to conduct a von Neumann stability analysis of the newly-obtained curl-free DG-like schemes. We use this stability analysis to find the maximum CFL number that is available at all orders up to the fourth order. We also use the stability analysis to show that our new class of DG-like schemes have superior wave propagation characteristics. The curl of a vector field only manifests itself in two or three dimensions. As a result, our von Neumann stability analysis is also two-dimensional. Because the stability analysis is multidimensional by necessity, it is not technically feasible to extend it to very high orders. The computer algebra systems that we use for this stability analysis cannot be pushed beyond fourth order. iii) Dumbser et al. [26] proposed PNPM schemes where all the modes that are up to Nth degree were evolved while all the higher modes up to Mth degree are reconstructed. The PNPM schemes had the great advantage that they displayed wave propagation properties almost as good as classical DG schemes of comparable order while permitting substantially larger timesteps. Balsara and Käppeli [10,11] found a similar trend in their study of divergence-preserving PNPM-like schemes. Therefore, the third goal of this paper is to analyze PNPM-like schemes that are globally curl-preserving. We intend to show that such PNPM-like schemes are competitive with their DG-like counterparts at comparable order of accuracy. We also intend to show that the maximum CFL number of PNPM-like schemes is much larger than that of DG-like schemes.
In this paper we analyze Eq. (2) with ( ) = 0 and ( , ) = 0 . A constant velocity " " is specified. The plan of the paper is as follows. Section 2 introduces a DG-like formulation for a curl-preserving model equation. Section 3 provides the essential ideas behind curlfree and curl-preserving reconstruction while pointing the reader to further literature. Section 4 presents the von Neumann stability analysis. Section 5 presents results from the von Neumann analysis of globally curl-free DG-like schemes. Section 6 presents some numerical results to show that the proposed schemes meet their order of accuracy. Section 7 presents conclusions.

DG-Like Formulation for the Curl-Preserving Model Equation
Let us consider Eq. (2) to understand its physics, thereby leading us to a curl-preserving DG-like formulation for that equation. Let us first consider Eq. (2) in its curl-free form (i.e., with ∇ × = 0 ) since the curl has to be kept mathematically zero in that limit. Since we are considering a minimalist system, we can also take ( ) = 0 because we are ignoring any further PDEs in the system. We write the vector components in two-dimensions as = (J x , J y ) and assume a constant velocity = v x , v y , so that the dot product becomes The equations that we have to solve can be written as Using the constraint in each of the two equations, they can be written as In other words, Eq. (4) tells us that we should be able to advect two components of a vector field. However, importantly, we should be able to carry out this advection of the vector field in a fashion that preserves the curl-free aspect of the vector field over each zone. We can only carry out such an advection if the component J x is collocated at the x-edges of the mesh and the component J y is collocated at the y-edges of the mesh. Furthermore, the potential, v x J x + v y J y , should be collocated at the vertices of the mesh, as shown in Fig. 1. Furthermore, Eq. (4) shows us that the potential should be obtained via multidimensional upwinding at the vertices of the mesh. For this simple example, multidimensional upwinding can be enforced visually depending on the direction of the velocity field. For example, Fig. 1 shows the multidimensionally upwinded potentials, when both components of the velocity are positive. On a mesh with zone sizes Δx and Δy in the x-and y-directions, the curl-free update equations at the first order (with positive velocity components) can be written as Fig. 1 The components of the curl-free vector field around the four zones centered around (i, j), (i − 1, j), (i − 1, j − 1), and (i, j − 1). A first order curl-free reconstruction is used in this figure; though the use of a higher order reconstruction of the vector field is also possible. The multidimensionally upwinded potentials at the vertices of the zone (i, j) are also shown for the situation where both components of the velocity are positive. Again, while the potentials are shown explicitly for the first order case in the figure, a higher order reconstruction will yield more accurate potentials at the vertices of the mesh Equation (5) shows us that the collocation described in Fig. 1 is crucial for maintaining a curl-free update and it will have to be built into our curl-free DG-like scheme. The use of identical potentials at each vertex of the mesh allows us to claim that the discrete version of the curl-free constraint is preserved forever for each and every zone of the mesh. In other words, the update equations are globally curl-free. Equation (5) also shows us the importance of multidimensional upwinding providing a unique potential at each vertex of the mesh. To prove that Eq. (6) follows from Eq. (5) at all orders of accuracy, please write out equations like Eq. Fig. 1. Please write the equations out in terms of the potentials i+1∕2,j+1∕2 , i−1∕2,j+1∕2 , i+1∕2,j−1∕2 , and i−1∕2,j−1∕2 shown at the vertices in Fig. 1. Notice that these potentials can be made as accurate as we desire by the use of a higher-order reconstruction of the vector field. The pairwise cancellation will show that Eq. (6) is satisfied at all orders-in other words, we have a mimetic scheme. More generally, when Eq. (2) is incorporated into a larger PDE system, the multidimensional Riemann solver (Balsara [3][4][5][6], Balsara et al. [8], Balsara and Dumbser [7], Balsara et al. [16], Balsara and Nkonga [14]) provides us with multidimensional upwinding consistent with the waves propagating in all directions.
The process of designing a higher-order DG-like scheme consists of endowing the primal variables with higher-order moments and then evolving those moments consistently with the governing equations. Consider the zone (i, j) in Fig. 1. We center the coordinates at the zone center so that the two-dimensional zone has extent −Δx∕2, Δx∕2 × −Δy∕2, Δy∕2 . At the right y-edge and the top x-edge of the zone we assert the moments up to fourth order as In a DG-like scheme, the modes in Eq. (7) become time-evolutionary. If only the linear part of Eq. (7) is retained, we get a second-order DG-like scheme; if the quadratic part of Eq. (7) is retained, we get a third-order DG-like scheme; if all the terms of Eq. (7) are retained, we get a fourth-order DG-like scheme. Equation (7) makes it easy to see the trial functions that are asserted in each of the mesh edges. We will use test functions identical to the trial functions. Let us first write the Galerkin projection in the abstract and then specialize it to Eq. (7). It is useful to realize that a traditional finite volume DG method is derived from a Gauss' law-based vector identity From the one-dimensional gradients involved in Eq. (3) we realize that the DG-like methods that we seek depend on the product rule for derivatives applied one-dimensionally. In other words, we rely on the identity In the above equation, " " is a test function at the mesh edges. Applying the above identity to the y-component of Eq.
(2) we can make the Galerkin projection in the y-edge of the mesh as follows: Here (y) is a test function at the y-edge of Fig. 1. Similarly, applying the above identity to the x-component of Eq. (2) we get There is not much going on in Eqs. (8) and (9) other than an integration by parts along with a one-dimensional Galerkin projection. However, in the next paragraph we interpret the above two equations to bring out the physics of the situation.
Equation (8) allows us to write the update equations for the evolutionary modes of J y (y, t) as

3
The angled brackets, 〈〉, represent line integrated averages of sufficiently high order within a y-edge (and later, similarly, for the x-edge). The potentials * * with the double star superscripts denote the potentials obtained at the mesh vertices using a multidimensional Riemann solver. The potentials * (y) with the single star superscripts denote potentials obtained by the application of one-dimensional Riemann solvers at the y-edge of the desired zone. These one-dimensional Riemann solvers may be invoked at multiple quadrature points at the y-edge so that the terms ⟨ * (y)⟩ and are accurately evaluated. The source terms have a similar interpretation so that the * and * variables in S y ( * , * ) are obtained from one-dimensional Riemann solvers. We see from Eq. (10a) that the update of the mean value, J y 0 (t) , will be curl-preserving and approach curl-free evolution in the limit where the source term tends to zero. Equations (10b)-(10d) show the same type of body terms that arise in a classical DG scheme due to the integration by parts with a test function with the key difference that they are now applied to the edges of the mesh. The terms involving (∇ × ) z in Eq. (10) also have a special interpretation. These terms are exactly zero when the evolution is curl-free. When the evolution is only curl-preserving, these terms will be proportional to the discrete circulation around the zone, but only if a curl-preserving reconstruction from Balsara et al. [16] is used. From each of the two sides of a two-dimensional mesh, we can obtain terms that provide (∇ × ) z . The resulting is therefore an arithmetic average of the curl evaluated from either side of that edge. We make analogous interpretations for the terms with (∇ × ) z in Eqs. (10b)-(10d).
Equation (9) allows us to write the update equations for the evolutionary modes of The interpretation of the terms in Eq. (11) mirrors that of Eq. (10). In Eq. (11), the angled brackets, 〈〉, represent line integrated averages of sufficiently high order within the x-edge.
In this work, we are interested in analyzing curl-free evolution, with the result that all terms with (∇ × ) z , S x ( * , * ) and S y ( * , * ) can be set to zero in Eqs. (10) and (11). Without the support of a larger PDE system, it is not possible to specify the source terms. Nevertheless, it is important for the reader to understand what a curl-preserving reconstruction is. We illustrate the same for the simplest of cases in the next section.

Curl-Preserving Reconstruction
Equations (10) and (11) show that we need a reconstruction strategy within a zone that matches the vector components and their higher moments in the edges of the mesh. This is needed because we want the scheme to be globally curl-preserving. Consequently, each edge, as seen by its abutting zones, will have the same component of the vector field as well as its higher moments. Equation (2), as well as Eqs. (10) and (11) show that when the discrete circulation evaluated around a zone is small, it should make proportionately small contributions to the edges via the (∇ × ) z -dependent terms. In other words, the curl of the reconstructed vector field should match the discrete circulation (evaluated around each zone) as well as its higher moments. Figure 2 shows us an example of how this works in two dimensions and at second order. In Balsara et al. [12] we have presented two and threedimensional versions of such a reconstruction strategy at several orders. Here we just show some details for the second-order case so that the reader may appreciate the core ideas as they are presented in one self-contained place.
Consider the vector field that is shown in Fig. 2. We consider the zone to be a unit square spanning −1∕2, 1∕2 × −1∕2, 1∕2 . The modes at the left and right y-edges are given by The modes in the bottom and top x-edges are given by We can write the equations of a vector field which matches the values of the components and their linear variation within each edge as follows: (12) J 1 y + Δ y J 1 y y and J 2 y + Δ y J 2 y y.
Note that the coefficients a yy and b xx are needed for curl constraint satisfaction. The discrete circulation within the zone of interest is given by J 1 with the result that we can obtain the higher modes of that variable up to linear variation and write it as y is zero, its variation will also be zero, so that we get Δ x R z and Δ y R z as zero values. If the discrete circulation y is small, its variation as represented by Δ x R z and Δ y R z will also be proportionately small. We want to fix the coefficients a yy and b xx so that the curl of the vector field in Eq. (14) exactly matches Eq. (15). This is obtained by setting This shows that the reconstructed vector field in Eq. (14) has been reconstructed in curl-preserving fashion. Therefore, the growth of the curl in Eqs. (10) and (11) is perfectly well-controlled and consistent with the PDE in Eq. (2). Furthermore, when the discrete circulation is exactly zero, the (∇ × ) z -dependent terms in Eqs. (10) and (11) contribute

Fig. 2
Collocation of vector components along the edges of a two-dimensional control volume. As evaluated over the edges of the square element, the discrete circulation is fully specified. The mean value of the vector components and their linear variation are shown along each edge, in keeping with a second order accurate reconstruction scheme. The reconstruction problem for a curl-preserving reconstruction consists of obtaining a polynomial-based vector field that matches the specified mean circulation in the zone while simultaneously matching the edge values within each zone absolutely nothing. In other words, the curl-free limit is exactly retrieved by our choice of reconstruction and discretization. We can use Eq. (14) along with Eq. (16) to write the reconstructed curl-preserving vector field in terms of an orthonormal set of modes that span the zone itself. Projecting Eq. (14) into an orthonormal basis made of tensor product Legendre polynomials, we get Notice that Eq. (17) shows us that there is a transcription from the modes that we use in a curl-preserving DG-like scheme to the modes that we would use for a finite volume-based DG scheme. We see that the modes of a curl-preserving DG-like scheme just combine differently, consistent with the constraints, to give us the modes in a traditional, finite volumebased DG scheme. This ensures that the order property is always retained. Note though that curl constraint-preservation results in some higher-order modes in Eq. (17) that would not be present in a classical second-order finite volume-based DG scheme. Therefore, the reverse transcription, i.e., going from the modes of a finite volume-based DG scheme of a certain order to the modes of a curl-preserving DG-like scheme of the same order, does not hold. At second order, one cannot make much from this transcription because all the coefficients in Eq. (17) are fully determined. However, as shown in Section II.4 of Balsara et al. [12], at fourth order and beyond, some of the modes of the higher-order curl-preserving reconstruction have to be obtained volumetrically while others are obtained from the edges of the mesh.

Von Neumann Stability Analysis of Curl-Free DG Schemes: Second-Order Example
The von Neumann stability analysis of a DG scheme can be performed in two different styles. The first is to convert the DG equations into a finite-difference-like form (Liu et al. [36], Zhang and Shu [48], Balsara and Käppeli [10]). The second approach is to identify the minimal number of modes, endow them with harmonic variation and then to directly carry out the stability analysis on the primal variables of the DG scheme (Balsara and Käppeli [11]). The latter approach works very well because it quickly allows us to identify the smallest number of variables that should be retained in a constraint-preserving DG scheme.
Here we describe the basic ingredients that go into carrying out a von Neumann stability analysis for a second-order accurate curl-free DG-like scheme. Please focus on Fig. 3. In the right y-edge of the central zone we identify the modes J y+ 0 (t) and J y+ y (t) as the mean y-component of the vector field and its linear variation in the y-direction. Similarly, in the top x-edge of the central zone we identify the modes J x+ 0 (t) and J x+ x (t) as the mean x-component of the vector field and its linear variation in the x-direction. In the spirit of a DG scheme, the modes are endowed with time-dependence. In the spirit of a harmonic variation, we assume rectangular zones of sizes Δx and Δy so that the Fourier modes vary as e −i(k x x+k y y) where the wave vector is given by k x , k y . In fact, we simplify even further by assuming square zones in most parts of this paper. Because we use Fourier modes in a von (17) Neumann stability analysis, the modes along the left y-edge are related to the modes along the right y-edge. Similarly, the modes along the bottom x-edge are related to the modes along the top x-edge. The relationship goes as follows Figure 3 shows even further inter-relationships among the modes that reside along the edges once the Fourier modal variation is assumed. At first blush it would seem that each zone in Fig. 3 has four independent pieces of information given by J y+ . However, because of the curl-free constraint, there are only three independent pieces of information. This becomes apparent when we use Eq. (18) to write the discrete curl-free condition in zone (i, j) of Fig. 3 as follows: As a result of Eq. (19), the only independent variables in the von Neumann stability analysis of a second-order accurate, curl-free DG-like scheme are J y+ 0 (t) , J y+ y (t) and J x+ x (t) . This simplifies the analysis considerably. Figure 3 shows how the facial modes at all the mesh faces are inter-related because of the Fourier modes and their spatial variation. As a result, the curl-free reconstruction of the vector field in each zone of Fig. 3 can be symbolically carried out using a computer algebra system. Equations (10a), (10b), and (11b) (in their curl-free forms) can then again be symbolically expressed using the same computer algebra system. As a result, we obtain expressions for the time rate of change of J y+ 0 (t) , J y+ y (t) , and J x+ x (t) that can be written in terms of J y+ 0 (t) , J y+ y (t) , and J x+ x (t) . In other words, we have reduced the problem of evaluating a single stage in the multistage RK-timestepping to the problem of obtaining a linear system of ODEs that look as follows: The nine coefficients in the matrix shown in Eq. (20) depend only on the wave numbers k x and k y , the velocities v x and v y , and the zone sizes Δx and Δy . They are explicitly given in Appendix A. We can also formally define the vector of unknowns as As a result, Eq. (20) can be formally written as where "A" is the 3 × 3 matrix shown in Eq. (20). We then discretize Eq. (20) in time with an explicit m-stage Runge-Kutta scheme having a timestep Δt of the form Here "I" is the identity matrix. The expressions for the coefficients i,k and i,k can be found in Gottlieb et al. [33] and also Spiteri and Ruuth [46,47]. Given the linearity of our DG scheme, we can write the time update as Here "G" is known as the amplification matrix of the scheme. It depends on the coefficients of the Runge-Kutta scheme, on the timestep Δt and the matrix "A" from Eq. (20). For the second-order SSP-RK scheme we can write the amplification matrix as Likewise, for the third-order SSP-RK scheme we can write the amplification matrix as  (24) = + Δt + Δt 2 2 2 + Δt 3 3! 3 .
In the next section, we will use this amplification matrix to devise our von Neumann stability analysis. This completes our description of the mathematics associated with the von Neumann stability analysis in second order. Higher orders can be done similarly.

Results from the von Neumann Stability Analysis of Globally Curl-Free DG-Like Schemes
By taking a close look at Eq. (3) we realize that the von Neumann stability analysis depends on the angle that the velocity vector v x , v y makes with respect to the x-axis of the mesh. Furthermore, Eqs. (18) and (19) show us that the von Neumann stability analysis also depends on the angle that the wave vector k x , k y makes with respect to the velocity vector v x , v y . For this reason, the stability analysis depends on multiple parameters. Besides, owing to the fact that the curl only manifests itself in two or more dimensions, it has to be multidimensional. For all of these reasons, we have only been able to carry out a von Neumann stability analysis for curl-free WENO-like, PNPM-like and DG-like schemes up to fourth order of accuracy. However, we realize that such a stability analysis that is done in two dimensions and for a full scheme can give us a wealth of information, and that information is catalogued in the ensuing two sub-sections.
The first insight that we would like to extract from such a stability analysis is the maximal CFL number for which the scheme is stable. For Eq. (3) the Fourier modes are indeed propagating with the velocity vector; therefore, the velocity vector sets the signal speed. For each choice of the spatial accuracy, we can choose a temporal accuracy for our SSP-RK scheme that is comparable or greater than the spatial accuracy. The upshot is that for each choice of the spatial and temporal accuracy, we can identify a maximal CFL number. Please realize that this involves sweeping through all velocities in two dimensions and for each choice of the velocity we have to sweep over all the wavenumbers that are permitted on the mesh. Stable CFL numbers are identified as the ones for which all possible wave vectors return an amplification matrix all of whose eigenvalues have an absolute value that is less than or equal to unity. Such a study of the maximal CFL number is documented in Sect. 5.1.
DG-like schemes can be very accurate, even when they are compared to their WENOlike counterparts. But we need to visually appreciate that. For that reason, we choose velocity vectors that make angles of 0°, 15°, 30°, and 45° to the mesh. For each of those velocity vectors, we sweep through all possible angles that the wave vector can make with respect to the velocity. This allows us to visualize the dissipation and dispersion errors of the schemes that we analyze. This information is shown in Sect. 5.2.

Maximal CFL Numbers from Stability Analysis
We identify the CFL number in each of the two directions by C x = v x Δt Δx and C y = v y Δt Δy . For each CFL number in either of the two directions, we sweep over all possible wave numbers k x Δx, k y Δy ∈ −π∕2, π∕2 × −π∕2, π∕2 . Figure 4 shows a colorized plot of the eigenvector of the amplification matrix with the largest absolute value for second, third, and fourth-order curl-free DG-like schemes. The white polygons in Fig. 4 identify the domain of stability for which the absolute value described above is less than or equal to unity. The white circles in Fig. 4 are the largest circles that can be inscribed in the polygons. The radii of those circles give us the largest effective CFL number that we should use for each of those DG-like schemes. Figure 4 is intended to give us a glimpse of the process to find the largest effective CFL number that we can find from our von Neumann stability analysis. Figure 4 corresponds to a situation where the order of temporal accuracy of our SSP-RK time stepping schemes indeed matched the spatial accuracy of the DG-like discretization. But we can use several possible SSP-RK schemes with each DG-like discretization, as long as the temporal accuracy is at least as large as the spatial accuracy. Table 1 shows the largest effective CFL number for a range of curl-free DG-like spatial discretizations and a range of temporal accuracies. In all cases, SSP-RK schemes were used for the temporal update. We see that for each scheme, an increasing temporal accuracy results in a larger effective CFL number, a result that conforms to the findings of Zhang and Shu [48] and Liu et al. [36].
Because we have erected the machinery of the von Neumann stability analysis, we can also use it to analyze the largest effective CFL number when other spatial discretizations are used. For example, when we retain only the time-evolution of the zeroth mode in Fig. 4 The domain of stability for a a second order in space and time curl-free DG-like scheme that uses SSK-RK2 timestepping, b a third order in space and time curl-free DG-like scheme that uses SSK-RK3 timestepping, and c a fourth order in space and time curl-free DG-like scheme that uses SSP-RK (5,4) timestepping. The CFL numbers in the x-and y-directions are denoted by C x and C y , and the color coding shows the absolute value of the largest eigenvalue of the amplification matrix. The white polygon identifies the full domain of stability and the white circle identifies the largest circle that can be fit within the domain of stability. The radius of the white circle, therefore, gives us the maximal CFL number Table 1 The largest effective CFL number for a range of curl-free DG-like spatial discretizations and a range of temporal accuracies Eqs. (10a) and (11a), we get a family of curl-free WENO-like schemes. For those schemes, all the higher modes, up to the desired order of accuracy, have to be reconstructed. Table 2 shows the largest effective CFL number for a range of curl-free WENO-like spatial discretizations and a range of temporal accuracies. In all cases, SSP-RK schemes were used for the temporal update. As expected, we see that the curl-free WENO-like schemes have CFL numbers much larger than the curl-free DG-like schemes. Table 3 shows the largest effective CFL number for a range of curl-free P1PM-like spatial discretizations and a range of temporal accuracies. In addition to retaining the time evolution of the modes in Eqs. (10a) and (11a), such schemes also evolve the first moments from Eqs. (10b) and (11b). As before, SSP-RK schemes were used for the temporal update. Comparing the CFL numbers from Table 3 to those from Tables 1 and 2, we see that curlfree P1PM-like schemes give us CFL numbers that are somewhat smaller than those of their WENO counterparts but substantially larger than their DG counterparts. We will further see in the next sub-section that P1PM-like schemes retain their first moments and that gives them dissipation and dispersion properties that are closer to their DG counterparts. The physical reason for that is because the linear mode retains most of the variation in the zone; consequently, much of the accuracy is retained. We, therefore, understand why curlfree P1PM-like schemes retain an important utilitarian position in the full range of schemes studied here.
This completes our study of the CFL number of curl-free DG-like schemes and their cousins.

Dissipation and Dispersion Properties of DG and PNPM Schemes
We now wish to study the dissipation and dispersion properties of the curl-free WENOlike, PNPM-like, and DG-like schemes. We will study these properties for second, third, and fourth order, so that we have a clear understanding of the improved wave propagation properties of these schemes with the increasing order. By the same token, we will also be able to inter-compare between the WENO-like, PNPM-like, and DG-like schemes. We expect that retaining more moments and evolving them consistently with the governing PDE, should give us schemes with improved wave propagation characteristics. For all the data shown in this sub-section, the temporal accuracy was made to match the spatial accuracy. While the von Neumann stability analysis for curl-free WENO-like schemes was already documented in Balsara et al. [12], we present it again here so that one can intercompare with the PNPM-like and DG-like schemes. The von Neumann stability analysis of the curl-free PNPM-like and DG-like schemes is being presented for the very first time here.
In each instance, we choose velocity vectors that make angles of 0°, 15°, 30°, and 45° to the mesh and use a CFL number that is 0.9 times the maximum shown in either Table 1 or Table 2 or Table 3. We then let the wave number k x , k y sweep through all angles, from Table 3 The largest effective CFL number for a range of curl-free P1PM-like spatial discretizations and a range of temporal accuracies P1P2 P1P3 SSP-RK3 0.390 3 -SSP-RK(5,4) 0.626 0 0.679 9 In all cases, SSP-RK schemes were used for the temporal update −180 • to +180 • relative to the velocity vector v x , v y . For each of those angles, we plot out 1 − |amplification factor| for the scheme. If this number is non-negative and close to zero, it indicates that the scheme has low dissipation. The phase of the amplification factor gives us a measure of the propagation speed of the waves. For each angle between the velocity vector and the wave number, we also plot out the error in the phase speed. As the wavelength increases relative to the zone size, we expect the schemes to propagate waves with increasing accuracy. As a result, we consider wavelengths that are 5 Δx , 10Δx , and 15 Δx . In the next eight figures that follow, wavelengths of 5 Δx are always shown with a blue color, wavelengths of 10Δx are always shown with a green color and wavelengths of 15 Δx are always shown with a red color. Figures 5, 6, and 7 show the wave propagation characteristics of curl-free WENO-like schemes at second, third and fourth order, respectively. We see that as we go from second to fourth order, the dissipation (as measured by 1 − |amplification factor| ) improves by an order of magnitude for each of the three wavelengths considered here. Similarly, as we go from second to fourth order, the phase error is also reduced by an order of magnitude. Table 4 shows the minimum of the absolute value of the amplification factor for all possible velocity directions and angles between the velocity and wave number vectors for curlfree WENO-like schemes when we have waves with the wavelength 5Δx, 10Δx, and 15Δx. In the same table, we also show the maximum phase error for the similar situation and for the same wavelengths. In other words, Table 4 was extracted from Figs. 5, 6, and 7 and allows us to quantify the most significant aspect of those figures. Table 4, therefore, allows us to make an important practical decision. Say we want to carry out a simulation with WENO-like schemes we want to meet a target set of dissipation and dispersion properties, Table 4 shows us what our options are. We may indeed choose a lower order scheme and use a lot of zones to cover the characteristic wavelength in the simulation. But we see that we can also choose a higher-order scheme and use fewer zones to cover the characteristic wavelength in the simulation.
There is no second-order P1P1 scheme, because such a scheme would be identical to a second-order DG scheme. However, our study of CFL numbers has shown us that curl-free third-order P1P2-like and fourth-order P1P3-like schemes still retain a very robust CFL number. We, therefore, want to know whether such schemes have superior wave propagation characteristics relative to the WENO-like schemes that we studied in the previous paragraph. Figures 8 and 9 show the wave propagation characteristics of curl-free P1P2-like Table 4 The minimum of the absolute value of the amplification factor for all possible velocity directions and all angles between the velocity and wave number vectors for curl-free WENO-like schemes when we have waves with wavelength 5Δx, 10Δx, and 15Δx Min of |amplification factor|   The wave propagation characteristics for curl-preserving second order WENO-like schemes. a-d one minus the absolute value of the amplification factor when the velocity vector makes angles of 0°, 15°, 30°, and 45° relative to the x-direction of the 2D mesh. e-h the phase error, again for the same angles. The 2D wave vector can make any angle relative to the 2D direction of velocity propagation, therefore, the amplitude and phase information are shown w.r.t. the angle made between the velocity direction and the direction of the wave vector. In each plot, the blue curve refers to waves that span 5 cells per wavelength; the green curve refers to waves that span 10 cells per wavelength; the red curve refers to waves that span 15 waves per wavelength Fig. 6 The wave propagation characteristics for curl-preserving third order WENO-like schemes. a-d one minus the absolute value of the amplification factor when the velocity vector makes angles of 0°, 15°, 30°, and 45° relative to the x-direction of the 2D mesh. e-h the phase error, again for the same angles. The 2D wave vector can make any angle relative to the 2D direction of velocity propagation, therefore, the amplitude and phase information are shown w.r.t. the angle made between the velocity direction and the direction of the wave vector. In each plot, the blue curve refers to waves that span 5 cells per wavelength; the green curve refers to waves that span 10 cells per wavelength; the red curve refers to waves that span 15 waves per wavelength 1 3 Fig. 7 The wave propagation characteristics for curl-preserving fourth order WENO-like schemes. a-d one minus the absolute value of the amplification factor when the velocity vector makes angles of 0°, 15°, 30° and 45°, relative to the x-direction of the 2D mesh. e-h the phase error, again for the same angles. The 2D wave vector can make any angle relative to the 2D direction of velocity propagation, therefore, the amplitude and phase information are shown w.r.t. the angle made between the velocity direction and the direction of the wave vector. In each plot, the blue curve refers to waves that span 5 cells per wavelength; the green curve refers to waves that span 10 cells per wavelength; the red curve refers to waves that span 15 waves per wavelength 1 3 and P1P3-like schemes at third and fourth order, respectively. We should, therefore, compare Fig. 8 to Fig. 6 because they both pertain to third-order schemes. Similarly, we should compare Fig. 9 to Fig. 7 because they both pertain to fourth-order schemes. The results are quite interesting. We see that Fig. 8 and Fig. 6 show comparable quality of wave propagation indicating that at third order the advantages are minimal. This lack of significant improvement might have to do with the fact that SSP-RK3 time stepping has excessive stabilization. Now let us turn to comparing Figs. 9 and 7. At fourth order, we do see that the P1P3-like scheme outperforms the WENO-O4 scheme by almost an order of magnitude. It shows the value of designing PNPM schemes as half-way houses between WENO and DG schemes. Table 5 shows the minimum of the absolute value of the amplification factor for all possible velocity directions and all angles between the velocity and wave number vectors for curl-free PNPM-like schemes when we have waves with wavelength 5Δx, 10Δx, and 15Δx. In the same table, we also show the maximum phase error for a similar situation and for the same wavelengths. In other words, Table 5 was extracted from Figs. 8 and 9 and allows us to quantify the most significant aspect of those figures. Again, Table 5 can help with practical decision-making. It shows us, for example, that fourth-order P1P3-like schemes do give us a substantial improvement over fourth-order WENO-like scheme while incurring only a modest increase in computational complexity.
While they have the smallest CFL numbers, DG-like schemes hold out the promise of almost spectral-like accuracy with increasing order of accuracy; for finite volume-based approaches this is now viewed as an accepted fact. We are now in a position to test that contention as it pertains to curl-free DG-like schemes. Figures 10, 11, and 12 show the wave propagation characteristics of curl-free DG-like schemes at second, third, and fourth order, respectively. With increasing order, the curl-free DG-like schemes do show significant improvement when inter-compared amongst themselves. Let us, therefore, compare across algorithms since we have all the data concatenated in one place. Figure 10 should be compared to Fig. 5. Figure 11 should be compared to Figs. 6 and 8. Figure 12 should be compared to Figs. 7 and 9. We see that the wave propagation characteristics of the secondorder DG-like scheme are entirely competitive with the wave propagation characteristics of the fourth-order WENO-like scheme. The fourth-order DG-like scheme is also somewhat superior to the fourth-order P1P3-like scheme, but please recall that this comes with a substantial increase in programming complexity and a decrease in the CFL number. Table 6 shows the minimum of the absolute value of the amplification factor for all possible velocity directions and all angles between the velocity and wave number vectors for curl-free Table 5 The minimum of the absolute value of the amplification factor for all possible velocity directions and all angles between the velocity and wave number vectors for curl-free PNPM-like schemes for waves with wavelength 5Δx, 10Δx, and 15Δx We also show the maximum phase error for the same wavelengths Min of |amplification factor|  Fig. 8 The wave propagation characteristics for curl-preserving third order P1P2-like schemes. a-d one minus the absolute value of the amplification factor when the velocity vector makes angles of 0°, 15°, 30°, and 45° relative to the x-direction of the 2D mesh. e-h the phase error, again for the same angles. The 2D wave vector can make any angle relative to the 2D direction of velocity propagation, therefore, the amplitude and phase information are shown w.r.t. the angle made between the velocity direction and the direction of the wave vector. In each plot, the blue curve refers to waves that span 5 cells per wavelength; the green curve refers to waves that span 10 cells per wavelength; the red curve refers to waves that span 15 waves per wavelength Fig. 9 The wave propagation characteristics for curl-preserving fourth order P1P3-like schemes. a-d one minus the absolute value of the amplification factor when the velocity vector makes angles of 0°, 15°, 30°, and 45° relative to the x-direction of the 2D mesh. e-h the phase error, again for the same angles. The 2D wave vector can make any angle relative to the 2D direction of velocity propagation, therefore, the amplitude and phase information are shown w.r.t. the angle made between the velocity direction and the direction of the wave vector. In each plot, the blue curve refers to waves that span 5 cells per wavelength; the green curve refers to waves that span 10 cells per wavelength; the red curve refers to waves that span 15 waves per wavelength DG-like schemes when we have waves with wavelength 5Δx, 10Δx, and 15Δx. In the same table, we also show the maximum phase error for a similar situation and for the same wavelengths. In other words, Table 6 was extracted from Figs. 10, 11, and 12 and allows us to quantify the most significant aspect of those figures. We see that the dissipation and dispersion characteristics of the curl-free DG-like schemes that we have designed are indeed excellent. Comparing Tables 5 and 6 we also see that the curl-free PNPM-like schemes are not far behind. Therefore, DG-like schemes are the go-to scheme when superlative performance is the only driving consideration. However, if one wants lower computational complexity and more robust timesteps, the PNPM-like schemes also present themselves as attractive choices. While this section has been focused on curl-free methods, we point out that curl-preserving methods only require a few additional terms in Eqs. (10) and (11) compared to curl-free schemes. Without an underlying fluid-dynamical PDE system that supplies additional terms like density, velocity, and temperature, it is not possible to obtain the source terms and other types of terms which would make Eq. (2) curl-preserving instead of curlfree. Furthermore, computer algebra systems are just not adept enough to support a more extensive stability analysis for larger PDE systems. For these reasons, the analysis presented here is focused on curl-free methods, but the insights developed here will extend to all curl-preserving methods.

Numerical Results
In this section, we present numerical experiments confirming that the developed curl-free DG schemes reach their expected design accuracies. Results for all the DG-like (PNPN) and PNPM (P0PN for WENO-like and P1PN for HWENO-like) up to fourth order are reported for two smooth test problems. All the tests are run with 95% of the maximal CFL number (Tables 1, 2, 3) of the respective scheme. Table 6 The minimum of the absolute value of the amplification factor for all possible velocity directions and all angles between the velocity and wave number vectors for curl-free DG-like schemes when we have waves with wavelength 5Δx, 10Δx, and 15Δx We also show the maximum phase error for the same wavelengths Min of |amplification factor| λ = 5Δx λ = 10Δx λ = 15Δx  . 10 The wave propagation characteristics for curl-preserving second order DG-like schemes. a-d one minus the absolute value of the amplification factor when the velocity vector makes angles of 0°, 15°, 30°, and 45°, relative to the x-direction of the 2D mesh. e-h the phase error, again for the same angles. The 2D wave vector can make any angle relative to the 2D direction of velocity propagation, therefore, the amplitude and phase information are shown w.r.t. the angle made between the velocity direction and the direction of the wave vector. In each plot, the blue curve refers to waves that span 5 cells per wavelength; the green curve refers to waves that span 10 cells per wavelength; the red curve refers to waves that span 15 waves per wavelength 1 3

Plane-Wave Test Problem
The first test problem is the propagation of a plane wave in a Cartesian domain of size −1∕2, 1∕2 2 with periodic boundaries. The plane wave is advected with veloc- ity v x = v y = 1 diagonally through the domain. The curl-free field is initialized from a potential where we set k x = k y =2π. The x-and y-field components of the curl-free field are then given by The setup is run up to time t f = 1 , by which time the plane wave has propagated once through the domain, and the accuracy of the schemes is evaluated. Moreover, we also show the preservation of the quadratic field energy (J 2 x + J 2 y )∕2 highlighting the dissipation characteristics of the schemes. Table 7 shows the L 1 and L ∞ errors at the final time for the PNPN (N = 1, 2, 3) DG-like schemes for resolutions from 8 × 8 up to 64 × 64. Table 8 shows the convergence analysis of the P0PN (N = 1, 2, 3) WENO-like schemes. Table 9 shows the convergence analysis of the P1PN (N = 1, 2, 3) HWENO-like schemes. The tables also catalogue the final quadratic energy as a fraction of the initial quadratic energy. We observe that all schemes reach their design accuracy. We also take note of the improved quadratic field energy preservation with increasing order of accuracy. Note that we did nothing special in the scheme to ensure that quadratic field energy is conserved; as a result, the rather good preservation of quadratic energy is entirely a consequence of the accuracy of the method. This is especially true for the DG-like schemes which preserve quadratic energy very well especially as the resolution is increased. The WENO-like schemes show slightly inferior energy preservation characteristics. However, the latter allow much larger time steps due to their larger allowed CFL numbers. The HWENO-like schemes show nearly the same quadratic energy preservation properties as the DG-like schemes and, furthermore, allow larger time steps similar to the WENO-like schemes. Consequently, we see that the HWENO-like schemes may be viewed as an efficient compromise between the extreme accuracy of the DG-like schemes and the much larger time steps of the WENO-like schemes.

Vortex Test Problem
The second test problem consists of the advection of a localized curl-free vortex simi-  The wave propagation characteristics for curl-preserving third order DG-like schemes. a-d one minus the absolute value of the amplification factor when the velocity vector makes angles of 0°, 15°, 30°, and 45° relative to the x-direction of the 2D mesh. e-h the phase error, again for the same angles. The 2D wave vector can make any angle relative to the 2D direction of velocity propagation, therefore, the amplitude and phase information are shown w.r.t. the angle made between the velocity direction and the direction of the wave vector. In each plot, the blue curve refers to waves that span 5 cells per wavelength; the green curve refers to waves that span 10 cells per wavelength; the red curve refers to waves that span 15 waves per wavelength Fig. 12 The wave propagation characteristics for curl-preserving fourth order DG-like schemes. a-d one minus the absolute value of the amplification factor when the velocity vector makes angles of 0°, 15°, 30°, and 45° relative to the x-direction of the 2D mesh. e-h the phase error, again for the same angles. The 2D wave vector can make any angle relative to the 2D direction of velocity propagation, therefore, the amplitude and phase information are shown w.r.t. the angle made between the velocity direction and the direction of the wave vector. In each plot, the blue curve refers to waves that span 5 cells per wavelength; the green curve refers to waves that span 10 cells per wavelength; the red curve refers to waves that span 15 waves per wavelength

3
The advection velocity is set to v x = v y = 1 . The problem is simulated for the time t f = 20 , by which point the vortex was advected once through the square domain in the diagonal direction till it returns to its initial position. At final time point, we measure the accuracy in the L 1 and L ∞ errors norms. Moreover, we also show the preservation of the quadratic field energy (J 2 x + J 2 y )∕2 highlighting the dissipation characteristics of the schemes. Table 10 shows the L 1 and L ∞ errors at the final time for the PNPN (N = 1, 2, 3) RKDGlike schemes for resolutions from 16 × 16 up to 256 × 256. Table 11 shows the convergence analysis of the P0PN (N = 1, 2, 3) WENO-like schemes. Table 12 shows the convergence analysis of the P1PN (N = 1, 2, 3) HWENO-like schemes. The tables also catalogue the final quadratic energy as a fraction of the initial quadratic energy. We observe that all schemes reach their design accuracy on the chosen mesh resolutions, even for this highly spatially localized vortex. Note that much of the field variation is confined around a circle with unit radius, corresponding to one-tenth of the computational domain. We find that the presented schemes concurrently have quadratic energy preservation with, nevertheless, excellent accuracy.
As before, we did nothing special in the scheme to ensure that the quadratic field energy is conserved. Consequently, the rather good preservation of quadratic energy is entirely a consequence of the accuracy of the method. For the quadratic field energy, we observe similar trends to the plane wave test problem. To further highlight the point, Fig. 13 shows the quadratic field energy preservation characteristics from Tables 10, 11, and 12, graphically. Each panel shows all the available schemes up to fourth order of accuracy. Figure 13a underlines that the curl-free DG-like schemes show excellent energy-preserving properties with increasing order of accuracy. In Fig. 13b, we observe again that the energy preserving properties of the curl-free WENO-like schemes are somewhat inferior in the preasymptotic regime. However, the improving trend with increasing order of accuracy is also clearly visible. In Fig. 13c, we see that the curl-free P1PN-like schemes share almost the same preservation properties as the curl-free DG-like schemes. The latter fact, and their substantially larger allowed CFL numbers (hence, time steps), highlight again that the P1PN-like schemes are very efficient curl constraint-preserving methods that share desirable qualities from both full DG-like and WENO-like schemes.
Lastly, we also catalogue that the schemes designed here are curl-preserving and can indeed reach the curl-free limit even when they are integrated for long simulation times. Dumbser et al. [28] have shown that if a classical higher-order Godunov scheme is applied to Eq. (2), the numerical instability manifests itself as an explosive increase in the discrete circulation when the simulation is run over long periods of time; see [28,Fig. 5]. Therefore, we would like to demonstrate that the discrete circulation is held down to machine precision when the simulation is run for a long period of time when we use the methods designed here. We would indeed like to go one step further and plot out the time series of the maximum pointwise error in the curl of J. In other words, we know that the vector field is a polynomial, see Eq. (14) or Eq. (17) for instance, so that we can evaluate the pointwise curl at any point within a zone (because the polynomials are differentiable). We then choose 2 × 2, 3 × 3, or 4 × 4 uniformly spaced points that are internal to each zone at second, third, and fourth orders respectively. We then evaluate the maximum of the absolute value of the curl at each and every internal point for all the zones on the mesh and we plot this maximum value as a function of time. Let us ask why this demonstration matters? We see from Eqs. (10) and (11) that in a curl-preserving scheme we will also need terms like ⟨ , and other terms like it, at the edges of the mesh. These have to be evaluated from either side of the edge that is being considered. Therefore, when we approach the curl-free limit, a curl-preserving scheme should naturally obtain curl-free evolution. This demonstration that the maximum pointwise error in the curl of J remains close to machine zero over long simulation times guarantees that such a limit is met.
To show that the curl remains close to machine zero at all points on the mesh even during long-time integration, we have run the vortex problem on a 64 × 64 zone mesh to a final time of 200. This time corresponds to the vortex making ten passages through the periodic computational domain. Figure 14 shows the maximum pointwise error of the curl of J as a function of time for a 64 × 64 zone run of the vortex problem. Figure 14a shows the evolution of the maximum pointwise curl as a function of time for the second, third, and fourth order curl-free DG-like schemes. Figure 14b shows the evolution of the maximum pointwise curl as a function of time for the third and fourth order curl-free P1PN-like schemes. Figure 14c shows the evolution of the maximum pointwise curl as a function of time for the second, third, and fourth order curl-free WENO-like schemes. The figure shows that all our curl-preserving schemes can preserve the curl constraint up to machine accuracy in simulations that are run for long integration times.

Conclusions
Novel classes of PDEs have recently emerged and the physics of the PDEs requires keeping strict control of the curl of one or more vector fields. The PDEs are hyperbolic systems of great interest to science and engineering. Many of the hyperbolic systems resulting from the GPR formulation for hyperelasticity and compressible multiphase flow with and without surface tension have curl-preserving update equations (Godunov and Romenski [32], Romenski [40], Romenski et al. [41], Peshkov and Romenski [37,38], Dumbser et al. [27,30,31], Schmidmayer et al. [42]). The equations of general relativity, when cast in the FO-CCZ4 formulation, also have such a structure (Alic et al. [1,2], Brown et al. [18], Dumbser et al. [28,29]). Similarly, it has recently become possible to recast Schrödinger's equation in first-order hyperbolic form, and the time-evolution of this important equation also has curl-preserving constraints (Dhaouadi et al. [25], Busto et al. [19]). Experience has shown that if nothing special is done to account for the curl-preserving vector field, it can blow up in a finite amount of simulation time (Dumbser et al. [28]).
Prior work has shown that classical zone-centered Godunov methods can be adapted to such systems only if a GLM-type cleaning approach is included to suppress the build up of circulation on the mesh (Dumbser et al. [28]). The two-fold problem with this approach is as follows: (i) we often get a very large system of Lagrange multipliers that are not part of the original PDE system and (ii) the signal speed with which the Lagrange multipliers have to be advected often exceeds the physical signal speed in the problem by a substantial margin. Another alternative is to solve an elliptical system at every timestep (Boscheri et al. [17]), which makes each timestep very expensive. In Balsara et al. [12] we first presented curl-preserving WENO-like methods that overcame both of the above-mentioned limitations. The methods were based on inventing a novel globally curl-preserving reconstruction strategy that reconstructs the vector field over the zone's volume using the components of the vector field that were collocated at the edges of the mesh. Non-linear hybridization, via WENO methods, was seamlessly built into the curl-preserving reconstruction strategy. These edge-centered components were updated using multidimensionally upwinded potentials that were collocated at the vertices of the mesh. Multidimensional Riemann solvers, designed by the first author, provided the requisite multidimensional upwinding.
The resulting highly stable finite volume-like schemes for curl-preserving systems had the following desirable properties. (i) They did not blow up even after very long integration times. (ii) They did not need GLM-style cleaning with very high signal speeds. (iii) They could operate with large explicit timesteps. (iv) They did not require the solution of an  elliptic system. And (v) they could be extended to higher orders while incorporating nonlinear hybridization using WENO-like methods. It is, therefore, desirable to invent DG-like and PNPM-like variants of these WENO-like schemes so that they can inherit the same desirable features-such a task is fulfilled in this paper.
Since we know that DG and PNPM schemes provide more accurate alternatives to WENO schemes, it becomes interesting to design curl-free and curl-preserving variants of the such schemes. In this paper, we present for the very first time, globally curl-preserving  DG-like and PNPM-like schemes that share the beneficial traits of the globally curl-preserving WENO-like schemes designed by Balsara et al. [12]. The higher moments of the vector components that live in the edges of the mesh are, therefore, endowed with timeevolution that is consistent with the governing equations. This is accomplished by making a Galerkin projection within each edge that results in a weak form of update equation  for the higher-order edge-centered moments. The update utilizes the multidimensionally upwinded potentials at the vertices of the mesh. Such update equations have been documented in Sect. 2 for the model Eq.
(2) and some nuances of the curl-preserving reconstruction, and how it relates to traditional DG schemes, are highlighted in Sect. 3. It is well-known that zone-centered DG schemes have wave propagation characteristics superior to zone-centered WENO schemes. Such a superior behavior can be revealed by conducting a von-Neumann stability analysis of either scheme and inter-comparing the results. It is, therefore, interesting to carry out a von Neumann stability analysis of our newly-developed globally curl-free and curl-preserving DG-like and PNPM-like schemes. In Sect. 4 we present details of our von Neumann stability analysis. To highlight the curl-free aspect of the evolution, the analysis must absolutely be done in two or more dimensions. Therefore, our analysis is two-dimensional by its very design. By pushing the capabilities of computer algebra systems to the limits, we have been able to extend this two-dimensional curl-preserving von Neumann stability analysis up to fourth order of accuracy.
Section 5 shows the results of this von Neumann stability analysis. We present such an analysis for globally curl-free WENO-like, PNPM-like and DG-like schemes to facilitate inter-comparison. In Sect. 5.1 the limiting CFL numbers for all these schemes are derived and documented in Tables 1, 2, and 3. We find, unsurprisingly, that WENO-like schemes offer the largest CFL numbers while DG-like schemes restrict us to substantially smaller CFL numbers. The PNPM-like schemes give us quite large CFL numbers at a muchreduced computational complexity. In Sect. 5.2 we document the dissipation and dispersion properties of the same three schemes. We do this for waves with wavelengths 5, 10, and 15 times the zone size. For each family of schemes, the dissipation and dispersion properties do indeed improve with increasing order of accuracy, as expected. This is shown in the figures associated with Sect. 5.2 and also in Tables 4, 5, and 6. We find that WENOlike schemes have dissipation and dispersion properties that are noticeably inferior to the DG-like schemes at the same order. However, we find that PNPM-like schemes have dissipation and dispersion properties that approach those of DG-like schemes at the same order while offering substantially larger CFL numbers. Section 6 presents numerical results, where we show that our methods meet their design accuracies. We also show that with increasing order of accuracy the methods become very good at preserving quadratic energy. This is a welcome result, because the methods were not intentionally designed to preserve quadratic energy; yet they seem to do a good job. When the evolution of the PDE is curl-free, our methods also hold down the discrete circulation to machine accuracy over long integration times. The importance of this fact in the design of curl-preserving schemes is also discussed.