Partial differential equations (PDE) with sign-changing coefficients are used to model physical phenomena of meta-materials, for example, wave transmission problems between classical materials and meta-materials. Therefore, there is an emerging interest in the development of numerical methods for such PDEs. However, the bilinear forms for problems with sign-changing coefficients are not coercive, so standard techniques relying on coercivity cannot be applied. The lack of coercivity indeed poses fundamental challenges to study the well-posedness of such PDEs, as well as the development of numerical methods.
A popular approach in the study of numerical methods for PDEs with sign-changing coefficients is the “-coercivity” approach. This approach utilizes a linear operator that recovers the coercivity of the bilinear form. The existence of such an operator is a natural assumption because the existence of is equivalent to the well-posedness of the PDE in the sense of Banach–Nečas–Babuška in the functional analysis framework . The application of the operator to develop numerical methods was first proposed in . Since then, the -coercivity approach has been used in the development of numerical methods for Poisson-type problems , Helmholtz-type problems , and Maxwell-type problems [2, 5, 4].
In numerical methods that use -coercivity, the operator cannot be used directly, because the domain and the range of are not the discrete spaces in general. For this reason, is introduced, a discrete approximation of , which has the coercivity recovery property in the discrete setting. The existence of such is non-trivial and difficult to guarantee because the form of is unknown in general. For general geometry and meshes sufficient conditions for the existence of are proposed in 
, but these include non-quantitative constants such as a norm of the discrete interpolation ofand the ratio of the positive and negative coefficients (the contrast). It is furthermore shown in  that locally symmetric meshes across the transmission interface, where the coefficient changes sign, improve the quality of numerical solutions. However, the generation of such meshes is a non-trivial restriction in mesh generation for complex interfaces.
Hybridization (or static condensation) was originally introduced for mixed finite element methods to reduce computational cost , and it was further analyzed in . More recently, hybridization was combined with the discontinuous Galerkin (DG) method. This resulted in the hybridizable discontinuous Galerkin (HDG) method which was shown to be more efficient than standard DG methods. The HDG method was systematically presented in  for elliptic partial differential equations. Since its introduction, it has been extended to many different problems. These include HDG methods for the Helmholtz problem, e.g., [12, 17] and the Maxwell problem, e.g., [6, 14].
In this paper, we develop a numerical method for Poisson-type problems with a sign-changing coefficient that avoids the -coercivity approach. The numerical method we present in this paper is always well-posed without any conditions on the domain geometry or the ratio between the negative and positive coefficients. This is achieved by employing two key ideas for the problem: a mixed formulation of the problem and using a hybridized discontinuous Galerkin (HDG) method with sign-changing stabilization parameters. We will see that the sign-changing stabilization parameter and the discontinuous test function space in discontinuous Galerkin methods allow us to obtain stable numerical methods without any non-quantitative assumptions on mesh size or coefficients.
We also carry out an error analysis of the HDG method for this problem. The error analysis of HDG methods typically uses the Aubin–Nitsche duality argument with the full elliptic regularity assumption. However, this assumption is not feasible in our setting, because the PDE is not elliptic. By revealing a stabilized saddle point problem structure of HDG methods, we circumvent this obstacle and we can show an error estimate without using the duality argument. We note that such analysis has been applied also to HDG methods for standard Poisson-type problems to avoid the full elliptic regularity assumption in the a priori error analysis (see).
The outline of this paper is as follows. In Section 2 we introduce the Poisson problem with a sign-changing coefficient. We introduce the HDG method and discuss its well-posedness in Section 3. An a priori error analysis is presented in Section 4. In Section 5 we show that a superconvergence result can be obtained when a suitable regularity assumption is available. Our analysis is verified by numerical examples in Section 6, while conclusions are drawn in Section 7.
2 The Poisson problem with a sign-changing coefficient
Let , , be a polygonal domain that is divided into two subdomains and such that and . The boundaries of , and are denoted by, respectively, , and . The interface separating the domains is denoted by and we define furthermore and . We assume that and are disjoint subsets of such that
. The outward unit normal vector field onis denoted by .
Let be a scalar function defined as
where and are constants. The contrast is defined as . Throughout this paper we assume that
We consider the following partial differential equation (PDE) for the scalar : equationparentequation
where and are given boundary data and is a given source term. This is not a standard second order PDE, because the sign of is indefinite. Let where is the standard Sobolev space of functions such that the -norm of their gradient is bounded. The variational formulation for the problem (when ) is given by: Find such that
3 The HDG method
Introducing the auxiliary variable , Eq. 3 can be written as a system of first-order equations: equationparentequation
In this section we introduce the HDG method for Eq. 5.
Let be a triangulation of and be the set of faces in the triangulation with co-dimension 1. By we denote the subset of such that its elements are in . is defined similarly. We assume that is a conforming triangulation with respect to , i.e., there is a subset of such that . For later reference we define and as the subsets of interior facets of such that the facets are in and , respectively. We also define as the facets on , i.e., . In summary, is a disjoint union of , , , . We also define as the subset of such that .
Let , denote the diameter of by , and let . For scalar functions and vector functions we denote, respectively, and . Furthermore, for two functions and with well-defined traces on , we define , where is the unit outward normal vector field on . Additionally, we define , , and , with similar definitions on and .
We use the standard notation for Sobolev spaces based on the -norm, i.e., , , denotes the Sobolev space based on the -norm with -differentiability on the domain . We refer to  for a rigorous definition of this space. The norm on is denoted by . When we drop the subscript .
To define the HDG method let
where and are the spaces of scalar and -dimensional vector valued polynomials of degree at most on a domain . We will use the following finite element spaces:
3.2 The discrete formulation
For a sufficiently regular solution of the mixed problem Eq. 5, and for the restriction of on , we can derive a system of variational equations from Eq. 5 with test functions in as equationparentequation
Here we use to denote a piecewise constant function adapted to ; in the integration on is the trace of . Our HDG method for Eq. 5 is the discrete counterpart of this system of variational equations, i.e., we find such that equationparentequation
Here equations Eq. 9a, Eq. 9b and the discrete equations Eq. 10a, Eq. 10b look inconsistent. This is because we impose the restriction on so that we can take the same trial and test function spaces in our numerical method. This choice will be useful to reveal a stabilized saddle point structure of our numerical method, and also allow us to obtain optimal error estimates without the Aubin–Nitsche duality argument. To clarify the consistency between Eq. 9a, Eq. 9b and Eq. 10a, Eq. 10b, we point out that Eq. 10a, Eq. 10b with are equivalent to the discretized forms of Eq. 9a, Eq. 9b with defined by
For later use we define
Then the sum of the left-hand sides of Eq. 10 can be rewritten as
For brevity, but without loss of generality, we assume and in the remainder of this paper.
In the following lemma we prove well-posedness of the HDG method Eq. 10.
Lemma 1 (Well-posedness)
Let . If on and on , there exists a unique solution to Eq. 10.
To show well-posedness of Eq. 10 assume that . We will show that , , vanish.
be the characteristic function that has the value 1 in domainand 0 elsewhere. Then let , , in Eq. 10, and add all the equations. The evaluation of Eq. 10 with the first components of these test functions gives
Similarly, the evaluation of Eq. 10 with the second components of the above test functions gives
Since and are positive and negative on and , respectively, we can conclude that on and on . This also implies that is continuous on .
Since and it follows from Eq. 10a that , i.e., is constant on . Then on because on .
We consider the following norms on and :
where is the absolute value of . Additionally, we define with norm
Furthermore, in the remainder of this paper we let be a constant independent of the mesh size .
There exists which depends only on and such that for there exists satisfying and with on .
We will use a result in [16, p.176] which claims that for any given and satisfying , there exists such that , on , and .
We first note that if has mean value zero, then there exists such that and by taking in the above result.
Consider now a decomposition with where is the Lebesgue measure of , and where is the component of with mean-value zero. It is clear that . As just noted, there exists such that and . Since is a polygonal/polyhedral domain, one can find a function whose support resides on , , and . We can find , a renormalization of , satisfying , and it holds that
with a constant that depends on and . Applying the result in [16, p.176] there exists such that , on , and . Then is a function in satisfying the desired conditions.
Theorem 1 (inf-sup condition)
Suppose that . If on and on , and for every , with positive constants and independent of , then there exists a positive constant independent of such that
We first prove the following weak inf-sup condition: There exist independent of such that
We show this by proving an equivalent condition, i.e., there exist independent of such that for any there exists such that and
By Lemma 2 there exists with such that and with depending only on and . We use to denote the element-wise orthogonal projection of into . It then holds that
for all . Then
where the third equality holds because and are single-valued on facets, on by the definition of , and on .
The element-wise trace inequality and the Cauchy-Schwarz inequality give
holds, so the weak inf-sup condition Eq. 15 follows.
4 Error analysis
In this section we present the a priori error estimates of our HDG method. For this we recall an interpolation result (cf. [10, Theorem 2.1]): For suppose that is a piecewise constant function which is nonnegative if and nonpositive if such that . Then there exists a bounded interpolation operator from to such that equationparentequation
for and , and
hold with . Here,
where is the face where is attained, and the implicit constants are independent of and . We remark that the proof in  is only for nonnegative . The proof, however, can easily be adapted to nonpositive , so we omit its proof here. We will also require , the facet-wise projection from to .
For brevity of notation, we will denote the difference between an unknown and its approximation by . It will be convenient to also split the error into interpolation and approximation errors: