# Complex Analytic Dependence on the Dielectric Permittivity in ENZ Materials:
The Photonic Doping Example

Robert V. Kohn and Raghavendra Venkatraman  Courant Institute of Mathematical
Sciences  kohn<EMAIL_ADDRESS>
###### Abstract.

Motivated by the physics literature on “photonic doping” of scatterers made
from “epsilon-near-zero” (ENZ) materials, we consider how the scattering of
time-harmonic TM electromagnetic waves by a cylindrical ENZ region
$\Omega\times\mathbb{R}$ is affected by the presence of a “dopant”
$D\subset\Omega$ in which the dielectric permittivity is not near zero.
Mathematically, this reduces to analysis of a 2D Helmholtz equation
$\mathrm{div}\,(a(x)\nabla u)+k^{2}u=f$ with a piecewise-constant, complex
valued coefficient $a$ that is nearly infinite (say $a=\frac{1}{\delta}$ with
$\delta\approx 0$) in $\Omega\setminus\overline{D}.$ We show (under suitable
hypotheses) that the solution $u$ depends analytically on $\delta$ near $0$,
and we give a simple PDE characterization of the terms in its Taylor
expansion. For the application to photonic doping, it is the leading-order
corrections in $\delta$ that are most interesting: they explain why photonic
doping is only mildly affected by the presence of losses, and why it is seen
even at frequencies where the dielectric permittivity is merely small. Equally
important: our results include a PDE characterization of the leading-order
electric field in the ENZ region as $\delta\to 0$, whereas the existing
literature on photonic doping provides only the leading-order magnetic field.

##  1\. Introduction

A body of literature has developed concerning the design of electromagnetic
devices using epsilon-near-zero (ENZ) materials [10, 13, 18]. A significant
part of it focuses on the transverse magnetic setting, where the magnetic
field has the form $\vec{H}(x)=(0,0,H(x_{1},x_{2}))$ and the time-harmonic
Maxwell’s equations reduce to the two-dimensional equation

$\displaystyle-\nabla\cdot\left(\frac{1}{\varepsilon}\nabla
H\right)-\omega^{2}\mu H=f.$ (1.1)

Here $\varepsilon=\varepsilon(x_{1},x_{2})$ is the dielectric permittivity and
$\mu=\mu(x_{1},x_{2})$ is the magnetic permeability, both of which are in
general complex-valued and frequency-dependent (i.e.,
$\varepsilon=\varepsilon(\omega)$ and $\mu=\mu(\omega)$, where $\omega$
denotes frequency). Within this time-harmonic setting, the electric field is
described by the gradient of the scalar function $H$ through
$\vec{E}=\frac{1}{i\omega\varepsilon}(-\partial_{x_{2}}H,\partial_{x_{1}}H,0).$
The function $f$ models a source that can be quite general.

It is easy to see from (1.1) what is special about ENZ materials: since it is
natural to expect that $\frac{1}{\varepsilon}\nabla H$ is bounded, we expect
$H$ to be nearly constant in a region where $\varepsilon$ is nearly $0$. In
the limit $\varepsilon\rightarrow 0$, we are not solving a PDE in the ENZ
region but rather selecting the constant value of $H$. This creates curious
effects; for example, parts of the ENZ region that are distant from one
another are nevertheless tightly coupled, since they share the same (constant)
value of $H$; moreover the shape of the ENZ region is less important than it
would be for a PDE, since all that is happening there is the selection of a
constant $H$. This effect has for example been used (i) to design entirely new
types of waveguides [1, 19, 20, 21, 24], and (ii) to show how “dopants” can
give an ENZ region an “effective magnetic permeability” $\mu_{\mathrm{eff}}$
that is entirely different from its true magnetic permeability $\mu$ [15, 20,
24].

The preceding discussion assumes that the PDE (1.1) has a meaningful limit as
$\varepsilon\rightarrow 0$ in the region occupied by the ENZ material. We
shall prove this, but we are mainly interested in the corrections that appear
when $\varepsilon$ is _nonzero_. Making a choice, we shall focus on the second
application mentioned above– photonic doping– though the strategy we follow
should also be applicable to waveguides.

To state the main results of this paper, we begin by describing the photonic
doping problem that we study. The scatterer has cross-section
$\Omega\subset\mathbb{R}^{2}$, with finitely many cylindrical dopant rods. We
let $D\subset\Omega$ denote the union of the dopants (see Fig. 1), and
consider the Helmholtz equation (1.1) with piecewise constant $\varepsilon$
and $\mu$, with $\varepsilon$ small in magnitude in the ENZ region
$\Omega\setminus\overline{D}$. Our analysis is sufficiently general to permit
$\varepsilon$ and $\mu$ taking distinct constant complex values in the
exterior $\mathbb{R}^{2}\setminus\overline{\Omega},$ the ENZ region
$\Omega\setminus\overline{D},$ and each of the dopant rods, respectively. So
long as the $|\varepsilon|\ll 1$ in the ENZ region, the mathematical analysis
of the problem is not complicated by the value of these other constants, as
well as the (finite) number of connected components of the dopant $D$.

Figure 1. The photonic doping problem

However, for simplicity, and in order to focus on the mathematical essence, we
make some reductions. We assume that there is a single dopant rod $D$ (see
Fig. 2) immersed in the ENZ region; we take $\varepsilon=1$ in the exterior;
and we take $\varepsilon=\delta\in\mathbb{C}\setminus\\{0\\},|\delta|\ll 1$ in
the ENZ region. Also, to simplify the presentation, we assume that
$\varepsilon=1$ in the dopant. Thus, we shall study the case

$\displaystyle\varepsilon_{\delta}(x):=\left\\{\begin{array}[]{cc}1&x\in\mathbb{R}^{2}\setminus\overline{\Omega},\\\
\delta&x\in\Omega\setminus\overline{D},\\\ 1&x\in D.\end{array}\right.$ (1.5)

We shall assume that the dopant is not resonant (see Section 2.2 for a precise
statement, and Section 5 for some discussion about what happens if the dopant
is resonant). Our analysis requires no other condition on $\delta;$ in
particular, it applies equally to the case $\mathrm{Im}(\delta)>0$
(corresponding to a lossy material) and the case $\mathrm{Im}(\delta)<0$
(corresponding to a material with gain) – provided, of course, that $|\delta|$
is sufficiently small. As for the value of $\mu$, we take it to be the same in
all regions (though the case where it takes different values in the three
regions would not be fundamentally different). For the exact assumptions on
the regularity of $\Omega$ and $D$, as well as our assumptions on the source,
we refer the reader to Section 2.

Our main result, stated in Theorem 2.2 below, asserts that the solution to
(1.1) depends analytically on the complex variable $\delta$ in a neighborhood
of $0$. Let us emphasize that the behavior of the magnetic field $\vec{H}$ in
the limit $\delta=0$ is correctly found (though not proved) in the physics
literature [15, 20]. As previously mentioned, this amounts to selecting the
correct constant value of the scalar function $H$. However, the physics
literature offers no analytical understanding about:

  1. (1)

the character and magnitude of corrections due to $\delta\neq 0,$ and

  2. (2)

the character and magnitude of the electric field $\vec{E}$, which, as we show
below, has a nonzero limit as $|\delta|\to 0.$

As corollaries to our analysis, we obtain simple and reasonably explicit
information on both points. Such information is physically relevant. To
explain why, we remind the reader that for real materials, the permittivity
$\varepsilon=\varepsilon^{\prime}+i\varepsilon^{\prime\prime}$ is a complex-
analytic function of the frequency $\omega$; moreover its imaginary part
represents losses (which are always present in real materials), so
$\varepsilon^{\prime\prime}(\omega)>0$ when $\omega$ is real and positive (see
e.g. [12] sections 77 and 80). Thus, while the notion of an epsilon-near-zero
material is a useful idealization, the physically relevant value of our
parameter $\delta$ has positive imaginary part, and the real part of $\delta$
is nonzero except at isolated frequencies. For the effect of “photonic doping”
to be robust with respect to material losses and frequency variation, the
corrections due to nonzero $\delta$ must be small. Our analyticity results
show that these corrections are of order $O(|\delta|)$.

Concerning (2), we recall that the Poynting vector $S$ is defined by
$S:=\frac{1}{2}\vec{E}\times\overline{\vec{H}}$. (Here, as usual, for any
complex number $z,$ $\overline{z}$ denotes its complex conjugate.) The real
part of $S$ is the time-averaged power flow, while the imaginary part is the
reactive power flow. It is of interest to understand the spatial structure of
these flows. Since the Poynting vector involves both the magnetic field
$\vec{H}$ and the electric field $\vec{E},$ this requires understanding both
fields in the limit $\delta\rightarrow 0$. For ENZ waveguides, this was the
focus of [14]; we are not aware of a similar study in the photonic doping
setting.

### Some related literature

We close this introduction with a brief discussion of some related literature.

  1. (1)

The paper [3] by O. Bruno considers the PDE $\mathrm{div}(a(x)\nabla u)=0$ for
coefficient fields $a$ that take two values: the conductivity of the
inclusions is $\delta\in\mathbb{C}\setminus\\{0\\},$ and that of the
background matrix is $1.$ It studies the effective conductivity as a function
of $\delta$ for both periodic and random coefficient fields $a,$ and provides
conditions on the geometry of the inclusions that ensure analyticity of the
effective conductivity near $\delta\approx 0$ and $\delta\approx\infty.$ Our
approach is not fundamentally different, but our geometry is more specific;
moreover, our focus is the Helmholtz equation.

  2. (2)

The papers [2] and [23] study how a periodic array of inclusions can lead to a
frequency-dependent effective $\mu$, in a suitable homogenized limit. In
photonic doping, an effective permeability $\mu_{\mathrm{eff}}$ emerges due to
the presence of the dopant (see Remark 2.3 below). However, our setting is
different from [2, 23] since we do not consider an array of inclusions, and in
particular, no homogenization procedure is being carried out. And yet our
setting shares some features with that of [2, 23], since here and also there,
the effective permeability is directly related to the solution of a Helmholtz
equation in the dopant with boundary condition $1.$

  3. (3)

In the present work, we assume that the ENZ region is a uniform
electromagnetic medium, characterized by two constants, its complex
permittivity and permeability. There are real materials and experimental
systems for which this is the case (see for example [8, 9, 18]). A different
but related body of literature is about metamaterials that achieve an
effective $\varepsilon$ that is close to zero at certain frequencies, as a
consequence of their specific microstructure (see [18], or for a recent
example [16]). In such systems, the robustness of the effect to the variation
of the frequency $\omega$ or the presence of loss depends not only on the
effects studied here, but also on how well their behavior is captured by the
effective $\varepsilon$.

  4. (4)

In our study of photonic doping the parameter $\delta$ tends to zero while the
geometry is fixed. It is of course interesting to ask whether other exotic
effects can be achieved by letting the microstructure vary as a physical
parameter approaches $0$ or $\infty$. The answer is yes; indeed, many papers
on electromagnetic metamaterials have this character. For divergence-form
equations of the form $\mathrm{div}\,(a(x)\nabla u)=0$ with $u$ scalar, the
full set of limits obtainable this way has been characterized [4]. The
analogous question for linear elasticity has also been addressed, and the set
of all possible limits is surprisingly rich [5]. We are not aware of analogous
results for Maxwell’s equations.

### Organization of the paper

After precisely formulating the problem and recording some preliminaries, we
state our main results in Section 2. This section also contains a roadmap of
the proof of the main theorem (see subsection 2.5). We introduce certain
operators that arise in describing higher order correctors in Section 3. The
proof of the main theorem and its corollaries are in Section 4, and we
conclude the paper in Section 5 by considering what happens when the dopant is
resonant (or nearly so).

##  2\. Formulation of the Problem and Main Results

###  2.1. Setting up the problem

We are given a bounded domain $\Omega\subset\mathbb{R}^{2},$ with Lipschitz
boundary (i.e., such that the boundary is locally the graph of a Lipschitz
continuous function), and an open, Lipschitz subdomain $D\Subset\Omega$ that
represents the dopant. Here and throughout the paper, we use the notation
$A\Subset B$ to denote that $\overline{A}\subset B$, where $\overline{A}$ is
the closure of $A$. The ENZ region fills the set
$\Omega\setminus\overline{D}.$ The unit outward pointing normal to
$\partial\Omega$ and $\partial D$ are denoted by $\nu_{\Omega}$ and $\nu_{D}$
respectively.

We let $\delta\in\mathbb{C}\setminus\\{0\\}$ be a complex number, and define
the spatially inhomogeneous but isotropic permittivity function as in (1.5).
Suppose $f$ is a source that is supported compactly outside
$\overline{\Omega}.$ To be precise, we suppose

$\displaystyle f\in L^{2}(K)$

for some compact set $K\Subset\mathbb{R}^{2}\setminus\overline{\Omega}$ (see
Figure 2).

$K$$\Omega$$D$$\nu_{\Omega}$$\nu_{D}$ Figure 2. Set up of the problem: $D$
represents the dopant rod immersed in the ENZ material, which occupies the
open set $\Omega\setminus\overline{D}$. The scattering response of the system
is excited through a source $f$ supported in a compact set
$K\subset\mathbb{R}^{2}\setminus\overline{\Omega}$. The unit outward pointing
normal to the sets $D$ and $\Omega$ are denoted respectively by $\nu_{D}$ and
$\nu_{\Omega}.$

Finally, defining the operator

$\displaystyle\mathcal{L}_{\delta}:=-\nabla\cdot\left(\frac{1}{\varepsilon_{\delta}(x)}\nabla\right)-k^{2},$

where $k$ is a nonzero complex number with $0\leqslant\mathrm{arg}\,k<\pi$, we
consider the problem

$\displaystyle\mathcal{L}_{\delta}u_{\delta}=f,\mbox{ in }\mathbb{R}^{2},$
(2.1) $\displaystyle\lim_{r\to\infty}\sqrt{r}\left(\frac{\partial}{\partial
r}-ik\right)u_{\delta}=0.$ (2.2)

Here, as usual, $r=|x|$ is the distance to the origin. The second condition is
the Sommerfeld radiation condition, which assures uniqueness by selecting a
solution that is “outgoing” at infinity. One interprets the PDE in (2.1) along
$\partial D$ and $\partial\Omega$ by requiring continuity of $u_{\delta}$ and
$\frac{1}{\varepsilon_{\delta}}\frac{\partial u_{\delta}}{\partial\nu}.$ Our
overall goal is to understand the dependence of the solution $u_{\delta}$ on
the parameter $\delta$ near $\delta=0$.

###  2.2. The nonresonance condition and three auxiliary functions

Our analysis requires that the dopant be non-resonant, in other words that
$k^{2}$ not be a Dirichlet eigenvalue of the Laplacian in $D$. Equivalently:

We assume that the Helmholtz operator $-\Delta-k^{2}$ (2.3) has trivial kernel
in $H^{1}_{0}(D)$.

By elliptic theory, this condition holds for all but a countable collection of
$k^{2}\in(0,\infty).$ (For some discussion about what happens when this
condition fails, see Section 5.)

The solutions of the following three auxiliary problems will play key roles in
our analysis:

Problem 1: Let $s:\mathbb{R}^{2}\setminus\overline{\Omega}\to\mathbb{R}$
denote the unique solution to

$\left\\{\begin{aligned} -\Delta s&=k^{2}s+f\quad&\mbox{ in
}\mathbb{R}^{2}\setminus\overline{\Omega},\\\ s&=0\quad&\mbox{ along
}\partial\Omega,\\\ s&\mbox{ satisfies the radiation condition
(\ref{e.bcinfty}) at infinity.}\end{aligned}\right.$ (2.4)

Problem 2: Let $\psi_{e}$ denote the unique solution to

$\left\\{\begin{aligned} -\Delta\psi_{e}&=k^{2}\psi_{e}\quad&\mbox{ in
}\mathbb{R}^{2}\setminus\overline{\Omega},\\\ \psi_{e}&=1\quad&\mbox{ along
}\partial\Omega,\\\ \psi_{e}&\mbox{ satisfies the radiation condition
(\ref{e.bcinfty}) at infinity. }\end{aligned}\right.$ (2.5)

Problem 3: Let $\psi_{d}$ denote the unique solution to

$\left\\{\begin{aligned} -\Delta\psi_{d}&=k^{2}\psi_{d}\quad&\mbox{ in
}{D},\\\ \psi_{d}&=1\quad&\mbox{ along }\partial D.\end{aligned}\right.$ (2.6)

The existence and uniqueness of solutions to Problems 1 and 2 is well-known,
see for example [17, Theorem 9.11]. The existence and uniqueness of
$\psi_{d}\in H^{1}(D)$ follows by Fredholm theory from our nonresonance
hypothesis (2.3).

Our analysis will require dividing by the complex number

$\beta:=k^{2}|\Omega\setminus\overline{D}|+\int_{\partial\Omega}\frac{\partial\psi_{e}}{\partial\nu_{\Omega}}\,d\mathcal{H}^{1}-\int_{\partial
D}\frac{\partial\psi_{d}}{\partial\nu_{D}}\,d\mathcal{H}^{1}.$ (2.7)

(Our sign conventions for the normal derivatives are shown in Fig. 2:
$\nu_{\Omega}$ denotes the outward pointing unit normal to $\partial\Omega,$
and $\nu_{D}$ denotes the outward pointing unit normal to $\partial D.$) It is
therefore important to know that this never vanishes.

######  Lemma 2.1.

The complex number $\beta$ defined by (2.7) is always nonzero.

###### Proof.

A result sometimes known as Rellich’s lemma says that if $u$ solves $\Delta
u+k^{2}u=0$ in $\mathbb{R}^{2}\setminus\overline{\Omega}$ and satisfies the
radiation condition (2.2) at infinity, then

$-2\,\mbox{Im}\left(k\int_{\partial\Omega}u\frac{\partial\overline{u}}{\partial\nu_{\Omega}}\,d\mathcal{H}^{1}\right)=\\\
2\,\mbox{Im}(k)\int_{\mathbb{R}^{2}\setminus\overline{\Omega}}(|k|^{2}|u|^{2}+|\nabla
u|^{2})\,dx+\lim_{R\rightarrow\infty}\int_{|x|=R}(|\partial u/\partial
r|^{2}+k^{2}|u|^{2})\,d\mathcal{H}^{1}$

(see e.g. equation (3.10) of [6]); moreover the second term on the right can
be written in terms of the far-field pattern of $u$, and it never vanishes
(unless $u=0$). Specializing to $u=\psi_{e}$, we have

$-2\,\mbox{Im}\left(k\int_{\partial\Omega}\frac{\partial\overline{\psi}_{e}}{\partial\nu_{\Omega}}\,d\mathcal{H}^{1}\right)=\\\
2\,\mbox{Im}(k)\int_{\mathbb{R}^{2}\setminus\overline{\Omega}}(|k|^{2}|\psi_{e}|^{2}+|\nabla\psi_{e}|^{2})\,dx+\lim_{R\rightarrow\infty}\int_{|x|=R}(|\partial\psi_{e}/\partial
r|^{2}+k^{2}|\psi_{e}|^{2})\,d\mathcal{H}^{1}.$ (2.8)

The function $\psi_{d}$ satisfies a similar relation: multiplying its PDE by
$\overline{\psi}_{d}$ and integrating by parts gives

$\int_{\partial
D}\frac{\partial\psi_{d}}{\partial\nu_{D}}\,d\mathcal{H}^{1}=\int_{D}(|\nabla\psi_{d}|^{2}-k^{2}|\psi_{d}|^{2})\,dx.$

Taking the complex conjugate then multiplying by $k$ and taking the imaginary
part gives

$\mbox{Im}\left(k\int_{\partial
D}\frac{\partial\overline{\psi}_{d}}{\partial\nu_{D}}\,d\mathcal{H}^{1}\right)=\mbox{Im}(k)\int_{D}(|\nabla\psi_{d}|^{2}+|k|^{2}|\psi_{d}|^{2})\,dx.$
(2.9)

We turn now to the assertion of the lemma. If $k$ is real then $\psi_{d}$ is
real, and it follows from (2.8) that the imaginary part of
$\int_{\partial\Omega}\partial\psi_{e}/\partial\nu_{\Omega}\,d\mathcal{H}^{1}$
is nonzero. Since the other terms in (2.7) are real, $\beta\neq 0$. If $k$ is
not real then the imaginary part of $k$ is positive (since
$0\leqslant\mathrm{arg}\,k<\pi$). From (2.7) we have

$k\overline{\beta}=\overline{k}\,|k|^{2}|\Omega\setminus\overline{D}|+k\int_{\partial\Omega}\frac{\partial\overline{\psi}_{e}}{\partial\nu_{\Omega}}\,d\mathcal{H}^{1}-k\int_{\partial
D}\frac{\partial\overline{\psi}_{d}}{\partial\nu_{D}}\,d\mathcal{H}^{1}.$

Equations (2.8) and (2.9) show that the imaginary part of the right hand side
is strictly negative, from which it follows that $\beta\neq 0$. ∎

###  2.3. Notation, conventions, and some function spaces

Unless otherwise specified, all spaces of functions in the sequel consist of
complex-valued functions. All PDEs are to be interpreted in the weak sense, by
testing with an $H^{1}$ function.

We will frequently need to consider the restriction to $\partial\Omega$ of a
function in $H^{1}(\Omega)$, so let us recall the facts we’ll be using. As the
domain $\Omega$ is bounded with Lipschitz boundary, each function $f\in
H^{1}(\Omega)$ admits a trace $\mathrm{tr}(f)\in L^{2}(\partial\Omega).$ This
trace map is a bounded linear functional on $H^{1}(\Omega)$, whose kernel is
the space $H^{1}_{0}(\Omega),$ and range is the fractional Sobolev space
$H^{1/2}(\partial\Omega)$ (see [17, Theorem 3.38]). The space
$H^{1/2}(\partial\Omega)$ is a Hilbert space with respect to the norm defined
variationally via: for any $g\in H^{1/2}(\partial\Omega),$

$\displaystyle\|g\|_{H^{1/2}(\partial\Omega)}:=\inf\left\\{\|\widetilde{G}\|_{H^{1}(\Omega)}:G\in
H^{1}(\Omega)\mbox{ and }\mathrm{tr}\,G=g\right\\}.$ (2.10)

For a proof of the assertion that this definition is equivalent to other
definitions (say via integral kernels or by locally flattening the boundary,
using a partition of unity, and appealing to a Fourier based definition within
each flattened patch) of the fractional Sobolev space
$H^{1/2}(\partial\Omega)$ within the setting of Lipschitz domains, we refer
the reader to [11, Section 5.1]. Finally, we define $H^{-1/2}(\partial\Omega)$
to be the dual space of $H^{1/2}(\partial\Omega).$ It is equivalently
characterized as the space of traces of normal components $V\cdot\nu_{\Omega}$
of measurable vector fields $V$ on $\Omega$ such that
$V\in(L^{2}(\Omega))^{2},$ and the divergence of $V$ is square integrable,
i.e., $\mathrm{div}\,V\in L^{2}(\Omega)$ (see [22, Chapter 20] for a proof).

An application of this latter equivalence that we will use repeatedly, is an
estimate of the normal derivative of $s$ defined via (2.4). We let $\Psi\in
C^{1}_{c}(\mathbb{R}^{2}\setminus\overline{\Omega})$ denote a cut off function
that satisfies $\Psi\equiv 1$ in a neighborhood of $\Omega.$ Considering the
vector field $V:=\Psi\nabla s,$ from (2.4) it is clear that
$V\in(L^{2}(\mathbb{R}^{2}\setminus\overline{\Omega}))^{2},$ with
$\mathrm{div}\,V\in L^{2}(\mathbb{R}^{2}\setminus\overline{\Omega}),$ has
compact support, and finally, satisfies $V\cdot\nu=\frac{\partial
s}{\partial\nu}.$ Therefore, by the above characterization of
$H^{-1/2}(\partial\Omega),$ it follows that

$\displaystyle\left\|\frac{\partial
s}{\partial\nu_{\Omega}}\right\|_{H^{-1/2}(\partial\Omega)}$
$\displaystyle=\|V\cdot\nu_{\Omega}\|_{H^{-1/2}(\partial\Omega)}$ (2.11)
$\displaystyle\leqslant C(\|V\|_{L^{2}}+\|\mathrm{div}\,V\|_{L^{2}})\leqslant
C\|f\|_{L^{2}(\mathbb{R}^{2}\setminus\overline{\Omega})},$

using the estimate

$\displaystyle\int_{\mathrm{spt}\,\Psi}|s|^{2}+|\nabla s|^{2}\,dx\leqslant
C\int_{\mathbb{R}^{2}\setminus\overline{\Omega}}|f|^{2}\,dx,$

which is proved in showing the existence and uniqueness of $s$ (see for
example [17]).

###  2.4. A summary of our main results

Our main theorem is

######  Theorem 2.2.

Suppose the dopant is non-resonant, in other words that (2.3) holds. Then
there exists $\Lambda>0$, depending only on $k$ and the sets $\Omega$ and $D$
(and independent of the source $f$ and the value of $\delta$) such that if
$|\delta|<\frac{1}{\Lambda},$ then the map $\delta\mapsto u_{\delta}\in
H^{1}_{loc}(\mathbb{R}^{2};\mathbb{C})$ is complex analytic. To be precise,
there exist $\\{v_{j}\\}_{j=0}^{\infty}\subset H^{1}_{loc}(\mathbb{R}^{2}),$
independent of $\delta$, such that

$\displaystyle u_{\delta}=\sum_{j=0}^{\infty}\delta^{j}v_{j}\quad\quad\mbox{
in }H^{1}_{loc}(\mathbb{R}^{2}),$ (2.12)

with $v_{j}$ satisfying: for each open set $U\subset\mathbb{R}^{2},$ there
exists a constant $K_{U}>0$ (independent of $\delta$), such that

$\displaystyle\|v_{j}\|_{H^{1}(U)}\leqslant
K_{U}\Lambda^{j}\left\|f\right\|_{L^{2}}.$

As will be clear from the proof, we obtain a complete characterization of all
the “correctors” $v_{j}$ that appear in the statement of Theorem 2.2. But it
is the leading-order correctors that are the most important, and the easiest
to describe. Therefore we devote the remainder of this section to a self-
contained description of the leading-order correctors, and statements of some
corollaries that they permit us to draw.

The zeroth order term in (2.12) is of course the easiest to describe: it is

$\displaystyle
v_{0}(x):=\left\\{\begin{array}[]{cc}c^{*}\psi_{e}+s&x\in\mathbb{R}^{2}\setminus\overline{\Omega}\\\
c^{*}&x\in\Omega\setminus\overline{D}\\\ c^{*}\psi_{d}&x\in
D,\end{array}\right.$

where

$c^{*}:=-\frac{1}{\beta}\int_{\partial\Omega}\frac{\partial
s}{\partial\nu_{\Omega}}\,d\mathcal{H}^{1},$ (2.13)

which is finite since $\beta\neq 0$ by Lemma 2.1.

######  Remark 2.3.

As noted in the introduction, $v_{0}$ is constant in the ENZ region. In the
physics literature, the term “effective permeability” $\mu_{\mathrm{eff}}$ is
used for the choice of $\mu$ that would have given the same value to $v_{0}$
in the ENZ region without any dopant. In our setting, the value of this choice
is

$\displaystyle\mu_{\mathrm{eff}}:=\frac{|\Omega\setminus\overline{D}|+\int_{D}\psi_{d}\,dx}{|\Omega|}\mu.$
(2.14)

Indeed, the preceding description of $\mu_{\mathrm{eff}}$ reduces to the value
such that

$\displaystyle\omega^{2}\mu_{\mathrm{eff}}|\Omega|+\int_{\partial\Omega}\frac{\partial\psi_{e}}{\partial\nu_{\Omega}}\,d\mathcal{H}^{1}=\omega^{2}\mu|\Omega\setminus\overline{D}|+\int_{\partial\Omega}\frac{\partial\psi_{e}}{\partial\nu_{\Omega}}\,d\mathcal{H}^{1}-\int_{\partial
D}\frac{\partial\psi_{d}}{\partial\nu_{D}}\,d\mathcal{H}^{1}.$ (2.15)

A bit of manipulation using Green’s theorem, the PDE solved by $\psi_{d}$, our
convention that $k^{2}=\omega^{2}\mu$, and

$-\int_{\partial
D}\frac{\partial\psi_{d}}{\partial\nu_{D}}\,d\mathcal{H}^{1}=k^{2}\int_{D}\psi_{d}\,dx,$

reveals that (2.14) and (2.15) are equivalent.

To identify the leading order electric field
$\frac{1}{i\omega\varepsilon(x)}(-\partial_{x_{2}}u_{\delta},\partial_{x_{1}}u_{\delta},0)$
we must discuss the $j=1$ term in (2.12). To do so, we introduce new
correctors. We let $\phi_{0}\in H^{1}(\Omega\setminus\overline{D})$ denote the
unique weak solution to the PDE

$\left\\{\begin{aligned} -\Delta\phi_{0}&=k^{2}c^{*}\quad\quad\mbox{ in
}\Omega\setminus\overline{D}\\\
\frac{\partial\phi_{0}}{\partial\nu_{\Omega}}&=c^{*}\frac{\partial\psi_{e}}{\partial\nu_{\Omega}}+\frac{\partial
s}{\partial\nu_{\Omega}}\quad\quad\mbox{ on }\partial\Omega\\\
\frac{\partial\phi_{0}}{\partial\nu_{D}}&=c^{*}\frac{\partial\psi_{d}}{\partial\nu_{D}}\quad\quad\mbox{
on }\partial D,\\\
\strokedint_{\Omega\setminus\overline{D}}&\phi_{0}\,dx=0;\end{aligned}\right.$
(2.16)

and we define the exterior corrector $\lambda_{0}\in
H^{1}(\mathbb{R}^{2}\setminus\overline{\Omega}),$

$\left\\{\begin{aligned} -\Delta\lambda_{0}&=k^{2}\lambda_{0}\quad\quad\mbox{
in }\mathbb{R}^{2}\setminus\overline{\Omega},\\\
\lambda_{0}&=\phi_{0}\quad\quad\mbox{ on }\partial\Omega,\\\
\lambda_{0}&\mbox{ satisfies the radiation condition (\ref{e.bcinfty}) at
infinity},\end{aligned}\right.$ (2.17)

and the dopant corrector $\chi_{0}\in H^{1}(D),$ given by

$\left\\{\begin{aligned} -\Delta\chi_{0}&=k^{2}\chi_{0}\quad\quad\mbox{ in
}D,\\\ \chi_{0}&=\phi_{0}\quad\quad\mbox{ on }\partial D.\end{aligned}\right.$
(2.18)

We emphasize that the corrector $\phi_{0}$ in the ENZ region
$\Omega\setminus\overline{D}$ solves a Poisson equation with constant right
hand side, not a Helmholtz equation. The choice of $c^{*}$ by (2.13) exactly
ensures that the problem (2.16) is well-posed; the problems defining
$(\lambda_{0},\chi_{0})$ are also well-posed, since we have assumed that the
dopant isn’t resonant. Define

$e_{0}:=-\frac{1}{\beta}\left(\int_{\partial\Omega}\frac{\partial\lambda_{0}}{\partial\nu_{\Omega}}\,d\mathcal{H}^{1}-\int_{\partial
D}\frac{\partial\chi_{0}}{\partial\nu_{D}}\,d\mathcal{H}^{1}\right),$ (2.19)

and set $c_{\delta}:=c^{*}+\delta e_{0}.$ The following result captures the
essential information available from the first two terms in (2.12).

######  Corollary 2.4.

Let $|\delta|<\frac{1}{\Lambda},$ where $\Lambda$ is as in Theorem 2.2, and
let $v_{\delta}\in H^{1}_{loc}(\mathbb{R}^{2})$ be defined via

$v_{\delta}(x):=\left\\{\begin{array}[]{ccc}&c_{\delta}\psi_{e}(x)+s(x)+\delta\lambda_{0}(x)&x\in\mathbb{R}^{2}\setminus\overline{\Omega}\\\
&c_{\delta}+\delta\phi_{0}&x\in\Omega\setminus\overline{D}\\\
&c_{\delta}\psi_{d}+\delta\chi_{0}&x\in D.\end{array}\right.$ (2.20)

Then, for any open set $U\subset\mathbb{R}^{2}$ with compact closure, there
exists a constant $C_{U}>0$ independent of $\delta$ such that

$\displaystyle\|u_{\delta}-v_{\delta}\|_{H^{1}(U)}\leqslant C_{U}\delta^{2}.$
(2.21)

Let us note that we obtain local in $H^{1}$ estimates only because of the
unboundedness of the domain. In particular, within the ENZ region
$\Omega\setminus\overline{D},$ the differences above are global, in the
following sense: the estimate (2.21) says that $u_{\delta}$ is close to being
a constant, where closeness is measured in terms of
$H^{1}(\Omega\setminus\overline{D})$.

######  Remark 2.5.

A consequence of Corollary 2.4 is that at leading order, the electric field
distribution in the ENZ region is given by

$\displaystyle\vec{E}:=\frac{1}{i\omega}(-\partial_{x_{2}}\phi_{0},\partial_{x_{1}}\phi_{0},0).$

Indeed, this follows, since within the TM setting of our paper, the magnetic
field is $\vec{H}=(0,0,u_{\delta}),$ and the electric field is
$\vec{E}=\frac{1}{i\omega\varepsilon_{\delta}(x)}(-\partial_{x_{2}}u_{\delta},\partial_{x_{1}}u_{\delta},0).$
Similarly, from Corollary 2.4, it is also clear that one can read off
information about the electric field distribution in the other regions.

Our next corollary generalizes the so-called “Ideal Fluid Analogy” introduced
in [14]. In that paper, the authors consider the setting in which both the
dielectric permittivity $\varepsilon_{\delta}$ and the magnetic permeability
($k^{2}=\omega^{2}\mu$ in the notation of the introduction), both vanish in
the ENZ region. As explained in the introduction, our next corollary clarifies
and complements the result in two ways: first, of course, it describes
corrections to the leading order theory recalled in the introduction due to
nonzero $\delta$ (e.g., via the presence of losses), and second, it is not
restricted to the case of vanishing magnetic permeability. In order to state
this corollary, we recall the notion of the Poynting vector, as the authors in
[14] formulate their ideal fluid analogy in terms of the Poynting vector. It
is defined as the function
$S_{\delta}:=\frac{1}{2}\vec{E}_{\delta}\times\overline{\vec{H}_{\delta}},$
with $\overline{z}$ denoting the complex conjugate of a complex number
$z\in\mathbb{C}$ . We have the following characterization of the limit of
$S_{\delta}$ as $\delta\to 0:$

######  Corollary 2.6.

In the ENZ set $\Omega\setminus\overline{D},$ the Poynting vector field
$S_{\delta}$ converges in $L^{2}(\Omega\setminus\overline{D})$ strongly to
${S}=\frac{\overline{c^{*}}}{2i\omega}\nabla\phi_{0}$, which satisfies the PDE
system

$\left\\{\begin{aligned}
\nabla\cdot{S}&=-\frac{1}{2i\omega}k^{2}|c^{*}|^{2}=\frac{i\omega\mu}{2}|c^{*}|^{2},\quad\quad\mbox{
in }\Omega\setminus\overline{D},\\\ \nabla\times{S}&=0,\quad\quad\mbox{ in
}\Omega\setminus\overline{D},\\\
\nu_{\Omega}\cdot{S}&=\frac{1}{2i\omega}\overline{c^{*}}\left(c^{*}\frac{\partial\psi_{e}}{\partial\nu_{\Omega}}+\frac{\partial
s}{\partial\nu_{\Omega}}\right),\quad\quad\mbox{ on }\partial\Omega,\\\
\nu_{D}\cdot{S}&=\frac{1}{2i\omega}|c^{*}|^{2}\frac{\partial\psi_{d}}{\partial\nu_{D}}\quad\quad\mbox{
on }\partial D.\end{aligned}\right.$ (2.22)

The “ideal fluid analogy” explored in [14] rests on the observation that if
$\mu=0$ in the ENZ region then $\nabla\cdot S=0$ and $\nabla\times S=0$ there
(so the real and imaginary parts of $S$ are the gradients of harmonic
functions). Corollary 2.6 provides a simple extension to our setting, which
applies for any value of $\mu.$ In particular, the real and imaginary parts of
$S$ in the ENZ region are the gradients of functions satisfying Poisson-type
equations of the form $\Delta W=\mbox{constant},$ where the right-hand-side is
the real or imaginary part of $\tfrac{1}{2}i\omega\mu|c^{*}|^{2}.$

###  2.5. Main ingredients of the proof

The main idea in the proof of Theorem 2.2 and its corollaries is to pursue a
perturbative ansatz of the form

$\displaystyle u_{\delta}=\sum_{j=0}^{\infty}\delta^{j}v_{j},$ (2.23)

and plug it into the PDE (2.1). When $|\delta|\ll 1$, we expect that
$u_{\delta}$ is essentially constant in the ENZ region. To this end, by
linearity of the PDE (2.1), one seeks

$\displaystyle v_{0}:=\left\\{\begin{array}[]{cc}c\psi_{e}+s&\mbox{ in
}\mathbb{R}^{2}\setminus\overline{\Omega}\\\ c&\mbox{ in
}\Omega\setminus\overline{D}\\\ c\psi_{d}&\mbox{ in }D,\end{array}\right.$

for an unknown constant $c,$ where the functions $s,\psi_{e},$ and $\psi_{d}$
are as defined via (2.4), (2.5) and (2.6). Observing that continuity of
$v_{0}$ across $\partial\Omega$ and $\partial D$ is enforced by construction,
it remains to determine the unknown constant $c.$ The value of this constant
is determined at the next order. Indeed, plugging in (2.23) into (2.1), and
matching the terms at $O(\delta)$ and enforcing continuity of
$\varepsilon_{\delta}^{-1}\frac{\partial u_{\delta}}{\partial\nu_{\Omega}}$
(respectively $\varepsilon_{\delta}^{-1}\frac{\partial
u_{\delta}}{\partial\nu_{D}})$ at the interface $\partial\Omega$ (respectively
$\partial D$), one arrives at the system (2.16) for
$\phi_{0}:=v_{1}|_{\Omega\setminus\overline{D}}$ in the ENZ region, and the
systems (2.17) and (2.18) respectively for the restrictions $\lambda_{0}$ and
$\chi_{0}$ of the corrector $v_{1}$ in the regions
$\mathbb{R}^{2}\setminus\overline{\Omega}$ and $D$ respectively. We note that
the PDE (2.16) is a Poisson equation, and the value $c=c^{*}$ is exactly the
unique choice of the constant $c$ that makes this problem consistent. The
choice of the boundary conditions for $\phi_{0}$ ensure continuity of
$\varepsilon_{\delta}^{-1}\frac{\partial v_{0}}{\partial\nu}$ across both the
interfaces, and those for $\lambda_{0}$ and $\chi_{0}$ yield continuity of
$v_{1}$ across the interfaces. The subsequent order correctors are then
defined inductively, where at each stage, the constant $c^{*}$ is corrected by
suitable constants $\delta^{j}e_{j},$ that render the Poisson problems in the
ENZ region consistent, providing the Dirichlet data for the exterior and
dopant Helmholtz problems.

The bulk of the proof lies in setting up the induction procedure above, and
arguing that the system for the infinite system of correctors does indeed
close, with bounds sufficient to ensure analyticity. The key point, as one
might naturally expect in the context of divergence form PDEs with piecewise
constant coefficients, is that the regions $\Omega\setminus\overline{D}$, $D$
and $\mathbb{R}^{2}\setminus\overline{\Omega}$ interact with each other
through their Dirichlet-to-Neumann maps.

##  3\. An Operator Theoretic Set-up

We introduce various operators that will be used in the proof of our main
result.

###  3.1. Dirichlet-to-Neumann operators outside the ENZ region

Let $\mathcal{D}^{e}_{k}$ and $\mathcal{D}^{d}_{k}$ denote the Dirichlet-to-
Neumann operators associated with the Helmholtz operator $-\Delta-k^{2}$ in
the exterior and in the dopant regions respectively. To be precise, given
$g_{e}\in H^{1/2}(\partial\Omega)$ (resp. $g_{d}\in H^{1/2}(D)$), define these
operators via

$\left\\{\begin{aligned} -\Delta u_{e}&=k^{2}u_{e}\quad\mbox{ in
}\mathbb{R}^{2}\setminus\overline{\Omega},\\\ u_{e}&=g_{e}\quad\mbox{ on
}\partial\Omega,\\\ u_{e}&\mbox{ satisfies the radiation condition
(\ref{e.bcinfty}) at infinity }\\\ \mathcal{D}^{e}_{k}(g_{e})&:=\frac{\partial
u_{e}}{\partial\nu_{\Omega}}\in H^{-1/2}(\partial\Omega),\end{aligned}\right.$

and

$\left\\{\begin{aligned} -\Delta u_{d}&=k^{2}u_{d}\quad\mbox{ in }D,\\\
u_{d}&=g_{d}\quad\mbox{ on }\partial D,\\\
\mathcal{D}^{d}_{k}(g_{d})&:=\frac{\partial u_{d}}{\partial\nu_{D}}\in
H^{-1/2}(\partial D),\end{aligned}\right.$

respectively. A basic fact that we will use from [6, 7, 17] is

######  Lemma 3.1.

Under the assumptions in Section 2.2, the operators $\mathcal{D}^{e}_{k}$ and
$\mathcal{D}^{d}_{k}$ define bounded linear operators: to be precise, there
exist constants $C_{k}^{e},C_{k}^{d}$ such that

$\displaystyle\|\mathcal{D}^{e}_{k}(g_{e})\|_{H^{-1/2}(\partial\Omega)}\leqslant
C_{k}^{e}\|g_{e}\|_{H^{1/2}(\partial\Omega)},$
$\displaystyle\|\mathcal{D}^{d}_{k}(g_{d})\|_{H^{-1/2}(\partial D)}\leqslant
C_{k}^{d}\|g_{d}\|_{H^{1/2}(\partial D)},$

With this definition at hand, by a convenient abuse of notation, we will write

$\displaystyle\int_{\partial\Omega}\mathcal{D}^{e}_{k}(g)\,d\mathcal{H}^{1}:=\langle\mathcal{D}^{e}_{k}(g),\mathbf{1}\rangle_{H^{-1/2}(\partial\Omega),H^{1/2}(\partial\Omega)},$

where $\mathbf{1}(x)\equiv 1$ for all $x\in\partial\Omega$ and, we use a
similar convention at $\partial D.$ By duality,

$\displaystyle\left|\int_{\partial\Omega}\mathcal{D}^{e}_{k}(g_{e})\,d\mathcal{H}^{1}\right|\leqslant
C_{k}^{e}\|g_{e}\|_{H^{1/2}(\partial\Omega)}\|\mathbf{1}\|_{H^{1/2}(\partial\Omega)},$
(3.1) $\displaystyle\left|\int_{\partial
D}\mathcal{D}^{d}_{k}(g_{d})\,d\mathcal{H}^{1}\right|\leqslant
C_{k}^{d}\|g_{d}\|_{H^{1/2}(\partial D)}\|\mathbf{1}\|_{H^{1/2}(\partial D)}.$

A tensorized version of $\mathcal{D}_{k}^{e},\mathcal{D}_{k}^{d}$ is easily
defined via: given $(g_{e},g_{d})\in H^{1/2}(\partial\Omega)\times
H^{1/2}(\partial D),$ set

$\displaystyle\mathcal{D}_{k}(g_{e},g_{d}):=(\mathcal{D}_{k}^{e}(g_{e}),\mathcal{D}_{k}^{d}(g_{d})).$

###  3.2. An averaging operator

For any $(h_{e},h_{d})\in H^{-1/2}(\partial\Omega)\times H^{-1/2}(\partial
D),$ we define

$\displaystyle\mathcal{A}(h_{e},h_{d}):=-\frac{1}{\beta}\left(\int_{\partial\Omega}h_{e}\,d\mathcal{H}^{1}-\int_{\partial
D}h_{d}\,d\mathcal{H}^{1}\right)\in\mathbb{C},$

recalling that $\beta\neq 0$ by Lemma 2.1.

######  Lemma 3.2.

The functional $\mathcal{A}:H^{1/2}(\partial\Omega)\times H^{1/2}(\partial
D)\to\mathbb{C}$ is a linear map that is bounded in the sense that there
exists a constant $C_{\mathcal{A}}>0$ with

$\displaystyle|\mathcal{A}(h_{e},h_{d})|\leqslant
C_{\mathcal{A}}(\|h_{e}\|_{H^{-1/2}(\partial\Omega)}+\|h_{d}\|_{H^{-1/2}(\partial
D)}).$ (3.2)

###### Proof.

The desired estimate follows easily from (3.1) and the triangle inequality. ∎

###  3.3. The solution operators

We introduce solution operators in the various regions. In order to define
this operator, we let $\underline{L}^{2}(\Omega\setminus\overline{D})$ denote
the set of square integrable functions in $\Omega\setminus\overline{D},$ whose
average is zero: precisely,

$\underline{L}^{2}(\Omega\setminus\overline{D}):=\left\\{f\in
L^{2}(\Omega\setminus\overline{D}):\int_{\Omega\setminus\overline{D}}f=0\right\\},$

and

$\displaystyle\underline{H}^{1}(\Omega\setminus\overline{D}):=\left\\{f\in
H^{1}(\Omega\setminus\overline{D}):\int_{\Omega\setminus\overline{D}}f=0\right\\}\,.$

We define the solution operator in the ENZ region, denoted by
$\mathcal{P}_{k}:\underline{L}^{2}(\Omega\setminus\overline{D})\times
H^{-1/2}(\partial\Omega)\times H^{-1/2}(\partial
D)\to\underline{H}^{1}(\Omega\setminus\overline{D})$ by

$\displaystyle\mathcal{P}_{k}(g,h_{e},h_{d}):=\hat{\phi},$

where $\hat{\phi}=\hat{\phi}(g,h_{e},h_{d})$ is the unique $H^{1}$ weak
solution to the boundary value problem

$\left\\{\begin{aligned}
-\Delta\hat{\phi}&=k^{2}\mathcal{A}(h_{e},h_{d})+k^{2}g,\quad\quad\mbox{ in
}\Omega\setminus\overline{D},\\\
\frac{\partial\hat{\phi}}{\partial\nu_{\Omega}}&=\mathcal{A}(h_{e},h_{d})\frac{\partial\psi_{e}}{\partial\nu}+h_{e}\quad\quad\mbox{
on }\partial\Omega,\\\
\frac{\partial\hat{\phi}}{\partial\nu_{D}}&=\mathcal{A}(h_{e},h_{d})\frac{\partial\psi_{d}}{\partial\nu_{D}}+h_{d},\quad\quad\mbox{
on }\partial D,\\\
\strokedint_{\Omega\setminus\overline{D}}&\hat{\phi}=0.\end{aligned}\right.$
(3.3)

The choice of the constant $\mathcal{A}(h_{e},h_{d})$ is exactly such that
this system is consistent.

In order to define correctors in the exterior and in the dopant region, we
define the Helmholtz solution operators
$\mathcal{H}_{k}^{e}:H^{1/2}(\partial\Omega)\to
H^{1}_{loc}(\mathbb{R}^{2}\setminus\overline{\Omega})$, and respectively,
$\mathcal{H}_{k}^{d}:H^{1/2}(\partial D)\to H^{1}(D),$ via

$\left\\{\begin{aligned} -\Delta v&=k^{2}v\quad\quad\mbox{ in
}\mathbb{R}^{2}\setminus\overline{\Omega}\\\ v&=g_{e}\in
H^{1/2}(\partial\Omega),\\\ v&\mbox{ satisfies the radiation condition
(\ref{e.bcinfty}) at infinity}\\\
\mathcal{H}_{k}^{e}&(g_{e}):=v,\end{aligned}\right.$ (3.4)

and respectively,

$\left\\{\begin{aligned} -\Delta v&=k^{2}v\quad\quad\mbox{ in }D\\\
v&=g_{d}\in H^{1/2}(\partial D),\\\
\mathcal{H}_{k}^{d}&(g_{d}):=v,\end{aligned}\right.$ (3.5)

That $\mathcal{H}^{d}_{k}$ is well-defined is guaranteed by our assumption
(2.3) that the dopant isn’t resonant. As usual, we tensorize this to define
$\mathcal{H}_{k}:H^{1/2}(\partial\Omega)\times H^{1/2}(\partial D)\to
H^{1}_{loc}(\mathbb{R}^{2}\setminus\overline{\Omega})\times H^{1}(D)$ via

$\displaystyle\mathcal{H}_{k}(g_{e},g_{d})=(\mathcal{H}_{k}^{e}(g_{e}),\mathcal{H}_{k}^{d}(g_{d})).$
(3.6)

######  Lemma 3.3.

The operator
$\mathcal{P}_{k}:\underline{L}^{2}(\Omega\setminus\overline{D})\times
H^{-1/2}(\partial\Omega)\times H^{-1/2}(\partial D)\to
H^{1}(\Omega\setminus\overline{D}),$ and
$\mathcal{H}_{k}:H^{1/2}(\partial\Omega)\times H^{1/2}(D)\to
H^{1}_{loc}(\mathbb{R}^{2}\setminus\overline{\Omega})\times H^{1}(D)$ are both
bounded linear operators between the respective Hilbert spaces.

###### Proof.

See [17, Theorem 8.19] for the assertion regarding $\mathcal{P}_{k}$, and [17,
Theorem 9.11 and page 286] for the statement regarding $\mathcal{H}_{k}.$ ∎

Denoting the boundary trace operator for the ENZ region,
$\Omega\setminus\overline{D}$ by $\mathcal{T},$ it then follows that the
composition
$\mathcal{T}\mathcal{P}_{k}:\underline{L}^{2}(\Omega\setminus\overline{D})\times
H^{-1/2}(\partial\Omega)\times H^{-1/2}(\partial D)\to
H^{1/2}(\partial\Omega)\times H^{1/2}(\partial D)$ is a bounded linear
operator. The main object of our proof is the iteration map
$\mathcal{I}_{k}:\underline{H}^{1}(\Omega\setminus\overline{D})\times
H^{-1/2}(\partial\Omega)\times H^{-1/2}(\partial
D)\to\underline{H}^{1}(\Omega\setminus\overline{D})\times
H^{-1/2}(\partial\Omega)\times H^{-1/2}(\partial D)$ defined via

$\displaystyle(g,h_{e},h_{d})\in\underline{H}^{1}(\Omega\setminus\overline{D})\times
H^{-1/2}(\partial\Omega)\times H^{-1/2}(\partial D)$ (3.7)
$\displaystyle\quad\quad\quad\quad\quad\mapsto\mathcal{I}_{k}(g,h_{e},h_{d}):=(\mathcal{P}_{k}(g,h_{e},h_{d}),\mathcal{D}_{k}\mathcal{T}\mathcal{P}_{k}(g,h_{e},h_{d})).$

In particular, $\mathcal{I}_{k}$ is a bounded linear operator.

Recall the number $c^{*}$ defined in (2.13), and the leading order correctors
$(\phi_{0},\lambda_{0},\chi_{0})$ defined in (2.16), (2.17), and (2.18)
respectively. We note that

$\displaystyle c^{*}=\mathcal{A}\left(\frac{\partial
s}{\partial\nu_{\Omega}},0\right),$ (3.8)

$\phi_{0}:=\mathcal{P}_{k}\left(0,\frac{\partial
s}{\partial\nu_{\Omega}},0\right),$ (3.9)

and that

$\displaystyle(\lambda_{0},\chi_{0})=\mathcal{H}_{k}\mathcal{T}\mathcal{P}_{k}\left(0,\frac{\partial
s}{\partial\nu_{\Omega}},0\right).$ (3.10)

Finally, we note

$\displaystyle\left(\frac{\partial\lambda_{0}}{\partial\nu_{\Omega}},\frac{\partial\chi_{0}}{\partial\nu_{D}}\right)=\mathcal{D}_{k}\mathcal{T}\mathcal{P}_{k}\left(0,\frac{\partial
s}{\partial\nu_{\Omega}},0\right),$ (3.11)

so that the number $e_{0}$ defined in (2.19) satisfies

$\displaystyle
e_{0}=\mathcal{A}\left(\mathcal{D}_{k}\mathcal{T}\mathcal{P}_{k}\left(0,\frac{\partial
s}{\partial\nu_{\Omega}},0\right)\right).$ (3.12)

##  4\. Proofs of the Main Theorem and its Corollaries

###### Proof of Theorem 2.2.

We proceed in several steps, deriving recursive definitions of higher order
correctors to all orders. Before setting up the induction step, it is
convenient to define

$\displaystyle\phi_{-1}\equiv 0,\,\lambda_{-1}:=s,\,\chi_{-1}\equiv 0,\,\mbox{
and }e_{-1}:=c^{*}=\mathcal{A}\left(\frac{\partial
s}{\partial\nu_{\Omega}},0\right),$ (4.1)

and recall the definitions for $(\phi_{0},\lambda_{0},\chi_{0},e_{0})$ from
(3.9)–(3.12).

1. In this step, we define a hierarchy of correctors. For each $j=1,2,\ldots,$ we inductively define,

$\displaystyle\phi_{j}$
$\displaystyle:=\mathcal{P}_{k}\left(\phi_{j-1},\frac{\partial\lambda_{j-1}}{\partial\nu_{\Omega}},\frac{\partial\chi_{j-1}}{\partial\nu_{D}}\right),$
(4.2) $\displaystyle(\lambda_{j},\chi_{j})$
$\displaystyle=\mathcal{H}_{k}\mathcal{T}\mathcal{P}_{k}\left(\phi_{j-1},\frac{\partial\lambda_{j-1}}{\partial\nu_{\Omega}},\frac{\partial\chi_{j-1}}{\partial\nu_{D}}\right),$

and so, observing that

$\left(\frac{\partial\lambda_{j}}{\partial\nu_{\Omega}},\frac{\partial\chi_{j}}{\partial\nu_{D}}\right)=\mathcal{D}_{k}\mathcal{T}\mathcal{P}_{k}\left(\phi_{j-1},\frac{\partial\lambda_{j-1}}{\partial\nu_{\Omega}},\frac{\partial\chi_{j-1}}{\partial\nu_{D}}\right),$
(4.3)

we set

$\displaystyle
e_{j}:=-\frac{1}{\beta}\left[\int_{\partial\Omega}\frac{\partial\lambda_{j}}{\partial\nu_{\Omega}}\,d\mathcal{H}^{1}-\int_{\partial
D}\frac{\partial\chi_{j}}{\partial\nu_{D}}\,d\mathcal{H}^{1}\right]=\mathcal{A}\left(\frac{\partial\lambda_{j}}{\partial\nu_{\Omega}},\frac{\partial\chi_{j}}{\partial\nu_{D}}\right),\quad\quad
j\in\mathbb{N},$ (4.4)

where we recall that, by (2.7), that $\beta\neq 0.$ Therefore, by definition
of $\mathcal{I}_{k}$, and the base case (4.1), it follows that

$\displaystyle\left(\phi_{j},\frac{\partial\lambda_{j}}{\partial\nu_{\Omega}},\frac{\partial\chi_{j}}{\partial\nu_{D}}\right)=\mathcal{I}_{k}^{j+1}\left(0,\frac{\partial
s}{\partial\nu_{\Omega}},0\right),\quad j\in\mathbb{N}.$ (4.5)

By induction, (3.11), and (4.3), and (4.5),

$\displaystyle
e_{j}=\mathcal{A}\left(\frac{\partial\lambda_{j}}{\partial\nu_{\Omega}},\frac{\partial\chi_{j}}{\partial\nu_{D}}\right)$
$\displaystyle\quad=\mathcal{A}\left(\mathcal{D}_{k}\mathcal{T}\mathcal{P}_{k}\left(\phi_{j-1},\frac{\partial\lambda_{j-1}}{\partial\nu_{\Omega}},\frac{\partial\chi_{j-1}}{\partial\nu_{D}}\right)\right)=\mathcal{A}\left(\mathcal{D}_{k}\mathcal{T}\mathcal{P}_{k}\mathcal{I}_{k}^{j}\left(0,\frac{\partial
s}{\partial\nu_{\Omega}},0\right)\right).$

Similarly, combining (4.2) and (4.5),

$\displaystyle\phi_{j}$
$\displaystyle=\mathcal{P}_{k}\mathcal{I}_{k}^{j}\left(0,\frac{\partial
s}{\partial\nu_{\Omega}},0\right),$ (4.6)
$\displaystyle(\lambda_{j},\chi_{j})$
$\displaystyle=\mathcal{H}_{k}\mathcal{T}\mathcal{P}_{k}\mathcal{I}_{k}^{j}\left(0,\frac{\partial
s}{\partial\nu_{\Omega}},0\right).$

As always, in the foregoing equations, an exponent on an operator means
repeated composition.

2. In this step, we obtain bounds on the correctors constructed in Step 1. By Lemmas 3.1 and (3.3), $\mathcal{I}_{k}$ defined in (3.7) is a bounded linear operator; we let $\Lambda$ denote its operator norm. By Lemmas 3.2 and (2.11) we find 111Here, for two sequences $\\{a_{j}\\}_{j},\\{b_{j}\\}_{j}\subset\mathbb{R},$ we use the notation $a_{j}\lesssim b_{j}$ to mean $a_{j}\leqslant Cb_{j},j\in\mathbb{N},$ for a universal constant $C$ that is independent of $j.$

$\displaystyle|e_{j}|\lesssim\Lambda^{j}\left\|\frac{\partial
s}{\partial\nu_{\Omega}}\right\|_{H^{-1/2}(\partial\Omega)}\lesssim\Lambda^{j}\|f\|_{L^{2}}.$
(4.7)

Similarly, from (4.6), and Lemmas 3.1, 3.3, and (2.11), we estimate

$\displaystyle\|\phi_{j}\|_{H^{1}(\Omega\setminus\overline{D})}\lesssim\Lambda^{j}\|f\|_{L^{2}},\quad\|\chi_{j}\|_{H^{1}(D)}\lesssim\Lambda^{j}\|f\|_{L^{2}},\quad\mbox{and}$
(4.8)
$\displaystyle\|\lambda_{j}\|_{H^{1}(U)}\lesssim\Lambda^{j}\|f\|_{L^{2}}$

for all open sets $U\subset\mathbb{R}^{2}\setminus\overline{\Omega},$ with
compact closure.

Let $\delta\in\mathbb{C}$ be such that $|\delta|<\frac{1}{\Lambda}.$ For such
$\delta,$ the operator $(I-\delta\mathcal{I}_{k})$ is invertible, and its
inverse defines a bounded linear operator on the space
$\underline{H}^{1}(\Omega\setminus\overline{D})\times
H^{-1/2}(\partial\Omega)\times H^{-1/2}(\partial D).$ It is a standard fact
that the associated Neumann series converges absolutely on compact subsets of
$\\{|\delta|\leqslant\frac{1}{\Lambda}\\}$ . Using the linearity of the
operators $\mathcal{A},\mathcal{D}_{k},\mathcal{T},\mathcal{P}_{k}$ and the
definition of the Neumann series, and (3.8), we find

$\displaystyle\hat{c}_{\delta}:=c^{*}+\sum_{j=0}^{\infty}\delta^{j+1}e_{j}=c^{*}+\sum_{j=0}^{\infty}\delta^{j+1}\mathcal{A}\left(\mathcal{D}_{k}\mathcal{T}\mathcal{P}_{k}\mathcal{I}_{k}^{j}\left(0,\frac{\partial
s}{\partial\nu_{\Omega}},0\right)\right)$ (4.9)
$\displaystyle\quad\quad=c^{*}+\delta\mathcal{A}\left(\mathcal{D}_{k}\mathcal{T}\mathcal{P}_{k}\sum_{j=0}^{\infty}\delta^{j}\mathcal{I}_{k}^{j}\left(0,\frac{\partial
s}{\partial\nu_{\Omega}},0\right)\right)$
$\displaystyle\quad\quad=c^{*}+\delta\mathcal{A}\mathcal{D}_{k}\mathcal{T}\mathcal{P}_{k}(I-\delta\mathcal{I}_{k})^{-1}\left(0,\frac{\partial
s}{\partial\nu_{\Omega}},0\right),$

satisfies $|\hat{c}_{\delta}|<\infty$ for all $|\delta|<\frac{1}{\Lambda}$.
Next, we define, and compute using (3.9) and (4.6) that

$\displaystyle\hat{\phi}_{\delta}:=\sum_{j=0}^{\infty}\delta^{j}\phi_{j}=\sum_{j=0}^{\infty}\delta^{j}\mathcal{P}_{k}\mathcal{I}_{k}^{j}\left(0,\frac{\partial
s}{\partial\nu_{\Omega}},0\right)$ (4.10)
$\displaystyle\quad=\mathcal{P}_{k}\left((I-\delta\mathcal{I}_{k})^{-1}\left(0,\frac{\partial
s}{\partial\nu_{\Omega}},0\right)\right)$

where, the convergence is understood in the sense of the Hilbert space
$H^{1}(\Omega\setminus\overline{D})$. Finally, we define, and compute using
(3.10) and (4.6) that

$\displaystyle(\hat{\lambda}_{\delta},\hat{\chi}_{\delta}):=\sum_{j=0}^{\infty}\delta^{j}(\lambda_{j},\chi_{j})$
(4.11)
$\displaystyle\quad\quad=\sum_{j=0}^{\infty}\delta^{j}\mathcal{H}_{k}\mathcal{T}\mathcal{P}_{k}\mathcal{I}_{k}^{j}\left(0,\frac{\partial
s}{\partial\nu_{\Omega}},0\right)$
$\displaystyle\quad\quad=\mathcal{H}_{k}\mathcal{T}\mathcal{P}_{k}\left((I-\delta\mathcal{I}_{k})^{-1}\left(0,\frac{\partial
s}{\partial\nu_{\Omega}},0\right)\right),$

where, once again, the convergence of the infinite sum is understood in the
sense of the Hilbert space
$H^{1}_{loc}(\mathbb{R}^{2}\setminus\overline{D})\times H^{1}(D).$

3. In this step we define a candidate for the solution $u_{\delta}$ of (2.1). Indeed, define

$\hat{u}_{\delta}(x):=\left\\{\begin{array}[]{cc}\hat{c}_{\delta}\psi_{e}+s+\delta\hat{\lambda}_{\delta}&x\in\mathbb{R}^{2}\setminus\overline{\Omega}\\\
\hat{c}_{\delta}+\delta\hat{\phi}_{\delta}&x\in\Omega\setminus\overline{D},\\\
\hat{c}_{\delta}\psi_{d}+\delta\hat{\chi}_{\delta}&x\in D.\end{array}\right.$
(4.12)

By construction, $\hat{u}_{\delta}$ and
$\frac{1}{\varepsilon_{\delta}}\frac{\partial u_{\delta}}{\partial\nu}$ are
continuous along $\partial\Omega$ and $\partial D,$ and in particular,
$\hat{u}_{\delta}\in H^{1}_{loc}(\mathbb{R}^{2})$ by Lemma 3.3. Using partial
sums, it is easily seen that by construction, $\hat{u}_{\delta}$ is a weak
solution to (2.1). As, by assumption, (2.1) is non-resonant,
$\hat{u}_{\delta}$ coincides with its unique solution.

4. We are ready to complete the proof of the theorem. Defining $v_{j}:=\lambda_{j}$ in the exterior, $:=e_{j}+\phi_{j}$ in the ENZ region, and $:=\chi_{j}$ in the dopant, the conclusion of the theorem follows using the estimates (4.7) and (LABEL:e.jcorrest) from Step 2. The bound on $v_{j}$ follows by combining the estimates in Step 2 with (2.11). ∎

Having completed the proof of our main result, we turn to the proofs of the
corollaries stated in Section 2.4.

###### Proof of Corollary 2.4.

From the proof of Theorem 2.2, we know that $u_{\delta}=\hat{u}_{\delta},$
defined in (4.12). It follows, therefore, that

$u_{\delta}(x)-v_{\delta}(x)=\left\\{\begin{array}[]{cc}\left(\sum_{j=2}^{\infty}\delta^{j}e_{j}\right)\psi_{e}+\sum_{j=2}^{\infty}\delta^{j}\lambda_{j},&x\in\mathbb{R}^{2}\setminus\overline{\Omega}\\\
\sum_{j=2}^{\infty}\delta^{j}(e_{j}+\phi_{j})&x\in\Omega\setminus\overline{D}\\\
\left(\sum_{j=2}^{\infty}\delta^{j}e_{j}\right)\psi_{d}+\sum_{j=2}^{\infty}\delta^{j}\chi_{j},&x\in
D.\end{array}\right.$

In the above, each of the summations is interpreted as a convergent series in
$H^{1}$ of the appropriate region (and $H^{1}_{loc}$ in the exterior). The
proof of the corollary is then immediate from (4.6),
$|\delta|<\frac{1}{\Lambda},$ and estimates (4.7) and (LABEL:e.jcorrest). ∎

###### Proof of Corollary 2.6.

Since the Poynting vector is defined by
$\vec{S}_{\delta}=\frac{1}{2}\vec{E}_{\delta}\times\overline{\vec{H}_{\delta}},$
the proof of this corollary is a short computation that combines
$\vec{H}_{\delta}=(0,0,u_{\delta}),$ and
$\vec{E}_{\delta}=\frac{1}{i\omega}(-\partial_{x_{2}}\phi_{0},\partial_{x_{1}}\phi_{0},0)+O(\delta)$
(see Remark 2.5), along with (2.21) to pass to the limit $\delta\to 0$ limit.
The boundary conditions for $S$ come from the Neumann boundary conditions for
$\phi_{0}$ in (2.16). ∎

##  5\. What if the Dopant is Resonant?

For any choice of the dopant geometry, there will be some choices of $k$ for
which the dopant is resonant, and it is natural to ask what happens then. The
answer one finds in [15] is that:

if the dopant resonance is excited then the magnetic field (5.1)
$\displaystyle\mbox{vanishes in the ENZ region and }\mu_{\rm eff}=\infty\,.$

This section starts by briefly discussing what goes wrong when the dopant is
resonant; then we shall recover (5.1) using our framework, and we’ll go beyond
it by identifying the leading-order electric field.

The dopant is resonant when $k^{2}$ is an eigenvalue of the Laplace operator
on $D$ with zero Dirichlet boundary conditions, i.e., there exist finitely
many (say $m\in\mathbb{N}$) eigenfunctions $\\{U_{j}\\}_{j=1}^{m}\subset
H^{1}_{0}(D)$ that are orthonormal in $L^{2}(D),$ and satisfy, for every
$j=1,\ldots,m,$

$\displaystyle-\Delta U_{j}$ $\displaystyle=$ $\displaystyle
k^{2}U_{j}\quad\mbox{in $D$}$ $\displaystyle U_{j}$ $\displaystyle=$
$\displaystyle 0\quad\quad\mbox{at $\partial D$}.$

The finiteness of $m$ follows from standard theory. For such $k$, our
framework breaks down in one of two mutually exclusive ways:

  1. (1)

Suppose first that the eigenfunctions $U_{j},j=1,\ldots,m$ each satisfy
$\int_{D}U_{j}=0,$ i.e., they each have mean zero in $D.$ In this case the
auxiliary function $\psi_{d}$ exists but it isn’t unique. Indeed, setting
$\Psi_{d}:=\psi_{d}-1,$ the existence of $\psi_{d}$ is equivalent to the
solvability of the PDE

$\displaystyle-\Delta\Psi_{d}=k^{2}(\Psi_{d}+1)\quad\quad\mbox{ in }D,$

for $\Psi_{d}\in H^{1}_{0}(D).$ By the Fredholm alternative, $\Psi_{d}$
exists, as the constant function $1$ is orthogonal in $L^{2}(D)$ to $U_{j},$
i.e., $\int_{D}U_{j}=0$ for each $j=1,\ldots,m.$ However it is nonunique,
since for any $\\{\alpha_{j}\\}_{j=1}^{m}\subset\mathbb{C},$ the choices

$\displaystyle\psi_{d}=\Psi_{d}+1=1+k^{2}\left[\sum_{\lambda\in\sigma(-\Delta)\setminus\\{k^{2}\\}}\frac{\int_{D}\phi_{\lambda}\,dx}{\lambda-k^{2}}\phi_{\lambda}\right]+\sum_{j=1}^{m}\alpha_{j}U_{j}$
(5.2)

solve (2.6). (Here $\sigma(-\Delta)$ denotes the set of Laplacian eigenvalues
with zero Dirichlet boundary data, and $\phi_{\lambda}$ is the eigenfunction
associated with eigenvalue $\lambda,$ normalized so that
$\|\phi_{\lambda}\|_{L^{2}(D)}=1.$)

Heuristically, we can say in this case that “the dopant resonance isn’t
excited,” since there is apparently nothing preventing us from choosing
$\alpha_{j}=0$ for $j=1,\ldots,m$. Moreover, since $\int U_{j}=0,$ the value
of $c^{*}$ (which feels $\psi_{d}$ only via $\int_{D}\psi_{d}$) is oblivious
to the choice of $\\{\alpha_{j}\\}$. However, the correction $\chi_{0}$ (see
(2.18)) depends on the choice of $\\{\alpha_{j}\\}$.

Let us record that, from (5.2),

$\displaystyle\int_{D}\psi_{d}=|D|+k^{2}\left[\sum_{\lambda\in\sigma(-\Delta)\setminus\\{k^{2}\\}}\frac{\bigl{(}\int_{D}\phi_{\lambda}\,dx\bigr{)}^{2}}{\lambda-k^{2}}\right].$
(5.3)

(This sum is convergent: by Weyl’s law, the $j$th eigenvalue $\lambda_{j}\sim
j^{2/d}=j,$ so the factors $\frac{1}{\lambda-k^{2}}$ share the decay of only
the harmonic series [which logarithmically blows up!]. However, the
convergence of the series defining $\int\psi_{d}$ is assisted by the fact that
the sequence of $L^{2}$-normalized eigenfunctions
$\\{\phi_{\lambda}\\}_{\lambda}$ converge weakly in $L^{2}$ to zero, so
$\int_{D}\phi_{\lambda}\,dx\to 0.$)

  2. (2)

Suppose instead that one of the eigenfunctions $U_{j}$ is such that
$\int_{D}U_{j}\neq 0.$ (This happens, for instance, if $k^{2}$ is the smallest
Dirichlet eigenvalue of $-\Delta$ in $D$ since the eigenspace is then simple
and the eigenfunction – taken to be real-valued – doesn’t change sign.) In
this case, rather than nonuniqueness the difficulty is nonexistence: the
boundary value problem (2.6) describing $\psi_{d}$ doesn’t have a solution.

For reasons that will become clear presently, we view this as the case when
“the dopant resonance is excited.”

To achieve a better understanding, it is natural to consider a limit in which
the resonant situation is approached via nonresonant ones. There are (at
least) two natural approaches:

  1. (i)

We can ask what happens as $k^{2}$ approaches an eigenvalue along the real
axis; this means choosing an eigenvalue $\lambda_{*}$ of $-\Delta$ in $D$ with
the Dirichlet boundary condition, and considering $k$ such that
$\lambda_{*}-k^{2}=\gamma,$ where $\gamma$ is real with $0<|\gamma|\ll 1.$

  2. (ii)

We can consider what happens when resonance is prevented by material losses.
This can be achieved, for example, by taking the magnetic permeability $\mu$
to have a positive imaginary part. This amounts in our framework to taking
$k^{2}=\lambda_{*}+i\gamma$, where $\lambda_{*}$ is an eigenvalue (as above)
and $0<\gamma\ll 1$.

Either way, our analysis applies to the limit $\delta\rightarrow 0$ when
$\gamma$ is held fixed. We are therefore able to consider the limiting
behavior when $\delta\to 0$ first, then $\gamma\to 0.$

Let us pursue this idea in case (2) (when $\psi_{d}$ doesn’t exist). Our
analysis encompasses both alternatives (i) and (ii) suggested above, and is in
fact more general. Fixing $\lambda_{*}\in\sigma(-\Delta),$ we consider any
$k\in\mathbb{C}$ with $0\leqslant\arg\,k<\pi$ such that $k^{2}$ is near
$\lambda_{*}$, say

$\lambda_{*}-k^{2}=\gamma\quad\mbox{ with }|\gamma|\mbox{ small but nonzero}.$

We shall identify the limiting behavior as $|\gamma|\to 0$ with

$0<|\delta|\ll|\gamma|\ll 1.$

To begin, we observe (by an obvious adjustment to (5.3)) that

$\displaystyle\int_{D}\psi_{d}=|D|+k^{2}\left[\sum_{\lambda\in\sigma(-\Delta)}\frac{\left(\int_{D}\phi_{\lambda}\,dx\right)^{2}}{\lambda-k^{2}}\right].$
(5.4)

We are assuming that at least one of the eigenfunctions associated with
$\lambda=\lambda_{*}$ has nonzero mean. It follows that the value of (5.4)
blows up as $|\gamma|\to 0$, though it is finite (and large, of order
$1/|\gamma|$) for fixed $\gamma\neq 0$. Therefore

$|\beta|\sim\left|\int_{\partial
D}\frac{\partial\psi_{d}}{\partial\nu_{D}}\,d\mathcal{H}^{1}\right|\sim\frac{1}{|\gamma|},$

and it follows that $c^{*}$ is small, and $|\mu_{\rm eff}|$ is large:

$|c^{*}|=O(|\gamma|)\quad\mbox{and}\quad\mu_{\rm eff}=O(|\gamma|^{-1})$ (5.5)

For the rest of this section it is convenient to track the dependence of
various quantities on $\gamma,$ and we do this by a subscript: for instance,
we specify the dependence of $c^{*}$ on $\gamma$ through the notation
$c^{*}_{\gamma}.$

Theorem 2.2 applies, and asserts that the leading order corrections in the ENZ
region, in the $\delta\to 0$ limit (for fixed $\gamma$) satisfy (2.16):

$\left\\{\begin{aligned}
&-\Delta\phi^{\gamma}_{0}=k_{\gamma}^{2}c_{\gamma}^{*}\quad\quad\mbox{ in
}\Omega\setminus\overline{D}\\\
&\frac{\partial\phi_{0}^{\gamma}}{\partial\nu_{\Omega}}=c_{\gamma}^{*}\frac{\partial\psi_{e}}{\partial\nu_{\Omega}}+\frac{\partial
s}{\partial\nu_{\Omega}}\quad\quad\mbox{ on }\partial\Omega\\\
&\frac{\partial\phi^{\gamma}_{0}}{\partial\nu_{D}}=c_{\gamma}^{*}\frac{\partial\psi^{\gamma}_{d}}{\partial\nu_{D}}\quad\quad\mbox{
on }\partial D,\\\
&\strokedint_{\Omega\setminus\overline{D}}\phi^{\gamma}_{0}\,dx=0.\end{aligned}\right.$
(5.6)

It is convenient to pass to the $\gamma\to 0^{+}$ limit by writing (5.6) in
weak form. Toward this end, testing (2.16) with $h\in
C^{\infty}(\Omega\setminus\overline{D}),$ we obtain

$\int_{\Omega\setminus\overline{D}}\nabla\phi_{0}^{\gamma}\cdot\nabla
h\,dx-\int_{\partial\Omega}h\left(c_{\gamma}^{*}\frac{\partial\psi_{e}}{\partial\nu_{\Omega}}+\frac{\partial
s}{\partial\nu_{\Omega}}\right)\,d\mathcal{H}^{1}-\int_{\partial
D}h\left(c_{\gamma}^{*}\frac{\partial\psi^{\gamma}_{d}}{\partial\nu_{D}}\right)\,d\mathcal{H}^{1}=\\\
k_{\gamma}^{2}c_{\gamma}^{*}\int_{\Omega\setminus\overline{D}}h\,dx.$ (5.7)

Arguing formally, if $\phi_{0}^{\gamma}$ converges weakly to
$\widehat{\phi_{0}}$ in $H^{1}(D)$, then since $k_{\gamma}^{2}\to\lambda_{*}$
, $\psi_{e},s$ are independent of $\gamma,$ and $c^{*}_{\gamma}\to 0$ as
$\gamma\to 0,$ it is clear how to pass to the limit in the first two terms on
the left, as well as the term on the right hand side of (5.7); it remains to
pass to the limit in the term $\int_{\partial
D}hc_{\gamma}^{*}\frac{\partial\psi_{d}^{\gamma}}{\partial\nu_{D}}\,d\mathcal{H}^{1}.$
The evaluation of this limit is hindered by the fact that $c_{\gamma}^{*}\to
0$ while, formally, $\psi_{d}^{\gamma}$ blows up. This can however be easily
analyzed by integration by parts:

$\displaystyle c_{\gamma}^{*}\int_{\partial
D}h\frac{\partial\psi_{d}^{\gamma}}{\partial\nu_{D}}\,d\mathcal{H}^{1}$
$\displaystyle=c_{\gamma}^{*}\int_{D}h\Delta\psi_{d}^{\gamma}\,dx+c_{\gamma}^{*}\int_{D}\nabla
h\cdot\nabla\psi_{d}^{\gamma}\,dx$ (5.8)
$\displaystyle=-c_{\gamma}^{*}k_{\gamma}^{2}\int_{D}h\psi_{d}^{\gamma}\,dx+c_{\gamma}^{*}\int_{D}\nabla
h\cdot\nabla\psi_{d}^{\gamma}\,dx.$

Let $\\{U_{j}\\}_{j=1}^{m}$ denote the eigenfunctions of $-\Delta$ in $D$ with
homogeneous Dirichlet boundary conditions associated with an eigenvalue
$\lambda_{*}$; here, $m$ denotes the multiplicity of $\lambda_{*},$ which is
finite by Fredholm theory. Then, using (5.2) with
$\lambda_{*}-k^{2}=\gamma\not\in\sigma(-\Delta),$ one easily finds that

$\displaystyle\lim_{\gamma\to
0}\gamma\beta_{\gamma}=\lambda_{*}^{2}\sum_{j=1}^{m}\left(\int_{D}U_{j}\,dx\right)^{2},$

and so, the limit

$\displaystyle\lim_{\gamma\to
0}\frac{c^{*}_{\gamma}}{\gamma}=-\frac{\int_{\partial\Omega}\frac{\partial
s}{\partial\nu_{\Omega}}\,d\mathcal{H}^{1}}{\gamma\beta_{\gamma}}\to-\frac{\int_{\partial\Omega}\frac{\partial
s}{\partial\nu_{\Omega}}\,d\mathcal{H}^{1}}{\lambda_{*}^{2}\sum_{j=1}^{m}\left(\int_{D}U_{j}\,dx\right)^{2}}=:\overline{C}.$
(5.9)

exists and is finite. Since $k_{\gamma}^{2}\to\lambda_{*}$ as $\gamma\to 0,$
using (5.9), it then follows that the right hand side of (5.8) converges to

$\displaystyle\overline{C}\sum_{j=1}^{m}-\lambda_{*}^{2}\left(\int_{D}U_{j}\,dx\right)\left(\int_{D}hU_{j}\,dx\right)+\lambda_{*}\left(\int_{D}U_{j}\,dx\right)\left(\int_{D}\nabla
h\cdot\nabla U_{j}\,dx\right),$

which, by another integration by parts, rewrites as

$\displaystyle\overline{C}\lambda_{*}\sum_{j=1}^{m}\left(\int_{D}U_{j}\,dx\right)\left(\int_{\partial
D}h\frac{\partial U_{j}}{\partial\nu_{D}}\,d\mathcal{H}^{1}\right).$ (5.10)

Therefore, sending $\gamma\to 0$ in (5.7), and using (5.10), we find that
$\widehat{\phi_{0}}$ is the unique solution to the problem

$\left\\{\begin{aligned} &-\Delta\widehat{\phi_{0}}=0\quad\quad\mbox{ in
}\Omega\setminus\overline{D}\\\
&\frac{\partial\widehat{\phi_{0}}}{\partial\nu_{\Omega}}=\frac{\partial
s}{\partial\nu_{\Omega}}\quad\quad\mbox{ on }\partial\Omega\\\
&\frac{\partial\widehat{\phi_{0}}}{\partial\nu_{D}}=\overline{C}\lambda_{*}\sum_{j=1}^{m}\left(\int_{D}U_{j}\,dx\right)\frac{\partial
U_{j}}{\partial\nu_{D}}\quad\quad\mbox{ on }\partial D,\\\
&\strokedint_{\Omega\setminus\overline{D}}\widehat{\phi_{0}}\,dx=0.\end{aligned}\right.$
(5.11)

In particular, the electric field distribution in the ENZ region, to leading
order, is still identified by this procedure. (We note that (5.11) does have a
solution, since by (5.9) the consistency condition for Laplace’s equation with
Neumann boundary conditions is satisfied.)

The preceding discussion applies only in the double limit when $\delta\to 0$
first then $\gamma\to 0$. It would be interesting to have a more complete
understanding, for example one that applies for $\delta$ in a neighborhood of
$0$ that doesn’t depend on $\gamma$. However this seems to require a method
different from that of the present paper.

## Acknowledgements

The authors warmly thank Nader Engheta for bringing this problem to our
attention, and for helpful discussions. This work was partially supported by
grants from the National Science Foundation (DMS-2009746, RVK) and the Simons
Foundation (award # 733694, RVK and RV). We are grateful to the anonymous
referee for helpful comments on our paper.

## References

  * [1] A. Alù, M. G. Silveirinha, and N. Engheta.  Transmission-line analysis of $\varepsilon$-near-zero–filled narrow channels.  Physical Review E, 78(1):016604, 2008.
  * [2] G. Bouchitté and D. Felbacq.  Homogenization near resonances and artificial magnetism from dielectrics.  C. R. Math. Acad. Sci. Paris, 339(5):377–382, 2004.
  * [3] O. P. Bruno.  The effective conductivity of strongly heterogeneous composites.  Proc. Roy. Soc. London Ser. A, 433(1888):353–381, 1991.
  * [4] M. Camar-Eddine and P. Seppecher.  Closure of the set of diffusion functionals with respect to the Mosco-convergence.  Mathematical Models and Methods in Applied Sciences, 12(8):1153–1176, 2002.
  * [5] M. Camar-Eddine and P. Seppecher.  Determination of the closure of the set of elasticity functionals.  Archive for Rational Mechanics and Analysis, 170(3):211–245, 2003.
  * [6] D. Colton and R. Kress.  Integral equation methods in scattering theory, volume 72 of Classics in Applied Mathematics.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2013.  Reprint of the 1983 original [ MR0700400].
  * [7] D. Colton and R. Kress.  Inverse acoustic and electromagnetic scattering theory, volume 93 of Applied Mathematical Sciences.  Springer, Cham, [2019] ©2019.  Fourth edition of [ MR1183732].
  * [8] B. Edwards, A. Alù, M. G. Silveirinha, and N. Engheta.  Reflectionless sharp bends and corners in waveguides using epsilon-near-zero effects.  Journal of Applied Physics, 105(4):044905, 2009.
  * [9] B. Edwards, A. Alù, M. E. Young, M. Silveirinha, and N. Engheta.  Experimental verification of epsilon-near-zero metamaterial coupling and energy squeezing using a microwave waveguide.  Physical Review Letters, 100(3):033903, 2008.
  * [10] N. Engheta.  Pursuing near-zero response.  Science, 340(6130):286–287, 2013.
  * [11] A. Kirsch and F. Hettlich.  The mathematical theory of time-harmonic Maxwell’s equations, volume 190 of Applied Mathematical Sciences.  Springer, Cham, 2015.
  * [12] L. D. Landau, E. M. Lifshitz, and L. P. Pitaevskii.  Electrodynamics of continuous media.  Elsevier, second edition, 1984.  Course of Theoretical Physics volume 8.
  * [13] I. Liberal and N. Engheta.  The rise of near-zero-index technologies.  Science, 358(6370):1540–1541, 2017.
  * [14] I. Liberal, M. Lobet, Y. Li, and N. Engheta.  Near-zero-index media as electromagnetic ideal fluids.  Proceedings of the National Academy of Sciences, 117(39):24050–24054, 2020.
  * [15] I. Liberal, A. M. Mahmoud, Y. Li, B. Edwards, and N. Engheta.  Photonic doping of epsilon-near-zero media.  Science, 355(6329):1058–1062, 2017.
  * [16] M. Maier, M. Mattheakis, E. Kaxiras, M. Luskin, and D. Margetis.  Homogenization of plasmonic crystals: seeking the epsilon-near-zero effect.  Proceedings of the Royal Society A, 475(2230):20190220, 2019.
  * [17] W. McLean.  Strongly elliptic systems and boundary integral equations.  Cambridge University Press, Cambridge, 2000.
  * [18] X. Niu, X. Hu, S. Chu, and Q. Gong.  Epsilon-near-zero photonics: a new platform for integrated devices.  Advanced Optical Materials, 6(10):1701292, 2018.
  * [19] M. Silveirinha and N. Engheta.  Tunneling of electromagnetic energy through subwavelength channels and bends using $\varepsilon$-near-zero materials.  Physical Review Letters, 97(15):157403, 2006.
  * [20] M. Silveirinha and N. Engheta.  Design of matched zero-index metamaterials using nonmagnetic inclusions in epsilon-near-zero media.  Physical Review B, 75(7):075119, 2007.
  * [21] M. G. Silveirinha and N. Engheta.  Theory of supercoupling, squeezing wave energy, and field confinement in narrow channels and tight bends using $\varepsilon$ near-zero metamaterials.  Physical Review B, 76(24):245109, 2007.
  * [22] L. Tartar.  An introduction to Sobolev spaces and interpolation spaces, volume 3 of Lecture Notes of the Unione Matematica Italiana.  Springer, Berlin; UMI, Bologna, 2007.
  * [23] V. V. Zhikov.  On an extension and an application of the two-scale convergence method.  Mat. Sb., 191(7):31–72, 2000.
  * [24] Z. Zhou, Y. Li, E. Nahvi, H. Li, Y. He, I. Liberal, and N. Engheta.  General impedance matching via doped epsilon-near-zero media.  Physical Review Applied, 13(3):034005, 2020.