Numerical methods for scattering problems in periodic waveguides

In this paper, we propose new numerical methods for scattering problems in periodic waveguides. Based on [20], the “physically meaningful” solution, which is obtained via the Limiting Absorption Principle (LAP) and is called an LAP solution, is written as an integral of quasi-periodic solutions on a contour. The definition of the contour depends both on the wavenumber and the periodic structure. The contour integral is then written as the combination of finite propagation modes and a contour integral on a small circle. Numerical methods are developed and based on the two representations. Compared with other numerical methods, we do not need the LAP process during numerical approximations, thus a standard error estimation is easily carried out. Based on this method, we also develop a numerical solver for halfguide problems. The method is based on the result that any LAP solution of a halfguide problem can be extended to the LAP solution of a fullguide problem. At the end of this paper, we also give some numerical results to show the efficiency of our numerical methods.


Introduction
Numerical simulation for scattering problems in periodic waveguides is an interesting topic in both mathematics and related technologies, due to its wide applications in optics, nanotechnology, etc. As the well-posedness of the scattering problems is not

R. Zhang
The scattering problem is modeled by the following equations: where f is a function in L 2 (Ω) with compact support, and k > 0 is the real wavenumber.

Remark 1
We can also consider different boundary conditions on ∂ Ω, e.g., the Dirichlet boundary condition or the Robin boundary condition, in the similar way. In this paper, we only want to take the Neumann boundary condition as an example for our methods.
To find a "physically meaningful" solution of (1)-(2), we introduce the well-known Limiting Absorption Principle (LAP). That is, given any > 0, consider the following damped Helmholtz equation: From the Lax-Milgram theorem, the problem is uniquely solvable in H 1 (Ω). Moreover, the solution u decays exponentially as x 1 → ±∞ (see [8]). The LAP assumes that u converges in H 1 loc (Ω) when → 0 + , and the limit, denoted by u, is the "physically meaningful" LAP solution. In the following, we introduce a new formulation for LAP solutions. Based on the new formulation, we introduce new numerical methods to compute the LAP solutions efficiently.
For simplicity, we introduce the following operator: Let the spectrum of A be denoted by σ (A), then the problem (1)-(2) is uniquely solvable in H 1 (Ω) if and only if k 2 / ∈ σ (A). To study the spectrum property of A which plays an important role in the well-posedness of the problem (1)-(2), we introduce the Floquet-Bloch theory. For details we refer to [2,8] and for more general results we refer to [25].

Floquet-Bloch theory
In this section, we introduce the Floquet-Bloch theory to study the spectrum σ (A). For simplicity, we introduce the following domains (see Fig. 1): Then Ω = j∈Z Ω j with the left and right boundaries j and j+1 . Let We also introduce the space of quasi-periodic functions. A function φ ∈ H 1 loc (Ω) is called z-quasi-periodic, if it satisfies for some fixed complex number z ∈ C. We define the subspace of H 1 (Ω 0 ) by: Then the functions in H 1 z (Ω 0 ) can be extended to z-quasi-periodic functions. Especially, when z = 1, all functions in H 1 1 (Ω 0 ) can be extended into periodic functions in H 1 loc (Ω). We also denote H 1 1 (Ω 0 ) by H 1 per (Ω 0 ). From the Floquet-Bloch theory, the spectrum of A is closely related to Bloch wave solutions. A Bloch wave solution is a non-trivial solution of (1)- (2) in H 1 z (Ω 0 ) with f = 0 for some z ∈ C. If a Bloch wave solution exists in H 1 z (Ω 0 ), z is called a Floquet-multiplier. Define the operator: where D(A, Ω 0 ) is defined in the same way as D(A, Ω), and Ω is replaced by the periodic cell Ω 0 . Note that, A z is self-adjoint with respect to the L 2 -space equipped with the weighted inner product (φ, ψ) L 2 ,q = Ω 0 qφψ dx . Let σ (A z ) be the spectrum of A z , then k 2 ∈ σ (A z ) if and only if z is a Floquet multiplier. Let F(k 2 ) be the collection of all Floquet multipliers with wavenumber k and UF(k 2 ) := F(k 2 ) ∩ S 1 (S 1 is the unit circle in C) be the set of all unit Floquet multipliers. In this paper, when the wavenumber is fixed, we write F instead of F(k 2 ) for simplicity. We list the properties of the Floquet multipliers [20], for more details we refer to [2,6,8,21,25]: -UF is finite.
z ∈ F if and only if z −1 ∈ F, thus z ∈ UF if and only if z = z −1 ∈ UF.
-F is a discrete set and the only finite accumulation point of F can be 0.
-F(k 2 ) depends continuously on k 2 . A classical result from the Floquet-Bloch theory also shows that (see [25]):
For any fixed n ∈ N, the graph (α, μ n (α)) : α ∈ (−π, π] is called a dispersion curve, and all dispersion curves compose a dispersion diagram. Following [5,6], we first show the dispersion diagrams for two different examples: 1. Example 1. q = 1 is a constant function, and its dispersion diagram is shown in Fig. 2 (left). The dispersion curve is given analytically: 2. Example 2. q = 9 in a disk B (0, 0.5), 0.3 and q = 1 outside the disk, and its dispersion diagram is shown in Fig. 2 (right)..
In the right picture of Fig. 2, there are "stop bands" (in red). When k 2 lies in the stop bands, the horizontal line with height k 2 has no intersection with any dispersion curves. This implies there is no propagating mode and the scattering problem (1)-(2) has a unique solution in H 1 (Ω). The rest of bands are called "pass bands", such as the whole domain in the left picture of Fig. 2 and the white region in the right picture of Fig. 2. When k 2 lies in a pass band, from (7), there is at least one α ∈ (−π, π] such that there is a non-trivial α-quasi-periodic function ψ satisfying A α ψ = k 2 ψ. Thus it is a Bloch wave solution and is called a propagating Floquet mode. The case when k 2 lies in a pass band is particularly interesting and challenging. Thus we discuss more details about this case. When k 2 lies in a pass band, there is at least one α ∈ (−π, π] such that k 2 ∈ σ (A α ). Then the set Thus UF can be written as: The points in P are divided into the following three classes: -When μ n (α) > 0, ψ n (·, α) propagates from the left to the right; -when μ n (α) < 0, ψ n (·, α) propagates from the right to the left; -when μ n (α) = 0, we can not decide the direction that ψ n (·, α) propagates.

Remark 2
It is possible that there are two (or more) different dispersion curves passing through the point (α, k 2 ). Suppose the elements in P have Q different values In this case, α j is treated as L j different elements in P, i.e., α j,1 = α j,2 = · · · = α j,L j = α j .
As UF is symmetric, P ± is also symmetric, i.e., α ∈ P + if and only if −α ∈ P − . For details see Theorem 4, [2].  As the limiting absorption principle fails when the set P 0 is not empty, we make the following assumption.
Assumption 1 Assume that in this paper, P 0 = ∅.
The assumption is reasonable as the set k > 0 : P 0 = ∅ is "sufficiently small", i.e., the set is countable with at most one accumulation point at ∞ (see Theorem 5,[2]).
In our later proof of the new integral representation of LAP solutions, we also have to avoid the cases when P + ∩ P − = ∅. Luckily, with the similar method used in the proof of Theorem 5 in [2], we can also prove that this set is discrete in the following lemma. For the proof we refer to Appendix.
We define three subsets of UF from the definition of P ± and P 0 by: Then UF = S 0 We also divide the set F\UF into the following two subsets: The Bloch wave solution corresponding to z ∈ RS is evanescent on the right, while the one corresponding to z ∈ L S is evanescent on the left. Moreover, let

Remark 3
We conclude the properties of the sets RS, L S and S 0 ± from the properties of F as follows: -From the symmetry of F, the sets S 0 + and S 0 − , RS and L S are symmetric, i.e., This also implies that z ∈ S + ⇐⇒ z −1 ∈ S − .
We omit the proof here. For details we refer to [20].

Quasi-periodic problems
From the last section, quasi-periodic problems are very important in the investigation of scattering problems in periodic domains. In this section, we consider the z-dependent cell problem: where z = 0 is a complex number and f z ∈ L 2 (Ω 0 ) depends analytically on z.
We are interested in the well-posedness of the problem (11)- (13), and also the dependence of the solution u z on the quasi-periodicity z. To this end, it is more convenient to study the problems in a fixed domain.
. Note that as z −x 1 is a multi-valued function, we require that z lies in the branch cutting off along the negative real axis (denoted by C × := {z ∈ C\{0} : −π < arg(z) ≤ π }, where arg(z) is the argument of the complex number z). From direct calculation, v z is the solution of the following problem: For each fixed z, the problem is equivalent to the following equation. For definitions of the corresponding operators/elements we refer to Sect. 4, [20].
This implies that when As K z is compact and depends analytically on z ∈ C × , I − K z is an analytic family of Fredholm operators. From the Analytic Fredholm Theory (see Theorem VI.14, [26]), [20], the set of poles of (I − K z ) −1 is exactly F. Thus v z , or equivalently u z = z x 1 v(x), depends analytically on z ∈ C\ F {0} and meromorphically on z ∈ C\{0}.

The Floquet-Bloch transform
The Floquet-Bloch transform is a very important tool in the analysis of scattering problems in PDEs in periodic structures, see [2][3][4].
For a function φ ∈ C ∞ 0 (Ω), define the Floquet-Bloch transform of φ by The transform is well-defined for any smooth function with compact support, and it can be extended to more general cases. Define the Region of Convergence (ROC) as the domain in C such that the series (15) converges. Note that the DOC may be empty for given function φ. When φ decays exponentially at the rate γ , i.e., there is a γ > 0 and C > 0 such that φ satisfies the Floquet-Bloch transform of φ is still well-defined, and the ROC is the annulus Moreover, from the perturbation theory, the function (Fφ)(z, ·) depends analytically on z ∈ T γ . It is also easy to check that the transformed function (Fφ)(z, ·) is quasiperiodic (i.e., it satisfies (5)). We conclude some mapping properties of the operator F in the following proposition.

Proposition 1
The operator F has the following properties when z lies on the unit circle S 1 (see [11,21]): -F is an isomorphism between H s (Ω) and L 2 (S 1 ; H s z (Ω 0 )) ( s ∈ R), where -Fφ depends analytically on z ∈ T γ , if and only if φ decays exponentially with the rate γ . -Given ψ(z, x) := (Fφ)(z, x) for some φ ∈ H s (Ω) and satisfies (16) for some γ > 0, the inverse operator F is given by:

Remark 4
From now on, we assume that supp( f ) ⊂ Ω 0 . The results are easily extended to cases when supp( f ) lies in larger bounded domains.
We define the Floquet-Bloch transform w(z, x) := (Fu)(z, x), then the transformed field w(z, ·) is well-defined and depends analytically on z ∈ T γ as a H 1 (Ω 0 )-function valued function. It is also easy to check that for any z ∈ T γ , w(z, ·) ∈ H 1 z (Ω 0 ) satisfies (11)- (13). Note that the source term in (11) is f , as From the inverse Floquet-Bloch transform and Cauchy integral theorem, the solution of the original problem can be represented as: where C is a rectifiable curve in T γ encircling 0. From the exponential decay of u, (Fu)(z, ·) exists and depends analytically on z ∈ T γ . On the other hand, from Sect. 3.2, when z ∈ C\F, the problem (11)-(13) is uniquely solvable in H 1 z (Ω 0 ), and w(z, ·) depends analytically on z ∈ C\(F {0}) and meromorphically on z ∈ C\{0}. Thus (Fu)(z, ·) is extended meromorphically in C\{0}. From the analytic continuation, we obtain the following result (see [20]).
The integral representation of u is obtained from Cauchy integral theorem.
is the solution of (11)-(13) for z ∈ C\F. Then the solution of (1)-(2) is written as where C ⊂ C is a counter-clockwise closed rectifiable path encircling all the points in RS(= S + ) and does not encircle any point in L S(= S − ).

The Limiting absorption principle (LAP)
In this section, we consider the case when k 2 ∈ σ (A) with the help of the limiting absorption principle. In the first two sections, we always assume that Assumption 2 holds. First we consider the damped Helmholtz equation (3)- (4). The corresponding z-quasi-periodic problem is formulated as: Similar to the definition of K z , we denote the operator with k 2 + i by K z . From Theorem 4 in [22], the poles of the operator (I − K z ) −1 depend continuously on . First, we study the asymptotic behaviour of distribution of the poles when > 0 is sufficiently small.

Distribution of poles of the damped Helmholtz equations
From the Floquet-Bloch theory (see [2,25]), k 2 ∈ σ (A) implies that UF = ∅. For any z ∈ UF, k 2 is an eigenvalue of A z . As both S 0 + and S 0 − are finite sets, and S 0 + and S 0 − , RS and L S are symmetric in the sense of (10), they are written as where z + j ( j = 1, . . . , Q) are Q different values on the unit circle, and so are z − for any integer j ∈ N. From Assumption 1 and 2, UF = S 0 Analogous to the case = 0, we define the following sets depending on : From the continuous dependence of > 0, we have the following properties of Z ± j ( ) when j = 1, 2, . . . , Q. For the proof we refer to Appendix in [8].

Lemma 3 For any j
Thus we conclude the behaviour of Z ± j ( ) for sufficiently small : -for any z + j ∈ S 0 + , the points Z + j, ( ) ∈ S 0 + ( ) converge to z + j from the interior of the unit circle; -for any z − j ∈ S 0 − , the points Z − j, ( ) ∈ S 0 − ( ) converge to z − j from the exterior of the unit circle.
To make it clear, we present a visualization of examples of the curves in Fig. 4. The red rectangles are points in S 0 − and the blue diamonds are points in S 0 + . The asymptotic behavior of Z ± j, ( ) as → 0 can be seen from the picture. For the points in RS( ) or L S( ), we estimate their distributions for sufficiently small > 0 in the following lemma (see Lemma 17,[20]).

Lemma 4 Suppose for some
From Lemmas 3 and 4, when > 0 is sufficiently small, the sets S 0 ± ( ) and RS( ), L S( ) have following properties: R. Zhang

Integral representation of the LAP solution
Now we are prepared to formulate the LAP solution of (1)-(2) when k 2 ∈ σ (A). From Theorem 2, the solution u (∀ > 0) of the damped problem (3)-(4) is represented by the curve integral: But when → 0 + , as the poles in S 0 ± ( ) approach S 1 , and the integral becomes singular. From Theorem 2, we are aimed at finding out a proper curve C to replace the unit circle such that the function w (z, x) is well-posed and converges uniformly with respect to z when → 0. From the properties of the sets S 0 ± ( ), RS( ) and L S( ), we define a convenient curve C 0 as follows.
Definition 1 Let the piecewise analytic curve C 0 be defined as the boundary of the following domain : where δ > 0 is sufficiently small such that the following conditions are satisfied: i.e., the balls do not contain any point in RS L S.
From Assumption 2 and Lemma 2, we can always find a proper parameter δ such that both conditions are satisfied. We refer to Fig. 4 for a visualization of C 0 for the example when n = 1 and k 2 = 3π 2 .
Thus from Lemmas 3 and 4, there is a constant C = C(δ) > 0 such that holds uniformly for any fixed sufficiently small > 0, where d(X , Y ) is the Hausdorff distance between two subsets X , Y ⊂ C. From the choice of the curve C 0 , the interior of the symmetric difference of B δ and B(0, 1) is As for sufficiently small > 0, (I − K z ) −1 depends analytically on z in the above domain, from Cauchy integral theorem, u has the equivalent formulation From the asymptotic behaviours of the poles Z ± j ( ) explained in Lemmas 3-4, let tends to 0 + , we get the exact formulation for the LAP solution: We also refer to Sect. 6.2, [20] for a similar approach. Moreover, for any n ∈ Z, there is a constant C = C(n) > 0 such that Moreover, from regularity of elliptic equations, the LAP solution u ∈ H 2 loc (Ω). We can also replace C 0 by any closed curve that lies in the neighbourhood of C 0 and enclose zero and all poles in S + , but no poles in S − (see Fig. 4 (right)). The left curve is C 0 , and the right is a choice of C. Thus Now we have already formulated the LAP solution directly from cell problems (19)- (21), which also provides a nice and clear formulation for the numerical scheme. However, we still need to know the dispersion diagram to find the correct curve C 0 (or C).

Formulation for more general case
The formulation (25) is based on Assumption 2, which is not necessary for the LAP process. In this section, we recall another formulation introduced in [20] which does not require it.
For fixed k > 0, we recall the sets P, P ± and P 0 (also the corresponding S, S ± and S 0 ). With Assumption 1, P 0 = ∅ thus P = P + ∪ P − (and also S = S + ∪ S − ). However, since P + ∩ P − may not be empty, we can not use the old notations any more. Instead, we simply let Let ψ j, (β, ·) be the eigenfunctions corresponding to the eigenvalues μ j, (β). We choose a τ > 0 such that the circle |z| = exp(−τ ) enclose all poles in RS, then from [20] with a similar complex contour integral with the LAP process, for all n ∈ Z:

Numerical scheme
In this section, we consider two numerical methods to approximate LAP solutions of (1)-(2). The first method (CCI-method) is based on the complex curve integral (25) and the second method (PM-method) is based on the propagating modes (28). Both cases involve numerical integral of solutions w(z, x) C 0 of the quasi-periodic problem (19)-(21) (or equivalently, the periodic problem (14)). As the quasi-periodic problems indexed by zs are well-posed, they can be solved by classic numerical methods.
We apply the finite element method to solve the periodic problem (14). Suppose that Ω 0 is covered by a family of regular and quasi-uniform meshes (see [32,33]) M h with the largest mesh size h 0 > 0. To construct periodic nodal functions, we suppose that the nodal points on the left and right boundaries have the same height. Omitting the nodal points on the right boundary, let φ ( ) M : 1 ≤ ≤ M be the piecewise linear nodal functions that equals to one at the -th nodal and zero at other nodal points, then it is easily extended into a globally continuous and periodic function in H 1 loc (Ω). M is the number of nodal points which do not lie on the right boundary. Define the discretization subspace by: Thus we are looking for the finite element solution to (14) with the expansion By Theorem 14 in [12] the finite element approximation is estimated as follows.
where C > 0 is a constant independent of z ∈ C 0 . The estimations are also true for the function w(z, In both methods, w h (z, ·) is the numerical approximation of the quasi-periodic solution w(z, ·).

CCI-method
In this subsection, we approximate the LAP solution numerically with the help of the curve integral (25). w(z, ·) is supposed to be known as it is easily computed by standard numerical methods. We can always choose different n's such that computations are carried out in different Ω n 's, but in this section n is taken to be 0 as an example. As the curve integral on S 1 is standard, we only consider the case k 2 ∈ σ (A), so the LAP solution u is written in the form of (25). Recall that the number of different elements in S 0 + is Q, and so is it in S 0 − . Thus UF = S 0 + S 0 − has 2Q elements when Assumption 2 holds. Especially, from Remark 3 and Assumption 2, neither 1 nor −1 lies in UF. Then the curve C 0 is parameterized as follows.

R. Zhang
By assumption, for any j = 1, . . . , 2Q, there is a sufficiently small δ > 0 such that the ball B exp(iε j ), δ does not contain any other poles except exp(iε j ) itself. Let Thus the curve segment of S 1 ∩ C 0 is parameterized as , δ lies either in the exterior or in the interior of the unit disk, depending on j . As the curve segment is part of the circle with center exp(i j ), by choosing proper angles θ Now the whole curve C 0 is parameterized piecewisely. Thus the representation of u in (25) is written as Each integrand depends smoothly on t. Thus we only need to consider the numerical integration where a < b and g depends smoothly on t ∈ [a, b] and g(t, ·) ∈ H 2 (Ω 0 ) for any fixed t.
For an efficient numerical integral, we adopt the method introduced in [14], which comes originally from Sect. 9.6, [29]. Let s(τ ) : (a, b) → [−π, π] be a function in C ∞ (a, b) and strictly monotonically increasing function and satisfies s(−π) = a; s(π ) = b; s ( ) (−π) = s ( ) (π ) = 0, = 1, 2, . . . , N 0 for some positive integer N 0 . Let t = s(τ ), the integral I(g) becomes Thus the new integrand h depends smoothly on τ and is extended to a periodic function with respect to τ . We approximate the integral I(g) by trapezoidal rule. Let [−π, π] be divided uniformly into N subintervals, and the grid points be The integral is approximated by Then we recall the result of Theorem 9.33 in [29], and obtain the error estimation of the integral via (31): We apply the quadrature rule (31) to approximate u(x): where s j be the smooth function from [ + j , − j+1 ] to [−π, π] and s j+2Q are the smooth function from [θ − j , θ + j ] to [−π, π]. With (32), we conclude that Finally, we replace w(z, x) by the finite element solution w h (z, x). Then the LAP solution is approximated by ) .
(35) Then we conclude the error between u N ,h and u in the domain Ω 0 .

where C is a constant that depends on Q and N 0 , but does not depend on N and h.
Proof From the representations of u N ,h and u, and also the results from (34), we have the following error estimation: The estimation of the H 1 −norm is also obtained in the same way: The proof is finished.

PM-method
In this subsection, we introduce another numerical method based on (28). The integration part in (28) can be approximated in the same way as the CCI method, thus we only need to deal with the term expanded by eigenfunctions. The first step is to find out explicit values of β j and corresponding eigenfunctions ψ j, (β j , ·). From the variational formulation (14) and by replacing log(z) by iα in (14), the problem is to find α ∈ (−π, π] such that there is a non- for any φ ∈ H 1 per (Ω 0 ). From Theorem 2.4, [3], it is written as a generalized eigenvalue problem: By solving this problem, we obtain the values of β j and the corresponding eigenfunctions ψ j, (β j , ·) at the same time.
Then the values of f , ψ j, (β j , ·) are easily calculated by numerical integral on triangular meshes. The second step is to evaluate (μ j, ) (β j ). A natural idea is to compute the derivative directly by a standard formula, such as the Newton's difference quotient or the symmetric difference quotient. However, since in these formulas, the denominators are always small and the numerators are not exact, the errors are expected to be larger. Due to this disadvantage, we use an alternative formula to compute (μ j, ) (β j ).
First we need the definition of the energy flux: The energy flux is closely related to the directions of propagating modes. When φ is propagating to the right, E(φ; φ) > 0; while when it is propagating to the left, E(φ; φ) < 0. Then we recall the following equation from the proof of Theorem 3 in [2]: 2k .

Fig. 5 Periodic waveguide
From the above formula the sign of μ j, (β j ) is exactly the same as that of the energy flux E(ψ j, (β j , ·), ψ j, (β j , ·)). Thus (28) is rewritten as: where j, is an indicator function with the value 1 if E(ψ j, , ψ j, ) > 0, and with the value 0 if E(ψ j, , ψ j, ) < 0. We will show numerical results for the PM method also for wavenumbers that do not satisfy Assumption 2.

Half-guide problems
The numerical method introduced in this paper is also extended to half-guide problems. Let the half waveguide be defined by Ω + := ∞ n=1 Ω n (see Fig. 5). Recall that Then we are looking for the LAP solution u + ∈ H 1 loc (Ω + ) such that it satisfies From [1], for any φ ∈ H 1/2 ( 1 ) the LAP solution exists for all positive wavenumbers satisfying Assumption 1 except for a discrete set. In this section, we introduce numerical methods to approximate those LAP solutions with Assumptions 1 and 2.
To solve this problem with the method based on the approximation (39), a two-step method is proposed. The first step is to determine a source term f which approximates the boundary condition on 1 . The second step uses the source term to compute u + , with the same method as the full-guide problem. As the second step is exactly the same as the full-guide problem, we only discuss the first step in this section.
Let A be the operator from Since implies that A is a compact operator with dense range. We are looking for an f such that As the equation is severely ill-posed, we apply the well-known Tikhonov regularization method to find the "best solution" of this equation. For details we refer to Chapter 2.3 in [34].

Numerical method
Now we conclude the numerical scheme for half-guide scattering problems as follows. We look for a best solution c such that c − φ has the minimum value.  (19)- (21) with the source term f M 0 γ . The LAP solution u is then obtained by either the CCI method or the PM method.

Solve the equation
From the convergence of the Tikhonov method and numerical solution of (39), this numerical scheme converges as M 0 → ∞, γ → 0, h → 0 and N → ∞.

Formulation for special wavenumbers
In previous sections, all discussions for the CCI method are based on Assumption 2, i.e., P + ∩ P − = ∅. As it is not a necessary condition for the limiting absorption principle, we extend the CCI method to the case P + ∩ P − = ∅ in this section. Assume that for the wavenumber k 0 > 0, Assumption 1 holds but Assumption 2 does not hold.
Following the notations in Sect. 5.3, P = {β 1 , . . . , β Q }. Since P + ∩ P − = ∅, assume that there is a positive integer 1 ≤ Q ≤ Q such that For a visualization please refer to the right picture in Fig. 6. From (28), the solution can be written as finite number of propagation modes and a curve integral. If for some j ∈ {1, 2, . . . , Q}, μ j, (β j ) > 0 for all = 1, 2, . . . , L j , then from [20], holds for all n ∈ N and sufficiently small δ j > 0. In this case, we can merge this term with the curve integral in (28). Define the following domain: where δ > 0 satisfies all the conditions in Definition 1. For the visualization we refer to the red curve on the right picture in Fig. 6. Then when n ∈ N, When j = 1, 2, . . . , Q , since there are analytic functions μ j, passing through (β j , k 2 0 ) with both positive and negative values of μ j, (β j ), the integral on a small circle with center e iβ j no longer works. Note that since μ j, depends analytically on α with non-vanishing derivatives at β j , the inverse function, denoted by β j, , exists and depends analytically on k 2 (see Chapter 16, [27]). From inverse function theorem, when μ j, (β j ) > 0, β j, (k 2 ) < β j when k 2 < k 2 0 and β j, (k 2 ) > β j when k 2 > k 2 0 . Similar for the case that μ j, (β j ) < 0. Moreover, when k 2 = k 2 0 , β j, are separated for different thus the integral representation works in this case.
Here the balls B ± j are defined by: Furthermore, we suppose that δ 0 is sufficiently small such that B ± j ∩ S = ∅. Please refer to the black circles in the right picture in Fig. 6. From above arguments, when For j = 1, 2, . . . , Q , define u + j,k 2 by: Then from a similar arguments in [20], holds uniformly when k 2 lies in a sufficiently small neighourhood of k 2 0 . Thus As this function depends real analytically on k 2 , the LAP solution is the limit of (43) as k 2 → k 2 0 , where the second term is defined in (42).

Numerical method
Now we introduce the numerical method based on (43). The first term with k 2 = k 2 0 can be evaluated directly by the numerical method introduced in (33). For the second term, motivated by [5], an interpolation technique with respect to real valued k 2 is introduced to approximate the exact value. We conclude the algorithm as follows.
1. Compute the first term in (43) with k 2 = k 2 0 from the method introduced in (35), denoted by u 0,N ,h with the same parameters N and h. The following error estimations can be obtained in the same way as in Theorem 4: Note that C does not depend on N and h, but it depends on k 2 and may blow up when k 2 is close to k 2 0 . This is due to the asymptotic behavior of the poles when k 2 → k 2 0 : the poles approach exp(iα j ) which lie on the integral contours B ± j . To obtain the error between u k 2 0 ,N ,h,M and u k 2 0 , we only need to estimate the error from the interpolation.
As u + j,k 2 depends real analytically on k 2 ∈ K , there is a point K * ∈ K such that and there is a C 0 > 0 such that u j H 1 (Ω 0 ) ≤ C 0 uniformly for = 0, 1, 2, . . . . When we approximate u + j,k 2 by interpolation, from standard error estimation of interpolation, where |K | = max x =y∈K |x − y|. Thus we can finally obtain the error estimation of the algorithm: where = 0, 1. However, in numerical results, the convergence rate of the third term is difficult to analyze. On one hand, to make sure that the Taylor expansion converges, we require that |K | is sufficiently small; on the other hand, the error becomes larger when k 2 → k 2 0 due to the pole at k 2 0 , as C = C(|K |) may blow up. Moreover, C 0 can be a large number and it is impossible to be evaluated.
The point a 0 = (0, 0.5), and ζ(t) is a C 8 -continuous function defined by First, we draw the dispersion diagram and the corresponding z-space, see Fig. 8. From the dispersion diagram, we find out two stop bands, i.e., (2.956 ± 0.01, 7.574 ± 0.01) and (13.41±0.01, 15.49±0.01). When k 2 = 5, it lies in the first stop band, thus there is no eigenvalue on S 1 in the z-space, i.e., Q = 0. In this case, u is represented by (18) with the integral on the unit circle S 1 . However, when k 2 = 17, there are two points lying on the dispersion curve with μ n (α) = 17. This implies that Q = 1 and α + 1 = 0.9576(±0.01) (blue diamond) and α − 1 = −0.9576(±0.01) (red square). Thus we design the integral curve in (25) as the red curve on the right. Moreover, we also check the condition number of the matrix obtained from the finite element discretization of (14), to find out a rough location of poles of the function w(z, ·) (see Fig. 9). The parameter N 0 is fixed to be 6 for all the numerical examples in this section.

Full-guide problems
In the first part of this section, we show some numerical examples for scattering problems in full-waveguide, when Assumption 1 and 2 are both satisfied. The compactly supported source term f is defined as follows: where a 0 = (0.1, 0.4) .
With these data, we calculate the value of u for different parameters. For the finite element method, we choose h = 0.0025, 0.005, 0.01, 0.02 and N = 8, 16, 32, 64 for k 2 = 5, N = 4, 8, 16, 32, 64 for k 2 = 17. We also compute "exact solutions" u exa from the finite element method introduced in [5,6] with h = 0.005 and the Lagrangian element. First we show the relative errors with different h and N for both cases, defined where u N ,h is the numerical solution with parameter N and h. Note that for k 2 = 5, the CCI-method and PM-method are the same, and the results are shown in Table 1.
For k 2 = 17, the results from the CCI-method are shown in Table 2 and from the PM-method are in 3. We can see that the relative error decays as N gets larger and h gets smaller. Note that when N is sufficiently large (e.g., N ≥ 32), the relative error does not decay when N gets larger, this implies that the error brought by N is relatively smaller, compared with the error from h. The decay rate of the CCI-method and the PM-method are relevant.
As the convergence rate with respect to h is classical, we are especially interested in that of N . We fix h = 0.01 for both cases, and compute the relative error between u N ,h and u 256,h for k 2 = 5, u 128,h for k 2 = 17. From the result in (32), the error is expected to decay at the rate of O(N −6 ). From the two pictures in Fig. 10, the convergence is even faster than expected. We also compute the energy fluxes of propagating Bloch wave solutions for this example. We approximate the integrals with the same method. Fix parameters r = 0.1 and N = 16, and the mesh size h = 0.005. We obtain the values From above results, we conclude that u + 1 is propagating to the right, and u − 1 is propagating to the left according to the energy criteria, thus it coincides with the analysis in Sect. 6.3.

Half-guide problems
In the second part of this section, we show some numerical examples for half guide problems. The boundary data is given by For all the examples, we approximate the source term by where γ is the regularization parameter. As the numerical results also depend on γ , we show different results with respect to different regularization parameters, i.e., R. Zhang    Table  4 and with γ = 10 −5 in Table 5. For k 2 = 17, we show results for N = 4, 8, 16 and h = 0.02, 0.01, 0.005 with γ = 10 −2 in Table 6, and with γ = 10 −5 in Table  7. We also give the contour map for the solution with k 2 = 17, γ = 10 −5 , N = 16 and h = 0.05 in Fig. 11. For all these cases, we can see that the error decays when N gets larger (especially when h is small) and h gets smaller (especially when N is large). However, the decaying rate slows down significantly when the parameter N becomes sufficiently large (e.g., N ≥ 16). This comes from the cut-off approximation of the series of f and the regularization process. We also notice that the relative errors corresponding to γ = 10 −2 is larger than that to γ = 10 −5 , which is also as expected.

Special wavenumbers
In this section, we show some numerical results when Assumption 2 is not satisfied, i.e., S + ∩ S − = ∅. From the dispersion diagram, i.e., Fig. 12, when k 2 = 9.6663(±0.01), Assumption 2 is not satisfied. However, we can not find the exact value of k 2 such that Assumption 2 is not satisfied, but a value very close to the exact one. Then the method introduced in Sect. 8 and the PM method is used for the numerical simulation. From the dispersion curve, Q = Q = 2, α + 1 = α − 2 = −1.8333(±0.01) and α + 2 = α − 1 = 1.8333(±0.01). Let the curve C 0 := z ∈ C : |z| = 0.8 , and B ± j ( j = 1, 2) be defined by (41) with δ 0 = 0.1. For the visualization of the points and curves we refer to Fig. 12.
For the CCI-method, we choose two different interpolation strategies to carry out the numerical approximation. For the first strategy, set M = 3, and   For the second one, M = 5 and We still use the result obtained by the method introduced in [5,6] to produce the "exact solution", and compute the relative errors with the parameters h = 0.02, 0.01, 0.005, 0.0025 and N = 16, 32, 64, 128. The results are shown in Table  8-9. In both tables, the relative errors with these two different strategies are similar, and the error decays when N gets larger and h gets smaller. These results show that the CCI-method for special numbers is convergent. We also used the PM-method to approximate the LAP solution with this special wavenumber. We choose parameters h = 0.006, 0.003, 0.0015, 0.00075 and N = 16, 32, 64, 128, and the results are shown in Table 10. From the results, the relative errors decay as N gets larger and h gets smaller, but the decay is not as stable as the CCI-method.
We also check the energy fluxes corresponding to the propagating modes. The approximation is carried out with the help of the second strategy, with parameters h = 0.005 and N = 64. The energy fluxes corresponding to u ± 1 and u ± 2 are evaluated as follows: From the values, u + 1 and u + 2 are propagating to the right, u − 1 and u − 2 are propagating to the left. This coincides with the results shown in Sect. 6.2.

Conclusion
Now we compare the two different methods -the CCI-method and the PM-method. The CCI-method depends on a simplified integral representation (25) for the LAP solution with Assumption 1, where a suitable complex integral curve is to be designed specially depending on the behaviour of the Floquet multipliers with respect to the absorbing parameter . Then the LAP solution is approximated by the sum of finite number of solutions of quasi-periodic problems, which are all well-posed. However, to know the behaviour of the poles, for example by producing a dispersion diagram, may take a relatively longer time. When Assumption 1 no longer holds, an interpolation technique is introduced to make the CCI-method suitable for this situation. This makes the CCI-method more complicated. But this method does not require very accurate evaluations of the exceptional values and Bloch wave solution. On the other hand, the PM-method is based on the curve integral with finite number of non-selfadjoint eigenvalue problem. Compared to the CCI-method, it does not depend on Assumption 1. We can decide if a Floquet mode is acceptable by the sign of the energy flux, so we do not need to know the dispersion curve in principle. However, as the non-selfadjoint eigenvalue problems have complex eigenvalues, sometimes it may not be easy to find out all real eigenvalues in (−π, π]. Moreover, the errors in computing eigenvalues, eigenfunctions and derivatives of the dispersion curves are relatively larger, thus we need a much finer mesh to have a solution with similar accuracy of solutions from the CCI-method. Note that for the full waveguide problem with k = √ 17, the results from the PM-method is relevant to the CCI-method. The reason may be the almost orthogonality between f and the eigenfunctions. We also compare the methods introduced in this paper with other methods. The computational complexity of both methods are equivalent to that of [14]. The problems are different, for the Floquet-Bloch transformed field has finite number of poles in this paper, but one or two branch cuts in [14]. The methods introduced in [8,9,15] are based on the numerical evaluation of the DtN maps, which are described by the quadratic characteristic equation. The evaluations are carried out by the iteration method based on the cell problem, thus this may involve many times of the solutions of quasi-periodic problem. Another interesting method is introduced in [18], by approximating the LAP solution by finite number of propagating modes and a truncated problem. Suppose the truncated problem exists in the cells from −N to N , and the degree of freedom is M in one cell, then the degree of freedom for the whole problem is greater than N M. Thus they have to solve a system of at least N M × N M. However, for our problem we only need to solve several times M × M linear system, which is much more efficient. Due to the super algebraic convergence, we do not need to solve the M × M linear system too many times. Thus our method is faster than the one introduced in that paper. then there is a sequence α n ∈ (−π, π] such that μ(α n ) = 0, ∀n ∈ N.
As μ i 0 and μ j 0 are both analytic functions, μ is analytic as well. Thus either μ is a constant function equals to 0, or α n = α 0 except for a finite number of n's.
For the second case, suppose there is an N >> 1 such that α n = α 0 for any n ≥ N , then μ i 0 (α n ) = μ j 0 (α n ) = k 2 n implies that k 2 n = k 2 0 for any n ≥ N . This contradicts with the monotone decreasing property. Thus k 2 0 can not be an accumulation point, the proof is finished.