Two scale Hardy space infinite elements for scalar waveguide problems

We consider the numerical solution of the Helmholtz equation in domains with one infinite cylindrical waveguide. Such problems exhibit wavenumbers on different scales in the vicinity of cut-off frequencies. This leads to performance issues for non-modal methods like the perfectly matched layer or the Hardy space infinite element method. To improve the latter, we propose a two scale Hardy space infinite element method which can be optimized for wavenumbers on two different scales. It is a tensor product Galerkin method and fits into existing analysis. Up to arbitrary small thresholds it converges exponentially with respect to the number of longitudinal unknowns in the waveguide. Numerical experiments support the theoretical error bounds.


Introduction
This paper deals with the efficient numerical solution of time-harmonic scattering and resonance problems which exhibit wavenumbers on different scales. Such phenomena might occur e.g. in nearly incompressible elastic materials. Another example are guided waves for which the geometry causes a dispersion effect, i.e. waves can propagate with different speed levels. At so-called cut-off frequencies (also referred to as Wood's anomalies) there exist wavenumbers which change from being genuine complex to being purely real, most often through the origin. This situation causes severe performance issues for non-modal numerical methods like the perfectly matched layer or the Hardy space infinite element method. In this paper we propose and analyze an improved variant of the latter.
For the sake of simplicity we restrict the presentation to a scalar Helmholtz equation in an infinite cylindrical waveguide with sound-hard boundary conditions. The extension to problems with several waveguides [18] or to cylindrical elastic waveguides [15] is straightforward as long as the underlying radiation condition remains valid.
Analytical studies of solutions to waveguide problems are typically based on the waveguide modes, see e.g. [24]. For a fixed frequency ω > 0 the waveguide modes satisfying the Helmholtz equation have longitudinal wavenumbers κ n (ω) = ± ω 2 − λ n , while λ n denote for n ∈ {0, 1, 2, . . . } the eigenvalues of the negative Laplace operator on the cross section of the waveguide. There exist finitely many propagating modes with real wavenumbers and infinitely many evanescent or exponentially increasing modes with non-real wavenumbers. In order to avoid unphysical solutions, exponentially growing as well as propagating modes with negative group velocity ∂ω/∂κ n are sorted out by the radiation condition.
Some numerical methods are directly related to the waveguide modes, see e.g. [22,26,27]. Since the wavenumbers and the waveguide modes depend non-linearly on the frequency, these approaches result for resonance problems into the solution of a non-standard eigenvalue problem. Frequency independent methods like the perfectly matched layer method [1,16,18] or the Hardy space infinite element method [12,17,18] lead to linear eigenvalue problems. But the performance of these methods deteriorates in the vicinity of cut-off frequencies, where an evanescent mode becomes a propagating mode.
For such frequencies ω with ω 2 ≈ λñ for oneñ, there exist wavenumbers on highly different scales, since κñ(ω) vanishes for ω 2 = λñ. This leads to a loss of performance e.g. for the perfectly matched layer method, where the damping profile is typically adjusted to one wavenumber. In [3] a hybrid perfectly matched layer/modal based method was proposed in order to reduce the error of the perfectly matched layer method arising from waveguide modes with small longitudinal wavenumbers. However, the method depends again non-linearly on ω 2 .
For the standard Hardy space infinite element method such problems exist as well, see [18] for scalar problems and [15] for elastic problems. We will use in this paper a kind of two scale Hardy space method, which allows to adjust the method to wavenumbers on two different scales. In [14] this method was introduced for problems where the signs of phase velocity ω/κ and group velocity ∂ω(κ)/∂κ differ. For such problems the use of the method introduced in [14] was essential in order to guarantee physically correct solutions. In the scalar waveguide problem under consideration in this paper, phase and group velocities always have the same sign. So in terms of the radiation condition there is no need for such a modified method.
Nevertheless, we can use it as a two scale Hardy space method leading to a significant improvement in the vicinity of cut-off frequencies.
The rest of the paper is organized as follows: in Section 2 we formulate the scattering as well as the resonance problem. This includes the underlying radiation condition, the pole condition, and the variational framework of the Hardy space infinite element method. Section 3 contains the two scale Hardy space infinite elements and their performance in the presence of wavenumbers on different scales. References to the convergence analysis for scattering as well as resonance problems are given in Section 4. In Section 4.3 we improve the previously derived convergence rate by focusing on the error in the area of interest. To this end we study the approximated Dirichlet-to-Neumann operator, which is implicitly determined by the Hardy space method. The numerical experiments in Section 5 support the convergence results. The paper closes with a short conclusion.

General setting
Let int ⊂ R d and ϒ ⊂ R d−1 with d = 2, 3 be bounded Lipschitz domains,Υ := {0} × ϒ ⊂ ∂ int and ext := R >0 × ϒ be an infinite cylindrical waveguide with int ∩ ext =Υ. For given frequency ω ∈ C we will look for solutions u to the Helmholtz equation Since is unbounded, (2.1) needs to be equipped with boundary and radiation conditions to obtain a meaningful problem.
Hence, it holds u n (ξ ) = c n exp(iκ n (ω)ξ ) + d n exp(−iκ n (ω)ξ ), ω 2 = λ n c n + d n ξ, ω 2 = λ n , ξ >0 (2.3) with wavenumbers κ n (w) = ω 2 − λ n and constants c n , d n ∈ C. For real ω ≥ 0 we use the definition κ n (ω) := i λ n − ω 2 , ω ∈ [0, √ λ n ] ω 2 − λ n , ω > √ λ n . (2.4) u n is called outgoing if d n = 0 for all n ∈ N 0 . This definition can be motivated by the limiting absorption principle or by energy flux arguments (see e.g. [18]). Since this definition of outgoing solutions directly uses the waveguide modes, we call it a modal radiation condition. A radiation condition for complex ω has to be a holomorphic extension of the previous radiation condition for ω ∈ R >0 \ n∈N 0 { √ λ n } in order to be meaningful. Thereby the precise definition of holomorphicity in this context is kept vague. Remark that (2.4) is itself not analytic for ω 2 = λ n , n ∈ N 0 , and hence a requirement to be a holomorphic extension of the whole positive half-axis is not meaningful.
A possible way to construct such a radiation condition is to holomorphicly extend The domain of holomorphicity of such a construction will always exclude branch cuts starting from each √ λ n , n ∈ N 0 . For fixed θ ∈ [0, 2π) we may use the definition Since this definition coincides with (2.4) if and only if {0, π} ∈ (θ − 2π, θ), we obtain a holomorphic extension of (2.4) for all θ ∈ (π, 2π).

Remark 2.1
For the following we need to assume, that the sequence (|c n |) n∈N 0 in (2.3) decays fast enough, which is guaranteed if u| ϒ is sufficiently smooth. E.g.
(|c n |) n∈N 0 decays super-algebraically, if we assume that u| ϒ ∈ H k (ϒ) for all k ∈ N. This is always the case, if the domain is sufficiently smooth aroundΥ, if there is no change of boundary conditions in the neighborhood ofΥ and if there are no sources in the neighborhood ofΥ. One practical remedy for problems not fulfilling these conditions would be to start the exterior domain not exactly at the beginning of the waveguide. Using Weyl asymptotics λ n = O(n 2 d−1 ) (see e.g. [25] and the references therein) on the interface ϒ ⊂ R d−1 and by definition of κ n with θ ∈ (π, 2π), (iκ n (ω)) = O(− √ λ n ) holds for sufficiently large n ∈ N. Hence, for ξ > 0 the sequence u n (ξ ) = c n exp(iκ n (ω)ξ ) with Fourier coefficients c n := u(0, •), ϕ n L 2 (ϒ) decays exponentially (d = 2) or at least super algebraically (d = 3) for n → ∞ . If a part of the waveguide with length L is added to the interior domain, the new constantsc n := u(L, •), ϕ n L 2 (ϒ) are given by c n exp(iκ n (ω)L) and decay sufficiently fast.
Two different θ 1 , θ 2 ∈ (π, 2π) lead to different functions u for a complex frequency ω, if the signs of κ θ 1 n (ω) and κ θ 2 n (ω) differ at least for oneñ ∈ N 0 . The choice of θ is implicitly realized by the complex scaling radiation condition as well as by the pole condition (see [18]). In this paper we employ the latter. Roughly speaking the pole condition specifies regions in the complex plane in which the Laplace transform of the solution should have no singularities. In full details and for scalar waveguide problems this radiation condition is presented in [18,Sec. 6.1]. We remark that we don't distinguish between the Laplace transform of a function and its holomorphic (or meromorphic) extension to guarantee a smooth presentation. The Laplace transform u n (s) := (Lu n )(s) = ∞ 0 u n (ξ ) exp(−sξ )dξ of u n is given bŷ Soû n are meromorphic functions. For κ n = 0 they have two poles of order one at ±iκ n . In the case of a cut-off frequency, i.e. if κ n = 0, these poles collapse to one pole of order two. Figure 1 shows the location of the poles in the complex plane for an ω with λ 2 < ω 2 < λ 3 . Now, let s 0 ∈ C be a complex parameter with negative real and positive imaginary part. We decompose C into the oriented straight line := −is 0 R and the half spaces below − and above + (see Fig. 1). The pole condition states, that for all η ∈ ϒ the Laplace transformû(s, η) := (Lu(•, η))(s) should be well defined and have no poles p n := iκ n in − ∪ . This radiation condition determines uniquely the sought holomorphic extension of the modal radiation condition using iκ n (ω) ∈ + as definition of κ n . More explicitly, this definition of κ n coincides with the choice θ = 2 arg(s 0 ) in (2.5).
Analogy to PML Perfectly matched layer methods for waveguide structures can be interpreted as a complex scaling of the longitudinal direction of the waveguide, see e.g. [1,7,16,18]. More precisely, a complex local coordinateξ(ξ) is introduced such that the complex scaled functionsû n := u n •ξ are exponentially decaying for ξ → ∞. For a linear complex scaling with σ ∈ C, complex scaling is a rotation and stretching of iκ n by σ (confer with Fig. 1). Using (iσ κ n (ω)) < 0 for all n ∈ N 0 as definition of κ n is equivalent to the choice θ = arg(1/σ 2 ) in (2.5). So the pole condition with parameter s 0 is equivalent to a linear complex scaling radiation condition with σ = 1/s 0 .
In order to turn the pole condition into a variational framework, we define the Hardy Those ω ∈ C belonging to the set defined in (2.6) lead to poles iκ n ∈ . This is in particular a problem for cut-off frequencies, i.e. if ω = √ λñ for oneñ ∈ N 0 . The Laplace transformûñ has in this case a pole at the origin and therefore does not belong to H − ( ) ⊂ L 2 ( ).
For the following theory we have to assume iκ n (ω) ∈ + , but we will study in details the performance of the presented method in the neighborhood of cut-off frequencies ω ≈ √ λ n . In this case the waveguide is a paradigm for problems with wavenumbers on different scales.

Hardy space variational framework
In [18,Sec. 6.2] the Hardy space infinite element method based on the pole condition is presented in detail for scalar waveguides. There, the Hardy space H − ( ) is mapped unitary to the Hardy space H + (S 1 ) of the complex unit disk using a Möbius transformation. The variational framework is then build in H + (S 1 ). Since we want to use later on a slightly different discretization, we reformulate the framework using directly the Hardy space H − ( ) of the complex half plane. In one dimension this approach is presented in [14,Sec. 3]. So all proofs of the following results can be found in one of these references.
Since the waveguide ext is the Cartesian product of R >0 with the interface ϒ, our variational framework is based on tensor products. More precisely, we use the transversal Hilbert spaces Y 1 := L 2 (ϒ) and Y 2 := H 1 (ϒ). These space are combined with longitudinal Hilbert spaces X 1 , X 2 to the exterior Hilbert space with scalar product In this way, the one-dimensional longitudinal direction of the cylindrical waveguide can be treated separately from the transversal directions. While the variational framework for the transversal direction is rather standard, the one for the longitudinal direction will be non-standard for the Hardy-space method.
Analogy to PML In [18] the variational framework of the Hardy space method is treated simultaneously with the one for a linear complex scaling method as in (2.8).
For the complex scaling method the longitudinal spaces are the standard Sobolev spaces X 1 := L 2 (R >0 ) and X 2 := H 1 (R >0 ). Hence, the exterior space V ext is isometrically isomorphic to The main ingredient to derive a variational formulation posed in a Hardy space is the integral transformation with the bilinear form This transformation holds for suitable functions f and g, in particular for f = u n (u n outgoing) and test functions g = v with v as in [14,Lemma 3.3]. To build a functional framework it is necessary to cope with the formal trace and differential operatorsf for suitable functionsf ∈ H − ( ). To this end we introduce the Hilbert space with scalar product and the two linear operators T m , T s : (2.14) These operators are bounded and injective. Furthermore, Lu n , Lv ∈ T m (X 2 ) and it holds In particular, using (2.7) we compute for η ∈ ϒ and with bilinear forms Note, that there is a slight change in notation: The boundary value u 0 := u| ϒ as well as v 0 , f 0 , g 0 are no longer constants in C but functions depending on the transversal variable η. Next we construct the canonical space on which the bilinear forms a ext and b ext are bounded. To this end let X 1 be the completion of X 2 with respect to the scalar product As mentioned in the beginning of this section, we use the standard Sobolev spaces with standard scalar products for the transversal directions and define the exterior space V ext by (2.10a). The scalar product (2.10b) of V ext takes the form (2.20) Here we have used the notation In order to complete the framework, using direct calculations it can be shown that L ξ u ∈ V ext for outgoing solutions to the Helmholtz equation. Moreover, the transformed test functions are dense in V ext by [14,Lemma 3.3]. Hence, we can drop the restrictions on the test functions since the bilinear forms are bounded on V ext . Finally, it can be shown that Analogy to PML For a linear complex scaling method (2.8) the exterior bilinear forms (2.18b) are Note, that this is still an unbounded perfectly matched layer.

Scattering and resonance problem
We are now ready for a precise formulation of the scattering and resonance problems. Let the boundary of the domain be the disjoined union of the boundary of the waveguide ∂ ext \Υ, a NeumannΥ N and a Dirichlet partΥ D of ∂ int \Υ. Then one possible setting of boundary conditions is with a given g ∈ H 1/2 (Υ D ) and the normal derivative ν. So in total for given ω > 0 with ω 2 = λ n for all n ∈ N 0 we are looking for a solution u ∈ H 1 loc ( ) 2 to the Helmholtz equation (2.1), which satisfies the radiation condition (2.9) and the boundary conditions (2.21).
To express the variational formulation we introduce the interior bilinear forms defined on with scalar product Since the scalar product does not depend on g, we simply use V as index for the scalar product. For simplicity, due to trΥ f = f 0 we denote occasionally elements of V g with f F and omit the explicit notation of f 0 . So finally we are looking for with bilinear forms Remark 2. 2 We use here a very simple setting. Roughly speaking all kinds of generalizations are possible, as long as they change only the interior Helmholtz problem and fit into the standard H 1 ( int ) framework. This includes inhomogeneous media, source terms and different kinds of inhomogeneous boundary conditions on ∂ int \Υ.
For the exterior problem we might use homogeneous Dirichlet boundary conditions on the waveguide boundary instead of (2.21a). This would just change the eigenvalues λ n to these of the Dirichlet-Laplacian. Also incoming waves entering via the waveguide ext can be incorporated using jumps in the Cauchy data atΥ.
Some of these generalizations can be found in [18]. Since they have no influence over the Hardy space infinite element method, we keep the presentation in this paper as simple as possible.
Based on the variational formulation (2.23) of Scattering problem (2.1), (2.21), (2.9) it is straightforward to formulate the variational formulation of the corresponding resonance problem: We are looking for resonances ω ∈ C (with positive real part) and non-trivial resonance function u U ∈ V 0 such that The well-posedness of (2.23) and (2.25) will be discussed in Section 4.1.

Hardy space infinite elements
We perform a standard Galerkin method using a parametrized family (V 0 . For resonance problems we solve for eigenpairs ω t , The convergence analysis of both problems is given in Section 4. In principle, for both problems the discretization error is bounded by a constant times the best approximation error (3.1) In this section we first present the discretizations and study the corresponding best approximation errors.
For the longitudinal part of the exterior space let X (N) ⊂ X 2 with increasing N ∈ N be a second family of (N In the next subsections we present different choices for these longitudinal subspaces and study how they perform in the vicinity of cut-off frequencies. Based on these spaces we employ subspaces for the exterior space. Finally, we obtain a family of finite dimensional subspaces of V defined by Since we employed the traces of the interior spaces V int 0,h as transversal part of V ext h,N , the coupling trΥ v = v 0 between exterior and interior spaces is straightforward. to g is necessary which introduces an additional error. Since the corresponding analysis is standard we henceforth always assume g ∈ trΥ N) in V is shown. Hence, the best approximation error (3.1) with t = (h, N) converges to zero if h → 0 and N → ∞.
Analogy to PML The same tensor product element can be performed for a complex scaling method, see [18]. Since the longitudinal space X 2 = H 1 (R >0 ) for the PML is unbounded, typically a truncation of R >0 to a bounded layer [0, L] is used. Using homogeneous Dirichlet boundary conditions at the end of the layer L leads to a projection of functions in X 2 to functions with compact support in [0, L]. Since such a projection converges point wise, the convergence theory in Section 4.2 applies to the truncated PML as well. Nevertheless, a discretization of the finite layer is needed additionally. So overall, we get a triple of discretization parameters Remark 3.1 In [2] an exact complex scaling is proposed, where a truncation of the infinite layer is not needed. This can be achieved by replacing the linear complex scaling (2.8) with a more general scaling of the formξ = σ (ξ) such that lim ξ →L (ξ) = ∞ for a finite L > 0. Please note, that such scalings have to be used very carefully for waveguide problems.
The same holds true for the real part of κ n if (ξ)/ (ξ ) → 0. Therefore, for resonance problems in waveguides real and imaginary part should be scaled at least asymptotically proportional.
To derive convergence rates, let (u, U ) ∈ V be a solution to (2.23) N) and Then by definition of · V ext via (2.20) and the triangle inequality it holds ).
( 3.4) Only the last term depends on both discretization parameters h and N. Thus we introduce a U (N) Recall that U n is a meromorphic function with one pole at iκ n due to (2.16).
The construction of suitable spacesX to minimize (3.5) will be the main goal of the following subsections. All other terms in (3.4) including the additional term can be treated with standard finite element interpolation arguments, see e.g. [6,Sec. 4.4] or [8, Theorem 3.2.1].

Two pole basis functions
The original Hardy space infinite element method [12,17,18] employs the Hardy space H + (S 1 ) of the complex unit sphere. Since F ∈ H + (S 1 ) is the boundary value of functions, which are holomorphic in the complex unit disk, it can be written as The monomials z → z j , j ∈ N 0 , build an orthogonal basis of H + (S 1 ) with respect to the standard L 2 (S 1 ) scalar product. Since the Möbius transform maps bijectively ∪ − to {z ∈ C : |z| ≤ 1} and H + (S 1 ) to H − ( ) respectively, it leads to the orthogonal basis of H − ( ), see [14,Lemma A.4]. Since all basis functions are meromorphic with one pole at s 0 , we call it a one pole basis. In [14,Sec. 4] this basis was generalized to a two pole basis. Let s 1 ∈ s 0 R >0 be a second parameter. Then Remark 3.2 In [14] this basis was used in order to construct a simple Riesz basis of Hardy spaces H − ( ) with curved depending on s 0 and s 1 . Such Hardy spaces are useful e.g. for cylindrical elastic waveguides [15], where the signs of phase velocity ω/κ and group velocity ∂ω(κ)/∂κ differ. In this paper we always assume s 1 ∈ s 0 R >0 and = −is 0 R. Thus, the two pole functions Depending on the two poles we define the function with which we will characterize convergence rates. Note, that this is the absolute value of the reciprocal of the corresponding function in [14,Sec. 4.2]. With the definitions of Section 2.1 it is straightforward to see, that We analyze in the following the best approximation error (3.5) using the first N basis functions the right hand side of (3.5) is due to (3.9) bounded by . (3.13) In order to study this bound, we need some knowledge about the wavenumbers κ n , the Fourier coefficients c n and the transversal eigenvalues λ n . For large n it holds So the exponential decay of g s 0 ,s 1 (iκ n ) N in the number of basis functions N decelerates for increasing n. Nevertheless, the overall performance of the method suffers only mildly from this degeneration of rates for large n. Using the presupposition of Rem. 2.1 on the decomposition of , there exist constants C 3 , C 4 > 0 such that |c n | ≤ C 3 exp −C 4 √ λ n for sufficiently large n ≥ n 0 ∈ N 0 . For each n ≥ n 0 the deteriorating rate g s 0 ,s 1 (iκ n ) N in the number of basis functions N is multiplied with constants (1+λ n )|c n | 2 1−g s 0 ,s 1 (iκ n ) ≤ C 3 exp −C 4 √ λ n , which decrease at least super-algebraically in n. By the integral test there exist constants C, c > 0 independent of n 0 such that Therefore, using g s 0 ,s 1 (iκ n ) N < 1 for n ≥ n 0 the best approximation error in the Hardy space, which is the square root of (3.13), is bounded for each sufficiently large Hence, the best approximation error of the Hardy space infinite elements decays, up to an arbitrary small threshold exp(−c λ n 0 ), exponentially for increasing number of degrees of freedom N. The rate depends on g s 0 ,s 1 (iκ n ) < 1 for the first n 0 poles iκ n of U . Since = {s : g s 0 ,s 1 (s) = 1}, the convergence rate degenerates for fixed s 0 and s 1 if iκ n ∈ + approach . In particular this would happen for ω → √ λñ for onẽ n ∈ N 0 , since then iκñ(ω) → 0 ∈ . In contrast to the degeneration of rates for large n, in this case the constant cñ might not be negligible. Thus, it is meaningful to adjust the scales of the parameters s 0 and s 1 to the scales of the (low indexed) wavenumbers.
Analogy to PML The one-pole basis corresponds more or less to a linear complex scaling with one constant damping parameter σ and a kind of optimized infinite element discretization of the infinite layer. The two-pole basis can be most likely compared to a perfectly matched layer method, where simultaneously two different dampings would be used.
Note, that for fixed σ the truncation error of the PML method converges exponentially for L → ∞ with rate exp( (iκ n σ )L). So the linear complex scaling suffers from wavenumbers on different scales as well. Moreover, (3.14) contains the whole discretization error of the Hardy space infinite element method. For the PML typically only the truncation errors decays exponentially, while the discretization error of the finite layer decays with the usual finite element approximation errors.
The dependency of |c n |g s 0 ,s 1 (iκ n (ω)) N/2 with respect to n and ω is studied in Section 5 with the help of numerical experiments. For the moment we focus on g s 0 ,s 1 (iκ n (ω)) with respect to s 0 and s 1 . In Fig. 2 the level sets of g s 0 ,s 1 are given for the one pole basis s 1 = s 0 and a two pole basis s 1 = s 0 /10. The crucial factors are the values of the level sets at the poles iκ n . Therefore, some poles iκ n are plotted for a two dimensional waveguide with thickness 1 and Neumann boundary conditions. The frequency ω is chosen in the vicinity of a cut-off frequency, i.e. there exists one relatively small wavenumber κ 2 (ω) ≈ 0.15.
In the detail enlargements of Fig. 2 it can be seen that the pole iκ 2 (ω) near the origin for the two pole basis lies on the isoline with value 0.5 and for the one pole basis on the one with value ≈ 0.85. So the convergence rate of the second mode is much better for the two pole basis than for the one pole basis. Just as important is the fact, that the other rates are only mildly worsened by the two pole basis. This would change if a one pole basis scaled with respect to iκ 2 (ω), e.g. with parameters s 0 = s 1 = −0.2 + 0.2i, would be used.
The expected overall convergence rate is studied in Fig. 3 with respect to different parameters s 0 and s 1 . The curve and therewith the Hardy space H − ( ) is kept fixed. Using the arguments of Rem. 2.1 the coefficients c n of the evanescent modes vanish if a large part of the waveguide is added to the interior domain. In this case the solution in the waveguide is a linear combination of the propagating modes only. Therefore, only the rates for the propagating modes with positive wavenumbers κ n (ω) are considered in Fig. 3.
Parameters on the diagonal of Fig. 3 lead to s 0 = s 1 , i.e. a one pole basis. It can be seen, that the two pole basis with s 0 = s 1 gives much better convergence rates. Note, that this finding depends on the wavenumbers and will change for frequencies ω, which are not in the vicinity of a √ λ n . There is one drawback using the two pole basis. Numerical studies in [14, Fig. 6] indicate, that the ratio C 2 /C 1 of the Riesz constants in (3.9) increase linearly with |s 0 |/|s 1 |. Hence, the condition number of the discretization matrix might grow linearly with this factor.

Matrix allocation
The bilinear forms in (2.24) consist of interior and exterior parts. The interior parts are standard bilinear forms defined on the common space H 1 ( int )×H 1 ( int ). Since the matrix allocation for those bilinear forms is standard we omit it and state solely the matrices corresponding to the exterior bilinear forms (2.18).
Recall, that we introduced in Section 3.1 the tensor product space (3.16) which can be computed as [14,Sec. 4.4] by means of the residue theorem, since we are dealing with meromorphic quadratically decaying integrands. The off-diagonal elements of (3.16) alternate and the diagonal elements are constant with the exception of the first element. Note, that all matrices A ext , B ext , A int , B int are independent of the frequency ω and the wavenumbers κ n (ω). Hence, a discretization of (2.25) using this basis leads to a linear Matrix eigenvalue problem.

Hardy space mode matching for scattering problems
The two pole basis does not depend on wavenumbers or waveguide modes, which is a big advantage over standard mode matching techniques. To use the modes in a numerical method, they have to be computed first. For more complicated cross-sections ϒ such computations have to be done again numerically with sufficient accuracy. Typically this requires to solve a transversal eigenvalue problem for each frequency ω. So numerical studies over a large range of frequencies might lead to significant computational costs for the precomputation of the modes. Moreover, resonance problems lead to non-linear eigenvalue problems, if the wavenumbers or waveguide modes are used in the numerical method.
In this subsection we propose a manipulation of the two pole basis using one critical wavenumber κñ(ω) with g s 0 ,s 1 (iκñ(ω)) ≈ 1 explicitly in one extra basis function. This allows for an extremely efficient treatment of wavenumbers with small absolute value, e.g. ω 2 ≈ λñ for oneñ ∈ N 0 . But of course this critical wavenumber has to be precomputed and depends on the frequency ω. So the following variant should not be used for resonance problems.
We call it mode matching, since our Galerkin basis will include the transformed longitudinal part of the critical waveguide mode. Note, that we do not include the transversal part of any waveguide mode. This part is handled by the finite element discretization trΥ V int . So the transversal matrices from Section 3.
Due to the symmetry of the bilinear forms the matrices M long , S long are symmetric as well. Note, that this modified basis needs to be used carefully. The approximation error of the extra basis functions N in the spaceX (N−1) is qualitatively determined by g s 0 ,s 1 (p) N−1 . So if N is sufficiently large such that g s 0 ,s 1 (p) N−1 is small, the basis functions inX (N) become almost linearly dependent. However in this case, the mode matched basis function is not needed anyway.
Analogy to PML The Hardy space method offers the possibility of this implicit mode matching using only the information of one critical wavenumber. The discretization matrices remain local on ϒ and no information of the transversal part of the modes is necessary. For the PML one hybrid mode matched variant is known [3], which exploits the bi orthogonality of the waveguide modes. So mode matching combined with a PML is not as trivial as for the Hardy space infinite element method.

Analytical framework
The suitable analytical framework for (2.23) and (2.25) is the theory of Fredholm operators (see e.g. [21]) and holomorphic Fredholm operator functions (see e.g. [10,11]). We first show how to embed the Hardy space formulations into this framework and then give the convergence analysis.
Note, that in [18] a complete convergence analysis for scattering problems and a rudimentary convergence analysis for resonance problems is developed. We rely on these results and combine them with the abstract T-coercivity framework of [13] in order to give a full convergence analysis for resonance problems as well. At the end of this section we further show how to improve the standard estimate of the discretization error by the best approximation error.

Well-posedness
To relate sesquilinear forms and subsequently linear operators to the bilinear forms a and b defined in (2.24) we introduce a conjugation on V, i.e. an anti-linear, self inverse and unitary operator from V to V. Let z → z denote the standard conjugation for complex numbers z ∈ C. For u U ∈ V, in particular for u ∈ H 1 ( int ) and U ∈ So the standard complex conjugation is not the desired involution.
To this end in [17, Sec. 2.3] a suitable conjugation C S 1 in the Hardy space H + (S 1 ) of the complex unit disk was introduced by C S 1 F (z) := U(z) for all F ∈ H + (S 1 ). Using the Möbius transform (3.6) leads to the conjugation C := for the bilinear form q defined in (2.12). So C turns the bilinear form q into a sesquilinear form on H − ( ). Overall, we define C : It can easily be checked that C is indeed a conjugation on V. Hence let A, B : V ext → V ext be the bounded linear operators defined through The main difficulty to obtain Fredholm properties for A −ω 2 B is that it is not a weakly coercive operator, i.e. it is not a compact perturbation of a strictly coercive operator J . Strictly coercive is thereby meant in the usual V for some α > 0 and all u U ∈ V. The notion of T-coercivity [4,5] proves useful in our context. It relies on an automorphism T on V such that T * (A −ω 2 B) (T * denotes the adjoint of T ) is weakly coercive. Note, that weak T -coercivity implies Fredholmness of the operator A −ω 2 B.
The construction of T is carried out in [18,Sec. 4 and 6]. Remark that T is denoted S therein. The restriction of T to V ext , the exterior problem in the waveguide, reflects basically the behavior of the waveguide modes. Inside the waveguide for the finitely many guided modes (ω 2 > λ n ) T acts as multiplication with a complex scalar, while for the infinitely many evanescent modes (ω 2 < λ n ) T acts as the identity. This basic concept works for complex frequencies as well if iκ n (ω) / ∈ for all n ∈ N 0 . The set {ω ∈ C : ∃ñ ∈ N 0 with iκñ(ω) ∈ } coincides with the branch cuts of κ n (·) introduced in Sec. 2.1 whereby θ = arg(s 2 0 ). This set forms the essential spectrum σ ess of the operator function ω → A −ω 2 B, i.e. the set of values for ω for which the operator is not Fredholm.
It follows that for ω ∈ C \ σ ess the operator A −ω 2 B is bijective up to a discrete set σ dis which consists of eigenvalues of finite multiplicity and has no accumulation point in C \ σ ess . Hence, the resonance problem (2.25) is meaningful and if ω > 0 is no eigenvalue, (2.23) is well posed in the sense of Hadamard.
With the substitution λ = ω 2 the essential spectrum of λ → A −λ B consists of infinitely many half lines starting from λ n , n ∈ N 0 with angles s 2 0 . The dependency of σ ess on the artificial parameter s 0 is striking. Though, the essential spectrum is of little interest. More important is the dependency of the discrete spectrum on s 0 . We have seen that for positive ω all s 0 with negative real and positive imaginary part lead to the same radiation condition. Hence these resonances do not depend on s 0 . All other resonances depend in principle on the artificial parameter s 0 . However, if for a parameter s 0 and a complex resonance ω a different parameters 0 constitutes the same radiation condition at ω, then the resonance is kept invariant under the parameter change.
The latter can also be argued by holomorphicity properties. Since all resonances on the real axis are independent of s 0 , the identity theorem for holomorphic functions yields: If ω is a resonance for one s (1) 0 , then it is also a resonance for a second s (2) 0 if there exists a path in the complex plane connecting the real axis with ω such that Analogy to PML In [18,Sec. 5.2] this dependency is studied in more details for complex scaling. As indicated in Section 2.1, for a linear complex scaling or perfectly matched layer method the radiation condition is equivalent to the pole condition. So all comments on the dependency of the Hardy space method parameter s 0 remain valid for a linear complex scaling with parameter σ = 1/s 0 .

Convergence analysis
We finish this section with a short sketch of the convergence analysis based on [18] and [13]. The latter work, which is based on [19,20], derives a condition to ensure convergence of Galerkin methods to weakly T -coercive source and eigenvalue problems. The obtained condition is the existence of bounded linear operators T h,N on and to V h,N with lim h→0,N→∞ Since all other basis functions in the tensor product discretization of Section 3.
By definition (2.20) of the norm in V ext and continuous embedding · X 1 ≤ C · X 2 , it is straightforward to compute Hence, using (3.12) for C N−1 it follows lim h→0,N→∞ C − C N h,N = 0 which suffices to regain the former results.
So overall, the framework of [13] applies to the Hardy space infinite element method. In particular, using [13, Theorem 3.1(5)] for an isolated resonance ω of (2.25) with multiplicity one and corresponding normed eigenfunction (u, U ) ∈ V, there exist a constant C > 0 and a sequence of discrete resonances ω h,N such that for all sufficiently small h > 0 and all sufficiently large N ∈ N 0 , whereby (u * , U * ) ∈ V is a normed eigenfunction of the resonance ω to the adjoint resonance problem. The adjoint problem is of similar kind and hence the approximation properties for (u * , U * ) are the same and inf (v h

Improved convergence rates for scattering problems
For scattering problems, [13, Theorem 3.1(3)] or [18, Theorem 3.2] yield stability of the discrete operators and the usual error bound by the best approximation error. However, for practical applications only the physical variable u ∈ V int and hence the error measured in the interior domain u − u h,N H 1 ( int ) is relevant which is also the computed quantity for common numerical studies. Indeed the error bound of u − u h,N H 1 ( int ) by the best approximation error inf can be improved with respect to the longitudinal degrees of freedom by a duality argument, i.e. g s 0 ,s 1 (iκ n ) N/2 in (3.14) can be improved to g s 0 ,s 1 (iκ n ) N . To this end we introduce Dirichlet-to-Neumann operators and study the effect of a semi-discretization of the exterior problem on the error u − u h,N H 1 ( int ) .
Let ω > 0, ω 2 = λ n , n ∈ N 0 be no resonance of (2.25). We consider the standard Helmholtz Dirichlet-to-Neumann operator in cylindrical domains DtN : iκ n c n ϕ n . (4.5) That way (2.23) is equivalent to find u ∈ V int g such that To perform our analysis we first derive an identity for κ n as follows. Using the decomposition of Section 2.1 into the waveguide modes, the exterior part of the bilinear form in (2.23) diagonalizes with bilinear forms on X 2 ×X 2 and q as defined in (2.12). Note, that for solutions (u, ∞ In other words the Hardy space variational formulation of the waveguide problem can be interpreted as a realization of the standard Dirichlet-to-Neumann operator (4.5) on the interface to the interior domain.
We study the effect of a Galerkin discretization X (N) ⊂ X 2 in a similar way by introducing a perturbed Dirichlet-to-Neumann operator As in (4.10) the constants κ (N) n can be related to the bilinear forms s n . Using the results in [18] the bilinear forms s n in (4.8) correspond to uniformly in n ∈ N 0 strictly coercive sesquilinear forms with continuity constants of magnitude 1 + |κ n | 2 . Hence, for all N ∈ N there exist unique solutions which is the discretized version of (4.9). Note, that we only seek U (N) n and the scalar component c n is a priori fixed, i.e. we are solving a discrete Dirichlet problem in the waveguide. So the effective Dirichlet-to-Neumann operator using this discretization of X 2 is given by (4.11) with constants The effect of this perturbed Dirichlet-to-Neumann operator can be studied using the solution u N ∈ V int g to (4.6) with DtN N instead of DtN. Using the results from Section 4.2 this solution is unique for sufficiently large N ∈ N. More precisely, we have shown stability in the sense that there exist N 0 ∈ N and C 1 > 0 such that for all w ∈ V int 0 and all N > N 0 it holds Using w = u − u N and (4.6) we obtain |c n | 2 κ n −κ (N) n 2 (1+λ n ) −1/2 1 2 (4.14) with H 1 ( int ) → H 1/2 (Υ)-continuity constant C 2 . Using continuity of s n in (4.10) and (4.13) leads to with C > 0 independent of n ∈ N 0 . Hence, using Céa's Lemma we obtain essentially the same bound as in each summand in (3.5) and later on in (3.13). Therefore, for each waveguide mode with index n ∈ N 0 we gain the already proven convergence g s 0 ,s 1 (iκ n ) N/2 . The main advantage of this Dirichlet-to-Neumann approach is, that we can improve this order by a duality argument. Note, that by (3.16) the bilinear forms s n are tridiagonal with respect to ( long j ) j ∈N 0 and thus κ n resp. κ . Hence, (4.15) can be replaced by with a different constant C > 0. The clue is now that the naive estimate can be improved by an Aubin-Nitsche trick for δ : for all W N n ∈X (N) . In order to estimate the last term, we need to calculate W n . If we replace (4.17) by thenW n = δ(iκ n − s 0 )F iκ n due to (4.10). A matrix representation of both problems can be formulated using the matrices M long inf and S long inf of (3.16), whereas (4.17) leads to the N×N sub matrix of the N 0 ×N 0 matrix representation of (4.19). Hence, it holds W n = α n s 0 ,s 1 0 + β n F iκ n with constants α n , β n ∈ C. So by the results of Section 3.2 and in particular with (3.12), from both two terms in (4.18) we obtain a convergence order g s 0 ,s 1 (iκ n ) N/2 leading to the improved error bound (4.20) Analogy to PML The technique to obtain a better convergence rate through the interpretation of the semi-approximation of the waveguide as a perturbed Dirichletto-Neumann operator works also for PML. In fact, the rate exp( (σ iκ n )L) obtained in [18] by the simple estimate differs from the rate exp(2 (σ iκ n )L) obtained in [1] by a perturbed Dirichlet-to-Neumann argument (L denotes the layer width).
Note, that for waveguides it is more favorable to chose a Neumann condition at the end of the perfectly matched layer, because it leads to the better estimate min{1, |κ n |} exp(2 (σ iκ n )L).

Numerical results
We study in this section the performance of the Hardy space infinite element method in the vicinity of cut-off frequencies. First we solve a scattering problem with analytically known exact solution. In the second part we give self convergence tests for a resonance problem. Both examples use a two dimensional cylindrical waveguide of thickness one with homogeneous Neumann boundary conditions at the waveguide boundary.
Cylindrical waveguides in R 3 would only complicate the finite element discretization of the interior problem and of the interface ϒ. The performance of the Hardy space method in longitudinal direction is independent of the shape and dimension of the interface. Therefore we refrained from numerical studies in three dimensions.
All computations were done within the high order finite element package Netgen/NgSolve [28,29] and the module [23].

Scattering problem
Let Hence at the Dirichlet boundary {−0.3} × ϒ we impose the datum 6 n=0 cos(nπη). Due to the small interior domain the constants c n := exp (0.3iκ n (ω)) for n = 0, . . . , 6 have absolute values larger than | exp (0.3iκ 6 (0)) | ≈ 0.003 for all ω > 0. Hence, since the first seven waveguide modes contribute to the error in the waveguide, one might set n 0 = 7 in (4.20). Increasing the interior domain would reduce the number of evanescent modes with remarkable contribution to the overall error. Note, that in Section 5.2 a resonance problem where all modes are present is solved, see Fig. 10.
The interior space V int is discretized by a finite element space with a very coarse triangular mesh (see Fig. 7c for the discretization of a similar domain) and uniform polynomial order 10. This leads to finite element approximation errors of magnitude 10 −10 for u. For the Hardy space method we use a fixed parameter s 0 = −2 + 2i and parameters s 1 = cs 0 with varying c ∈ {0.001, 0.01, 0.1, 1}. Additionally, we use the mode matching of Section 3.4 with s 1 = s 0 . We measure the error between the exact solution and the discrete solution Figure 4 shows the error u − u h,N H 1 ( int ) for different frequencies ω ∈ [π, 3π ] for fixed number of degrees of freedom N = 20 in longitudinal direction. We observe that in the frequency region away from the cut-off frequencies {0, π, 2π, 3π, . . . } the error increases with decreasing c, but stays bounded away from 1. This was to be expected since for these frequencies ω all n = 0, . . . , 6 yield |κ n (ω)| 0 and the decisive factor g s 0 ,cs 0 (iκ n (ω)) (defined in (3.10)) increases with decreasing c. In the limit c = 0 it holds g s 0 ,0 (iκ n (ω)) = |iκ n (ω) − s 0 |/|iκ n (ω) + s 0 | < 1.
The error in the neighborhood of the cut-off frequencies, which we aimed to improve, is studied in Fig. 5 in logarithms scale with respect to the distance δ ω of ω to 2π = √ λ 2 . Mode matching is by far the best choice, but of course knowledge of the value κ 2 (ω) is needed. If the wavenumber is not available, the two-pole basis still leads to small errors using a proper scaling of s 1 with respect to the approximate scale of |κ 2 (ω)|. The value c = 0.01 seems to be almost optimal for δ ω ∈ [10 −6 , 10 −2 ]. For smaller distances, smaller values of c are needed to minimize the approximation error arising from g s 0 ,s 1 (iκ 2 (ω)) N . For larger distances δ ω > 10 −2 the standard one pole basis seems to be the best choice. Note, that in Fig. 2 the distance δ ω of ω to √ λ 2 was approximately 10 −3 . Figure 6a shows the convergence with respect to the number of longitudinal unknowns N for a frequency ω in the vicinity of √ λ 2 . In order to study the different terms in (4.20) some of the results of Fig. 6a are compared to |c n |g s 0 ,s 1 (iκ n ) N for n = 0, 2, 6. The values of g s 0 ,s 1 (iκ n ) are given in Table 1 and the constants c n are stated after (5.1). In particular the results for the mode matched variant are compared in Fig. 6c, the one pole basis (c = 1) in Fig. 6d, and the two pole basis with c = 0.05 in Fig. 6b.
Let us start with Fig. 6c and the results of the mode matched variant of the one pole Hardy space elements. Since the Laplace transformed longitudinal part of the second mode is handled explicitly in this variant of the method, the corresponding term with n = 2 does not appear in (4.20). From Fig. 6c we deduce that it has indeed no effect on the computational error. Since c 7 = 0, we expect a convergence with factor g s 0 ,s 1 (iκ 6 ), which is the maximal contributing entry of the first row of Table 1. This expectation is met by the results in Fig. 6c for N ≥ 10. Nevertheless, |c 6 | ≈ 0.003 |c 0 | = |c 1 | = 1. Therefore, for small values of N we see the faster convergence of the zeroth mode.
For the one pole basis Fig. 6d shows, that the error in the second mode dominates the errors of the other modes as expected by the values in Table 1. Only for small N it seems that the zeroth mode dominates the error. But clearly the performance of the one pole method suffers from the poor convergence of the second mode. Figure 6b supports the use of a properly scaled two pole method. Asymptotically dominant is the error in the 6th waveguide mode. For small N again the error in the zeroth mode dominates. Due to the scaling the convergence rate of the error in the second mode is so much improved, that it is negligible for the overall performance.
We finish this subsection with some final remarks on the index n 0 in (4.20). From Fig. 6c we may deduce, that the finite element error of the interior domain is approximately 10 −10 . If this is the given tolerance for the Hardy space method, n 0 should The mode matched variant with s 1 = s 0 is is labeled "special". Additionally, in (b)-(d) the computational error is compared to |c n |g s 0 ,s 1 (iκ n ) N for n ∈ {0, 2, 6} with g s 0 ,s 1 (iκ n ) given in Table 1 and constants c n stated after (5.1) be chosen such that C exp(−c λ n 0 ) in (4.20) is below this tolerance. Since in this case, the constants c n for all n ≥ n 0 are so small, that these modes will lead to errors below the given tolerance for any N ∈ N. Table 1 Expected convergence rates g s 0 ,s 1 (iκ n ) N compared to the numerical rates identified by the results of Fig. 6a. The convergence rate of the mode matching variant is approx. 0.6413 N g s 0 ,s 1 (iκ 0 ) g s 0 ,s 1 (iκ 1 ) g s 0 ,s 1 (iκ 2 ) g s 0 ,s 1 (iκ 3 ) g s 0 ,s 1 (iκ 4 ) g s 0 ,s 1 (iκ 5 ) g s 0 ,s 1 (iκ 6 ) g s 0 ,s 1 (iκ 7

Resonance problem
For resonance problems the frequency is part of the sought unknowns and therefore mode matching as in Section 3.4 would lead to a highly non-linear eigenvalue problem. Therefore we use at this point only the two pole basis of Section 3.2. The waveguide ext is exactly the same as in the preceding subsection, but we change the interior domain int ⊂ [−0.3, 0] × ϒ (see Fig. 7). We chose a parameter . So we obtain parameter dependent resonances ω α of (2.25) using homogeneous Neumann boundary conditions on the whole boundary ∂ . Figure 7 shows for three different values of α the absolute value of a resonance function corresponding to the resonance closest to 2π = √ λ 2 . They belong to the blue thick curve of resonances ω α in Fig. 8 which converges to 2π for α → 0. The thick red and green curves of resonances converge to π and 3π respectively. The computed eigenvalues in light colors build a discretization of an essential spectrum and correspond to wavenumbers with iκ n (ω) ∈ (see [18] for the details). So these numerical eigenvalues are discretization effects and must not be mistaken to be approximations of resonances of (2.25).
Since the blue resonances ω α of Fig. 8 converge to √ λ 2 for α → 0, the absolute value of iκ 2 (ω α ) becomes small (see the caption of Fig. 7 for squared magnitudes). In this case we are in the vicinity of a cut-off frequency and therefore expect the two pole basis to be superior to the one pole basis. Figure 9 supports this expectation. With decreasing distance of ω to √ λ 2 the two pole basis with smaller values for c becomes more and more efficient. The optimal values of c corresponds almost exactly to the values indicated by Fig. 5 for the scattering problem. Note, that Fig. 10 indicates that the main contribution to the error arises from the second waveguide mode. This effect becomes more and more dominant with decreasing distance to √ λ 2 . Nevertheless, all waveguide modes are present. Thus for the one pole basis adjusting s 0 only to the critical wavenumber κ 2 (ω) would lead to large errors for the other waveguide modes. So the two scale Hardy space method is indeed necessary for an efficient method for wavenumbers on different scales.

Conclusion
In [14] a two pole Hardy space infinite element method was developed to treat backward wave phenomena which occur e.g. in elastic waveguides [15]. Different to that goal, we proposed in the present work to apply a certain variant of the two pole elements also to systems which allow for wavenumbers at different scales. As model problem we considered time-harmonic acoustic cylindrical waveguides in the neighborhood of cut-off frequencies. If the method parameters s 0 and s 1 of the infinite elements are adjusted to the two different scales of the (low indexed) wavenumbers, common performance is recovered. In particular, up to arbitrary small thresholds exponential convergence with respect to the number of longitudinal unknowns is obtained. The method is independent of the frequency and therefore well suited for the numerical solution of resonance problems. Additionally, for scattering problems where the vanishing wavenumber is explicitly known, we generalized the mode matched variant proposed in [18] to the two pole elements. The transversal part of the mode itself is not needed, only the wavenumber enters the method and the discretization matrix stays local at the cross-section of the waveguide.
We reviewed the analysis of [18] within the T-coercivity framework of [13] which facilitates the proofs and extends the results of [18]. Through this approach the generalization of the convergence results to the two pole and mode matched elements becomes rather straightforward. We further derived an improved error estimate (compared to the best approximation error) w.r.t. the longitudinal degrees of freedom. Hence we support the claimed convergence and performance of our proposed method with a sound analytical framework.
In Section 5 we gave different numerical experiments which validate our theoretical results and guide potential users of the method with the optimal choice of parameters. An obvious extension of the method would be to use the mode matched variant with more than one extra waveguide mode or to use a kind of multi scale Hardy space infinite element method with more than two poles. In principle this seems to be possible, but all discretizations of the Hardy space have to be chosen not only to optimize the best approximation error, but also to ensure moderate condition numbers of the discretization matrices. So a detailed analysis as for the two pole elements is necessary in order to construct well suited basis functions.