Mathematical analysis of subwavelength resonant acoustic scattering in multi-layered high-contrast structures

We consider the subwavelength resonant phenomenon for the multi-layered high-contrast (MLHC) metamaterials. Assume that the MLHC metamaterials consist of two types of materials: the host matrix material and the high-contrast material. We write

$$D=\cup_{j=1}^{N}D_{j}$$

to denote the entire resonator-nested, where $D_{j}$ is the bounded doubly-connected domain lying between the interior boundary $\Gamma_{j}^{-}$ and the exterior boundary $\Gamma_{j}^{+}$, $j=1,2,\ldots,N$. Each $\Gamma_{j}^{-}$ surrounds $\Gamma_{j+1}^{+}$ and $\Gamma_{j}^{+}$ surrounds $\Gamma_{j}^{-}$ $(j=1,2,\ldots,N-1)$. Denote $D_{j}^{\prime}$ the region of the gap between the $\Gamma_{j}^{-}$ and $\Gamma_{j+1}^{+}$, $j=1,2,\ldots,N-1$. Let $D_{0}^{\prime}$ and $D_{N}^{\prime}$ be the unbounded domain with boundary $\Gamma_{1}^{+}$ and the bounded domain with boundary $\Gamma_{N}^{-}$, respectively. Then, the host matrix material can be written by

$$\mathbb{R}^{3}\setminus{D}=\cup_{j=0}^{N}D^{\prime}_{j}.$$

The configuration of the considered metamaterial is characterized by the density $\rho(x)$ and the bulk modulus $\kappa(x)$ which are given by

$$\rho(x)=\left\{\begin{array}[]{ll}\rho_{\mathrm{r}},&x\in D_{j},\;j=1,2,\ldots,N,\\ \rho,&x\in D_{j}^{\prime},\;j=0,1,\ldots,N,\end{array}\right.\;\;\mbox{ and }\;\;\kappa(x)=\left\{\begin{array}[]{ll}\kappa_{\mathrm{r}},&x\in D_{j},\;j=1,2,\ldots,N,\\ \kappa,&x\in D_{j}^{\prime},\;j=0,1,\ldots,N.\end{array}\right.$$

(2.1)

The wave speeds and wavenumbers of resonators and host matrix materials,respectively, are given by

$$v_{\mathrm{r}}=\sqrt{\frac{\kappa_{\mathrm{r}}}{\rho_{\mathrm{r}}}},\quad v=\sqrt{\frac{\kappa}{\rho}},\quad k_{\mathrm{r}}=\frac{\omega}{v_{\mathrm{r}}},\quad k=\frac{\omega}{v}.$$

(2.2)

We also introduce two dimensionless contrast parameters:

$$\delta=\frac{\rho_{\mathrm{r}}}{\rho},\,\,\tau=\frac{k_{\mathrm{r}}}{k}=\frac{v}{v_{\mathrm{r}}}=\sqrt{\frac{\rho_{\mathrm{r}}\kappa}{\rho\kappa_{\mathrm{r}}}}.$$

(2.3)

We will assume that $v=O(1)$, $v_{\mathrm{r}}=O(1)$, and $\tau=O(1)$; meanwhile $\delta\ll 1.$ This high-contrast assumption is the cause of the underlying system’s subwavelength resonant response and will be at the center of our subsequent analysis.

We will consider the scattering of a time-harmonic acoustic wave by this MLHC structure. This is described by the Helmholtz system of equations

$$\begin{cases}\displaystyle\Delta u+\frac{\omega^{2}}{v^{2}}u=0&\text{in }\mathbb{R}^{3}\setminus D,\\ \displaystyle\Delta u+\frac{\omega^{2}}{v_{\mathrm{r}}^{2}}u=0&\text{in }D,\\ \vskip 3.0pt plus 1.0pt minus 1.0pt\cr\displaystyle u|_{+}-u|_{-}=0&\text{on }\Gamma^{\pm}_{j},\;j=1,2,\ldots,N,\\ \vskip 3.0pt plus 1.0pt minus 1.0pt\cr\displaystyle\delta\frac{\partial u}{\partial\nu}|_{+}=\frac{\partial u}{\partial\nu}|_{-}&\text{on }\Gamma^{+}_{j},\;j=1,2,\ldots,N,\\ \vskip 3.0pt plus 1.0pt minus 1.0pt\cr\displaystyle\frac{\partial u}{\partial\nu}|_{+}=\delta\frac{\partial u}{\partial\nu}|_{-}&\text{on }\Gamma^{-}_{j},\;j=1,2,\ldots,N,\\ \vskip 3.0pt plus 1.0pt minus 1.0pt\cr\displaystyle u^{s}:=u-u^{in}&\mbox{satisfies the Sommerfeld radiation condition,}\end{cases}$$

(2.4)

where $u^{in}$ is the incoming wave, and the notation $\nu$ denotes the outward normal on $\Gamma_{j}^{\pm}$. By the Sommerfeld radiation condition, the scattered wave $u^{s}$ satisfies

$$\left(\frac{\partial}{\partial|x|}-\mathrm{i}k_{0}\right)u^{s}=O(|x|^{-2})\quad\mbox{ as }|x|\to+\infty.$$

(2.5)

In this subsection, we briefly introduce the boundary layer potential operators and then establish the representation formula of the solution (2.4), and we also refer to [15, 29] for more relevant discussions on layer potential techniques.

Let $G_{k}$ be the outgoing fundamental solution to the PDO $\Delta+k^{2}$ in $\mathbb{R}^{3}$, which is given by

$$G_{k}(x)=-\frac{e^{\mathrm{i}k|x|}}{4\pi|x|}.$$

(2.6)

Let $\Omega$ be a bounded domain with a $C^{1,\eta}\,(0<\eta<1)$ boundary $\Gamma$. The single layer potential $\mathcal{S}_{\Gamma}^{k}$ associated with wavenumber $k$ can defined by

$$\mathcal{S}_{{\Gamma}}^{k}[\phi](x)=\int_{{\Gamma}}G_{k}(x-y)\phi(y)~{}\mathrm{d}\sigma(y),\quad x\in\mathbb{R}^{3},$$

(2.7)

where $\phi\in L^{2}(\Gamma)$ is the density function. There hold the following jump relations on the surface [18]

$$\frac{\partial}{\partial\nu}\mathcal{S}^{k}_{\Gamma}[\phi]\Big{|}_{\pm}=\left(\pm\frac{1}{2}I+\mathcal{K}_{{\Gamma}}^{k,*}\right)[\phi]\quad\mbox{on }{\Gamma},$$

(2.8)

where the subscripts $+$ and $-$ denote evaluation from outside and inside the boundary $\Gamma$, respectively, and $\mathcal{K}_{\Gamma}^{k,*}$ is the Neumann-Poincaré (NP) operator defined by

$$\mathcal{K}_{\Gamma}^{k,*}[\phi](x)=\mbox{p.v.}\;\int_{{\Gamma}}\frac{\partial G_{k}(x-y)}{\partial\nu_{x}}\phi(y)~{}\mathrm{d}\sigma(y),\quad x\in\Gamma.$$

Here p.v. stands for the Cauchy principle value. In what follows, we denote by $\mathcal{S}_{\Gamma}$ and $\mathcal{K}_{\Gamma}^{*}$ be the single-layer and Neumann-Poincaré operators $\mathcal{S}_{\Gamma}^{k}$ and $\mathcal{K}_{\Gamma}^{k,*}$, by formally taking $k=0$ respectively.

With the help of the layer potentials in subsection 2.2, the solution to the Helmholtz system (2.4) can be represented by (cf. [32])

$$u(x)=\begin{cases}\displaystyle u^{in}+\mathcal{S}_{\Gamma_{1}^{+}}^{k}[\psi_{1}^{+}](x),&\quad x\in D_{0}^{\prime},\\ \vskip 3.0pt plus 1.0pt minus 1.0pt\cr\displaystyle\mathcal{S}_{{\Gamma_{j}^{+}}}^{k_{\mathrm{r}}}[\phi_{j}^{+}](x)+\mathcal{S}_{{\Gamma_{j}^{-}}}^{k_{\mathrm{r}}}[\psi_{j}^{-}](x),&\quad x\in{D_{j}},\;j=1,2,\ldots,N,\\ \vskip 3.0pt plus 1.0pt minus 1.0pt\cr\displaystyle\mathcal{S}_{{\Gamma_{j}^{-}}}^{k}[\phi_{j}^{-}](x)+\mathcal{S}_{{\Gamma_{j+1}^{+}}}^{k}[\psi_{j+1}^{+}](x),&\quad x\in{D_{j}^{\prime}},\;j=1,2,\ldots,N-1,\\ \vskip 3.0pt plus 1.0pt minus 1.0pt\cr\displaystyle\mathcal{S}_{{\Gamma_{N}^{-}}}^{k}[\phi_{N}^{-}](x),&\quad x\in{D^{\prime}_{N}},\end{cases}$$

(2.9)

where $\psi_{j}^{\pm},\phi_{j}^{\pm}\in L^{2}(\Gamma_{j}^{\pm})$, $j=1,2,\ldots,N.$ Using the transmission conditions in (2.4), and the jump relations for the single layer potentials, we can obtain that $\psi_{j}^{\pm}$ and $\phi_{j}^{\pm}$, $j=1,2,\ldots,N,$ satisfy the following system of boundary integral equations:

$$\mathcal{A}(\omega,\delta)[\Psi]=(u^{in},\delta\frac{\partial u^{in}}{\partial\nu},0,0,\ldots,0)^{T}.$$

(2.10)

Here

$$\Psi=(\psi_{1}^{+},\phi_{1}^{+},\psi_{1}^{-},\phi_{1}^{-},\psi_{2}^{+},\phi_{2}^{+},\psi_{2}^{-},\phi_{2}^{-},\ldots,\psi_{N}^{+},\phi_{N}^{+},\psi_{N}^{-},\phi_{N}^{-})^{T},$$

and the $4N$-by-$4N$ matrix type operator $\mathcal{A}(\omega,\delta)$ has the block tridiagonal form

$$\begin{split}\mathcal{A}(\omega,\delta):&=\begin{pmatrix}\mathcal{M}_{\Gamma_{1}^{+}}&\mathcal{R}_{\Gamma_{1}^{+},\Gamma_{1}^{-}}&&&&\\ \mathcal{L}_{\Gamma_{1}^{-},\Gamma_{1}^{+}}&\mathcal{M}_{\Gamma_{1}^{-}}&\mathcal{R}_{\Gamma_{1}^{-},\Gamma_{2}^{+}}&&&\\ &\mathcal{L}_{\Gamma_{2}^{+},\Gamma_{1}^{-}}&\mathcal{M}_{\Gamma_{2}^{+}}&\mathcal{R}_{\Gamma_{2}^{+},\Gamma_{2}^{-}}&&\\ &&\mathcal{L}_{\Gamma_{2}^{-},\Gamma_{2}^{+}}&\mathcal{M}_{\Gamma_{2}^{-}}&\mathcal{R}_{\Gamma_{2}^{-},\Gamma_{3}^{+}}&\\ &&&\ddots&\ddots&\ddots&\\ &&&&\mathcal{L}_{\Gamma_{N}^{+},\Gamma_{N-1}^{-}}&\mathcal{M}_{\Gamma_{N}^{+}}&\mathcal{R}_{\Gamma_{N}^{+},\Gamma_{N}^{-}}\\ &&&&&\mathcal{L}_{\Gamma_{N}^{-},\Gamma_{N}^{+}}&\mathcal{M}_{\Gamma_{N}^{-}}\end{pmatrix}\end{split}$$

(2.11)

where $\mathcal{M}_{\Gamma_{i}^{\pm}}$ is the self-interaction of the exterior and the interior boundaries for the $i$-th resonator, respectively, defined by

$$\quad\mathcal{M}_{\Gamma_{i}^{+}}=\begin{pmatrix}-\mathcal{S}^{k}_{\Gamma_{i}^{+}}&\mathcal{S}^{k_{\mathrm{r}}}_{\Gamma_{i}^{+}}\\ -\delta(\frac{1}{2}I+\mathcal{K}_{\Gamma_{i}^{+}}^{k,*})&-\frac{1}{2}I+\mathcal{K}_{\Gamma_{i}^{+}}^{k_{\mathrm{r}},*}\end{pmatrix},\quad\mathcal{M}_{\Gamma_{i}^{-}}=\begin{pmatrix}-\mathcal{S}^{k_{\mathrm{r}}}_{\Gamma_{i}^{-}}&\mathcal{S}^{k}_{\Gamma_{i}^{-}}\\ -(\frac{1}{2}I+\mathcal{K}_{\Gamma_{i}^{-}}^{k_{\mathrm{r}},*})&\delta(-\frac{1}{2}I+\mathcal{K}_{\Gamma_{i}^{-}}^{k,*})\end{pmatrix},$$

(2.12)

$\mathcal{L}_{\Gamma_{i}^{-},\Gamma_{i}^{+}}$ and $\mathcal{R}_{\Gamma_{i}^{-},\Gamma_{i+1}^{+}}$ encode the effect of the exterior boundaries of the $i$-th resonator and the $(i+1)$-th resonator, respectively, on the interior boundary of the $i$-th resonator, as defined by

$$\mathcal{L}_{\Gamma_{i}^{-},\Gamma_{i}^{+}}=\begin{pmatrix}\displaystyle 0&-\mathcal{S}^{k_{\mathrm{r}}}_{\Gamma_{i}^{-},\Gamma_{i}^{+}}\\ \vskip 3.0pt plus 1.0pt minus 1.0pt\cr\displaystyle 0&-\mathcal{K}_{\Gamma_{i}^{-},\Gamma_{i}^{+}}^{k_{\mathrm{r}},*}\end{pmatrix},\quad\mathcal{R}_{\Gamma_{i}^{-},\Gamma_{i+1}^{+}}=\begin{pmatrix}\displaystyle\mathcal{S}^{k}_{\Gamma_{i}^{-},\Gamma_{i+1}^{+}}&0\\ \vskip 3.0pt plus 1.0pt minus 1.0pt\cr\displaystyle\delta\mathcal{K}_{\Gamma_{i}^{-},\Gamma_{i+1}^{+}}^{k,*}&0\end{pmatrix},$$

(2.13)

$\mathcal{L}_{\Gamma_{i}^{+},\Gamma_{i-1}^{-}}$ and $\mathcal{R}_{\Gamma_{i}^{+},\Gamma_{i}^{-}}$ encode the effect of the interior boundaries of the $(i-1)$-th resonator and the $i$-th resonator, respectively, on the exterior boundary of the $i$-th resonator, as defined by

$$\mathcal{L}_{\Gamma_{i}^{+},\Gamma_{i-1}^{-}}=\begin{pmatrix}0&-\mathcal{S}^{k}_{\Gamma_{i}^{+},\Gamma_{i-1}^{-}}\\ \vskip 3.0pt plus 1.0pt minus 1.0pt\cr 0&-\delta\mathcal{K}_{\Gamma_{i}^{+},\Gamma_{i-1}^{-}}^{k,*}\end{pmatrix},\quad\mathcal{R}_{{\Gamma_{i}^{+},\Gamma_{i}^{-}}}=\begin{pmatrix}\displaystyle\mathcal{S}^{k_{\mathrm{r}}}_{\Gamma_{i}^{+},\Gamma_{i}^{-}}&0\\ \vskip 3.0pt plus 1.0pt minus 1.0pt\cr\displaystyle\mathcal{K}_{\Gamma_{i}^{+},\Gamma_{i}^{-}}^{k_{\mathrm{r}},*}&0\end{pmatrix}.$$

(2.14)

Here, we introduced the operators $\mathcal{S}^{k}_{\Gamma_{i},\Gamma_{j}}:L^{2}(\Gamma_{j})\to L^{2}(\Gamma_{i})$ and $\mathcal{K}_{\Gamma_{i},\Gamma_{j}}^{k,*}:L^{2}(\Gamma_{j})\to L^{2}(\Gamma_{i})$ as defined respectively by

$$\mathcal{S}^{k}_{\Gamma_{i},\Gamma_{j}}[\varphi]=\mathcal{S}^{k}_{\Gamma_{j}}[\varphi]\big{|}_{\Gamma_{i}}\;\mbox{ and }\;\mathcal{K}_{\Gamma_{i},\Gamma_{j}}^{k,*}[\varphi]=\frac{\partial}{\partial\nu_{i}}\mathcal{S}^{k}_{\Gamma_{j}}[\varphi]\big{|}_{\Gamma_{i}},\mbox{ for }\forall\varphi\in L^{2}(\Gamma_{j}).$$

In order to understand the resonant behavior of the $N$-layer nested scatterers, we shall give the definition of the subwavelength resonant frequencies and resonant modes of the system based on the high contrast $\delta$ between the materials.

When the material parameters are real, it is easy to see that $\overline{\mathcal{A}(\omega,\delta)}=\mathcal{A}(-\overline{\omega},\delta)$, from which we can see that the subwavelength resonant frequencies will be symmetric about the imaginary axis.

Based on Lemma 2.1, we will subsequently state results only for the subwavelength resonant frequencies with positive real parts. The existence of subwavelength resonant frequencies is given by the following Theorem, which was proved in [32] by using the generalized Rouché theorem [15].

In this subsection, we recall the first characterization of resonances based on spherical wave expansions. By rather lengthy and tedious calculations, the exact formulas for the eigenfrequencies of single-resonator and dual-resonator, as well as numerical computations for any finite multi-layer nested resonators, have been derived and presented in [32].

Let $j_{n}(t)$ and $h^{(1)}_{n}(t)$ be the spherical Bessel and Hankel functions of the first kind of order $n$, respectively, and $Y_{n}(\hat{x})$ be the spherical harmonics. By using spherical coordinates, the solution $u$ to (2.4), with the material parameters given in (2.1)–(2.3), can be written as

$$u(x)=\begin{cases}\displaystyle u^{in}+\sum_{n=0}^{+\infty}a^{+}_{1,n}h^{(1)}_{n}(kr)Y_{n},&\quad x\in D_{0}^{\prime},\\ \vskip 3.0pt plus 1.0pt minus 1.0pt\cr\displaystyle\sum_{n=0}^{+\infty}\left(b^{+}_{j,n}j_{n}(k_{\mathrm{r}}r)Y_{n}+a^{-}_{j,n}h^{(1)}_{n}(k_{\mathrm{r}}r)Y_{n}\right),&\quad x\in{D_{j}},\;j=1,2,\ldots,N,\\ \vskip 3.0pt plus 1.0pt minus 1.0pt\cr\displaystyle\sum_{n=0}^{+\infty}\left(b^{-}_{j,n}j_{n}(kr)Y_{n}+a^{+}_{j+1,n}h^{(1)}_{n}(kr)Y_{n}\right),&\quad x\in{D_{j}^{\prime}},\;j=1,2,\ldots,N,\end{cases}$$

(3.4)

where $a^{+}_{N+1,n}=0$. By using the transmission conditions across $\Gamma^{\pm}_{j}$, $j=1,2,\ldots,N$, we find that the constants satisfy

$$\bm{A}_{(n)}(\omega,\delta)[\bm{a}^{\pm}]=(u^{in}|_{\Gamma_{1}^{+}},\delta\frac{\partial u^{in}}{\partial\nu}|_{\Gamma_{1}^{+}},0,0,\cdots,0)^{T},$$

for all $n\in\mathbb{N}\cup\{0\}$. Here

$$\bm{a}^{\pm}=\left(a^{+}_{1,n},b^{+}_{1,n},a^{-}_{1,n},b^{-}_{1,n},\cdots,a^{+}_{N,n},b^{+}_{N,n},a^{-}_{N,n},b^{-}_{N,n}\right)^{T}$$

and the $4N$-by-$4N$ matrix type operator $\bm{A}_{(n)}(\omega,\delta)$ has the block tridiagonal form

$$\bm{A}_{(n)}(\omega,\delta):=\begin{pmatrix}{M}^{+}_{1,n}&{R}^{+,-}_{1,1,n}&&&&&\\ \vskip 3.0pt plus 1.0pt minus 1.0pt\cr{L}^{-,+}_{1,1,n}&{M}^{-}_{1,n}&{R}^{-,+}_{1,2,n}&&&&\\ \vskip 3.0pt plus 1.0pt minus 1.0pt\cr&{L}^{+,-}_{2,1,n}&{M}^{+}_{2,n}&{R}^{+,-}_{2,2,n}&&&\\ \vskip 3.0pt plus 1.0pt minus 1.0pt\cr&&{L}^{-,+}_{2,2,n}&{M}^{-}_{2,n}&{R}^{-,+}_{2,3,n}&&&\\ &&&\ddots&\ddots&\ddots&\\ &&&&{L}^{+,-}_{N,N-1,n}&{M}^{+}_{N,n}&{R}^{+,-}_{N,N,n}\\ \vskip 3.0pt plus 1.0pt minus 1.0pt\cr&&&&&{L}^{-,+}_{N,N,n}&{M}^{-}_{N,n}\end{pmatrix},$$

(3.5)

where

$${M}_{i,n}^{+}=\begin{pmatrix}-h^{(1)}_{n}(kr_{i}^{+})&j_{n}(k_{\mathrm{r}}r_{i}^{+})\\ \vskip 3.0pt plus 1.0pt minus 1.0pt\cr-\delta h^{(1)^{\prime}}_{n}(kr_{i}^{+})&\tau j^{\prime}_{n}(k_{\mathrm{r}}r_{i}^{+})\end{pmatrix},\quad{M}_{i,n}^{-}=\begin{pmatrix}-h^{(1)}_{n}(k_{\mathrm{r}}r_{i}^{-})&j_{n}(kr_{i}^{-})\\ \vskip 3.0pt plus 1.0pt minus 1.0pt\cr-\tau h^{(1)^{\prime}}_{n}(k_{\mathrm{r}}r_{i}^{-})&\delta j^{\prime}_{n}(kr_{i}^{-})\end{pmatrix},$$

(3.6)

$${L}^{-,+}_{i,i,n}=\begin{pmatrix}0&-j_{n}(k_{\mathrm{r}}r_{i}^{-})\\ \vskip 3.0pt plus 1.0pt minus 1.0pt\cr 0&-\tau j^{\prime}_{n}(k_{\mathrm{r}}r_{i}^{-})\end{pmatrix},\quad{R}^{+,-}_{i,i,n}=\begin{pmatrix}h^{(1)}_{n}(k_{\mathrm{r}}r_{i}^{+})&0\\ \vskip 3.0pt plus 1.0pt minus 1.0pt\cr\tau h^{(1)^{\prime}}_{n}(k_{\mathrm{r}}r_{i}^{+})&0\end{pmatrix},$$

and

$${L}^{+,-}_{i,i-1,n}=\begin{pmatrix}0&-j_{n}(kr_{i}^{+})\\ \vskip 3.0pt plus 1.0pt minus 1.0pt\cr 0&-\delta j^{\prime}_{n}(kr_{i}^{+})\end{pmatrix},\quad{R}^{-,+}_{i,i+1,n}=\begin{pmatrix}h^{(1)}_{n}(kr_{i}^{-})&0\\ \vskip 3.0pt plus 1.0pt minus 1.0pt\cr\delta h^{(1)^{\prime}}_{n}(kr_{i}^{-})&0\end{pmatrix}.$$

It is important to emphasize that, from a physical perspective, we are concerned with the resonance of nested materials, which corresponds to the system’s lowest resonant frequency. At this frequency, the system exhibits a factor corresponding to $n=0$, as the lowest resonance is characterized by the fewest number of oscillations [32, 13]. It can be seen in the proofs of [32, Theorems 4.1-4.2] that the primary reason for mode splitting lies in the fact that as the number of nested resonators increases, the degree of the corresponding characteristic polynomial also increases, while the type of resonance (which consists solely of monopolar resonances) remains unchanged. The mathematical justification for the monopole approximation is given in Section 4. Consequently, the $4N$-by-$4N$ matrix

$$\bm{A}(\omega,\delta):=\bm{A}_{(0)}(\omega,\delta)$$

(3.7)

becomes singular.

Although Proposition 3.1 provides a complete characterization of the solution to (2.4), it does not directly facilitate the prediction of scattering resonances. This is primarily due to the computational complexity involved in deriving the asymptotic expansions of the determinant $\det(\bm{A}(\omega,\delta))$ of the $4N$-by-$4N$ matrix $\bm{A}(\omega,\delta)$ by using the following asymptotic expansions:

$$j_{0}(t)=\frac{\sin t}{t}=1-\frac{t^{2}}{6}+\frac{t^{4}}{120}+O(t^{6}),$$

(3.8)

$$h^{(1)}_{0}(t)=\frac{\sin t}{t}-\mathrm{i}\frac{\cos t}{t}=1-\frac{t^{2}}{6}+\frac{t^{4}}{120}+\mathrm{i}(-\frac{1}{t}+\frac{t}{2}-\frac{t^{3}}{24})+O(t^{5}),$$

(3.9)

for $t\ll 1$. We remark in Subsection 3.2 that the scattering problem (2.4) can be reformulated as a $2N\times 2N$ linear system based on the Dirichlet-to-Neumann (DtN) map, which gives a second characterization of the resonances.

In this section, we characterize the DtN map of the Helmholtz problem (2.4) in MLHC metamaterial $D=\cup_{j=1}^{N}D_{j}$, with the aim of applying the DtN approach developed in [41] to analyze subwavelength resonances.

The following lemma provides an explicit expression for the solution to exterior Dirichlet problems on $\mathbb{R}^{3}\setminus D$.

We give the exact matrix representation of the DtN map $\mathcal{T}^{k}$ in the following proposition.