Discontinuous Galerkin method for the Helmholtz transmission problem in two-level homogeneous media

In this paper, the discontinuous Galerkin (DG) method is developed and analyzed for solving the Helmholtz transmission problem (HTP) with the first order absorbing boundary condition in two-level homogeneous media. This whole domain is separated into two disjoint subdomains by an interface, where two types of transmission conditions are provided. The application of the DG method to the HTP gives the discrete formulation. A rigorous theoretical analysis demonstrates that the discrete formulation can retain absolute stability without any mesh constraint. We prove that the errors in \begin{document}$ H^{1} $\end{document} and \begin{document}$ L^{2} $\end{document} norms are bounded by \begin{document}$ C_{1}kh + C_{2}k^{4}h^{2} $\end{document} and \begin{document}$ C_{1}kh^{2} + C_{2}k^{3}h^{2} $\end{document} , respectively, where \begin{document}$ C_1 $\end{document} and \begin{document}$ C_2 $\end{document} are positive constants independent of the wave number \begin{document}$ k $\end{document} and the mesh size \begin{document}$ h $\end{document} . Numerical experiments are conducted to verify the accuracy of the theoretical results and the efficiency of the numerical method.

1. Introduction. The acoustic scattering problems are often involved in many practical applications, such as biomedical imaging [10], non-destructive testing [30], geophysics and radar detecting [13], fuzzy recognition and image processing. These problems describe a general physical process where the acoustic waves are forced to deviate from the original direction of propagation due to localized non-uniformities in the medium through which they pass. In fact, such a non-uniformity is caused by an impenetrable or penetrable scatterer. So far, there have been a large number of works devoted to investigating the two kinds of acoustic scattering problems in theoretical analysis and numerical approximation.
On one hand, the scattering of time-harmonic acoustic waves by an impenetrable bounded obstacle can be described as the Helmholtz problem. Lots of numerical methods have been widely used to discretize the Helmholtz equation with various types of boundary conditions. In general, the numerical performance of any finite element solution of the Helmholtz problem depends significantly on the wave number k. When k is very large, the mesh size h has to be sufficiently small for the scheme to resolve the oscillations. In [28], Ihlenburg and Babuška showed that the H 1 -error bound for the finite element solution contains a pollution term that is related to the loss of stability with large wave numbers, and then, Babuška et al. [7] addressed the question of whether it is possible to reduce the pollution effect. In order to eliminate or substantially reduce the pollution errors, various numerical methods have been developed in the literature for solving the Helmholtz problem. Feng and Wu [22] developed some interior penalty discontinuous Galerkin methods and proved the proposed discontinuous Galerkin methods are stable without any mesh constraint. We refer the reader to [1,2,3,6,9,11,12,16,18,23,25,29,35,36,39] and the references therein for a detailed account on other numerical approaches for the Helmholtz problem.
On the other hand, the scattering of time-harmonic acoustic waves by a penetrable bounded obstacle implies that the acoustic waves can propagate through the surface of the obstacle. This is generally called Helmholtz transmission problem (HTP), that consists of a pair of Helmholtz equations coupled with transmission conditions under different wave numbers [33]. The so-called transmission conditions express the continuity of the medium and the equilibrium of the forces acting on it. In fact, the HTP is an unbounded region problem. In the present numerical investigation, there exists several techniques of numerical approximation for the HTP. One of them is the coupling of the boundary element method (BEM) and finite element method (FEM), which decomposes the unbounded domain by introducing an artificial boundary containing the obstacle inside. Then the BEM is to solve the exterior scattering problem outside the artificial boundary while the FEM is to solve the interior one [17,26,27,31,32,33]. Concerning the coupling of BEM and FEM, Johnson and Nédélec firstly derived the significant result addressing the theoretical justification in [31]. By the other technique, one can enforce some artificial boundary condition on the auxiliary boundary to reduce the original problem to a boundary value problem, for instance, the Dirichlet to Neumann mapping [24] and the first order absorbing boundary condition. In this paper, we consider the HTP with the first order absorbing boundary condition using the discontinuous Galerkin method.
By investigating the existing references of the coupling of the BEM and FEM for HTP, we note that the interface is always assumed to be sufficiently smooth in R d . However, the requirement of smooth interface is often too strict to be satisfied in practical obstacle. Therefore, it is necessary and valuable to extend the smooth interface to the case of Lipschitz continuous interface. The discontinuous Galerkin (DG) methods have received a lot of attention and undergone intensive studies by many researchers [4,5,8,14,15,19,21,37,38]. In [15], the authors concluded that the DG methods have the following main advantages: highly parallelizable; well suited to dealing with complicated geometries; easily handled adaptivity strategies. It is also well known that the trial and test spaces of DG methods are very easy to construct and they can naturally handle inhomogeneous boundary conditions. Motivated by above discussions, we aim to investigate the Helmholtz transmission problem in two-level homogeneous medium and develop the DG method for the HTP. The focus of the paper is to establish the rigorous stability and error analysis. Compared with the work in [22], we extend the Helmholtz problem of impenetrable obstacle to the HTP. Since the Helmholtz operator is not a coercive elliptic operator, the difficult part of the theoretical analysis is to establish the stability estimates for the numerical solution. For this purpose, we will make use of a local version of the Rellich identity and follow the stability analysis for the PDE solutions given in [18,22,25].
The remainder of the paper is organized as follows. We first describe the classical Helmholtz transmission problem in Section 2. In Section 3, the DG method for the Helmholtz transmission problem is formulated. The stability result for the DG method is established in Section 4. In Section 5, the error estimates of the discrete solutions to the HTP in H 1 and L 2 norms are strictly derived. In Section 6, we present some numerical experiments to gauge the theoretical estimates and to test the performance of the proposed DG method.
2. Statement of the Helmholtz transmission problem. In this section, we introduce the Helmholtz transmission problem (HTP) in two-level homogeneous media. Let Ω ⊂ R d , d = 2, 3, be an open, bounded polygonal or polyhedral domain. As shown in Figure 1, the interface Γ separates the whole domain Ω into two connected disjoint subdomains Ω 1 and Ω 2 with the corresponding wave numbers k 1 and k 2 , respectively. Denote by Γ 0 the outer boundary of Ω. Thus, we can express that Ω = Ω 1 ∪ Ω 2 ∪ Γ, Ω 1 ∩ Ω 2 = ∅ and ∂Ω 1 = Γ 0 ∪ Γ, ∂Ω 2 = Γ. The unknown u 1 (respectively u 2 ) stands for the scattered field (resp. transmitted field) with given wave number k 1 (resp. k 2 ) in the domain Ω 1 (resp. Ω 2 ). Then, the Helmholtz transmission problem in two-level homogeneous media reads: where i = √ −1 denotes the imaginary unit and µ = . The normal vector n on Γ is oriented from Ω 1 to Ω 2 . On Γ 0 , n is taken to be the unit outward normal vector. In this paper, we consider the case of k 2 < k 1 . All terms f 1 , f 2 , g, g 1 and g 2 on the right-hand side are given functions. The equation (3) is known as the first order absorbing boundary condition [20] on Γ 0 . The equations (4) and (5) are the transmission conditions on Γ.
We now recall the definition of star-shaped domains [22].
where n Ω is the unit outward normal vector to ∂Ω. If c Ω is positive, Ω is said to be strictly star-shaped.
Throughout the paper, we assume that Ω and Ω 2 are strictly star-shaped domains. Under these assumptions the following stability estimate for the problem (1) -(5) was proved in [34].
Theorem 2.2. Assume that Ω 2 is star-shaped, k 2 < k 1 and g 1 = g 2 = 0. Then the solution of (1) -(5) satisfies where C is a constant independent of k, but dependent on the measures of Ω 1 and Ω 2 .
3. Discontinuous Galerkin approximation. To formulate the discontinuous Galerkin (DG) method, we first need to introduce some notations. The standard space, norm and inner product notations are adopted. In particular for any region Q, (·, ·) Q and ·, · ∂Q denote the L 2 -inner product on complex-valued spaces L 2 (Q) and L 2 (∂Q), respectively. When Q stands for the whole domain Ω, we drop the subscript in the inner products. Throughout the paper, C is used to denote a generic positive constant which is independent of the wave number k and the mesh size h. We use the expression A B to denote that A ≤ CB. We also use A B to denote that A and B are equivalent. In general, Re v and Im v denote the real and imaginary parts of v, respectively. Let T h be a family of quasi-uniform triangulations of the domain Ω which are aligned with the interface Γ, with characteristic mesh size h > 0. Assume that any element K ∈ T h cannot be cut by Γ. Then, the partition T h can be naturally divided into two sets denoted by T h,1 = T h ∩Ω 1 and T h,2 = T h ∩Ω 2 . Denote by h K the diameter of any K ∈ T h . We also define ∂K = set of all boundary edges of element K, If two elements K and K are neighbors and share one common side e, there are two traces of v along e. We define a jump and an average for v, and assume that the normal vector n e is oriented from K to K : For every e ∈ E B 0 , n e be the unit outward normal vector to Γ 0 , set [v] = {v} = v| e . If e ∈ E Γ , n e is oriented from Ω 1 to Ω 2 across the interface Γ, set . We define the broken Sobolev space equipped with the broken gradient seminorm: We introduce three bilinear forms J 0 (u, v), J 1 (u, v), L 1 (u, v) : H × H → R that penalize the jumps of the function values, the normal derivatives and the tangential derivatives values: where τ e be the unit tangential vector to the edge e and the parameters γ 0,e , γ 1,e and β 1,e are the corresponding penalty parameters. Let We now define the sesquilinear form a h (·, ·) and linear form l(·) : For any v ∈ H, we define mesh dependent norms/seminorms Then the variational formulation of the problem (1) -(5) is finding u ∈ H and For any K ∈ T h , let P 1 (K) denote the set of all linear polynomials on K. We define the DG approximation spaces V h and V h,Γ as 4. Stability estimate. In this section, we focus on the stability estimate for the scheme (18). The crucial ingredients of the analysis are to use a special test function (18) and to use the Rellich identity on the element [18]. The following lemma establishes some integral identities which play an important role in our analysis.
where x Ω denotes the point in the star-shaped domain definition for Ω andv is the complex conjugate of v.
Remark 2. The identity (20) can be viewed as a local version of the Rellich identity for the Laplacian ∆ [18]. Since V h ⊂ H, the identity (19), (20), (21) also hold for The next lemmas will give some estimates for the stability analysis.

Lemma 4.2.
For any u h ∈ V h , the following estimates hold where The real and imaginary parts of Therefore, by taking the real and imaginary parts of the equations (24), we get (22) and (23). (18). Then Proof.
Then, by multiplying k 2 and adding 2k 2 u h 2 L 2 (Ω) on the above equality, we obtain where v h = α · ∇u h , u h is a linear polynomial on K, and hence, v h ∈ V h . Testing (18) by v h ∈ V h and taking the real part of the resulting equation, we get Using the identity |a| 2 − |b| 2 According to (20) in Lemma 4.1, it yields The facts that u h is piecewise linear and Plugging (30), (31), (32), (33) and (34) into (29), we obtain (28).
Next, we will derive the stability estimate for scheme (18).
Then there exists a positive constant C sta such that where Proof. First, we rewrite (28) as following Re(1 − µ) ∂u h ∂n e , u h e + I 1 + I 2 + I 3 + I 4 + I 5 + I 6 + I 7 + I 8 .
Next, we bound each term on the right hand side of (37). For an edge e ∈ E I , let K e and K e denote the two elements in T h that share e. For an edge e ∈ E Γ , let K e and K e denote the two elements that share e and K e ∈ T h,1 , K e ∈ T h,2 . For an edge e ∈ E B 0 , let K e denote the element that has e as an edge and K e = ∅. Applying the Hölder inequality, the trace and Young inequalities, we obtain In view of (21) in Lemma 4.1 and the Young inequality, we have Since u h is piecewise linear on K, ∇u h is piecewise constant vector on K and ∇u h L 2 (e) ≤ Ch Note ∂v h ∂n e K = ∂u h ∂n e K , we have Using the Hölder and trace inequalities, we obtain Based on the inverse inequality and the following identities v h = α · ∇u h , α · ∇u h = α · n e ∂u h ∂n e + α · τ e ∂u h ∂τ e , we obtain .
The continuous and coercive properties for the sesquilinear form a h (·, ·) are given in the following lemma.
In addition, for any 0 < ε < 1, there exists a positive constant C such that Proof. From (13), the definition of a h (·, ·), for any v ∈ H, w ∈ H ∩ H 1 Γ (Ω), we have, Therefore, we obtain From (26) and (27), we have Using the Cauchy inequality, Young inequality and trace inequality, we obtain In addition, there exist a constant C on each element K such that Plugging (54) -(56) into (53), we get . Therefore, we have The proof is completed.
The Lemma 4.5 directly implies the well-posedness of the scheme (18).

Elliptic projection.
In this subsection, we introduce an elliptic projection and derive an error estimate for the projection. For any w ∈ H ∩ H 1 Γ (Ω), we define its elliptic projection w h ∈ V h,Γ by Letting u ∈ H be the solution of the problem (1) -(5) and u h ∈ V h be its elliptic projection. Then, the definition (58) implies the following Galerkin orthogonality Lemma 5.1. Suppose u ∈ H ∩ H 1 Γ (Ω) is the solution of the problem (1) - (5). Then there hold the following estimates: where λ = 1 + 1 γ0 . Proof. Letting u h be the P 1 -conforming finite element interpolation of u on the mesh T h and denote η h = u h − u h . From η h + u − u h = u − u h and the Galerkin orthogonality (59), we get Assuming ε = 1 2 , C > 1 2 in (49) and following (62), we get which together with the regularity for the solution (cf. Theorem 1 in [18] and Section 3 in [34]) and estimates of finite element interpolation (cf. [11,14]) gives The above inequality together with the fact that u − u h = u − u h − η h implies the estimate (60). Next, we use the Aubin-Nitsche's duality argument [11,14] to show (61). Consider the following auxiliary problem It shows that ω satisfies Testing equation (63) by u − u h and using (59), we get Then, from (60), (67) and the above inequality, we obtain (61). The proof is completed.

Error estimates.
In this subsection, we shall derive error estimates for DG approximation (18). This goal will be achieved by making use of the stability estimate in Theorem 4.4 and the projection error estimate in Lemma 5.1. Letting u ∈ H and u h ∈ V h be the solution of (1) -(5) and of (18), respectively. Define the error function e h = u − u h . Subtracting (18) from (17), the following error equation holds Letting u h be the elliptic projection of u. Defining η = u − u h , ξ = u h − u h , we have e h = η − ξ. From (59) and (68), we get The above equation implies that ξ ∈ V h is the solution of (18) with the source form f = −k 2 η and g = 0. Then by Theorem 4.4 and Lemma 5.1, we can obtain the following lemma.
Lemma 5.2. The difference ξ between the discrete solution u h and its elliptic projection u h satisfies the following estimate where C sta is defined in Theorem 4.4 and λ = 1 + 1 γ0 . Combining Lemma 5.1 and Lemma 5.2, we obtain the main result in this paper. Theorem 5.3. Let u ∈ H and u h ∈ V h be the unique solutions of (1) -(5) and of (18), respectively. Then there exists two positive constants C 1 and C 2 such that the following error estimates hold where C sta is defined in Theorem 4.4 and λ = 1 + 1 γ0 . Proof. It follows from Lemma 5.1 and Lemma 5.2 and the triangle inequality that The proof is completed.
Remark 3. The estimates of Theorem 5.3 hold for any h > 0. There are so-called preasymptotic error estimates [28] (i.e., for the mesh in the regime k 2 h ≥ 1). 6. Numerical experiments. In this section, we will provide some numerical experiments to gauge the theoretical results and to test the performance of the DG method for solving the problem (1) -(5) in two dimensions. To this end, the following exact solutions are given where J ν (z)(ν = 0, 1) are Bessel functions of the first kind. The corresponding f 1 and f 2 are where r = x 2 + y 2 . The datum g, g 1 , g 2 on the right-hand side can be determined by the equations (3) -(5). We consider the domain Ω as the unit regular hexagon centered at the origin, Ω 2 as a smaller regular hexagon whose edge length is half of Ω and Ω 1 = Ω \ Ω 2 . For any positive integer m, let T 1/2 m denote the regular triangulation that consists of 6 · 2 2m congruent and equilateral triangles of size h = 1/2 m . The triangulation of Ω 2 consists of 3 · 2 2m−1 triangles. See Figure 2 for a sample triangulation T 1/8 of Ω 1 and Ω 2 , respectively. The mesh generation and all computations are conducted in the MATLAB c environment.
We next test the effect of the large wave number for the error of the DG solution. Figure 7 plots the relative error of the DG solution in H 1 -seminorm of different k 1 , k 2 . In this figure, it can be seen that for a large wave number, the relative error remains to be unchanged for large h, then decays with h being decreasing and converges at the order of O(h). This indicates the presence of the pollution effect.  Set the penalty parameters as γ 0,e = 100, γ 1,e = 0.01 + 0.07i, β 1,e = 1 and k 1 = 100, k 2 = 90, h = 1/64. In Figure 8 and Figure 9, we plot the DG numerical solution u h,1 (u h,2 ) and the exact solution u 1 (u 2 ), respectively. It shows that the DG solution has the same shape as the exact solution although its amplitude is not very accurate. However, Figure 10 and Figure 11 show that when we set k 1 = 100, k 2 = 50, the DG solution u h,1 have different shape and amplitude with the exact solution u 1 .
Remark 4. When we design the numerical tests, we found that if k 1 and k 2 have bigger difference, the error will be bigger. Therefore, we set that the values of k 1 and k 2 are near in the numerical tests. And we will research that how the difference between k 1 and k 2 has affected the errors in our future work.
The traces of the DG solution with penalty parameters γ 0,e = 100, γ 1,e = −0.04+ 0.08i, β 1,e = 1 in xz-plane display the damping for waves. In Figure 12, we display   This indicates an excellent convergence of the DG method when the mesh is refined. Furthermore, we can obtain a fast convergence by using high order polynomials with  mesh refinement in the DG method. Figure 13 displays the traces plot for the case of k 1 = 100, k 2 = 90 and h = 1/64, 1/128 when piecewise quadratic polynomials are employed in the DG method. It shows that the computational result captures both the fast oscillation and the magnitude of the exact solution very well.
Appendix: The derivation processes of the DG formulation. In this section we present the detailed derivation processes of the DG formulation.
Testing (1), (2) by v ∈ H ∩ H 1 Γ (Ω) and using integration by parts, we get Then the variational formulation of the problem (1) -(5) is: finding u ∈ H and u 1 − u 2 = g 1 on Γ, such that where a h (·, ·) and l(·) are defined in (13) and (14), respectively. The three terms J 0 (·, ·), J 1 (·, ·), L 1 (·, ·) are the penalty terms, which penalize the jump of the function values, the jump of the normal derivatives and the jump of the tangential derivatives, respectively. The term J 0 is necessary to maintain the stability of the discontinuous Galerkin methods. In Remark 3.1(e) of the reference [22], the authors have demonstrated that without J 1 term and L 1 term in a h (·, ·) the methods of this paper are still stable and convergent, but under a more stringent mesh constraint. What's more, the above derivation processes of the formulation show that a h (·, ·) + ik ·, · Γ0 is a consistent discretization for −∆, i.e., J 0 (u, v) = J 1 (u, v) = L 1 (u, v) = 0 for all u ∈ H 2 (Ω), v ∈ H ∩ H 1 Γ (Ω).