Stabilization in relation to wavenumber in HDG methods

Simulation of wave propagation through complex media relies on proper understanding of the properties of numerical methods when the wavenumber is real and complex. Numerical methods of the Hybrid Discontinuous Galerkin (HDG) type are considered for simulating waves that satisfy the Helmholtz and Maxwell equations. It is shown that these methods, when wrongly used, give rise to singular systems for complex wavenumbers. A sufficient condition on the HDG stabilization parameter for guaranteeing unique solvability of the numerical HDG system, both for Helmholtz and Maxwell systems, is obtained for complex wavenumbers. For real wavenumbers, results from a dispersion analysis are presented. An asymptotic expansion of the dispersion relation, as the number of mesh elements per wave increase, reveal that some choices of the stabilization parameter are better than others. To summarize the findings, there are values of the HDG stabilization parameter that will cause the HDG method to fail for complex wavenumbers. However, this failure is remedied if the real part of the stabilization parameter has the opposite sign of the imaginary part of the wavenumber. When the wavenumber is real, values of the stabilization parameter that asymptotically minimize the HDG wavenumber errors are found on the imaginary axis. Finally, a dispersion analysis of the mixed hybrid Raviart-Thomas method showed that its wavenumber errors are an order smaller than those of the HDG method.


Background
Wave propagation through complex structures, composed of both propagating and absorbing media, are routinely simulated using numerical methods. Among the various numerical methods used, the Hybrid Discontinuous Galerkin (HDG) method has emerged as an attractive choice for such simulations. The easy passage to high order using interface unknowns, condensation of all interior variables, availability of error estimators and adaptive algorithms, are some of the reasons for the adoption of HDG methods.
It is important to design numerical methods that remain stable as the wavenumber varies in the complex plane. For example, in applications like computational lithography, one finds absorbing materials with complex refractive index in parts of the domain of simulation. Other examples are furnished by meta-materials. A separate and important reason for requiring such stability emerges in the computation of resonances by iterative searches in the complex plane. It is common for such iterative algorithms to solve a source problem with a complex wavenumber as its current iterate. Within such algorithms, if the HDG method is used for discretizing the source problem, it is imperative that the method remains stable for all complex wavenumbers.
One focus of this study is on complex wavenumber cases in acoustics and electromagnetics, motivated by the above-mentioned examples. Ever since the invention of the HDG method in [4], it has been further developed and extended to other problems in many works (so many so that it is now impractical to list all references on the subject here). Of particular interest to us are works that applied HDG ideas to wave propagation problems such as [6,8,10,11,12,13,14]. We will make detailed comparisons with some of these works in a later section. However, none of these references address the stability issues for complex wavenumber cases. While the choice of the HDG stabilization parameter in the real wave number case can be safely modeled after the well-known choices for elliptic problems [5], the complex wave number case is essentially different. This will be clear right away from a few elementary calculations in the next section, which show that the standard prescriptions of stabilization parameters are not always appropriate for the complex wave number case. This then raises further questions on how the HDG stabilization parameter should be chosen in relation to the wavenumber, which are addressed in later sections.
Another focus of this study is on the difference in speeds of the computed and the exact wave, in the case of real wavenumbers. By means of a dispersion analysis, one can compute the discrete wavenumber of a wave-like solution computed by the HDG method, for any given exact wavenumber. An extensive bibliography on dispersion analyses for the standard finite element method can be obtained from [1,7]. For nonstandard finite element methods however, dispersion analysis is not so common [9], and for the HDG method, it does not yet exist. We will show that useful insights into the HDG method can be obtained by a dispersion analysis. In multiple dimensions, the discrete wavenumber depends on the propagation angle. Analytic computation of the dispersion relation is feasible in the lowest order case. We are thus able to study the influence of the stabilization parameter on the discrete wavenumber and offer recommendations on choosing good stabilization parameters. The optimal stabilization parameter values are found not to depend on the wavenumber. In the higher order case, since analytic calculations pose difficulties, we conduct a dispersion analysis numerically.
We begin, in the next section, by describing the HDG methods. We set the stage for this study by showing that the commonly chosen HDG stabilization parameter values for elliptic problems are not appropriate for all complex wavenumbers. In the subsequent section, we discover a constraint on the stabilization parameter, dependent on the wavenumber, that guarantees unique solvability of both the global and the local HDG problems. Afterward, we perform a dispersion analysis for both the HDG method and a mixed method and discuss the results.

Methods of the HDG type
We borrow the basic methodology for constructing HDG methods from [4] and apply it to the time-harmonic Helmholtz and Maxwell equations (written as first order systems). While doing so, we set up the notations used throughout, compare the formulation we use with other existing works, and show that for complex wavenumbers there are stabilization parameters that will cause the HDG method to fail.
Undesirable HDG stabilization parameters for the Helmholtz system. We begin by considering the lowest order HDG system for Helmholtz equation. Let k be a complex number. Consider the Helmholtz system on Ω ⊂ R 2 with homogeneous Dirichlet boundary conditions,î where f ∈ L 2 (Ω). Let T h denote a square or triangular mesh of disjoint elements K, so Ω = ∪ K∈T h K, and let F h denote the collection of edges. The HDG method produces an approximation (⃗ u, φ,φ) to the exact solution ( ⃗ U, Φ,Φ), whereΦ denotes the trace of Φ on the collection of element boundaries ∂T h . The HDG solution with polynomial spaces V (K), W (K), and M (F ) specified differently depending on element type: Triangles Squares Here, for a given domain D, P p (D) denotes polynomials of degree at most p, and Q p (D) denotes polynomials of degree at most p in each variable. The HDG solution satisfies
One of the main reasons to use an HDG method is that all interior unknowns (⃗ u, φ) can be eliminated to get a global system for solely the interface unknowns (φ). This is possible whenever the local system is uniquely solvable. (For details on this elimination and other perspectives on HDG methods, see [4].) In the lowest order (p = 0) case, on a square element K of side length h, if we use a basis in the following order then the element matrix for the system (4) is This shows that if then M is singular, and so the HDG method will fail. The usual recipe of choosing τ = 1 is therefore inappropriate when k is complex valued.
Intermediate case of the 2D Maxwell system. It is an interesting exercise to consider the 2D Maxwell system before going to the full 3D case. In fact, the HDG method for the 2D Maxwell system can be determined from the HDG method for the 2D Helmholtz system. The 2D Maxwell system iŝ where J ∈ L 2 (Ω), and the scalar curl ∇× ⋅ and the vector curl ⃗ ∇× ⋅ are defined by Here This also shows that the HDG method for Helmholtz equation should yield an HDG method for the 2D Maxwell system. We thus conclude that there exist stabilization parameters that will cause the HDG system for 2D Maxwell system to fail.
To examine this 2D HDG method, if we let ⃗ H and E denote the HDG approximations for R⃗ r and E, respectively, then the HDG system (2) with ⃗ u and φ replaced by −R ⃗ H and E, respectively, gives where ⃗ t = R⃗ n the tangent vector, and we have used the 2D cross product defined by In particular, the numerical flux prescription (3) implies where R ⃗ H denotes the numerical trace of R ⃗ H. We rewrite this in terms of ⃗ H and E, to obtainĤ One may rewrite this again, asĤ This expression is notable because it will help us consistently transition the numerical flux prescription from the Helmholtz to the full 3D Maxwell case discussed next. A comparison of this formula with those in the existing literature is included in Table 1.
The 3D Maxwell system. Consider the 3D Maxwell system on Ω ⊂ R 3 with a perfect electrically conducting boundary condition, where ⃗ J ∈ (L 2 (Ω)) 3 . For this problem, T h denotes a cubic or tetrahedral mesh, and F h denotes the collection of mesh faces. The HDG method approximates the exact solution The discrete spaces are defined by with polynomial spaces Y (K) and N (F ) specified by: Tetrahedra Cubes Our HDG method for (8) is where, in analogy with (7), we now set numerical flux bŷ where ( ⃗ E −Ê) t denotes the tangential component, or equivalentlŷ Note that the 2D system (6) is obtained from the 3D Maxwell system (8) by assuming symmetry in x 3 -direction. Hence, for consistency between 2D and 3D formulations, we should have the same form for the numerical flux prescriptions in 2D and 3D. The HDG method is then equivalently written as For comparison with other existing formulations, see Table 1.
Again, let us look at the solvability of the local element problem In the lowest order (p = 0) case, on a cube element K of side length h, if we use a basis in the following order then the 6 × 6 element matrix for the system (11) is where I 3 denotes the 3 × 3 identity matrix. Again, exactly as in the Helmholtz casecf.
then the local static condensation required in the HDG method will fail in the Maxwell case also.
Behavior on tetrahedral meshes. For the lowest order (p = 0) case on a tetrahedral element, just as for the cube element described above, there are bad stabilization parameter values. Consider, for example, the tetrahedral element of size h defined by with a basis ordered as in (12). The element matrix for the system (11) is then We immediately see that the rows become linearly dependent if Hence, for τ = −îkh (3 √ 3 + 6), the HDG method will fail on tetrahedral meshes. For orders p ⩾ 1, the element matrices are too complex to find bad parameter values so simply. Instead, we experiment numerically. Setting τ = −î, which is equivalent to the choice made in [14] (see Table 1), we compute the smallest singular value of the element matrix M (the matrix of the left hand side of (11) with K set by (14)) for a range of normalized wavenumbers kh. Figures 1a and 1b show that, for orders p = 1 and p = 2, there are values of kh for which τ = −î results in a singular value very close to zero. Taking a closer look at the first nonzero local minimum in Figure 1a, we find that the local matrix corresponding to normalized wavenumber kh ≈ 7.49 has an estimated condition number exceeding 3.9 × 10 15 , i.e., for all practical purposes, the element matrix is singular. To illustrate how a different choice of stabilization parameter τ can affect the conditioning of the element matrix, Figures 1c and 1d show the smallest singular values for the same range of kh, but with τ = 1. Clearly the latter choice of τ is better than the former.
From another perspective, Figure 1e shows the smallest singular value of the element matrix as τ is varied in the complex plane, while fixing kh to 1. Figure 1f is similar except that we fixed kh to the value discussed above, approximately 7.49. In both cases,

Results on unisolvent stabilization
We now turn to the question of how we can choose a value for the stabilization parameter τ that will guarantee that the local matrices are not singular. The answer, given by a condition on τ , surprisingly also guarantees that the global condensed HDG matrix is nonsingular. These results are based on a tenuous stability inherited from the fact nonzero polynomials are never waves, stated precisely in the ensuing lemma. Then we give the condition on τ that guarantees unisolvency, and before concluding the section, present some caveats on relying solely on this tenuous stability.
As is standard in all HDG methods, the unique solvability of the element problem allows the formulation of a condensed global problem that involves only the interface unknowns. We introduce the following notation to describe the condensed systems. First, for Maxwell's equations, for any . If all the sources in (10) vanish, then the condensed global problem forÊ ∈ N h takes the form where By following a standard procedure [4] we can express a(⋅, ⋅) explicitly as follows: Here we have used the complex conjugate of (15b) with ⃗ w = ⃗ H Λ , along with the definition ofĤ Λ , and then used (15a).
Similarly, for the Helmholtz equation, let (⃗ u η , φ η ) ∈ V h × W h denote the fields such that, for all K ∈ T h , the functions (⃗ u η K , φ η K ) solve the element problem (4) for given dataφ = η. If the sources in (2) vanish, then the condensed global problem forφ ∈ M h is written as where the form is found, as before, by the standard procedure: The sesquilinear forms a(⋅, ⋅) and b(⋅, ⋅) are used in the main result, which gives sufficient conditions for the solvability of the local problems (11), (4) and the global problems (16), (17). Before proceeding to the main result, we give a simple lemma, which roughly speaking, says that nontrivial harmonic waves are not polynomials.
Proof. We use a contradiction argument. If E ≡ ⃗ 0, then we may assume without loss of generality that at least one of the components of ⃗ E is a polynomial of degree exactly p. But this contradicts k 2 ⃗ E = ⃗ ∇×( ⃗ ∇× ⃗ E) because all components of ⃗ ∇×( ⃗ ∇× ⃗ E) are polynomials of degree at most p − 2. Hence ⃗ E ≡ ⃗ 0. An analogous argument can be used for the Helmholtz case as well.
Then, in the Maxwell case, the local element problem (11) and the condensed global problem (16) are both unisolvent. Under the same condition, in the Helmholtz case, the local element problem (4) and the condensed global problem (17) are also unisolvent.
Proof. We first prove the theorem for the local problem for Maxwell's equations. Assume (18) holds and setÊ = ⃗ 0 in the local problem (11). Unisolvency will follow by showing that ⃗ E and ⃗ H must equal ⃗ 0. Choosing ⃗ v = ⃗ E, and ⃗ w = ⃗ H, then subtracting (11b) from (11a), we get Under condition (18b), we immediately have that the fields ⃗ E and ⃗ H are zero on K. Otherwise, (18a) implies ⃗ E × ⃗ n ∂K = 0, and then (11) giveŝ By Lemma 1 this equation has no nontrivial solutions in the space Y (K). Thus, the element problem for Maxwell's equations is unisolvent. We prove that the global problem for Maxwell's equations is unisolvent by showing thatÊ = ⃗ 0 is the unique solution of equation (16). This is done in a manner almost identical to what was done above for the local problem: First, set η =Ê in equation (16) and take the real part to get This immediately shows that if condition (18b) holds, then the fields ⃗ E and ⃗ H are zero on Ω ⊂ R 3 and the proof is finished. In the case of condition (18a), we have ⃗ n×(Ê − ⃗ E ∂K ) = ⃗ 0 for all K. Using equations (10), this yields so Lemma 1 proves that the fields on element interiors are zero, which in turn implieŝ E = ⃗ 0 also. Thus, the theorem holds for the Maxwell case. The proof for the Helmholtz case is entirely analogous.
Note that even with Dirichlet boundary conditions and real k, the theorem asserts the existence of a unique solution for the Helmholtz equation. However, the exact Helmholtz problem (1) is well-known to be not uniquely solvable when k is set to one of an infinite sequence of real resonance values. The fact that the discrete system is uniquely solvable even when the exact system is not, suggests the presence of artificial dissipation in HDG methods. We will investigate this issue more thoroughly in the next section.
However, we do not advocate relying on this discrete unisolvency near a resonance where the original boundary value problem is not uniquely solvable. The discrete matrix, although invertible, can be poorly conditioned near these resonances. Consider, for example, the Helmholtz equation on the unit square with Dirichlet boundary conditions. The first resonance occurs at k = π √ 2. In Figure 2, we plot the condition number σ max σ min of the condensed HDG matrix for a range of wavenumbers near the resonance k = π √ 2, using a small fixed mesh of mesh size h = 1 4, and a value of τ = 1 that satisfies (18). We observe that although the condition number remains finite, as predicted by the theorem, it peaks near the resonance for both the p = 0 and the p = 1 cases. We also observe that a parameter setting of τ = −î that does not satisfy the conditions of the theorem produce much larger condition numbers, e.g., the condition numbers that are orders of magnitude greater than 10 10 (off axis limits of Figure 2b) for k near the resonance were obtained for p = 1 and τ = −î. To summarize the caveat, even though the condition number is always bounded for values of τ that satisfy (18), it may still be practically infeasible to find the HDG solution near a resonance.

Results of dispersion analysis for real wavenumbers
When the wavenumber k is complex, we have seen that it is important to choose the stabilization parameter τ such that (18b) holds. We have also seen that when k is real, the stability obtained by (18a) is so tenuous that it is of negligible practical value. For real wavenumbers, it is safer to rely on stability of the (un-discretized) boundary value problem, rather than the stability obtained by a choice of τ .
The focus of this section is on real k and the Helmholtz equation (1). In this case, having already separated the issue of stability from the choice of τ , we are now free to optimize the choice of τ for other goals. By means of a dispersion analysis, we now proceed to show that some values of τ are better than others for minimizing discrepancies in wavespeed. Since dispersion analyses are limited to the study of propagation of plane waves (that solve the Helmholtz equation), we will not explicitly consider the Maxwell HDG system in this section. However, since we have written the Helmholtz and Maxwell system consistently with respect to the stabilization parameter (see the transition from (3) to (9) via (7)), we anticipate our results for the 2D Helmholtz case to be useful for the Maxwell case also.
The dispersion relation in the one-dimensional case. Consider the HDG method (2) in the lowest order (p = 0) case in one dimension (1D) -after appropriately interpreting the boundary terms in (2). We follow the techniques of [1] for performing a dispersion analysis. Using a basis on a segment of size h in this order u 1 = 1, φ 1 = 1,φ 1 = 1,φ 2 = 1, the HDG element matrix takes the form M = M 11 M 12 M 21 M 22 where The Schur complement for the two endpoint basis functions {φ 1 ,φ 2 } is then Applying this matrix on an infinite uniform grid (of elements of size h), we obtain the stencil at an arbitrary point. Ifψ j denotes the solution (trace) value at the jth point (j ∈ Z), then the jth equation reads In a dispersion analysis, we are interested in how this equation propagates plane waves on the infinite uniform grid. Hence, substitutingψ j = exp(îk h jh), we get the following dispersion relation for the unknown discrete wavenumber k h : Simplifying, .
This is the dispersion relation for the HDG method in the lowest order case in one dimension. Even when τ and k are real, the argument of the arccosine is not. Hence in general, indicating the presence of artificial dissipation in HDG methods. Note however that if τ is purely imaginary and kh is sufficiently small, (20) implies that Im(k h ) = 0. Let us now study the case of small kh (i.e., large number of elements per wavelength). As kh → 0, using the approximation cos −1 (1 − x 2 2) ≈ x + x 3 24 + ⋯ valid for small x, and simplifying (20), we obtain Comparing this with the discrete dispersion relation of the standard finite element method in one space dimension (see [1]), namely k h h − kh ≈ O((kh) 3 ), we find that wavespeed discrepancies from the HDG method can be larger depending on the value of τ . In particular, we conclude that if we choose τ = ±î, then the error k h h − kh in both methods are of the same order O((kh) 3 ).
Before concluding this discussion of the one-dimensional case, we note an alternate form of the dispersion relation suitable for comparison with later formulas. Using the half-angle formula, equation (20) can be rewritten as where c = cos(k h h 2).
Lowest order two-dimensional case. In the 2D case, we use an infinite grid of square elements of side length h. The HDG element matrix associated to the lowest order (p = 0) case of (2) is now larger, but the Schur complement obtained after condensing out all interior degrees of freedom is only 4 × 4 because there is one degree of freedom per edge. Note that horizontal and vertical edges represent two distinct types of degrees of freedom, as shown in Figures 5a and 5b. Hence there are two types of stencils.
For conducting dispersion analysis with multiple stencils, we follow the approach in [7] (described more generally in the next subsection). Accordingly, let C 1 and C 2 denote the infinite set of stencil centers for the two types of stencils present in our case. Then, we get an infinite system of equations for the unknown solution (numerical trace) valueŝ ψ 1,⃗ p 1 andψ 2,⃗ p 2 at all ⃗ p 1 ∈ C 1 and ⃗ p 2 ∈ C 2 , respectively. We are interested in how this infinite system propagates plane wave solutions in every angle θ. Therefore, with the ansatzψ j,⃗ p j = a j exp(î⃗ κ h ⋅ ⃗ p j ) for constants a j (j = 1 and 2), where the discrete wave vector is given by we proceed to find the relation between the discrete wavenumber k h and the exact wavenumber k. Substituting the ansatz into the infinite system of equations and simplifying, we obtain a 2 × 2 system F [ a 1 a 2 ] = 0 where F = 2 khτ 2 c 1 c 2 d 1 (4 τ +îkh) + 2 khτ 2 c 1 2 d 2 (4 τ +îkh) + 2 khτ 2 c 2 2 2 khτ 2 c 1 c 2 and, for j = 1, 2, Hence the 2D dispersion relation relating k h to k in the HDG method is To formally compare this to the 1D dispersion relation, consider these two sufficient conditions for det(F ) = 0 to hold: where k 1 = k cos θ and k 2 = k sin θ. (Indeed, multiplying (26) j by d j+1 (j = 1) or d j−1 (j = 2) and summing over j = 1, 2, one obtains a multiple of det(F ).) The equations in (26) can be simplified to which are relations that have a form similar to the 1D relation (23). Hence we use asymptotic expansions of arccosine for small kh, similar to the ones used in the 1D case, to obtain an expansion for k h j , for j = 1, 2. The final step in the calculation is the use of the simple identity Simplifying the above-mentioned expansions for each term on the right hand side above, we obtain as kh → 0. Thus, we conclude that if we want dispersion errors to be O((kh) 3 ), then we must choose a prescription that is not very useful in practice because it depends on the propagation angle θ. However, we can obtain a more practically useful condition by setting τ to be the constant value that best approximates ± 1 2î cos(4θ) + 3 for all 0 ⩽ θ ⩽ π 2, namely τ = ±î These values of τ asymptotically minimize errors in discrete wavenumber over all angles for the lowest order 2D HDG method. Note that for any purely imaginary τ , (27) implies that k h j is real if kh is sufficiently small, so Im(k h ) = 0, thus eliminating artificial dissipation. We now report results of numerical computation of k h = k h (θ) by directly applying a nonlinear solver to the 2D dispersion relation (25) (for a set of propagation angles θ). The obtained values of the real part Re k h (θ) are plotted in Figure 3a, for a few fixed values of τ . The discrepancy between the exact and discrete curves quantifies the difference in the wave speeds for the computed and the exact wave. Next, analyzing the computed k h (θ) for values of τ on a uniform grid in the complex plane, we found that the values of τ that minimize kh − k h (θ)h are purely imaginary. As shown in Figure 4, these τ -values approach the asymptotic values determined analytically in (30). A second validation of our analysis is performed by considering the maximum error over all θ for each value of τ and then determining the practically optimal value of τ . The results, given in Table 2, show that the optimal τ values do approach the analytically  kh Optimal τ , Optimal τ , Table 2. Numerically found values of τ that minimize kh − k h (θ)h for all θ in the p = 0 case. determined value (see (31)) of ±î √ 3 2 ≈ ±0.866î. Further numerical results for the p = 0 case are presented together with a higher order case in the next subsection.
Higher order case. To go beyond the p = 0 case, we extend a technique of [7] (as in [9]). Using a higher order HDG stencil, we want to obtain an analogue of (25), which can be numerically solved for the discrete wavenumber k h = k h (θ). The accompanying dispersive, dissipative, and total errors are defined respectively by Again, we consider an infinite lattice of h × h square elements with the ansatz that the HDG degrees of freedom interpolate a plane wave traveling in the θ direction with wavenumber k h . The lowest order and next higher order HDG stencils are compared in Figure 5. Note that the figure only shows the interactions of the degrees of freedom   corresponding to theφ variable-the only degrees of freedom involved after elimination of the ⃗ u and φ degrees of freedom via static condensation. The lowest order method has two node types (shown in Figures 5a-5b), while the first order method has four node types (shown in Figures 5c-5f). For a method with S distinct node types, denote the solution value at a node of the s th type, 1 ⩽ s ⩽ S, located at ⃗ lh ∈ R 2 , by ψ s, ⃗ l . With our ansatz that these solution values interpolate a plane wave, we have for some constants a s . Now, to develop notation to express each stencil's equation, we fix a stencil within the lattice. Suppose that it corresponds to a node of the t th type, 1 ⩽ t ⩽ S, that is located at ⃗ h. For 1 ⩽ s ⩽ S, define J t,s = { ⃗ l ∈ R 2 ∶ a node of type s is located at (⃗  + ⃗ l)h} and, for ⃗ l ∈ J t,s , denote the stencil coefficient of the node at location (⃗  + ⃗ l)h by D t,s, ⃗ l . The stencil coefficient is the linear combination of the condensed local matrix entries that would likewise appear in the global matrix of equation (17 This is the equation that we solve to determine k h for a given θ for any order (cf. (25)).
Results of the dispersion analysis are shown in Figures 3 and 6. These figures combine the results from previously discussed p = 0 case and the p = 1 cases to facilitate comparison. Here, we set k = 1 and h = π 4, i.e., 8 elements per wavelength. Figure 6 shows the dispersive, dissipative, and total errors for various values of τ ∈ C. For both the lowest order and first order cases, although the dispersive error is minimized at a value of τ having nonzero real part, the total error is minimized at a purely imaginary value of τ . This is attributed to the small dissipative errors for such τ . Specifically, the total error is minimized when τ = 0.87î in the p = 1 case. This is close to the optimal value of τ found (both analytically and numerically) for p = 0. This value of τ reduces the total wavenumber error by 90% in the p = 1 case, relative to the total error when using τ = 1.  dissip , and total error total 1 total for various τ ∈ C. Here, k = 1, h = π 4, and 1 disp , 1 dissip and 1 total denote the errors when τ = 1, respectively.
Comparison with dispersion relation for the Hybrid Raviart-Thomas method. The HRT (Hybrid Raviart-Thomas) method is a classical mixed method [2,3,15] which has a similar stencil pattern, but uses different spaces. Namely, the HRT method for the Helmholtz equation is defined by exactly the same equations as (2) but with these choices of spaces on square elements: V (K) = Q p+1,p (K) × Q p,p+1 (K), W (K) = Q p (K), and M (F ) = P p (F ). Here Q l,m (K) denotes the space of polynomials which are of degree at most l in the first coordinate and of degree at most m in the second coordinate. The general method of dispersion analysis described in the previous subsection can be applied for the HRT method. We proceed to describe our new findings, which in the lowest order case includes an exact dispersion relation for the HRT method.
In the p = 0 case, after statically condensing the element matrices and following the procedure leading to (25), we find that the discrete wavenumber k h for the HRT method satisfies the 2D dispersion relation where c j , as defined in (24), depends on k h j , which in turn depends on k h . Similar to the HDG case, we now observe that the two equations are sufficient conditions for (36) to hold. Indeed, if l j is the left hand side above, then l 1 (2c 2 2 + 1) + l 2 (2c 2 1 + 1) equals the left hand side of (36). The equations of (37) can immediately be solved: Hence, using (28) and simplifying using the same type of asymptotic expansions as the ones we previously used, we obtain as kh → 0. Comparing with (29), we find that in the lowest order case, the HRT method has an error in wavenumber that is asymptotically one order smaller than the HDG method for any propagation angle, irrespective of the value of τ . To conclude this discussion, we report the results from numerically solving the nonlinear solution (36) for k h (θ) for an equidistributed set of propagation angles θ. We have also calculated the analogue of (36) for the p = 1 case (following the procedure described in the previous subsection). Recall the dispersive, dissipative, and total errors in the wavenumbers, as defined in (33). After scaling by the mesh size h, these errors for both the HDG and the HRT methods are graphed in Figure 7 for p = 0 and p = 1. We find that the dispersive errors decrease at the same order for the HRT method and the HDG method with τ = 1. While (38) suggests that the dissipative errors for the HRT method should be of higher order, our numerical results found them to be zero (up to machine accuracy). The dissipative errors also quickly fell to machine zero for the HDG method with the previously discussed "best" value of τ =î √ 3 2, as seen from Figure 7.  Figure 7. Convergence rates as kh → 0

Conclusions
These are the findings in this paper: (1) There are values of stabilization parameters τ that will cause the HDG method to fail in time-harmonic electromagnetic and acoustic simulations using complex wavenumbers. (See equation (5) et seq.) (2) If the wavenumber k is complex, then choosing τ so that Re(τ ) Im(k) ⩽ 0 guarantees that the HDG method is uniquely solvable. (See Theorem 1.) (3) If the wavenumber k is real, then even when the exact wave problem is not wellposed (such as at a resonance), the HDG method remains uniquely solvable when Re(τ ) ≠ 0. However, in such cases, we found the discrete stability to be tenuous. (See Figure 2 and accompanying discussion.) (4) For real wavenumbers k, we found that the HDG method introduces small amounts of artificial dissipation (see equation (21)) in general. However, when τ is purely imaginary and kh is sufficiently small, artificial dissipation is eliminated (see equation (32)). In 1D, the optimal values of τ that asymptotically minimize the total error in the wavenumber (that quantifies dissipative and dispersive errors together) are τ = ±î (see equation (22)). (5) In 2D, for real wavenumbers k, the best values of τ are dependent on the propagation angle. Overall, values of τ that asymptotically minimize the error in the discrete wavenumber (considering all angles) is τ = ±î √ 3 2 (per equation (31)). While dispersive errors dominate the total error for τ =î √ 3 2, dissipative errors dominate when τ = 1 (see Figure 7). (6) The HRT method, in both the numerical results and the theoretical asymptotic expansions, gave a total error in the discrete wavenumber that is asymptotically one order smaller than the HDG method. (See (38) and Figure 7.)