Numerical methods for scattering problems in periodic waveguides

Abstract

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.

1 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 always true, the Limiting Absorption Principle (LAP) is commonly used to find out the “physically meaningful” solution. That is, the “physically meaningful solution” is assumed to be the limit of a family of solutions with absorbing material, when the positive absorption parameter tends to 0. In this paper, the limit is called the LAP solution. It has been proved that LAP solutions exist for planar waveguides filled up with periodic material, and we refer to [1,2,3,4] for different proofs.

In recent years, some numerical methods have been developed to solve this kind of problems based on the LAP. As the problem with absorbing medium is uniquely solvable and the unique solution decays exponentially along the periodic waveguide, the solution is easily approximated by problems defined in bounded domains and then computed by standard numerical methods. We refer to [6, 23] for a numerical method with the help of a so-called recursive doubling process to approximate Robin-to-Robin maps on left-and right-boundaries of a periodic cell, and the original solution is approximated by the one with a sufficiently small positive absorption parameter. Based on this idea, the solution is also approximated by an extrapolation technique with respect to small positive absorption parameters in [5], and the method is also extended to scattering problems in locally perturbed periodic layers in [7]. On the other hand, the LAP has also been applied to pick up “proper” propagation modes. With the asymptotic behaviour of propagating modes when the absorption parameter tends to 0, it is possible to determine whether certain modes propagate to the left or to the right. With these propagating modes, Dirichlet-to-Neumann maps on the boundaries of periodic cells are approximated and the solution for the whole problem is computed. For details we refer to [8, 15]. For other numerical methods we also refer to [16,17,18].

Recently, the Floquet-Bloch transform has been applied to both theoretical analysis and numerical simulation for scattering problems in periodic structures. We refer to [9, 10, 19] for scattering problems with (perturbed) periodic media, to [11,12,13,14] for problems with periodic surfaces, and to [2,3,4] for periodic waveguides. From the Floquet-Bloch theory, the unique solution of the periodic waveguide problem with absorption is written as a contour integral on the unit circle, where the integrand is a family of quasi-periodic solutions depending analytically on the quasi-periodicities. When the absorbing parameter tends to 0, the integrand may become irregular as poles may approach the unit circle during this process. Based on [1, 8, 15], we replace the unit circle by a small modification of it. For the choice of the new curve, we require that the integral does not change for any sufficiently small absorbing parameter, and the integrand does not have any poles on the new curve when the parameter is 0. We also show that when the wavenumber satisfies certain conditions, such a curve always exists. Luckily, we have proved that the wavenumbers such that the conditions are not satisfied compose a discrete subset of $\mathbb {R}_+$. Thus for any positive wavenumber except for this discrete set, we can write the LAP solution as a contour integral on a closed curve, where the integrand is a family of quasi-periodic solutions of cell problems. As the quasi-periodic cell problems are classical, the numerical simulation of LAP solutions is also easily carried out. The curve can be easily chosen as a piecewise analytic one, as there are only finite number of eigenvalues on the unit circle. Finally, as the solution depends piecewise analytically on the curve, a high order numerical method is designed based on the contour integral. Moreover, we can also extend this method to more general wavenumbers, with an interpolation technique inspired by the paper [5].

The numerical method is also extended to halfguide problems. From [20], any LAP solution of a halfguide problem is (not uniquely) extended to an LAP solution of a fullguide problem. Thus the numerical method can be designed as two steps. The first step is to find out the source term for the fullguide problems such that its solution approximates that of the halfguide problem, for given boundary data. We compute the boundary values for basis functions in the space where the source term lies in, and then find out the corresponding coefficients by the least square method that approximates the boundary data. As the choice of source terms is not unique and the problem is severely ill-posed, a Tikhonov regularization technique is adopted. With this source term, we go to the second step, i.e., to solve the fullguide problem with the method developed in this paper and approximate the solution of the halfguide problem. Moreover, the methods can also be extended to more general wavenumbers in the same way.

The rest of this paper is organized as follows. We recall the mathematical model of the scattering problem in the second section, and introduce the Floquet-Bloch theory in the third section. Then we consider the quasi-periodic cell problems in the fourth section. In Sect. 5, we apply the Floquet-Bloch transform to obtain a simplified integral representation of the LAP solutions. In Sect. 6, we develop a numerical method for LAP solutions, based on the integral representation. Then we extend the method to halfguide problems in Sect. 7 and more general wavenumbers in Sect. 8. In the last section, we present numerical results to show the efficiency of our algorithms.

2 The mathematical model of direct scattering problems

Let $\varOmega =\mathbb {R}\times (0,1)\subset \mathbb {R}^2$ be a closed waveguide (see Fig. 1) with boundary

$$\begin{aligned} \partial \varOmega =\{x\in \mathbb {R}^2:\,x_2=0\text { or }x_2=1\}. \end{aligned}$$

The boundary is assumed to be impenetrable. Suppose $\varOmega$ is filled up with periodic material with a real-valued refractive index q satisfying:

$$\begin{aligned} q(x_1+1,x_2)=q(x_1,x_2),\quad q\ge c_0>0 \quad \forall \,x\in \varOmega . \end{aligned}$$

Moreover, we require that $q\in L^\infty (\varOmega )$.

Fig. 1

Periodic waveguide

The scattering problem is modeled by the following equations:

$$\begin{aligned} \varDelta u+k^2 q u= & {} f\quad \text { in }\varOmega ; \end{aligned}$$

(1)

$$\begin{aligned} \frac{\partial u}{\partial x_2}= & {} 0\quad \text { on }\partial \, \varOmega , \end{aligned}$$

(2)

where f is a function in $L^2(\varOmega )$ with compact support, and $k>0$ is the real wavenumber.

Remark 1

We can also consider different boundary conditions on $\partial \,\varOmega$, 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 $\epsilon >0$, consider the following damped Helmholtz equation:

$$\begin{aligned} \varDelta u_\epsilon +(k^2+\mathrm {i}\epsilon )q u_\epsilon= & {} f\quad \text { in }\varOmega ; \end{aligned}$$

(3)

$$\begin{aligned} \frac{\partial u_\epsilon }{\partial {x_2}}= & {} 0\quad \text { on }\partial \,\varOmega . \end{aligned}$$

(4)

From the Lax-Milgram theorem, the problem is uniquely solvable in $H^1(\varOmega )$. Moreover, the solution $u_\epsilon$ decays exponentially as $x_1\rightarrow \pm \infty$ (see [8]). The LAP assumes that $u_\epsilon$ converges in $H^1_{loc}(\varOmega )$ when $\epsilon \rightarrow 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:

$$\begin{aligned} {\mathcal {A}}u:=-q^{-1} \varDelta u {\text { with } D({\mathcal {A}},\varOmega ):=\Big \{u\in H^1(\varOmega ):\,\varDelta u\in L^2(\varOmega ),\,\left. \frac{\partial u}{\partial x_2}\right| _{\partial \,\varOmega }=0\Big \}}. \end{aligned}$$

Let the spectrum of ${\mathcal {A}}$ be denoted by $\sigma ({\mathcal {A}})$, then the problem (1)–(2) is uniquely solvable in $H^1(\varOmega )$ if and only if $k^2\notin \sigma ({\mathcal {A}})$. To study the spectrum property of ${\mathcal {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].

3 Floquet-Bloch theory and quasi-periodic problems

3.1 Floquet-Bloch theory

In this section, we introduce the Floquet-Bloch theory to study the spectrum $\sigma (A)$. For simplicity, we introduce the following domains (see Fig. 1):

$$\begin{aligned} \varOmega _j:= (-1/2+j,1/2+j]\times (0,1),\quad \Gamma _j=\{-1/2+j\}\times (0,1),\quad j\in \mathbb {Z}. \end{aligned}$$

Then $\varOmega =\bigcup _{j\in \mathbb {Z}}\varOmega _j$ with the left and right boundaries $\Gamma _j$ and $\Gamma _{j+1}$. Let

$$\begin{aligned} \partial \,\varOmega _j:=\partial \,\varOmega \cap \overline{\varOmega _j}=\Big \{x\in \mathbb {R}^2:\,-1/2+j<x_1\le 1/2+j,\,x_2=0\text { or }1\Big \}. \end{aligned}$$

We also introduce the space of quasi-periodic functions. A function $\phi \in H^1_{loc}(\varOmega )$ is called z-quasi-periodic, if it satisfies

$$\begin{aligned} \left. \phi \right| _{\Gamma _{j+1}}=z\, \left. \phi \right| _{\Gamma _j} ,\quad \forall j\in \mathbb {Z}\end{aligned}$$

(5)

for some fixed complex number $z\in \mathbb {C}$. We define the subspace of $H^1(\varOmega _0)$ by:

$$\begin{aligned} H_z^1(\varOmega _0):=\Big \{\phi \in H^1(\varOmega _0):\,\phi \text { satisfies } (5) \hbox { for }j=0\Big \}. \end{aligned}$$

Then the functions in $H_z^1(\varOmega _0)$ can be extended to z-quasi-periodic functions. Especially, when $z=1$, all functions in $H^1_1(\varOmega _0)$ can be extended into periodic functions in $H^1_{loc}(\varOmega )$. We also denote $H^1_1(\varOmega _0)$ by $H^1_{per}(\varOmega _0)$.

From the Floquet-Bloch theory, the spectrum of ${\mathcal {A}}$ is closely related to Bloch wave solutions. A Bloch wave solution is a non-trivial solution of (1)–(2) in $H^1_z(\varOmega _0)$ with $f=0$ for some $z\in \mathbb {C}$. If a Bloch wave solution exists in $H^1_z(\varOmega _0)$, z is called a Floquet- multiplier. Define the operator:

$$\begin{aligned} {\mathcal {A}}_z u=-q^{-1}\,\varDelta u\,\text { {with domain} }D_z({\mathcal {A}},\varOmega _0):=D({\mathcal {A}},\varOmega _0)\cap H_z^1(\varOmega _0), \end{aligned}$$

(6)

where $D({\mathcal {A}},\varOmega _0)$ is defined in the same way as $D({\mathcal {A}},\varOmega )$, and $\varOmega$ is replaced by the periodic cell $\varOmega _0$. Note that, ${\mathcal {A}}_z$ is self-adjoint with respect to the $L^2$-space equipped with the weighted inner product $(\phi ,\psi )_{L^2,q}=\int _{\varOmega _0}q\phi \overline{\psi }\,\mathrm {d}x \,$. Let $\sigma ({\mathcal {A}}_z)$ be the spectrum of ${\mathcal {A}}_z$, then $k^2\in \sigma (A_z)$ if and only if z is a Floquet multiplier.

Let $\mathbb F(k^2)$ be the collection of all Floquet multipliers with wavenumber k and ${\mathbb U\mathbb F}(k^2):=\mathbb {F}(k^2)\cap \mathbb {S}^1$ ($\mathbb {S}^1$ is the unit circle in $\mathbb {C}$) be the set of all unit Floquet multipliers. In this paper, when the wavenumber is fixed, we write $\mathbb {F}$ instead of $\mathbb 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]:

• $\mathbb {UF}$ is finite.

• $z\in {\mathbb F}$ if and only if $z^{-1}\in {\mathbb F}$, thus $z\in {\mathbb {UF}}$ if and only if $\overline{z}=z^{-1}\in {\mathbb {UF}}$.

• $\mathbb {F}$ is a discrete set and the only finite accumulation point of $\mathbb {F}$ can be 0.

• $\mathbb F(k^2)$ depends continuously on $k^2$.

Fig. 2

Dispersion diagram. Left: Case I; Right: Case II

A classical result from the Floquet-Bloch theory also shows that (see [25]):

$$\begin{aligned} \sigma ({\mathcal {A}})=\bigcup _{|z|=1}\sigma ({\mathcal {A}}_z). \end{aligned}$$

(7)

Thus, it is particularly important to study the spectrum of ${\mathcal {A}}_z$ when $|z|=1$. For simplicity, let $\alpha =-\mathrm {i}\log (z)$ where $\alpha \in (-\pi ,\pi ]$. We replace ${\mathcal {A}}_z$ by ${\mathcal {A}}_\alpha$ in the rest of this section, by abuse of notation. Then (5) becomes

$$\begin{aligned} \left. u\right| _{\Gamma _{j+1}}=\exp (\mathrm {i}\alpha ) \left. u\right| _{\Gamma _j}, \quad \forall j\in \mathbb {Z}. \end{aligned}$$

(8)

Denote the spectrum of ${\mathcal {A}}_\alpha$ by $\sigma ({\mathcal {A}}_\alpha )$. As ${\mathcal {A}}_\alpha$ is self-adjoint, $\sigma ({\mathcal {A}}_\alpha )$ is a discrete subset of $(0,\infty )$. By rearranging the order of the points in $\sigma ({\mathcal {A}}_\alpha )$ properly, we obtain a family of analytic functions $\{\mu _n(\alpha ):\,n\in \mathbb {N}\}$ and $\{\psi _n(\cdot ,\alpha ):\,n\in \mathbb {N}\}$:

$$\begin{aligned} {\mathcal {A}}_\alpha \psi _n(\cdot ,\alpha )=\mu _n(\alpha )\psi _n(\cdot ,\alpha ),\quad \sigma ({\mathcal {A}}_\alpha )=\bigcup _{n\in \mathbb {N}}\{\mu _n(\alpha )\}. \end{aligned}$$

Thus $\sigma ({\mathcal {A}})=\bigcup _{n\in \mathbb {N}}\bigcup _{\alpha \in (-\pi ,\pi ]}\big \{\mu _n(\alpha )\big \}$. Both $\mu _n(\alpha )$ and $\psi _n(\cdot ,\alpha )$ are extended into analytic functions in $\alpha$ in a sufficiently small neighbourhood of $(-\pi ,\pi )\times \{0\}$.

For any fixed $n\in \mathbb {N}$, the graph $\big \{(\alpha ,\mu _n(\alpha )):\,\alpha \in (-\pi ,\pi ]\big \}$ 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:

  $$\begin{aligned} \mu _{jm}(\alpha )=j^2\pi ^2+(\alpha +2\pi m)^2,\quad j\in \mathbb {N},\,m\in \mathbb {Z}. \end{aligned}$$

• 2. Example 2. $q=9$ in a disk $B\big ((0,0.5),0.3\big )$ 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(\varOmega )$. 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 $\alpha \in (-\pi ,\pi ]$ such that there is a non-trivial $\alpha$-quasi-periodic function $\psi$ satisfying ${\mathcal {A}}_\alpha \psi =k^2\psi$. 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 $\alpha \in (-\pi ,\pi ]$ such that $k^2\in \sigma (A_\alpha )$. Then the set

$$\begin{aligned} P :=\left\{ \alpha \in (-\pi ,\pi ]:\,\exists \, n\in \mathbb {N},\,\mathrm{s.t.,}\, \mu _n(\alpha )=k^2\right\} \ne \emptyset . \end{aligned}$$

Thus $\mathbb {UF}$ can be written as:

$$\begin{aligned} \mathbb {UF}=\left\{ \exp (\mathrm {i}\alpha ):\,\alpha \in P\right\} . \end{aligned}$$

The points in P are divided into the following three classes:

• When $\mu '_n(\alpha )>0$, $\psi _n(\cdot ,\alpha )$ propagates from the left to the right;

• when $\mu '_n(\alpha )<0$, $\psi _n(\cdot ,\alpha )$ propagates from the right to the left;

• when $\mu '_n(\alpha )=0$, we can not decide the direction that $\psi _n(\cdot ,\alpha )$ propagates.

For physical interpretations of the Floquet modes $\psi _n(\cdot ,\alpha )$ we refer to Remark 4, [2]. Based on the above classification, we define the following three sets:

$$\begin{aligned} P_\pm:= & {} \left\{ \alpha \in (-\pi ,\pi ]:\,\exists \, n\in \mathbb {N}\text { s.t., }\mu _n(\alpha )=k^2\text { and }\pm \mu '_n(\alpha )>0\right\} ;\\ P_0:= & {} \left\{ \alpha \in (-\pi ,\pi ]:\,\exists \, n\in \mathbb {N}\text { s.t., }\mu _n(\alpha )=k^2\text { and }\mu '_n(\alpha )=0\right\} . \end{aligned}$$

Then $P =P_+ \bigcup P_- \bigcup P_0$.

Remark 2

It is possible that there are two (or more) different dispersion curves passing through the point $(\alpha ,k^2)$. Suppose the elements in P have Q different values $\alpha _1,\dots ,\alpha _Q$. For any $j=1,2,\dots ,Q$, there are $L_j$ ($L_j\ge 1$) different dispersion curves $\mu _{j,\ell }$ passing through $(\alpha ,k^2)$ ($\ell =1,2,\dots ,L_j$) such that

$$\begin{aligned} \mu _{j,1}(\alpha _j)=\mu _{j,2}(\alpha _j)=\cdots \mu _{j,L_j}(\alpha _j)=k^2. \end{aligned}$$

In this case, $\alpha _j$ is treated as $L_j$ different elements in P, i.e.,

$$\begin{aligned} \alpha _{j,1}=\alpha _{j,2}=\cdots =\alpha _{j,L_j}=\alpha _j. \end{aligned}$$

As $\mathbb {UF}$ is symmetric, $P_\pm$ is also symmetric, i.e., $\alpha \in P_+$ if and only if $-\alpha \in 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 =\emptyset$.

The assumption is reasonable as the set $\big \{k>0:\,P_0 \ne \emptyset \big \}$ is “sufficiently small”, i.e., the set is countable with at most one accumulation point at $\infty$ (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_+ \cap P_- \ne \emptyset$. 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.

Lemma 1

The set $\big \{k\in \mathbb {R}_+:\,P_+ \cap P_- \ne \emptyset \big \}$ is countable, and its only accumulation point is $\infty$.

Assumption 2

Until Sect. 7, we assume that k satisfies $P_+ \cap P_- =\emptyset$.

We define three subsets of $\mathbb {UF}$ from the definition of $P_\pm$ and $P_0$ by:

$$\begin{aligned} S_\pm ^0:=\left\{ z=\exp (\mathrm {i}\alpha ):\,\alpha \in P_\pm \right\} ,\quad S_0^0:=\left\{ z=\exp (\mathrm {i}\alpha ):\,\alpha \in P_0 \right\} . \end{aligned}$$

(9)

Then $\mathbb {UF}=S_+^0\bigcup S_-^0\bigcup S_0^0$.

Fig. 3

Example for $n=1$ and $k^2=3\pi ^2$. $\mathbb {UF}$ in $\alpha$-space and z-space. Red squares denotes the points in $P_-$ ($\S _-^0$), while blue diamonds denotes the points in $P_+$ ($\S _+^0$)

We also divide the set $\mathbb {F}{\setminus }\mathbb {UF}$ into the following two subsets:

$$\begin{aligned} RS:=\left\{ z\in \mathbb {F}:\,|z|<1\right\} ;\quad LS:=\left\{ z\in \mathbb {F}:\,|z|>1\right\} . \end{aligned}$$

The Bloch wave solution corresponding to $z\in RS$ is evanescent on the right, while the one corresponding to $z\in LS$ is evanescent on the left. Moreover, let

$$\begin{aligned} S_+:=S_+^0\bigcup RS,\quad S_-:=S_-^0\bigcup LS. \end{aligned}$$

Remark 3

We conclude the properties of the sets RS, LS and $S_\pm ^0$ from the properties of $\mathbb {F}$ as follows:

• From the symmetry of $\mathbb {F}$, the sets $S_+^0$ and $S_-^0$, RS and LS are symmetric, i.e.,

  $$\begin{aligned} z\in S_+^0\iff z^{-1}=\overline{z}\in S_-^0;\quad z\in RS\iff z^{-1}\in LS. \end{aligned}$$

  (10)

  This also implies that

  $$\begin{aligned} z\in S_+\iff z^{-1}\in S_-. \end{aligned}$$

• When Assumption 1 is satisfied, $S_0^0=\emptyset$, thus $S_+^0\bigcup S_-^0=\mathbb {UF}$.

• When Assumption 2 is satisfied, $S_+^0\cap S_-^0=\emptyset$.

• From (10), if $\pm 1\in S_+^0$ (or $\pm 1\in S_-^0$), then $\pm 1\in S_-^0\cap S_+^0$. If Assumption 2 is satisfied, $S_+^0\cap S_-^0=\emptyset$, then $\pm 1\notin \mathbb {UF}{\setminus } S_0^0$. If Assumption 1 is also satisfied, $S_0^0=\emptyset$ implies that $\pm 1\notin \mathbb {F}$.

Lemma 2

Let $k>0$. There is a $\tau >0$ such that $RS\in B(0,\exp (-\tau ))$ and $LS\in \mathbb {C}{\setminus }\overline{B(0,\exp (\tau ))}$.

We omit the proof here. For details we refer to [20].

3.2 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:

$$\begin{aligned} \varDelta u_z+k^2 qu_z= & {} f_z\quad \text { in }\varOmega _0; \end{aligned}$$

(11)

$$\begin{aligned} \frac{\partial u_z}{\partial {x_2}}= & {} 0\quad \text { on }\partial \varOmega \cap \overline{\varOmega _0}; \end{aligned}$$

(12)

$$\begin{aligned} u|_{\Gamma _1}-z\, u|_{\Gamma _0}= & {} 0; \end{aligned}$$

(13)

where $z\ne 0$ is a complex number and $f_z\in L^2(\varOmega _0)$ depends analytically on z.

We are interested in the well-posedness of the problem (1112)–(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. Let $v_z=z^{-x_1}u_z(x)$, then $v_z\in H_{per}^1(\varOmega _0)$. 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 $\mathbb {C}_\times :=\{z\in \mathbb {C}{\setminus }\{0\}:\,-\pi <\arg (z)\le \pi \}$, where $\arg (z)$ is the argument of the complex number z). From direct calculation, $v_z$ is the solution of the following problem:

$$\begin{aligned} \varDelta v_z+2\log (z)\frac{\partial v_z}{\partial x_1}+(k^2q+\log ^2(z))v_z= & {} z^{-x_1} f_z\quad \text { in }\varOmega _0;\\ \frac{\partial v_z}{\partial \nu }= & {} 0\quad \text { on }\partial \,{\varOmega _0};\\ v_z|_{\Gamma _1}- v_z|_{\Gamma _0}= & {} 0. \end{aligned}$$

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].

$$\begin{aligned} (I-{\mathcal {K}}_z)v_z=z^{-x_1}\widetilde{f}_z. \end{aligned}$$

This implies that when $I-{\mathcal {K}}_z$ is invertible, then

$$\begin{aligned} v_z=\left( I-{\mathcal {K}}_z\right) ^{-1}z^{-x_1}\widetilde{f}_z. \end{aligned}$$

(14)

As ${\mathcal {K}}_z$ is compact and depends analytically on $z\in \mathbb {C}_\times$, $I-{\mathcal {K}}_z$ is an analytic family of Fredholm operators. From the Analytic Fredholm Theory (see Theorem VI.14, [26]), $(I-{\mathcal {K}}_z)^{-1}$ depends meromorphically on z in $\mathbb {C}{\setminus }\{0\}\}$. From Sect. 4 in [20], the set of poles of $\left( I-{\mathcal {K}}_z\right) ^{-1}$ is exactly $\mathbb {F}$. Thus $v_z$, or equivalently $u_z=z^{x_1}v(x)$, depends analytically on $z\in \mathbb {C}{\setminus }\!\left( \mathbb {F}\bigcup \{0\}\right)$ and meromorphically on $z\in \mathbb {C}{\setminus }\{0\}$.

4 The Floquet-Bloch transform and its application

4.1 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 $\phi \in C_0^\infty (\varOmega )$, define the Floquet-Bloch transform of $\phi$ by

$$\begin{aligned} (\mathcal {F}\phi )(z,x):=\sum _{n=-\infty }^\infty \phi (x_1+n,x_2)z^{-n},\quad x\in \varOmega _0,\,z\in \mathbb {C}. \end{aligned}$$

(15)

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 $\mathbb {C}$ such that the series (15) converges. Note that the DOC may be empty for given function $\phi$. When $\phi$ decays exponentially at the rate $\gamma$, i.e., there is a $\gamma >0$ and $C>0$ such that $\phi$ satisfies

$$\begin{aligned} |\phi (x_1,x_2)|\le C\exp (-\gamma |x_1|),\quad \forall \,x\in \varOmega , \end{aligned}$$

(16)

the Floquet-Bloch transform of $\phi$ is still well-defined, and the ROC is the annulus

$$\begin{aligned} T_\gamma =\big \{z\in \mathbb {C}:\,\exp (-\gamma )<|z|<\exp (\gamma )\big \}. \end{aligned}$$

Moreover, from the perturbation theory, the function $(\mathcal {F}\phi )(z,\cdot )$ depends analytically on $z\in T_\gamma$. It is also easy to check that the transformed function $(\mathcal {F}\phi )(z,\cdot )$ is quasi-periodic (i.e., it satisfies (5)). We conclude some mapping properties of the operator $\mathcal {F}$ in the following proposition.

Proposition 1

The operator $\mathcal {F}$ has the following properties when z lies on the unit circle $\mathbb {S}^1$ (see [11, 21]):

• $\mathcal {F}$ is an isomorphism between $H^s(\varOmega )$ and $L^2(\mathbb {S}^1;H^s_z(\varOmega _0))$ ( $s\in \mathbb {R}$), where

  $$\begin{aligned} L^2(\mathbb {S}^1;H^s_z(\varOmega _0)):=\left\{ \phi \in \mathcal {D}'(\mathbb {S}^1\times \varOmega _0):\,\left[ \int _{\mathbb {S}^1}\left\| \phi (z,\cdot )\right\| ^2_{H^s_z(\varOmega _0)}\,\mathrm {d}z \,\right] ^{1/2}<\infty \right\} . \end{aligned}$$

• $\mathcal {F}\phi$ depends analytically on $z\in T_\gamma$, if and only if $\phi$ decays exponentially with the rate $\gamma$.

• Given $\psi (z,x):=(\mathcal {F}\phi )(z,x)$ for some $\phi \in H^s(\varOmega )$ and satisfies (16) for some $\gamma >0$, the inverse operator $\mathcal {F}$ is given by:

  $$\begin{aligned} (\mathcal {F}^{-1}\psi )(x_1+n,x_2)=\frac{1}{2\pi \mathrm {i}}\oint _{\mathbb {S}^1}\psi (z,x)z^{n-1}\,\mathrm {d}z \,. \end{aligned}$$

  (17)

4.2 Application of the Floquet-Bloch transform

In this section, we apply the Floquet-Bloch transform $\mathcal {F}$ to the scattering problem (1)–(2), when $k^2\notin \sigma (A)$. We are particularly interested in the case:

• for $k>0$, $k^2\notin \sigma (A)$, i.e., $k^2$ lies in a stop band; or

• $k^2$ is no longer real, i.e., $k^2=k^2_0+\mathrm {i}\epsilon$ for some fixed $k_0>0$ and $\epsilon >0$.

When either of the two conditions is satisfied, the problem (1)–(2) is uniquely solvable in $H^1(\varOmega )$. Moreover, u decays exponentially at the infinity, i.e., u satisfies (16) for some $C>0$ and $\gamma >0$ (see [24]).

Remark 4

From now on, we assume that $\mathrm{supp}(f)\subset \varOmega _0$. The results are easily extended to cases when $\mathrm{supp}(f)$ lies in larger bounded domains.

We define the Floquet-Bloch transform $w(z,x):=(\mathcal {F}u)(z,x)$, then the transformed field $w(z,\cdot )$ is well-defined and depends analytically on $z\in T_\gamma$ as a $H^1(\varOmega _0)$-function valued function. It is also easy to check that for any $z\in T_\gamma$, $w(z,\cdot )\in H_z^1(\varOmega _0)$ satisfies (1112)–(13). Note that the source term in (1112) is f, as $(\mathcal {F}f )(z,x)=f(x)$ for any $z\in T_\gamma$.

From the inverse Floquet-Bloch transform and Cauchy integral theorem, the solution of the original problem can be represented as:

$$\begin{aligned} \begin{aligned} u(x_1+n,x_2)= (\mathcal {F}^{-1} w)(x_1+n,x_2)&=\frac{1}{2\pi \mathrm {i}}\oint _{\mathbb {S}^1} w(z,x)z^{n-1}\,\mathrm {d}z \,\\&=\frac{1}{2\pi \mathrm {i}}\oint _{\mathcal {C}} w(z,x)z^{n-1}\,\mathrm {d}z \,, \end{aligned} \end{aligned}$$

where $\mathcal {C}$ is a rectifiable curve in $T_\gamma$ encircling 0.

From the exponential decay of u, $(\mathcal {F}u)(z,\cdot )$ exists and depends analytically on $z\in T_\gamma$. On the other hand, from Sect. 3.2, when $z\in \mathbb {C}{\setminus }\mathbb {F}$, the problem (1112)–(13) is uniquely solvable in $H_z^1(\varOmega _0)$, and $w(z,\cdot )$ depends analytically on $z\in \mathbb {C}{\setminus }(\mathbb {F}\bigcup \{0\})$ and meromorphically on $z\in \mathbb {C}{\setminus }\{0\}$. Thus $(\mathcal {F}u)(z,\cdot )$ is extended meromorphically in $\mathbb {C}{\setminus }\{0\}$. From the analytic continuation, we obtain the following result (see [20]).

Theorem 1

When $k^2\notin \sigma ({\mathcal {A}})$, the Floquet-Bloch transformed field $(\mathcal {F}u)(z,x)$ is extended to an analytic function in $\mathbb {C}{\setminus }(\mathbb {F}\bigcup \{0\})$ and a meromorphic function in $\mathbb {C}{\setminus }\{0\}$ by the solution $w(z,\cdot )$ of (1112)–(13).

The integral representation of u is obtained from Cauchy integral theorem.

Theorem 2

Suppose $k^2\notin \sigma ({\mathcal {A}})$. $w(z,\cdot )$ is the solution of (1112)–(13) for $z\in \mathbb {C}{\setminus }\mathbb {F}$. Then the solution of (1)–(2) is written as

$$\begin{aligned} u(x_1+n,x_2)= (\mathcal {F}^{-1} w)(x_1+n,x_2)=\frac{1}{2\pi i}\oint _{\mathcal {C}} w(z,x)z^{n-1}\,\mathrm {d}z \,, \end{aligned}$$

(18)

where $\mathcal {C}\subset \mathbb {C}$ is a counter-clockwise closed rectifiable path encircling all the points in $RS(=S_+)$ and does not encircle any point in $LS(=S_-)$.

5 The Limiting absorption principle (LAP)

In this section, we consider the case when $k^2\in \sigma ({\mathcal {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:

$$\begin{aligned} \varDelta w_\epsilon (z,\cdot )+(k^2+\mathrm {i}\epsilon )qw_\epsilon (z,\cdot )= & {} f\quad \text { in }\varOmega _0; \end{aligned}$$

(19)

$$\begin{aligned} \frac{\partial w_\epsilon (z,\cdot )}{\partial {x_2}}= & {} 0\quad \text { on }\partial \,{\varOmega _0};\end{aligned}$$

(20)

$$\begin{aligned} w_\epsilon (z,\cdot )|_{\Gamma _1}-z\,w_\epsilon (z,\cdot )|_{\Gamma _0}= & {} 0. \end{aligned}$$

(21)

Similar to the definition of ${\mathcal {K}}_z$, we denote the operator with $k^2+\mathrm {i}\epsilon$ by ${\mathcal {K}}_z^\epsilon$. From Theorem 4 in [22], the poles of the operator $(I-{\mathcal {K}}_z^{\epsilon })^{-1}$ depend continuously on $\epsilon$. First, we study the asymptotic behaviour of distribution of the poles when $\epsilon >0$ is sufficiently small.

5.1 Distribution of poles of the damped Helmholtz equations

From the Floquet-Bloch theory (see [2, 25]), $k^2\in \sigma ({\mathcal {A}})$ implies that $\mathbb {UF}\ne \emptyset$. For any $z\in \mathbb {UF}$, $k^2$ is an eigenvalue of ${\mathcal {A}}_z$. As both $S_+^0$ and $S_-^0$ are finite sets, and $S_+^0$ and $S_-^0$, RS and LS are symmetric in the sense of (10), they are written as

$$\begin{aligned}&S_+^0=\{z_{j,\ell }^+:\,j=1,\dots ,Q;\,\ell =1,\dots ,L_j\};\,\, RS=\{z_{Q+1}^+,z_{Q+2}^+,\dots \}; \end{aligned}$$

(22)

$$\begin{aligned}&S_-^0=\{z_{j,\ell }^-:\,j=1,\dots ,Q;\,\ell =1,\dots ,L_j\};\,\, LS=\{z_{Q+1}^-,z_{Q+2}^-,\dots \}; \end{aligned}$$

(23)

where $z_j^+$ ($j=1,\dots ,Q$) are Q different values on the unit circle, and so are $z_j^-$ ($j=1,\dots ,Q$). Moreover, $z_{j,1}^\pm =\dots =z_{j,L_j}^\pm =z_j^\pm$ for $j=1,\dots ,Q$, and $z_j^+=\left( z_j^-\right) ^{-1}$ for any integer $j\in \mathbb {N}$. From Assumption 1 and 2, $\mathbb {UF}=S_+^0\bigcup S_-^0$, $\mathbb {F}=S_+^0\bigcup S_-^0\bigcup RS\bigcup LS$; moreover, $S_+^0\cap S_-^0=\emptyset$.

From the continuous dependence of poles, for any $z_{j,\ell }^\pm$ or $z_j^\pm \in \mathbb {F}$ with $j\in \mathbb {N}$, there is a continuous function $Z_{j,\ell }^\pm (\epsilon )$ or $Z_j^\pm (\epsilon )$, such that $\{Z_{j,\ell }^+(\epsilon ),Z_{j,\ell }^-(\epsilon ):\,j=1,\dots ,Q;\,\ell =1,\dots ,L_j\}\bigcup \{Z_j^+(\epsilon ),Z_j^-(\epsilon ):\,j\ge Q+1\}$ are exactly the set of all poles with respect to $k^2+\mathrm {i}\epsilon$. Moreover, $\lim _{\epsilon \rightarrow 0}Z^\pm _{j,\ell }(\epsilon )=Z^\pm _{j,\ell }(0)=z_j^\pm$ for $j=1,\dots ,Q$ and $\ell =1,\dots ,L_j$ and $\lim _{\epsilon \rightarrow 0}Z^\pm _{j}(\epsilon )=Z^\pm _{j}(0)=z_j^\pm$ for $j\ge Q+1$. Analogous to the case $\epsilon =0$, we define the following sets depending on $\epsilon$:

$$\begin{aligned}&S_+^0(\epsilon )=\{Z_{j,\ell }^+(\epsilon ):\,j=1,\dots ,Q;\,\ell =1,\dots ,L_j\};\\&S_-^0(\epsilon )=\{Z_{j,\ell }^-(\epsilon ):\,j=1,\dots ,Q;\,\ell =1,\dots ,L_j\};\\&LS(\epsilon )=\{Z_{Q+1}^-(\epsilon ),Z_{Q+2}^-(\epsilon ),\dots \};\\&RS(\epsilon )=\{Z_{Q+1}^+(\epsilon ),Z_{Q+2}^+(\epsilon ),\dots \}. \end{aligned}$$

From the continuous dependence of $\epsilon >0$, we have the following properties of $Z_j^\pm (\epsilon )$ when $j=1,2,\dots ,Q$. For the proof we refer to Appendix in [8].

Lemma 3

For any $j=1,2,\dots ,Q$, when $\epsilon >0$ is sufficiently small, the functions satisfy $|Z_{j,\ell }^+(\epsilon )|<1$ and $|Z_{j,\ell }^-(\epsilon )|>1$.

Thus we conclude the behaviour of $Z_j^\pm (\epsilon )$ for sufficiently small $\epsilon$:

• for any $z_j^+\in S^0_+$, the points $Z_{j,\ell }^+(\epsilon )\in S_+^0(\epsilon )$ converge to $z_j^+$ from the interior of the unit circle;

• for any $z_j^-\in S^0_-$, the points $Z_{j,\ell }^-(\epsilon )\in S_-^0(\epsilon )$ 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,\ell }^\pm (\epsilon )$ as $\epsilon \rightarrow 0$ can be seen from the picture.

For the points in $RS(\epsilon )$ or $LS(\epsilon )$, we estimate their distributions for sufficiently small $\epsilon >0$ in the following lemma (see Lemma 17, [20]).

Fig. 4

Left: the curve $\mathcal {C}_0$; Right: a choice of $\mathcal {C}$. Red solid curves: $Z_j^\pm (\epsilon )$; grey dashed curves: the unit circle. The red rectangles are points in $S^0_-$ and the blue diamonds are points in $S^0_+$. The arrows show the direction of the poles when $\epsilon \rightarrow 0$

Lemma 4

Suppose for some $\tau >0$, $RS\subset B(0,\exp (-\tau ))$ and $LS\subset \mathbb {C}{\setminus }\overline{B(0,\exp (\tau ))}$. For any $\tau _1\in (0,\tau )$, there exists $\epsilon _0>0$ such that $Z_j^+(\epsilon )\in B(0,\exp (-\tau _1))$ and $Z_j^-(\epsilon )\in \mathbb {C}{\setminus }\overline{B(0,\exp (\tau _1))}$ for any $j\ge Q+1$ and $\epsilon \in (0,\epsilon _0)$.

From Lemmas 3 and 4, when $\epsilon >0$ is sufficiently small, the sets $S_\pm ^0(\epsilon )$ and $RS(\epsilon )$, $LS(\epsilon )$ have following properties:

$$\begin{aligned} S_+^0(\epsilon )\bigcup RS(\epsilon )\subset B(0,1);\quad S_-^0(\epsilon )\bigcup LS(\epsilon )\subset \mathbb {C}{\setminus }\overline{B(0,1)}. \end{aligned}$$

(24)

5.2 Integral representation of the LAP solution

Now we are prepared to formulate the LAP solution of (1)–(2) when $k^2\in \sigma ({\mathcal {A}})$. From Theorem 2, the solution $u_\epsilon$ ($\forall \,\epsilon >0$) of the damped problem (3)–(4) is represented by the curve integral:

$$\begin{aligned} u_\epsilon (x_1+n,x_2)=\frac{1}{2\pi \mathrm {i}}\oint _{\mathbb {S}^1}w_\epsilon (z,x)z^{n-1}\,\mathrm {d}z \,. \end{aligned}$$

But when $\epsilon \rightarrow 0^+$, as the poles in $S_\pm ^0(\epsilon )$ approach $\mathbb {S}^1$, and the integral becomes singular. From Theorem 2, we are aimed at finding out a proper curve $\mathcal {C}$ to replace the unit circle such that the function $w_\epsilon (z,x)$ is well-posed and converges uniformly with respect to z when $\epsilon \rightarrow 0$. From the properties of the sets $S_\pm ^0(\epsilon )$, $RS(\epsilon )$ and $LS(\epsilon )$, we define a convenient curve $\mathcal {C}_0$ as follows.

Definition 1

Let the piecewise analytic curve $\mathcal {C}_0$ be defined as the boundary of the following domain :

$$\begin{aligned} B_\delta :=B(0,1)\bigcup \left[ \bigcup _{n=1}^Q B(z_j^+,\delta )\right] \!{\setminus }\!\left[ \bigcup _{n=1}^Q \overline{B(z_j^-,\delta )}\right] , \end{aligned}$$

where $\delta >0$ is sufficiently small such that the following conditions are satisfied:

• For any $j\ne \ell$, $B(z_j^\pm ,\delta )\cap B(z_\ell ^\pm ,\delta )=\emptyset$.

• $\left[ \bigcup _{n=1}^Q B(z_j^+,\delta )\right] \bigcup \left[ \bigcup _{n=1}^Q B(z_j^-,\delta )\right] \subset T_\tau$, i.e., the balls do not contain any point in $RS\bigcup LS$.

From Assumption 2 and Lemma 2, we can always find a proper parameter $\delta$ such that both conditions are satisfied. We refer to Fig. 4 for a visualization of $\mathcal {C}_0$ for the example when $n=1$ and $k^2=3\pi ^2$.

Thus from Lemmas 3 and 4, there is a constant $C=C(\delta )>0$ such that

$$\begin{aligned} d\left( \mathcal {C}_0,\,RS(\epsilon )\bigcup LS(\epsilon )\bigcup S_+^0(\epsilon )\bigcup S_-^0(\epsilon )\right) >C(\delta ) \end{aligned}$$

holds uniformly for any fixed sufficiently small $\epsilon >0$, where d(X, Y) is the Hausdorff distance between two subsets $X,\,Y\subset \mathbb {C}$. From the choice of the curve $\mathcal {C}_0$, the interior of the symmetric difference of $B_\delta$ and B(0, 1) is

$$\begin{aligned} \Big [\bigcup _{n=1}^Q\left( B(z_n^+,\delta ){\setminus }\overline{B(0,1)}\right) \Big ]\bigcup \Big [\bigcup _{n=1}^Q\left( B(z_n^-,\delta )\cap {B(0,1)}\right) \Big ]. \end{aligned}$$

As for sufficiently small $\epsilon >0$, $(I-{\mathcal {K}}_z^\epsilon )^{-1}$ depends analytically on z in the above domain, from Cauchy integral theorem, $u_\epsilon$ has the equivalent formulation

$$\begin{aligned} u_\epsilon (x_1+n,x_2)=\frac{1}{2\pi \mathrm {i}}\oint _{\mathcal {C}_0}w_\epsilon (z,x)z^{n-1}\,\mathrm {d}z \,. \end{aligned}$$

From the asymptotic behaviours of the poles $Z_j^\pm (\epsilon )$ explained in Lemmas 3–4, let $\epsilon$ tends to $0^+$, we get the exact formulation for the LAP solution:

$$\begin{aligned} u(x_1+n,x_2)=\frac{1}{2\pi \mathrm {i}}\oint _{\mathcal {C}_0}w(z,x)z^{n-1}\,\mathrm {d}z \,. \end{aligned}$$

(25)

We also refer to Sect. 6.2, [20] for a similar approach. Moreover, for any $n\in \mathbb {Z}$, there is a constant $C=C(n)>0$ such that

$$\begin{aligned} \Vert u\Vert _{H^1(\varOmega _n)}\le C \Vert f\Vert _{L^2(\varOmega _0)}. \end{aligned}$$

(26)

Moreover, from regularity of elliptic equations, the LAP solution $u\in H^2_{loc}(\varOmega )$.

We can also replace $\mathcal {C}_0$ by any closed curve that lies in the neighbourhood of $\mathcal {C}_0$ and enclose zero and all poles in $S_+$, but no poles in $S_-$ (see Fig. 4 (right)). The left curve is $\mathcal {C}_0$, and the right is a choice of $\mathcal {C}$. Thus

$$\begin{aligned} u(x_1+n,x_2)=\frac{1}{2\pi \mathrm {i}}\oint _{\mathcal {C}}w(z,x)z^{n-1}\,\mathrm {d}z \,. \end{aligned}$$

(27)

Now we have already formulated the LAP solution directly from cell problems (1920)–(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 $\mathcal {C}_0$ (or $\mathcal {C}$).

5.3 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_\pm$ and $P_0$ (also the corresponding S, $S_\pm$ and $S_0$). With Assumption 1, $P_0=\emptyset$ thus $P=P_+\cup P_-$ (and also $S=S_+\cup S_-$). However, since $P_+\cap P_-$ may not be empty, we can not use the old notations any more. Instead, we simply let $P=\{\beta _1,\beta _2,\dots ,\beta _Q\}$ and each $\beta _j$ ($j=1,2,\dots ,Q$) is related to $L_j$ analytic functions $\mu _{j,\ell }$ ($\ell =1,2,\dots ,L_J$) which define $L_j$ dispersion diagrams. Thus if $\mu '_{j,\ell }(\beta _j)>0(<0)$ then $\beta _j\in P_+$ ($P_-$).

Let $\psi _{j,\ell }(\beta ,\cdot )$ be the eigenfunctions corresponding to the eigenvalues $\mu _{j,\ell }(\beta )$. We choose a $\tau >0$ such that the circle $|z|=\exp (-\tau )$ enclose all poles in RS, then from [20] with a similar complex contour integral with the LAP process, for all $n\in \mathbb {Z}$:

$$\begin{aligned} u(x_1+n,x_2)={\left\{ \begin{array}{ll} \displaystyle \begin{aligned}&{} \mathrm {i}\left[ \sum _{j=1}^Q\sum _{(\mu _{j,\ell })'(\beta _j)>0}\frac{\exp (\mathrm {i}n\beta _j)\left<f,\psi _{j,\ell }(\beta _j,\cdot )\right>}{(\mu _{j,\ell })'(\beta _j)}\psi _{j,\ell }(\beta _j,\cdot )\right] \\ {} &{}\,+ \frac{1}{2\pi \mathrm {i}}\oint _{|z|=\exp (-\tau )} w(z,x)z^{n-1}\,\mathrm {d}z \,,\qquad \qquad \qquad \qquad n\ge 0; \end{aligned}\\ \\ \displaystyle \begin{aligned} &{} \mathrm {i}\left[ \sum _{j=1}^Q\sum _{(\mu _{j,\ell })'(\beta _j)<0}\frac{\exp (\mathrm {i}n\beta _j)\left<f,\psi _{j,\ell }(\beta _j,\cdot )\right>}{(\mu _{j,\ell })'(\beta _j)}\psi _{j,\ell }(\beta _j,\cdot )\right] \\ {} &{}\,+\frac{1}{2\pi \mathrm {i}}\oint _{|z|=\exp (\tau )} w(z,x)z^{n-1}\,\mathrm {d}z \,,\qquad \qquad \qquad \qquad \quad n< 0.\end{aligned} \end{array}\right. } \end{aligned}$$

(28)

6 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)\big |_{\mathcal {C}_0}$ of the quasi-periodic problem (1920)–(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 $\varOmega _0$ is covered by a family of regular and quasi-uniform meshes (see [32, 33]) $\mathcal {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 $\left\{ \phi _M^{(\ell )}:\,1\le \ell \le M\right\}$ be the piecewise linear nodal functions that equals to one at the $\ell$-th nodal and zero at other nodal points, then it is easily extended into a globally continuous and periodic function in $H^1_{loc}(\varOmega )$. M is the number of nodal points which do not lie on the right boundary. Define the discretization subspace by:

$$\begin{aligned} V_h:=\mathrm{span}\left\{ \phi _M^{(1)},\phi _M^{(2)},\dots ,\phi _M^{(M)}\right\} \subset H^1_{per}(\varOmega _0). \end{aligned}$$

Thus we are looking for the finite element solution to (14) with the expansion

$$\begin{aligned} v_z^h(x)=\sum _{\ell =1}^M c_z^\ell \phi _m^{(\ell )}(x) \end{aligned}$$

satisfies

$$\begin{aligned} a_z(v_z^h,\phi _h)=-\int _{\varOmega _0}z^{-x_1}f\overline{\phi _h}\,\mathrm {d}x \,,\quad \forall \,\phi _h\in V_h. \end{aligned}$$

By Theorem 14 in [12] the finite element approximation is estimated as follows.

Theorem 3

Suppose $f\in L^2(\varOmega _0)$ and $q\in W^{1,\infty }(\varOmega _0)$, the solution $v_z\in H^2(\varOmega _0)$ and the error between $v_z^h$ and $v_z$ is bounded by

$$\begin{aligned} \Vert v_z^h-v_z\Vert _{L^2(\varOmega _0)}\le Ch^2\Vert f\Vert _{L^2(\varOmega _0)},\quad \Vert v_z^h-v_z\Vert _{H^1_{per}(\varOmega _0)}\le Ch\Vert f\Vert _{L^2(\varOmega _0)} \end{aligned}$$

where $C>0$ is a constant independent of $z\in \mathcal {C}_0$. The estimations are also true for the function $w(z,x)=(\mathcal {F}u)(z,x)$ defined in Sect. 4.2 when $z\in \mathcal {C}_0$, i.e.,

$$\begin{aligned}&\Vert w_h(z,\cdot )-w(z,\cdot )\Vert _{L^2(\varOmega _0)}\le Ch^2\Vert f\Vert _{L^2(\varOmega _0)};\\&\Vert w_h(z,\cdot )-w(z,\cdot )\Vert _{H^1_{per}(\varOmega _0)}\le Ch\Vert f\Vert _{L^2(\varOmega _0)}. \end{aligned}$$

In both methods, $w_h(z,\cdot )$ is the numerical approximation of the quasi-periodic solution $w(z,\cdot )$.

6.1 CCI-method

In this subsection, we approximate the LAP solution numerically with the help of the curve integral (25). $w(z,\cdot )$ 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 $\varOmega _n$’s, but in this section n is taken to be 0 as an example. As the curve integral on $\mathbb {S}^1$ is standard, we only consider the case $k^2\in \sigma (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 $\mathbb {UF}=S_+^0\bigcup S_-^0$ has 2Q elements when Assumption 2 holds. Especially, from Remark 3 and Assumption 2, neither 1 nor $-1$ lies in $\mathbb {UF}$. Then the curve $\mathcal {C}_0$ is parameterized as follows.

Suppose $\mathbb {UF}=\left\{ \varepsilon _j:\, j=1,2,\dots ,2Q\right\}$ and

$$\begin{aligned} -\pi<\varepsilon _1<\varepsilon _2<\cdots<\varepsilon _{2Q-1}<\varepsilon _{2Q}<\pi . \end{aligned}$$

By assumption, for any $j=1,\dots ,2Q$, there is a sufficiently small $\delta >0$ such that the ball $B\big (\exp (\mathrm {i}\varepsilon _j),\delta \big )$ does not contain any other poles except $\exp (\mathrm {i}\varepsilon _j)$ itself. Let $\varepsilon _j^-<\varepsilon _j^+$ be two angles such that $\left\{ \exp (\mathrm {i}\varepsilon _j^+),\exp (\mathrm {i}\varepsilon _j^-)\right\} =\partial B\big (\exp (\mathrm {i}\varepsilon _j),\delta \big )\cap \mathbb {S}^1$. For sufficiently small $\delta >0$,

$$\begin{aligned} -\pi<\,\varepsilon _1^-<\varepsilon _1<\varepsilon _1^+\,<\,\varepsilon _2^-<\varepsilon _2<\varepsilon _2^+\,<\,\cdots \,<\,\varepsilon _{2Q}^-<\varepsilon _{2Q}<\varepsilon _{2Q}^+\,<\,\pi . \end{aligned}$$

(29)

Thus the curve segment of $\mathbb {S}^1\cap \mathcal {C}_0$ is parameterized as

$$\begin{aligned} \mathcal {C}_j:=\left\{ \exp (\mathrm {i}\varepsilon ):\,\varepsilon \in [\varepsilon _j^+,\varepsilon _{j+1}^-]\right\} ,\quad j=1,2,\dots ,2Q, \end{aligned}$$

where $\varepsilon _{2Q+1}^-:=\varepsilon _1^-+2\pi$.

For any $j\in \{1,2,\dots ,2Q\}$, $\mathcal {C}_0\cap \partial B\big (\exp (\mathrm {i}\varepsilon _j),\delta \big )$ lies either in the exterior or in the interior of the unit disk, depending on $\epsilon _j$. As the curve segment is part of the circle with center $\exp (\mathrm {i}\epsilon _j)$, by choosing proper angles $\theta _j^-<\theta _j^+<\theta _j^-+2\pi$. Thus the curve segment $\mathcal {C}_0\cap \partial B(\exp (\mathrm {i}\varepsilon _j),\delta )$ is parameterized as

$$\begin{aligned} \mathcal {D}_{j}=\big \{\exp (\mathrm {i}\varepsilon _j)+\delta \exp (\mathrm {i}\varepsilon ),\quad \theta \in (\theta _j^-,\theta _j^+)\big \},\quad j=1,2,\dots ,2Q. \end{aligned}$$

Now the whole curve $\mathcal {C}_0$ is parameterized piecewisely. Thus the representation of u in (25) is written as

$$\begin{aligned} u(x)&=\frac{1}{2\pi \mathrm {i}}\sum _{j=1}^{2Q}\int _{\mathcal {C}_j}w(z,x)z^{-1}\,\mathrm {d}z \,+\frac{1}{2\pi \mathrm {i}}\sum _{j=1}^{2Q}\int _{\mathcal {D}_j}w(z,x)z^{-1}\,\mathrm {d}z \,\\&=\frac{1}{2\pi }\sum _{j=1}^{2Q}\int _{\epsilon _j^+}^{\epsilon _{j+1}^-}w(\exp (\mathrm {i}t),x)\,\mathrm {d}t \,\\&\quad +\frac{1}{2\pi }\sum _{j=1}^{2Q}\int _{\theta _j^-}^{\theta _j^+}w(\exp (\mathrm {i}\alpha _j)+\delta \exp (\mathrm {i}t),x)\frac{\delta \exp (\mathrm {i}t)}{\exp (\mathrm {i}\alpha _j)+\delta \exp (\mathrm {i}t)}\,\mathrm {d}t \, \end{aligned}$$

Each integrand depends smoothly on t. Thus we only need to consider the numerical integration

$$\begin{aligned} {\mathcal {I}}(g):=\int _a^b g(t,x)\,\mathrm {d}t \,, \end{aligned}$$

where $a<b$ and g depends smoothly on $t\in [a,b]$ and $g(t,\cdot )\in H^2(\varOmega _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(\tau ):\,(a,b)\rightarrow [-\pi ,\pi ]$ be a function in $C^\infty (a,b)$ and strictly monotonically increasing function and satisfies

$$\begin{aligned} s(-\pi )=a;\,s(\pi )=b;\,s^{(\ell )}(-\pi )=s^{(\ell )}(\pi )=0,\,\ell =1,2,\dots ,N_0 \end{aligned}$$

(30)

for some positive integer $N_0$. Let $t=s(\tau )$, the integral ${\mathcal {I}}(g)$ becomes

$$\begin{aligned} {\mathcal {I}}(g)=\int _{-\pi }^{\pi } g(s(\tau ),x)s'(\tau )\,\mathrm {d}\tau \,:=\int _{-\pi }^{\pi } h(\tau ,x)\,\mathrm {d}\tau \,. \end{aligned}$$

Thus the new integrand h depends smoothly on $\tau$ and is extended to a periodic function with respect to $\tau$.

We approximate the integral ${\mathcal {I}}(g)$ by trapezoidal rule. Let $[-\pi ,\pi ]$ be divided uniformly into N subintervals, and the grid points be

$$\begin{aligned} t_j=-\pi +\frac{2\pi }{N}j,\quad j=1,\dots ,N. \end{aligned}$$

The integral is approximated by

$$\begin{aligned} {\mathcal {I}}_N(g)=\frac{2\pi }{N}\sum _{\ell =1}^N h(t_\ell ,x)=\frac{2\pi }{N}\sum _{\ell =1}^N g(s(t_\ell ),x)s'(t_\ell ). \end{aligned}$$

(31)

Then we recall the result of Theorem 9.33 in [29], and obtain the error estimation of the integral via (31):

$$\begin{aligned} \left\| {\mathcal {I}}(g)-{\mathcal {I}}_N(g)\right\| _{H^2(\varOmega _0)}\le CN^{-N_0}\Vert s\Vert _{C^{N_0}([-\pi ,\pi ];H^2(\varOmega _0))}. \end{aligned}$$

(32)

We apply the quadrature rule (31) to approximate u(x):

$$\begin{aligned} u_N(x)=&\frac{2\pi }{N}\sum _{j=1}^{2Q}\left[ \sum _{\ell =1}^N w(\exp (\mathrm {i}s_j(t_\ell )),x)s'_j(t_\ell )\right. \nonumber \\&\left. +\sum _{\ell =1}^N w(\exp (\mathrm {i}\alpha _j)+\delta \exp (\mathrm {i}s_{j+2Q}(t_\ell )),x)\frac{\delta s'_{j+2Q}(t_\ell )\exp (\mathrm {i}s_{j+2Q}(t_\ell ))}{\exp (\mathrm {i}\alpha _j)+\delta \exp (\mathrm {i}s_{j+2Q}(t_\ell ))}\right] ,\nonumber \\ \end{aligned}$$

(33)

where $s_j$ be the smooth function from $[\epsilon _j^+,\epsilon _{j+1}^-]$ to $[-\pi ,\pi ]$ and $s_{j+2Q}$ are the smooth function from $[\theta _j^-,\theta _{j}^+]$ to $[-\pi ,\pi ]$. With (32), we conclude that

$$\begin{aligned} \left\| u(x)-u_N(x)\right\| _{H^2(\varOmega _0)}\le CN^{-N_0+1/2}\Vert w\Vert _{C^{N_0}([-\pi ,\pi ];H^2(\varOmega _0))}. \end{aligned}$$

(34)

Finally, we replace w(z, x) by the finite element solution $w_h(z,x)$. Then the LAP solution is approximated by

$$\begin{aligned} \begin{aligned} u_{N,h}(x)&=\frac{1}{N}\sum _{j=1}^{2Q}\left[ \sum _{\ell =1}^N w_h(\exp (\mathrm {i}s_j(t_\ell )),x)s'_j(t_\ell )\right. \\&\quad \left. +\sum _{\ell =1}^N w_h(\exp (\mathrm {i}\epsilon _j)+\delta \exp (\mathrm {i}s_j(t_\ell )),x)\frac{\delta s'_{j+2Q}(t_\ell )\exp (\mathrm {i}s_j(t_\ell ))}{\exp (\mathrm {i}\alpha _j)+\delta \exp (\mathrm {i}s_j(t_\ell ))}\right] . \end{aligned} \end{aligned}$$

(35)

Then we conclude the error between $u_{N,h}$ and u in the domain $\varOmega _0$.

Theorem 4

Let $u_{N,h}$ be the numerical LAP solution comes from the finite element method and the integral approximation (35). Then the error is bounded by

$$\begin{aligned} \begin{aligned}&\Vert u_{N,h}-u\Vert _{L^2(\varOmega _0)}\le C\left( N^{-N_0}+h^2\right) \Vert f\Vert _{L^2(\varOmega _0)};\\ {}&\Vert u_{N,h}-u\Vert _{H^1(\varOmega _0)}\le C\left( N^{-N_0}+h\right) \Vert f\Vert _{L^2(\varOmega _0)}, \end{aligned} \end{aligned}$$

(36)

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:

$$\begin{aligned} \begin{aligned}&\!\Vert u_{N,h}-u\Vert _{L^2(\varOmega _0)}\\&\le \Vert u_{N,h}-u_N\Vert _{L^2(\varOmega _0)}+\Vert u_N-u\Vert _{L^2(\varOmega _0)}\\&\le \left\| \frac{1}{N}\sum _{j=1}^{2Q}\left[ \sum _{\ell =1}^N (w_h-w)(\exp (\mathrm {i}s_j(t_\ell )),x)s'_j(t_\ell )\right. \right. \\&\quad \left. \left. +\sum _{\ell =1}^N (w-w_h)(\exp (\mathrm {i}\epsilon _j)+\delta \exp (\mathrm {i}s_j(t_\ell )),x)\frac{\delta s'_{j+2Q}(t_\ell )\exp (\mathrm {i}s_j(t_\ell ))}{\exp (\mathrm {i}\alpha _j)+\delta \exp (\mathrm {i}s_j(t_\ell ))}\right] \right\| _{L^2(\varOmega _0)}\\&\qquad +CN^{-N_0}\Vert f\Vert _{H^1(\varOmega _0)}\\&\le \frac{C}{N}\sum _{j=1}^{2Q}\sum _{\ell =1}^N\left[ \Vert (w-w_h)(\exp (\mathrm {i}\epsilon _j)+\delta \exp (\mathrm {i}s_j(t_\ell )),\cdot )\Vert _{L^2(\varOmega _0)}\right. \\&\qquad \left. +\Vert (w-w_h)(\exp (\mathrm {i}\epsilon _j)+\delta \exp (\mathrm {i}s_j(t_\ell )),\cdot )\Vert _{L^2(\varOmega _0)}\right] +CN^{-N_0}\Vert f\Vert _{H^1(\varOmega _0)}\\&\le C\left( h^2+N^{-N_0}\right) \Vert f\Vert _{H^1(\varOmega _0)}. \end{aligned} \end{aligned}$$

The estimation of the $H^1-$norm is also obtained in the same way:

$$\begin{aligned} \Vert u_{N,h}-u\Vert _{H^1(\varOmega _0)}\le C\left( h+N^{-N_0}\right) \Vert f\Vert _{H^1(\varOmega _0)}. \end{aligned}$$

The proof is finished.

6.2 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 $\beta _j$ and corresponding eigenfunctions $\psi _{j,\ell }(\beta _j,\cdot )$. From the variational formulation (14) and by replacing $\log (z)$ by $\mathrm {i}\alpha$ in (14), the problem is to find $\alpha \in (-\pi ,\pi ]$ such that there is a non-trival $v\in H^1_{per}(\varOmega _0)$ solving

$$\begin{aligned} \int _{\varOmega _0}\left[ \nabla v\cdot \nabla \overline{\phi }+\mathrm {i}\alpha \left( v\frac{\partial \overline{\phi }}{\partial x_1}-\frac{\partial v}{\partial x_1}\overline{\phi }\right) -(k^2 q-\alpha ^2)v\overline{\phi }\right] \,\mathrm {d}x \, =0, \end{aligned}$$

for any $\phi \in H^1_{per}(\varOmega _0)$. From Theorem 2.4, [3], it is written as a generalized eigenvalue problem:

$$\begin{aligned} \left( \begin{matrix} I-K_0 &{} 0 \\ 0 &{} I \end{matrix} \right) \left( \begin{matrix} u_1 \\ u_2 \end{matrix} \right) =-\alpha \left( \begin{matrix} B &{} C^{1/2} \\ -C^{1/2} &{} 0 \end{matrix} \right) \left( \begin{matrix} u_1 \\ u_2 \end{matrix} \right) . \end{aligned}$$

(37)

By solving this problem, we obtain the values of $\beta _j$ and the corresponding eigenfunctions $\psi _{j,\ell }(\beta _j,\cdot )$ at the same time. Then the values of $\left<f,\psi _{j,\ell }(\beta _j,\cdot )\right>$ are easily calculated by numerical integral on triangular meshes.

The second step is to evaluate $(\mu _{j,\ell })'(\beta _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 $(\mu _{j,\ell })'(\beta _j)$.

First we need the definition of the energy flux:

$$\begin{aligned} \mathcal {E}(\phi ,\phi )=4k\mathfrak {I}\int _{\Gamma _j}\frac{\partial \phi }{\partial x_1}\overline{\phi }\,\mathrm {d}s \,. \end{aligned}$$

The energy flux is closely related to the directions of propagating modes. When $\phi$ is propagating to the right, $\mathcal {E}(\phi ; \phi ) > 0$; while when it is propagating to the left, $\mathcal {E}(\phi ; \phi ) < 0$. Then we recall the following equation from the proof of Theorem 3 in [2]:

$$\begin{aligned}&\mu '_{j,\ell }(\beta _j)\int _{\varOmega _0}q(x)\psi _{j,\ell }(\beta _j,x)\overline{\psi _{j,\ell }(\beta _j,x)}\,\mathrm {d}x \,\\&\quad \quad =\mathfrak {I}\int _{\Gamma _j}\frac{\partial \phi }{\partial x_1}\overline{\phi }\,\mathrm {d}s \,=\frac{E(\psi _{j,\ell }(\beta _j,\cdot ),\psi _{j,\ell }(\beta _j,\cdot ))}{2 k}. \end{aligned}$$

Thus we have the following formula to compute the value of $\mu '_{j,\ell }(\beta _j)$:

$$\begin{aligned} \mu '_{j,\ell }(\beta _j)=\frac{E(\psi _{j,\ell }(\beta _j,\cdot ),\psi _{j,\ell }(\beta _j,\cdot ))}{2k\int _{\varOmega _0}q(x)\psi _{j,\ell }(\beta _j,x)\overline{\psi _{j,\ell }(\beta _j,x)}\,\mathrm {d}x \,}. \end{aligned}$$

From the above formula the sign of $\mu '_{j,\ell }(\beta _j)$ is exactly the same as that of the energy flux $E(\psi _{j,\ell }(\beta _j,\cdot ),\psi _{j,\ell }(\beta _j,\cdot ))$.

Thus (28) is rewritten as:

$$\begin{aligned} u(x_1+n,x_2)={\left\{ \begin{array}{ll} \displaystyle \begin{aligned}&{} \mathrm {i}\left[ \sum _{j=1}^Q\sum _{\ell =1}^{L_j}\mathbb {1}_{j,\ell } \frac{\exp (\mathrm {i}n\beta _j)\left<f,\psi _{j,\ell }(\beta _j,\cdot )\right>}{(\mu _{j,\ell })'(\beta _j)}\psi _{j,\ell }(\alpha _j,\cdot )\right] \\ {} &{}\,+ \frac{1}{2\pi \mathrm {i}}\oint _{|z|=\exp (-\tau )} w(z,x)z^{n-1}\,\mathrm {d}z \,,\qquad \qquad \qquad \qquad n\ge 0; \end{aligned}\\ \\ \displaystyle \begin{aligned}&{} \mathrm {i}\left[ \sum _{j=1}^Q\sum _{\ell =1}^{L_j}\left( 1-\mathbb {1}_{j,\ell }\right) \frac{\exp (\mathrm {i}n\alpha _j)\left<f,\psi _{j,\ell }(\beta _j,\cdot )\right>}{(\mu _{j,\ell })'(\beta _j)}\psi _{j,\ell }(\beta _j,\cdot )\right] \\ {} &{}\,+\frac{1}{2\pi \mathrm {i}}\oint _{|z|=\exp (\tau )} w(z,x)z^{n-1}\,\mathrm {d}z \,,\qquad \qquad \qquad \qquad \quad n< 0;\end{aligned} \end{array}\right. } \end{aligned}$$

where $\mathbb {1}_{j,\ell }$ is an indicator function with the value 1 if $\mathcal {E}(\psi _{j,\ell },\psi _{j,\ell })>0$, and with the value 0 if $\mathcal {E}(\psi _{j,\ell },\psi _{j,\ell })<0$. We will show numerical results for the PM method also for wavenumbers that do not satisfy Assumption 2.

7 Half-guide problems

The numerical method introduced in this paper is also extended to half-guide problems. Let the half waveguide be defined by $\varOmega _+:=\bigcup _{n=1}^\infty \varOmega _n$ (see Fig. 5). Recall that

$$\begin{aligned} \varOmega _0=(-1/2,1/2]\times (0,1);\quad \varOmega _j=\varOmega _0+(j,0)^\top . \end{aligned}$$

Then we are looking for the LAP solution $u_+\in H_{loc}^1(\varOmega _+)$ such that it satisfies

$$\begin{aligned} \varDelta u_++k^2 qu_+=0\quad \text { in }\varOmega _+;\quad u_+=\phi \quad \text { on }\Gamma _1. \end{aligned}$$

(38)

From [1], for any $\phi \in H^{1/2}(\Gamma _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.

Fig. 5

Periodic waveguide

7.1 Approximation of LAP solutions

First, we define the following space:

$$\begin{aligned} \mathcal {U}:=\left\{ u(f)\big |_{\Gamma _1}:\,u(f)\in H^1_{loc}(\varOmega ) \text { is the LAP solution of}\, (1)-(2)\hbox { with }f\in L^2(\varOmega )\right\} . \end{aligned}$$

From the analysis in [20], $\overline{\mathcal {U}}=H^{1/2}(\Gamma _1)$. For any fixed $\phi \in H^{1/2}(\Gamma _1)$, given any sufficiently small $\epsilon _0>0$, there is an $\varOmega _0$-supported function f such that $\left\| u(f)\big |_{\Gamma _1}-\phi \right\| _{H^{1/2}(\Gamma _1)}\le \epsilon _0$. Then $u_+$ is approximated by:

$$\begin{aligned} \widetilde{u}_+(x_1+n,x_2)=\frac{1}{2\pi \mathrm {i}}\oint _{\mathcal {C}_0}w(z,x)z^{n-1}\,\mathrm {d}z \,,\quad x\in \varOmega _0,\,n\in \mathbb {N}\end{aligned}$$

(39)

where $w(z,\cdot )$ solves (1112)–(13) with the right hand side f, and $\mathcal {C}_0$ is defined by Definition 1.

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 $\Gamma _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 $L^2(\varOmega _0)$ to $H^{1/2}(\Gamma _1)$ defined by

$$\begin{aligned} f\mapsto \left. \left( \frac{1}{2\pi \mathrm {i}}\oint _{\mathcal {C}_0}w(z,x)z^{-1}\,\mathrm {d}z \,\right) \right| _{\Gamma _1}. \end{aligned}$$

(40)

Since $u(f)\in H^2(\varOmega _0)$, $Af=u(f)\in H^{3/2}(\Gamma _1)\subset H^{1/2}(\Gamma _1)$. The compact embedding of $H^{3/2}(\Gamma _1)$ in $H^{1/2}(\Gamma _1)$ implies that A is a compact operator with dense range. We are looking for an f such that

$$\begin{aligned} Af=\phi \quad \text { on }\Gamma _1. \end{aligned}$$

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].

7.2 Numerical method

Now we conclude the numerical scheme for half-guide scattering problems as follows.

• 1. Choose a family of basis functions for the space $L^2(\varOmega _0)$ and denote it by $\{\phi _\ell \}_{\ell =1}^\infty$. Fix a sufficiently large integer $M_0>0$ and construct the matrix

  $$\begin{aligned} \Phi :=\left( A\phi _1,A\phi _2,\dots ,A\phi _{M_0}\right) . \end{aligned}$$

  We look for a best solution c such that $\Phi c-\phi$ has the minimum value.

• 2. Fix a small regularization parameter $\gamma >0$, compute the vector

  $$\begin{aligned} c(\gamma )=\left[ \Phi ^*\Phi +\gamma I\right] ^{-1}\left( \Phi ^*\phi \right) . \end{aligned}$$

• 3. Approximate the function f by $f_\gamma ^{M_0}$ with coefficients $c(\gamma )$:

  $$\begin{aligned} f_\gamma ^{M_0}(x)=\sum _{\ell =1}^{M_0}c_\ell (\gamma )\phi _\ell (x). \end{aligned}$$

• 4. Solve the equation (1920)–(21) with the source term $f_\gamma ^{M_0}$. The LAP solution u is then obtained by either the CCI method or the PM method.

From the convergence of the Tikhonov method and numerical solution of (39), this numerical scheme converges as $M_0\rightarrow \infty$, $\gamma \rightarrow 0$, $h\rightarrow 0$ and $N\rightarrow \infty$.

8 Special wavenumbers

8.1 Formulation for special wavenumbers

In previous sections, all discussions for the CCI method are based on Assumption 2, i.e., $P_+\cap P_-=\emptyset$. As it is not a necessary condition for the limiting absorption principle, we extend the CCI method to the case $P_+\cap P_-\ne \emptyset$ 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=\{\beta _1,\dots ,\beta _Q\}$. Since $P_+\cap P_-\ne \emptyset$, assume that there is a positive integer $1\le Q'\le Q$ such that

$$\begin{aligned} P_+\cap P_-=\{\beta _1,\dots ,\beta _{Q'}\}. \end{aligned}$$

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\in \{1,2,\dots ,Q\}$, $\mu '_{j,\ell }(\beta _j)>0$ for all $\ell =1,2,\dots ,L_j$, then from [20],

$$\begin{aligned} \mathrm {i}\sum _{(\mu _{j,\ell })'(\beta _j)>0}\frac{\exp (\mathrm {i}n\beta _j)\left<f,\psi _{j,\ell }(\beta _j,\cdot )\right>}{(\mu _{j,\ell })'(\beta _j)}\psi _{j,\ell }(\beta _j,\cdot ) =\frac{1}{2\pi \mathrm {i}}\int _{\partial B(e^{\mathrm {i}\beta _j},\delta _j)}w(z,x)z^{n-1}\,\mathrm {d}z \, \end{aligned}$$

holds for all $n\in \mathbb {N}$ and sufficiently small $\delta _j>0$. In this case, we can merge this term with the curve integral in (28). Define the following domain:

$$\begin{aligned} B_{\tau _0,\delta }=B(0,\exp (-\tau _0))\bigcup \left[ \bigcup _{\beta _j\in P^+{\setminus } P^-} B(e^{\mathrm {i}\beta _j},\delta )\right] {\setminus } \left[ \bigcup _{\beta _j\in P^-{\setminus } P^+} \overline{B(e^{\mathrm {i}\beta _j},\delta )}\right] ,\end{aligned}$$

where $\delta >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\in \mathbb {N}$,

$$\begin{aligned} u(x_1+n,x_2)=&\mathrm {i}\left[ \sum _{j=1}^{Q'}\sum _{(\mu _{j,\ell })'(\beta _j)>0}\frac{\exp (\mathrm {i}n\beta _j)\left<f,\psi _{j,\ell }(\beta _j,\cdot )\right>}{(\mu _{j,\ell })'(\beta _j)}\psi _{j,\ell }(\beta _j,\cdot )\right] \\ {}&+ \frac{1}{2\pi \mathrm {i}}\oint _{\partial B_{\tau _0,\delta }} w(z,x)z^{n-1}\,\mathrm {d}z \,, \end{aligned}$$

When $j=1,2,\dots ,Q'$, since there are analytic functions $\mu _{j,\ell }$ passing through $(\beta _j,k_0^2)$ with both positive and negative values of $\mu '_{j,\ell }(\beta _j)$, the integral on a small circle with center $e^{\mathrm {i}\beta _j}$ no longer works. Note that since $\mu _{j,\ell }$ depends analytically on $\alpha$ with non-vanishing derivatives at $\beta _j$, the inverse function, denoted by $\beta _{j,\ell }$, exists and depends analytically on $k^2$ (see Chapter 16, [27]). From inverse function theorem, when $\mu '_{j,\ell }(\beta _j)>0$, $\beta _{j,\ell }(k^2)<\beta _j$ when $k^2<k_0^2$ and $\beta _{j,\ell }(k^2)>\beta _j$ when $k^2>k_0^2$. Similar for the case that $\mu '_{j,\ell }(\beta _j)<0$. Moreover, when $k^2\ne k_0^2$, $\beta _{j,\ell }$ are separated for different $\ell$ thus the integral representation works in this case.

Fig. 6

Example for $n=1$ and $k^2_0\approx 61.685$. There are four unit Floquet multipliers, corresponding to $\alpha \approx \pm 1.5708,\,\pm 0.9151$. From the definition of $S_0^\pm$, $P'=2$ and $P=3$, the green circles are $z_1\approx \exp (1.5708\mathrm {i})$ and $z_2\approx \exp (-1.5708\mathrm {i})$, the red square is $z_3^-\approx \exp (- 0.9151\mathrm {i})$ and the blue diamond is $z_3^+\approx \exp ( 0.9151\mathrm {i})$. The red curve is $\mathcal {C}_0(\tau _0)$, the balls with black boundaries are $B_j^\pm$ where $j=1,2$

Here the balls $B_j^\pm$ are defined by:

$$\begin{aligned} B_j^\pm :=B\left( \exp \left[ \mathrm {i}(\beta _j\pm \delta _0)\right] ,2\sin \left( \frac{\delta _0}{2}\right) \right) . \end{aligned}$$

(41)

Furthermore, we suppose that $\delta _0$ is sufficiently small such that $\overline{B_j^\pm }\cap S=\emptyset$. Please refer to the black circles in the right picture in Fig. 6. From above arguments, when $\mu '_{j,\ell }(\beta _j)>0$,

$$\begin{aligned} e^{\mathrm {i}\beta _{j,\ell }(k^2)}\in B_j^-\text { when }k^2<k_0^2;\quad e^{\mathrm {i}\beta _{j,\ell }(k^2)}\in B_j^+\text { when }k^2>k_0^2; \end{aligned}$$

when $\mu '_{j,\ell }(\beta _j)<0$,

$$\begin{aligned} e^{\mathrm {i}\beta _{j,\ell }(k^2)}\in B_j^+\text { when }k^2<k_0^2;\quad e^{\mathrm {i}\beta _{j,\ell }(k^2)}\in B_j^-\text { when }k^2>k_0^2; \end{aligned}$$

For $j=1,2,\dots ,Q'$, define $u_{j,k^2}^+$ by:

$$\begin{aligned} u_{j,k^2}^+(x):={\left\{ \begin{array}{ll} \frac{1}{2\pi \mathrm {i}}\int _{\partial B_j^-}w_{k^2}(z,x)z^{-1}\,\mathrm {d}z \,\quad \text { when }k^2<k_0^2;\\ \frac{1}{2\pi \mathrm {i}}\int _{\partial B_j^+}w_{k^2}(z,x)z^{-1}\,\mathrm {d}z \,\quad \text { when }k^2>k_0^2. \end{array}\right. } \end{aligned}$$

(42)

Then from a similar arguments in [20],

$$\begin{aligned} u_{j,k^2}^+=\mathrm {i}\left[ \sum _{(\mu _{j,\ell })'(\beta _j)>0}\frac{\exp (\mathrm {i}n\beta _j)\left<f,\psi _{j,\ell }(\beta _j,\cdot )\right>}{(\mu _{j,\ell })'(\beta _j)}\psi _{j,\ell }(\beta _j,\cdot )\right] \end{aligned}$$

holds uniformly when $k^2$ lies in a sufficiently small neighourhood of $k_0^2$. Thus

$$\begin{aligned} u_{k^2}(x_1+n,x_2)= \frac{1}{2\pi \mathrm {i}}\oint _{\partial B_{\tau _0,\delta }} w_{k^2}(z,x)z^{n-1}\,\mathrm {d}z \,+\sum _{j=1}^{Q'} u_{j,k^2}^+. \end{aligned}$$

(43)

As this function depends real analytically on $k^2$, the LAP solution is the limit of (43) as $k^2\rightarrow k_0^2$, where the second term is defined in (42).

8.2 Numerical method

Now we introduce the numerical method based on (43). The first term with $k^2=k_0^2$ 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_0^2$ from the method introduced in (35), denoted by $u_{0,N,h}$ with the same parameters N and h.

• 2. Fix M ($M\ge 2$) different wavenumbers $k_m^2\in K{\setminus }\{k_0^2\}$, $m=1,2,\dots ,M$, evaluate the values of $u_{k_m^2,j}^+$ for $j=1,2,\dots ,Q'$ in (43). The values are approximated by (42) with N points, and denoted by $u_{j,k_m^2,N,h}^+$.

• 3. With these M points, we approximate the value at $k_0^2$ by interpolation from $u_{j,k_m^2,N,h}^+$ where $m=1,2,\dots ,M$.

• 4. $u_{k_0^2,N,h}$ is then approximated by

  $$\begin{aligned} {u}_{k^2_0,N,M,h}(x_1+n,x_2)={u}_{0,N,h}+\sum _{j=1}^{Q'} {u}_{j,k_0^2,N,M,h}^+. \end{aligned}$$

The following error estimations can be obtained in the same way as in Theorem 4:

$$\begin{aligned} \begin{aligned} \Vert u_{0,N,h}-u_0\Vert _{H^\ell (\varOmega _0)},\,&\left\| u_{j,k^2,N,h}^+-u_{j,k^2}^+\right\| _{H^\ell (\varOmega _0)}\\ {}&\le C\left( N^{-N_0}+h^{2-\ell }\right) \Vert f\Vert _{L^2(\varOmega _0)},\quad \ell =0,1.\end{aligned} \end{aligned}$$

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_0^2$. This is due to the asymptotic behavior of the poles when $k^2\rightarrow k_0^2$: the poles approach $\exp (\mathrm {i}\alpha _j)$ which lie on the integral contours $B_j^\pm$. To obtain the error between $u_{k^2_0,N,h,M}$ and $u_{k_0^2}$, we only need to estimate the error from the interpolation. As $u_{j,k^2}^+$ depends real analytically on $k^2\in K$, there is a point $K^*\in K$ such that

$$\begin{aligned} u_{j,k^2}^+=\sum _{\ell =0}^\infty u_\ell ^j (k^2-K^*)^\ell , \end{aligned}$$

and there is a $C_0>0$ such that $\left\| u_\ell ^j\right\| _{H^1(\varOmega _0)}\le C_0^\ell$ uniformly for $\ell =0,1,2,\dots$. When we approximate $u_{j,k^2}^+$ by interpolation, from standard error estimation of interpolation,

$$\begin{aligned} \left\| {u}_{j,k_0^2,N,M,h}^+-{u}_{j,k_0^2}^+\right\| _{H^1(\varOmega _0)}\le C_0^{M} |K|^{M}, \end{aligned}$$

where $|K|=\max _{x\ne y\in K}|x-y|$. Thus we can finally obtain the error estimation of the algorithm:

$$\begin{aligned} \left\| {u}_{k^2_0,N,M,h}-{u}_{k^2_0}\right\| _{H^\ell (\varOmega _0)}\le C\left( N^{-N_0}+h^{2-\ell }+C_0^M|K|^{M}\right) \Vert f\Vert _{L^2(\varOmega _0)}, \end{aligned}$$

(44)

where $\ell =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 $\left\| u_{j,k^2,N,{M},h}^+-u_{j,k^2}^+\right\| _{H^\ell (\varOmega _0)}$ becomes larger when $k^2\rightarrow k_0^2$ due to the pole at $k_0^2$, as $C=C(|K|)$ may blow up. Moreover, $C_0$ can be a large number and it is impossible to be evaluated.

9 Numerical results

In this section, we present some numerical examples to show the efficiency of both the CCI-method and the PM-method. For all the examples, the function q is chosen as follows (see Fig. 7):

$$\begin{aligned} q(x)={\left\{ \begin{array}{ll} 1,\quad |x-a_0|>0.3;\\ 9,\quad 0.1<|x-a_0|<0.3;\\ 1+8\zeta (|x-a_0|),\quad \text {otherwise.} \end{array}\right. } \end{aligned}$$

The point $a_0=(0,0.5)$, and $\zeta (t)$ is a $C^8$-continuous function defined by

$$\begin{aligned} \zeta (t)={\left\{ \begin{array}{ll} 1,\quad t\le a;\\ 0,\quad t\ge b;\\ 1-\left[ \int _{\tau =a}^b (\tau -a)^4(\tau -b)^4\,\mathrm {d}\tau \,\right] ^{-1}\left[ \int _{\tau =a}^t (\tau -a)^4(\tau -b)^4\,\mathrm {d}\tau \,\right] ,\, a<t<b \end{array}\right. } \end{aligned}$$

with $a=0.1$, $b=0.3$.

Fig. 7

Left: surface plot of q; right: contour plot of q

Fig. 8

Left: dispersion diagram; right: z-space

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\pm 0.01,7.574\pm 0.01)$ and $(13.41\pm 0.01,15.49\pm 0.01)$. When $k^2=5$, it lies in the first stop band, thus there is no eigenvalue on $\mathbb {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 $\mathbb {S}^1$. However, when $k^2=17$, there are two points lying on the dispersion curve with $\mu _n(\alpha )=17$. This implies that $Q=1$ and $\alpha _1^+=0.9576(\pm 0.01)$ (blue diamond) and $\alpha _1^-=-0.9576(\pm 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,\cdot )$ (see Fig. 9). The parameter $N_0$ is fixed to be 6 for all the numerical examples in this section.

Fig. 9

Condition number on the complex plane

9.1 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:

$$\begin{aligned} f(x)={\left\{ \begin{array}{ll} 0,\quad |x-a_0|>0.3;\\ 3\cos (2\pi x_1)\sin (2\pi x_2),\quad 0.1<|x-a_0|<0.3;\\ 3\zeta (|x-a_0|)\cos (2\pi x_1)\sin (2\pi x_2),\quad \text {otherwise;} \end{array}\right. } \end{aligned}$$

where $a_0=(0.1,0.4)^\top$.

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 by

$$\begin{aligned} err_{N,h}=\frac{\Vert u_{N,h}-u_{exa}\Vert _{L^2(\varOmega _0)}}{\Vert u_{exa}\Vert _{L^2(\varOmega _0)}}, \end{aligned}$$

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\ge 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.

Table 1 Relative $L^2$-errors for $k^2=5$

Table 2 Relative $L^2$-errors for $k^2=17$-CCI-method

Table 3 Relative $L^2$-errors for $k^2=17$-PM-method

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.

Fig. 10

Left: $k^2=5$; Right: $k^2=17$

We also compute the energy fluxes of propagating Bloch wave solutions for this example. We approximate the integrals

$$\begin{aligned} u_1^\pm =\frac{1}{2\pi \mathrm {i}}\int _{\left| z-\exp (\mathrm {i}\alpha _1^\pm )\right| =r}w(z,x)z^{-1}\,\mathrm {d}z \, \end{aligned}$$

with the same method. Fix parameters $r=0.1$ and $N=16$, and the mesh size $h=0.005$. We obtain the values

$$\begin{aligned} \mathcal {E}(u_1^+,u_1^+)\approx 0.00115,\quad \mathcal {E}(u_1^-,u_1^-)\approx -0.00115. \end{aligned}$$

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.

9.2 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

$$\begin{aligned} \phi (x_1,x_2)=\sin x_2+0.5 x_2^2+\exp (\mathrm {i}x_1)\cos (2 x_2). \end{aligned}$$

For all the examples, we approximate the source term by

$$\begin{aligned} f_\gamma (x_1,x_2)\approx \sum _{j=-20}^{20}\sum _{\ell =0}^{10}\widehat{f}_{j,\ell }(\gamma )\exp (2\mathrm {i}\pi j x_1)\cos (\pi \ell x_2), \end{aligned}$$

where $\gamma$ is the regularization parameter. As the numerical results also depend on $\gamma$, we show different results with respect to different regularization parameters, i.e., $\gamma =10^{-2}$ and $10^{-5}$. The “exact solutions” also come from the method introduced in [5, 6], and the relative error $err_{N,h}$ is defined in the same way. For both cases, the numerical solutions are computed by the CCI-method. We show the results for $k^2=5$ with parameters $N=8,16,32$ and $h=0.02,0.01,0.005$ with $\gamma =10^{-2}$ in Table 4 and with $\gamma =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 $\gamma =10^{-2}$ in Table 6, and with $\gamma =10^{-5}$ in Table 7. We also give the contour map for the solution with $k^2=17$, $\gamma =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\ge 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 $\gamma =10^{-2}$ is larger than that to $\gamma =10^{-5}$, which is also as expected.

Fig. 11

Contour map of the solution with $k^2=17$. Left: real part of the solution; Right: imaginary part of the solution

Table 4 Relative $L^2$-errors for $k^2=5$, $\gamma =10^{-2}$

Table 5 Relative $L^2$-errors for $k^2=5$, $\gamma =10^{-5}$

Table 6 Relative $L^2$-errors for $k^2=17$, $\gamma =10^{-2}$

Table 7 Relative $L^2$-errors for $k^2=17$, $\gamma =10^{-5}$

9.3 Special wavenumbers

In this section, we show some numerical results when Assumption 2 is not satisfied, i.e., $S_+\cap S_-\ne \emptyset$. From the dispersion diagram, i.e., Fig. 12, when $k^2=9.6663(\pm 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$, $\alpha _1^+=\alpha _2^-=-1.8333(\pm 0.01)$ and $\alpha _2^+=\alpha _1^-=1.8333(\pm 0.01)$. Let the curve $\mathcal {C}_0:=\big \{z\in \mathbb {C}:\,|z|=0.8\big \}$, and $B_j^\pm$ ($j=1,2$) be defined by (41) with $\delta _0=0.1$. For the visualization of the points and curves we refer to Fig. 12.

Fig. 12

$k^2=9.6663(\pm 0.01)$. Left: dispersion diagram ; right: z-space

For the CCI-method, we choose two different interpolation strategies to carry out the numerical approximation. For the first strategy, set $M=3$, and

$$\begin{aligned} k_1^2=k^2-0.1,\,k_2^2=k^2+0.1,\,k_3^2=k^2+0.2. \end{aligned}$$

For the second one, $M=5$ and

$$\begin{aligned} k_1^2=k^2-0.2,\,k_2^2=k^2-0.1,\,k_3^2=k^2+0.1,\,k_4^2=k^2+0.2,\,k_5^2=k^2+0.3.\end{aligned}$$

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.

Table 8 Relative $L^2$-errors for $k^2=9.6663(\pm 0.01)$: CCI-method-strategy 1

Table 9 Relative $L^2$-errors for $k^2=9.6663(\pm 0.01)$: CCI-method-strategy 2

Table 10 Relative $L^2$-errors for $k^2=9.6663(\pm 0.01)$: PM-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^\pm$ and $u_2^\pm$ are evaluated as follows:

$$\begin{aligned} \begin{aligned}&\mathcal {E}(u_1^+,u_1^+)\approx 0.00235,\quad \mathcal {E}(u_1^-,u_1^-)\approx -0.00235,\\&\mathcal {E}(u_2^+,u_2^+)\approx 0.0428,\quad \mathcal {E}(u_2^-,u_2^-)\approx -0.0428. \end{aligned} \end{aligned}$$

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.

9.4 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 $\epsilon$. 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 $(-\pi ,\pi ]$. 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=\sqrt{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 NM. Thus they have to solve a system of at least $NM\times NM$. However, for our problem we only need to solve several times $M\times M$ linear system, which is much more efficient. Due to the super algebraic convergence, we do not need to solve the $M\times M$ linear system too many times. Thus our method is faster than the one introduced in that paper.

Appendix

Proof of Lemma 1

Assume that the set has a bounded accumulation point $k_0\in \mathbb {R}_+$, i.e., there is a sequence $k_n$ such that

$$\begin{aligned} P_+(k_n)\cap P_-(k_n)\ne \emptyset ,\quad \lim _{n\rightarrow \infty }k_n=k_0. \end{aligned}$$

Thus the sequence has a strictly monotones subsequence which also converges to $k_0$. Without loss of generality, we assume that the subsequence is monotonically decreasing and is still denoted by $k_n$, i.e.,

$$\begin{aligned} k_0^2<\cdots<k_n^2<k_{n-1}^2<\cdots <k_1^2. \end{aligned}$$

This implies that for any $n\in \mathbb {N}$, there is a pair $(i_n\,,j_n)\in \mathbb {N}\times \mathbb {N}$ such that

$$\begin{aligned} \exists \, \alpha _n\in (-\pi ,\pi ],\,s.t.,\,\mu _{i_n}(\alpha _n)=\mu _{j_n}(\alpha _n)=k^2_n \end{aligned}$$

that satisfies

$$\begin{aligned} \mu _{i_n}'(\alpha _n)>0,\quad \mu _{j_n}'(\alpha _n)<0. \end{aligned}$$

As $\lim _{n\rightarrow \infty }\mu _n(\alpha )=\infty$ for any $\alpha \in (-\pi ,\pi ]$, there should be a subsequence of pairs $(i_n,\,j_n)$ such that $i_n=i_0$ and $j_n=j_0$ where $i_0$ and $j_0$ are two constant positive integers. Still denote the subsequence of $k_n^2$ by $k_n^2$, then for any $n\in \mathbb {N}$, there is an $\alpha _n\in (-\pi ,\pi ]$ such that

$$\begin{aligned} \mu _{i_0}(\alpha _n)=\mu _{j_0}(\alpha _n)=k^2_n,\quad \mu '_{i_0}(\alpha _n)>0,\,\mu '_{j_0}(\alpha _n)<0. \end{aligned}$$

Define the function

$$\begin{aligned} \mu (\alpha ):=\mu _{i_0}(\alpha )-\mu _{j_0}(\alpha ), \end{aligned}$$

then there is a sequence $\alpha _n\in (-\pi ,\pi ]$ such that

$$\begin{aligned} \mu (\alpha _n)=0,\quad \forall n\in \mathbb {N}. \end{aligned}$$

As $\mu _{i_0}$ and $\mu _{j_0}$ are both analytic functions, $\mu$ is analytic as well. Thus either $\mu$ is a constant function equals to 0, or $\alpha _n=\alpha _0$ except for a finite number of n’s.

For the first case, $\mu '(\alpha )=0$ for any $\alpha _n$, which contradicts with

$$\begin{aligned} \mu '(\alpha _n)=\mu '_{i_0}(\alpha _n)-\mu '_{j_0}(\alpha _n)>0. \end{aligned}$$

For the second case, suppose there is an $N>>1$ such that $\alpha _n=\alpha _0$ for any $n\ge N$, then $\mu _{i_0}(\alpha _n)=\mu _{j_0}(\alpha _n)=k^2_n$ implies that $k^2_n=k^2_0$ for any $n\ge N$. This contradicts with the monotone decreasing property. Thus $k_0^2$ can not be an accumulation point, the proof is finished.