\begin{scope}[xshift = 1.3cm]
\draw (-0.2,0) circle (8pt);
\end{scope};
\end{scope};

\begin{scope}[xshift=1.85cm]
\draw[<->, opacity=0.5] (0.2,0) -- (1.1,0) node[pos=0.5, yshift=-7pt,]{$l$};
\draw (0.2,0) circle (8pt);
%\draw[fill] (0.2,0) circle (0.5pt);
\begin{scope}[xshift = 1.3cm]
\draw (-0.2,0) circle (8pt);
\end{scope};
\end{scope};

\begin{scope}[xshift=4.05cm]
\draw (0.2,0) circle (8pt);
%\draw[<->, opacity=0.5] (0.2,0) -- (1.1,0) node[pos=0.5, yshift=-5pt,]{$d$};
%	\draw[opacity=0.5] (0.65,1) -- (0.65,-1);
\begin{scope}[xshift = 1.3cm]
\draw (-0.2,0) circle (8pt);
\end{scope};
\end{scope};

\begin{scope}[xshift=6.25cm]
\draw (0.2,0) circle (8pt);
%	\draw[opacity=0.5] (0.65,1) -- (0.65,-1) node[above left]{$\phi^z = \pi$};
%\draw[<->, opacity=0.5] (0.2,0) -- (1.1,0) node[pos=0.5, yshift=-5pt,]{$d$};
\begin{scope}[xshift = 1.3cm]
\draw (-0.2,0) circle (8pt);
\end{scope};
\end{scope};

\begin{scope}[xshift=-4.05cm]
\draw (0.2,0) circle (8pt);
%\draw[<->, opacity=0.5] (0.2,0) -- (1.1,0) node[pos=0.5, yshift=-5pt,]{$d$};
\begin{scope}[xshift = 1.3cm]
\draw (-0.2,0) circle (8pt);
\end{scope};
\end{scope};

\begin{scope}[xshift=-6.25cm]
\draw (0.2,0) circle (8pt);
%	\draw[opacity=0.5] (-0.5,1) -- (-0.5,-1);
%	\draw[opacity=0.5]  (1.8,1) -- (1.8,-1) node[above left]{$\phi^z = 0$};
%\draw[<->, opacity=0.5] (0.2,0) -- (1.1,0) node[pos=0.5, yshift=-5pt,]{$d$};
\begin{scope}[xshift = 1.3cm]
\draw (-0.2,0) circle (8pt);
\end{scope};
\end{scope};

%	\begin{scope}[xshift=-6.25cm]
%	\draw [decorate,opacity=0.5,decoration={brace,amplitude=10pt},xshift=-4pt,yshift=0pt]
%	(-0.25,0.5) -- (1.8,0.5) node [black,midway,xshift=-0.6cm]{};
%	\end{scope}

\begin{scope}[xshift=-4.05cm]
\draw [decorate,opacity=0.5,decoration={brace,amplitude=10pt},xshift=-4pt,yshift=0pt]
(-0.25,0.5) -- (1.8,0.5) node [black,midway,xshift=-0.6cm]{};
\end{scope}

\begin{scope}[xshift=-1.85cm]
\draw [decorate,opacity=0.5,decoration={brace,amplitude=10pt},xshift=-4pt,yshift=0pt]
(-0.25,0.5) -- (1.8,0.5) node [black,midway,xshift=-0.6cm]{};
\end{scope}

\begin{scope}[xshift=0.72cm]
\draw [decorate,opacity=0.5,decoration={brace,amplitude=10pt},xshift=-4pt,yshift=0pt]
(-0.25,0.5) -- (1.8,0.5) node [black,midway,xshift=-0.6cm]{};
\end{scope}

\begin{scope}[xshift=2.92cm]
\draw [decorate,opacity=0.5,decoration={brace,amplitude=10pt},xshift=-4pt,yshift=0pt]
(-0.25,0.5) -- (1.8,0.5) node [black,midway,xshift=-0.6cm]{};
\end{scope}

%	\begin{scope}[xshift=5.12cm]
%	\draw [decorate,opacity=0.5,decoration={brace,amplitude=10pt},xshift=-4pt,yshift=0pt]
%	(-0.25,0.5) -- (1.8,0.5) node [black,midway,xshift=-0.6cm]{};
%	\end{scope}

\begin{scope}[xshift=-4.05cm-4pt]
\draw[opacity=0.5]    (0.775,0.9) to[out=90,in=-100] (1.54,1.5);
\end{scope}

\begin{scope}[xshift=-1.85cm-4pt]
\draw[opacity=0.5]    (0.775,0.9) to[out=100,in=-80] (-0.64,1.5);
\end{scope}

%	\begin{scope}[xshift=-1.85cm-4pt]
%	\draw[opacity=0.5]    (0.775,0.9) to[out=100,in=-45] (-2.25,1.5);
%	\end{scope}

\begin{scope}[xshift=0.72cm-4pt]
\draw[opacity=0.5]    (0.775,0.9) to[out=90,in=-100] (2.04,1.5);
\end{scope}

\begin{scope}[xshift=2.92cm-4pt]
\draw[opacity=0.5]    (0.775,0.9) to[out=100,in=-80] (-0.14,1.5);
\end{scope}

\begin{scope}[xshift=-2.4cm]
\draw [decorate,opacity=0.5,decoration={brace,amplitude=10pt,mirror},xshift=-4pt,yshift=0pt]
(-1.125,1.9) -- (0.925,1.9) node [black,midway,xshift=-0.6cm]{};
\begin{scope}[xshift=-4pt]
\draw[opacity=0.5, dotted]    (-1.25,2) -- (-1.25,2.7);
\draw[opacity=0.5, dotted]    (1.05,2) -- (1.05,2.7);
\end{scope}
\begin{scope}[xshift = -4pt-0.75cm]
\draw[opacity=0.5] (0.2,2.3) circle (8pt);
\node[opacity=0.5] at (0.65,2.9) {$\varphi_j^z = 0$};
%	\draw[opacity=0.5] (0.65,1) -- (0.65,-1) node[above left]{$\phi^z = \pi$};
%\draw[<->, opacity=0.5] (0.2,0) -- (1.1,0) node[pos=0.5, yshift=-5pt,]{$d$};
\begin{scope}[xshift = 1.3cm]
\draw[opacity=0.5] (-0.2,2.3) circle (8pt);
\end{scope}
\end{scope}
\end{scope}

\begin{scope}[xshift=2.87cm]
\draw [decorate,opacity=0.5,decoration={brace,amplitude=10pt,mirror},xshift=-4pt,yshift=0pt]
(-1.125,1.9) -- (0.925,1.9) node [black,midway,xshift=-0.6cm]{};
\begin{scope}[xshift=-4pt]
\draw[opacity=0.5, dotted]    (-1.25,2) -- (-1.25,2.7);
\draw[opacity=0.5, dotted]    (1.05,2) -- (1.05,2.7);
\end{scope}
\begin{scope}[xshift = -4pt-0.75cm]
\draw[opacity=0.5] (-0.05,2.3) circle (8pt);
\node[opacity=0.5] at (0.65,2.9) {$\varphi_j^z = \pi$};
%	\draw[opacity=0.5] (0.65,1) -- (0.65,-1) node[above left]{$\phi^z = \pi$};
%\draw[<->, opacity=0.5] (0.2,0) -- (1.1,0) node[pos=0.5, yshift=-5pt,]{$d$};
\begin{scope}[xshift = 1.3cm]
\draw[opacity=0.5] (0.05,2.3) circle (8pt);
\end{scope}
\end{scope}
\end{scope}
\end{tikzpicture}
\caption{Two-dimensional cross-section of a finite dimer chain with 13 resonators, heuristically showing how to identify unit cells with different Zak phases on either side of the edge.} \label{fig:finite_SSH}
\end{figure}

We model wave propagation in the crystal $D$ by the Helmholtz problem
\begin{equation} \label{eq:scattering_finite}
\left\{
\begin{array} {ll}
\ds \Delta {u}+ \omega^2 {u}  = 0 \quad &\text{in } \R^3 \setminus \p D, \\
\nm
\ds  {u}|_{+} -{u}|_{-}  =0  \quad &\text{on } \partial D, \\
\nm
\ds  \delta \frac{\partial {u}}{\partial \nu} \bigg|_{+} - \frac{\partial {u}}{\partial \nu} \bigg|_{-} =0 \quad& \text{on } \partial D, \\
\nm
\ds |x| \left(\tfrac{\p}{\p|x|}-\i \omega\right)u \to 0
&\text{as } {|x|} \rightarrow \infty.
\end{array}
\right.
\end{equation}
\Cref{fig:perturbation1} shows a comparison between the band-gap frequencies of a chain with a topological defect (as in \Cref{fig:finite_SSH}), compared against a chain with a point defect (analogously as in \Cref{sec:defects}).  Not only is the topological band-gap frequency further away from the edges of the band gap, but it also exhibits a lower variance when random Gaussian errors is imposed on the resonator locations. For large error, we see that the band-gap frequency of the topological defect chain (\Cref{fig:nondilute}) is much more robust than the band-gap frequency of the point-defect chain.

In \Cref{rmk:tightbind} we observed that, under the dilute assumption, the capacitance formulation is analogous to the tight-binding model commonly utilized in studies of quantum-mechanical systems. The tight-binding model is often coupled with a nearest-neighbour approximation, whereby long-range interactions are neglected. \Cref{fig:tight_bind} shows the bulk and band-gap frequencies computed using such nearest-neighbour approximation. Compared to \Cref{fig:dilute}, where all interactions are taken into account, we observe that the nearest-neighbour approximation is not accurate. This discrepancy shows a fundamental difference between topological edge modes in the setting of classical waves and quantum-mechanical waves. Moreover, we see that the band-gap frequency in \Cref{fig:tight_bind} is unaffected by the error. This can be seen as a consequence of chiral symmetry, which is not present without the nearest-neighbour approximation.

\begin{figure}[p]
\begin{center}
\begin{subfigure}[b]{0.45\linewidth}
\includegraphics[height=5.0cm]{protected_41_nondilute}
\caption{Dimer chain with 41 resonators, separation distances $d=3, d' = 6$.} \label{fig:nondilute}
\end{subfigure}
\hspace{10pt}
\begin{subfigure}[b]{0.45\linewidth}
\includegraphics[height=4.1cm]{pointdefect_41v2}
\caption{Point-defect chain with 41 resonators, separation distance $d=12$ and defect radius $R_d = 0.99$.} \label{fig:point_defect_errors}
\end{subfigure}
\caption{Simulation of band-gap frequency (red) and bulk frequencies (black) of a topological defect (a) and point defect (b) chain, with Gaussian $\mathcal{N}(0,\sigma^2)$ errors added to the resonator positions. The standard deviation $\sigma$ is expressed as a percentage of the average resonator separation. %In all cases, the resonator radius was $R=1$.
\label{fig:perturbation1}
\end{center}
\end{figure}

\begin{figure}[p]
\begin{center}
\begin{subfigure}[b]{0.45\linewidth}
\includegraphics[height=5.0cm]{protected_41}
\caption{Dilute dimer chain with 41 resonators, separation distances $d=12, d' = 42$.} \label{fig:dilute}
\end{subfigure}
\hspace{10pt}
\begin{subfigure}[b]{0.45\linewidth}
\includegraphics[height=5.0cm]{tightbinding}
\caption{Nearest-neighbour approximation for the dilute dimer chain from (a).% with 41 resonators, separation distances $d=12, d' = 42$.
} \label{fig:tight_bind}
\end{subfigure}
\caption{Similar simulations as in \Cref{fig:perturbation1}, but for a dilute chain. (a): fully-coupled simulations, where all interactions are taken into account. (b): nearest-neighbour approximation, where long-range interactions are neglected. The nearest-neighbour approximation provides a chirally symmetric problem, in which the centre frequency is preserved. This approximation is not accurate, and the centre frequency will have a small but non-zero variance corresponding to small but non-zero long-range interactions.
\label{fig:perturbation2}
\end{center}
\end{figure}

\subsubsection{Edge modes in a dislocated chain}
In \Cref{sec:bandinv} we saw that the structure under consideration has different topological properties in the different cases $l>l'$ and $l<l'$. These two regimes only differ by a choice of unit cell: if we shift the unit cell $Y$ by $L/2$ the regimes will be swapped. Correspondingly, edges between the different topological regimes are created by introducing a shift between the unit cells on either side of the edge. This is exemplified in \Cref{fig:finite_SSH}. Another example of this principle is to translate, or \emph{dislocate}, half the chain by some distance $d$ as illustrated in \Cref{fig:infinite}. In this case, we consider a chain of resonators given by
$$\Dc_\mathrm{disloc} = \Bigg(\bigcup_{m=-\infty}^{-1} D+(mL,0,0)\Bigg)\cup\Bigg(\bigcup_{m=0}^\infty D+(mL+d,0,0)\Bigg),$$
where $D$ is the resonator dimer defined as in the beginning of this section.

Qualitatively, it is straightforward to understand how the dislocation will affect the existence and behaviour of band-gap frequencies. As $d$ increases from $0$, we will detune the dimer coupling, which is responsible for the band gap, and band-gap frequencies will therefore appear from each edge of the band gap. On the other hand, when $d$ is very large the two half-space chains will decouple, and the bulk-boundary correspondence suggests that each of these half-chains will support a single band-gap frequency. We therefore expect only a single band-gap frequency in this limit. As $d$ varies between $0$ and $∞$ these two frequencies hybridize, and will together cover the whole band gap. This is sketched in \Cref{fig:gapcrossing}.

We now set out to prove that the picture outlined above is correct. Similarly to \Cref{sec:defects}, we will apply a fictitious source method to model the dislocated structure $_disloc$ in terms of the original structure $$ along with fictitious sources on the dislocated resonators. The analysis is separated into three parts, depending on the value of $d$:
\begin{itemize}
\item $d\ll1$. In this case, we can use asymptotic expansions in terms of $d$ to prove that there is a band-gap frequency emerging from each edge of the band gap.
\item $d=mL$ for some $m>0$. This particular dislocation is equivalent to removing $m$ dimers from $\Dc$. This observation simplifies the problem and allows explicit computations of the band-gap frequencies in terms of the eigenvalue problem of certain Toeplitz matrices.
\item $d>d_0$, where $d_0$ is the width of one resonator. Due to technical reasons, we restrict the remaining values of $d$, so that the dislocated and the original resonators are not overlapping.
\end{itemize}
\begin{figure}
\centering
\begin{tikzpicture}[scale=0.5]
\begin{scope}[shift={(1.5,2)}]
\draw (-6,0) circle (0.3);
\draw (-5,0) circle (0.3);
\draw (-3,0) circle (0.3);
\draw (-2,0) circle (0.3);
\draw (0,0) circle (0.3);
\draw (1,0) circle (0.3);
\draw (3,0) circle (0.3);
\draw (4,0) circle (0.3);
\draw (6,0) circle (0.3);
\draw (7,0) circle (0.3);
\draw (9,0) circle (0.3);
\draw (10,0) circle (0.3);
\draw (12,0) circle (0.3);
\draw (13,0) circle (0.3);
\node at (14.5,0) {$\dots$};
\node at (-7.5,0) {$\dots$};
\draw[dashed,opacity=0.5,red] (3.5,-0.5) -- (3.5,0.6);
\draw[-Latex,opacity=0.5,red] (3.5,-0.6) to[out=-95,in=10] (2.3,-1.6);
\draw[-Latex,opacity=0.5,red] (3.5,-0.6) to[out=-85,in=170] (4.7,-1.6);
\end{scope}
\draw (-6,0) circle (0.3);
\draw (-5,0) circle (0.3);
\draw (-3,0) circle (0.3);
\draw (-2,0) circle (0.3);
\draw (0,0) circle (0.3);
\draw (1,0) circle (0.3);
\draw (3,0) circle (0.3);
\draw[dashed,opacity=0.5,red] (3.5,-0.5) -- (3.5,0.5);
\draw[dashed,opacity=0.5,red] (6.5,-0.5) -- (6.5,0.5);
\draw (7,0) circle (0.3);
\draw (9,0) circle (0.3);
\draw (10,0) circle (0.3);
\draw (12,0) circle (0.3);
\draw (13,0) circle (0.3);
\draw (15,0) circle (0.3);
\draw (16,0) circle (0.3);
\node at (17.5,0) {$\dots$};
\node at (-7.5,0) {$\dots$};
\draw[<->,opacity=0.5,red] (3.5,-0.6) -- (6.5,-0.6) node[pos=0.5,below]{\small $d$};
%	\draw[dotted] (-6.5,-1.5) rectangle (3.5,1.5);
%	\draw[dotted] (6.5,-1.5) rectangle (16.5,1.5);
% 	\draw[<->] (3,0.5) -- (7,0.5) node[pos=0.5, yshift=7pt]{$l+d$};
% 	\draw[<->] (0,0.5) -- (1,0.5) node[pos=0.5, yshift=7pt]{$l$};
% 	\draw[<->] (0,-0.5) -- (3,-0.5) node[pos=0.5, yshift=-7pt]{$L$};
% 	\draw[<->] (9,0.5) -- (10,0.5) node[pos=0.5, yshift=7pt]{$l$};
% 	\draw[<->] (7,-0.5) -- (10,-0.5) node[pos=0.5, yshift=-7pt]{$L$};
%	\draw[<->, opacity=0.5] (0,-0.8) -- (5,-0.8) node[pos=0.5, yshift=-7pt,]{$L$};
\end{tikzpicture}
\caption{We start with an array of pairs of subwavelength resonators, known to have a subwavelength band gap. A dislocation (with size $d>0$) is introduced to create band-gap frequencies.} \label{fig:infinite}
\end{figure}

\begin{figure}[htb]
\centering
\begin{tikzpicture}
\draw[line width=0.5mm,red] plot [smooth] coordinates {(0,0.5) (1,0.8) (5.5,1.05) (7.7,1.1)};
\draw[line width=0.5mm,red] plot [smooth] coordinates {(0,2.1) (1,1.3) (4,1.12) (7.7,1.1)};
\draw[gray!60!white,fill=gray!60!white] (0,0) -- (7.7,0) -- (7.7,0.5) --  (0,0.5);
% 	\draw plot [smooth] coordinates {(7.7,0.5) (4,0.55) (2,0.5) (0,0.5)};
\draw (0,0.5) -- (7.7,0.5);
\draw[gray!60!white,fill=gray!60!white] (0,2.6) -- (7.7,2.6) -- (7.7,2.1) -- (0,2.1);
% 	\draw plot [smooth] coordinates {(7.7,2.1) (4,2.1) (2,2.05) (0,2)} -- (0,2.5);
\draw (7.7,2.1) -- (0,2.1);
\draw[->,line width=0.4mm] (0,0) -- (8,0) node[pos=0.95, yshift=-8pt]{$d$};
\draw[->,line width=0.4mm] (0,0) -- (0,2.85) node[pos=0.85, xshift=-8pt]{$\omega$};
\draw[<-,opacity=0.5] (3.5,0.95) to[out=80,in=180] (4,1.7);
\node[right,opacity=0.5] at (4,1.7) {\small band-gap frequencies};
\draw[<-,opacity=0.5] (3.4,1.15) to[out=80,in=180] (4,1.7);
\draw [decorate,opacity=0.5,decoration={brace,amplitude=10pt}]
(7.8,2.1) -- (7.8,0.5) node [midway,right,xshift=3mm]{\small band gap};
\node at (3.8,2.35){\color{white}\small essential spectrum};
\node at (3.8,0.25){\color{white}\small essential spectrum};
\draw [decorate,opacity=0.5,decoration={brace,amplitude=10pt}]
(-0.6,0) -- (-0.6,2.6) node [midway,xshift=-4mm,align=right,left]{\small subwavelength\\ regime};
\end{tikzpicture}
\caption{As the dislocation size $d$ increases from zero, a band-gap frequency appears from each edge of the subwavelength band gap. These two frequencies converge to a single value within the subwavelength band gap as $d\to\infty$.} \label{fig:gapcrossing}
\end{figure}

With these ideas at hand, the following two theorems were proved in [7]. The first result, valid for small $d$, shows the emergence of a band-gap frequency from each edge of the band gap.
\begin{thm}\label{prop:smalld}
Assume that $D_1$ and $D_2$ are strictly convex. For small enough $d$ and $\delta$, and in the case $l_0 > 1/2$, there are two band-gap frequencies $\omega_1(d), \omega_2(d)$ such that $\omega_j(d) \rightarrow \omega_j^\diamond, j=1,2$ as $d\rightarrow 0$. In the case $l_0 < 1/2$, there are no band-gap frequencies as $d, \delta \rightarrow 0$.
\end{thm}
\begin{thm} \label{thm:main}
Assume that the resonators are in the dilute regime specified by \eqref{eq:dilute_SSH} and that $l_0>1/2$. Then, for small enough $\delta$ and $\epsilon$, there exists some $d_0=O(\epsilon)$ such that there are two band-gap frequencies $\omega_1(d)$ and $\omega_2(d)$
for all $d\in [d_0,\infty)$, both of which converge to the same value $\omega_\infty$ as $d\rightarrow \infty$.
\end{thm}
In particular, \Cref{thm:main} states that the band-gap frequencies will cover an interval $ℐ:=[ω_1(d_0), ω_2(d_0)]$ inside the band gap, and therefore allows us to fine-tune the system to achieve optimal robustness. In \Cref{fig:eigenmodes} we see numerically computed eigenmodes, clearly showing how the edge modes of the ``half system'' hybridize to create the two band-gap frequencies.

\begin{figure}[htb]
\begin{minipage}[b]{0.48\linewidth}
\vspace{0pt}
\centering
\includegraphics[width=\linewidth,trim={1cm 1cm 1cm 0},clip]{30mode}
\end{minipage}
\begin{minipage}[b]{0.48\linewidth}
\centering
\includegraphics[width=\linewidth,trim={1cm 1cm 1cm 5cm},clip]{halfmode}
\vspace{10pt}
\end{minipage}
\caption{\textit{Left:} The two edge modes for an array of 42 spherical resonators of radius $1$. Here, we simulate an array with parameters $L=9$, $l=6$, $d=30$ and $\delta=1/7000$. \textit{Right:} For comparison, the edge mode of the corresponding `half system' is shown, which can be thought of as the $d=\infty$ case. In both cases, the eigenmodes are shown directly above the corresponding system of resonators.} \label{fig:eigenmodes}
\end{figure}

\subsubsection{Non-Hermitian band inversion and edge modes}\label{sec:pt_top}
In the previous sections, we assumed that the wave speeds satisfy $v=v_1=v_2=1$, and only considered topological phenomena originating from the geometry of the resonator chain. In this section we do the opposite, namely we study topologically protected modes in structures where the geometry is periodic, and where topological edges are introduced in terms of the resonator wave speeds. Most importantly, we allow the wave speeds to be complex, giving a non-Hermitian capacitance formulation similarly as in \Cref{sec:exceptional}. Details of this analysis are found in [18].

We let $D$ be the dimer $D=D_1∪D_2$ as defined previously in \Cref{sec:bandinv}. For $i=1,2$ and $m∈$, we introduce the notation
$$D_i^m := D_i + (mL,0,0), \qquad k_i^m := \frac{\omega}{v_i^m}.$$
Here, $v_i^m$ denotes the wave speed inside $D_i^m$. We will assume that the contrast parameter is given by $δ∈$ for all resonators. We then consider
\begin{equation} \label{eq:scattering_pttop}
\left\{
\begin{array} {ll}
\ds \Delta {u}+ k^2 {u}  = 0 & \text{in } \R^3 \setminus \Dc, \\
\nm
\ds \Delta {u}+ (k_i^m)^2 {u}  = 0 & \text{in } D_i^m, \\
\nm
\ds  {u}|_{+} -{u}|_{-}  =0  & \text{on } \p \Dc, \\
\nm
\ds  \delta \frac{\partial {u}}{\partial \nu} \bigg|_{+} - \frac{\partial {u}}{\partial \nu} \bigg|_{-} =0 & \text{on } \p \Dc, \\
\nm
\ds u(x_1,x_2,x_3) & \text{satisfies the outgoing radiation condition as } \sqrt{x_2^2+x_3^2} \rightarrow \infty.
\end{array}
\right.
\end{equation}
In the general setting, this equation is not periodic and cannot be approached using Floquet-Bloch theory. For prescribed values $v_1,v_2∈$ we will consider two different cases, depending on $v_i^m$:
\begin{align}\label{eq:vperiodic}
\text{Periodic structure:} &  \qquad v_1^m = v_1, \ m\in \Z, & &v_2^m = v_2,  \ m\in \Z,
\\ \label{eq:vdefect}
\text{Defect structure:} & \qquad v_1^m = \begin{cases}v_1, & m \leq 0, \\ v_2, &m > 0,\end{cases} & &v_2^m = \begin{cases}v_2, & m \leq 0, \\ v_1, & m > 0.\end{cases}
\end{align}
We emphasize that complex values of $v_i$ correspond to non-Hermitian structures. In the case $v_1 = v_2$, the periodic structure is $$-symmetric.

We begin by studying the periodic structure. In \Cref{sec:bandinv} we described the concept of band inversion in the Hermitian case. In the present case, since the quasiperiodic generalized capacitance matrix $^α$ is non-Hermitian, the band inversion phenomenon is formulated in a slightly different manner. For $j = 1,2,$ we let $𝐮_j$ and $𝐰_j$ denote a bi-orthogonal system of eigenvectors of $^α$. In other words, $𝐮_j$ and $𝐰_j$ are eigenvectors of $^α$ and $(^α)^*$, respectively, and satisfy $𝐰_i·𝐮_j = δ_ij$.

We let the right Bloch eigenmodes $u_j^α$  be the modes of \eqref{eq:scattering_pttop}, and the left Bloch eigenmodes  $w_j^α$ be the modes of \eqref{eq:scattering_pttop} with wave speeds given by $v_1$ and $v_2$, respectively. From \Cref{prop:eigenvector_quasi}, we then have
\begin{align*}
u_j^\alpha &= \mathbf{u}_{j}\cdot  \textbf{S}_D^{\alpha,k}(x) + O(\delta^{1/2}), \\
w_j^\alpha &= \mathbf{w}_{j}\cdot\textbf{S}_D^{\alpha,k}(x) + O(\delta^{1/2}).
\end{align*}

\begin{defn}[non-Hermitian Zak phase] \label{defn:zakNH}
For a non-degenerate band $\omega_j^\alpha$ of the periodic structure \eqref{eq:vperiodic}, we let $u_j^\alpha$ and $w_j^\alpha$ be a family of normalized right, respectively left, eigenmodes which depend continuously on $\alpha$.
We define the (non-Hermitian) Zak phase $\varphi_j^{\mathrm{zak}}$ by
\begin{align*}
\varphi_j^{\mathrm{zak}} &:= \frac{\i}{2} \int_{Y^*}\left( \Big\langle w_j^\alpha, \frac{\p u_j^\alpha  }{\p \alpha}\Big\rangle + \Big\langle u_j^\alpha, \frac{\p  w_j^\alpha}{\p \alpha} \Big\rangle \right) \dx \alpha.
\end{align*}
\end{defn}
In the case $v_1 = v_2∈$, this definition coincides with \Cref{defn:zak}, and we therefore choose the same notation for these two definitions. In the sequel, we will occasionally write $φ_j^zak(v_1, v_2)$ to denote the Zak phase defined with wave speed $v_1$ inside $D_1$ and $v_2$ inside $D_2$.

In the Hermitian case, a non-zero Zak phase is equivalent to an inverted band structure. The fact that the Hermitian Zak phase is quantized originates from the fact that the eigenmodes are purely monopole and dipole modes at $α= 0$ and $α= π/L$. Unlike the Hermitian case, the non-Hermitian Zak phase is not quantized. A non-integer value of the Zak phase can be attributed to a ``partial'' band inversion, where the eigenmodes are expressed as (complex) linear combinations of monopole and dipole modes, which swap as $α$ traverses $Y^*$.

For the next result, we will assume that the Hermitian counterpart of the structure is topologically trivial. More precisely, we assume
\begin{equation}\label{eq:trivial}
\varphi_j^{\mathrm{zak}}\big(\Re(v_1), \Re(v_2)\big) = 0.
\end{equation}
Then the following result from [18] holds.
\begin{prop} \label{prop:zak}
Assume that the chain is periodic and topologically trivial, in other words that it satisfies \eqref{eq:vperiodic} and \eqref{eq:trivial}. Then we have
$$\varphi_j^{\mathrm{zak}}(v_1, v_2) = - \varphi_j^{\mathrm{zak}}(v_2, v_1) + O(\delta) \quad \text{and} \quad \varphi_j^{\mathrm{zak}}(\overline{v_1}, \overline{v_2}) = \varphi_j^{\mathrm{zak}}(v_1, v_2) + O(\delta).$$
In particular, if $v_2 = \overline{v_1}$, we have $\varphi_j^{\mathrm{zak}}(v_1, \overline{v_1}) = O(\delta).$
\end{prop}
\begin{remark}
\Cref{prop:zak} provides intuition on how to create structures supporting edge modes. \Cref{prop:zak} shows that distinct Zak phases can, in general, be achieved by swapping $v_1$ and $v_2$ while keeping the geometry fixed. Therefore, the defect specified in \eqref{eq:vdefect} results in Zak phases with opposite signs on the different sides of the edge. The reason we assume \eqref{eq:trivial} is to emphasize that distinct Zak phases can originate as a pure non-Hermitian effect, which disappears in the Hermitian limit as $\mathrm{Im}(v_1), \mathrm{Im}(v_2) \rightarrow 0$.
\end{remark}

We now turn to the analysis of the defect problem, specified in \eqref{eq:vdefect} and illustrated in \Cref{fig:edge}.
\begin{figure}[tbh]
\centering
\begin{tikzpicture}[scale=1.5]
\pgfmathsetmacro{\rb}{0.25pt}
\pgfmathsetmacro{\rs}{0.2pt}
\coordinate (a) at (0.25,0);
\coordinate (b) at (1.05,0);

\draw (-0.5,0.85) -- (-0.5,-1);
\def\Done{ plot [smooth cycle] coordinates {($(a)+(210:\rb)$) ($(a)+(270:\rs)$) ($(a)+(330:\rb)$) ($(a)+(30:\rs)$) ($(a)+(90:\rb)$) ($(a)+(150:\rs)$)
\def\Dtwo{ plot [smooth cycle] coordinates {($(b)+(30:\rb)$) ($(b)+(90:\rs)$) ($(b)+(150:\rb)$) ($(b)+(210:\rs)$) ($(b)+(270:\rb)$) ($(b)+(330:\rs)$) }};
\draw \Done;
\draw \Dtwo;
\pattern[pattern=north east lines, opacity=0.45, pattern color = red] \Done;
\pattern[pattern=crosshatch dots, opacity=0.45, pattern color = blue] \Dtwo;
\draw (a) node{$v_2$};
\draw (b) node{$v_1$};
\draw (0.65,0.9) node{$m=1$};

\begin{scope}[xshift=-2.3cm]
\coordinate (a) at (0.25,0);
\coordinate (b) at (1.05,0);
\draw[dashed, opacity=0.5] (-0.5,0.85) -- (-0.5,-1);
\draw \Done;
\draw \Dtwo;
\pattern[pattern=north east lines, opacity=0.45, pattern color = red] \Dtwo;
\pattern[pattern=crosshatch dots, opacity=0.45, pattern color = blue] \Done;
\draw (a) node{$v_1$};
\draw (b) node{$v_2$};
\draw (0.65,0.9) node{$m=0$};
\end{scope}
\begin{scope}[xshift=-4.6cm]
\coordinate (a) at (0.25,0);
\coordinate (b) at (1.05,0);
\draw \Done;
\draw \Dtwo;
\pattern[pattern=north east lines, opacity=0.45, pattern color = red] \Dtwo;
\pattern[pattern=crosshatch dots, opacity=0.45, pattern color = blue] \Done;
\draw (a) node{$v_1$};
\draw (b) node{$v_2$};
\draw (0.65,0.9) node{$m=-1$};
\begin{scope}[xshift = 1.2cm]
\draw (-1.6,0) node{$\cdots$};
\end{scope};
\end{scope}
\begin{scope}[xshift=2.3cm]
\coordinate (a) at (0.25,0);
\coordinate (b) at (1.05,0);
\draw[dashed, opacity=0.5] (-0.5,0.85) -- (-0.5,-1);
\draw \Done;
\draw \Dtwo;
\pattern[pattern=north east lines, opacity=0.45, pattern color = red] \Done;
\pattern[pattern=crosshatch dots, opacity=0.45, pattern color = blue] \Dtwo;
\draw (a) node{$v_2$};
\draw (b) node{$v_1$};
\draw (0.65,0.9) node{$m=2$};
\draw (1.7,0) node{$\cdots$};
\end{scope}
\end{tikzpicture}
\caption[Single]{Illustration of the edge. The special case $v_1 = \overline{v_2}$ corresponds to local $\P \T$ symmetry. Legend:
\raisebox{-2pt}{\tikz{
\draw[pattern=crosshatch dots, opacity=0.45, pattern color = blue] (0,0) circle (5pt);
\draw (0,0) circle (5pt);}} material parameter $v_1$,
\raisebox{-2pt}{\tikz{
\draw[pattern=north east lines, opacity=0.45, pattern color = red] (0,0) circle (5pt);
\draw (0,0) circle (5pt);}} material parameter $v_2$. \label{fig:edge} }
\end{figure}
Since the geometry of this defect chain is periodic, we can utilize a Floquet-Bloch approach to derive a capacitance matrix characterization of any localized mode. We define the matrices
$$A = \begin{pmatrix} 1 & b \\[0.3em] b & 1\end{pmatrix}, \quad B = |D_1|\begin{pmatrix} \ds \left(\delta v_2^2\right)^{-1} & \ds b\left(\delta v_1^2\right)^{-1}  \\[0.3em] \ds b\left(\delta v_1^2\right)^{-1}  & \ds \left(\delta v_2^2\right)^{-1} \end{pmatrix},$$
for some parameter $b∈$, which can be interpreted as the decay of the localized mode between two consecutive resonators. The corresponding capacitance matrix, whose eigenvalue problem provides a discrete approximation to the localized mode and corresponding frequency, is given by
$$\C_\mathrm{edge} = B^{-1}C^\alpha A.$$
The eigenvalues of this matrix depend on $α∈Y^*$. In order for localized modes to exist, there must be an eigenvalue $μ= μ_j^α(b_0)$ of $_edge$ which is constant in $α$ for some $b=b_0$.
We can then compute $b$ as $b=b_±$, where
\begin{equation} \label{eq:b}
b_\pm = \frac{1}{2}\left(l\left(1 - \frac{v_1^2}{v_2^2}\right) \pm \sqrt{l^2\left(1 - \frac{v_1^2}{v_2^2}\right)^2 + \frac{4v_1^2}{v_2^2}}\right), \qquad l = \frac{\lambda_2+\lambda_1}{\lambda_2-\lambda_1},
\end{equation}
with $λ_1 = C_11^π/L + C_12^π/L$ and $λ_2 = 2C_11^0$. Depending on the values of $v_1$ and $v_2$, we either have $|b_+| = |b_-| = 1$ (in which case there are no localized modes) or that $|b_±|<1$ while $|b_∓| >1$. Based on these ideas, we can prove the following result [18].
\begin{thm} \label{thm:nonH}
Suppose that the array of resonators has a defect in the material parameters specified by \eqref{eq:vdefect}. Then, for small $\delta$,
\begin{itemize}
\item if $v_1 = \overline{v_2}$ with $\mathrm{Im}(v_1)^2 \leq \frac{\mathrm{Re}(v_1)^2}{l^2-1}$ (unbroken $\P\T$-symmetry), the structure does not support simple localized modes in the subwavelength regime.
\item if $v_1 = \overline{v_2}$ with $\mathrm{Im}(v_1)^2 > \frac{\mathrm{Re}(v_1)^2}{l^2-1}$ (broken $\P\T$-symmetry) or if $v_1 \neq \overline{v_2}$ (no $\P\T$-symmetry), the frequency $\omega$ of a simple localized mode in the subwavelength regime must satisfy
$$\omega = \sqrt{\mu_j^\alpha(b_0)} + O(\delta).$$
Here, $b_0$ is the value of $b$ specified by \eqref{eq:b} satisfying $|b_0| < 1$.
\end{itemize}
\end{thm}
\begin{remark}
\Cref{thm:nonH} provides a characterization of the possible frequency $\omega$ and decay length $b$ of simple localized modes in the subwavelength regime. In order to prove that such mode indeed exists, we would have to prove that one eigenvalue $\mu_j^\alpha(b_0)$ is indeed constant in $\alpha$. Analytically, this is obscured by the fact that the capacitance coefficients have a complicated dependency on $\alpha$. Numerically, however, the eigenvalues $\mu$ are straightforward to compute, providing convincing evidence for such modes to exist. Moreover, the localized modes can easily be computed in large finite chains (see \Cref{fig:edgemodes}), demonstrating excellent agreement between the numerical and theoretically predicted values.
\end{remark}

\begin{figure}[htb]
\begin{subfigure}[b]{0.44\linewidth}
\centering
\includegraphics[width=0.95\linewidth]{modekl}
\caption{Localized mode in the case of small imaginary part.} \label{fig:model}
\end{subfigure}\hfill
\begin{subfigure}[b]{0.44\linewidth}
\centering
\includegraphics[width=0.95\linewidth]{modek}
\caption{Localized mode in the case of large imaginary part.} \label{fig:mode}
\end{subfigure}
\caption{Plots of localized modes in finite, large, arrays of resonators. Observe the different $x$-axis scales in the two subfigures, demonstrating significantly different degrees of localization. The localized modes have purely non-Hermitian origin and disappear in the Hermitian limit when $\Im(v_i)\to0$. }  \label{fig:edgemodes}
\end{figure}

\begin{figure}[h]
\begin{minipage}{0.15\linewidth}
\centering\small  	Hermitian:
\end{minipage}\hfill	\begin{minipage}{0.8\linewidth}
\begin{center}
\begin{tikzpicture}[scale=1.9]
\draw[->] (-0.8,0) -- (0.8,0) node[below]{\footnotesize$\alpha$};
\draw [->] (0,-0.1) -- (0,1) node[left]{\footnotesize$\omega$};

\draw plot[smooth] coordinates {(-0.8,0.7) (-0.55,0.6) (-0.2,0.75) (0.0,0.8) (0.2,0.75)  (0.55,0.6) (0.8,0.7)};
\draw plot[smooth] coordinates {(0.0,0.0) (0.3,0.35)  (0.55,0.48) (0.8,0.4)};
\draw plot[smooth] coordinates {(0.0,0.0) (-0.3,0.35)  (-0.55,0.48) (-0.8,0.4)};
\draw (0,-0.23) node{\footnotesize $l_0 < 1/2$};
\draw (0,-0.4) node{\footnotesize No band inversion};
\begin{scope}[xshift=2.25cm]
\coordinate (d1) at (0.55,0.55);
\coordinate (d2) at (-0.55,0.55);

\draw[->] (-0.8,0) -- (0.8,0) node[below]{\footnotesize$\alpha$};
\draw [->] (0,-0.1) -- (0,1) node[left]{\footnotesize$\omega$};

\draw plot[smooth] coordinates {(-0.8,0.4) (d2) (-0.2,0.75) (0.0,0.8) (0.2,0.75)  (d1) (0.8,0.4)};
\draw plot[smooth] coordinates {(0.0,0.0) (0.3,0.35)  (d1) (0.8,0.7)};
\draw plot[smooth] coordinates {(0.0,0.0) (-0.3,0.35)  (d2) (-0.8,0.7)};
\draw (0,-0.4) node{\footnotesize Dirac cone degeneracy};
\draw (0,-0.23) node{\footnotesize $l_0 = 1/2$};
\end{scope}

\begin{scope}[xshift=4.5cm]
\draw[->] (-0.8,0) -- (0.8,0) node[below]{\footnotesize$\alpha$};
\draw [->] (0,-0.1) -- (0,1) node[left]{\footnotesize$\omega$};

\draw plot[smooth] coordinates {(-0.8,0.7) (-0.55,0.6) (-0.2,0.75) (0.0,0.8) (0.2,0.75)  (0.55,0.6) (0.8,0.7)};
\draw plot[smooth] coordinates {(0.0,0.0) (0.3,0.35)  (0.55,0.48) (0.8,0.4)};
\draw plot[smooth] coordinates {(0.0,0.0) (-0.3,0.35)  (-0.55,0.48) (-0.8,0.4)};

\draw (0,-0.23) node{\footnotesize$l_0 > 1/2$};
\draw (0,-0.4) node{\footnotesize Full band inversion};
\end{scope}
\end{tikzpicture}
\end{center}
\end{minipage}
\begin{minipage}{0.15\linewidth}
\centering	 \small Non-Hermitian:
\end{minipage}	\hfill
\begin{minipage}{0.8\linewidth}
\begin{subfigure}{0.33\linewidth}
\begin{center}
\begin{tikzpicture}
\draw (0,0) node{\includegraphics[clip,width=0.75\linewidth]{bandnpt.png}};
\draw[->] (-1.65,-0.56) -- (1.65,-0.56) node[below]{\footnotesize $\alpha$};
\draw[->] (0,-0.71) -- (0,1.2) node[right]{\footnotesize $\omega$};
\draw (0,-1.4) node{\footnotesize Partial inversion};
\draw (0,-1.1) node{\footnotesize $|v_1| < |v_2|$};
\end{tikzpicture}
\end{center}
\end{subfigure}\begin{subfigure}{0.33\linewidth}
\vspace{4pt}
\begin{center}
\begin{tikzpicture}
\draw (0,0) node{\includegraphics[width=0.75\linewidth]{bandpt.png}};
\draw[->] (-1.65,-0.5) -- (1.65,-0.5) node[below]{\footnotesize $\alpha$};
\draw[->] (0,-0.65) -- (0,1.2) node[right]{\footnotesize $\omega$};
\draw (0,-1.4) node{\footnotesize Exceptional point};
\draw (0,-1.1) node{\footnotesize $v_1 = \overline{v_2}$};
\end{tikzpicture}
\end{center}
\end{subfigure}\begin{subfigure}{0.33\linewidth}
\vspace{-4pt}
\begin{center}
\begin{tikzpicture}
\draw (0,0) node{\includegraphics[width=0.75\linewidth]{bandnpt.png}};
\draw[->] (-1.65,-0.56) -- (1.65,-0.56) node[below]{\footnotesize $\alpha$};
\draw[->] (0,-0.71) -- (0,1.2) node[left]{\footnotesize $\omega$};
\draw (0,-1.4) node{\footnotesize Partial inversion};
\draw (0,-1.1) node{\footnotesize $|v_1| > |v_2|$};
\draw (0.85,1.2) node{\includegraphics[width=0.33\linewidth]{L.png}};
\end{tikzpicture}
\end{center}
\end{subfigure}
\end{minipage}
\caption{Comparison between the topological phase transitions studied in \Cref{sec:bandinv} and \Cref{sec:pt_top}. In the Hermitian case studied in \Cref{sec:bandinv}, the Dirac cone found in the symmetric case $l_0 = 1/2$ can open into topologically distinct band gaps. In the non-Hermitian case, the degeneracy corresponds to an exceptional point. This exceptional point can open into separable bands which have distinct, albeit non-quantized, Zak phases.}\label{fig:phasetransition}
\end{figure}
In both the Hermitian and non-Hermitian cases, there are degeneracies associated to symmetric structures. Breaking the symmetry can open the degeneracy into topologically distinct band gaps as illustrated in \Cref{fig:phasetransition}. In the Hermitian case, the degeneracy is a linear intersection known as a \emph{Dirac cone}. In the non-Hermitian, $$-symmetric case, there can be exceptional point degeneracies which open when the $$-symmetry is detuned. In both the Hermitian and the non-Hermitian case, edge modes are created when the two topologically different phases are joined along an interface.

\subsection{Bound states in the continuum and Fano resonances}
Sometimes, localized modes can exist in periodic structures without a defect. This typically happens when the structure has certain symmetries, resulting in resonant modes in the radiation continuum whose far-field radiation vanishes. Such states are known as \emph{bound states in the continuum}.
\begin{defn}
A resonant mode $u_n^\alpha$ of \eqref{eq:scattering_quasi} is said to be a bound state in the continuum if the corresponding resonant frequency $\omega^\alpha_n$ is real, satisfies $|\alpha| < \omega^\alpha_n/v$ and the mode satisfies,
$$u_n^\alpha(x_l,x_0) = O(e^{-K|x_0|}), \quad |x_0| \to \infty, \quad K>0.$$
\end{defn}
As we shall see, there are two symmetry conditions required to achieve bound states in the continuum in the subwavelength resonator arrays. The first condition is a symmetry condition of the structure $D$, while the second condition is that $α= 0$, which corresponds to modes that radiate perpendicularly to the structure. We then have the following result from [8].
\begin{thm}
Assume that $d-d_l = 1$ and that $D = D_1\cup D_2 \subset Y$. We assume that $D$ satisfies the symmetry conditions $$\P D_1 = D_2, \quad \P_0 D_1 = D_1 \quad \P_0 D_2 = D_2,$$
and $\delta_1v_2^2 = \delta_2 v_2^2 \in \R$, where $\P,\P_0:\R^d\to\R^d$ are the parity operators
$$\P(x) = -x\quad\text{and}\quad \P_0(x_l,x_0) = (x_l,-x_0).$$
Moreover, we assume that $\alpha = 0$. Then, for small enough $\delta$, the second resonant mode $u_2^0$ is a bound state in the continuum.
\end{thm}
Interestingly, we can never have the first resonant mode $u_1^0$ as a bound state in the continuum. This is due to the fact that the first mode corresponds to the broad, low-frequency response of the screen which is not due to the local resonance of the structure.

An interesting problem is now to describe the behaviour when the symmetry conditions are no longer satisfied. In particular, when $_0 D_i ≠D_i$, we can use \Cref{thm:higher_periodic} to conclude that $ω_2^α$ has a small but non-zero imaginary part. As we will see next, this causes an interesting \emph{Fano-type} transmission scattering.

Given a unit vector $=(w_l,w_0)∈^d$ with $w_0>0$ we define the wave vectors
$$\k_+ = \frac{\omega}{v}\begin{pmatrix}w_l\\w_0\end{pmatrix}, \qquad \k_- = \frac{\omega}{v}\begin{pmatrix}w_l\\-w_0\end{pmatrix}.$$
We also let $α= ω/vw_l$. We now assume that
$$\uin(x) = c_1 e^{\i \k_-\cdot x} + c_2 e^{\i\k_+\cdot x},$$
and seek the behaviour of the solution $u$ of \eqref{eq:scattering_p}. In the first radiation continuum, the scattered field $u-$ consists of a single propagating mode as $|x_0| →∞$.
We will write $ f∼g$ to denote that two functions $f,g$ are equal up to exponentially decaying factors, in the sense that there is some constant $K>0$ such that
$$|f(x_l,x_0) - g(x_l,x_0)| = O(e^{-Kx_0}) \text{ as } x_0\rightarrow \infty.$$ We therefore have
\begin{equation}\label{eq:utot}
u \sim \begin{cases}c_1 e^{\i\k_-\cdot x} + d_1 e^{\i\k_+\cdot x}, & x_0\rightarrow \infty, \\ c_2 e^{\i\k_+\cdot x} + d_2 e^{\i\k_-\cdot x}, & x_0\rightarrow -\infty,\end{cases}
\end{equation}
\begin{equation} \label{eq:S_def}
\begin{pmatrix} d_1 \\ d_2 \end{pmatrix} =  S\begin{pmatrix} c_1 \\ c_2 \end{pmatrix}, \qquad S = \begin{pmatrix} r_+ & t_- \\ t_+ & r_- \end{pmatrix}.
\end{equation}
$S$ is known as the scattering matrix.
The reflection and transmission coefficients $r_+, t_+$ are the coefficients of the outgoing part of the field in the case $(x) =  e^ı_̨-·x$, \ie{} when the incident field is a plane wave from the positive $x_0$ direction (and reversely for $r_-, t_-$). The following theorem was proved in [8].
\begin{thm} \label{thm:phano}
Assume that  $\P D = D$ and that $0\leq \omega \leq K\sqrt{\delta}$ for some constant $K$. Let $c=\int_{\D}y_0\psi_1^0(t)\de \sigma(y)$ and assume that $c \neq 0$. Then we have the following asymptotic expansion of the scattering matrix as $\delta \to 0$
\begin{equation}
S = \frac{\omega_1^\alpha}{\omega_1^\alpha-\omega}\begin{pmatrix} 1&1\\1&1\end{pmatrix} +\frac{2\i\omega\Im(\omega_2^\alpha)}{(\omega_2^\alpha)^2-\omega^2}\begin{pmatrix} 1&-1\\-1&1\end{pmatrix} - \begin{pmatrix} 1&0\\0&1\end{pmatrix} + O(\delta^{1/2}),
\end{equation}
where the error term is uniform with respect to $\omega$.
\end{thm}
The scattering matrix $S$ contains two transmission peaks originating from the resonances $ω_1^α$ and $ω_2^α$. Since the imaginary part of $ω_2^α$ is very small, this corresponds to a sharp peak which will interfere with the broader peak associated to $ω_1^α$ to create an asymmetric, Fano-type, transmission peak as illustrated in \Cref{fig:fano}.

\begin{figure}[tbh]
\begin{subfigure}[b]{0.48\linewidth}
\centering
\begin{tikzpicture}[scale=0.4]
\begin{scope}[rotate=40]
\draw (0.8,0) circle (0.5);
\draw (-0.8,0) circle (0.5);
\end{scope}
\coordinate (oo) at (220:1.8);
\coordinate (a) at ($(0:3)+(oo)$);
\coordinate (b) at ($(40:3)+(oo)$);
\coordinate (o) at ($(0,0)+(oo)$);
\draw[opacity=0.7] (a) -- (o) -- (b);
\pic [draw, opacity=0.7, ->, "$\theta$", angle eccentricity=1.3,angle radius=0.7cm] {angle = a--o--b};
\node at (7,0){$\cdots$};
\node at (-7,0){$\cdots$};
\draw[->,thick] (-2,-4) -- (-1,-2) node[pos=0.7,below right]{$\uin$};
\draw[dashed] (-2,3) -- (-2,-2.5) (2,3) -- (2,-2.5);
\draw[<->] (-2,2.5) -- (2,2.5) node[pos=0.5,above]{$L$};
\begin{scope}[xshift=4cm,rotate=40]
\draw (0.8,0) circle (0.5);
\draw (-0.8,0) circle (0.5);
\end{scope}
\begin{scope}[xshift=-4cm,rotate=40]
\draw (0.8,0) circle (0.5);
\draw (-0.8,0) circle (0.5);
\end{scope}
\end{tikzpicture}
\vspace{0.5cm}
\caption{Sketch of a screen of resonators with an incident plane wave $\uin$. In this case, we have resonators arranged in a $\mathcal{P}$-symmetric dimer that is inclined at an angle of $\theta$ to the plane of the metascreen.} \label{fig:2D}
\end{subfigure}\hfill
\begin{subfigure}[b]{0.48\linewidth}
\begin{center}
\includegraphics[width=1\linewidth]{fano0T}
\end{center}
\caption{Transmission spectrum for the system in (a) with $\theta = 0.025\pi$. } \label{fig:fanoT}
\end{subfigure}
\caption{Sketch of a resonator screen (a) and corresponding transmission spectrum (b). If $\theta=0$, the resonant modes at $\alpha=0$ are bound states in the continuum and the transmission spectrum does not show a resonant peak at $\omega \approx 3$. For small but non-zero $\theta$, a characteristic Fano-type transmission peak appears at the resonant frequency.} \label{fig:fano}
\end{figure}

\subsection{Extraordinary transmission and unidirectional reflection}
We now investigate further scattering phenomena of resonator screens. In particular, we will extend the analysis in the previous section (valid in the Hermitian case when $v_i∈$) to the non-Hermitian case with a $$-symmetric screen. We assume that
$$\delta_1v_1^2 = a+\i b, \quad \delta_2v_2^2 = a-\i b,$$
in other words that $δ_1v_1^2 = (δ_2v_2^2 )$. Then we can apply \Cref{thm:higher_periodic} to conclude that
$$\omega_2^\alpha = \sqrt{\frac{2aC_{11}^0}{|D_1|}} + \frac{\i k_0}{4C_{11}^0}\left(\frac{b^2|Y_l|}{a^2} - \frac{c^2}{|Y_l|}\right) + O(\delta^{3/2}),$$
where $k_0$ and $c$ are defined as in the previous section. The following theorem, proved in [11], describes the scattering behaviour of this $$-symmetrical screen of resonators.
\begin{thm} \label{thm:unidir}
Assume that $D$ is $\P\T$-symmetric, so that $\P D_1 = D_2$ while $\delta_1v_1^2 = a+\i b$ and $\delta_2v_2^2 = a-\i b$ for $a>0,b\geq0$. Moreover, assume that $\P_3D_1=D_2$. Let $c=\int_{\D}y_0\psi_1^0(t)\de \sigma(y)$ and $\omega_* = \sqrt{\frac{2aC_{11}^0}{|D_1|}}$. Assume that $b|Y_l| \neq a|c| $ and that $\omega\in \R$ is in the subwavelength regime such that $\omega-\omega_* = O(\delta)$.	We then have the following asymptotic expansion of the scattering matrix:
$$S =\frac{2\i \omega\Im(\omega_2^\alpha )}{(\omega_2^\alpha)^2 - \omega^2}\begin{pmatrix} 1&-1\\-1&1\end{pmatrix} + \frac{2k_0 b}{a|D_1|\big((\omega_2^\alpha)^2 - \omega^2\big)}\begin{pmatrix} - ac & \i b|Y_l| \\ \i b|Y_l| & ac \end{pmatrix} - \begin{pmatrix} 1&0\\0&1\end{pmatrix} + O(\delta^{1/2}),$$
where the error term is uniform with respect to $\omega$ in a neighbourhood of $\omega_*$.	In particular, we have
$$r_\pm =- \frac{\omega_*^2-\omega^2 \pm \frac{2 k_0 bc}{|D_1|}}{(\omega_2^\alpha)^2 - \omega^2}+ O(\delta^{1/2}),$$
and, at leading order, $r_+$ and $r_-$ vanish at $\omega = \omega_+$ and $\omega = \omega_-$, respectively, which are given by
$$\omega_+^2 = \omega_*^2 +  \frac{2 k_0 bc}{|D_1|}, \quad \omega_-^2 = \omega_*^2-\frac{2 k_0 bc}{|D_1|}.$$
\end{thm}
Comparing \Cref{thm:phano} (which is valid only for $b=0$) and \Cref{thm:unidir}, we see that there is an extra term in \Cref{thm:unidir} corresponding to the non-zero gain and loss $b≠0$. This term is responsible for the approximate zeros of $r_+$ and $r_-$ leading to unidirectional reflection. This is numerically verified in \Cref{fig:unidir}, where we observe that the two reflectances vanish on different sides of the critical frequency $ω_*$. Moreover, at least formally, we see that when $(ω_2^α)=0$, the singularity of the scattering matrix $S$ will not vanish in the case $b≠0$. In the Hermitian case, real resonances correspond to bound states in the continuum, which decouple from the far-field and therefore cannot be excited by incoming waves. In the non-Hermitian case, however, we can have real resonances with modes which are excited by incoming waves. Such resonances correspond to extraordinary transmission, where the transmitted field is greatly amplified. This amplification, which is impossible in the Hermitian case due to energy conservation, is possible due to the energy input in the non-Hermitian case [11].

\begin{figure}[h]
\begin{subfigure}{0.35\linewidth}
\begin{center}
\begin{tikzpicture}[scale=0.15]
\draw (0,-1.5) circle (1);
\draw (0,1.5) circle (1);
\node at (0,-1.5){\small $-$};
\node at (0,1.5){\small $+$};
\node at (11,0){$\cdots$};
\node at (-11,0){$\cdots$};
\draw[->,thick] (-4,6) -- (-1,3) node[pos=0.7,xshift=4pt,yshift=8pt]{\small $\uin$};
\draw[->,thick] (1,3) -- (4,6) node[pos=0.7,right]{\small $R_+$};
\draw[->,thick] (2.5,-3) -- (5.5,-6) node[pos=0.7,right]{\small $T_+$};
\begin{scope}[xshift=4cm]
\draw (0,-1.5) circle (1);
\draw (0,1.5) circle (1);
\node at (0,-1.5){\small $-$};
\node at (0,1.5){\small $+$};
\end{scope}
\begin{scope}[xshift=-4cm]
\draw (0,-1.5) circle (1);
\draw (0,1.5) circle (1);
\node at (0,-1.5){\small $-$};
\node at (0,1.5){\small $+$};
\end{scope}
\begin{scope}[xshift=8cm]
\draw (0,-1.5) circle (1);
\draw (0,1.5) circle (1);
\node at (0,-1.5){\small $-$};
\node at (0,1.5){\small $+$};
\end{scope}
\begin{scope}[xshift=-8cm]
\draw (0,-1.5) circle (1);
\draw (0,1.5) circle (1);
\node at (0,-1.5){\small $-$};
\node at (0,1.5){\small $+$};
\end{scope}
\begin{scope}[yshift=-16cm]
\draw (0,-1.5) circle (1);
\draw (0,1.5) circle (1);
\node at (0,-1.5){\small $-$};
\node at (0,1.5){\small $+$};
\node at (11,0){$\cdots$};
\node at (-11,0){$\cdots$};
\draw[->,thick] (-4,-6) -- (-1,-3) node[pos=0.7,left]{\small $\uin$};
\draw[->,thick] (2.5,3) -- (5.5,6) node[pos=0.7,right]{\small $T_-$};
\draw[->,thick] (1,-3) -- (4,-6) node[pos=0.7,right]{\small $R_-$};
\begin{scope}[xshift=4cm]
\draw (0,-1.5) circle (1);
\draw (0,1.5) circle (1);
\node at (0,-1.5){\small $-$};
\node at (0,1.5){\small $+$};
\end{scope}
\begin{scope}[xshift=-4cm]
\draw (0,-1.5) circle (1);
\draw (0,1.5) circle (1);
\node at (0,-1.5){\small $-$};
\node at (0,1.5){\small $+$};
\end{scope}
\begin{scope}[xshift=8cm]
\draw (0,-1.5) circle (1);
\draw (0,1.5) circle (1);
\node at (0,-1.5){\small $-$};
\node at (0,1.5){\small $+$};
\end{scope}
\begin{scope}[xshift=-8cm]
\draw (0,-1.5) circle (1);
\draw (0,1.5) circle (1);
\node at (0,-1.5){\small $-$};
\node at (0,1.5){\small $+$};
\end{scope}
\end{scope}
\end{tikzpicture}
\end{center}
\end{subfigure}
\begin{subfigure}{0.6\linewidth}
\begin{center}
\includegraphics[width=\linewidth]{scattering2}
\end{center}
\end{subfigure}
\caption{Plot of the transmittance $T_\pm = |t_\pm|^2$ (blue) and reflectance $R_\pm = |r_\pm|^2$ (red) as functions of the frequency. Due to reciprocity we have $T_+ = T_-$. However, we observe that $R_+\neq R_-$. In particular, there are points when $R_\pm$ vanish while $R_\mp$ is non-zero, known as unidirectional reflection. Moreover, in this non-Hermitian case the transmission coefficients are not bounded by unity and can attain large peak values, known as extraordinary transmission.} \label{fig:unidir}
\end{figure}

\subsection{Time-modulated metamaterials}\label{sec:time}
We have now seen a variety of phenomena that can occur when the material parameters vary periodically in the spatial variable $x$. Mathematically, we can treat the time variable $t$ in a similar fashion. If the material parameters depend periodically on $t$, there can be conceptually similar phenomena which nevertheless have fundamentally different physical implications.

The Helmholtz equation we have studied so far originates from the scalar wave equation when posed in the frequency domain, and is valid only when the material parameters are constant in time. To study the time-dependent case, we return to the wave equation
\begin{equation}\label{eq:wave}
\left(\frac{\p }{\p t } \frac{1}{\kappa(x,t)} \frac{\p}{\p t} - \nabla \cdot \frac{1}{\rho(x,t)} \nabla\right) u(x,t) = 0, \quad x\in \R^d, t\in \R.
\end{equation}
Here, $κ(x,t)$ and $ρ(x,t)$ are the material parameters. We consider the case of a  finite collection of resonators $D=D_1∪...∪D_N$ in $d=3$ spatial dimensions. We assume that the modulation is only performed inside the resonators, so that
\begin{equation} \label{eq:resonatormod}
\kappa(x,t) = \begin{cases}
\kappa_0, & x \in \R^3 \setminus \overline{D}, \\  \kappa_i(t), & x\in D_i,
\end{cases}, \qquad \rho(x,t) = \begin{cases}
\rho_0, & x \in \R^3 \setminus \overline{D}, \\  \rho_i(t), & x\in D_i. \end{cases}
\end{equation}
We assume that $κ_i ∈C^1()$ for each $i=1,...,N$. We will consider the case when the modulation of $1/ρ$ and $1/κ$ consist of a finite Fourier series with a large number of Fourier coefficients:
$$\frac{1}{\rho_i(t)} = \sum_{n = -M}^M r_{i,n} e^{\i n \Omega t}, \qquad \frac{1}{\kappa_i(t)} = \sum_{n = -M}^M k_{i,n} e^{\i n \Omega t},$$
for some $M∈$ satisfying $M = O(δ^-γ/2)$ for some $0<γ<1$. Moreover, we assume that they are periodic in $t$ with period $T$ and frequency $Ω= 2π/T$.

As before we define the (time-dependent) contrast parameter and wave speeds as
$$\delta_i(t) = \frac{\rho_i(t)}{\rho_0}, \quad v_i(t) = \sqrt{\frac{\kappa_i(t)}{\rho_i(t)}}, \quad v_0 = \sqrt{\frac{\kappa_0}{\rho_0}}, $$
and assume that
$$\delta_i(t) = O(\delta), \quad v_i(t) = O(1) \quad v = O(1), \quad \text{for all } \ t\in \R,$$
for $i=1,...,N$, where $δ≪1$.

The notion of frequency is slightly altered in this time-modulated setting. Since the wave equation \eqref{eq:wave} is periodic in $t$, we can apply the Floquet transform in $t$ and obtain the differential problem
\begin{equation} \label{eq:wave_transf}
\begin{cases}\ \ds \left(\frac{\p }{\p t } \frac{1}{\kappa(x,t)} \frac{\p}{\p t} - \nabla \cdot \frac{1}{\rho(x,t)} \nabla\right) u(x,t) = 0,\\[0.3em]
\	u(x,t)e^{-\i \omega t} \text{ is $T$-periodic in $t$}.
\end{cases}
\end{equation}
It is apparent that $ω$ plays the role of the quasiperiodicity in the Floquet-Bloch theory, and we will refer to $ω$ as a quasifrequency. Since $ω$ is defined modulo $Ω$, we define the time-Brillouin zone $Y_t^* := ℂ / (Ω)$. Observe that we allow complex quasifrequencies; real $ω$ correspond to solutions which are bounded in $t$, whereas non-real $ω$ correspond to exponentially increasing or decaying solutions.

Due to the periodic nature of $Y_t^*$, the usual definition of subwavelength frequencies does not apply to quasifrequencies. For example, in the particular case when $Ω= O(δ^1/2)$ (which will be of interest later on), the whole Brillouin zone scales as $O(δ^1/2)$, meaning that \emph{all} quasifrequencies tend to zero as $δ→0$. In order to distinguish between these quasifrequencies, we introduce the following definition.
\begin{defn}[Subwavelength quasifrequency] \label{def:sub}
A quasifrequency $\omega = \omega(\delta) \in Y^*_t$ of \eqref{eq:wave_transf} is said to be a subwavelength quasifrequency if there is a corresponding Bloch solution $u(x,t)$, depending continuously on $\delta$, which can be written as
$$u(x,t)= e^{\i \omega t}\sum_{n = -\infty}^\infty v_n(x)e^{\i n\Omega t},$$
$$\omega \rightarrow 0 \ \text{and} \ M\Omega \rightarrow 0 \ \text{as} \ \delta \to 0,$$
for some integer-valued function $M=M(\delta)$ such that, as $\delta \to 0$, we have
$$\sum_{n = -\infty}^\infty \|v_n\|_{L^2(D)} = \sum_{n = -M}^M \|v_n\|_{L^2(D)} + o(1).$$
\end{defn}
The following theorem, proved in [19], gives the capacitance matrix approximation to the subwavelength quasifrequencies as $δ→0$.
\begin{thm}
Assume that the material parameters are given by \eqref{eq:resonatormod}. Then, as $\delta \to 0$, the quasifrequencies $\omega \in Y^*_t$ to the wave equation \eqref{eq:wave} in the subwavelength regime are, to leading order, given by the quasifrequencies of the system of ordinary differential equations for $y_i(t)$,
\begin{equation}\label{eq:C_ODE}
\sum_{j=1}^N C_{ij} y_j(t) = -\frac{|D_i|}{\delta_i(t)}\frac{\dx}{\dx t}\left(\frac{1}{\delta_iv_i^2}\frac{\dx (\delta_iy_i)}{\dx t}\right),
\end{equation}
for $i=1,...,N$.
\end{thm}
The matrix appearing in the left-hand side of \eqref{eq:C_ODE} is the capacitance matrix. We can rewrite \eqref{eq:C_ODE} into the following system of Hill equations
\begin{equation}	\label{eq:hill}
\Psi''(t)+ M(t)\Psi(t)=0,
\end{equation}
where the vector $Ψ$ and the matrix $M$ are defined as
$$\Psi(t) = \left(\frac{\sqrt{\delta_i(t)}}{v_i(t)}y_i(t)\right)_{i=1}^N, \quad M(t) = W_1(t)C W_2(t) + W_3(t),$$
with $W_1, W_2$ and $W_3$ being the diagonal matrices with diagonal entries
$$\left(W_1\right)_{ii} = \frac{v_i\delta_i^{3/2}}{|D_i|}, \qquad \left(W_2\right)_{ii} =\frac{v_i}{\sqrt{\delta_i}}, \qquad \left(W_3\right)_{ii} = \frac{\sqrt{\delta_i}v_i}{2}\frac{\dx }{\dx t}\frac{1}{(\delta_iv_i^2)^{3/2}}\frac{\dx (\delta_iv_i^2)}{\dx t},$$
for $i=1,...,N$.
\begin{figure}
\centering
\begin{subfigure}[b]{0.4\linewidth}
\centering
\begin{tikzpicture}
\draw (-1.5,0) circle (1);
\draw (1.5,0) circle (1);
\node at (-1.5,0){$D_1$};
\node at (1.5,0){$D_2$};
\node at (-1.6,-1.5){\small $\kappa_1(t) =1+\epsilon\sin(\Omega t)$};
\node at (1.6,-1.5){\small $\kappa_2(t) = 1-\epsilon\sin(\Omega t)$};
\end{tikzpicture}
\vspace{0.5cm}
\caption{A $\mathcal{PT}$-symmetric pair of spherical resonators.}
\end{subfigure}
\hspace{0.1cm}
\begin{subfigure}[b]{0.55\linewidth}
\includegraphics[width=\linewidth]{EPtime}
\caption{The two subwavelength resonant frequencies.}
\end{subfigure}
\caption{Similarly as in the non-Hermitian case demonstrated in \Cref{fig:EP2}, there can be exceptional points in the time-modulated case. The two subwavelength quasifrequencies of a pair of $\mathcal{PT}$-symmetric, time-modulated resonators can be approximated through the capacitance formulation, and an asymptotic exceptional point occurs at $\epsilon \approx0.3$.} \label{fig:EPtime}
\end{figure}
The time-modulated case shares many similarities with the non-Hermitian case (with complex parameters) studied before. Both these cases have energy input and output to the system, and not surprisingly we can find exceptional points in the time-modulated systems. \Cref{fig:EPtime} shows the emergence of an exceptional point in a dimer structure, similar to \Cref{fig:EP2} but instead due to the time-modulation.
\subsection{Near-zero metamaterials}
In \Cref{sec:doubleneg} we saw the emergence of exotic parameter values, namely negative effective parameters, due to the small-scale structure of the metamaterial. In the present section we will observe effective material parameters which are close to zero, in which case the Helmholtz equation reduces to the Laplace equation and wave propagation occurs without phase change (corresponding to ``infinite'' phase velocity).

\subsubsection{Near-zero refractive index in honeycomb crystals}
\begin{figure}[htb]
\begin{subfigure}[b]{0.5\linewidth}
\centering
\begin{tikzpicture}[scale=1]
\begin{scope}[xshift=-5cm,scale=1]
\coordinate (a) at (1,{1/sqrt(3)});
\coordinate (b) at (1,{-1/sqrt(3)});
\pgfmathsetmacro{\rb}{0.25pt}
\pgfmathsetmacro{\rs}{0.2pt}

\draw (0,0) -- (1,{1/sqrt(3)}) node[xshift=2pt,yshift=3pt]{\footnotesize$Y$} -- (2,0) -- (1,{-1/sqrt(3)}) -- cycle;
\begin{scope}[xshift = 1.33333cm]
\draw plot [smooth cycle] coordinates {(0:\rb) (60:\rs) (120:\rb) (180:\rs) (240:\rb) (300:\rs) };
\draw (0,0) node{\footnotesize $D_2$};
\end{scope}
\begin{scope}[xshift = 0.666667cm, rotate=60]
\draw plot [smooth cycle] coordinates {(0:\rb) (60:\rs) (120:\rb) (180:\rs) (240:\rb) (300:\rs) };
\draw (0,0) node{\footnotesize $D_1$};
\end{scope}

\draw[opacity=0.2] ({2/3},0) -- ({4/3},0)
($0.5*(1,{1/sqrt(3)})$) -- ({2/3},0)
($0.5*(1,{-1/sqrt(3)})$) -- ({2/3},0)
($(1,{1/sqrt(3)})+0.5*(1,{-1/sqrt(3)})$) -- ({4/3},0)
($0.5*(1,{1/sqrt(3)})+(1,{-1/sqrt(3)})$) -- ({4/3},0);

\begin{scope}[shift = (a)]
\begin{scope}[xshift = 1.33333cm]
\draw plot [smooth cycle] coordinates {(0:\rb) (60:\rs) (120:\rb) (180:\rs) (240:\rb) (300:\rs) };
\end{scope}
\begin{scope}[xshift = 0.666667cm, rotate=60]
\draw plot [smooth cycle] coordinates {(0:\rb) (60:\rs) (120:\rb) (180:\rs) (240:\rb) (300:\rs) };
\end{scope}
\draw[opacity=0.2] ({2/3},0) -- ({4/3},0)
($0.5*(1,{1/sqrt(3)})$) -- ({2/3},0)
($0.5*(1,{-1/sqrt(3)})$) -- ({2/3},0)
($(1,{1/sqrt(3)})+0.5*(1,{-1/sqrt(3)})$) -- ({4/3},0)
($0.5*(1,{1/sqrt(3)})+(1,{-1/sqrt(3)})$) -- ({4/3},0);
\end{scope}
\begin{scope}[shift = (b)]
\begin{scope}[xshift = 1.33333cm]
\draw plot [smooth cycle] coordinates {(0:\rb) (60:\rs) (120:\rb) (180:\rs) (240:\rb) (300:\rs) };
\end{scope}
\begin{scope}[xshift = 0.666667cm, rotate=60]
\draw plot [smooth cycle] coordinates {(0:\rb) (60:\rs) (120:\rb) (180:\rs) (240:\rb) (300:\rs) };
\end{scope}
\draw[opacity=0.2] ({2/3},0) -- ({4/3},0)
($0.5*(1,{1/sqrt(3)})$) -- ({2/3},0)
($0.5*(1,{-1/sqrt(3)})$) -- ({2/3},0)
($(1,{1/sqrt(3)})+0.5*(1,{-1/sqrt(3)})$) -- ({4/3},0)
($0.5*(1,{1/sqrt(3)})+(1,{-1/sqrt(3)})$) -- ({4/3},0);
\end{scope}
\begin{scope}[shift = ($-1*(a)$)]
\begin{scope}[xshift = 1.33333cm]
\draw plot [smooth cycle] coordinates {(0:\rb) (60:\rs) (120:\rb) (180:\rs) (240:\rb) (300:\rs) };
\end{scope}
\begin{scope}[xshift = 0.666667cm, rotate=60]
\draw plot [smooth cycle] coordinates {(0:\rb) (60:\rs) (120:\rb) (180:\rs) (240:\rb) (300:\rs) };
\end{scope}
\draw[opacity=0.2] ({2/3},0) -- ({4/3},0)
($0.5*(1,{1/sqrt(3)})$) -- ({2/3},0)
($0.5*(1,{-1/sqrt(3)})$) -- ({2/3},0)
($(1,{1/sqrt(3)})+0.5*(1,{-1/sqrt(3)})$) -- ({4/3},0)
($0.5*(1,{1/sqrt(3)})+(1,{-1/sqrt(3)})$) -- ({4/3},0);
\end{scope}
\begin{scope}[shift = ($-1*(b)$)]
\begin{scope}[xshift = 1.33333cm]
\draw plot [smooth cycle] coordinates {(0:\rb) (60:\rs) (120:\rb) (180:\rs) (240:\rb) (300:\rs) };
\end{scope}
\begin{scope}[xshift = 0.666667cm, rotate=60]
\draw plot [smooth cycle] coordinates {(0:\rb) (60:\rs) (120:\rb) (180:\rs) (240:\rb) (300:\rs) };
\end{scope}
\draw[opacity=0.2] ({2/3},0) -- ({4/3},0)
($0.5*(1,{1/sqrt(3)})$) -- ({2/3},0)
($0.5*(1,{-1/sqrt(3)})$) -- ({2/3},0)
($(1,{1/sqrt(3)})+0.5*(1,{-1/sqrt(3)})$) -- ({4/3},0)
($0.5*(1,{1/sqrt(3)})+(1,{-1/sqrt(3)})$) -- ({4/3},0);
\end{scope}
\begin{scope}[shift = ($(a)+(b)$)]
\begin{scope}[xshift = 1.33333cm]
\draw plot [smooth cycle] coordinates {(0:\rb) (60:\rs) (120:\rb) (180:\rs) (240:\rb) (300:\rs) };
\end{scope}
\begin{scope}[xshift = 0.666667cm, rotate=60]
\draw plot [smooth cycle] coordinates {(0:\rb) (60:\rs) (120:\rb) (180:\rs) (240:\rb) (300:\rs) };
\end{scope}
\draw[opacity=0.2] ({2/3},0) -- ({4/3},0)
($0.5*(1,{1/sqrt(3)})$) -- ({2/3},0)
($0.5*(1,{-1/sqrt(3)})$) -- ({2/3},0)
($(1,{1/sqrt(3)})+0.5*(1,{-1/sqrt(3)})$) -- ({4/3},0)
($0.5*(1,{1/sqrt(3)})+(1,{-1/sqrt(3)})$) -- ({4/3},0);
\end{scope}
\begin{scope}[shift = ($-1*(a)-(b)$)]
\begin{scope}[xshift = 1.33333cm]
\draw plot [smooth cycle] coordinates {(0:\rb) (60:\rs) (120:\rb) (180:\rs) (240:\rb) (300:\rs) };
\end{scope}
\begin{scope}[xshift = 0.666667cm, rotate=60]
\draw plot [smooth cycle] coordinates {(0:\rb) (60:\rs) (120:\rb) (180:\rs) (240:\rb) (300:\rs) };
\end{scope}
\draw[opacity=0.2] ({2/3},0) -- ({4/3},0)
($0.5*(1,{1/sqrt(3)})$) -- ({2/3},0)
($0.5*(1,{-1/sqrt(3)})$) -- ({2/3},0)
($(1,{1/sqrt(3)})+0.5*(1,{-1/sqrt(3)})$) -- ({4/3},0)
($0.5*(1,{1/sqrt(3)})+(1,{-1/sqrt(3)})$) -- ({4/3},0);
\end{scope}
\begin{scope}[shift = ($(a)-(b)$)]
\begin{scope}[xshift = 1.33333cm]
\draw plot [smooth cycle] coordinates {(0:\rb) (60:\rs) (120:\rb) (180:\rs) (240:\rb) (300:\rs) };
\end{scope}
\begin{scope}[xshift = 0.666667cm, rotate=60]
\draw plot [smooth cycle] coordinates {(0:\rb) (60:\rs) (120:\rb) (180:\rs) (240:\rb) (300:\rs) };
\end{scope}
\draw[opacity=0.2] ({2/3},0) -- ({4/3},0)
($0.5*(1,{1/sqrt(3)})$) -- ({2/3},0)
($0.5*(1,{-1/sqrt(3)})$) -- ({2/3},0)
($(1,{1/sqrt(3)})+0.5*(1,{-1/sqrt(3)})$) -- ({4/3},0)
($0.5*(1,{1/sqrt(3)})+(1,{-1/sqrt(3)})$) -- ({4/3},0);
\end{scope}
\begin{scope}[shift = ($-1*(a)+(b)$)]
\begin{scope}[xshift = 1.33333cm]
\draw plot [smooth cycle] coordinates {(0:\rb) (60:\rs) (120:\rb) (180:\rs) (240:\rb) (300:\rs) };
\end{scope}
\begin{scope}[xshift = 0.666667cm, rotate=60]
\draw plot [smooth cycle] coordinates {(0:\rb) (60:\rs) (120:\rb) (180:\rs) (240:\rb) (300:\rs) };
\end{scope}
\draw[opacity=0.2] ({2/3},0) -- ({4/3},0)
($0.5*(1,{1/sqrt(3)})$) -- ({2/3},0)
($0.5*(1,{-1/sqrt(3)})$) -- ({2/3},0)
($(1,{1/sqrt(3)})+0.5*(1,{-1/sqrt(3)})$) -- ({4/3},0)
($0.5*(1,{1/sqrt(3)})+(1,{-1/sqrt(3)})$) -- ({4/3},0);
\end{scope}
\end{scope}
\end{tikzpicture}
\end{subfigure}
\begin{subfigure}[b]{0.4\linewidth}
\centering
\begin{tikzpicture}[scale=1.5]
\coordinate (a) at ({1/sqrt(3)},1);
\coordinate (b) at ({1/sqrt(3)},-1);
\coordinate (c) at ({2/sqrt(3)},0);
\coordinate (K1) at ({1/sqrt(3)},{1/3});
\coordinate (K2) at ({1/sqrt(3)},{-1/3});
\coordinate (K3) at (0,{-2/3});
\coordinate (K4) at ({-1/sqrt(3)},{-1/3});
\coordinate (K5) at ({-1/sqrt(3)},{1/3});
\coordinate (K6) at (0,{2/3});

\draw[fill] (K1) circle(1pt) node[xshift=6pt,yshift=4pt]{$\alpha_1^*$};
\draw[fill] (K2) circle(1pt) node[xshift=6pt,yshift=-4pt]{$\alpha_2^*$};
\draw[fill] (0,0) circle(1pt) node[left]{$\Gamma$};

\draw (K1) -- (K2) -- (K3) -- (K4) -- (K5) node[above]{$Y^*$} -- (K6) -- cycle;
\end{tikzpicture}
\vspace{10pt}
\end{subfigure}
\caption{Honeycomb crystal and corresponding Brillouin zone.}\label{fig:honeycomb}
\end{figure}

We consider a two-dimensional infinite honeycomb crystal in two dimensions depicted in Figure \ref{fig:honeycomb}. We let the lattice $Λ$ be generated by the lattice vectors
$$ l_1 = L\left( \frac{\sqrt{3}}{2}, \frac{1}{2} \right),~~l_2 = L\left( \frac{\sqrt{3}}{2}, -\frac{1}{2}\right).$$
We assume that each unit cell contains two resonators, $D=D_1∪D_2$, such that each resonator is invariant under rotation by $2π/3$ and so that $D$ is invariant under rotation by $π$.

The dual lattice  $Λ^*$ is generated by $α_1$ and $α_2$ given by
$$ \alpha_1 = \frac{2\pi}{L}\left( \frac{1}{\sqrt{3}}, 1\right),~~\alpha_2 = \frac{2\pi}{L}\left(\frac{1}{\sqrt{3}}, -1 \right).$$
The points $$\alpha_1^*= \frac{2\alpha_1+\alpha_2}{3}, \quad \alpha^*_2 = \frac{\alpha_1+2\alpha_2}{3},$$ in the Brillouin zone are called \emph{Dirac points}. Next, we will study the band functions and Bloch modes around these points. For simplicity, we only consider the analysis around the Dirac point $α_* := α_1^*$, the other point having a similar behaviour.

\begin{figure}[tbh]
\centering
\begin{subfigure}[b]{0.48\linewidth}
\centering
\includegraphics[height=5.5cm]{hl_ndilute_zoom}
\caption{Subwavelength band structure of a honeycomb structure, exhibiting a Dirac cone at $\alpha = K$.}
\label{fig:cone}
\end{subfigure}\hfill
\begin{subfigure}[b]{0.48\linewidth}
\centering
\includegraphics[height=6cm]{Sr1}
\caption{Small-scale behaviour of the eigenmodes at the Dirac point.}
\label{fig:Smr}
\end{subfigure}
\caption{The honeycomb structure exhibits a Dirac cone at the corner of the Brillouin zone. Corresponding eigenmodes are rapidly oscillating, and is periodic across one hexagon in the honeycomb structure.} \label{fig:smallscale}
\end{figure}

\begin{figure}[tbh]
\begin{subfigure}[b]{0.47\linewidth}
\hspace{-10pt}
\includegraphics[width=\linewidth]{h1D1}
\caption{Large-scale behaviour of the eigenmodes close to the Dirac point. }
\label{fig:h1D}
\end{subfigure}\hfill
\begin{subfigure}[b]{0.49\linewidth}
\includegraphics[width=\linewidth]{linear}
\caption{Spatial frequency $f$ of the envelopes as a function of the frequency shift $\varepsilon=\omega-\omega_*$.}
\end{subfigure}
\caption{For frequencies close to the Dirac frequency, the eigenfunctions oscillate on two distinct scales. The spatial frequency of the large-scale oscillations vanish at the centre of the Dirac cone, corresponding to zero refractive index.} \label{fig:largescale}
\end{figure}

At $α= α_*$, the generalized capacitance matrix $^α$ has an eigenvalue of multiplicity 2: $λ_1^α_* = λ_2^α_*$. The next result shows that this asymptotic degeneracy is in fact an exact degeneracy, and moreover that the band functions intersect in a Dirac cone at this point (see \Cref{fig:cone}) [14].
\begin{thm}\label{thm:honeycomb}
For $\alpha$ close to $\alpha_*$ and $\delta$ small enough, the first two band functions form a Dirac cone, \ie{},
\begin{equation} \label{eq:dirac}
\begin{matrix}
\ds \omega_1^\alpha = \omega_*- \mu|\alpha - \alpha_*| \big[ 1+ O(|\alpha-\alpha_*|) \big], \\[0.5em]
\ds \omega_2^\alpha = \omega_*+ \mu|\alpha - \alpha_*| \big[ 1+ O(|\alpha-\alpha_*|) \big],
\end{matrix}
\end{equation}
where $\omega_*$ and $\mu$ are independent of $\alpha$ and satisfy
$$\omega_*= \sqrt{\lambda_1^{\alpha_*}} + O(\delta) \quad \text{and} \quad \mu =  |c|\sqrt\delta\mu_0 + O(\delta),  \quad \mu_0=\frac{1}{2}\sqrt{\frac{v_\mathrm{r}^2 }{|D_1|C_{11}^{\alpha_*}}}, \quad c = \left|\frac{\partial C_{12}^{\alpha}}{\partial \alpha_1}\Big|_{\alpha=\alpha_*}\right|,$$
as $\delta \rightarrow 0$. Moreover, the error term $O(|\alpha-\alpha_*|)$ in \eqref{eq:dirac} is uniform in $\delta$.
\end{thm}
Next, we will see that the near-zero property follows as a direct consequence of the Dirac cone. We consider a homogenization setting close to the Dirac points. We rescale the unit cell by replacing $Y$ with $sY$ for some small $s>0$. To fix the order of the resonant frequencies, \ie{} $ω_i^α= O(1)$, we assume that $δ= O(s^2)$ as $s→0$. We then have the following result [20].
\begin{thm} \label{thm:nearzero}
For frequencies $\omega$ close to the Dirac frequency $\omega_*$, namely, $\omega-\omega_* = \beta \sqrt\delta$, the following asymptotic behaviour of the Bloch eigenfunction $u^{\alpha_*/s + \tilde{\alpha}}_s$ holds:
u_{s}^{\alpha_*/s+\tilde{\alpha}}(x) = \begin{bmatrix}A e^{\i\tilde{\alpha}\cdot x}\\ B e^{\i\tilde{\alpha}\cdot x}\end{bmatrix}\cdot \textbf{S}_D^{\alpha_*,k}\left(\frac{x}{s}\right) + O(s),
where the macroscopic field $[\tilde{u}_{1}, \tilde{u}_{2}]^T:=[A e^{\i\tilde{\alpha}\cdot x}, B e^{\i\tilde{\alpha}\cdot x}]^T$ satisfies the two-dimensional Dirac equation
\mu_0
\begin{bmatrix}
0 & (-c\i)( \p_1 - \i \p_2)
\\
(-\overline{c}\i)( \p_1 + \i \p_2) & 0
\end{bmatrix}
\begin{bmatrix}
\tilde{u}_{1}
\\
\tilde{u}_{2}
\end{bmatrix}
= \frac{\omega-\omega_*}{\sqrt\delta}
\begin{bmatrix}
\tilde{u}_{1}
\\
\tilde{u}_{2}
\end{bmatrix}.
\end{thm}
The system of Dirac equations can be considered as a homogenized equation for the honeycomb structure. Each $ũ_j$ satisfy the Helmholtz equation
\begin{equation} \label{eq:hom}
\Delta \tilde{u}_j + \frac{(\omega-\omega_*)^2}{\mu^2} \tilde{u}_j = 0.
\end{equation}
In particular, at $ω=ω_*$ this effective equation reduces to the Laplace equation corresponding to effective zero refractive index. Equation \eqref{eq:hom} describes the large-scale behaviour of the eigenmodes, illustrated in \Cref{fig:largescale}. We emphasize that in addition to this large-scale behaviour, there will be small-scale oscillations described by the functions $S_D^α_*,k$ as illustrated in \Cref{fig:smallscale}.

\subsubsection{Double-near zero in time-modulated materials}

\begin{figure}[h]
\begin{center}
\includegraphics[width=0.75\linewidth]{honeycomb}
\end{center}
\caption{Illustration of the trimer honeycomb lattice, with phase-shifted time-modulations inside the trimers.} \label{fig:trihoney}
\end{figure}

\begin{figure}[h]
\begin{center}
\begin{subfigure}[t]{0.4\linewidth}
\vspace{0pt}
\begin{center}
\begin{tikzpicture}
\draw (0.25,1.38) node{\includegraphics[width=1\linewidth]{HLresmod_a}};
%\draw[draw=none,fill,opacity=0.3] (0,0) -- (1.2,3.2) -- (-1.2,3.2) -- cycle;
%\draw (0,3) node{\tiny first radiation};
%\draw (0,2.8) node{\tiny continuum};
\end{tikzpicture}
\end{center}
\caption{Static band structure}
\end{subfigure}
\hspace{20pt}
\begin{subfigure}[t]{0.4\linewidth}
\vspace{0pt}
\begin{center}
\begin{tikzpicture}
\draw (0.25,1.38) node{\includegraphics[width=1\linewidth]{HLresmod_b}};
%\draw[draw=none,fill,opacity=0.3] (0,0) -- (1.2,3.2) -- (-1.2,3.2) -- cycle;
%\draw (0,3) node{\tiny first radiation};
%\draw (0,2.8) node{\tiny continuum};
\end{tikzpicture}
\end{center}
\caption{Modulated band structure}
\end{subfigure}
\end{center}
\caption{In time-modulated structures, Dirac cones may appear at the origin of the Brillouin zone. Such points are associated to double near-zero materials, enabling wave transmission without phase changes and with strong interactions with incoming waves.}\label{fig:diracGamma}
\end{figure}

The Dirac cone observed in the previous section is located at the corner of the Brillouin zone, which means that corresponding modes do not lie in the radiation continuum. For physical structures of finite dimensions, these waves will be confined to the structure and will not interact with incoming waves. Physically, this can be seen as an impedance mismatch, leading to poor transmission, due to the fact that this structure corresponds to a \emph{single}-near zero material, where only one of the material parameters vanish.

In order to overcome this low transmission, it is desirable to create structures with Dirac cones at the origin of the Brillouin zone. Corresponding modes lie in the radiation continuum, and will therefore interact with incoming waves. In order to achieve this, we study a time-modulated honeycomb structure as illustrated in \Cref{fig:trihoney}. Using the theory from \Cref{sec:time} we can compute the (quasi-)band structure associated to this material. In \Cref{fig:diracGamma} we see the static (folded) band structure, with six bands in the subwavelength regime. For a particular modulation strength, we see that a Dirac cone degeneracy may appear at the origin of the Brillouin zone, enabling a double-near zero material around this point.

\section{Concluding remarks}

In this review, we have studied several Helmholtz scattering problems posed in the subwavelength regime and repeatedly encountered the concept of capacitance. We have studied the mathematical properties of the generalized capacitance matrix, in both finite and infinite, periodic settings. We were then able to use the capacitance formulation to study several different interesting physical phenomena in the field of subwavelength metamaterials. This formulation emerged from a functional analytic approach, where the generalized capacitance matrix describes the perturbation of the kernel of a non-linear integral operator for asymptotically small parameter values. Similar approaches apply to a plethora of other subwavelength scattering problems, including high-contrast dielectric particles, plasmonic nanoparticles and Helmholtz resonators [3, 24, 25, 26, 27]. Suitable capacitance formulations can thereby be used to characterize a wide range of subwavelength resonance phenomena.

Approximating classical wave systems in terms of generalized capacitance matrices shares similarities with the tight-binding approximation that is used widely in quantum theory, in the sense that both formulations provide a discrete approximation to a continuous differential problem. As we have seen, there are nevertheless fundamental differences, for instance due to the strong interactions between the subwavelength resonators. As observed in \Cref{rmk:tightbind}, the correspondence between the capacitance formulation and the tight-binding approximation holds only in the case of dilute resonators. Even in the dilute case, long-range interactions between subwavelength resonators are relatively strong and nearest-neighbour approximations are not generally appropriate. The strength of the capacitance formulation is that the capacitance matrix accounts for these strong interactions, thereby providing a unified mathematical model to study challenging problems in subwavelength physics.

\begin{comment}
\newpage
\section*{PLAN}

\begin{enumerate}
\item Introduction
\begin{enumerate}
\item Subwavelength metamaterials
\begin{itemize}
\item summarize the applications
\item metamaterials, cochlea, sensing, feasibility of air bubbles... \cite{ammari2019subwavelength, ammari2020mimicking, leroy2009design, leroy2016design}
\end{itemize}
\item Mathematical background
\begin{itemize}
\item Summary of other (possibly less rigorous or less useful) methods (e.g. scattering matrix, tight-binding approximations, coupled-mode approximations, perturbation theory, homogenization), why our novel approach is better
\end{itemize}
\item General functional analytic approach for subwavelength resonators
\begin{itemize}
\item Give an overview in terms of pole-pencil decompositions etc
\item Review similar functional analytic methods for various settings
\item \emph{e.g.}, Helmholtz resonators, plasmonic particles [25, 26], dielectric resonators,...
\item ``In this paper, we will explore these ideas in detail through a  case study of a specific example..."
\end{itemize}
\item Capacitance coefficients
\begin{itemize}
\item Some background on their use in electrostatics \cite{lekner2011capacitance, smolic2020capacitance}
\item Some other papers giving an overview of the capacitance matrix \cite{diaz2011positivity, smolic2020capacitance, herrera2008geometrical}
\end{itemize}
\end{enumerate}

\item Finite systems
\begin{enumerate}
\item Main results (for 3d)
\item Understanding capacitance coefficients
\begin{itemize}
\item Key properties (symmetric, positivity,...)
\item Basis of eigenvectors (in Hermitian case)
\item Dilute approximations [10]
\item Bispherical coordinates for two spheres \cite{ammari2020close,ammari2019double,lekner2011capacitance}
\end{itemize}
\item close-to-touching resonators [12]
\item higher-order approximations
\begin{itemize}
\item e.g. we can find the resonances to higher-order (the imaginary part) [13, 5]
\end{itemize}
\item Two-dimensional models
\item Numerical results
\begin{itemize}
\item multipole expansions (and exact numerical approach), using capacitance matrix as initial guess
\item computing the capacitance matrix numerically
\end{itemize}
\end{enumerate}

\item Periodic structures
\begin{enumerate}
\item Floquet theory
\item Quasi-periodic layer potentials (alpha non zero and for alpha near zero) (add the expansion of G for alpha near zero from the exceptional point paper).
\item Quasi-periodic and periodic capacitance matrices
\item Quasi-periodic non-Hermitian capacitance matrix
\item Bloch eigenmodes and eigenfrequencies for square lattices versus honeycomb lattices.
%				\item non-Hermitian systems
%				\begin{itemize}
%					\item Unidirectional metascreens [11]
%					\item Edge modes [18]
%				\end{itemize}
\item higher-order approximations
\begin{itemize}
\item e.g. we can find the scattered field to higher-order [11]
\end{itemize}
\end{enumerate}
\item Applications of the capacitance matrix in subwavelength physics
\begin{enumerate}
\item Effective meduim theory, double-negative metamaterials [17],
\item Graded metamaterials (cochlea)
\item Sensitivity enhancement using high-order exceptional points
\item Defects in periodic structures (point and line defect, method of fictitious sources  [15, 21])
\item Topological properties of an SSH chain
\begin{itemize}
\item topological phases
\item Hermitian band inversion and localized modes
\item Non-Hermitian band inversion and localized modes
\item Robust wave guides [10, 7, 18]
\end{itemize}
\item BICs and Fano resonances for infinite periodic structure
\item Extraordinary transmission and unidirectional reflection (from the exceptional point paper)
\item Time-modulated metamaterials

\end{enumerate}
\item Concluding remarks
\begin{enumerate}
\item the concept of capacitance matrix  applies to other subwavelength resonance phenomena (plasmonics, high-contrast dielectric particles, Helmholtz resonators);
\item to compare  with quantum theory and Hamiltonian systems
\begin{itemize}
\item slow decay so nearest neighbour approximation not suitable
\item $\C$ takes full interactions into account (``strong interactions'')
\end{itemize}
\end{enumerate}
\end{enumerate}
\end{comment}

\appendix
\section{Abstract capacitance matrix}\label{app:abstract}
In this section, we describe how the capacitance formulation emerges from the structure of a general integral operator $$, which describes a subwavelength resonance problem. We work with the functional analytic approach described in \Cref{sec:approach}. At $δ=0$, we assume that the operator $(ω,0)$ has a characteristic value $ω=0$ of multiplicity $2N$, admitting the following pole-pencil decomposition:
\begin{equation}\label{eq:polepencil}
\A(\omega,0)^{-1} = \frac{K}{\omega^2} + \Rc(\omega), \quad \text{ for } \quad K = \sum_{i=1}^N \langle \Phi_i, \cdot\rangle \Psi_i,
\end{equation}
where $((0,0)) = span{Ψ_j}$, $(^*(0,0)) = span{Φ_j}$ and $$ is holomorphic for $ω$ in a neighbourhood of $0$.
Moreover, we assume that $(ω,δ)$, for small but non-zero $δ$, satisfies
$$\A(\omega,\delta) = \A(\omega,0) + \L(\omega,\delta),$$
for some operator $Ł$ satisfying (in corresponding operator norm) $Ł = O(δ)$ uniformly for $ω$ in a neighbourhood of $0$.

In this abstract setting, we can derive a capacitance formulation of the subwavelength resonances. We solve the equation $(ω,δ)Φ= 0$. Multiplying with $(ω,0)^-1$, we have
\begin{align*}
0 &= \A(\omega,0)^{-1}\A(\omega,\delta)\Phi = \A(\omega,0)^{-1}\left(\A(\omega,0) + \L\right)\Phi = \left(I + \frac{K\L}{\omega^2} + \Rc \L\right)\Phi.
\end{align*}
Defining $ℬ(ω,δ) = ω^2(ω) Ł(ω,δ)$ yields
\begin{equation}\label{eq:eigen}
\left(\omega^2I + K\L + \mathcal{B}\right)\Phi = 0.
\end{equation}
The characteristic values are therefore determined by \eqref{eq:eigen}, which is in general a non-linear eigenvalue problem since $Ł$ and $ℬ$ depend on $ω$. If we consider the subwavelength resonances, we have $ℬ = O(ω^2δ)$ uniformly for $ω$ and $δ$ around $0$. Similarly, we have $Ł= Ł_0 + Ł̂$, where $Ł̂=O(ωδ)$. Therefore, the subwavelength resonances are approximated by the eigenvalues of the finite-rank operator $-KŁ_0$ whose restriction to $((0,0))$ is given by the generalized capacitance matrix:
$$\C_{ij} =  -\langle \Phi_i, \L_0\Psi_j \rangle.$$
Then the characteristic values satisfy
$$\omega_n = \pm \sqrt{\lambda_n} + O(\delta),$$
where $λ_n$ are the eigenvalues of $$.

As an example of this formulation, for a finite collection of $N$ resonators in $d=3$ (as considered in \Cref{sec:finite}), the operator $$ is given by
\begin{equation*}
\A(\omega,\delta)= \begin{pmatrix}
\widetilde{\S}_D^\omega & -\S_D^k \\
-\frac{1}{2}I+\widetilde{\K}_D^{\omega,*} & -\widetilde\delta \left(\frac{1}{2}I+\K_D^{k,*}\right)
\end{pmatrix}.
\end{equation*}
With above notation, we then have that
\begin{equation*}
\A(\omega,0)= \begin{pmatrix}
{\S}_D^0 & -\S_D^0 \\
-\frac{1}{2}I+{\K}_D^{0,*} & 0
\end{pmatrix}, \quad \L(\omega,\delta) = \begin{pmatrix}
0 & 0 \\
0 & -\widetilde\delta \left(\frac{1}{2}I+\K_D^{k,*}\right)
\end{pmatrix}, \quad \L_0(\delta) = \begin{pmatrix}
0 & 0 \\
0 & -\widetilde\delta \left(\frac{1}{2}I+\K_D^{0,*}\right)
\end{pmatrix}.
\end{equation*}
Moreover, it can be shown that $(ω,0)^-1$ satisfies \eqref{eq:polepencil} where
$$\Phi_i = -\frac{v_i^2}{|D_i|}\begin{pmatrix} 0 \\ \chi_{\p D_i} \end{pmatrix},\qquad \Psi_j = \begin{pmatrix} \psi_i \\ \psi_i \end{pmatrix}, \quad \psi_i = (\S_D^0)^{-1}[\chi_{\p D_i}].$$
From this, it is straightforward to compute $_ij=δ_iv_i^2/|D_i|⟨χ_D_i, ψ_j⟩$, as defined in \eqref{eq:GCM} in \Cref{defn:GCM}.

\bibliographystyle{abbrv}
\begin{thebibliography}{10}

[1]
A.~Al{\`u} and N.~Engheta.
\newblock Achieving transparency with plasmonic and metamaterial coatings.
\newblock {\em Phys. Rev. E}, 72(1):016623, 2005.

[2]
H.~Ammari and G.~Bao.
\newblock Mathematical modeling of near-field optics.
\newblock {\em J. Math. Anal. Appl.}, 265(2):430 -- 446, 2002.

[3]
H.~Ammari, A.~Dabrowski, B.~Fitzpatrick, P.~Millien, and M.~Sini.
\newblock Subwavelength resonant dielectric nanoparticles with high refractive
\newblock {\em Math. Meth. Appl. Sci.}, 42(18):6567--6579, 2019.

[4]
H.~Ammari and B.~Davies.
\newblock A fully coupled subwavelength resonance approach to filtering
auditory signals.
\newblock {\em Proc. R. Soc. A}, 475(2228):20190049, 2019.

[5]
H.~Ammari and B.~Davies.
\newblock A biomimetic basis for auditory processing and the perception of
natural sounds.
\newblock {\em arXiv preprint arXiv:2005.12794}, 2020.

[6]
H.~Ammari and B.~Davies.
\newblock Mimicking the active cochlea with a fluid-coupled array of
subwavelength {Hopf} resonators.
\newblock {\em Proc. R. Soc. A}, 476(2234):20190870, 2020.

[7]
H.~Ammari, B.~Davies, and E.~O. Hiltunen.
\newblock Robust edge modes in dislocated systems of subwavelength resonators.
\newblock {\em arXiv preprint arXiv:2001.10455}, 2020.

[8]
H.~Ammari, B.~Davies, E.~O. Hiltunen, H.~Lee, and S.~Yu.
\newblock Bound states in the continuum and {Fano} resonances in subwavelength
resonator arrays.
\newblock {\em arXiv preprint arXiv:2103.12470}, 2021.

[9]
H.~Ammari, B.~Davies, E.~O. Hiltunen, H.~Lee, and S.~Yu.
\newblock High-order exceptional points and enhanced sensing in subwavelength
resonator arrays.
\newblock {\em Stud. Appl. Math.}, 146(2):440--462, 2021.

[10]
H.~Ammari, B.~Davies, E.~O. Hiltunen, and S.~Yu.
\newblock Topologically protected edge modes in one-dimensional chains of
subwavelength resonators.
\newblock {\em J. Math. Pures Appl.}, 144:17--49, 2020.

[11]
H.~Ammari, B.~Davies, H.~Lee, E.~O. Hiltunen, and S.~Yu.
\newblock Exceptional points in parity--time-symmetric subwavelength
\newblock {\em arXiv preprint arXiv:2003.07796}, 2020.

[12]
H.~Ammari, B.~Davies, and S.~Yu.
\newblock Close-to-touching acoustic subwavelength resonators: eigenfrequency
separation and gradient blow-up.
\newblock {\em Multiscale Model. Simul.}, 18(3):1299–1317, 2020.

[13]
H.~Ammari, B.~Fitzpatrick, D.~Gontier, H.~Lee, and H.~Zhang.
\newblock Minnaert resonances for acoustic waves in bubbly media.
\newblock {\em Ann. I. H. Poincar{\'e}-An.}, 35(7):1975--1998, 2018.

[14]
H.~Ammari, B.~Fitzpatrick, E.~O. Hiltunen, H.~Lee, and S.~Yu.
\newblock Honeycomb-lattice {M}innaert bubbles.
\newblock {\em SIAM J. Math. Anal.}, 52(6):5441--5466, 2020.

[15]
H.~Ammari, B.~Fitzpatrick, E.~O. Hiltunen, and S.~Yu.
\newblock Subwavelength localized modes for acoustic waves in bubbly crystals
with a defect.
\newblock {\em SIAM J. Appl. Math.}, 78(6):3316--3335, 2018.

[16]
H.~Ammari, B.~Fitzpatrick, H.~Kang, M.~Ruiz, S.~Yu, and H.~Zhang.
\newblock {\em Mathematical and Computational Methods in Photonics and
Phononics}, volume 235 of {\em Mathematical Surveys and Monographs}.
\newblock American Mathematical Society, Providence, 2018.

[17]
H.~Ammari, B.~Fitzpatrick, H.~Lee, S.~Yu, and H.~Zhang.
\newblock Double-negative acoustic metamaterials.
\newblock {\em Quart. Appl. Math.}, 77(4):767--791, 2019.

[18]
H.~Ammari and E.~O. Hiltunen.
\newblock Edge modes in active systems of subwavelength resonators.
\newblock {\em arXiv preprint arXiv:2006.05719}, 2020.

[19]
H.~Ammari and E.~O. Hiltunen.
\newblock Time-dependent high-contrast subwavelength resonators.
\newblock {\em arXiv preprint arXiv:2012.10274}, 2020.

[20]
H.~Ammari, E.~O. Hiltunen, and S.~Yu.
\newblock A high-frequency homogenization approach near the {D}irac points in
bubbly honeycomb crystals.
\newblock {\em Arch. Ration. Mech. Anal.}, 238(3):1559--1583, 2020.

[21]
H.~Ammari, E.~O. Hiltunen, and S.~Yu.
\newblock Subwavelength guided modes for acoustic waves in bubbly crystals with
a line defect.
\newblock {\em To appear in J. Eur. Math. Soc.}, 2021.

[22]
H.~Ammari, H.~Kang, and H.~Lee.
\newblock {\em Layer Potential Techniques in Spectral Analysis}, volume 153 of
{\em Mathematical Surveys and Monographs}.
\newblock American Mathematical Society, Providence, 2009.

[23]
H.~Ammari, H.~Lee, and H.~Zhang.
\newblock Bloch waves in bubbly crystal near the first band gap: a
high-frequency homogenization approach.
\newblock {\em SIAM J. Math. Anal.}, 51(1):45--59, 2019.

[24]
H.~Ammari, B.~Li, and J.~Zou.
\newblock Mathematical analysis of electromagnetic scattering by dielectric
nanoparticles with high refractive indices.
\newblock {\em arXiv preprint arXiv:2003.10223}, 2020.

[25]
H.~Ammari, P.~Millien, M.~Ruiz, and H.~Zhang.
\newblock Mathematical analysis of plasmonic nanoparticles: the scalar case.
\newblock {\em Arch. Rational Mech. Anal.}, 224(2):597--658, 2017.

[26]
H.~Ammari, M.~Ruiz, S.~Yu, and H.~Zhang.
\newblock Mathematical analysis of plasmonic resonances for nanoparticles: the
full {Maxwell} equations.
\newblock {\em J. Differ. Equations}, 261(6):3615--3669, 2016.

[27]
H.~Ammari and H.~Zhang.
\newblock A mathematical theory of super-resolution by using a system of
sub-wavelength {Helmholtz} resonators.
\newblock {\em Commun. Math. Phys.}, 337(1):379--428, 2015.

[28]
H.~Ammari and H.~Zhang.
\newblock Effective medium theory for acoustic waves in bubbly fluids near
{Minnaert} resonant frequency.
\newblock {\em SIAM J. Math. Anal.}, 49(4):3252--3276, 2017.

[29]
G.~Bao, D.~C. Dobson, and J.~A. Cox.
\newblock Mathematical studies in rigorous grating theory.
\newblock {\em J. Opt. Soc. Amer.}, A 12(5):1029--1042, 1995.

[30]
L.~G. Bennetts, M.~A. Peter, and R.~V. Craster.
\newblock Graded resonator arrays for spatial frequency separation and
amplification of water waves.
\newblock {\em J. Fluid Mech.}, 854:R4, 2018.

[31]
A.-S. Bonnet-Bendhia and F.~Starling.
\newblock Guided waves by electromagnetic gratings and non-uniqueness examples
for the diffraction problem.
\newblock {\em Math. Method. Appl. Sci.}, 17(5):305--338, 1994.

[32]
L.~C. Botten, M.~Cadilhac, G.~Derrick, D.~Maystre, R.~McPhedran,
M.~Nevi{\`e}re, and P.~Vincent.
\newblock {\em Electromagnetic Theory of Gratings}, volume~22.
\newblock Springer Science \& Business Media, 2013.

[33]
\newblock On the scattered field generated by a ball inhomogeneity of constant
\newblock {\em Asymptotic Anal.}, 77(3-4):197--246, 2012.

[34]
D.~Colton and R.~Kress.
\newblock {\em Integral Equation Methods in Scattering Theory}.
\newblock Wiley, New York, 1983.

[35]
M.~Costabel, E.~Darrigrand, and E.-H. Kon{\'e}.
\newblock Volume and surface integral equations for electromagnetic scattering
by a dielectric body.
\newblock {\em J. Comput. Appl. Math.}, 234(6):1817--1825, 2010.

[36]
M.~Costabel, E.~Darrigrand, and H.~Sakly.
\newblock Volume integral equations for electromagnetic scattering in two
\newblock {\em Comput. Math. Appl.}, 70(8):2087--2101, 2015.

[37]
R.~A. Diaz and W.~J. Herrera.
\newblock The positivity and other properties of the matrix of capacitance:
Physical and mathematical implications.
\newblock {\em J. Electrostat.}, 69(6):587--595, 2011.

[38]
\newblock The bulk-edge correspondence for continuous honeycomb lattices.
\newblock {\em Commun. Part. Diff. Eq.}, 44(12):1406--1430, 2019.

[39]
S.~Dyatlov and M.~Zworski.
\newblock {\em Mathematical Theory of Scattering Resonances}, volume 200 of
{\em Graduate Studies in Mathematics}.
\newblock American Mathematical Society, Providence, RI, 2019.

[40]
C.~L. Fefferman, J.~P. Lee-Thorp, and M.~I. Weinstein.
\newblock Topologically protected states in one-dimensional continuous systems
and {D}irac points.
\newblock {\em Proc. Natl. Acad. Sci. USA}, 111(24):8759--8763, 2014.

[41]
C.~L. Fefferman, J.~P. Lee-Thorp, and M.~I. Weinstein.
\newblock Honeycomb {S}chrödinger operators in the strong binding regime.
\newblock {\em Commun. Pure Appl. Math.}, 71(6):1178--1270, 2018.

[42]
I.~Gohberg and E.~Sigal.
\newblock An operator generalization of the logarithmic residue theorem and the
theorem of {R}ouch\'{e}.
\newblock {\em Sb. Math.}, 13(4):603--625, 1971.

[43]
\newblock Resonance theory in atom-surface scattering.
\newblock {\em Comm. Math. Phys.}, 126(2):263--290, 1989.

[44]
D.~J. Hansen, C.~Poignard, and M.~S. Vogelius.
\newblock Asymptotically precise norm estimates of scattering from a small
circular inhomogeneity.
\newblock {\em Applicable Anal.}, 86(4):433--458, 2007.

[45]
H.~Hodaei, A.~U. Hassan, S.~Wittek, H.~Garcia-Gracia, R.~El-Ganainy, D.~N.
Christodoulides, and M.~Khajavikhan.
\newblock Enhanced sensitivity at higher-order exceptional points.
\newblock {\em Nature}, 548(7666):187--191, 2017.

[46]
M.~S. Jang and H.~Atwater.
\newblock Plasmonic rainbow trapping structures for light localization and
spectrum splitting.
\newblock {\em Phys. Rev. Lett.}, 107(20):207401, 2011.

[47]
N.~Jim{\'e}nez, V.~Romero-Garcia, V.~Pagneux, and J.-P. Groby.
\newblock Rainbow-trapping absorbers: Broadband, perfect and asymmetric sound
absorption by subwavelength panels for transmission problems.
\newblock {\em Sci. Rep.}, 7(1):13595, 2017.

[48]
M.~Kadic, G.~W. Milton, M.~van Hecke, and M.~Wegener.
\newblock {3D} metamaterials.
\newblock {\em Nat. Rev. Phys.}, 1(3):198--210, 2019.

[49]
A.~Karlos and S.~J. Elliott.
\newblock Cochlea-inspired design of an acoustic rainbow sensor with a smoothly
varying frequency response.
\newblock {\em Sci. Rep.}, 10(1):1--11, 2020.

[50]
\newblock {\em {Floquet Theory for Partial Differential Equations}}.
\newblock Number~60 in Operator Theory: Advances and Applications.
{Birkh\"auser Verlag}, Basel, 1993.

[51]
F.~Lemoult, N.~Kaina, M.~Fink, and G.~Lerosey.
\newblock Soda cans metamaterial: A subwavelength-scaled phononic crystal.
\newblock {\em Crystals}, 6(7), 2016.

[52]
J.~Lin and H.~Zhang.
\newblock Mathematical theory for topological photonic materials in one
\newblock {\em arXiv preprint arXiv:2101.05966}, 2021.

[53]
P.~A. Martin.
\newblock {\em Multiple scattering: interaction of time-harmonic waves with N
\newblock Number 107. Cambridge University Press, 2006.

[54]
J.~C. Maxwell.
\newblock {\em A Treatise on Electricity and Magnetism}, volume~1.
\newblock Oxford: Clarendon Press, 1873.

[55]
T.~Meklachi, S.~Moskow, and J.~C. Schotland.
\newblock Asymptotic analysis of resonances of small volume high contrast
linear and nonlinear scatterers.
\newblock {\em J. Math. Phys.}, 59(8):083502, 2018.

[56]
G.~W. Milton and N.-A.~P. Nicorovici.
\newblock On the cloaking effects associated with anomalous localized
\newblock {\em Proc. R. Soc. A.}, 462(2074):3027--3059, 2006.

[57]
\newblock {O}n musical air-bubbles and the sounds of running water.
\newblock {\em Philos. Mag.}, 16(104):235--248, 1933.

[58]
J.~B. Pendry.
\newblock Negative refraction makes a perfect lens.
\newblock {\em Phys. Rev. Lett.}, 85(18):3966, 2000.

[59]
M.~Rupin, G.~Lerosey, J.~de~Rosny, and F.~Lemoult.
\newblock Mimicking the cochlea with an active acoustic metamaterial.
\newblock {\em New J. Phys.}, 21:093012, 2019.

[60]
D.~R. Smith, J.~B. Pendry, and M.~C. Wiltshire.
\newblock Metamaterials and negative refractive index.
\newblock {\em Science}, 305(5685):788--792, 2004.

[61]
W.~P. Su, J.~R. Schrieffer, and A.~J. Heeger.
\newblock Solitons in polyacetylene.
\newblock {\em Phys. Rev. Lett.}, 42:1698--1701, Jun 1979.

[62]
K.~L. Tsakmakidis, A.~D. Boardman, and O.~Hess.
\newblock ‘{Trapped} rainbow’ storage of light in metamaterials.
\newblock {\em Nature}, 450:397--401, 2007.

[63]
V.~G. Veselago.
\newblock The electrodynamics of substances with simultaneously negative values
of $\epsilon$ and $\mu$.
\newblock {\em Sov. Phys. Usp.}, 10(4):509--514, 1968.

[64]
F.~Vollmer, S.~Arnold, and D.~Keng.
\newblock Single virus detection from the reactive shift of a
whispering-gallery mode.
\newblock {\em Proc. Natl. Acad. Sci. U.S.A.}, 105(52):20701--20704, 2008.

[65]
P.~R. Wallace.
\newblock The band theory of graphite.
\newblock {\em Phys. Rev.}, 71(9):622, 1947.

[66]
\newblock Sensors operating at exceptional points: general theory.
\newblock {\em Phys. Rev. A}, 93(3):033809, 2016.

[67]
S.~Yves, R.~Fleury, T.~Berthelot, M.~Fink, F.~Lemoult, and G.~Lerosey.
\newblock Crystalline metamaterials for topological properties at subwavelength
\newblock {\em Nat. Commun.}, 8(1):16023, Jul 2017.

[68]
S.~Zhang, X.~Zhang, L.~Jin, and Z.~Song.
\newblock High-order exceptional points in supersymmetric arrays.
\newblock {\em Phys. Rev. A}, 101(3):033820, 2020.

[69]
J.~Zhu, Y.~Chen, X.~Zhu, F.~J. Garcia-Vidal, X.~Yin, W.~Zhang, and X.~Zhang.
\newblock Acoustic rainbow trapping.
\newblock {\em Sci. Rep.}, 3:1728, 2013.

\end{thebibliography}

\end{document}