# A highly efficient and accurate numerical method for the electromagnetic
scattering problem with rectangular cavities

Peijun Li  Department of Mathematics, Purdue University, West Lafayette,
Indiana 47907, USA<EMAIL_ADDRESS>and  Xiaokai Yuan  School of
Mathematics, Jilin University, Changchun 130012, Jilin, China
<EMAIL_ADDRESS>
###### Abstract.

This paper presents a robust numerical solution to the electromagnetic
scattering problem involving multiple multi-layered cavities in both
transverse magnetic and electric polarizations. A transparent boundary
condition is introduced at the open aperture of the cavity to transform the
problem from an unbounded domain into that of bounded cavities. By employing
Fourier series expansion of the solution, we reduce the original boundary
value problem to a two-point boundary value problem, represented as an
ordinary differential equation for the Fourier coefficients. The analytical
derivation of the connection formula for the solution enables us to construct
a small-scale system that includes solely the Fourier coefficients on the
aperture, streamlining the solving process. Furthermore, we propose accurate
numerical quadrature formulas designed to efficiently handle the weakly
singular integrals that arise in the transparent boundary conditions. To
demonstrate the effectiveness and versatility of our proposed method, a series
of numerical experiments are conducted.

###### Key words and phrases:

Electromagnetic scattering, the Helmholtz equation, cavity scattering problem,
transparent boundary condition, hypersingular and weakly singular integrals,
numerical quadratures.

###### 2010 Mathematics Subject Classification:

78A25, 78M25, 35Q60, 65N80

The first author is supported in part by the NSF grant DMS-2208256. The second
author is partially supported by the NSFC grants 12201245 and 12171017.

##  1\. Introduction

Electromagnetic cavity scattering problems find significant applications in
various fields. For instance, in radar and remote sensing, a comprehensive
understanding of cavity scattering is crucial for radar systems utilized in
target detection, identification, and tracking [4]. In wireless communication
systems, cavities can arise from surrounding structures or obstacles, and
analyzing cavity scattering aids in predicting signal propagation,
interference, and overall system performance [17]. Furthermore, cavity
scattering plays a significant role in the behavior of metamaterials and
photonic devices, influencing their unique electromagnetic properties and
guiding applications in areas like superlensing, cloaking, and wave
manipulation [14, 26]. Rigorous analysis and accurate computation of cavity
scattering are vital for technological progress and effectively addressing
real-world challenges.

The crucial industrial and military applications of cavity scattering problems
have made them a focal point of interest for both engineering and mathematical
communities. In the engineering community, researchers have initiated the
investigation of electromagnetic scattering by cavities filled with penetrable
materials [15, 24]. The well-posedness of cavity scattering problems has been
rigorously analyzed using integral equation methods or variational approaches,
with detailed studies available in [1, 2, 21, 27]. For rectangular shaped
cavity scattering, a refined stability estimate with an explicit wavenumber
dependence has been derived [6, 7]. Lately, the exploration of subwavelength
enhancement has emerged as an important theme in mathematical research [22,
23]. In [3, 8, 14], the field enhancement is explored for both single and
double rectangle cavities under various boundary conditions. These studies aim
to provide a more profound mathematical understanding of the subwavelength
enhancement phenomenon. We refer to [20] for a survey of recent developments
in mathematical modeling and analysis of cavity scattering problems.

Numerical methods have been extensively studied for the solution of
electromagnetic cavity scattering problems. Since the problem is formulated in
an unbounded domain, it is essential to employ an appropriate artificial
boundary condition to reformulate it into a bounded domain. Several approaches
for introducing artificial boundary conditions include using the Green’s
function method on the aperture [13, 19], employing a perfectly matched layer
in the exterior of the cavity [11], applying the Fourier transform on the
ground [2], constructing transparent boundary conditions by utilizing the
Fourier series expansion on the semi-circle over the cavity [31], or adopting
the method of boundary integral equations [18]. In [28, 29, 5], numerical
quadrature formulas were developed to discretize the hypersingular integrals
in the method of Green’s function. The wave field inside the cavity was
approximated using a second-order finite difference scheme. The approach was
further extended to a fourth-order scheme in [32]. In [11] and [31], an
adaptive finite element method was developed, combining perfectly matched
layer and transparent boundary conditions to handle the possible singularity
of the solution, respectively. For problems involving scattering by multiple
cavities, the Gauss–Seidel technique or a preconditioned iterative method was
employed to accelerate computation [12, 21, 30, 33]. Notably, in the methods
mentioned above, the model equation is discretized within the entire cavity.
When the cavity possesses a rectangular shape, the field inside can be
approximated using its Fourier series expansion, which reduces the scattering
problem to one-dimensional ordinary differential equations for the Fourier
coefficients. In [9], these ordinary differential equations are discretized
using a second-order finite difference scheme. The Fourier coefficients inside
the cavity are expressed in terms of the Fourier coefficients on the aperture
through Gaussian elimination. Finally, a linear system on the aperture is
obtained by applying a finite difference scheme, which provides an effective
approach to solve the scattering problem for rectangular cavities.

This paper presents a highly efficient and accurate numerical method for
solving the electromagnetic scattering problem involving multiple multi-
layered rectangular cavities. We assume that these cavities are embedded in
the ground, with their apertures aligned with the ground, and their interiors
filled possibly with multiple layered media. To tackle this problem, we follow
the approach presented in [9], where the field inside the cavity is expanded
using its Fourier series, and the Helmholtz equation is reduced to ordinary
differential equations. By studying the equations and the transmission
conditions across each layered medium, we establish a connection formula that
links the Fourier coefficients of the solution in each layer to the Fourier
coefficients of the solution on the aperture. This connection formula, along
with the transparent boundary condition on the aperture, leads to a small
$N$-by-$N$ linear system, where $N$ represents the Fourier truncation number.
To efficiently generate the linear system, we design an alternative
transparent boundary condition, which only involves weakly singular integrals.
However, due to the singularity of the Hankel function, direct application of
high-order quadrature formulas is not feasible. To address this issue, we
utilize the power series of Bessel functions to deduce a recursive formula,
enhancing the regularity of the integrand function. This enables us to adopt
high-order Gaussian quadratures. Once the system is solved, the field inside
the cavity can be obtained immediately using the connection formula.

The proposed method offers the advantage of significantly reducing memory and
computational costs. As we only need to solve the system on the aperture and
store the connection formula of the Fourier coefficients, the required
computational resources are dramatically reduced. A series of numerical
experiments is conducted to demonstrate the efficiency and versatility of our
proposed method. It proves to be efficient and accurate in handling the cavity
scattering problem in both transverse magnetic (TM) and transverse electric
(TE) polarizations.

The structure of this paper is as follows. Section 2 focuses on the model
formulation, wherein the two fundamental polarizations are introduced. In
Sections 3 and 4, we derive the connection formula for scattering in TM and TE
polarization, respectively. This includes scenarios with a single empty
cavity, a single multi-layered cavity, and multiple cavities filled with
multi-layered media. Section 5 is dedicated to deriving an alternative
artificial boundary condition and proposing quadrature formulas for the
involved weakly singular integrals. In Section 6, we provide numerical
examples to demonstrate the features of the proposed method. The paper
concludes with overall reflections and avenues for future research in Section
7.

##  2\. Problem formulation

We examine the electromagnetic scattering by rectangular cavities situated
within an unbounded ground plane. Given the time dependence of the
electromagnetic field as $e^{-{\rm i}\omega t}$, with $\omega>0$ denoting the
angular frequency, the wave propagation obeys the time-harmonic Maxwell’s
equations:

$\nabla\times\boldsymbol{E}={\rm
i}\omega\mu\boldsymbol{H},\quad\nabla\times\boldsymbol{H}=-{\rm
i}\omega\epsilon\boldsymbol{E}+\sigma\boldsymbol{E},$ (2.1)

where $\boldsymbol{E}$ and $\boldsymbol{H}$ represent the electric field and
the magnetic field, respectively, $\mu$ is the magnetic permeability,
$\epsilon$ denotes the electric permittivity, and $\sigma$ stands for the
electrical conductivity. We assume that the medium is non-magnetic, implying a
constant magnetic permeability $\mu$ throughout. However, the electric
permittivity $\epsilon$ and electrical conductivity $\sigma$ are allowed to
vary as spatial functions.

In this work, we focus on the electromagnetic scattering problem in TM and TE
polarizations. Let $\Omega\subset\mathbb{R}^{2}$ represent the cross-section
of the $z$-invariant cavity. Its boundary is denoted by
$\partial\Omega=\Gamma_{c}\cup\Gamma_{g}$, where $\Gamma_{c}$ represents the
boundary of the cavity, including the vertical walls and the horizontal
bottom, while $\Gamma_{g}$ denotes the infinite ground plane. The aperture of
the cavity, aligned with $\Gamma_{g}$, is denoted by $\Gamma$. The cavity may
be filled vertically with a layered inhomogeneous medium. The problem geometry
is illustrated in Figure 1.

Figure 1. Geometry of multiple multi-layered rectangular cavities.

In TM polarization, the incident and total electric fields are perpendicular
to the magnetic field and take the form
$\boldsymbol{E}^{i}=(0,0,u^{i}),\boldsymbol{E}=(0,0,u)$. It can be verified
from (2.1) that $u$ satisfies the two-dimensional Helmholtz equation

$\Delta u+\kappa^{2}u=0\quad\text{in}~{}\mathbb{R}^{2}_{+}\cup\Omega,$ (2.2)

where the wavenumber $\kappa=(\omega^{2}\epsilon\mu+{\rm
i}\omega\mu\sigma)^{1/2}$ with $\Im\kappa\geq 0$. For TE polarization, the
magnetic field can be represented as $\boldsymbol{H}=(0,0,u)$. Similarly, we
can show from (2.1) that $u$ satisfies the two-dimensional generalized
Helmholtz equation

$\nabla\cdot(\kappa^{-2}\nabla
u)+u=0\quad\text{in}~{}\mathbb{R}^{2}_{+}\cup\Omega.$ (2.3)

By considering that both the ground plane and the cavity boundary exhibit
perfect electric conductivity (PEC), we have

$\nu\times\boldsymbol{E}=0\quad\text{on}~{}\Gamma_{g}\cup\Gamma_{c},$ (2.4)

where $\nu$ denotes the unit normal vector to the surfaces $\Gamma_{g}$ and
$\Gamma_{c}$. Under TM polarization, the boundary condition (2.4) simplifies
to

$u=0\quad\text{on}~{}\Gamma_{g}\cup\Gamma_{c}.$ (2.5)

In TE polarization, the boundary condition (2.4) is equivalent to

$\partial_{\nu}u=0\quad\text{on}~{}\Gamma_{g}\cup\Gamma_{c}.$ (2.6)

When the medium in the upper half space is homogeneous and isotropic, it can
be characterized by a constant wavenumber denoted as $\kappa_{0}$. Let the
cavity be illuminated from above by a time-harmonic plane wave

$u^{i}(x,y)=e^{{\rm i}(\alpha x-\beta
y)},\quad\alpha=\kappa_{0}\sin\theta,\quad\beta=\kappa_{0}\cos\theta,$

where $\theta\in(-\pi/2,\pi/2)$ is the incident angle. The scattered field
$u^{s}$ in TM polarization can be described as $u^{s}=u-u^{i}+u^{r}$, and in
TE polarization, it is given by $u^{s}=u-u^{i}-u^{r}$, where $u^{r}=e^{{\rm
i}(\alpha x+\beta y)}$ is referred to as the reflection field. In both cases,
the scattered field satisfies the Sommerfeld radiation condition:

$\lim\limits_{r\rightarrow\infty}\sqrt{r}\left(\partial_{r}u^{s}-{\rm
i}\kappa_{0}u^{s}\right)=0,\quad r=\sqrt{x^{2}+y^{2}}.$ (2.7)

Based on (2.7), the transparent boundary conditions (TBC) can be formulated on
the aperture [32]. In TM polarization, the TBC can be expressed as

$\partial_{y}u=I_{\rm TM}(u)+f,\quad f=-2{\rm i}\beta e^{{\rm i}\alpha x},$
(2.8)

where $I_{\rm TM}$ represents the Dirichlet-to-Neumann (DtN) operator, while
the associated hypersingular integral is defined in the sense of the Hadamard
finite-part, and its expression is given by

$I_{\rm TM}(u)(x)=\frac{{\rm
i}\kappa_{0}}{2}\int_{\Gamma}\frac{1}{|x-x^{\prime}|}H_{1}^{(1)}(\kappa_{0}|x-x^{\prime}|)u(x^{\prime},0){\rm
d}x^{\prime}.$ (2.9)

Here, $H_{1}^{(1)}$ represents the Hankel function of the first kind with
order one. In TE polarization, the TBC is described by

$u=I_{\rm TE}(u)+g,\quad g=2e^{{\rm i}\alpha x},$ (2.10)

where the Neumann-to-Dirichlet (NtD) operator $I_{\rm TE}$ involves a weakly
singular integral and is defined by

$I_{\rm TE}(u)(x)=-\frac{{\rm
i}}{2}\int_{\Gamma}H_{0}^{(1)}\left(\kappa_{0}|x-x^{\prime}|\right)\partial_{y^{\prime}}u(x^{\prime},0){\rm
d}x^{\prime}.$ (2.11)

Here, $H_{0}^{(1)}$ is the Hankel function of the first kind with order zero.

It is important to note that the quantities $u$ and $\partial_{y}u$ appearing
in (2.8) and (2.10) are to be interpreted as being evaluated from the upper
side of $\Gamma$, indicated more precisely as $u(x,0^{+})$ and
$\partial_{y}u(x,0^{+})$, respectively.

##  3\. TM polarization

In this section, we focus on the boundary value problem (2.2), (2.5), (2.8) in
TM polarization. To present our findings clearly, we start by considering the
scattering phenomena from a single empty cavity. Subsequently, we extend our
investigation to the case of a single multi-layered cavity. Finally, we
explore the general scenario of multiple cavities, with each cavity being
filled possibly with a multi-layered medium.

###  3.1. A single empty cavity

Assuming that the empty cavity is filled with a homogeneous medium
characterized by the wavenumber $\kappa=\kappa_{0}$. Let
$\Gamma=[a,b]\times\\{0\\}$ denote the aperture of the cavity, where $a$ and
$b$ are the coordinates along the $x$-axis, and the width and depth of the
cavity are represented by $w=b-a$ and $h$, respectively.

Because of the boundary condition (2.5), the total field inside the cavity can
be approximated using sine functions, i.e., it admits the Fourier series
expansion

$u(x,y)=\sum\limits_{n=1}^{N}u^{(n)}(y)\sin\frac{n\pi(x-a)}{w},$ (3.1)

where $N$ is a positive integer that controls the accuracy of the numerical
solution. Substituting (3.1) into (2.2), we obtain the second order ordinary
differential equations for the Fourier coefficients

$\displaystyle u^{(n)^{\prime\prime}}(y)+\beta_{n}^{2}u^{(n)}(y)=0,\quad
y\in(-h,0),$

where
$\beta_{n}=\big{(}\kappa^{2}-\left(\frac{n\pi}{w}\right)^{2}\big{)}^{1/2},\Im\beta_{n}\geq
0$ for $n=1,\dots,N.$

Let $u^{(n)},n=1,\dots,N$, be the Fourier coefficients of the solution on
$\Gamma$. Since $u$ vanishes on the bottom of the cavity, we have from (3.1)
that $u^{(n)}(-h)=0,n=1,\dots,N$. Therefore, the $n$-th order Fourier
coefficient satisfies a two-point boundary value problem

$\left\\{\begin{aligned}
&u^{(n)^{\prime\prime}}(y)+\beta_{n}^{2}u^{(n)}(y)=0,\quad-h<y<0,\\\
&u^{(n)}(0)=u^{(n)},\quad u^{(n)}(-h)=0.\end{aligned}\right.$ (3.2)

A straightforward calculation shows that the solution to (3.2) is given by

$u^{(n)}(y)=\left\\{\begin{aligned} &\frac{1}{1-e^{2{\rm
i}\beta_{n}h}}\big{(}e^{-{\rm i}\beta_{n}y}-e^{2{\rm i}\beta_{n}h}e^{{\rm
i}\beta_{n}y}\big{)}u^{(n)}&\quad\text{if }~{}\beta_{n}\neq 0,\\\
&\Big{(}1+\frac{1}{h}y\Big{)}u^{(n)}&\quad\text{if
}~{}\beta_{n}=0.\end{aligned}\right.$ (3.3)

Let $n_{0}$ denote the term that could lead to $\beta_{n}=0$. Combining (3.1)
and (3.3) yields that the expression for $u(x,y)$ is given as follows:

$\displaystyle u(x,y)$ $\displaystyle=\sum\limits_{n=1,n\neq
n_{0}}^{N}\frac{1}{1-e^{2{\rm i}\beta_{n}h}}\big{(}e^{-{\rm
i}\beta_{n}y}-e^{2{\rm i}\beta_{n}h}e^{{\rm
i}\beta_{n}y}\big{)}\sin\frac{n\pi(x-a)}{w}u^{(n)}$
$\displaystyle\quad+\Big{(}1+\frac{1}{h}y\Big{)}\sin\frac{n_{0}\pi(x-a)}{w}u^{(n_{0})}.$
(3.4)

Taking the normal derivative of (3.1) from the lower side of $\Gamma$, we
obtain

$\displaystyle\partial_{y}u(x,0^{-})=\sum\limits_{n=1}^{N}s^{(n)}u^{(n)}\sin\frac{n\pi(x-a)}{w},$
(3.5)

where $s^{(n_{0})}=1/h$ and $s^{(n)}=-{\rm i}\beta_{n}(1+e^{2{\rm
i}\beta_{n}h})/(1-e^{2{\rm i}\beta_{n}h})$ for $n\neq n_{0}$. Using the
transmission conditions on $\Gamma$:
$u(x,0^{+})=u(x,0^{-}),\partial_{y}u(x,0^{+})=\partial_{y}u(x,0^{-}),$ and
substituting (3.5) into the TBC (2.8) leads to

$\displaystyle\sum\limits_{n=1}^{N}s^{(n)}u^{(n)}\sin\frac{n\pi(x-a)}{w}=\sum\limits_{n=1}^{N}u^{(n)}I_{\rm
TM}\Big{(}\sin\frac{n\pi(x-a)}{w}\Big{)}+f.$ (3.6)

Multiplying both sides of (3.6) by $\sin\frac{m\pi(x-a)}{w}$, where
$m=1,...,N$, and integrating over $\Gamma$, we obtain

$D_{\rm TM}U=M_{\rm TM}U+F,$ (3.7)

where $U=[u^{(1)},u^{(2)},\dots,u^{(N)}]^{\top}$, the $m$-th entry of the
vector $F$ is defined by

$F(m)=-2{\rm i}\beta\int_{0}^{w}e^{{\rm i}\alpha(x+a)}\sin\frac{m\pi x}{w}{\rm
d}x,$

$D_{\rm TM}$ is a diagonal matrix with its diagonal entries $D_{\rm
TM}(m,m)=ws^{(m)}/2,$ and the elements of the matrix $M_{\rm TM}$, with
dimensions $N\times N$, are expressed as

$M_{\rm TM}(m,n)=\int_{\Gamma}\sin\frac{m\pi(x-a)}{w}I_{\rm
TM}\Big{(}\sin\frac{n\pi(x^{\prime}-a)}{w}\Big{)}{\rm d}x.$ (3.8)

Once the linear system (3.7) is solved, the solution inside the cavity can be
explicitly computed using the expansion (3.1). However, the matrix entries
(3.8) contain hypersingular integrals, which make it challenging to devise a
high-order quadrature. In Section 5, we will reformulate (2.8) to an
alternative form that incorporates weakly singular integrals. Additionally, we
will introduce high-order quadrature formulas to effectively tackle this
concern.

###  3.2. A single multi-layered cavity

Consider a cavity situated at $[a,b]$ and filled with an $L$-layered
inhomogeneous medium, where $L\geq 2$ is an integer. Let $w=b-a$ and $h$ be
the width and depth of the cavity, respetively. Denote the sequence of values
$-h=y_{L}<y_{L-1}<\cdots<y_{1}<y_{0}=0$. Consider the interfaces between
different media represented by $\Gamma_{l}=[a,b]\times\left\\{y_{l}\right\\}$,
where $l=1,...,L-1$. The aperture and bottom of the cavity are denoted by
$\Gamma_{0}$ and $\Gamma_{L}$, respectively. For each $l=1,...,L$, let
$\Omega_{l}$ be the region between $\Gamma_{l-1}$ and $\Gamma_{l}$. It is
assumed that the medium in each layer is homogeneous and characterized by a
constant wavenumber $\kappa_{l}$.

Let $u_{l}(x,y)$ denote the total field within the $l$-th layer, which
satisfies

$\Delta u_{l}(x,y)+\kappa_{l}^{2}u_{l}(x,y)=0\quad\text{in }~{}\Omega_{l},$
(3.9)

along with the transmission conditions on the interface $\Gamma_{l}$:

$u_{l}(x,y_{l})=u_{l+1}(x,y_{l}),\quad\partial_{y}u_{l}(x,y_{l})=\partial_{y}u_{l+1}(x,y_{l}).$
(3.10)

The solution of (3.9) can be approximated by

$u_{l}(x,y)=\sum\limits_{n=1}^{N}u_{l}^{(n)}(y)\sin\frac{n\pi(x-a)}{w}.$
(3.11)

Consider $u_{l}^{(n)},n=1,\dots,N$, for $l=0,1,\dots,L$, as the Fourier
coefficients of the solution on $\Gamma_{l}$. Note from (2.5) that
$u_{L}^{(n)}=0,n=1,\dots,N$. For $l=1,\dots,L$, substituting (3.11) into (3.9)
reveals that the Fourier coefficients satisfy

$\left\\{\begin{aligned}
&u_{l}^{(n)^{\prime\prime}}(y)+\big{(}\beta_{l}^{(n)}\big{)}^{2}u_{l}^{(n)}(y)=0,\quad
y_{l}<y<y_{l-1},\\\ &u_{l}^{(n)}(y_{l})=u_{l}^{(n)},\quad
u_{l}^{(n)}(y_{l-1})=u_{l-1}^{(n)},\end{aligned}\right.$ (3.12)

where
$\beta_{l}^{(n)}=\big{(}\kappa_{l}^{2}-\left(\frac{n\pi}{w}\right)^{2}\big{)}^{1/2}$
with $\Im\beta_{l}^{(n)}\geq 0.$

A straightforward computation shows that the solution of (3.12) is given by

$u_{l}^{(n)}(y)=\left\\{\begin{aligned}
&\frac{1}{\zeta_{l}}\Big{[}\big{(}e^{{\rm
i}\beta_{l}^{(n)}\left(y-y_{l-1}\right)}-e^{-{\rm
i}\beta_{l}^{(n)}\left(y-y_{l-1}\right)}\big{)}u_{l}^{(n)}\\\ &\hskip
28.45274pt-\big{(}e^{{\rm i}\beta_{l}^{(n)}\left(y-y_{l}\right)}-e^{-{\rm
i}\beta_{l}^{(n)}\left(y-y_{l}\right)}\big{)}u_{l-1}^{(n)}\Big{]}&\quad\text{if
}~{}\beta_{l}^{(n)}\neq 0,\\\
&\frac{1}{h_{l}}\Big{[}\big{(}u_{l}^{(n)}-u_{l-1}^{(n)}\big{)}y+u^{(n)}_{l-1}y_{l}-u^{(n)}_{l}y_{l-1}\Big{]}&\quad\text{if
}~{}\beta_{l}^{(n)}=0,\end{aligned}\right.$ (3.13)

where $h_{l}=y_{l}-y_{l-1},\zeta_{l}=e^{{\rm i}\beta_{l}^{(n)}h_{l}}-e^{-{\rm
i}\beta_{l}^{(n)}h_{l}}$.

Taking the derivatives of (3.13) from the upper and lower sides of
$\Gamma_{l}$, and using (3.10)–(3.11), we deduce a symmetric tri-diagonal
linear system for cases where $L\geq 3$, which is referred to as the
connection formula and is presented as

$D_{\rm TM}^{(n)}\boldsymbol{u}^{(n)}=\boldsymbol{b}_{\rm TM}^{(n)},$ (3.14)

where
$\boldsymbol{u}^{(n)}=[u_{1}^{(n)},u_{2}^{(n)},\dots,u_{L-1}^{(n)}]^{\top}$,
$\boldsymbol{b}_{\rm TM}^{(n)}=[-a_{1}^{(n)}u_{0}^{(n)},0,\cdots,0]^{\top}$,
and the diagonal and sub-diagonal entries of the matrix $D_{\rm TM}^{(n)}$ are
given by

$D_{\rm TM}^{(n)}(l,l)=b_{l}^{(n)}+b_{l+1}^{(n)},\quad D_{\rm
TM}^{(n)}(l,l+1)=a_{l+1}^{(n)}.$

Here

$\displaystyle a_{l}^{(n)}$ $\displaystyle=\left\\{\begin{aligned}
&-\frac{1}{h_{l}}\quad&\text{if }~{}\beta_{l}^{(n)}=0,\\\ &-2\frac{{\rm
i}\beta_{l}^{(n)}}{\zeta_{l}}\quad&\text{if }~{}\beta_{l}^{(n)}\neq
0,\end{aligned}\right.$ (3.15) $\displaystyle b_{l}^{(n)}$
$\displaystyle=\left\\{\begin{aligned} &\frac{1}{h_{l}}\quad&\text{if
}~{}\beta_{l}^{(n)}=0,\\\ &\frac{{\rm
i}\beta_{l}^{(n)}}{\zeta_{l}}\big{(}e^{{\rm i}\beta_{l}^{(n)}h_{l}}+e^{-{\rm
i}\beta_{l}^{(n)}h_{l}}\big{)}\quad&\text{if }~{}\beta_{l}^{(n)}\neq
0.\end{aligned}\right.$ (3.16)

When $L=2$, the connection formula (3.14) is simply reduced to

$\big{(}b_{1}^{(n)}+b_{2}^{(n)}\big{)}u_{1}^{(n)}=-a_{1}^{(n)}u_{0}^{(n)}.$
(3.17)

By solving (3.14) or (3.17), we can retrieve the Fourier coefficients of the
solution on the interfaces ${\boldsymbol{u}}^{(n)}$ in terms of the Fourier
coefficients of the solution on the aperture $u_{0}^{(n)}$. Since
$u_{0}^{(n)}$ is unknown, in practice, we can first solve (3.14) or (3.17)
with the right-hand side $\hat{\boldsymbol{b}}^{(n)}=[1,0,..,0]^{\top}$, and
denote the corresponding solution as
$\hat{\boldsymbol{u}}^{(n)}=[\hat{u}_{1}^{(n)},\hat{u}_{2}^{(n)},\dots,\hat{u}_{L-1}^{(n)}]^{\top}$.
Then, the solution ${\boldsymbol{u}}^{(n)}$ can be obtained by multiplying the
factor $-a_{1}^{(n)}u_{0}^{(n)}$ to $\hat{\boldsymbol{u}}^{(n)}$.

Next, we discuss the approach to determining $u_{0}^{(n)},n=1,\dots,N$. For
the case when $l=1$, we take the normal derivative of (3.11) on $\Gamma$ from
below, leading to

$\partial_{y}u_{1}(x,0^{-})=\sum\limits_{n=1}^{N}{\hat{s}}^{(n)}u_{0}^{(n)}\sin\frac{n\pi(x-a)}{w},$
(3.18)

where ${\hat{s}}^{(n)}=-b_{1}^{(n)}+(a_{1}^{(n)})^{2}\hat{u}_{1}^{(n)}$. We
mention that the normal derivative can also be described in the form of (3.18)
with ${\hat{s}}^{(n)}=-b_{1}^{(n)}$ for the case where $L=1$.

Utilizing the transmission conditions on $\Gamma$:
$u(x,0^{+})=u_{1}(x,0^{-}),\partial_{y}u(x,0^{+})=\partial_{y}u_{1}(x,0^{-})$,
and replacing (3.5) with (3.18) and following the same discussion as that in
Subsection 3.1, we arrive at an analogous formulation of the system depicted
in (3.7): $D_{\rm TM}U=M_{\rm TM}U+F$, except that the diagonal matrix $D_{\rm
TM}$ is now defined as $D_{\rm TM}(m,m)=w{\hat{s}}^{(m)}/2,m=1,\dots,N.$.
After the system is solved, the solution in each layer can be computed
explicitly by (3.11).

###  3.3. Multiple multi-layered cavities

Assume that the ground contains a total of $K$ cavities, and we denote the
aperture of the $k$-th cavity as
$\Gamma_{k}=[a_{k},b_{k}]\times\left\\{0\right\\}$, with a width of
$w_{k}=b_{k}-a_{k}$. Inside the $k$-th cavity, it is filled with an
inhomogeneous medium consisting of $L_{k}$ layers. The interfaces between
different layers are denoted by
$\Gamma_{k,l}=[a_{k},b_{k}]\times\left\\{y_{l}^{k}\right\\}$, where
$l=1,..,L_{k}-1$. The aperture and the bottom of the $k$-th cavity are denoted
by $\Gamma_{k,0}$ and $\Gamma_{k,L_{k}}$, respectively. In the $l$-th layer of
the $k$-th cavity, we assume that the medium is homogeneous and characterized
by a constant wavenumber $\kappa_{k,l}$.

Denote by $u(x,y;k)$ the total field restricted in the $k$-th cavity. The DtN
operator (2.9) can be written as

$I_{\rm TM}(u)(x)=\frac{{\rm
i}\kappa_{0}}{2}\sum\limits_{k=1}^{K}\int_{\Gamma_{k}}\frac{1}{|x-x^{\prime}|}H_{1}^{(1)}(\kappa_{0}|x-x^{\prime}|)u(x^{\prime},0;k){\rm
d}x^{\prime}.$

The solution in the $l$-th layer of the $k$-th cavity can be approximated by

$u_{l}(x,y;k)=\sum\limits_{n=1}^{N}u_{l}^{(n)}(y;k)\sin\frac{n\pi(x-a_{k})}{w_{k}},\quad
y\in(y_{l}^{k},y_{l-1}^{k}).$ (3.19)

It can be verified that the Fourier coefficients of (3.19) satisfy the same
system as (3.12), and the solution is given by

$u_{l}^{(n)}(y;k)=\left\\{\begin{aligned}
&\frac{1}{\zeta^{k}_{l}}\Big{[}\big{(}e^{{\rm
i}\beta_{l,k}^{(n)}(y-y^{k}_{l-1})}-e^{-{\rm
i}\beta_{l,k}^{(n)}(y-y^{k}_{l-1})}\big{)}u_{l,k}^{(n)}\\\ &\hskip
28.45274pt-\big{(}e^{{\rm i}\beta_{l,k}^{(n)}(y-y^{k}_{l})}-e^{-{\rm
i}\beta_{l,k}^{(n)}(y-y^{k}_{l})}\big{)}u_{l-1,k}^{(n)}\Big{]}&\quad\text{if
}~{}\beta_{l,k}^{(n)}\neq 0,\\\
&\frac{1}{h_{l}^{k}}\Big{[}\big{(}u_{l,k}^{(n)}-u_{l-1,k}^{(n)}\big{)}y+y^{k}_{l}u^{(n)}_{l-1,k}-y^{k}_{l-1}u^{(n)}_{l,k}\Big{]}&\quad\text{if
}~{}\beta_{l,k}^{(n)}=0,\end{aligned}\right.$ (3.20)

where $u^{(n)}_{l,k},n=1,\dots,N,l=0,1,\dots,L,k=1,\dots,K$ are the Fourier
coefficients of the solution on $\Gamma_{k,l}$,
$h^{k}_{l}=y^{k}_{l}-y^{k}_{l-1},\zeta^{k}_{l}=e^{{\rm
i}\beta_{l,k}^{(n)}h^{k}_{l}}-e^{-{\rm i}\beta_{l,k}^{(n)}h^{k}_{l}}$, and
$\beta_{l,k}^{(n)}=\big{(}\kappa_{k,l}^{2}-\big{(}\frac{n\pi}{w_{k}}\big{)}^{2}\big{)}^{1/2},\Im\beta_{l,k}^{(n)}\geq
0.$

Continuing in a manner analogous to the discourse presented in Subsection 3.2,
we can deduce an equivalent connection formula, similar to the one provided in
(3.14), with the inclusion of an additional subscript $k$. For $l=1$, the
normal derivative of (3.19) on the aperture of the $k$-th cavity is given by

$\partial_{y}u_{1}(x,0^{-};k)=\sum\limits_{n=1}^{N}s_{k}^{(n)}u_{0,k}^{(n)}\sin\frac{n\pi(x-a_{k})}{w_{k}},$
(3.21)

where $s_{k}^{(n)}=-b_{1,k}^{(n)}+(a_{1,k}^{(n)})^{2}\hat{u}_{1,k}^{(n)}$,
with the coefficients $a_{l,k}^{(n)}$ and $b_{l,k}^{(n)}$ defined in
(3.15)–(3.16), with an additional subscript $k$. When $L_{k}=1$, the
coefficient is $s_{k}^{(n)}=-b_{1,k}^{(n)}$.

Likewise, using the transmission conditions:
$u(x,0^{+})=u_{1}(x,0^{-};k),\partial_{y}u(x,0^{+})=\partial_{y}u_{1}(x,0^{-};k),x\in(a_{k},b_{k})$,
upon substituting (3.5) with (3.21), and proceeding with a comparable
discussion as presented in Subsection 3.1, we can deduce the linear system
governing the Fourier coefficients of the solution on the apertures:

${\mathbb{D}}_{\rm TM}{\mathbb{U}}={\mathbb{M}}_{\rm
TM}{\mathbb{U}}+{\mathbb{F}},$ (3.22)

where

${\mathbb{U}}=\begin{bmatrix}U_{1}\\\ U_{2}\\\ \vdots\\\
U_{K}\end{bmatrix},\quad{\mathbb{F}}=\begin{bmatrix}F_{1}\\\ F_{2}\\\
\vdots\\\ F_{K}\end{bmatrix},\quad{\mathbb{D}}_{\rm
TM}=\begin{bmatrix}D^{(1)}_{\rm TM}&0&\cdots&0\\\ 0&D^{(2)}_{\rm
TM}&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\ 0&0&\cdots&D^{(K)}_{\rm
TM}\end{bmatrix}$

with $U_{k}=[u_{0,k}^{(1)},u_{0,k}^{(2)},\dots,u_{0,k}^{(N)}]^{\top}$,
$D^{(k)}_{\rm TM}=\frac{w_{k}}{2}{\rm diag}(s_{k}^{(1)},\dots,s_{k}^{(N)})$,
the elements of the $k$-th block of the vector ${\mathbb{F}}$ being given by

$F_{k}(m)=-2{\rm i}\beta\int_{0}^{w_{k}}e^{{\rm
i}\alpha(x+a_{k})}\sin\frac{m\pi x}{w_{k}}{\rm d}x,$

and the elements within the block matrix ${\mathbb{M}}_{\text{TM}}$ consist of
hypersingular integrals, which will be transformed into forms involving weakly
singular integrals. An examination of this transformation process will be
presented in Section 5.

Once the system (3.22) is solved, the Fourier coefficients on the inner
interfaces of the $k$-th cavity can be obtained through the connection
formula. The total field in the $l$-th layer of the $k$-th cavity can be
computed using (3.19)–(3.20).

##  4\. TE polarization

In this section, we provide a concise overview of the numerical solution for
the boundary value problem (2.3), (2.6), (2.10) in TE polarization, as its
fundamental concept aligns closely with that of the TM polarization.

###  4.1. A single empty cavity

Consider the scattering by a single empty cavity. Due to the boundary
condition (2.6), the total field can be approximated by the Fourier series of
cosine functions:

$u(x,y)=\sum\limits_{n=0}^{N}u^{(n)}(y)\cos\frac{n\pi(x-a)}{w}.$ (4.1)

Denote by $u^{(n)},n=0,1,\dots,N,$ the Fourier coefficients of the solution on
$\Gamma$. Substituting (4.1) into (2.3) and (2.6) yields that the Fourier
coefficients satisfy

$\left\\{\begin{aligned}
&u^{(n)^{\prime\prime}}(y)+\beta_{n}^{2}u^{(n)}(y)=0,\quad-h<y<0,\\\
&u^{(n)}(0)=u^{(n)},\quad u^{(n)^{\prime}}(-h)=0,\end{aligned}\right.$

which has a unique solution given by

$u^{(n)}(y)=\left\\{\begin{aligned} &\frac{1}{1+e^{2{\rm
i}\beta_{n}h}}\big{(}e^{-{\rm i}\beta_{n}y}+e^{2{\rm i}\beta_{n}h}e^{{\rm
i}\beta_{n}y}\big{)}u^{(n)}\quad&\text{if }~{}\beta_{n}\neq 0,\\\
&u^{(n)}\quad&\text{if }~{}\beta_{n}=0.\end{aligned}\right.$ (4.2)

Taking the normal derivative of (4.1) on $\Gamma$ from below and using (4.2)
gives

$\partial_{y}u(x,0^{-})=\sum\limits_{n=0}^{N}t^{(n)}u^{(n)}\cos\frac{n\pi(x-a)}{w},$
(4.3)

where $t^{(n_{0})}=0$ and $t^{(n)}={\rm i}\beta_{n}(e^{2{\rm
i}\beta_{n}h}-1)/(1+e^{2{\rm i}\beta_{n}h})$ for $n\neq n_{0}$.

As the cavity is assumed to be empty, and the wavenumber within the cavity
remains consistent with that in free space, the following transmission
conditions are imposed on $\Gamma$:
$u(x,0^{+})=u(x,0^{-}),\partial_{y}u(x,0^{+})=\partial_{y}u(x,0^{-})$.
Substituting (4.3) into (2.10) leads to

$\displaystyle\sum\limits_{n=0}^{N}u^{(n)}(0)\cos\frac{n\pi(x-a)}{w}$ (4.4)
$\displaystyle=-\frac{\rm
i}{2}\int_{\Gamma}H_{0}^{(1)}\left(\kappa|x-x^{\prime}|\right)\bigg{[}\sum\limits_{n=0}^{N}t^{(n)}u^{(n)}(0)\cos\frac{n\pi(x^{\prime}-a)}{w}\bigg{]}{\rm
d}x^{\prime}+g(x).$

Multiplying both sides of (4.4) with $\cos\frac{m\pi(x-a)}{w},m=0,1,\dots,N$,
and integrating on the aperture, we obtain

$D_{\rm TE}U=\tilde{M}_{\rm TE}U+G,\quad\tilde{M}_{\rm TE}(m,n)=t^{(n)}M_{\rm
TE}(m,n),$ (4.5)

where $U=[u^{(0)},u^{(1)},\dots,u^{(N)}]^{\top}$, $D_{\rm TE}$ is a diagonal
matrix whose diagonal entry is $\frac{w}{2}$ unless for $D_{\rm TE}(0,0)=w$,
and the $m$-th component of vector $G$ is

$G(m)=2\int_{\Gamma}e^{{\rm i}\alpha x}\cos\frac{m\pi(x-a)}{w}{\rm
d}x=2e^{{\rm i}\alpha a}\int_{0}^{w}e^{{\rm i}\alpha x}\cos\frac{m\pi
x}{w}{\rm d}x,$ (4.6)

and the elements of the matrix $M_{\rm TE}$, by a change of variables, are
expressed by

$M_{\rm TE}(m,n)=-\frac{\rm
i}{2}\left(\frac{w}{2\pi}\right)^{2}\int_{0}^{2\pi}\int_{0}^{2\pi}H_{0}^{(1)}\Big{(}\frac{\kappa_{0}w}{2\pi}|x-x^{\prime}|\Big{)}\cos\frac{nx^{\prime}}{2}\cos\frac{mx}{2}{\rm
d}x^{\prime}{\rm d}x.$ (4.7)

The entries of $M_{\rm TE}$ contain weakly singular integrals, and an
efficient quadrature formula for evaluating them is proposed in Section 5.
After solving (4.5), the total field in the cavity can be computed explicitly
by using (4.1)–(4.2).

###  4.2. A single multi-layered cavity

Denote by $u_{l}$ the total field in the $l$-th layer, which can be
approximated by

$u_{l}(x,y)=\sum\limits_{n=0}^{N}u_{l}^{(n)}(y)\cos\frac{n\pi(x-a)}{w}.$ (4.8)

A straightforward calculation shows that the Fourier coefficients satisfy
exactly the same system (3.12) in the $l$-th layer. Moreover, the transmission
conditions across the interface $\Gamma_{l}$ are

$u_{l}(x,y_{l})=u_{l+1}(x,y_{l}),\quad\frac{1}{\kappa_{l}^{2}}\partial_{y}u_{l}(x,y_{l})=\frac{1}{\kappa_{l+1}^{2}}\partial_{y}u_{l+1}(x,y_{l}),\quad
l=1,\dots,L-1,$

which give, after using the orthogonality properties of the cosine functions,
that

$\frac{1}{\kappa_{l}^{2}}\big{(}a_{l}^{(n)}u^{(n)}_{l-1}+b_{l}^{(n)}u^{(n)}_{l}\big{)}=-\frac{1}{\kappa_{l+1}^{2}}\big{(}b_{l+1}^{(n)}u^{(n)}_{l}+a_{l+1}^{(n)}u^{(n)}_{l+1}\big{)},\quad
l=1,\dots,L-1,$

where $a^{(n)}_{l},b^{(n)}_{l}$ for $l=1,\dots,L$ are defined in
(3.15)–(3.16). On the bottom of the cavity, we deduce from the boundary
condition (2.6) that

$a_{L}^{(n)}u^{(n)}_{L-1}+b_{L}^{(n)}u^{(n)}_{L}=0.$ (4.9)

Hence we can obtain a similar tri-diagonal system for the connection formula

$D_{\rm TE}^{(n)}\boldsymbol{u}^{(n)}=\boldsymbol{b}_{\rm TE}^{(n)},$ (4.10)

where
$\boldsymbol{u}^{(n)}=[u_{1}^{(n)},u_{2}^{(n)},\dots,u_{L}^{(n)}]^{\top}$,
$\boldsymbol{b}_{\rm
TE}^{(n)}=[-\frac{1}{\kappa_{1}^{2}}a_{1}^{(n)}u^{(n)}_{0},0,\dots,0]^{\top}$,
and the entries of the symmetric tri-diagonal matrix $D_{\rm TE}^{(n)}$ are
specified by

$D_{\rm
TE}^{(n)}(l,l)=\frac{b_{l}^{(n)}}{\kappa_{l}^{2}}+\frac{b_{l+1}^{(n)}}{\kappa_{l+1}^{2}},\quad
D^{(n)}(l,l+1)=\frac{a_{l+1}^{(n)}}{\kappa_{l+1}^{2}}$ (4.11)

with $b_{L+1}^{(n)}/\kappa_{L+1}^{2}=0$. In the case of $L=1$, the connection
formula (4.11) is reduced to (4.9).

Denote by
$\hat{\boldsymbol{u}}^{(n)}=[\hat{u}_{1}^{(n)},\hat{u}_{2}^{(n)},\dots,\hat{u}_{L}^{(n)}]^{\top}$
the solution of (4.10) with the right-hand side given by
$\hat{\boldsymbol{b}}^{(n)}=[1,0,...,0]^{\top},$ then we have
${\boldsymbol{u}}^{(n)}=-\frac{1}{\kappa_{1}^{2}}a_{1}^{(n)}u^{(n)}_{0}\hat{\boldsymbol{u}}^{(n)}.$
Taking the derivative of (4.8) on the aperture from below for $l=1$, and using
the transmission conditions $u(x,0^{+})=u_{1}(x,0^{-})$ and
$\frac{1}{\kappa_{0}^{2}}\partial_{y}u(x,0^{+})=\frac{1}{\kappa_{1}^{2}}\partial_{y}u_{1}(x,0^{-})$,
we obtain

$\displaystyle\partial_{y}u(x,0^{+})=\left(\frac{\kappa_{0}}{\kappa_{1}}\right)^{2}\partial_{y}u_{1}(x,0^{-})=\sum\limits_{n=0}^{N}{\hat{t}}^{(n)}u_{0}^{(n)}\cos\frac{\pi
n(x-a)}{w},$ (4.12)

where
${\hat{t}}^{(n)}=\big{(}\frac{\kappa_{0}}{\kappa_{1}}\big{)}^{2}\big{[}\frac{1}{\kappa_{1}^{2}}(a_{1}^{(n)})^{2}\hat{u}^{(n)}_{1}-b_{1}^{(n)}\big{]}$.
Substituting (4.12) into (2.10), multiplying both sides by
$\cos\frac{m\pi(x-a)}{w}$, and integrating it on the aperture, we obtain

$\hat{D}_{\rm TE}U=\hat{M}_{\rm TE}U+G,$ (4.13)

where $U=[u^{(0)},u^{(1)},\dots,u^{(N)}]^{\top}$, the entries of $G$ are given
in (4.6), the diagonal matrix $\hat{D}_{\rm TE}$ and the matrix $\hat{M}_{\rm
TE}$ are given by

$\hat{D}_{\rm TE}(m,m)=\left\\{\begin{aligned} &w\quad&\text{if }~{}m=0,\\\
&\frac{w}{2}\quad&\text{if }~{}m\neq 0,\end{aligned}\right.\qquad\hat{M}_{\rm
TE}(m,n)={\hat{t}}^{(n)}M_{\rm TE}(m,n)$

with the elements of $M_{\rm TE}$ being given in (4.7). After the system
(4.13) is solved, the solution in each layer can be computed explicitly by
(4.8).

###  4.3. Multiple multi-layered cavities

Let $u(x,y;k)$ denote the total field of the $k$-th cavity. The NtD operator
(2.10) on the $k$-th aperture can be reformulated as

$\displaystyle I_{\rm TE}(u)(x)=-\frac{\rm
i}{2}\sum\limits_{k=1}^{K}\int_{\Gamma_{k}}H_{0}^{(1)}(\kappa_{0}|x-x^{\prime}|)\partial_{y^{\prime}}u(x^{\prime},0;k){\rm
d}x^{\prime}.$ (4.14)

The solution $u(x,y;k)$ in the $l$-th layer can be approximated by

$u_{l}(x,y;k)=\sum\limits_{n=0}^{N}u_{l}^{(n)}(y;k)\cos\frac{n\pi(x-a_{k})}{w_{k}},\quad
y\in(y_{l}^{k},y_{l-1}^{k}).$ (4.15)

Following the same discussion in Section 4.2, we find the following expression
within the $k$-th cavity:

$\displaystyle\partial_{y}u(x,0^{+};k)=\left(\frac{\kappa_{0}}{\kappa_{k,1}}\right)^{2}\partial_{y}u_{1}(x,0^{-};k)=\sum\limits_{n=0}^{N}t_{k}^{(n)}u_{0,k}^{(n)}\cos\frac{n\pi(x-a_{k})}{w_{k}}.$

where
$t_{k}^{(n)}=\big{(}\frac{\kappa_{0}}{\kappa_{k,1}}\big{)}^{2}\big{(}-b_{1,k}^{(n)}+\frac{1}{\kappa_{k,1}^{2}}(a_{1,k}^{(n)})^{2}\hat{u}_{1,k}\big{)}$.

By substituting the above equation into (2.10) and (4.14), we deduce

${\mathbb{D}}_{\rm TE}{\mathbb{U}}={\mathbb{M}}_{\rm
TE}{\mathbb{U}}+{\mathbb{G}},$ (4.16)

where

${\mathbb{U}}=\begin{bmatrix}U_{1}\\\ U_{2}\\\ \vdots\\\
U_{K}\end{bmatrix},\quad{\mathbb{G}}=\begin{bmatrix}G_{1}\\\ G_{2}\\\
\vdots\\\ G_{K}\end{bmatrix},\quad{\mathbb{D}}_{\rm
TE}=\begin{bmatrix}D^{(1)}_{\rm TE}&0&\cdots&0\\\ 0&D^{(2)}_{\rm
TE}&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\ 0&0&\cdots&D^{(K)}_{\rm
TE}\end{bmatrix}$

with $U_{k}=[u_{0,k}^{(0)},u_{0,k}^{(1)},\dots,u_{0,k}^{(N)}]^{\top}$,
$D^{(k)}_{\rm TE}=w_{k}$ for $k=1$ and $w_{k}/2$ otherwise, the elements of
the $k$-th block of the vector ${\mathbb{G}}$ being defined as

$G_{k}(m)=2\int_{0}^{w_{k}}e^{{\rm i}\alpha(x+a_{k})}\cos\frac{m\pi
x}{w_{k}}{\rm d}x,$

and the block form of the matrix ${\mathbb{M}}_{\rm TE}$ is

${\mathbb{M}}_{\rm TE}=\begin{bmatrix}M_{1,1}&M_{1,2}&\cdots&M_{1,K}\\\
M_{2,1}&M_{2,2}&\cdots&M_{2,K}\\\ \vdots&\vdots&\ddots&\vdots\\\
M_{K,1}&M_{K,2}&\cdots&M_{K,K}\end{bmatrix}$

with the entries being given by

$M_{k,j}(m,n)=-\frac{{\rm
i}t^{(n)}_{j}}{2}\int_{0}^{w_{k}}\int_{0}^{w_{j}}H_{0}^{(1)}(\kappa_{0}|x+a_{k}-x^{\prime}-a_{j}|)\cos\frac{n\pi
x^{\prime}}{w_{j}}\cos\frac{m\pi x}{w_{k}}{\rm d}x^{\prime}{\rm d}x.$

Once the system (4.16) is solved, the Fourier coefficients on the inner
interfaces of the $k$-th cavity can be obtained using the connection formula.
Subsequently, the total field in the $l$-th layer of the $k$-th cavity can be
expressed by (4.15), where the coefficients are defined according to (3.13).

##  5\. Numerical quadratures

Since the first-order Hankel function $H_{1}^{(1)}(t)$ is singular at $t=0$,
the DtN operator (2.9) for the TM polarization is defined in the sense of the
Hadamard finite-part. Developing a high-order numerical integral quadrature
formula is challenging for hypersingular integrals. In this section, we
present an alternative transparent boundary condition (TBC) containing weakly
singular integrals and propose a high-order numerical quadrature rule for
these specific cases.

###  5.1. An alternative transparent boundary condition

We investigate the transparent boundary condition (2.8) in TM polarization and
explore the application of numerical quadrature techniques to handle the
weakly singular integrals involved in the transparent boundary conditions.

######  Lemma 5.1.

The transparent boundary condition (2.8) can be reformulated as follows:

$\displaystyle\partial_{x_{2}}u(x_{1},0)$ $\displaystyle=\frac{{\rm
i}\kappa_{0}^{2}}{2}\int_{\Gamma}u(y_{1},0)H_{0}^{(1)}(\kappa_{0}|x_{1}-y_{1}|){\rm
d}y_{1}$ $\displaystyle\quad+\frac{\rm
i}{2}\partial_{x_{1}}\left[\int_{\Gamma}\partial_{y_{1}}u(y_{1},0)H_{0}^{(1)}(\kappa_{0}|x_{1}-y_{1}|){\rm
d}y_{1}\right]-2{\rm i}\beta e^{{\rm i}\alpha x_{1}}.$ (5.1)

###### Proof.

Let $G_{\rm TM}(\boldsymbol{x},\boldsymbol{y})$ represent the Dirichlet
Green’s function in the upper half-plane for TM polarization. Specifically, it
is given by

$G_{\rm TM}(\boldsymbol{x};\boldsymbol{y})=\frac{\rm
i}{4}H_{0}^{(1)}(\kappa_{0}|\boldsymbol{x}-\boldsymbol{y}|)-\frac{\rm
i}{4}H_{0}^{(1)}(\kappa_{0}|\boldsymbol{x}-\boldsymbol{y}^{\prime}|),$

where $\boldsymbol{x}=(x_{1},x_{2})$, $\boldsymbol{y}=(y_{1},y_{2})$, and
$\boldsymbol{y}^{\prime}=(y_{1},-y_{2})$. It follows from the Sommerfeld
radiation condition (2.7) and Green’s integral formula that for any
$(x_{1},x_{2})\in\mathbb{R}_{+}^{2}$, the scattered field $u^{s}(x_{1},x_{2})$
can be expressed as

$\displaystyle
u^{s}(x_{1},x_{2})=-\int_{\Gamma_{g}\cup\Gamma}\partial_{y_{2}}u^{s}(y_{1},0)G_{\rm
TM}(x_{1},x_{2};y_{1},0){\rm d}y_{1}$
$\displaystyle+\int_{\Gamma_{g}\cup\Gamma}u^{s}(y_{1},0)\partial_{y_{2}}G_{\rm
TM}(x_{1},x_{2};y_{1},0){\rm d}y_{1}.$

Based on the homogeneous Dirichlet boundary condition (2.5), the total field
satisfies

$\displaystyle u(x_{1},x_{2})-u^{b}(x_{1},x_{2})$
$\displaystyle=-\int_{\Gamma_{g}\cup\Gamma}\Big{[}\partial_{y_{2}}(u(y_{1},0)-u^{b}(y_{1},0))G_{\rm
TM}(x_{1},x_{2};y_{1},0)$
$\displaystyle\quad-(u(y_{1},0)-u^{b}(y_{1},0))\partial_{y_{2}}G_{\rm
TM}(x_{1},x_{2};y_{1},0)\Big{]}{\rm d}y_{1}$
$\displaystyle=\int_{\Gamma}u(y_{1},0)\partial_{y_{2}}G_{\rm
TM}(x_{1},x_{2};y_{1},0){\rm d}y_{1}.$ (5.2)

where $u^{b}=u^{i}-u^{r}=e^{{\rm i}(\alpha x_{1}-\beta x_{2})}-e^{{\rm
i}(\alpha x_{1}+\beta x_{2})}$ is referred to the background field.

Define an auxiliary function as the Dirichlet Green’s function in the upper
half-plane for TE polarization:

$G_{\rm TE}(\boldsymbol{x},\boldsymbol{y})=\frac{\rm
i}{4}H_{0}^{(1)}(\kappa_{0}|\boldsymbol{x}-\boldsymbol{y}|)+\frac{\rm
i}{4}H_{0}^{(1)}(\kappa_{0}|\boldsymbol{x}-\boldsymbol{y}^{\prime}|).$

Noting $\partial_{y_{2}}G_{\rm
TM}(x_{1},x_{2};y_{1},y_{2})=-\partial_{x_{2}}G_{\rm
TE}(x_{1},x_{2};y_{1},y_{2})$, we have from (5.1) that

$u(x_{1},x_{2})-u^{b}(x_{1},x_{2})=-\int_{\Gamma}u(y_{1},0)\partial_{x_{2}}G_{\rm
TE}(x_{1},x_{2};y_{1},0){\rm d}y_{1}.$ (5.3)

Taking the partial derivative of (5.3) with respect to $x_{2}$, we obtain

$\displaystyle\partial_{x_{2}}u(x_{1},x_{2})-\partial_{x_{2}}u^{b}(x_{1},x_{2})=-\int_{\Gamma}u(y_{1},0)\partial^{2}_{x_{2}}G_{\rm
TE}(x_{1},x_{2};y_{1},0){\rm d}y_{1}$
$\displaystyle=\kappa_{0}^{2}\int_{\Gamma}u(y_{1},0)G_{\rm
TE}(x_{1},x_{2};y_{1},0){\rm d}y_{1}$ $\displaystyle\quad-\frac{\rm
i}{4}\int_{\Gamma}u(y_{1},0)\left[\partial_{x_{1}y_{1}}^{2}H_{0}^{(1)}(\kappa_{0}|\boldsymbol{x}-\boldsymbol{y}|)+\partial_{x_{1}y_{1}}^{2}H_{0}^{(1)}(\kappa_{0}|\boldsymbol{x}-\boldsymbol{y}^{\prime}|)\right]{\rm
d}y_{1}$ $\displaystyle=\kappa_{0}^{2}\int_{\Gamma}u(y_{1},0)G_{\rm
TE}(x_{1},x_{2};y_{1},0){\rm d}y_{1}$ $\displaystyle\quad+\frac{\rm
i}{4}\int_{\Gamma}\partial_{y_{1}}u(y_{1},0)\left[\partial_{x_{1}}H_{0}^{(1)}(\kappa_{0}|\boldsymbol{x}-\boldsymbol{y}|)+\partial_{x_{1}}H_{0}^{(1)}(\kappa_{0}|\boldsymbol{x}-\boldsymbol{y}^{\prime}|)\right]{\rm
d}y_{1},$ (5.4)

where we have employed integration by parts and utilized the following facts:

$\displaystyle\partial_{x_{2}}^{2}G_{\rm TE}(\boldsymbol{x};\boldsymbol{y})$
$\displaystyle=-\partial_{x_{1}}^{2}G_{\rm
TE}(\boldsymbol{x};\boldsymbol{y})-\kappa_{0}^{2}G_{\rm
TE}(\boldsymbol{x};\boldsymbol{y}),$ $\displaystyle\partial_{x_{1}}^{2}G_{\rm
TE}(\boldsymbol{x};\boldsymbol{y})$ $\displaystyle=\frac{\rm
i}{4}\left[\partial_{x_{1}}^{2}H_{0}^{(1)}(\kappa_{0}|\boldsymbol{x}-\boldsymbol{y}|)+\partial_{x_{1}}^{2}H_{0}^{(1)}(\kappa_{0}|\boldsymbol{x}-\boldsymbol{y}^{\prime}|)\right]$
$\displaystyle=-\frac{\rm
i}{4}\left[\partial_{x_{1}y_{1}}^{2}H_{0}^{(1)}(\kappa_{0}|\boldsymbol{x}-\boldsymbol{y}|)+\partial_{x_{1}y_{1}}^{2}H_{0}^{(1)}(\kappa_{0}|\boldsymbol{x}-\boldsymbol{y}^{\prime}|)\right].$

Taking $x_{2}\rightarrow 0^{+}$ in (5.1), and using the continuity of the
total field and the single layer potential across the aperture $\Gamma$, we
deduce

$\displaystyle\partial_{x_{2}}u(x_{1},0)-\partial_{x_{2}}u^{b}(x_{1},0)$
$\displaystyle=\kappa_{0}^{2}\int_{\Gamma}u(y_{1},0)G_{\rm
TE}(x_{1},0;y_{1},0){\rm d}y_{1}+\frac{\rm
i}{2}\int_{\Gamma}\partial_{y_{1}}u(y_{1},0)\partial_{x_{1}}H_{0}^{(1)}(\kappa_{0}|x_{1}-y_{1}|){\rm
d}y_{1}$ $\displaystyle=\frac{{\rm
i}\kappa_{0}^{2}}{2}\int_{\Gamma}u(y_{1},0)H_{0}^{(1)}(\kappa_{0}|x_{1}-y_{1}|){\rm
d}y_{1}+\frac{\rm
i}{2}\partial_{x_{1}}\Big{[}\int_{\Gamma}\partial_{y_{1}}u(y_{1},0)H_{0}^{(1)}(\kappa_{0}|x_{1}-y_{1}|){\rm
d}y_{1}\Big{]},$

which completes the proof by noting $\partial_{x_{2}}u^{b}(x_{1},0)=-2{\rm
i}\beta e^{{\rm i}\alpha x_{1}}$. ∎

Substituting the expansion (3.1) into (5.1), we get

$\displaystyle\partial_{x_{2}}u(x_{1},0)$ $\displaystyle=\frac{{\rm
i}\kappa_{0}^{2}}{2}\int_{0}^{w}\bigg{(}\sum\limits_{n=1}^{N}u^{(n)}\sin\frac{n\pi
y_{1}}{w}H_{0}^{(1)}(\kappa_{0}|x_{1}-y_{1}-a|)\bigg{)}{\rm d}y_{1}$
$\displaystyle\quad+\frac{\rm
i}{2}\partial_{x_{1}}\bigg{[}\int_{0}^{w}\sum\limits_{n=1}^{N}u^{(n)}\frac{n\pi}{w}\cos\frac{n\pi
y_{1}}{w}H_{0}^{(1)}(\kappa_{0}|x_{1}-y_{1}-a|){\rm d}y_{1}\bigg{]}-2{\rm
i}\beta e^{{\rm i}\alpha x_{1}}.$ (5.5)

By multiplying both sides of (5.1) by $\sin\frac{m\pi(x_{1}-a)}{w}$ and
integrating over $\Gamma$, we obtain the equivalent entries of the matrix
$M_{\rm TM}$ in (3.8):

$\displaystyle M_{\rm TM}(m,n)$ $\displaystyle=\frac{{\rm
i}\kappa_{0}^{2}}{2}\int_{0}^{w}\left[\int_{0}^{w}\sin\frac{n\pi
y}{w}H_{0}^{(1)}\left(\kappa_{0}|x-y|\right){\rm d}y\right]\sin\frac{m\pi
x}{w}{\rm d}x$ $\displaystyle\quad-\frac{{\rm
i}mn\pi^{2}}{2w^{2}}\int_{0}^{w}\left[\int_{0}^{w}\cos\frac{n\pi
y}{w}H_{0}^{(1)}(\kappa_{0}|x-y|){\rm d}y\right]\cos\frac{m\pi x}{w}{\rm d}x.$
(5.6)

For multiple multi-layered cavities, following the same discussion, we can
derive the TBC on the $k$-th cavity as

$\displaystyle\partial_{x_{2}}u(x_{1},0;k)$ $\displaystyle=\frac{{\rm
i}\kappa_{0}^{2}}{2}\sum\limits_{j=1}^{K}\int_{\Gamma_{j}}u(y_{1},0;j)H_{0}^{(1)}(\kappa_{0}|x_{1}-y_{1}|){\rm
d}y_{1}$ $\displaystyle\quad+\frac{\rm
i}{2}\sum\limits_{j=1}^{K}\partial_{x_{1}}\left[\int_{\Gamma_{j}}\partial_{y_{1}}u(y_{1},0;j)H_{0}^{(1)}(\kappa_{0}|x_{1}-y_{1}|){\rm
d}y_{1}\right]-2{\rm i}\beta e^{{\rm i}\alpha x_{1}}.$ (5.7)

Substituting (3.19) into (5.1) and using a change of variables, we deduce

$\displaystyle\partial_{x_{2}}u(x_{1},0;k)=\frac{{\rm
i}\kappa_{0}^{2}}{2}\sum\limits_{j=1}^{K}\int_{0}^{w_{j}}\left(\sum\limits_{n=1}^{N}u^{(n)}_{0,j}\sin\frac{n\pi
y_{1}}{w_{j}}H_{0}^{(1)}(\kappa_{0}|x_{1}-y_{1}-a_{j}|)\right){\rm d}y_{1}$
$\displaystyle+\frac{\rm
i}{2}\sum\limits_{j=1}^{K}\partial_{x_{1}}\left[\int_{0}^{w_{j}}\sum\limits_{n=1}^{N}u^{(n)}_{0,j}\frac{n\pi}{w_{j}}\cos\frac{n\pi
y_{1}}{w_{j}}H_{0}^{(1)}(\kappa_{0}|x_{1}-y_{1}-a_{j}|){\rm
d}y_{1}\right]-2{\rm i}\beta e^{{\rm i}\alpha x_{1}}.$ (5.8)

By multiplying both sides of (5.1) by $\sin\frac{\pi m(x_{1}-a_{k})}{w_{k}}$
and integrating over $\Gamma_{k}$, we can obtain a new form of the matrix
${\mathbb{M}}_{\rm TM}$ in (3.22), which only contains weakly singular
integrals and has entries respresented as follows:

$\displaystyle M_{k,j}(m,n)$ $\displaystyle=\frac{{\rm
i}\kappa_{0}^{2}}{2}\int_{0}^{w_{k}}\int_{0}^{w_{j}}\sin\frac{n\pi
y}{w_{j}}H_{0}^{(1)}\left(\kappa_{0}|x+a_{k}-y-a_{j}|\right)\sin\frac{m\pi
x}{w_{k}}{\rm d}y{\rm d}x$ $\displaystyle\quad-\frac{{\rm
i}mn\pi^{2}}{2w_{j}w_{k}}\int_{0}^{w_{k}}\int_{0}^{w_{j}}\cos\frac{n\pi
y}{w_{j}}H_{0}^{(1)}\left(\kappa_{0}|x+a_{k}-y-a_{j}|\right)\cos\frac{m\pi
x}{w_{k}}{\rm d}y{\rm d}x,$

where $k,j=1,\dots,K$, and $m,n=1,\dots,N$.

###  5.2. Numerical quadratures

The matrices ${\mathbb{M}}_{\rm TM}$ and ${\mathbb{M}}_{\rm TE}$ contain
weakly singular integrals in the following forms subsequent to the use of a
change of variables:

$\displaystyle\int_{0}^{2\pi}\int_{0}^{2\pi}\sin\frac{ns}{2}H_{0}^{(1)}\left(\frac{\kappa_{0}w}{2\pi}|s-t|\right)\sin\frac{mt}{2}{\rm
d}s{\rm d}t,$ (5.9)
$\displaystyle\int_{0}^{2\pi}\int_{0}^{2\pi}\cos\frac{ns}{2}H_{0}^{(1)}\left(\frac{\kappa_{0}w}{2\pi}|s-t|\right)\cos\frac{mt}{2}{\rm
d}s{\rm d}t.$

In this section, we develop an efficient approach to accurately evaluate
(5.9).

######  Lemma 5.2.

If $m+n$ is odd, then

$\displaystyle\int_{0}^{2\pi}\left(\int_{0}^{2\pi}H_{0}^{(1)}\left(\frac{\kappa_{0}w}{2\pi}|s-t|\right)\sin\frac{ns}{2}{\rm
d}s\right)\sin\frac{mt}{2}{\rm d}t$ $\displaystyle=0,$ (5.10)
$\displaystyle\int_{0}^{2\pi}\left(\int_{0}^{2\pi}H_{0}^{(1)}\left(\frac{\kappa_{0}w}{2\pi}|s-t|\right)\cos\frac{ns}{2}{\rm
d}s\right)\cos\frac{mt}{2}{\rm d}t$ $\displaystyle=0.$ (5.11)

###### Proof.

We only show the proof of (5.10) since the proof is similar to (5.11). For any
$t\in[0,2\pi]$, define

$F(t)=\int_{0}^{2\pi}H_{0}^{(1)}\left(\frac{\kappa_{0}w}{2\pi}|s-t|\right)\sin\frac{ns}{2}{\rm
d}s.$

Let $\xi=2\pi-s$, we have from a simple calculation that

$\displaystyle F(2\pi-t)$
$\displaystyle=\int_{0}^{2\pi}H_{0}^{(1)}\left(\frac{\kappa_{0}w}{2\pi}|s-2\pi+t|\right)\sin\frac{ns}{2}{\rm
d}s$
$\displaystyle=\int_{0}^{2\pi}H_{0}^{(1)}\left(\frac{\kappa_{0}w}{2\pi}|2\pi-\xi-2\pi+t|\right)\sin\frac{n}{2}(2\pi-\xi){\rm
d}\xi=(-1)^{(n+1)}F(t).$

A straightforward calculation shows that

$\displaystyle\int_{0}^{2\pi}\left(\int_{0}^{2\pi}H_{0}^{(1)}\left(\frac{\kappa_{0}w}{2\pi}|s-t|\right)\sin\frac{ns}{2}{\rm
d}s\right)\sin\frac{mt}{2}{\rm d}t$
$\displaystyle=\int_{0}^{\pi}F(t)\sin\frac{mt}{2}{\rm
d}t+\int_{\pi}^{2\pi}F(t)\sin\frac{mt}{2}{\rm d}t$
$\displaystyle=\int_{0}^{\pi}F(t)\sin\frac{mt}{2}{\rm
d}t+(-1)^{(m+n)}\int_{0}^{\pi}F(t)\sin\frac{mt}{2}{\rm d}t=0,$

if $m+n$ is odd, which completes the proof. ∎

Next, we address numerical quadrature in scenarios where the sum of $m$ and
$n$ results in an even number. The established technique outlined in [10,
Section 3.5] for computing weakly singular integrals associated with the
zeroth-order Hankel function necessitates the integrand function to be
$2\pi$-periodic. Regrettably, this condition cannot be met in our current
context.

It is clear to note from (5.9) that

$\displaystyle\int_{0}^{2\pi}\int_{0}^{2\pi}\sin\frac{ns}{2}H_{0}^{(1)}\Big{(}\frac{\kappa_{0}w}{2\pi}|s-t|\Big{)}\sin\frac{mt}{2}{\rm
d}s{\rm d}t$
$\displaystyle=\int_{0}^{2\pi}\int_{0}^{2\pi}\sin\frac{ns}{2}\left[H_{0}^{(1)}\Big{(}\frac{\kappa_{0}w}{2\pi}|s-t|\Big{)}-\frac{2{\rm
i}}{\pi}J_{0}\Big{(}\frac{\kappa_{0}w}{2\pi}|s-t|\Big{)}\ln|s-t|\right]\sin\frac{mt}{2}{\rm
d}s{\rm d}t$ $\displaystyle\quad+\frac{2{\rm
i}}{\pi}\int_{0}^{2\pi}\int_{0}^{2\pi}J_{0}\Big{(}\frac{\kappa_{0}w}{2\pi}|s-t|\Big{)}\ln|s-t|\sin\frac{ns}{2}\sin\frac{mt}{2}{\rm
d}s{\rm d}t.$ (5.12)

Recalling the power series of the Bessel function (cf. [25])

$\displaystyle H_{0}^{(1)}(z)$ $\displaystyle=J_{0}(z)+{\rm i}Y_{0}(z)$
$\displaystyle=J_{0}(z)+\frac{2{\rm
i}}{\pi}\Big{(}\ln\frac{z}{2}+\gamma\Big{)}J_{0}(z)+\frac{2{\rm
i}}{\pi}\sum\limits_{k=1}^{\infty}(-1)^{k+1}\bigg{(}\sum\limits_{j=1}^{k}\frac{1}{j}\bigg{)}\frac{\left(\frac{1}{4}z^{2}\right)^{k}}{(k!)^{2}},$

we know that the first term of (5.2) is analytic. Moreover, we have

$\lim\limits_{t\rightarrow
s}H_{0}^{(1)}\Big{(}\frac{\kappa_{0}w}{2\pi}|s-t|\Big{)}-\frac{2{\rm
i}}{\pi}J_{0}\Big{(}\frac{\kappa_{0}w}{2\pi}|s-t|\Big{)}\ln|s-t|=1+\frac{2{\rm
i}}{\pi}\gamma+\frac{2{\rm i}}{\pi}\ln\frac{\kappa_{0}w}{4\pi}.$

The integrand function of the second term is weakly singular, and directly
applying a high order quadrature rule is not feasible. Our idea is to
repeatedly apply integration by parts to increase the regularity of the
integrand function.

Specifically, it follows from the expansion of Bessel function (cf. [25]) that

$\displaystyle
J_{0}\Big{(}\frac{\kappa_{0}w}{2\pi}|s-t|\Big{)}=\sum\limits_{k=0}^{\infty}(-1)^{k}\frac{1}{(k!)^{2}}\Big{[}\frac{\kappa_{0}w}{4\pi}(s-t)\Big{]}^{2k}=\sum\limits_{k=0}^{K}(-1)^{k}\frac{1}{(k!)^{2}}\Big{[}\frac{\kappa_{0}w}{4\pi}(s-t)\Big{]}^{2k}$
$\displaystyle+\bigg{[}J_{0}\Big{(}\frac{\kappa_{0}w}{2\pi}|s-t|\Big{)}-\sum\limits_{k=0}^{K}(-1)^{k}\frac{1}{(k!)^{2}}\Big{(}\frac{\kappa_{0}w}{4\pi}(s-t)\Big{)}^{2k}\bigg{]}.$
(5.13)

Since the second term has higher regularity of $C^{2K+2}$, substituting (5.2)
into (5.2), we only need to handle the integral

$S_{k}(n,m)=\int_{0}^{2\pi}\int_{0}^{2\pi}(t-s)^{k-1}\ln|t-s|\sin\frac{mt}{2}\sin\frac{ns}{2}{\rm
d}s{\rm d}t.$

Through straightforward calculations and integration by parts, we deduce

$\displaystyle
S_{k}(n,m)=-\frac{1}{k(k+1)}\left(\frac{m}{2}\right)^{2}S_{k+2}(n,m)+\frac{1}{k(k+1)^{2}}\left(\frac{m}{2}\right)^{2}T^{(1)}_{k}(n,m)$
$\displaystyle+\frac{m}{2k^{2}}T_{k}^{(2)}(n,m)-\frac{m}{2k(k+1)^{2}}\left[(-1)^{m+n}-(-1)^{k}\right]T_{k}^{(3)}(n)$
$\displaystyle+\frac{m}{2k(k+1)}\left[(-1)^{m+n}-(-1)^{k}\right]W_{k}(n),$
(5.14)

where

$\displaystyle
T^{(1)}_{k}(n,m)=\int_{0}^{2\pi}\int_{0}^{2\pi}\sin\frac{ns}{2}(t-s)^{k+1}\sin\frac{mt}{2}{\rm
d}s{\rm d}t,$ $\displaystyle
T_{k}^{(2)}(n,m)=\int_{0}^{2\pi}\int_{0}^{2\pi}\sin\frac{ns}{2}(t-s)^{k}\cos\frac{mt}{2}{\rm
d}s{\rm d}t,$ $\displaystyle
T_{k}^{(3)}(n)=\int_{0}^{2\pi}\sin\frac{ns}{2}s^{k+1}{\rm d}s,\quad
W_{k}(n)=\int_{0}^{2\pi}\sin\frac{ns}{2}s^{k+1}\ln s\,{\rm d}s.$

The regularity of $W_{k}$ can be increased by applying the same idea.
Explicitly, we have

$\displaystyle
W_{k}(n)=A_{k}(n)-\left(\frac{n}{2}\right)^{2}\frac{1}{(k+2)(k+3)}W_{k+2}(n),$

where

$\displaystyle
A_{k}(n)=\frac{n}{2}\frac{1}{k+2}\left(\frac{1}{k+2}+\frac{1}{k+3}\right)\int_{0}^{2\pi}s^{k+2}\cos\frac{ns}{2}$
$\displaystyle-\frac{(-1)^{n}n}{2(k+2)(k+3)}\left(2\pi\right)^{k+3}\ln(2\pi).$

The corresponding parts of the integrals involving cosine functions can be
evaluated similarly as follows:

$\displaystyle P_{k}(n,m)$
$\displaystyle=\int_{0}^{2\pi}\int_{0}^{2\pi}(t-s)^{k-1}\ln|t-s|\cos\frac{mt}{2}\cos\frac{ns}{2}{\rm
d}s{\rm d}t$
$\displaystyle=-\left(\frac{m}{2}\right)^{2}\frac{1}{k(k+1)}P_{k+2}(n,m)+\left(\frac{m}{2}\right)^{2}\frac{1}{k(k+1)^{2}}U_{k}^{(1)}(n,m)-\frac{m}{2k^{2}}U_{k}^{(2)}(n,m)$
$\displaystyle\quad-\left[(-1)^{(m+n)}-(-1)^{k}\right]\frac{1}{k^{2}}U_{k}^{(3)}(n)+\left[(-1)^{(m+n)}-(-1)^{k}\right]\frac{1}{k}X_{k}(n),$
(5.15)

where

$\displaystyle
U_{k}^{(1)}(n,m)=\int_{0}^{2\pi}\int_{0}^{2\pi}\cos\frac{ns}{2}(t-s)^{k+1}\cos\frac{mt}{2}{\rm
d}s{\rm d}t,$ $\displaystyle
U_{k}^{(2)}(n,m)=\int_{0}^{2\pi}\int_{0}^{2\pi}\cos\frac{ns}{2}(t-s)^{k}\sin\frac{mt}{2}{\rm
d}s{\rm d}t,$ $\displaystyle
U_{k}^{(3)}(n)=\int_{0}^{2\pi}\cos\frac{ns}{2}s^{k}{\rm d}s,\quad
X_{k}(n)=\int_{0}^{2\pi}\cos\frac{ns}{2}s^{k}\ln s{\rm d}s.$

Moreover, the function $X_{k}$ satisfies the recurrence relation

$X_{k}(n)=Y_{k}(n)-\left(\frac{n}{2}\right)^{2}\frac{1}{(k+1)(k+2)}X_{k+2}(n),$

where

$\displaystyle
Y_{k}(n)=-\frac{n}{2}\frac{1}{k+1}\left(\frac{1}{k+1}+\frac{1}{k+2}\right)\int_{0}^{2\pi}s^{k+1}\sin\frac{ns}{2}{\rm
d}s$
$\displaystyle+\frac{(-1)^{n}(2\pi)^{k+1}}{\left(k+1\right)^{2}}\left[(k+1)\ln(2\pi)-1\right].$

By using the recurrence relations (5.2) and (5.2), we can iteratively enhance
the regularity of the integrand functions in (5.2) until they meet the
requirements of integral quadrature formulas, such as the Gaussian quadrature
formula.

##  6\. Numerical Experiments

In this section, we present a series of numerical examples to provide
compelling evidence regarding the efficacy of the method in precisely
analyzing electromagnetic scattering phenomena linked to rectangular cavities.

###  6.1. Order of accuracy

First, we assess the accuracy of the proposed method, demonstrating its
precision and reliability in handling the cavity scattering problems in both
TM and TE polarizations. In this experiment, we aim to evaluate the convergent
order of the proposed method in the $L^{2}$-norm. Due to the unavailability of
an analytic solution, we conduct a comparison of the results obtained on a
finer mesh.

The cavity is situated in the interval $[-0.5,0.5]$ with a depth of $h=1.5$.
It is subjected to illumination by a plane wave with a wavenumber of
$\kappa_{0}=1.5$ and an incidence angle of $\theta=\pi/9$. For the numerical
evaluation of (5.9), we opt for the composite 4-point Gaussian quadrature
formula as a representative numerical integration technique. We set the
truncation number to $N=30$. It is worth mentioning that by employing higher-
order numerical quadrature formulas, even higher orders of convergence can be
achieved. Figure 2 displays the results for both TM and TE polarizations. The
$x$-axis represents the mesh size, while the $y$-axis corresponds to the error
measured in the $L^{2}$-norm. The red dashed line depicts the convergent order
of the proposed method. For comparison purposes, we also include the solid
blue curve, which represents the theoretical eighth-order convergence. The
results clearly demonstrate that the proposed method exhibits an accuracy with
an order of $O(h^{8})$, thereby affirming its effectiveness and reliability in
numerical simulations.

Figure 2. Example 6.1: (Left) Order of convergence in TM polarization; (Right)
Ordr of convergence in TE polarization.

###  6.2. Radar cross-section

Next, we explore the application of our proposed method in radar cross-section
(RCS) computations, which highlights its capability to accurately predict the
scattering behavior of practical radar systems.

In two dimensions, the RCS is defined by (cf. [16])

$\sigma(\varphi):=\lim\limits_{r=\rightarrow\infty}2\pi
r\frac{|u^{s}(x,y;\varphi)|^{2}}{|u^{i}(x,y;\theta)|^{2}},$

where $u^{i}$ is the incident field, $u^{s}$ denotes the scattered field, and
$\theta$ and $\varphi$ represent the incident and observation angles,
respectively. When $\theta$ and $\varphi$ are equal, $\sigma$ is referred to
as the backscatter RCS, and it is defined as

$\text{Backscatter RCS}(\sigma)(\varphi)=10\log_{10}\sigma(\varphi){\rm dB}.$

In TM polarization, with measurements taken on the aperture of the cavity, the
RCS can be expressed as follows (cf. [31]):

$\sigma(\varphi)=\kappa_{0}\Big{|}\sin\varphi\int_{\Gamma}u^{s}(x,0)e^{{\rm
i}\kappa_{0}\cos\varphi x}{\rm d}x\Big{|}^{2}.$

In this experiment, we replicate a benchmark example documented in [16].
Specifically, we consider the backscatter RCS of a single rectangular cavity
in TM polarization. The cavity has a width $w=\lambda$ and a depth
$h=\lambda/4$, where $\lambda=2\pi/\kappa_{0}$ represents the wavelength in
free space. For our analysis, we set $\kappa_{0}=32\pi$ and $N=150$, along
with an incident angle of $\theta=\pi/3$. The obtained numerical backscatter
RCS results are shown in Figure 3. To provide a basis for comparison, we also
present the results obtained using the adaptive finite element TBC method
[31]. The adaptive TBC method is terminated when the total number of nodal
points reaches 15000. In the left part of Figure 3, we plot the results for
the empty cavity. The solid red line represents the backscatter RCS computed
by our proposed method, while the blue circle points represent the results
obtained from the adaptive TBC method. In the right part of Figure 3, we
consider the cavity filled with a lossy medium characterized by an electric
permittivity of $\epsilon=4+{\text{i}}$ and a magnetic permeability of
$\mu=1$. Similar to the empty cavity case, we depict the results of the
proposed method using a solid red line and the adaptive TBC method using blue
circle points. The results clearly demonstrate the consistency between the
proposed method and the adaptive TBC method. However, it is worth noting that
the proposed method only requires solving a small scale system with 151
unknowns, making it significantly more efficient in accurately computing the
backscatter RCS of the cavity compared to the adaptive TBC method. This
highlights the efficiency and accuracy of the proposed method in tackling the
backscatter RCS computation for the considered cavity configuration.

Figure 3. Example 6.2: (Left) The backscatter RCS of an empty cavity; (Right)
The backscatter RCS of a lossy cavity.

###  6.3. Subwavelength enhancement

In this experiment, we explore the efficacy of our proposed method in
analyzing field enhancement phenomena in subwavelength structures. By focusing
on TE polarization, we can demonstrate potential applications of the method in
nanophotonics and metamaterial research.

First, we consider the case of a single cavity with the following dimensions:
cavity width $w=50\,\mu m$ and cavity depth $h=1\,cm$. The cavity is
illuminated from above at an incident angle of $\theta=\pi/6$. To study the
behavior of the electric field enhancement, we vary the wavelength in a range
from $0.2\,cm$ to $62\,cm$. Figure 4 presents the plot of the electric field
enhancement factor $Q_{E}$ against the wavenumber $\kappa$, where $Q_{E}$ is
defined as the ratio of the $L^{2}$-norm of the electric field $u$ to the
$L^{2}$-norm of the incident electric field $u^{i}$, both integrated over the
cavity domain $D$, i.e.,

$Q_{E}=\frac{\|u\|_{L^{2}(D)}}{\|u^{i}\|_{L^{2}(D)}}.$

As depicted in the left part of Figure 4, it can be observed that the
enhancement factors exhibit peaks at resonant frequencies located near
$\kappa=\pi/2+n\pi$, where $n$ is an integer. These numerical findings are
consistent with the theoretical results reported in [14], highlighting the
capability of our proposed method in capturing the electric field enhancement
phenomena in the subwavelength cavity structure.

In the case of two cavities, we position them in the ground with a width of
$w=0.2\,\mu m$, a depth of $h=1.5\,\mu m$, and a distance of $d=0.5\,\mu m$
between them. The two cavities are illuminated from above by a plane wave with
an incident angle of $\theta=-\pi/9$. For our analysis, we vary the wavelength
in the range from $500\,cm$ to $2500\,cm$. As shown in the right part of
Figure 4, the blue dashed curve represents the enhancement factor for a single
cavity with the same width and depth, exhibiting a resonance frequency around
$\kappa_{0}=9000$. On the other hand, the red dashed line and black solid line
in the plot represent the enhancement factors for the left and right cavity,
respectively. It is evident from the results that the type of enhancement for
the two cavities is different. The left cavity exhibits a single peak near
$\kappa_{0}$, while the right cavity shows two resonance frequencies around
$\kappa_{0}$. The enhancement first peaks to the left of $\kappa_{0}$ and then
decreases dramatically. After that, the enhancement achieves a second peak to
the right of $\kappa_{0}$. These results demonstrate that the enhancement
caused by the two cavities exhibits both an antisymmetric mode and a
symmetrical mode. The numerical findings are in line with the theoretical
results reported in [3] and the experimental observations conducted in [26].

Figure 4. Example 6.3: (Left) The electric field enhancement factor by a
single cavity; (Right) The electric field enhancement factor by two cavities.

###  6.4. Multiple multi-layered cavities

Finally, we present a challenging scenario involving the scattering by
multiple cavities, which demonstrates the suitability of our proposed method
in effectively handling complex configurations of multiple multi-layered
cavities.

In this experiment, we consider both TM and TE polarizataions for a
configuration comprising three cavities, situated at the following locations:
$[-0.6,-0.1]\times[0,-0.1]\cup[0,0.2]\times[0,-0.5]\cup[0.3,0.6]\times[0,-0.3]$.
The first cavity is assumed to be empty, while the second cavity is filled
with a three-layered medium. This medium is characterized by wavenumbers
$\kappa=\pi,2\pi$, and $10\pi$, and it is separated at the inner interfaces
$\left\\{y=-\frac{1}{6}\right\\}$ and $\left\\{y=-\frac{1}{3}\right\\}$. The
third cavity is filled with a two-layered medium, characterized by wavenumbers
$\kappa=1+0.5{\text{i}}$ and $\kappa=0.5$, and it is separated at the middle
of the cavity. In order to evaluate the accuracy of our proposed method, we
once again conduct a comparison with the results attained through the adaptive
finite element TBC method [31]. The comparison is presented in Figures 5 and
6. In both figures, the red solid lines and blue circle points depict the
magnitudes of the total electric field on the diagonal of the left, middle,
and right cavities, respectively, obtained using our proposed method and the
adaptive TBC method. The comparison evidently demonstrates that the proposed
method produces accurate outcomes, as they closely align with the results
obtained from the adaptive TBC method. Nevertheless, it is important to note
that the proposed method necessitates solving a substantially smaller system.

Figure 5. Example 6.4: The total field on the diagonal of left, middle, and
right cavity in TM polarization.

Figure 6. Example 6.4: The total field on the diagonal of left, middle, and
right cavity in TE polarization.

##  7\. Conclusion

This paper introduces a highly efficient and accurate numerical method for
addressing electromagnetic scattering problems involving rectangular cavities.
By using the Fourier series expansion, the original boundary value problem is
transformed into one-dimensional ordinary differential equations for the
Fourier coefficients. A connection formula is established to link the Fourier
coefficients in each layer to those on the aperture of the cavity. This
approach enables us to solve the system solely on the aperture and store the
connection formula of the Fourier coefficients, significantly reducing
computational resources.

Furthermore, we propose an alternative TBC on the aperture, which involves
only weakly singular integrals. To handle the singularity of the Hankel
function, we utilize the power series of Bessel functions to deduce a
recursive formula, enhancing the smoothness of the integrand function.
Consequently, high-order Gaussian quadratures can be employed, enhancing the
efficiency of the method.

The numerical results demonstrate the effectiveness of our approach in
accurately solving scattering problems with rectangular cavities. This method
proves to be a valuable tool with diverse practical applications in radar
systems, wireless communications, metamaterials, photonic devices, and more.
Future research could explore extending this method to the three-dimensional
Maxwell’s equations and other geometries and optimizing its efficiency for
even more complex scattering scenarios.

## References

  * [1] H. Ammari, G. Bao, and A. Wood, An integral equation method for the electromagnetic scattering from cavities, Math. Meth. Appl. Sci., 23 (2000), 1057–1072.
  * [2] H. Ammari, G. Bao, and A. Wood, Analysis of the electromagnetic scattering from a cavity, Jpn. J. Indus. Appl. Math., 19 (2001), 301–308.
  * [3] J. Babadjian, E. Bonnetier, and F. Triki, Enhancement of electromagnetic fields caused by interacting subwavelength cavities, Multiscale Model. Simul., 8 (2010), 1383–1418.
  * [4] G. Bao and J. Lai, Optimal shape design of a cavity for radar cross section reduction, SIAM J. Control Optim., 52 (2014), 2122–2140.
  * [5] G. Bao and W. Sun, A fast algorithm for the electromagmetic scattering from a large cavity, SIAM J. Sci. Comput., 27 (2005), 553–574.
  * [6] G. Bao and K. Yun, Stability for the electromagnetic scattering from large cavities, Arch. Ration. Mech. Anal., 3 (2016), 1003–1044.
  * [7] G. Bao, K. Yun, and Z. Zhou, Stability of the scattering from a large electromagnetic cavity in two dimensions, SIAM J. Math. Anal., 44 (2012), 383–404.
  * [8] E. Bonnetier and F. Triki, Asymptotic of the Green function for the diffraction by a perfectly conducting plane perturbed by a sub-wavelength rectangular cavity, Math. Methods Appl. Sci., 33 (2010), 772–798.
  * [9] Y. Chen, X. Jiang, P. Li, and J. Lai, A fast algorithm for the electromagnetic scattering from a large rectangular cavity in three dimensions, J. Comput. Phys., 437 (2021), 110331.
  * [10] D. Colton and R. Kress, Inverse Acoustic and Electromagnetic Scattering Theory, Springer, New York, 2013.
  * [11] Y. Chen, P. Li, and X. Yuan, An adaptive finite element PML method for the open cavity scattering problems, Commun. Comput. Phys., 29 (2021), 1505–1540.
  * [12] K. Du, A composite preconditioner for the electromagnetic scattering from a large cavity, J. Comput. Phys., 230 (2011), 8089–8108.
  * [13] K. Du, W. Sun, and X. Zhang, Arbitrary high-order $C^{0}$ tensor product Galerkin finite element methods for the electromagnetic from a large cavity, J. Comput. Phys., 242 (2013), 181–195.
  * [14] Y. Gao, P. Li, and X. Yuan, Electromagnetic field enhancement in a subwavelength rectangular open cavity, Photons, Phonons, and PDEs in Partial Differential Equations and Applications, ed. by H. Ammari, H. Kang, and H. Zhang, Springer, 2:55 (2021), 1–51.
  * [15] J.-M. Jin, Electromagnetic scattering from large, deep, and arbitrarily-shaped open cavities, Electromagn., 1 (1998), 3–34.
  * [16] J.-M. Jin, The Finite Element Method in Electromagnetics , Wiley & Son, New York, 2002.
  * [17] J.-M. Jin and J. L. Volakis, A hybrid finite element method for scattering and radiation by micro strip patch antennas and arrays residing in a cavity, IEEE Trans. Antennas Propag., 39 (1991), 1598–1604.
  * [18] J. Lai, S. Ambikasaran, and L. Greengard, A fast direct solver for high frequency scattering from a large cavity in two dimensions, SIAM J. Sci. Comput., 36 (2014), B887–B903.
  * [19] H. Li, H. Ma, and W. Sun, Legendre spectral Galerkin method for electromagnetic scattering from large cavities, SIAM J. Numer. Anal., 51 (2013), 253–276.
  * [20] P. Li, A survey of open cavity scattering problems, J. Comp. Math., 36 (2018), 1–16.
  * [21] P. Li and A. Wood, A two-dimensional Helmholtz equation solution for the multiple cavity scattering problem, J. Comput. Phys., 240 (2013), 100–120.
  * [22] J. Lin and F. Reitich, Electromagnetic field enhancement in small gaps: a rigorous mathematical theory, SIAM J. Appl. Math., 75 (2015), 2290–2310.
  * [23] J. Lin and H. Zhang, Scattering and field enhancement of a perfect conducting narrow slit, SIAM J. Appl. Math., 77 (2017), 951–976.
  * [24] J. Liu and J.-M. Jin, A special higher order finite-element method for scattering by deep cavities, IEEE Trans. Antennas Propag., 5 (2000), 694–703.
  * [25] F. Olver, D. Lozier, R. Biosvert, and C. Clark, NIST Handbook of Mathematical Functions, Cambridge University Press, New York, 2010.
  * [26] J. Le Perchec, P. Quémerais, A. Barbara, and T. López-Ríos, Controlling strong electromagnetic fields at a subwavelength scale, Phys. Rev. Lett., 3 (2006), 1–5.
  * [27] T. Van and A. Wood, Finite element analysis for 2-D cavity problem, IEEE Trans. Antennas Propag., 51 (2003), 1–8.
  * [28] Y. Wang, K. Du, and W. Sun, A second-order method for the electromagnetic scattering from a large cavity, Numer. Math. Theor. Meth. Appl., 1 (2008), 357–382.
  * [29] J. Wu, Y. Wang, W. Li, and W. Sun, Toeplitz-type approximations to the Hadamard integral operators and their applications in electromagnetic cavity problems, Appl. Numer. Math., 58 (2008), 101–121.
  * [30] X. Wu and W. Zheng, An adaptive perfectly matched layer method for multiple cavity scattering problems, Commun. Comput. Phys., 19 (2016), 534–558.
  * [31] X. Yuan, G. Bao, and P. Li, An adaptive finite element DtN method for the open cavity scattering problems, CSIAM Trans. Appl. Math., 1 (2020), 316–345.
  * [32] M. Zhao, Z. Qiao, and T. Tang, A fast high order method for electromagnetic scattering by large open cavities, J. Comput. Math., 29(2011), 287–304.
  * [33] M. Zhao and N. Zhu, A fast preconditioned iterative method for the electromagnetic scattering by multiple cavities with high wave numbers, J. Comput. Phys., 398 (2019), 108826.