††thanks: The work of Liqun Cao and Yongwei Zhang is supported by National
Natural Science Foundation of China (grant 11971030, 11571353, 91330202) and
(grant 11871441, 11671369), respectively.

# Efficient multiscale algorithms for simulating nonlocal optical response of
metallic nanostructure arrays

Yongwei Zhang School of Mathematics and Statistics, Zhengzhou University,
Zhengzhou 450001, China (zhangyongwei@lsec.cc.ac.cn). Chupeng Ma  Institute
for Applied Mathematics and Interdisciplinary Center for Scientific Computing,
Heidelberg University, Im Neuenheimer Feld 205, Heidelberg 69120, Germany;
(machupeng@lsec.cc.ac.cn). Liqun Cao Corresponding author. LSEC, NCMIS,
University of Chinese Academy of Sciences, Institute of Computational
Mathematics and Scientific/Engineering Computing, Academy of Mathematics and
Systems Science, Chinese Academy of Sciences, Beijing 100190, China
(clq@lsec.cc.ac.cn). Dongyang Shi School of Mathematics and Statistics,
Zhengzhou University, Zhengzhou 450001, China (shi_dy@zzu.edu.cn).

###### Abstract

In this paper, we consider numerical simulations of the nonlocal optical
response of metallic nanostructure arrays inside a dielectric host, which is
of particular interest to the nanoplasmonics community due to many unusual
properties and potential applications. Mathematically, it is described by
Maxwell’s equations with discontinuous coefficients coupled with a set of
Helmholtz-type equations defined only on the domains of metallic
nanostructures. To solve this challenging problem, we develop an efficient
multiscale method consisting of three steps. First, we extend the system into
the domain occupied by the dielectric medium in a novel way and result in a
coupled system with rapidly oscillating coefficients. A rigorous analysis of
the error between the solutions of the original system and the extended system
is given. Second, we derive the homogenized system and define the multiscale
approximate solution for the extended system by using the multiscale
asymptotic method. Third, to fix the inaccuracy of the multiscale asymptotic
method inside the metallic nanostructures, we solve the original system in
each metallic nanostructure separately with boundary conditions given by the
multiscale approximate solution. A fast algorithm based on the $LU$
decomposition is proposed for solving the resulting linear systems. By
applying the multiscale method, we obtain the results that are in good
agreement with those obtained by solving the original system directly at a
much lower computational cost. Numerical examples are provided to validate the
efficiency and accuracy of the proposed method.

###### keywords:

nonlocal optical response, metallic nanostructure arrays, multiscale
asymptotic method, finite element method.

###### AMS:

65N30, 65N55, 65F10, 65Y05

##  1 Introduction

Metallic nanostructures have attracted great interest of engineers and
scientists in the last two decades due to the phenomenon of localized surface
plasmon resonance (LSPR) [16]. When a metallic nanoparticle is irradiated by
light at optical frequencies, the light wave couples with the collective
oscillation of conduction band electrons at the metal surface and excites
LSPR, leading to many unusual and fascinating physical properties. It enables
the confinement of light at the nanoscale, the enormous local fields
enhancement of the incident light, and the squeezing of light beyond the
diffraction limit. As a powerful platform for manipulation of light at the
nanoscale, metallic nanostructures have found a wide range of applications in
various fields, such as near-field scanning microscopy, ultrasensitive sensing
and detection, solar cells, and related studies have formed a fast-growing and
highly interdisciplinary field – Nanoplasmonics [13].

In order to understand and make use of the phenomenon of LSPR, an appropriate
modeling for describing the interaction between light and metal nanostructures
is required. The most commonly used model for describing the optical response
of metals is the Drude model (later extended into the Lorentz–Drude model)
[5]. This is a classical model based on the kinetic theory of electrons in a
metal which assumes that electrons in the metal do not interact with each
other and move against a positively charged ionic background due to external
fields. The Drude model has achieved a great success in explaining many
physical properties of bulk metals. However, as the size of metals becomes
much smaller than the wavelength of the incident light, nonlocal effects due
to the Pauli exclusion principle become important and the Drude model fails to
explain experimentally observable phenomena. To overcome this, some more
sophisticated material models, such as a nonlinear hydrodynamic model which
couples macroscopic Maxwell’s equations for electromagnetic fields with
hydrodynamic equations for free electrons in metals, are used [6]. The linear-
response approximation of the hydrodynamic model yields the nonlocal
hydrodynamic Drude (NHD) model, which has been a popular model in the
computational study of metallic nanostructures due to its ability to describe
the nonlocal optical response of metal nanoparticles with a high accuracy and
a low computational effort [4, 17, 22]. Mathematically, the NHD model is a
coupled system of Maxwell’s equations for the electric field and a Helmholtz-
type equation for the polarization current.

Much effort has been devoted to solving the NHD model numerically to simulate
the nonlocal optical properties of metallic nanostructures. In [9], a Galerkin
finite element method based on the Nédélec elements is proposed to solve the
NHD model to simulate complex-shaped nano-plasmonic scatters. In [26], a
computational scheme within the framework of a boundary integral equation and
a method of moments algorithm is developed for the NHD model to predict the
interaction of light with metallic nanoparticles. The well-posedness and the
convergence of finite element method for the NHD model are rigorously proved
in [12]. The hybrid discontinuous Galerkin methods for the NHD model have also
been considered [11, 19, 20, 23]. For other more numerical methods for this
coupled system, we refer the reader to [8, 10, 27] and references therein.

It is worth noting that almost all the existing numerical studies on the NHD
model simulate only one or several metallic nanostructures. Compared with
individual and small clusters of nanoparticles, when arranged into
multidimensional arrays, metallic nanoparticles strongly interact with each
other such that the collective properties can be rationally designed by
changing the inter-particle spacing, array period, and overall composition,
providing unparalleled opportunities for realizing materials with interesting
and unusual photonic and metamaterial properties [18, 24]. Indeed, metallic
nanoparticle arrays have found applications in plasmonic solar cells [14] and
optical sensing [7] and the literature of related studies is fast growing in
recent years. However, due to the huge computational cost, numerical
simulations of optical properties of metallic nanoparticle arrays are very
limited, especially in the case of three-dimensional arrays. A simple and
widely used method for describing the macroscopic electromagnetic properties
of composite materials is given by the Maxwell–Garnett effective medium theory
(EMT) [3] in which the effective permittivity is determined by the volume
fraction of the inclusions, the permittivity of the inclusions, and the
background material. However, in general, the Maxwell–Garnett formula is only
valid at low volume fraction. In the case of plasmon resonance, it has been
reported that the Maxwell–Garnett approach is correct only at volume fraction
of the inclusions lower than $10^{-5}$ [1]. Furthermore, in the case
considered in this paper, due to the spatially-dispersive permittivity of the
metallic nanostructures, the application of the Maxwell–Garnett formula will
lead to a much more complicated effective system. Therefore, in order to
simulate the nonlocal optical response of realistic metallic nanoparticle
arrays, some accurate and efficient methods beyond the Maxwell–Garnett formula
are necessary.

In this paper, we consider numerical simulations of the nonlocal optical
response of metallic nanostructure arrays embedded in a dielectric medium with
the NHD model. Mathematically, it involves solving a coupled system of
Maxwell’s equations defined on the whole domain and a set of Helmholtz-type
equations defined in each metallic nanostructure separately. It is a
challenging task to solve this system numerically, mainly due to two reasons.
First, the coefficients of Maxwell’s equations are rapidly oscillating, whose
numerical solution is computationally expensive (or even impossible) in
practice with a standard method, especially when there are a huge number of
metallic nanostructures. Second, in contrast to many other coupled systems,
this system is coupled in a novel way in which an equation defined on the
whole domain is coupled with another equation defined on many unconnected
domains, making it difficult to apply the classical multiscale asymptotic
method directly. In order to solve this system efficiently, some special
strategies to deal with such a coupling are necessary.

In this paper, we develop an efficient multiscale method to solve this coupled
system, which consists of three steps. First, we extend the equation satisfied
by the polarization current into the domain occupied by the dielectric medium
by setting some coefficients of the equation sufficiently large outside the
metallic nanostructures. The error between the solutions of the original
system and the extended system is rigorously proved. In this way, we obtain a
coupled system defined on the whole domain of the scatter with rapidly
oscillating coefficients. Second, we use the multiscale asymptotic method to
derive the homogenized system and define the multiscale approximate solution
for the extended system. The coefficients of the homogenized system are
constant, making it much easier to solve. In view of the fact that the
polarization current has little effect on the homogenized coupled system in
practice, we propose to solve the homogenized Maxwell’s equations rather than
the homogenized coupled system to further reduce the computational burden.
Numerical results show that the multiscale asymptotic method produces
satisfactory results outside the metallic nanostructures while it fails to
capture the oscillations of the electric field inside the metallic
nanostructures. Third, in order to remedy the failure of the multiscale
asymptotic method inside the metallic nanostructures, we solve the original
system in each metallic nanostructure separately with boundary conditions
given by the multiscale approximate solution. Since now the system is only
solved in each metallic nanostructure with constant coefficients, the
computational cost is much lower. In addition, by taking advantage of the
periodic arrangement of metallic nanostructures, we propose a fast algorithm
based on the $LU$ decomposition for solving the resulting linear systems to
further reduce the computational cost and make the algorithm applicable to the
system with a huge number of metallic nanostructures. By these three steps, we
get the results agreeing well with those obtained from solving the original
system directly with a much lower computational effort. The multiscale method
developed in this paper might be able to be generalized to simulate the
optical properties of composite materials with the non-dispersive media matrix
and the dispersive media inclusions, such as the Drude, Lorentz–Drude and
Debye media.

The rest of this paper is organized as follows. In section 2, we briefly
introduce the NHD model and describe the problem considered in this paper. The
extended system is given at the end of this section. In section 3, we derive
the homogenized system and define the multiscale approximate solution for the
extended system by using the technique of multiscale asymptotic expansion. In
section 4, we describe the multiscale approaches and the associated numerical
algorithms in detail. In section 5, we present some numerical experiments to
demonstrate the accuracy and efficiency of our method.

##  2 Problem Formulation

Throughout this paper, we use $\mathbf{curl}$, $\mathbf{div}$, and
$\mathbf{grad}$ to denote the curl, divergence, and gradient operators,
respectively.

###  2.1 Nonlocal hydrodynamic Drude model

The NHD model is derived via the linearization of the hydrodynamic model [4,
17]. Within the hydrodynamic model, the free electrons in a metal are
considered as a charged and compressible fluid which is described in terms of
the charge density $n(\mathbf{x},t)$, the electron fluid velocity ${\bf
v}(\mathbf{x},t)$, and the electron pressure $p(\mathbf{x},t)$. Under the
influence of macroscopic electromagnetic fields, the motion of the electron
fluid is determined by the Euler equation

${\displaystyle m_{e}(\partial_{t}+{\bf v}\cdot\mathbf{grad}+\gamma){\bf
v}=-e({\bf E}+\mu_{0}{\bf v}\times{\bf H})-\frac{\mathbf{grad}\;p}{n}},$ (1)

along with the continuity equation

${\displaystyle\partial_{t}n+\mathbf{div}(n{\bf v})=0,}$ (2)

where $-e$ is the electron charge, $m_{e}$ is the effective electron mass, and
$\gamma>0$ is the damping constant. The term $-e({\bf E}+\mu_{0}{\bf
v}\times{\bf H})$ represents the Lorentz force (assuming non-magnetic
materials) in which the electric field $\mathbf{E}$ and the magnetic field
$\mathbf{H}$ satisfy the macroscopic Maxwell’s equations

$\left\\{\begin{array}[]{lll}{\displaystyle\mathbf{curl}\;{\bf
E}=-{\mu_{0}}\partial_{t}{\bf H},}\\\\[5.69054pt]
{\displaystyle\mathbf{curl}\;{\bf
H}=\varepsilon_{0}\varepsilon_{\infty}\partial_{t}{\bf E}+{\bf
J}.}\end{array}\right.$ (3)

Here $\mu_{0}$ and ${\varepsilon_{0}}$ are the magnetic permeability and
electric permittivity of vacuum, respectively, and $\varepsilon_{\infty}$ is
the relative electric permittivity of the metal that takes into account the
polarization of bound electrons. The hydrodynamic equations (1)-(2) are
coupled with Maxwell’s equations (3) via the polarization current density
$\mathbf{J}$ of the free electrons:

$\mathbf{J}=-en\mathbf{v}.$ (4)

The above equations (1)-(4) constitute the hydrodynamic model containing
nonlinear and nonlocal effects. In order to get a simplified system of
equations more suitable for numerical computation, we linearize the
hydrodynamic model in a perturbative manner. The physical fields are expanded
into a non-oscillating term (e.g. the constant equilibrium electron density
$n_{0}$) and a small first-order dynamic term as follows:

$n({\bf x},t)\approx n_{0}+n_{1}({\bf x},t),\quad{\bf v}({\bf x},t)\approx{\bf
v}_{0}+{\bf v}_{1}({\bf x},t).$ (5)

Similar expansions can be written for the electric and magnetic fields. Since
in the absence of external fields ${\bf v}={\bf v}_{0}={\bf 0}$, the two
nonlinear (higher-order) terms ${\bf v}\cdot{\bf grad}\,{\bf v}$ and ${\bf
v}\times{\bf H}$ in (1) are dropped out due to linearization. By using the
Thomas-Fermi theory of metals [17], the electron pressure $p(\mathbf{x},t)$ is
given by

$p(\mathbf{x},t)=(3\pi^{2})^{2/3}\frac{\hbar^{2}}{5m_{e}}n({\bf x},t)^{5/3}.$
(6)

In view of (5) and (6), the term $\displaystyle\frac{{\bf grad}\,p}{n}$ in (1)
can be linearized as

$\frac{\mathbf{grad}\;p}{n}\approx
m_{e}\beta^{2}\frac{\mathbf{grad}\;n}{n_{0}},$ (7)

where
$\displaystyle\beta^{2}=\frac{(3\pi)^{2/3}}{3}\frac{\hbar^{2}}{m_{e}^{2}}n_{0}^{2/3}$
is a parameter representing the nonlocality. Using the approximation (7) and
neglecting the two nonlinear terms in (1) due to linearization, we arrive at
the following linearized Euler equation

$\partial_{t}{\bf v}=\frac{-e}{m_{e}}{\bf E}-\gamma{\bf
v}-\beta^{2}\frac{\mathbf{grad}\;n}{n_{0}}$ (8)

and the linearized continuity equation

$\partial_{t}n+n_{0}\mathbf{div}\;{\bf v}=0.$ (9)

Differentiating (8) with respect to time $t$, inserting the linearized current
density ${\bf J}\approx-en_{0}{\bf v}$ and making use of (9), we come to

$\partial_{tt}{\bf J}+\gamma\partial_{t}{\bf
J}-\beta^{2}\mathbf{grad}\,(\mathbf{div}\,{\bf
J})-\omega_{p}^{2}\varepsilon_{0}\partial_{t}{\bf E}=0,$ (10)

where $\omega_{p}=\sqrt{n_{0}e^{2}/(m_{e}\varepsilon_{0})}$ is the plasma
frequency. Combining (3) and (10), we obtain the NHD model for metals

$\left\\{\begin{array}[]{@{}l@{}}{\displaystyle\mathbf{curl}\;{\bf
E}=-\mu_{0}\partial_{t}{\bf H},}\\\\[5.69054pt]
{\displaystyle\mathbf{curl}\;{\bf
H}=\varepsilon_{0}\varepsilon_{\infty}\partial_{t}{\bf E}+{\bf
J},}\\\\[5.69054pt] {\displaystyle\partial_{tt}{\bf J}+\gamma\partial_{t}{\bf
J}-\beta^{2}\mathbf{grad}\,(\mathbf{div}\,{\bf
J})-\omega_{p}^{2}\varepsilon_{0}\partial_{t}{\bf E}=0.}\end{array}\right.$
(11)

Replacing $\partial_{t}$ with $-{\rm i}\omega$ in (11) by Fourier
transformation in the time domain, where ${\rm i}$ is the imaginary unit and
$\omega$ is the angular frequency, and eliminating the magnetic field ${\bf
H}$, we obtain the NHD model in the frequency domain

$\left\\{\begin{array}[]{@{}l@{}}{\displaystyle\mathbf{curl}\;(\mu_{0}^{-1}\mathbf{curl}\;{\bf
E})-\varepsilon_{0}\varepsilon_{\infty}\omega^{2}{\bf E}={\rm i}\omega{\bf
J},}\\\\[5.69054pt] {\displaystyle\omega(\omega+{\rm i}\gamma){\bf
J}+\beta^{2}\mathbf{grad}\,(\mathbf{div}\,{\bf J})={\rm
i}\omega\omega^{2}_{p}\varepsilon_{0}{\bf E}.}\end{array}\right.$ (12)

######  Remark 2.1.

By Fourier transformation in the space domain, we replace the operator ${\bf
grad}\,{\bf div}$ with $-{\bf k}^{2}$ in the second equation of (12) and get

${\bf J}=\frac{{\rm i}\omega\omega^{2}_{p}\varepsilon_{0}\,{\bf
E}}{{\omega(\omega+{\rm i}\gamma)-\beta^{2}{\bf k}^{2}}},$ (13)

which leads to the spatially-dispersive (relative) permittivity for the metal

$\epsilon(\omega,{\bf
k})=\epsilon_{\infty}-\frac{\omega_{p}^{2}}{\omega(\omega+{\rm
i}\gamma)-\beta^{2}{\bf k}^{2}}.$ (14)

The parameter $\beta$ represents the level of nonlocality. As
$\beta\rightarrow 0$, the NHD model reduces to the classical local-response
Drude model.

###  2.2 Electromagnetic scattering by metallic nanostructure arrays

In this paper, we consider the problem as shown in Figure 2.1 where the
metallic nanostructure arrays embedded in a dielectric medium are illuminated
by an incident plane wave. The interaction of light with metallic
nanostructure arrays is described by the NHD model

$\left\\{\begin{array}[]{@{}l@{}}{\displaystyle\mathbf{curl}\,(\mu_{\eta}^{-1}\,({\bf
x})\mathbf{curl}\,{\bf E}_{\eta})-\varepsilon_{\eta}({\bf x})\omega^{2}\,{\bf
E}_{\eta}-{\mathrm{{i}}}\omega{\bf J}_{\eta}={\bf 0},\,\;\quad{\rm
in}\;\;\Omega_{s},}\\\\[5.69054pt]
{\displaystyle\omega(\omega+{\mathrm{{i}}}\gamma){\bf
J}_{\eta}+\beta^{2}\mathbf{grad}\,(\mathbf{div}\,{\bf J}_{\eta})-{\rm
i}\omega\omega^{2}_{p}\varepsilon_{0}{\bf E}_{\eta}={\bf 0},\;\quad{\rm
in}\;\;\Omega_{s,\eta}.}\end{array}\right.$ (15)

Here $\Omega_{s}\subset\mathbb{R}^{3}$ is a bounded Lipschitz domain occupied
by the scatter (the metallic nanostructures and the dielectric medium).
$\Omega_{s,\eta}=\cup_{k=1}^{N}\Omega_{\eta}^{k}$ denotes the domain occupied
by the metallic nanostructures, where $\Omega_{\eta}^{k}$ is the domain
occupied by the $k$-th nanostructure and $N$ is the number of metallic
nanostructures. $\eta$ is the relative size of the periodic microstructure and
we denote by $Y=(0,l_{1})\times(0,l_{2})\times(0,l_{3})$ the rescaled
reference cell, where $l_{1}$, $l_{2}$, and $l_{3}$ are positive constants.
$\mu_{\eta}(\mathbf{x})$ and $\varepsilon_{\eta}(\mathbf{x})$ are the magnetic
permeability and electric permittivity, respectively. We assume that
$\mu_{\eta}(\mathbf{x})=\mu(\mathbf{x}/\eta)$,
$\varepsilon_{\eta}(\mathbf{x})=\varepsilon(\mathbf{x}/\eta)$, where $\mu({\bf
y})$ and $\varepsilon({\bf y})$ satisfy

$(\mathbf{A}_{1}).$ $\mu({\bf y})$ and $\varepsilon({\bf y})$ are $Y$-periodic
in $\mathbb{R}^{3}$, i.e.,

$\mu({\bf y}+kl_{i}\mathbf{e}_{i})=\mu({\bf y}),\quad\varepsilon({\bf
y}+kl_{i}\mathbf{e}_{i})=\varepsilon({\bf y}),\quad\forall
k\in\mathbb{Z},\quad i=1,2,3,\quad{\bf y}\in\mathbb{R}^{3},$

and $\\{\mathbf{e}_{1},\,\mathbf{e}_{2},\,\mathbf{e}_{3}\\}$ is the canonical
basis of $\mathbb{R}^{3}$.

$(\mathbf{A}_{2}).$ There exists positive constants $\mu_{1}$, $\mu_{2}$,
$\varepsilon_{1}$, and $\varepsilon_{2}$ such that $\mu({\bf y})$ and
$\varepsilon({\bf y})$ satisfy

$0<\mu_{1}\leq\mu({\bf y})\leq\mu_{2}<\infty,\quad
0<\varepsilon_{1}\leq\varepsilon({\bf
y})\leq\varepsilon_{2}<\infty,\quad\forall{\bf y}\in Y.$

######  Remark 2.2.

In the simple case as illustrated in Fig 2.1, where there is only one metallic
nanostructure inside the reference cell, $\mu({\bf y})$ and $\varepsilon({\bf
y})$ are piecewise constant functions given by

$\displaystyle\mu({\bf
y})=\left\\{\begin{array}[]{lll}\mu_{m},\,\quad\mathbf{y}\in
Y_{1}\\\\[8.53581pt] \mu_{d},\,\,\,\quad\mathbf{y}\in
Y_{2},\end{array}\right.\quad\varepsilon({\bf
y})=\left\\{\begin{array}[]{l}\varepsilon_{m},\,\quad\mathbf{y}\in
Y_{1}\\\\[8.53581pt] \varepsilon_{d},\,\;\quad\mathbf{y}\in
Y_{2},\end{array}\right.$

where $Y_{1}\subset Y$ and $Y_{2}\subset Y$ are the domains occupied by the
metallic nanostructure and the dielectric medium, respectively,
$\mu_{m}\,$($\varepsilon_{m}$) and $\mu_{d}\,$($\varepsilon_{d}$) are the
magnetic permeability (electric permittivity) of the metal and dielectric
medium, respectively.

Plane waveexcitation Fig. 1: Scattering setting: an incident electromagnetic
wave is incident on metallic nanostructure arrays embedded in a dielectric
medium. We assume that the dielectric medium is surrounded by free space. The
rescaled reference cell is plotted on the right.

Note that the polarization current $\mathbf{J}_{\eta}$ only exists in the
metallic nanostructures. For the well-posedness of the problem, we need to
impose some appropriate boundary conditions on $\partial\Omega_{s,\eta}$ for
$\mathbf{J}_{\eta}$. In this paper, we use the hard-wall boundary conditions

${\displaystyle{\bf n}\cdot{\bf J}_{\eta}=0,\quad{\rm
on}\;\;\partial\Omega_{s,\eta},}$ (16)

which implies that the electrons are confined within the metal and spill-out
of electrons outside the metal is neglected.

We assume that the scatter is surrounded by free space in which the (total)
electric field satisfies the following Maxwell’s equations

${\displaystyle\mathbf{curl}\,(\mu_{0}^{-1}\,\mathbf{curl}\,\mathbf{E})-\varepsilon_{0}\omega^{2}\,\mathbf{E}=\mathbf{0}}$
(17)

and the scattered field $\mathbf{E}_{s}=\mathbf{E}-\mathbf{E}^{inc}$ satisfies
the Silver-Müller radiation condition

${\displaystyle(\mu_{0}^{-1}\,\mathbf{curl}\,{\bf E}_{s})\times\hat{\bf
x}-{\mathrm{{i}}}\omega{\bf E}_{s}=o(1/|{\mathbf{x}}|),\quad{\rm
as}\quad|\mathbf{x}|\rightarrow\infty}.$ (18)

Here $\mathbf{E}^{inc}$ denotes the electric field of the incident wave. In
addition, we need to impose the transmission conditions across the boundary
$\partial\Omega_{s}$ for the electric field

$\mathbf{n}\times\mathbf{E}_{\eta}=\mathbf{n}\times\mathbf{E},\quad\mathbf{n}\times(\mu_{\eta}^{-1}\mathbf{curl}\,\mathbf{E}_{\eta})=\mathbf{n}\times(\mu_{0}^{-1}\mathbf{curl}\,\mathbf{E}).$
(19)

The scattering problem (15)-(19) is defined on the whole space. To make this
unbounded electromagnetic problem suitable for numerical computations, it is
necessary to truncate the infinite domain. Different techniques can be used,
such as absorbing boundary conditions, perfectly matched layer (PML), and
boundary integral equation techniques. In this paper, we apply the following
absorbing boundary condition on the boundary of a much larger domain $\Omega$
that contains the domain $\Omega_{s}$ of the scatter

${\displaystyle(\mu_{0}^{-1}\,\mathbf{curl}\,{\bf E})\times{\bf
n}-{\mathrm{{i}}}\omega({\bf n}\times{\bf E})\times{\bf n}={\bf g},\quad{\rm
on}\;\;\partial\Omega,}$ (20)

where ${\bf g}({\bf x})=(\mu_{0}^{-1}\,\mathbf{curl}\,{\bf E}^{inc})\times{\bf
n}-{\mathrm{{i}}}\omega({\bf n}\times{\bf E}^{inc})\times{\bf n}$. (20) is the
first-order approximation of the Silver-Müller radiation condition [21]. It is
important to note that the multiscale method developed in this paper is
independent of the choice of domain truncation techniques.

Combining (15)-(20), Maxwell’s equations coupled with the NHD model which
describe electromagnetic scattering by metallic nanostructure arrays embedded
in a dielectric host are given by

$\left\\{\begin{array}[]{@{}l@{}}{\displaystyle\mathbf{curl}\,(\mu_{\eta}^{-1}\,({\bf
x})\mathbf{curl}\,{\bf E}_{\eta})-\varepsilon_{\eta}({\bf x})\omega^{2}\,{\bf
E}_{\eta}-{\mathrm{{i}}}\omega{\bf J}_{\eta}={\bf 0},\quad\;{\rm
in}\;\;\Omega_{s},}\\\\[5.69054pt]
{\displaystyle\omega(\omega+{\mathrm{{i}}}\gamma){\bf
J}_{\eta}+\beta^{2}\mathbf{grad}\,(\mathbf{div}\,{\bf J}_{\eta})-{\rm
i}\omega\omega^{2}_{p}\varepsilon_{0}{\bf E}_{\eta}={\bf 0},\;\quad{\rm
in}\;\Omega_{s,\eta},}\\\\[5.69054pt]
{\displaystyle\mathbf{curl}\,(\mu_{0}^{-1}\,\mathbf{curl}\,\mathbf{E})-\varepsilon_{0}\omega^{2}\,\mathbf{E}=\mathbf{0},\qquad\qquad\qquad\quad\;\;\quad{\rm
in}\;\;\Omega/\Omega_{s},}\\\\[5.69054pt]
{\displaystyle\mathbf{n}\times(\mu_{\eta}^{-1}\mathbf{curl}\,\mathbf{E}_{\eta})=\mathbf{n}\times(\mu_{0}^{-1}\mathbf{curl}\,\mathbf{E}),\qquad{\rm
on}\;\partial\Omega_{s},}\\\\[5.69054pt]
{\displaystyle\mathbf{n}\times\mathbf{E}_{\eta}=\mathbf{n}\times\mathbf{E},\qquad\qquad\quad\qquad\qquad\;\qquad{\rm
on}\;\partial\Omega_{s},}\\\\[5.69054pt] {\displaystyle{\bf n}\cdot{\bf
J}_{\eta}=0,\qquad\qquad\;\;\,\qquad\qquad\qquad\;\qquad\quad{\rm
on}\;\;\partial\Omega_{s,\eta},}\\\\[5.69054pt]
{\displaystyle(\mu_{0}^{-1}\,\mathbf{curl}\,{\bf E})\times{\bf
n}-{\mathrm{{i}}}\omega({\bf n}\times{\bf E})\times{\bf n}={\bf
g},\;\;\quad{\rm on}\;\;\partial\Omega.}\end{array}\right.$ (21)

Due to the rapidly oscillating coefficients of Maxwell’s equations, it is
usually very difficult to solve the coupled system (21) directly by a standard
numerical method. In general, this difficulty can be overcome by using the
homogenization and multiscale asymptotic methods. However, an interesting
feature of this coupled system is that the equation satisfied by the
polarization current is only defined in the metallic nanostructures while
couples with Maxwell’s equations defined on the whole domain, which prevents
the application of the usual multiscale asymptotic method. To tackle this
problem, we extend the equation satisfied by the polarization current outside
the metallic nanostructures in a novel way.

First, we introduce the $Y$-periodic functions $\gamma_{\lambda}(\mathbf{y})$
and $\beta^{2}_{\lambda}(\mathbf{y})$ which are given by

$\displaystyle\gamma_{\lambda}(\mathbf{y})=\left\\{\begin{array}[]{lll}\gamma,\,\,\quad\mathbf{y}\in
Y_{1}\\\\[8.53581pt] \lambda,\,\,\,\quad\mathbf{y}\in
Y_{2},\end{array}\right.\quad\beta^{2}_{\lambda}(\mathbf{y})=\left\\{\begin{array}[]{l}\beta^{2},\,\quad\mathbf{y}\in
Y_{1}\\\\[8.53581pt] \lambda,\,\,\,\;\quad\mathbf{y}\in
Y_{2},\end{array}\right.$ (26)

in the reference cell, where $Y_{1}\subset Y$ and $Y_{2}\subset Y$ are the
domains occupied by metallic nanostructures and the dielectric medium,
respectively, and $\lambda>0$ is the extension parameter which is assumed to
be sufficiently large.

Next we set $\gamma_{\eta,\lambda}(\mathbf{x})=\gamma_{\lambda}({\bf x}/\eta)$
and
$\beta^{2}_{\eta,\lambda}(\mathbf{x})=\beta^{2}_{\lambda}(\mathbf{x}/\eta)$,
which satisfy

$\displaystyle\gamma_{\eta,\lambda}(\mathbf{x})=\left\\{\begin{array}[]{l}\gamma,\,\quad\mathbf{x}\in\Omega_{s,\eta}\\\\[8.53581pt]
\lambda,\,\quad\mathbf{x}\in\Omega_{s}/\Omega_{s,\eta},\end{array}\right.\quad\beta^{2}_{\eta,\lambda}(\mathbf{x})=\left\\{\begin{array}[]{l}\beta^{2},\,\quad\mathbf{x}\in\Omega_{s,\eta}\\\\[8.53581pt]
\lambda,\,\;\,\quad\mathbf{x}\in\Omega_{s}/\Omega_{s,\eta},\end{array}\right.$
(31)

respectively.

Finally, we extend the equation satisfied by the polarization current into the
whole of $\Omega_{s}$ by replacing the constants $\gamma$ and $\beta^{2}$ in
the equation by $\gamma_{\eta,\lambda}(\mathbf{x})$ and
$\beta^{2}_{\eta,\lambda}(\mathbf{x})$, respectively. We denote by
$({\mathbf{J}}_{\eta,\lambda},{\mathbf{E}}_{\eta,\lambda},{\mathbf{E}}_{\lambda})$
the extended solution, which satisfies the extended system

$\left\\{\begin{array}[]{@{}l@{}}{\displaystyle\mathbf{curl}\,(\mu_{\eta}^{-1}\mathbf{curl}\,{\bf
E}_{\eta,\lambda})-\varepsilon_{\eta}\omega^{2}\,{\bf
E}_{\eta,\lambda}-{\mathrm{{i}}}\omega{\bf J}_{\eta,\lambda}={\bf
0},\qquad\quad\qquad\quad\;{\rm in}\;\;\Omega_{s},}\\\\[5.69054pt]
{\displaystyle\omega(\omega+{\mathrm{{i}}}\gamma_{\eta,\lambda}){\bf
J}_{\eta,\lambda}+\mathbf{grad}\,(\beta^{2}_{\eta,\lambda}\mathbf{div}\,{\bf
J}_{\eta,\lambda})-{\rm i}\omega\omega^{2}_{p}\varepsilon_{0}{\bf
E}_{\eta,\lambda}={\bf 0},\;\quad{\rm in}\;\;\Omega_{s},}\\\\[5.69054pt]
{\displaystyle\mathbf{curl}\,(\mu_{0}^{-1}\,\mathbf{curl}\,{\mathbf{E}}_{\lambda})-\varepsilon_{0}\omega^{2}\,{\mathbf{E}}_{\lambda}=\mathbf{0},\qquad\qquad\quad\qquad\qquad\quad\;\;\quad{\rm
in}\;\;\Omega/\Omega_{s},}\\\\[5.69054pt]
{\displaystyle\mathbf{n}\times(\mu_{\eta}^{-1}\mathbf{curl}\,{\mathbf{E}}_{\eta,\lambda})=\mathbf{n}\times(\mu_{0}^{-1}\mathbf{curl}\,{\mathbf{E}}_{\lambda}),\,\qquad{\rm
on}\;\partial\Omega_{s},}\\\\[5.69054pt]
{\displaystyle\mathbf{n}\times{\mathbf{E}}_{\eta,\lambda}=\mathbf{n}\times{\mathbf{E}}_{\lambda},\qquad\qquad\quad\qquad\qquad\;\qquad{\rm
on}\;\partial\Omega_{s},}\\\\[5.69054pt] {\displaystyle{\bf n}\cdot{\bf
J}_{\eta,\lambda}=0,\qquad\quad\qquad\;\,\qquad\qquad\qquad\;\qquad\quad{\rm
on}\;\,\partial\Omega_{s},}\\\\[5.69054pt]
{\displaystyle(\mu_{0}^{-1}\,\mathbf{curl}\,{\bf E}_{\lambda})\times{\bf
n}-{\mathrm{{i}}}\omega({\bf n}\times{\bf E}_{\lambda})\times{\bf n}={\bf
g},\;\;\;\quad{\rm on}\;\;\partial\Omega.}\end{array}\right.$ (32)

We have the following error estimate of the extension:

######  Theorem 1.

Let $({\mathbf{J}}_{\eta},{\mathbf{E}}_{\eta},{\mathbf{E}})$ and
$({\mathbf{J}}_{\eta,\lambda},{\mathbf{E}}_{\eta,\lambda},{\mathbf{E}}_{\lambda})$
be the solutions of the original system (21) and the extended system (32),
respectively. Assume that $\mathbf{g}\in\mathbf{L}^{2}(\partial\Omega)$ and
$\Omega$, $\Omega_{s}$, and $\Omega_{s,\eta}$ are bounded, simply-connected,
Lipschitz domains in $\mathbb{R}^{3}$ with $\bar{\Omega}_{s}\subset\Omega$ and
$\bar{\Omega}_{s,\eta}\subset\Omega_{s}$. It holds for any $\lambda>0$,

$\begin{array}[]{lll}{\displaystyle\|{\mathbf{J}}_{\eta}-{\mathbf{J}}_{\eta,\lambda}\|_{\mathbf{H}(\mathbf{div};\Omega_{s,\eta})}+\|{\mathbf{E}}_{\eta}-{\mathbf{E}}_{\eta,\lambda}\|_{\mathbf{H}(\mathbf{curl};\Omega_{s})}}\\\\[5.69054pt]
{\displaystyle\qquad\quad+\,\|{\mathbf{E}}-{\mathbf{E}}_{\lambda}\|_{\mathbf{H}(\mathbf{curl};\Omega/\Omega_{s})}\leq{C}/{\lambda^{\frac{1}{2}}},}\\\\[5.69054pt]
\end{array}$ (33)

where $C$ is a positive constant independent of $\lambda$ but might depend on
$\eta$.

######  Remark 2.3.

Theorem 1 indicates that the error between the original system and the
extended system decreases with increasing extension parameter $\lambda$.

The proof of Theorem 1 is given in the Appendix.

##  3 Homogenization and Multiscale asymptotic method

In this section, we derive the homogenized system and define the multiscale
approximate solutions for the extended system (32) by using the technique of
multiscale asymptotic expansion. For brevity, we omit the dependence of the
electric field and the polarization current on the parameter $\lambda$ and
write $({\bf E}_{\eta,\lambda},{\bf J}_{\eta,\lambda})$ as $({\bf
E}_{\eta},{\bf J}_{\eta})$ in what follows. We consider the homogenization of
the following equations

$\left\\{\begin{array}[]{@{}l@{}}{\displaystyle\mathbf{curl}\,(\mu_{\eta}^{-1}\mathbf{curl}\,{\bf
E}_{\eta})-\varepsilon_{\eta}\omega^{2}{\bf E}_{\eta}={\mathrm{{i}}}\omega{\bf
J}_{\eta},}\\\\[5.69054pt] {\displaystyle\gamma^{\ast}_{\eta}{\bf
J}_{\eta}+\mathbf{grad}\,(\beta^{\ast}_{\eta}\mathbf{div}\,{\bf
J}_{\eta})={\rm i}\omega{\bf E}_{\eta},}\end{array}\right.$ (34)

where

$\mu_{\eta}=\mu({\bf x}/\eta),\quad\varepsilon_{\eta}=\varepsilon({\bf
x}/\eta),\quad\gamma^{\ast}_{\eta}=\frac{\omega(\omega+{\rm i}\gamma({\bf
x}/\eta))}{\omega_{p}^{2}\varepsilon_{0}},\quad\beta^{\ast}_{\eta}=\frac{\beta^{2}({\bf
x}/\eta)}{\omega_{p}^{2}\varepsilon_{0}},$ (35)

and $\mu({\bf y})$, $\varepsilon({\bf y})$, $\gamma({\bf y})$, $\beta({\bf
y})$ are $Y$-periodic functions. We want to study the behaviour of
$({\bf{E}}_{\eta},\bf J_{\eta})$ as $\eta\to 0$.

For the sake of simplicity, we introduce the operators

$\begin{array}[]{@{}c@{}}\mathcal{A}^{\eta}={\mathbf{curl}}\,\mu_{\eta}^{-1}\,{\mathbf{curl}}-\varepsilon_{\eta}\omega^{2},\quad\mathcal{B}^{\eta}={\mathbf{grad}}\,\beta^{\ast}_{\eta}\,{\mathbf{div}}+\gamma^{\ast}_{\eta}.\end{array}$

We look for a formal asymptotic expansion of the form

$\begin{array}[]{@{}l@{}}{\bf{E}}_{\eta}={\bf{E}}_{0}({\bf x},{\bf
x}/\eta)+\eta{\bf{E}}_{1}({\bf x},{\bf x}/\eta)+\eta^{2}{\bf{E}}_{2}({\bf
x},{\bf x}/\eta)+...,\\\\[5.69054pt] {\bf{J}}_{\eta}={\bf{J}}_{0}({\bf x},{\bf
x}/\eta)+\eta{\bf{J}}_{1}({\bf x},{\bf x}/\eta)+\eta^{2}{\bf{J}}_{2}({\bf
x},{\bf x}/\eta)+...,\end{array}$ (36)

where for $k=0,1,\cdots$, ${\bf{E}}_{k}({\bf x},{\bf y})$ and
${\bf{J}}_{k}({\bf x},{\bf y})$ are functions of both variables ${\bf x}$ and
${\bf y}$ and $Y$-periodic with respect to ${\bf y}$.

Let ${\bm{\Psi}}={\bm{\Psi}}({\bf x},{\bf y})\in\mathbb{R}^{3}$ be a function
depending on two variables of $\mathbb{R}^{3}$ and denote by
${\bm{\Psi}}_{\eta}={\bm{\Psi}}({\bf x},{\bf x}/\eta)$. By using the chain
rule, we have

$\begin{array}[]{lll}{\displaystyle\mathbf{curl}\,{\bm{\Psi}}_{\eta}({\bf
x})=\big{[}\mathbf{curl}_{\bf x}\,{\bm{\Psi}}+\frac{1}{\eta}\mathbf{curl}_{\bf
y}\,{\bm{\Psi}}\big{]}({\bf x},\frac{\bf x}{\eta})}\\\\[5.69054pt]
{\displaystyle\mathbf{div}\,{\bm{\Psi}}_{\eta}({\bf
x})=\big{[}\mathbf{div}_{\bf x}\,{\bm{\Psi}}+\frac{1}{\eta}\mathbf{div}_{\bf
y}\,{\bm{\Psi}}\big{]}({\bf x},\frac{\bf x}{\eta})}.\end{array}$

With this in mind, we can write $\mathcal{A}^{\eta}{\bm{\Psi}}_{\eta}$ and
$\mathcal{B}^{\eta}{\bm{\Psi}}_{\eta}$ as follows:

$\begin{array}[]{lll}{\displaystyle\mathcal{A}^{\eta}{\bm{\Psi}}_{\eta}({\bf
x})=\big{[}\big{(}\eta^{-2}\mathcal{A}_{1}+\eta^{-1}\mathcal{A}_{2}+\mathcal{A}_{3}\big{)}{\bm{\Psi}}\big{]}({\bf
x},\frac{\bf x}{\eta}),}\\\\[8.53581pt]
{\displaystyle\mathcal{B}^{\eta}{\bm{\Psi}}_{\eta}({\bf
x})=\big{[}\big{(}\eta^{-2}\mathcal{B}_{1}+\eta^{-1}\mathcal{B}_{2}+\mathcal{B}_{3}\big{)}{\bm{\Psi}}\big{]}({\bf
x},\frac{\bf x}{\eta}),}\end{array}$ (37)

where

$\left\\{\begin{array}[]{@{}l@{}}{\displaystyle\mathcal{A}_{1}=\mathbf{curl}_{\bf
y}\,\mu^{-1}({\bf y})\,\mathbf{curl}_{\bf y}}\\\\[5.69054pt]
{\displaystyle\mathcal{A}_{2}=\mathbf{curl}_{\bf x}\,\mu^{-1}({\bf
y})\,\mathbf{curl}_{\bf y}+\mathbf{curl}_{\bf y}\,\mu^{-1}({\bf
y})\,\mathbf{curl}_{\bf x}}\\\\[5.69054pt]
{\displaystyle\mathcal{A}_{3}=\mathbf{curl}_{\bf x}\,\mu^{-1}({\bf
y})\,\mathbf{curl}_{\bf x}-\varepsilon({\bf y})\omega^{2},}\end{array}\right.$
(38)

and

$\left\\{\begin{array}[]{@{}l@{}}{\displaystyle\mathcal{B}_{1}=\mathbf{grad}_{\bf
y}\,\beta^{\ast}({\bf y})\,\mathbf{div}_{\bf y}}\\\\[5.69054pt]
{\displaystyle\mathcal{B}_{2}=\mathbf{grad}_{\bf x}\,\beta^{\ast}({\bf
y})\,\mathbf{div}_{\bf y}+\mathbf{grad}_{\bf y}\,\beta^{\ast}({\bf
y})\,\mathbf{div}_{\bf x}}\\\\[5.69054pt]
{\displaystyle\mathcal{B}_{3}=\mathbf{grad}_{\bf x}\,\beta^{\ast}({\bf
y})\,\mathbf{div}_{\bf x}+\gamma^{\ast}({\bf y}).}\end{array}\right.$ (39)

Next substituting the expansions (36) into the equations (34), using (37)-(39)
and equating the power-like terms of $\eta$, we get the following equations

$\begin{array}[]{@{}l@{}}{\displaystyle\eta^{-2}:\qquad\mathcal{A}_{1}{\bf{E}}_{0}=\bf{0}}\\\\[5.69054pt]
{\displaystyle\eta^{-1}:\qquad\mathcal{A}_{1}{\bf{E}}_{1}+\mathcal{A}_{2}{\bf{E}}_{0}=\bf{0}}\\\\[5.69054pt]
{\displaystyle\;\,\eta^{0}:\qquad\mathcal{A}_{1}{\bf{E}}_{2}+\mathcal{A}_{2}{\bf{E}}_{1}+\mathcal{A}_{3}{\bf{E}}_{0}={\rm
i}\omega{\bf{J}}_{0},}\end{array}$ (40)

and

$\begin{array}[]{@{}l@{}}{\displaystyle\eta^{-2}:\qquad\mathcal{B}_{1}{\bf{J}}_{0}=\bf{0}}\\\\[5.69054pt]
{\displaystyle\eta^{-1}:\qquad\mathcal{B}_{1}{\bf{J}}_{1}+\mathcal{B}_{2}{\bf{J}}_{0}=\bf{0}}\\\\[5.69054pt]
{\displaystyle\;\,\eta^{0}:\qquad\mathcal{B}_{1}{\bf{J}}_{2}+\mathcal{B}_{2}{\bf{J}}_{1}+\mathcal{B}_{3}{\bf{J}}_{0}={\rm
i}\omega{\bf{E}}_{0}.}\end{array}$ (41)

Since ${\bf x}$ and ${\bf y}$ are considered as independent variables,
equations (40)-(41) are PDEs in ${\bf y}$ while ${\bf x}$ plays the role of a
parameter.

Note that

$\begin{array}[]{lll}{\displaystyle(\mathcal{A}_{1}{\bf
E}_{0},{\bf{E}}_{0})_{Y}=(\mu^{-1}({\bf y})\mathbf{curl}_{\bf
y}\,{\bf{E}}_{0},\mathbf{curl}_{\bf y}\,{\bf{E}}_{0})_{Y},}\\\\[5.69054pt]
{\displaystyle(\mathcal{B}_{1}{\bf J}_{0},{\bf{J}}_{0})_{Y}=(\beta^{\ast}({\bf
y})\mathbf{div}_{\bf y}\,{\bf{J}}_{0},\mathbf{div}_{\bf
y}\,{\bf{J}}_{0})_{Y}.}\end{array}$

Consequently, the first equations of (40) and (41) are equivalent to

$\begin{array}[]{@{}l@{}}{\displaystyle\mathbf{curl}_{\bf
y}\,{\bf{E}}_{0}={\bf 0},\quad\mathbf{div}_{\bf y}\,{\bf
J}_{0}}=\bf{0},\end{array}$ (42)

respectively. Using (42), we can reduce the second equations of (40) and (41)
to

$\begin{array}[]{@{}l@{}}{\displaystyle\mathbf{curl}_{\bf y}\,(\mu^{-1}({\bf
y})\,\mathbf{curl}_{\bf y}\,{\bf{E}}_{1}+\mu^{-1}({\bf y})\mathbf{curl}_{\bf
x}\,{\bf{E}}_{0})={\bf 0},}\\\\[5.69054pt] {\displaystyle\mathbf{grad}_{\bf
y}\,(\beta^{\ast}({\bf y})\,\mathbf{div}_{\bf
y}\,{\bf{J}}_{1}+\beta^{\ast}({\bf
y})\mathbf{div}_{\bf{x}}\,{\bf{J}}_{0})=\bf{0}.}\end{array}$ (43)

Applying $\mathbf{div}_{\bf y}$ and $\mathbf{curl}_{\bf y}$ to the third
equations of (40) and (41), respectively, we have

$\begin{array}[]{@{}l@{}}{\displaystyle\mathbf{div}_{\bf
y}\,\mathbf{curl}_{\bf x}(\mu^{-1}({\bf y})\,\mathbf{curl}_{\bf
y}{\bf{E}}_{1}+\mu^{-1}({\bf y})\,\mathbf{curl}_{\bf
x}{\bf{E}}_{0})-\omega^{2}\mathbf{div}_{\bf y}(\varepsilon({\bf
y}){\bf{E}}_{0})}\\\\[5.69054pt] {\displaystyle\qquad={\rm
i}\omega\,\mathbf{div}_{\bf y}{\bf{J}}_{0}={\bf 0},}\\\\[5.69054pt]
{\displaystyle\mathbf{curl}_{\bf y}\,\mathbf{grad}_{\bf x}(\beta^{\ast}({\bf
y})\,\mathbf{div}_{\bf y}{\bf{J}}_{1}+\beta^{\ast}({\bf y})\,\mathbf{div}_{\bf
x}{\bf{J}}_{0})+\mathbf{curl}_{\bf y}\big{(}\gamma^{\ast}({\bf
y})\,{\bf{J}}_{0}\big{)}}\\\\[5.69054pt] {\displaystyle\qquad={\rm
i}\omega\,\mathbf{curl}_{\bf y}{\bf{E}}_{0}={\bf 0},}\end{array}$ (44)

where we have used (42). Since $\mathbf{div}_{\bf y}\,\mathbf{curl}_{\bf
x}=-\mathbf{div}_{\bf{x}}\,\mathbf{curl}_{\bf y}$ and $\mathbf{curl}_{\bf
y}\,\mathbf{grad}_{\bf x}=-\mathbf{curl}_{\bf x}\,\mathbf{grad}_{\bf y}$,
combining (43) and (44), we obtain

$\begin{array}[]{@{}l@{}}\mathbf{div}_{\bf y}(\varepsilon({\bf
y}){\bf{E}}_{0})=\mathbf{curl}_{\bf y}\big{(}\gamma^{\ast}({\bf
y})\,{\bf{J}}_{0}\big{)}={\bf 0}.\end{array}$ (45)

Denote by

$\begin{array}[]{@{}l@{}}{\displaystyle{\bf{w}}_{1}=\mu^{-1}({\bf
y})\,\mathbf{curl}_{\bf y}\,{\bf{E}}_{1}+\mu^{-1}({\bf y})\mathbf{curl}_{\bf
x}\,{\bf{E}}_{0},}\\\\[5.69054pt] {\displaystyle{w}_{2}=\beta^{\ast}({\bf
y})\,\mathbf{div}_{\bf y}\,{\bf{J}}_{1}+\beta^{\ast}({\bf
y})\mathbf{div}_{\bf{x}}\,{\bf{J}}_{0}},\end{array}$ (46)

and

$\begin{array}[]{@{}l@{}}\mathcal{M}_{Y}({\bf{w}}_{1})=\tilde{{\bf{w}}}_{1},\,\quad\mathcal{M}_{Y}({{w}}_{2})=\tilde{{{w}}}_{2},\\\\[5.69054pt]
\mathcal{M}_{Y}({\bf{E}}_{0})=\tilde{{\bf{E}}}_{0},\,\quad\mathcal{M}_{Y}({\bf{J}}_{0})=\tilde{{\bf{J}}}_{0},\end{array}$
(47)

where $\mathcal{M}_{Y}(f)=\frac{1}{|Y|}\int_{Y}f({\bf y})\,d{\bf y}$ denotes
the mean value of $f$ over the reference cell $Y$.

It follows from (42)-(43) and the definition of ${\bf{w}}_{1},{{w}}_{2}$ that

$\begin{array}[]{@{}l@{}}\mathbf{curl}_{\bf y}\,{\bf{w}}_{1}={\bf
0},\,\quad\mathbf{grad}_{\bf y}\,{w}_{2}={\bf 0},\\\\[5.69054pt]
\mathbf{div}_{\bf y}\,(\mu({\bf y}){\bf{w}}_{1})=\mathbf{div}_{\bf
y}\,\mathbf{curl}_{\bf
x}\,{\bf{E}}_{0}=-\mathbf{div}_{\bf{x}}\,\mathbf{curl}_{\bf
y}\,{\bf{E}}_{0}={\bf 0}.\end{array}$ (48)

Therefore, we see that $\mathbf{curl}_{\bf y}\,({\bf{w}}_{1}-\tilde{{\bf
w}}_{1})={\bf 0}$ and $\mathcal{M}_{Y}({\bf{w}}_{1}-\tilde{\bf
w}_{1})={\bm{0}}$, which imply that there exists a $Y$-periodic function
${\psi}({\bf x},{\bf y})$ such that

$\begin{array}[]{@{}l@{}}{\bf{w}}_{1}-\tilde{{\bf w}}_{1}=\mathbf{grad}_{\bf
y}\,{\psi}.\end{array}$ (49)

Similarly, in view of (42), there exists $Y$-periodic functions ${\phi}({\bf
x},{\bf y})$ and ${\bm{\varphi}}({\bf x},{\bf y})$ such that

$\begin{array}[]{@{}l@{}}{\bf{E}}_{0}-\tilde{{\bf E}}_{0}=\mathbf{grad}_{\bf
y}\,{\phi},\quad{\bf{J}}_{0}-\tilde{{\bf J}}_{0}=\mathbf{curl}_{\bf
y}\,{\bm{\varphi}}.\end{array}$ (50)

In addition, since ${\bf grad}_{\bf y}\,w_{2}={\bf 0}$, we see that
$w_{2}=\tilde{w}_{2}$. Substituting (49) into the last equation of (48) gives

$\begin{array}[]{@{}l@{}}{\displaystyle\mathbf{div}_{\bf y}\,\big{(}\mu({\bf
y})\,\mathbf{grad}_{\bf y}\,{\psi}\big{)}=-\mathbf{div}_{\bf y}\,(\mu({\bf
y})\tilde{{\bf w}}_{1})=-\sum_{i=1}^{3}(\tilde{{\bf
w}}_{1})_{i}\,\mathbf{div}_{\bf y}\,(\mu({\bf y}){\bf e}_{i})},\end{array}$
(51)

where $(\tilde{{\bf w}}_{1})_{i}$ denotes the $i$-th component of $\tilde{{\bf
w}}_{1}$ and $\\{{\bf e}_{1},{\bf e}_{2},{\bf
e}_{3}\\}=\\{(1,0,0)^{T},(0,1,0)^{T},(0,0,1)^{T}\\}$ is the canonical basis of
$\mathbb{R}^{3}$. Similarly, we deduce from (50) and (45) that

$\begin{array}[]{@{}l@{}}{\displaystyle\mathbf{div}_{\bf
y}\,\big{(}\varepsilon({\bf y})\,\mathbf{grad}_{\bf
y}\,{\phi}\big{)}=-\mathbf{div}_{\bf y}\,(\varepsilon({\bf y})\tilde{{\bf
E}}_{0})=-\sum_{i=1}^{3}(\tilde{{\bf E}}_{0})_{i}\,\mathbf{div}_{\bf
y}\,(\varepsilon({\bf y}){\bf e}_{i}),}\\\\[11.38109pt]
{\displaystyle\mathbf{curl}_{\bf y}\big{(}\gamma^{\ast}({\bf
y})\mathbf{curl}_{\bf y}\,{\bm{\varphi}}\big{)}=-\mathbf{curl}_{\bf
y}\,\big{(}\gamma^{\ast}({\bf y})\tilde{{\bf
J}}_{0}\big{)}=-\sum_{i=1}^{3}(\tilde{{\bf J}}_{0})_{i}\,\mathbf{curl}_{\bf
y}\,\big{(}\gamma^{\ast}({\bf y})\mathbf{e}_{i}\big{)}.}\end{array}$ (52)

To proceed further, we introduce the scalar-valued cell functions
$\theta_{i}^{\mu}({\bf y})$, $\theta_{i}^{\varepsilon}({\bf y})$, and the
vector-valued cell functions $\Theta_{i}^{\gamma}({\bf y})$, $i=1,2,3$, which
are the solutions of the following problems:

$\left\\{\begin{array}[]{@{}l@{}}{\displaystyle\mathbf{div}_{\bf y}\,(\mu({\bf
y})\,\mathbf{grad}_{\bf y}\,\theta_{i}^{\mu})=-\mathbf{div}_{\bf y}\,(\mu({\bf
y}){\bf e}_{i}),\quad{\rm in}\;\;Y}\\\\[5.69054pt]
{\displaystyle\theta_{i}^{\mu}\quad Y\,{\rm-}\,{\rm
periodic},}\end{array}\right.$ (53)

$\left\\{\begin{array}[]{@{}l@{}}{\displaystyle\mathbf{div}_{\bf
y}\,(\varepsilon({\bf y})\,\mathbf{grad}_{\bf
y}\,\theta_{i}^{\varepsilon})=-\mathbf{div}_{\bf y}\,(\varepsilon({\bf y}){\bf
e}_{i}),\quad{\rm in}\;\;Y}\\\\[5.69054pt]
{\displaystyle\theta_{i}^{\varepsilon}\quad Y\,{\rm-}\,{\rm
periodic},}\end{array}\right.$ (54)

and

$\left\\{\begin{array}[]{@{}l@{}}{\displaystyle\mathbf{curl}_{\bf
y}\,(\gamma^{\ast}({\bf y})\,\mathbf{curl}_{\bf
y}\,\Theta_{i}^{\gamma})=-\mathbf{curl}_{\bf y}\,(\gamma^{\ast}({\bf y}){\bf
e}_{i}),\quad{\rm in}\;\;Y}\\\\[5.69054pt] {\displaystyle\mathbf{div}_{\bf
y}\,\Theta_{i}^{\gamma}=0,\quad{\rm in}\;\;Y}\\\\[5.69054pt]
{\displaystyle\Theta_{i}^{\gamma}\quad Y\,{\rm-}\,{\rm
periodic},}\end{array}\right.$ (55)

respectively.

Next we define the vectors
$\bm{\theta}^{\mu}=(\theta_{1}^{\mu},\theta_{2}^{\mu},\theta_{3}^{\mu})$,
$\bm{\theta}^{\eta}=(\theta_{1}^{\eta},\theta_{2}^{\eta},\theta_{3}^{\eta})$,
and the matrix
$\bm{\Theta}^{\gamma}=\\{\Theta_{1}^{\gamma},\Theta_{2}^{\gamma},\Theta_{3}^{\gamma}\\}$.
Applying the superposition principle to the equations (51) and (52), we can
write $\psi$, $\phi$, and ${\bm{\varphi}}$ as

$\begin{array}[]{@{}l@{}}{\displaystyle\psi=\sum_{i=1}^{3}(\tilde{{\bf
w}}_{1})_{i}\,\theta_{i}^{\mu}=\bm{\theta}^{\mu}\cdot\tilde{{\bf
w}}_{1},\quad\phi=\sum_{i=1}^{3}(\tilde{{\bf
E}}_{0})_{i}\,\theta_{i}^{\mu}=\bm{\theta}^{\varepsilon}\cdot\tilde{{\bf
E}}_{0}},\\\\[5.69054pt]
{\displaystyle\qquad\qquad\quad{\bm{\varphi}}=\sum_{i=1}^{3}(\tilde{{\bf
J}}_{0})_{i}\,\Theta_{i}^{\gamma}=\bm{\Theta}^{\gamma}\cdot\tilde{{\bf
J}}_{0},}\end{array}$ (56)

respectively. Combining (49), (50), and (56), we obtain

$\begin{array}[]{@{}l@{}}{\bf{w}}_{1}=(I+\mathbf{grad}_{\bf
y}\,{\bm{\theta}^{\mu}})\tilde{{\bf
w}}_{1},\quad{\bf{E}}_{0}=(I+\mathbf{grad}_{\bf
y}\,{\bm{\theta}}^{\varepsilon})\tilde{{\bf E}}_{0},\\\\[8.53581pt]
{\bf\qquad\qquad\qquad{J}}_{0}=(I+\mathbf{curl}_{\bf
y}\,{\bm{\Theta}^{\gamma}})\tilde{{\bf J}}_{0}.\end{array}$ (57)

Integrating the two equations in (46) over the reference cell $Y$ after
multiplying them by $\mu({\bf y})$ and $1/\beta^{\ast}({\bf y})$,
respectively, it follows that

$\begin{array}[]{@{}l@{}}\mathbf{curl}_{\bf x}\,\tilde{{\bf
E}}_{0}=\mathcal{M}_{Y}(\mu({\bf
y}){\bf{w}}_{1})=\mathcal{M}_{Y}\big{(}\mu({\bf y})(I+\mathbf{grad}_{\bf
y}\,{\bm{\theta}}^{\mu})\big{)}\tilde{{\bf w}}_{1},\\\\[5.69054pt]
\mathbf{div}_{\bf{x}}\,\tilde{{\bf
J}}_{0}=\mathcal{M}_{Y}(\frac{1}{\beta^{\ast}({\bf
y})}w_{2})=\mathcal{M}_{Y}(\frac{1}{\beta^{\ast}({\bf
y})})\tilde{w}_{2},\end{array}$ (58)

where we have used (57) and the fact that $w_{2}=\tilde{w}_{2}$. Similarly,
applying $\mathcal{M}_{Y}$ to the third equations of (40) and (41),
respectively, recalling the definition of ${\bf w}_{1}$ and $w_{2}$, and using
(57), we arrive at

$\begin{array}[]{@{}l@{}}\mathbf{curl}_{\bf
x}\,\tilde{\bf{w}}_{1}-\omega^{2}\mathcal{M}_{Y}\big{(}\varepsilon({\bf
y})(I+\mathbf{grad}\,{\bm{\theta}^{\eta}})\big{)}\tilde{{\bf E}}_{0}={\rm
i}\omega\tilde{{\bf J}}_{0},\\\\[5.69054pt]
\mathbf{grad}_{\bf{x}}\,\tilde{{w}}_{2}+\mathcal{M}_{Y}\big{(}\gamma^{\ast}({\bf
y})(I+\mathbf{curl}_{\bf y}\,{\bm{\Theta}^{\gamma}})\big{)}\tilde{{\bf
J}}_{0}={\rm i}\omega\tilde{{\bf E}}_{0}.\end{array}$ (59)

Substituting (58) into (59), we obtain the homogenized system associated to
(34)

$\left\\{\begin{array}[]{@{}l@{}}\mathbf{curl}_{\bf
x}\,\big{(}\widehat{\mu}^{-1}\,\mathbf{curl}_{\bf x}\,\tilde{{\bf
E}}_{0}\big{)}-\omega^{2}\widehat{\varepsilon}\,\tilde{{\bf E}}_{0}={\rm
i}\omega\tilde{{\bf J}}_{0},\\\\[5.69054pt]
\mathbf{grad}_{\bf{x}}\,\big{(}\widehat{\beta^{\ast}}\,\mathbf{div}_{\bf{x}}\,\tilde{{\bf
J}}_{0}\big{)}+\widehat{\gamma^{\ast}}\,\tilde{{\bf J}}_{0}={\rm
i}\omega\tilde{{\bf E}}_{0},\end{array}\right.$ (60)

where the homogenized coefficients are given by

$\begin{array}[]{@{}l@{}}\widehat{\mu}=\mathcal{M}_{Y}(\mu({\bf
y})(I+\mathbf{grad}_{\bf
y}\,{\bm{\theta}}^{\mu})),\quad\widehat{\varepsilon}=\mathcal{M}_{Y}(\varepsilon({\bf
y})(I+\mathbf{grad}_{\bf y}\,{\bm{\theta}}^{\varepsilon})),\\\\[8.53581pt]
\widehat{\gamma^{\ast}}=\mathcal{M}_{Y}(\gamma^{\ast}({\bf
y})(I+\mathbf{curl}_{\bf
y}\,{\bm{\Theta}^{\gamma}})),\quad\widehat{\beta^{\ast}}=\mathcal{M}_{Y}(\frac{1}{\beta^{\ast}({\bf
y})})^{-1},\end{array}$ (61)

respectively, and the cell functions
$\bm{\theta}^{\mu}=(\theta_{1}^{\mu},\theta_{2}^{\mu},\theta_{3}^{\mu})$,
$\bm{\theta}^{\eta}=(\theta_{1}^{\eta},\theta_{2}^{\eta},\theta_{3}^{\eta})$,
and
$\bm{\Theta}^{\gamma}=\\{\Theta_{1}^{\gamma},\Theta_{2}^{\gamma},\Theta_{3}^{\gamma}\\}$
are defined in (53)-(55), respectively. Note that although $\mu({\bf y})$,
$\varepsilon({\bf y})$, and $\gamma^{\ast}({\bf y})$ are scalar-valued
functions, in general, the homogenized coefficients $\widehat{\mu}$,
$\widehat{\varepsilon}$, and $\widehat{\gamma^{\ast}}$ are non-diagonal
matrices.

In addition, we define the multiscale approximate solution for
$({\bf{E}}_{\eta},\bf J_{\eta})$ as follows:

$\left\\{\begin{array}[]{@{}l@{}}\tilde{\bf{E}}_{0,\eta}=(I+\mathbf{grad}_{\bf
y}\,{\bm{\theta}^{\varepsilon}})\,\tilde{{\bf E}}_{0},\\\\[5.69054pt]
\tilde{\bf{J}}_{0,\eta}=(I+\mathbf{curl}_{\bf
y}\,{\bm{\Theta}^{\gamma}})\,\tilde{{\bf J}}_{0}.\end{array}\right.$ (62)

It is worth pointing out that the homogenized system and the multiscale
approximate solution are only defined in the domain of the scatter. The
homogenized system (60) is coupled to Maxwell’s equations defined in the
exterior domain with constant coefficients and forms the homogenized coupled
system

$\left\\{\begin{array}[]{@{}l@{}}{\displaystyle\mathbf{curl}\,\big{(}\widehat{\mu}^{-1}\,\mathbf{curl}\,\tilde{{\bf
E}}_{0}\big{)}-\omega^{2}\widehat{\varepsilon}\,\tilde{{\bf E}}_{0}-{\rm
i}\omega\tilde{{\bf J}}_{0}={\bf 0},\;\quad\;{\rm
in}\;\;\Omega_{s},}\\\\[5.69054pt]
{\displaystyle\mathbf{grad}\,\big{(}\widehat{\beta^{\ast}}\,\mathbf{div}\,\tilde{{\bf
J}}_{0}\big{)}+\widehat{\gamma^{\ast}}\,\tilde{{\bf J}}_{0}-{\rm
i}\omega\tilde{{\bf E}}_{0}={\bf 0},\qquad\quad{\rm
in}\;\;\Omega_{s},}\\\\[5.69054pt]
{\displaystyle\mathbf{curl}\,(\mu_{0}^{-1}\,\mathbf{curl}\,\mathbf{E})-\varepsilon_{0}\omega^{2}\,\mathbf{E}=\mathbf{0},\qquad\qquad\;\;\quad{\rm
in}\;\;\Omega/\Omega_{s},}\\\\[5.69054pt]
{\displaystyle\mathbf{n}\times(\widehat{\mu}^{-1}\mathbf{curl}\,\tilde{\mathbf{E}}_{0})=\mathbf{n}\times(\mu_{0}^{-1}\mathbf{curl}\,\mathbf{E}),\,\qquad{\rm
on}\;\;\partial\Omega_{s},}\\\\[5.69054pt]
{\displaystyle\mathbf{n}\times\tilde{\mathbf{E}}_{0}=\mathbf{n}\times\mathbf{E},\qquad\qquad\quad\qquad\qquad\;\qquad{\rm
on}\;\;\partial\Omega_{s},}\\\\[5.69054pt] {\displaystyle{\bf
n}\cdot\tilde{\bf
J}_{0}=0,\qquad\qquad\;\;\,\qquad\qquad\qquad\;\qquad\quad{\rm
on}\;\;\partial\Omega_{s},}\\\\[5.69054pt]
{\displaystyle(\mu_{0}^{-1}\,\mathbf{curl}\,{\bf E})\times{\bf
n}-{\mathrm{{i}}}\omega({\bf n}\times{\bf E})\times{\bf n}={\bf
g},\;\;\quad{\rm on}\;\;\partial\Omega.}\end{array}\right.$ (63)

######  Remark 3.1.

The coefficients of the homogenized system are constant in the domain of the
scatter, making the numerical solution of the homogenized problem much easier
than the original problem. Physically, the homogenization procedure is
equivalent to replacing the heterogeneous scatter comprising of metallic
nanostructures and the dielectric medium by an averaging homogeneous scatter
as illustrated in Fig 3.1.

For simplicity, we set

$\displaystyle\mu_{\rm
eff}=\left\\{\begin{array}[]{l}\widehat{\mu},\,\;\;\quad\rm{in}\;\;\Omega_{s}\\\\[8.53581pt]
\mu_{0},\,\,\quad\rm{in}\;\;\Omega/\Omega_{s},\end{array}\right.\quad\varepsilon_{\rm
eff}=\left\\{\begin{array}[]{l}\widehat{\varepsilon},\,\;\,\;\quad\rm{in}\;\;\Omega_{s}\\\\[8.53581pt]
\varepsilon_{0},\,\;\,\quad\rm{in}\;\;\Omega/\Omega_{s},\end{array}\right.$
(68)

and denote by $\tilde{{\mathbf{E}}}_{0}$ the electric field on the whole
domain $\Omega$. Then, the homogenized coupled system (63) can be rewritten as

$\left\\{\begin{array}[]{@{}l@{}}{\displaystyle\mathbf{curl}\,\big{(}{\mu}_{\rm
eff}^{-1}\,\mathbf{curl}\,\tilde{{\bf
E}}_{0}\big{)}-\omega^{2}{\varepsilon}_{\rm eff}\,\tilde{{\bf E}}_{0}-{\rm
i}\omega\tilde{{\bf J}}_{0}={\bf 0},\;\;\quad\;{\rm
in}\;\;\Omega,}\\\\[8.53581pt]
{\displaystyle\mathbf{grad}\,\big{(}\widehat{\beta^{\ast}}\,\mathbf{div}\,\tilde{{\bf
J}}_{0}\big{)}+\widehat{\gamma^{\ast}}\,\tilde{{\bf J}}_{0}-{\rm
i}\omega\tilde{{\bf E}}_{0}={\bf 0},\,\quad\qquad\quad{\rm
in}\;\;\Omega_{s},}\\\\[5.69054pt] {\displaystyle{\bf n}\cdot\tilde{\bf
J}_{0}=0,\qquad\qquad\quad\;\,\qquad\qquad\qquad\;\qquad\quad{\rm
on}\;\;\partial\Omega_{s},}\\\\[5.69054pt] {\displaystyle(\mu_{\rm
eff}^{-1}\,\mathbf{curl}\,\tilde{\bf E}_{0})\times{\bf
n}-{\mathrm{{i}}}\omega({\bf n}\times\tilde{\bf E}_{0})\times{\bf n}={\bf
g},\;\;\quad{\rm on}\;\;\partial\Omega.}\end{array}\right.$ (69)

Plane waveexcitationHomogenizationPlane waveexcitation Fig. 2: Left: the
original problem; right: the homogenized problem.

##  4 Multiscale approaches and the associated numerical algorithms

Based on the discussion in the previous sections, we have the multiscale
approach for solving the problem (21) consisting of the following steps:

Step 1. Choose a sufficiently large parameter $\lambda>0$, set the
coefficients $\gamma_{\eta,\lambda}$ and $\beta^{2}_{\eta,\lambda}$ according
to (31), and extend the equation satisfied by the polarization current into
the domain occupied by the dielectric medium.

Step 2. Solve the equations (53)-(55) in the reference cell $Y$ to get the
cell functions $\theta_{i}^{\mu}({\bf y})$, $\theta_{i}^{\varepsilon}({\bf
y})$, $\Theta_{i}^{\gamma}({\bf y})$, $i=1,2,3$, and compute the homogenized
coefficients $\widehat{\mu}$, $\widehat{\varepsilon}$,
$\widehat{\gamma^{\ast}}$, and $\widehat{\beta^{\ast}}$ given by (61).

Step 3. Solve the homogenized coupled system (69) in the domain $\Omega$ to
get the homogenized solution ($\tilde{\bf E}_{0}$, $\tilde{\bf J}_{0}$).

Step 4. Compute the multiscale approximate solution (62) by adding some
correctors to the homogenized solution.

The proposed multiscale approach reduces the direct solution of the
challenging problem (21) to the solution of several cell problems and a
homogenized problem with constant coefficients, which is capable of saving
much computational cost especially when there exists a huge number of metallic
nanostructures (the parameter $\eta$ is very small).

Numerical results in section 5 show that the electric field calculated by the
multiscale approach is in good agreement with the reference solution
(numerical solution of the original system on a very fine mesh) outside the
metallic nanostructures. However, it is much less accurate inside the metallic
nanostructures. The main reason for the failure of the multiscale approach
inside the metallic nanostructures is that due to the large extension
parameter $\lambda$, the imaginary part of the homogenized coefficient
$\widehat{\gamma^{\ast}}$ is also very large, making the polarization current
rather small and thus the electric field loses the nonlocal information. In
fact, we have

######  Lemma 2.

Assume that there exists a positive constant $\alpha$ such that the
homogenized coefficients $\widehat{\gamma^{\ast}}$ satisfies

${\rm
Im}(\widehat{\gamma^{\ast}}){\bm{\xi}}\cdot{\bm{\xi}}\geq\alpha|{\bm{\xi}}|^{2},\quad\forall{\bm{\xi}}\in\mathbb{R}^{3}.$
(70)

Then the solution of the homogenized coupled system (69) satisfies

$\|\tilde{\bf J}_{0}\|_{\mathbf{L}^{2}(\Omega_{s})}\leq\frac{C}{\alpha},$ (71)

where $C$ is a positive constant independent of $\alpha$.

The proof of this lemma is very similar to that of Lemma 5 and thus we omit it
here.

######  Lemma 3.

Let $(\tilde{\bf E}_{0},\tilde{\bf J}_{0})$ be the solution of (69) and
$\mathbf{E}$ be the solution of the Maxwell’s equations

$\left\\{\begin{array}[]{@{}l@{}}{\displaystyle\mathbf{curl}\,\big{(}{\mu}_{\rm
eff}^{-1}\,\mathbf{curl}\,{{\bf E}}\big{)}-\omega^{2}{\varepsilon}_{\rm
eff}\,{{\bf E}}={\bf 0},\;\quad\;{\rm in}\;\;\Omega,}\\\\[8.53581pt]
{\displaystyle(\mu_{\rm eff}^{-1}\,\mathbf{curl}\,{\bf E})\times{\bf
n}-{\mathrm{{i}}}\omega({\bf n}\times{\bf E})\times{\bf n}={\bf
g},\;\;\quad{\rm on}\;\;\partial\Omega.}\end{array}\right.$ (72)

Under the assumption of Lemma 2, we have the following estimate

$\|\tilde{\mathbf{E}}_{0}-\mathbf{E}\|_{\mathbf{H}(\mathbf{curl};\Omega)}\leq\frac{C}{\alpha}.$
(73)

###### Proof.

Let $\mathbf{e}=\tilde{\mathbf{E}}_{0}-\mathbf{E}$. By subtracting the first
equation of (69) from (72), we see that $\mathbf{e}$ satisfies

$\left\\{\begin{array}[]{@{}l@{}}{\displaystyle\mathbf{curl}\,\big{(}{\mu}_{\rm
eff}^{-1}\,\mathbf{curl}\,{{\bf e}}\big{)}-\omega^{2}{\varepsilon}_{\rm
eff}\,{{\bf e}}={\rm i}\omega\tilde{\bf J}_{0},\;\quad\;{\rm
in}\;\;\Omega,}\\\\[8.53581pt] {\displaystyle(\mu_{\rm
eff}^{-1}\,\mathbf{curl}\,{\bf e})\times{\bf n}-{\mathrm{{i}}}\omega({\bf
n}\times{\bf e})\times{\bf n}={\bf 0},\;\;\quad{\rm
on}\;\;\partial\Omega.}\end{array}\right.$ (74)

From Theorem 4.17 of [15], we get the estimate

$\|\mathbf{e}\|_{\mathbf{H}(\mathbf{curl};\Omega)}\leq C\|\tilde{\bf
J}_{0}\|_{\mathbf{L}^{2}(\Omega_{s})},$ (75)

which yields the desired estimate (73) by using (71). ∎

######  Remark 4.1.

The above two lemmas show that as the imaginary part of the homogenized
coefficient $\widehat{\gamma^{\ast}}$ becomes larger, the influence of the
polarization current on the electric field becomes smaller. If $\alpha$ tends
to infinity, the coupled system (69) will reduce to the pure Maxwell’s
equations and the electric field will lose the nonlocal information inside the
metallic nanostructures. Although we are not able to give an estimate of the
lower bound $\alpha$ in terms of the extension parameter $\lambda$, we have
tested some practical examples in section 5 and verified numerically that in
these examples, $\alpha$ tends to infinity as $\lambda\rightarrow\infty$.

Motivated by the numerical results of the previous multiscale approach and
Lemma 2 and 3, we propose the following modified multiscale approach.

Modified Multiscale Approach:

Step 1. Solve the equations (53)-(54) in the reference cell $Y$ to get the
cell functions $\theta_{i}^{\mu}({\bf y})$, $\theta_{i}^{\varepsilon}({\bf
y})$, and compute the homogenized coefficients $\widehat{\mu}$ and
$\widehat{\varepsilon}$ given by (61).

Step 2. Solve the following homogenized Maxwell’s equation on the whole domain

$\left\\{\begin{array}[]{@{}l@{}}{\displaystyle\mathbf{curl}\,\big{(}{\mu}_{\rm
eff}^{-1}\,\mathbf{curl}\,{{\bf E}}^{0}\big{)}-\omega^{2}{\varepsilon}_{\rm
eff}\,{{\bf E}}^{0}={\bf 0},\;\quad\;{\rm in}\;\;\Omega,}\\\\[8.53581pt]
{\displaystyle(\mu_{\rm eff}^{-1}\,\mathbf{curl}\,{\bf E}^{0})\times{\bf
n}-{\mathrm{{i}}}\omega({\bf n}\times{\bf E}^{0})\times{\bf n}={\bf
g},\;\;\quad{\rm on}\;\;\partial\Omega.}\end{array}\right.$ (76)

Step 3. Compute the multiscale approximate solution
${\bf{E}}^{0}_{\eta}=(I+\mathbf{grad}_{\bf
y}\,{\bm{\theta}^{\varepsilon}})\,{{\bf E}}^{0}$ and then find the modified
multiscale approximate solutions ${\bf{E}}^{M}_{\eta}$ for the electric field
such that ${\bf{E}}^{M}_{\eta}={\bf{E}}^{0}_{\eta}$ outside the metallic
nanostructures and ${\bf{E}}^{M}_{\eta}$ satisfies

$\left\\{\begin{array}[]{@{}l@{}}{\displaystyle\mathbf{curl}\,(\mu_{k}^{-1}\mathbf{curl}\,{\bf
E}^{M}_{\eta})-\varepsilon_{k}\,\omega^{2}\,{\bf
E}^{M}_{\eta}-{\mathrm{{i}}}\omega{\bf J}^{M}_{\eta}={\bf 0},\;\qquad\quad{\rm
in}\;\;\Omega^{k}_{\eta},}\\\\[5.69054pt]
{\displaystyle\frac{\omega(\omega+{\mathrm{{i}}}\gamma)}{\omega^{2}_{p}\varepsilon_{0}}{\bf
J}^{M}_{\eta}+\frac{\beta^{2}}{\omega^{2}_{p}\varepsilon_{0}}\mathbf{grad}\,(\mathbf{div}\,{\bf
J}^{M}_{\eta})-{\mathrm{{i}}}\omega{\bf E}^{M}_{\eta}={\bf 0},\quad{\rm
in}\;\;\Omega^{k}_{\eta},}\\\\[11.38109pt] {\displaystyle{\bf n}\cdot{\bf
J}^{M}_{\eta}=0,\;\;\;\;\qquad\qquad{\rm
on}\;\;\partial\Omega^{k}_{\eta},}\\\\[5.69054pt] {\displaystyle{\bf
n}\times{\bf E}^{M}_{\eta}={\bf n}\times{\bf E}_{\eta}^{0},\quad\quad{\rm
on}\;\;\partial\Omega^{k}_{\eta},}\end{array}\right.$ (77)

for $k=1,\cdots,N$, where $\Omega_{\eta}^{k}$ is the domain occupied by the
$k$-th metallic nanostructure, $N$ is the number of the metallic
nanostructures, $\mu_{k}=\mu_{\eta}|_{\Omega_{\eta}^{k}}$, and
$\varepsilon_{k}=\varepsilon_{\eta}|_{\Omega_{\eta}^{k}}$.

There are two main differences between the modified multiscale approach and
the previous one. First, in our modified approach, we only solve the
homogenized Maxwell’s equation (76) instead of the homogenized coupled system
(69) on the whole domain. This is based on the observation that the
homogenized coupled system (69) is close to the pure Maxwell’s equations if we
choose a sufficiently large extension parameter $\lambda$. By doing so, we can
reduce much computational cost since we don’t need to solve the cell equation
(55), either. Second, in our modified approach, we modify the multiscale
approximate solution ${\bf{E}}^{0}_{\eta}$ inside the metallic nanostructures
by solving the original coupled system with boundary conditions given by
${\bf{E}}^{0}_{\eta}$ in each metallic nanostructure. This is based on the
observation in numerical experiments that ${\bf{E}}^{0}_{\eta}$ computed by
the previous multiscale approach agrees well with the reference solution
outside the metallic nanostructures while substantially deviates from the
reference solution inside the metallic nanostructures. Note that (77) is
solved in each metallic nanostructure separately, where the coefficients
$\mu_{k}$, $\varepsilon_{k}$, $\gamma$ and $\beta$ are constant. Therefore,
the computational burden of solving (77) is much lower than solving the
coupled system on the whole domain. In addition, in view of the periodic
arrangement of metallic nanostructures within the dielectric medium, we have a
fast algorithm for solving (77) which is described in detail below.

Denoting by
$\mathbf{E}^{M,k}_{\eta}=\mathbf{E}^{M}_{\eta}|_{\Omega_{\eta}^{k}}$,
$\mathbf{J}^{M,k}_{\eta}=\mathbf{J}^{M}_{\eta}|_{\Omega_{\eta}^{k}}$, and
assuming that
$\mathbf{E}^{M,k}_{\eta}=\widetilde{\mathbf{E}}^{M,k}_{\eta}+{\bf{E}}^{0}_{\eta}|_{\Omega_{\eta}^{k}}$,
then the problem (77) is equivalent to finding
$\widetilde{\mathbf{E}}^{M,k}_{\eta}$ such that

$\left\\{\begin{array}[]{@{}l@{}}{\displaystyle\mathbf{curl}\,(\mu_{k}^{-1}\mathbf{curl}\,\widetilde{\bf
E}^{M,k}_{\eta})-\varepsilon_{k}\,\omega^{2}\,\widetilde{\bf
E}^{M,k}_{\eta}-{\mathrm{{i}}}\omega{\bf
J}^{M}_{\eta}=-\mathbf{curl}\,(\mu_{k}^{-1}\mathbf{curl}\,{\bf
E}^{0}_{\eta})+\varepsilon_{k}\,\omega^{2}\,{\bf E}^{0}_{\eta},\quad{\rm
in}\;\;\Omega^{k}_{\eta},}\\\\[5.69054pt]
{\displaystyle\frac{\omega(\omega+{\mathrm{{i}}}\gamma)}{\omega^{2}_{p}\varepsilon_{0}}{\bf
J}^{M,k}_{\eta}+\frac{\beta^{2}}{\omega^{2}_{p}\varepsilon_{0}}\mathbf{grad}\,(\mathbf{div}\,{\bf
J}^{M,k}_{\eta})-{\mathrm{{i}}}\omega\widetilde{\bf
E}^{M,k}_{\eta}={\mathrm{{i}}}\omega{\bf E}^{0}_{\eta},\quad{\rm
in}\;\;\Omega^{k}_{\eta},}\\\\[11.38109pt] {\displaystyle{\bf n}\cdot{\bf
J}^{M,k}_{\eta}=0,\;\;\;\qquad{\rm
on}\;\;\partial\Omega^{k}_{\eta},}\\\\[5.69054pt] {\displaystyle{\bf
n}\times\widetilde{\bf E}^{M,k}_{\eta}={\bf 0},\quad\quad{\rm
on}\;\;\partial\Omega^{k}_{\eta}.}\end{array}\right.$ (78)

Assume that
${\bf{E}}^{0}_{\eta}|_{\Omega_{\eta}^{k}}\in\mathbf{H}(\mathbf{curl};\Omega_{\eta}^{k})$
for $k=1,\cdots,N$. The variational formulation of (78) is: Find
$\widetilde{\mathbf{E}}^{M,k}_{\eta}\in\mathbf{H}_{0}(\mathbf{curl};\Omega_{\eta}^{k})$,
${\mathbf{J}}^{M,k}_{\eta}\in\mathbf{H}_{0}(\mathbf{div};\Omega_{\eta}^{k})$
such that

$\left\\{\begin{array}[]{@{}l@{}}{\displaystyle\mu_{k}^{-1}(\mathbf{curl}\,\widetilde{\bf
E}^{M,k}_{\eta},\,\mathbf{curl}\,\mathbf{u})-\varepsilon_{k}\,\omega^{2}(\widetilde{\bf
E}^{M,k}_{\eta},\,\mathbf{u})-{\mathrm{{i}}}\omega({\bf
J}^{M,k}_{\eta},\,\mathbf{u})}\\\\[5.69054pt]
{\displaystyle\quad=-\mu_{k}^{-1}(\mathbf{curl}\,{\bf
E}^{0}_{\eta},\,\mathbf{curl}\,\mathbf{u})+\varepsilon_{k}\,\omega^{2}({\bf
E}^{0}_{\eta},\,\mathbf{u}),\qquad\forall\mathbf{u}\in\mathbf{H}_{0}(\mathbf{curl};\Omega_{\eta}^{k})}\\\\[5.69054pt]
{\displaystyle\frac{\omega(\omega+{\mathrm{{i}}}\gamma)}{\omega^{2}_{p}\varepsilon_{0}}({\bf
J}^{M,k}_{\eta},\,\mathbf{w})-\frac{\beta^{2}}{\omega^{2}_{p}\varepsilon_{0}}(\mathbf{div}\,{\bf
J}^{M,k}_{\eta},\,\mathbf{div}\,\mathbf{w})-{\mathrm{{i}}}\omega(\widetilde{\bf
E}^{M,k}_{\eta},\,\mathbf{w})}\\\\[11.38109pt]
{\displaystyle\quad={\mathrm{{i}}}\omega({\bf
E}^{0}_{\eta},\,\mathbf{w}),\qquad\forall\mathbf{w}\in\mathbf{H}_{0}(\mathbf{div};\Omega_{\eta}^{k}).}\end{array}\right.$
(79)

Let $\mathcal{T}^{k}_{h}$ be a uniform triangulation of $\Omega_{\eta}^{k}$
into tetrahedrons of maximal diameter $h$. We define the lowest order Nédélec
$\mathbf{H}_{0}(\mathbf{curl};\Omega_{\eta}^{k})$-conforming and
Raviart–Thomas $\mathbf{H}_{0}(\mathbf{div};\Omega_{\eta}^{k})$-conforming
finite element spaces as follows [15]:

$\begin{array}[]{lll}{\displaystyle{X}^{k}_{h}=\\{{\bf
u}_{h}\in\mathbf{H}_{0}(\mathbf{curl};\Omega_{\eta}^{k}):\;{\bf
u}_{h}|_{T}=\mathbf{a}_{T}+\mathbf{b}_{T}\times\mathbf{x}\;\;{\rm
with}\;\mathbf{a}_{T},\mathbf{b}_{T}\in\mathbb{R}^{3},\;\;\forall
T\in\mathcal{T}^{k}_{h}\\},}\\\\[5.69054pt] {\displaystyle{Y}^{k}_{h}=\\{{\bf
u}_{h}\in\mathbf{H}_{0}(\mathbf{div};\Omega_{\eta}^{k}):\;{\bf
u}_{h}|_{T}={a}_{T}\mathbf{x}+\mathbf{b}_{T}\;\;{\rm
with}\;{a}_{T}\in\mathbb{R},\mathbf{b}_{T}\in\mathbb{R}^{3},\;\;\forall
T\in\mathcal{T}^{k}_{h}\\}}.\end{array}$

The finite element approximation of (79) is formulated as follows: find
$\widetilde{\mathbf{E}}^{M,k}_{\eta,h}\in{X}^{k}_{h}$,
${\mathbf{J}}^{M,k}_{\eta,h}\in{Y}^{k}_{h}$ such that

$\left\\{\begin{array}[]{@{}l@{}}{\displaystyle\mu_{k}^{-1}(\mathbf{curl}\,\widetilde{\bf
E}^{M,k}_{\eta,h},\,\mathbf{curl}\,\mathbf{u}_{h})-\varepsilon_{k}\,\omega^{2}(\widetilde{\bf
E}^{M,k}_{\eta,h},\,\mathbf{u}_{h})-{\mathrm{{i}}}\omega({\bf
J}^{M,k}_{\eta,h},\,\mathbf{u}_{h})}\\\\[5.69054pt]
{\displaystyle\quad=-\mu_{k}^{-1}(\mathbf{curl}\,{\bf
E}^{0}_{\eta},\,\mathbf{curl}\,\mathbf{u}_{h})+\varepsilon_{k}\,\omega^{2}({\bf
E}^{0}_{\eta},\,\mathbf{u}_{h}),\qquad\forall\mathbf{u}_{h}\in{X}^{k}_{h},}\\\\[5.69054pt]
{\displaystyle\frac{\omega(\omega+{\mathrm{{i}}}\gamma)}{\omega^{2}_{p}\varepsilon_{0}}({\bf
J}^{M,k}_{\eta,h},\,\mathbf{w}_{h})-\frac{\beta^{2}}{\omega^{2}_{p}\varepsilon_{0}}(\mathbf{div}\,{\bf
J}^{M,k}_{\eta,h},\,\mathbf{div}\,\mathbf{w}_{h})-{\mathrm{{i}}}\omega(\widetilde{\bf
E}^{M,k}_{\eta,h},\,\mathbf{w}_{h})}\\\\[11.38109pt]
{\displaystyle\quad={\mathrm{{i}}}\omega({\bf
E}^{0}_{\eta},\,\mathbf{w}_{h}),\qquad\forall\mathbf{w}_{h}\in{Y}^{k}_{h}.}\end{array}\right.$
(80)

The system (80) can be rewritten as the following algebraic equation

$\mathbb{A}_{k}\mathbf{x}_{k}=\mathbf{b}_{k},$ (81)

where $\mathbb{A}_{k}$ is the matrix, $\mathbf{x}_{k}$ and $\mathbf{b}_{k}$
are the vectors representing the unknowns and the right-hand term,
respectively.

If there is only one metallic nanostructure in the reference cell as shown in
Fig 2.1, then the domains $\Omega_{\eta}^{k}$ $(k=1,\cdots,N)$ are essentially
the same (translation-invariant). Therefore, we can generate the same mesh
$\mathcal{T}_{h}$ for all the domains $\Omega_{\eta}^{k}$, which results in
the same matrix $\mathbb{A}$ for the algebraic equation (81) due to the
translation-invariant property of the $\mathbf{curl}$ and $\mathbf{div}$
operators and the fact that the coefficients $\mu_{k}$ and $\varepsilon_{k}$
take the same value in all the domains, respectively. In this case, we need to
solve a set of linear equations

$\mathbb{A}\mathbf{x}_{k}=\mathbf{b}_{k},\qquad k=1,\cdots,N$ (82)

with the same matrix and many different right-hand vectors. To solve this
problem, we can perform an $LU$ decomposition of the matrix $\mathbb{A}$ once
and then solve the triangular matrices for the different $\mathbf{b}_{k}$.
This algorithm can be easily extended to the case that there exist several
metallic nanostructures in the reference cell. In this way, we can make the
computational burden of solving (77) controllable and insensitive to the
number of the metallic nanostructures.

Finally, we briefly discuss the numerical methods for solving the scalar cell
equations (53)-(54) and the Maxwell’s equations (76). Note that the cell
equations are the same as for the scalar elliptic operators
$-\mathbf{div}\,\mu_{\eta}\,\mathbf{grad}$ and
$-\mathbf{div}\,\varepsilon_{\eta}\,\mathbf{grad}$ and they can be solved
numerically by the adaptive finite element method as discussed in [25]. As for
the Maxwell’s equations (76), a Galerkin finite element method based on the
Nédélec elements can be used and the resulting linear system can be solved by
a multigrid preconditioned Minres algorithm. For more details, we refer the
reader to [2, 15].

######  Remark 4.2.

Although we don’t use the technique of extension and the homogenized equation
for the polarization current explicitly in the modified multiscale approach,
they motivate and justify the modified approach. In fact, as indicated in
Lemma 2 and 3, it is because the imaginary part of the homogenized coefficient
$\widehat{\gamma^{\ast}}$ is sufficiently large due to the large extension
parameter $\lambda$ that we can replace the homogenized coupled system (69) by
the homogenized Maxwell’s equations (76) in our modified approach. In
addition, the multiscale asymptotic method discussed in section 3 can be used
for the homogenization of composite materials in which the inclusions and the
matrix are both metallic nanostructures with different physical parameters.

##  5 Numerical examples

(a) (b)

Fig. 3: (a) periodic gold nanoparticle arrays embedded in a dielectric medium.
(b) The rescaled reference cell $Q$.

In this section, we present several numerical experiments to validate the
effectiveness of our method.

######  Example 5.1.

In this example, we consider the scattering of a plane wave by periodic gold
nanoparticle arrays embedded in a dielectric medium as shown in Fig. 5.1.
There are 4 nanospheres with a radius of 2 nm in each direction and the inter-
particle distance is 1 nm. The nanoparticle arrays are irradiated by a plane
wave $\mathbf{E}^{inc}={\rm exp}({\mathrm{{i}}}\omega y)\mathbf{e}_{x}$
propagating in the $y-$direction. The Silver–Müller absorbing boundary
condition is set on the boundary of a sphere of radius 40 nm to truncate the
computational domain. We use the same material parameters for the gold
nanoparticle as in [17], which are summarized in Table 1. For the dielectric
medium, we consider the following two cases:

######  Case 5.1.

Silicon dioxide ($SiO_{2}$):
$\quad\mu=\mu_{0},\quad\varepsilon=3.9\,\varepsilon_{0}$;

######  Case 5.2.

water ($H_{2}O$) : $\quad\mu=\mu_{0},\quad\varepsilon=80\,\varepsilon_{0}$.

Table 1: Physical parameters of the NHD model. $\epsilon_{0}$ and $\mu_{0}$ are the electric permittivity and magnetic permeability of free space, respectively. $\omega_{p}$ | $\gamma$ | $\beta$ | $\epsilon$ | $\mu$
---|---|---|---|---
$1.37\times 10^{16}$ rad/s |  $1.08\times 10^{14}$ rad/s |  $1.08\times 10^{6}$ m/s | $9.5\,\epsilon_{0}$ | $\mu_{0}$

First, we solve the original system (21) and the extended system (32) directly
on a very fine mesh and test the error between their solutions. Set

$err=\|{\bf E}_{\eta}-{\bf
E}_{\eta,\lambda}\|_{\mathbf{H}_{T}(\mathbf{curl},\Omega)}+\|{\bf
J}_{\eta,\lambda}-{\bf J}_{\eta}\|_{\mathbf{H}(\mathbf{div},\Omega_{s})}$

where $({\bf E}_{\eta},{\bf J}_{\eta})$ and $({\bf E}_{\eta,\lambda},{\bf
J}_{\eta,\lambda})$ are the numerical solutions of the original system and the
extended system, respectively. In Fig. 5.2, we display the errors in
logarithmic scale as a function of $\lambda/\gamma$, which allows us to
visualize the convergence rates as the slopes of the curves. From Fig. 5.2, we
observe that the extension method has the first-order convergence rate, which
is faster than the result given in Theorem 1. A sharper error estimate for the
extension method will be investigated in our future work.

(a)

(b)

Fig. 4: Convergence rates of the extension (logarithmic scale). (a): Case 5.1;
(b): Case 5.2.

Next, we give some numerical results to demonstrate the efficiency and
accuracy of the proposed multiscale approach. Since in general it is
impossible to get the analytic solution $({\bf E}_{\eta},{\bf J}_{\eta})$ of
the original problem (21), in order to show the accuracy of our method, we
regard the numerical solution of the original system on a very fine mesh as
the reference solution.

(a)

(b)

Fig. 5: Display of the electric field on the line $x=y=z$ calculated by the
original multiscale method. (a): The real part of the $y$ component of the
electric field in Case 5.1 at $\omega=0.75\omega_{p}$; (b): The real part of
the $x$ component of the electric field in Case 5.2 at
$\omega=0.48\omega_{p}$.

We first solve the problem by the original multiscale approach presented at
the beginning of section 4. The numerical results for the electric field on
the line $x=y=z$ in Case 5.1 and 5.2 are plotted in Fig. 5 (a)-(b). We clearly
see that the multiscale approximate solution agrees well with the reference
solution outside the metallic nanostructures while it fails to capture the
oscillations of the electric field inside the metallic nanostructures.

To illustrate why the original multiscale approach fails inside the metallic
nanostructures, we compute the minimal eigenvalue (denoted by $\alpha$) of the
imaginary part of the homogenized coefficient $\widehat{\gamma^{\ast}}$ by
solving the cell problem (55) numerically. The values of $\alpha/\gamma$ are
displayed in Fig. 5.4 as a function of $\lambda/\gamma$, from which we see
that the minimal eigenvalue $\alpha$ tends to infinity as the extension
parameter $\lambda\rightarrow\infty$. This observation combined with Lemma 3
show that the electric field computed by the original multiscale approach with
a large extension parameter loses the nonlocal information, leading to the
failure of the original multiscale approach.

Fig. 6: The minimal eigenvalue $\alpha$ of $\widehat{\gamma^{\ast}}$
(logarithmic scale).

Next we use the modified multiscale approach to solve the problem. The
computational costs of solving the problem directly and solving by the
modified multiscale approach for Case 5.1 at $\omega=0.75\omega_{p}$ are given
in Table 2, which clearly show that our method is capable to reduce the
computational cost dramatically.

Table 2: Comparison of computational costs for Case 5.1 at $\omega=0.75\omega_{p}$. | original problem | cell problem | homogenized problem | modified problem
---|---|---|---|---
Elements | 2166403 | 68053 | 514247 | 12999
Dof | 31032404 | 13749 | 1231616 | 197766
CPU Time (s) | 3881.9739 | 1.2144 | 58.3542 | 23.4486

Without confusion, we denote by $({\bf E}_{\eta},{\bf J}_{\eta})$, ${\bf
E}^{0}$, ${\bf E}^{0}_{\eta}$, $({\bf E}^{M}_{\eta},{\bf J}^{M}_{\eta})$ the
reference solution of the original problem, the homogenized solution for the
electric field by solving the homogenized Maxwell equations (76) numerically,
the multiscale approximate solution for the electric field, and the modified
multiscale approximate solution by solving (77) inside each metallic
nanostructure numerically, respectively. The numerical errors of the modified
multiscale approach for Case 5.1 and 5.2 are displayed in Table 3 and 4,
respectively. Here we abbreviate $\|\cdot\|_{\mathbf{H}_{T}({\bf
curl};\,\Omega)}$ and $\|\cdot\|_{\mathbf{H}({\bf div};\,\Omega_{s})}$ as
$\|\cdot\|_{\Omega}$ and $\|\cdot\|_{\Omega_{s}}$ for brevity.

Table 3: Comparison of the numerical errors in Case 5.1. $\omega/\omega_{p}$ | $\frac{\|{\bf E}^{0}-\mathbf{E}_{\eta}\|_{\Omega}}{\|\mathbf{E}_{\eta}\|_{\Omega}}$ | $\frac{\|{\bf E}^{0}_{\eta}-\mathbf{E}_{\eta}\|_{\Omega}}{\|\mathbf{E}_{\eta}\|_{\Omega}}$ | $\frac{\|{\bf E}^{M}_{\eta}-\mathbf{E}_{\eta}\|_{\Omega}}{\|\mathbf{E}_{\eta}\|_{\Omega}}$ | $\frac{\|\mathbf{J}^{M}_{\eta}-\mathbf{J}_{\eta}\|_{\Omega_{s}}}{\|\mathbf{J}_{\eta}\|_{\Omega_{s}}}$
---|---|---|---|---
0.55 | 0.021147 | 0.021449 | 0.017761 | 0.184561
0.65 | 0.019662 | 0.015671 | 0.014059 | 0.178703
0.75 | 0.020467 | 0.013447 | 0.012417 | 0.146900
Table 4: Comparison of the numerical errors in Case 5.2. $\omega/\omega_{p}$ | $\frac{\|{\bf E}^{0}-\mathbf{E}_{\eta}\|_{\Omega}}{\|\mathbf{E}_{\eta}\|_{\Omega}}$ | $\frac{\|{\bf E}^{0}_{\eta}-\mathbf{E}_{\eta}\|_{\Omega}}{\|\mathbf{E}_{\eta}\|_{\Omega}}$ | $\frac{\|{\bf E}^{M}_{\eta}-\mathbf{E}_{\eta}\|_{\Omega}}{\|\mathbf{E}_{\eta}\|_{\Omega}}$ | $\frac{\|\mathbf{J}^{M}_{\eta}-\mathbf{J}_{\eta}\|_{\Omega_{s}}}{\|\mathbf{J}_{\eta}\|_{\Omega_{s}}}$
---|---|---|---|---
0.48 | 0.028744 | 0.028645 | 0.008425 | 0.091273
0.58 | 0.035826 | 0.034270 | 0.009712 | 0.046189
0.68 | 0.034353 | 0.024464 | 0.022433 | 0.140880

To show the accuracy of the modified multiscale method, in Fig. 7 we compare
the modified multiscale approximate solution with the reference solution for
the electric field in the same setting as Fig. 5. We find that the modified
multiscale approximate solution for the electric field is in good agreement
with the reference solution in the whole computational domain.

(a)

(b)

Fig. 7: Display of the electric field on the line $x=y=z$ calculated by the
modified multiscale method. (a): The real part of the $y$ component of the
electric field in Case 5.1 at $\omega=0.75\omega_{p}$; (b): The real part of
the $x$ component of the electric field in Case 5.2 at
$\omega=0.48\omega_{p}$.

Fig. 8 (a)-(d) display the numerical results for the electric field based on
the modified multiscale approach on the intersection $y=1.5$ nm at
$\omega=0.75\omega_{p}$ in Case 5.1.

In Fig. 9 (a)-(d), we compare the multiscale approximate solution with the
reference solution for the polarization current on the intersection $y=1.5$ nm
in Case 5.1 and 5.2.

(a) (b) (c) (d)

Fig. 8: The $x$-component of the electric field $\mathbf{E}$ on the
intersection $y=1.5$ nm at $\omega=0.75\omega_{p}$ in Case 5.1: (a) the
homogenized solution $|\mathbf{E}^{0}_{x}|$; (b) the multiscale approximate
solution $|\mathbf{E}^{0}_{\eta,x}|$; (c) the modified multiscale approximate
solution $|\mathbf{E}^{M}_{\eta,x}|$; (d) the reference solution
$|{\mathbf{E}}_{\eta,x}|$.

(a) (b) (a) (b)

Fig. 9: The $x$-component of the polarization current $\mathbf{J}$ on the
intersection $y=1.5$ nm: (a) the modified multiscale approximate solution
$|\mathbf{J}^{M}_{\eta,x}|$ in Case 5.1 at $\omega=0.75\omega_{p}$; (b) the
reference solution $|\mathbf{J}_{\eta,x}|$ in Case 5.1 at
$\omega=0.75\omega_{p}$; (c) the modified multiscale approximate solution
$|\mathbf{J}^{M}_{\eta,x}|$ in Case 5.2 at $\omega=0.48\omega_{p}$; (d) the
reference solution $|\mathbf{J}_{\eta,x}|$ in Case 5.2 at
$\omega=0.48\omega_{p}$.

##  6 Conclusions

We have presented a novel multiscale approach to simulate the nonlocal optical
response of metallic nanostructure arrays embedded in a dielectric medium by
solving a nonlocal hydrodynamic Drude (NHD) model with rapidly oscillating
coefficients. Our approach contains two key steps. First, we solve a
homogenized problem on the whole domain, which is much easier to solve than
the original system since the coefficients are constant. Next, we solve the
original system in each metallic nanostructure separately to capture the
nonlocal information of the electric field. A fast algorithm based on the $LU$
decomposition is proposed to solve the resulting linear system. In this way,
we make the numerical results of our approach agree well with the reference
solution on the whole domain. More importantly, the computational burden of
our approach is much lower than that of solving the original system directly.
Numerical experiments are presented to demonstrate the accuracy and efficiency
of the proposed multiscale approach.

Based on the work in this paper, in the near future we will study applications
of the multiscale approach to the time-dependent NHD model to simulate the
transient nonlocal optical response of metallic nanostructure arrays.
Furthermore, from the view of theoretical analysis, another focus of future
works is to prove the convergence of the homogenization method developed in
this paper.

##  Appendix A Proof of the error bound for the extension

We provide a proof of Theorem 1. To begin with, we introduce some notations
used in the proof. Let ${L}^{p}(\Omega)$ and
$\mathbf{L}^{p}(\Omega)=[L^{p}(\Omega)]^{3}$ be the Lebesgue spaces of
complex-valued functions and vector-valued functions, respectively. $L^{2}$
inner-products in $L^{2}(\Omega)$ and $\mathbf{L}^{2}(\Omega)$ are denoted by
$(\cdot,\cdot)$ without ambiguity. To avoid confusion, we use
$(\cdot,\cdot)_{s}$ and $(\cdot,\cdot)_{s,\eta}$ to denote the $L^{2}$ inner-
products in $L^{2}(\Omega_{s})$ ($\mathbf{L}^{2}(\Omega_{s})$) and
$L^{2}(\Omega_{s,\eta})$ ($\mathbf{L}^{2}(\Omega_{s,\eta})$), respectively.

We define

$\begin{array}[]{lll}{\displaystyle\mathbf{H}(\mathbf{curl};\Omega)=\\{\mathbf{u}\in\mathbf{L}^{2}(\Omega)\,|\,\,\mathbf{curl}\,\mathbf{u}\in\mathbf{L}^{2}(\Omega)\\},}\\\\[5.69054pt]
{\displaystyle\mathbf{H}(\mathbf{div};\Omega)=\\{\mathbf{u}\in\mathbf{L}^{2}(\Omega)\,|\,\,\mathbf{div}\,\mathbf{u}\in
L^{2}(\Omega)\\},}\\\\[5.69054pt] \end{array}$ (83)

which are equipped with the norms

$\begin{array}[]{lll}{\displaystyle\|\mathbf{u}\|_{\mathbf{H}(\mathbf{curl};\Omega)}=\|\mathbf{u}\|_{\mathbf{L}^{2}(\Omega)}+\|\mathbf{curl}\,\mathbf{u}\|_{\mathbf{L}^{2}(\Omega)},}\\\\[5.69054pt]
{\displaystyle\|\mathbf{u}\|_{\mathbf{H}(\mathbf{div};\Omega)}=\|\mathbf{u}\|_{\mathbf{L}^{2}(\Omega)}+\|\mathbf{div}\,\mathbf{u}\|_{L^{2}(\Omega)}.}\end{array}$

In addition,

$\begin{array}[]{lll}{\displaystyle\mathbf{H}_{T}(\mathbf{curl};\Omega)=\\{\mathbf{u}\in\mathbf{H}(\mathbf{curl};\Omega)\,|\,\,\mathbf{u}_{T}=({\bf
n}\times{\bf u})\times{\bf n}\in{\bf L}^{2}(\partial\Omega)\,\,\,{\rm
on}\,\,\partial\Omega\\},}\\\\[5.69054pt]
{\displaystyle\mathbf{H}_{0}(\mathbf{curl};\Omega)=\\{\mathbf{u}\in\mathbf{H}(\mathbf{curl};\Omega)\,|\,\,\mathbf{u}\times\mathbf{n}={\bf
0}\,\,\,{\rm on}\,\,\partial\Omega\\},}\\\\[5.69054pt]
{\displaystyle\mathbf{H}_{0}(\mathbf{div};\Omega)=\\{\mathbf{u}\in\mathbf{H}(\mathbf{div};\Omega)\,|\,\,\mathbf{u}\cdot\mathbf{n}=0\,\,\,{\rm
on}\,\,\partial\Omega\\}.}\\\\[5.69054pt] \end{array}$ (84)

Functions in $\mathbf{H}_{T}(\mathbf{curl};\Omega)$ are equipped with the norm

$\|\mathbf{u}\|_{\mathbf{H}_{T}(\mathbf{curl};\Omega)}=\|\mathbf{u}\|_{\mathbf{L}^{2}(\Omega)}+\|\mathbf{curl}\,\mathbf{u}\|_{\mathbf{L}^{2}(\Omega)}+\|\mathbf{u}_{T}\|_{\mathbf{L}^{2}(\partial\Omega)}.$
(85)

For the sake of convenience, we denote by

$\begin{array}[]{lll}{\displaystyle
X(\Omega)=\mathbf{H}_{T}(\mathbf{curl};\Omega),\;\qquad\|\mathbf{u}\|_{X(\Omega)}=\|\mathbf{u}\|_{\mathbf{H}_{T}(\mathbf{curl};\Omega)},}\\\\[5.69054pt]
{\displaystyle
Y(\Omega_{s})=\mathbf{H}_{0}(\mathbf{div};\Omega_{s}),\;\qquad\|\mathbf{v}\|_{Y(\Omega_{s})}=\|\mathbf{v}\|_{\mathbf{H}(\mathbf{div};\Omega_{s})},}\\\\[5.69054pt]
{\displaystyle
Y(\Omega_{s,\eta})=\mathbf{H}_{0}(\mathbf{div};\Omega_{s,\eta}),\quad\|\mathbf{v}\|_{Y(\Omega_{s,\eta})}=\|\mathbf{v}\|_{\mathbf{H}(\mathbf{div};\Omega_{s,\eta})}.}\end{array}$
(86)

The weak formulation of the original system (21) reads as follows: Find
$(\mathbf{E}_{\eta},{\bf J}_{\eta})\in X(\Omega)\times Y(\Omega_{s,\eta})$,
such that the equations

$\left\\{\begin{array}[]{@{}l@{}}{\displaystyle\big{(}\tilde{\mu}_{\eta}^{-1}\mathbf{curl}\,{\bf
E}_{\eta},\,\mathbf{curl}\,{\bf u})-\omega^{2}(\tilde{\varepsilon}_{\eta}{\bf
E}_{\eta},\,{\bf u})-{\mathrm{{i}}}\omega\langle{\bf E}_{\eta,T},{\bf
u}_{T}\rangle=\langle{\bf g},{\bf
u}_{T}\rangle+{\mathrm{{i}}}\omega\big{(}{\bf J}_{\eta},\,{\bf
u}\big{)}_{s,\eta}}\\\\[5.69054pt]
{\displaystyle\big{(}\beta^{2}\mathbf{div}\,{\bf
J}_{\eta},\,\mathbf{div}\,{\bf
v}\big{)}_{s,\eta}-(\omega(\omega+{\mathrm{{i}}}\gamma){\bf J}_{\eta},\,{\bf
v})_{s,\eta}=-{\mathrm{{i}}}\omega\omega^{2}_{p}\varepsilon_{0}\big{(}{\bf
E}_{\eta},\,{\bf v}\big{)}_{s,\eta}}\end{array}\right.$ (87)

hold for each $({\bf u},{\bf v})\in X(\Omega)\times Y(\Omega_{s,\eta})$, where
$\tilde{\mu}_{\eta}$ and $\tilde{\varepsilon}_{\eta}$ are given by

$\tilde{{\mu}}_{\eta}=\left\\{\begin{array}[]{ll}{{\mu}_{\eta}},&{\rm
in}\,\,\;\Omega_{s}\\\\[5.69054pt] {\mu_{0}},&{\rm
in}\,\,\;\Omega/{\Omega}_{s},\end{array}\right.\qquad\tilde{{\varepsilon}}_{\eta}=\left\\{\begin{array}[]{ll}{{\varepsilon}_{\eta}},&{\rm
in}\,\,\;\Omega_{s}\\\\[5.69054pt] {\varepsilon_{0}},&{\rm
in}\,\,\;\Omega/{\Omega}_{s},\end{array}\right.$ (88)

${\bf u}_{T}=({\bf n}\times{\bf u})\times{\bf n}$ with ${\bf n}$ being the
unit outward normal to $\partial\Omega$, and $\langle\cdot,\cdot\rangle$
denotes the $L^{2}$ inner product in $\mathbf{L}^{2}(\partial\Omega)$. Note
that here we use $\mathbf{E}_{\lambda}$ to denote the electric field on the
whole domain for convenience.

Similarly, we can give the weak formulation for the extended system (32). Find
$(\mathbf{E}_{\eta,\lambda},{\bf J}_{\eta,\lambda})\in X(\Omega)\times
Y(\Omega_{s})$, such that the equations

$\left\\{\begin{array}[]{@{}l@{}}{\displaystyle\big{(}\tilde{\mu}_{\eta}^{-1}\mathbf{curl}\,{\bf
E}_{\eta,\lambda},\,\mathbf{curl}\,{\bf
u})-\omega^{2}(\tilde{\varepsilon}_{\eta}{\bf E}_{\eta,\lambda},\,{\bf
u})-{\mathrm{{i}}}\omega\langle{\bf E}_{\eta,\lambda,T},{\bf
u}_{T}\rangle=\langle{\bf g},{\bf
u}_{T}\rangle+{\mathrm{{i}}}\omega\big{(}{\bf J}_{\eta,\lambda},\,{\bf
u}\big{)}_{s}}\\\\[5.69054pt]
{\displaystyle\big{(}\beta^{2}_{\eta,\lambda}\mathbf{div}\,{\bf
J}_{\eta,\lambda},\,\mathbf{div}\,{\bf
v}\big{)}_{s}-(\omega(\omega+{\mathrm{{i}}}\gamma_{\eta,\lambda}){\bf
J}_{\eta,\lambda},\,{\bf
v})_{s}=-{\mathrm{{i}}}\omega\omega^{2}_{p}\varepsilon_{0}\big{(}{\bf
E}_{\eta,\lambda},\,{\bf v}\big{)}_{s}}\end{array}\right.$ (89)

hold for each $({\bf u},{\bf v})\in X(\Omega)\times Y(\Omega_{s})$. In the
same spirit, here $\mathbf{E}_{\eta,\lambda}$ denotes the electric field on
the whole domain.

The following lemma gives the well-posedness of the equations (87) and (89).

######  Lemma 4.

[12, Theorem 3.1] Let $\Omega$, $\Omega_{s}$, and $\Omega_{s,\eta}$ be the
bounded, simply-connected, Lipschitz domains in $\mathbb{R}^{3}$ with
$\bar{\Omega}_{s}\subset\Omega$ and $\bar{\Omega}_{s,\eta}\subset\Omega_{s}$.
Assumes that $\mathbf{g}\in\mathbf{L}^{2}(\partial\Omega)$. There exist a
unique solution $({\bf E}_{\eta},\,{\bf J}_{\eta})\in X(\Omega)\times
Y(\Omega_{s,\eta})$ to the equations (87) satisfying

$\|{\bf E}_{\eta}\|_{X(\Omega)}+\|{\bf J}_{\eta}\|_{Y(\Omega_{s,\eta})}\leq
C\|{\bf g}\|_{{\bf L}^{2}(\partial\Omega)},$ (90)

where $C$ is a positive constant. Similarly, there exist a unique solution
$({\bf E}_{\eta,\lambda},\,{\bf J}_{\eta,\lambda})\in X(\Omega)\times
Y(\Omega_{s})$ to the equations (89) satisfying

$\|{\bf E}_{\eta,\lambda}\|_{X(\Omega)}+\|{\bf
J}_{\eta,\lambda}\|_{Y(\Omega_{s})}\leq C\|{\bf g}\|_{{\bf
L}^{2}(\partial\Omega)}.$ (91)

In order to estimate the error between $({\bf E}_{\eta},{\bf{J}}_{\eta})$ and
$({\bf E}_{\eta,\lambda},{\bf{J}}_{\eta,\lambda})$, we need a decay estimate
for the extended polarization current ${\bf J}_{\eta,\lambda}$ in
$\Omega_{s}/\Omega_{s,\eta}$ with respect to $\lambda$.

######  Lemma 5.

Let $({\bf E}_{\eta,\lambda},{\bf J}_{\eta,\lambda})$ be the solution of the
extended system (89). There exists a constant $C>0$ independent of $\lambda$
such that

$\|{\bf J}_{\eta,\lambda}\|_{\mathbf{L}^{2}(\Omega_{s}/\Omega_{s,\eta})}\leq
C/{\lambda}^{1/2},\,\quad\|\mathbf{div}\,{\bf
J}_{\eta,\lambda}\|_{\mathbf{L}^{2}(\Omega_{s}/\Omega_{s,\eta})}\leq
C/{\lambda}^{1/2}.$ (92)

###### Proof.

By choosing the test function $({\bf u},{\bf v})=({\bf E}_{\eta,\lambda},{\bf
J}_{\eta,\lambda})$ in (89), and taking the imaginary part of the two
equations of (89), we obtain

$\begin{array}[]{@{}l@{}}{\displaystyle\omega\|{\bf
E}_{\eta,\lambda,T}\|^{2}_{\mathbf{L}^{2}(\partial\Omega)}=-\omega\,{\rm
Re}\big{(}{\bf J}_{\eta,\lambda},\,{\bf E}_{\eta,\lambda}\big{)}_{s}-{\rm
Im}\langle{\bf g},{\bf E}_{\eta,\lambda,T}\rangle,}\\\\[5.69054pt]
{\displaystyle\big{(}\gamma_{\eta,\lambda}{\bf J}_{\eta,\lambda},\,{\bf
J}_{\eta,\lambda}\big{)}_{s}=\omega^{2}_{p}\varepsilon_{0}\,{\rm
Re}\big{(}{\bf E}_{\eta,\lambda},\,{\bf
J}_{\eta,\lambda}\big{)}_{s},}\end{array}$ (93)

which implies that

$\begin{array}[]{lll}{\displaystyle\big{(}\gamma_{\eta,\lambda}{\bf
J}_{\eta,\lambda},\,{\bf
J}_{\eta,\lambda}\big{)}_{s}+\omega^{2}_{p}\varepsilon_{0}\|{\bf
E}_{\eta,\lambda,T}\|^{2}_{\mathbf{L}^{2}(\partial\Omega)}\leq\frac{\omega^{2}_{p}\varepsilon_{0}}{\omega}|\langle{\bf
g},{\bf E}_{\eta,\lambda,T}\rangle|}\\\\[5.69054pt]
{\displaystyle\qquad\leq\frac{\omega^{2}_{p}\varepsilon_{0}}{2}\|{\bf
E}_{\eta,\lambda,T}\|^{2}_{\mathbf{L}^{2}(\partial\Omega)}+C\|{\bf
g}\|^{2}_{\mathbf{L}^{2}(\partial\Omega)}.}\end{array}$ (94)

Since $\gamma_{\eta,\lambda}|_{\Omega_{s}/\Omega_{s,\eta}}=\lambda$, from (94)
we deduce

$\lambda\|{\bf
J}_{\eta,\lambda}\|^{2}_{\mathbf{L}^{2}(\Omega_{s}/\Omega_{s,\eta})}\leq\big{(}\gamma_{\eta,\lambda}{\bf
J}_{\eta,\lambda},\,{\bf J}_{\eta,\lambda}\big{)}_{s}\leq C\|{\bf
g}\|^{2}_{\mathbf{L}^{2}(\partial\Omega)},$ (95)

which gives the first inequality in (92). By applying Theorem 4.17 of [15] to
the first equation of (89), we further have

$\|{\bf E}_{\eta,\lambda}\|_{X(\Omega)}\leq C(\|{\bf
g}\|_{\mathbf{L}^{2}(\partial\Omega)}+\|{\bf
J}_{\eta,\lambda}\|_{\mathbf{L}^{2}(\Omega_{s})})\leq C\|{\bf
g}\|_{\mathbf{L}^{2}(\partial\Omega).}$ (96)

Next we choose ${\bf v}={\bf J}_{\eta,\lambda}$ in the second equation of (89)
and take the real part of the equation to get

$\|\beta_{\eta,\lambda}\mathbf{div}\,{\bf J}_{\eta,\lambda}\|^{2}_{{\bf
L}^{2}(\Omega_{s})}=\omega^{2}\|{\bf J}_{\eta,\lambda}\|^{2}_{{\bf
L}^{2}(\Omega_{s})}+\omega\omega_{p}^{2}\varepsilon_{0}{\rm Im}\,\big{(}{\bf
E}_{\eta,\lambda},\,{\bf J}_{\eta,\lambda}\big{)}_{s}\leq C\|{\bf
g}\|^{2}_{\mathbf{L}^{2}(\partial\Omega)},$ (97)

where we have used (96). Then the second inequality in (92) follows
immediately from (97) and the fact that
$\beta_{\eta,\lambda}|_{\Omega_{s}/\Omega_{s,\eta}}=\lambda^{\frac{1}{2}}$. ∎

Next we construct a function $\mathbf{w}\in{\bf H}(\rm{\bf
div},\Omega_{s,\eta})$ that has the same normal trace on
$\partial\Omega_{s,\eta}$ as ${\bf J}_{\eta,\lambda}$.

######  Lemma 6.

Let $({\bf E}_{\eta,\lambda},{\bf J}_{\eta,\lambda})$ be the solution of the
extended system (89). Then there exists a $\mathbf{w}\in{\bf H}(\rm{\bf
div},\Omega_{s,\eta})$ such that
$\mathbf{n}\cdot\mathbf{w}=\mathbf{n}\cdot{\bf J}_{\eta,\lambda}$ on
$\partial\Omega_{s,\eta}$, and

$\|{\bf w}\|_{{\bf H}(\rm{\bf div},\Omega_{s,\eta})}\leq C/\lambda^{1/2},$
(98)

where $C$ is a positive constant independent of $\lambda$.

###### Proof.

By Lemma 5 and [15, Theorem 3.24], we find

$\|\mathbf{n}\cdot{\bf
J}_{\eta,\lambda}\|_{{H}^{-\frac{1}{2}}(\partial\Omega_{s,\eta})}\leq C\|{\bf
J}_{\eta,\lambda}\|_{{\bf H}({\rm{\bf div}},\Omega_{s}/\Omega_{s,\eta})}\leq
C/\lambda^{1/2}.$ (99)

Let $\phi\in H^{1}(\Omega_{s,\eta})$ be the solution of the Neumann problem

$(\mathbf{grad}\,\phi,\mathbf{grad}\,\psi)_{s,\eta}+(\phi,\psi)_{s,\eta}=<\mathbf{n}\cdot{\bf
J}_{\eta,\lambda},\,\psi>_{{H}^{-\frac{1}{2}}(\partial\Omega_{s,\eta}),{H}^{\frac{1}{2}}(\partial\Omega_{s,\eta})}$
(100)

for all $\psi\in H^{1}(\Omega_{s,\eta})$. Taking $\phi$ as the test function
in (100) gives

$\begin{array}[]{lll}{\displaystyle\|\phi\|_{{H}^{1}(\Omega_{s,\eta})}^{2}\leq\,|<\mathbf{n}\cdot{\bf
J}_{\eta,\lambda},\,\phi>_{{H}^{-\frac{1}{2}}(\partial\Omega_{s,\eta}),{H}^{\frac{1}{2}}(\partial\Omega_{s,\eta})}|}\\\\[5.69054pt]
{\displaystyle\quad\leq\|\mathbf{n}\cdot{\bf
J}_{\eta,\lambda}\|_{{H}^{-\frac{1}{2}}(\partial\Omega_{s,\eta})}\|\phi\|_{{H}^{\frac{1}{2}}(\partial\Omega_{s,\eta})}\leq
C\|\mathbf{n}\cdot{\bf
J}_{\eta,\lambda}\|_{{H}^{-\frac{1}{2}}(\partial\Omega_{s,\eta})}\|\phi\|_{{H}^{1}(\Omega_{s,\eta})},}\end{array}$

which yields

$\|\phi\|_{{H}^{1}(\Omega_{s,\eta})}\leq C/\lambda^{1/2}$ (101)

by a use of (100).

Let ${\bf w}=\mathbf{grad}\,\phi\in{\bf L}^{2}(\Omega_{s,\eta})$. Taking
$\psi\in\mathcal{C}_{0}^{\infty}(\Omega_{s,\eta})$ in (100), we see that

$({\bf w},\mathbf{grad}\,\psi)_{s,\eta}+(\phi,\psi)_{s,\eta}=0$ (102)

for all $\psi\in\mathcal{C}_{0}^{\infty}(\Omega_{s,\eta})$, which implies that
$\mathbf{div}\,{\bf w}=\phi\in L^{2}(\Omega_{s,\eta})$. Thus ${\bf w}\in{\bf
H}(\rm{\bf div},\Omega_{s,\eta})$ satisfies
$\mathbf{n}\cdot\mathbf{w}=\mathbf{n}\cdot{\bf J}_{\eta,\lambda}$ on
$\partial\Omega_{s,\eta}$, and

$\|{\bf w}\|_{{\bf H}(\rm{\bf
div},\Omega_{s,\eta})}=\|\phi\|_{{H}^{1}(\Omega_{s,\eta})}\leq
C/\lambda^{1/2},$ (103)

by using (101). ∎

Now we turn to the proof of Theorem 1.

Let ${\bf e}={\bf E}_{\eta,\lambda}-{\bf E}_{\eta}$ and ${\bf j}={\bf
J}_{\eta,\lambda}-{\bf J}_{\eta}-{\bf w}$. Note that
$\mathbf{n}\cdot\mathbf{j}=0$ on $\partial\Omega_{s,\eta}$ and thus
$\mathbf{j}\in Y(\Omega_{s,\eta})$. By subtracting (87) from (89) and
recalling that $\gamma_{\eta,\lambda}|_{\Omega_{s,\eta}}=\gamma$ and
$\beta^{2}_{\eta,\lambda}|_{\Omega_{s,\eta}}=\beta^{2}$, we see that $({\bf
e},\,{\bf j})$ satisfies

$\left\\{\begin{array}[]{@{}l@{}}{\displaystyle\big{(}\tilde{\mu}_{\eta}^{-1}\mathbf{curl}\,{\bf
e},\,\mathbf{curl}\,{\bf u})-\omega^{2}(\tilde{\varepsilon}_{\eta}{\bf
e},\,{\bf u})-{\mathrm{{i}}}\omega\langle{\bf e}_{T},{\bf
u}_{T}\rangle={\mathrm{{i}}}\omega\big{(}{\bf j},\,{\bf
u}\big{)}_{s,\eta}}\\\\[5.69054pt]
{\displaystyle\qquad\qquad\qquad\qquad\qquad+\,{\mathrm{{i}}}\omega\big{(}{\bf{\bf{w}}},\,{\bf
u}\big{)}_{s,\eta}+\mathrm{i}\omega\big{(}\mathbf{J}_{\eta,\lambda},\,\mathbf{u}\big{)}_{\Omega_{s}/\Omega_{s,\eta}},}\\\\[8.53581pt]
{\displaystyle\beta^{2}\big{(}\mathbf{div}\,{\bf j},\,\mathbf{div}\,{\bf
v}\big{)}_{s,\eta}-\omega(\omega+{\mathrm{{i}}}\gamma)({\bf j},\,{\bf
v})_{s,\eta}=-{\mathrm{{i}}}\omega\omega_{p}^{2}\varepsilon_{0}\big{(}{\bf
e},\,{\bf v}\big{)}_{s,\eta}}\\\\[5.69054pt]
{\displaystyle\quad\quad\quad\quad\quad\quad\quad-\beta^{2}\big{(}\mathbf{div}\,{\bf
w},\,\mathbf{div}\,{\bf
v}\big{)}_{s,\eta}+\omega(\omega+{\mathrm{{i}}}\gamma)({\bf w},\,{\bf
v})_{s,\eta}}\end{array}\right.$ (104)

for all $({\bf u},{\bf v})\in X(\Omega)\times Y(\Omega_{s,\eta})$. We will
prove that $({\bf e},\,{\bf j})$ satisfies the estimate

$\|{\bf e}\|_{X(\Omega)}+\|{\bf j}\|_{Y(\Omega_{s,\eta})}\leq
C/{\lambda}^{1/2},$ (105)

where $C>0$ is independent of $\lambda$. Theorem 1 follows from (105) and
Lemma 6.

To prove (105), we first choose
$(\mathbf{u},\,\mathbf{v})=(\mathbf{e},\,\mathbf{j})$ in (104) and take the
imaginary part of the equations to get

$\begin{array}[]{lll}{\displaystyle\|\mathbf{e}_{T}\|^{2}_{\mathbf{L}^{2}(\partial\Omega)}=-{\rm
Re}\,(\mathbf{j},\,\mathbf{e})_{s,\eta}-{\rm
Re}\,(\mathbf{w},\,\mathbf{e})_{s,\eta}-\mathrm{Re}\,(\mathbf{J}_{\eta,\lambda},\,\mathbf{e})_{\Omega_{s}/\Omega_{s,\eta}},}\\\\[5.69054pt]
{\displaystyle\gamma\|\mathbf{j}\|^{2}_{\mathbf{L}^{2}(\Omega_{s,\eta})}=\omega^{2}_{p}\varepsilon_{0}\,\mathrm{Re}\,(\mathbf{e},\,\mathbf{j})_{s,\eta}-\gamma\,\mathrm{Re}\,(\mathbf{w},\,\mathbf{j})_{s,\eta}.}\end{array}$
(106)

Combining the two equations in (106) and using Lemma 5 and 6, we obtain

$\begin{array}[]{lll}{\displaystyle\omega^{2}_{p}\varepsilon_{0}\|\mathbf{e}_{T}\|^{2}_{\mathbf{L}^{2}(\partial\Omega)}+\gamma\|\mathbf{j}\|^{2}_{\mathbf{L}^{2}(\Omega_{s,\eta})}}\\\\[5.69054pt]
{\displaystyle\leq\omega^{2}_{p}\varepsilon_{0}|(\mathbf{w},\,\mathbf{e})_{s,\eta}|+\omega^{2}_{p}\varepsilon_{0}|(\mathbf{J}_{\eta,\lambda},\,\mathbf{e})_{\Omega_{s}/\Omega_{s,\eta}}|+\gamma|(\mathbf{w},\,\mathbf{j})_{s,\eta}|}\\\\[5.69054pt]
{\displaystyle\leq
C\big{(}\|\mathbf{w}\|_{\mathbf{L}^{2}(\Omega_{s,\eta})}+\|\mathbf{J}_{\eta,\lambda}\|_{\mathbf{L}^{2}(\Omega_{s}/\Omega_{s,\eta})}\big{)}\|\mathbf{e}\|_{\mathbf{L}^{2}(\Omega)}+\frac{\gamma}{2}\|\mathbf{j}\|^{2}_{\mathbf{L}^{2}(\Omega_{s,\eta})}+\frac{\gamma}{2}\|\mathbf{w}\|^{2}_{\mathbf{L}^{2}(\Omega_{s,\eta})}}\\\\[8.53581pt]
{\displaystyle\leq{C}{\lambda^{-\frac{1}{2}}}\|\mathbf{e}\|_{\mathbf{L}^{2}(\Omega)}+\frac{\gamma}{2}\|\mathbf{j}\|^{2}_{\mathbf{L}^{2}(\Omega_{s,\eta})}+{C}{\lambda^{-1}},}\end{array}$

which implies that

$\|\mathbf{j}\|_{\mathbf{L}^{2}(\Omega_{s,\eta})}\leq{C}{\lambda^{-\frac{1}{4}}}\|\mathbf{e}\|^{\frac{1}{2}}_{\mathbf{L}^{2}(\Omega)}+{C}{\lambda^{-\frac{1}{2}}}.$
(107)

Now we consider the first equation of (104) separately in which $\mathbf{j}$
and $\mathbf{w}$ are considered as given sources. By employing Theorem 4.17 of
[15] again, we find

$\|{\bf e}\|_{X(\Omega)}\leq
C\big{(}\|\mathbf{j}\|_{\mathbf{L}^{2}(\Omega_{s,\eta})}+\|\mathbf{w}\|_{\mathbf{L}^{2}(\Omega_{s,\eta})}+\|\mathbf{J}_{\eta,\lambda}\|_{\mathbf{L}^{2}(\Omega_{s}/\Omega_{s,\eta})}\big{)},$

and consequently, with a use of (107), Lemma 5, and Lemma 6, we come to

$\|{\bf
e}\|_{X(\Omega)}\leq{C}{\lambda^{-\frac{1}{4}}}\|\mathbf{e}\|^{\frac{1}{2}}_{\mathbf{L}^{2}(\Omega)}+{C}{\lambda^{-\frac{1}{2}}}\leq\frac{1}{2}\|\mathbf{e}\|_{\mathbf{L}^{2}(\Omega)}+{C}{\lambda^{-\frac{1}{2}}},$

which yields that

$\|{\bf e}\|_{X(\Omega)}\leq{C}{\lambda^{-\frac{1}{2}}}.$ (108)

Next we consider the second equation of (104) separately with given sources
$\mathbf{e}$ and $\mathbf{w}$. By applying a similar argument as in the proofs
of Theorem 4.17 in [15] and Theorem 3.1 in [12], we get the following estimate

$\|{\bf j}\|_{Y(\Omega_{s,\eta})}\leq
C\big{(}\|\mathbf{e}\|_{\mathbf{L}^{2}(\Omega_{s,\eta})}+\|{\bf w}\|_{{\bf
H}(\rm{\bf div},\Omega_{s,\eta})}\big{)}.$

With Lemma 6 and (108), it follows that

$\|{\bf j}\|_{Y(\Omega_{s,\eta})}\leq C{\lambda^{-\frac{1}{2}}}.$ (109)

Combining (108) and (109) gives (105) and we complete the proof of Theorem 1.

## References

  * [1] B. Belyaev and V. Tyurnev, Electrodynamic calculation of effective electromagnetic parameters of a dielectric medium with metallic nanoparticles of a given size, Journal of Experimental and Theoretical Physics, 127 (2018), pp. 608–619.
  * [2] Z. Chen, L. Wang, and W. Zheng, An adaptive multilevel method for time-harmonic maxwell equations with singularities, SIAM Journal on Scientific Computing, 29 (2007), pp. 118–138.
  * [3] T. C. Choy, Effective medium theory: principles and applications, vol. 165, Oxford University Press, 2015.
  * [4] C. Ciraci, J. B. Pendry, and D. R. Smith, Hydrodynamic model for plasmonics: a macroscopic approach to a microscopic problem, Chem. phys. Chem., 14 (2013), pp. 1109–1116.
  * [5] P. Drude, On the electron theory of metals, Annals of physics, 306 (1900), pp. 566–613.
  * [6] A. Eguiluz, S. Ying, and J. Quinn, Influence of the electron density profile on surface plasmons in a hydrodynamic model, Physical Review B, 11 (1975), p. 2118.
  * [7] S. Enoch, R. Quidant, and G. Badenes, Optical sensing based on plasmon coupling in nanoparticle arrays, Optics express, 12 (2004), pp. 3422–3427.
  * [8] Y. Eremin, A. Doicu, and T. Wriedt, The numerical scheme of the discrete sources method to analyze 3D plasmonic nanostructures accounting for the non-local effect, J. Comput. Phys., 388 (2019), pp. 357–370.
  * [9] K. R. Hiremath, L. Zschiedrich, and F. Schmidt, Numerical solution of nonlocal hydrodynamic drude model for arbitrary shaped nano-plasmonic structures using Nédélec finite elements, J. Comput. Phys., 231 (2012), pp. 5890–5896.
  * [10] Y. Huang, J. Li, and W. Yang, Theoretical and numerical analysis of a non-local dispersion model for light interaction with metallic nanostructures, Comput. Math. Appl., 72 (2016), pp. 921–932.
  * [11] L. Li, S. Lanteri, N. A. Mortensen, and M. Wubs, A hybridizable discontinuous Galerkin method for solving nonlocal optical response models, Comput. Phys. Commun., 219 (2017), pp. 99–107.
  * [12] C. Ma, Y. Zhang, and J. Zou, Mathematical and numerical analysis of a nonlocal drude model in nanoplasmonics, arXiv preprint arXiv:1906.04790, (2019).
  * [13] S. A. Maier, Plasmonics: fundamentals and applications, Springer Science & Business Media, 2007.
  * [14] S. Mokkapati, F. Beck, A. Polman, and K. Catchpole, Designing periodic arrays of metal nanoparticles for light-trapping applications in solar cells, Applied Physics Letters, 95 (2009), p. 053115.
  * [15] P. Monk, Finite element methods for Maxwell’s equations, Oxford University Press, 2003.
  * [16] J. Pitarke, V. Silkin, E. Chulkov, and P. Echenique, Theory of surface plasmons and surface-plasmon polaritons, Reports on progress in physics, 70 (2006), p. 1.
  * [17] S. Raza, S. I. Bozhevolnyi, M. Wubs, and N. A. Mortensen, Nonlocal optical response in metallic nanostructures, Journal of Physics: Condensed Matter, 27 (2015), p. 183204.
  * [18] M. B. Ross, C. A. Mirkin, and G. C. Schatz, Optical properties of one-, two-, and three-dimensional arrays of plasmonic nanostructures, The journal of physical chemistry C, 120 (2016), pp. 816–830.
  * [19] N. Schmitt, C. Scheid, S. Lanteri, A. Moreau, and J. Viquerat, A DGTD method for the numerical modeling of the interaction of light with nanometer scale metallic structures taking into account non-local dispersion effects, J. Comput. Phys., 316 (2016), pp. 396–415.
  * [20] N. Schmitt, C. Scheid, J. Viquerat, and S. Lanteri, Simulation of three-dimensional nanoscale light interaction with spatially dispersive metals using a high order curvilinear DGTD method, J. Comput. Phys., 373 (2018), pp. 210–229.
  * [21] B. Stupfel and R. Mittra, A theoretical study of numerical absorbing boundary conditions, IEEE transactions on antennas and propagation, 43 (1995), pp. 478–487.
  * [22] G. Toscano, S. Raza, A.-P. Jauho, N. A. Mortensen, and M. Wubs, Modified field enhancement and extinction by plasmonic nanowire dimers due to nonlocal response, Optics express, 20 (2012), pp. 4176–4188.
  * [23] F. Vidal-Codina, N. C. Nguyen, S.-H. Oh, and J. Peraire, A hybridizable discontinuous Galerkin method for computing nonlocal electromagnetic effects in three-dimensional metallic nanostructures, J. Comput. Phys., 355 (2018), pp. 548–565.
  * [24] W. Wang, M. Ramezani, A. I. Väkeväinen, P. Törmä, J. G. Rivas, and T. W. Odom, The rich photonic world of plasmonic nanoparticle arrays, Materials today, 21 (2018), pp. 303–314.
  * [25] L. Zhang, L. Cao, and J. Luo, Multiscale analysis and computation for a stationary schrödinger–poisson system in heterogeneous nanostructures, Multiscale Modeling & Simulation, 12 (2014), pp. 1561–1591.
  * [26] X. Zheng, M. Kupresak, R. Mittra, and G. A. Vandenbosch, A boundary integral equation scheme for simulating the nonlocal hydrodynamic response of metallic antennas at deep-nanometer scales, IEEE Trans. Antennas Propag., 66 (2018), pp. 4759–4771.
  * [27] X. Zheng, M. Kupresak, V. V. Moshchalkov, R. Mittra, and G. A. Vandenbosch, A potential-based formalism for modeling local and hydrodynamic nonlocal responses from plasmonic waveguides, IEEE Transactions on Antennas and Propagation, 67 (2019), pp. 3948–3960.