Nanoplasmonics is an active research field concerned with the study of optical properties of metallic nanostructures [4, 20]. The interaction of metallic nanostrucutes with light at optical frequencies produces the excitation of (localized) surface plasmons, i.e., the collective oscillation of conduction band electrons at the metal surface, leading to many unusual and fascinating properties. It enables the confinement of light at the nanoscale vicinity of metal surfaces, the enormous local fields enhancement of the incident wave, and the squeezing of light beyond the diffraction limit, providing unparalleled means for manipulation of light at the nanoscale. As a result, nanoplasmonics has found applications in different fields, such as near-field scanning microscopy , ultrasensitive sensing and detection , and plasmonic waveguilding . To be able to understand and make use of nanoplasmonic phenomena, an appropriate modeling for describing the optical response of metallic nanostructures is required.
Historically, the interaction of light with metals has been described by the classical theory of light-matter interaction in which light and matter (the free electrons) are described by Maxwell’s equations and Newtonian mechanics, respectively. The most widely used model to describe the optical response of metals is Drude’s model . In this model, the collective oscillation of conduction electrons in metals subject to driving optical fields is analysed within the framework of local-response approximation (LRA), where the material response occurs only in the point of space of the perturbation and there is no response at even short distances. Drude’s model has achieved a great success in the modeling of bulk metals. However, as the size of metallic structures shrinks down to nanometer scales, nonlocal interaction effects between electrons become predominant and Drude’s model is inadequate to explain experimentally observable phenomena.
In view of the limitation of Drude’s model, the nonlocal response theories for metallic nanostructures have gained considerable interest and attention in the past decade and some improved models have been proposed, including the hydrodynamic model , the nonlocal hydrodynamic Drude (NHD) model , and the generalized nonlocal optical response (GNOR) model . In the hydrodynamic model, the free electrons in metals are modeled as a charged fluid and described by hydrodynamic equations of Euler-type. The NHD and GNOR models are derived from the hydrodynamic model by the linear-response approximation. All these models are coupled to Maxwell’s equations (in both frequency domain and time domain) and form coupled systems of PDEs.
Due to a relatively simple form and the successful interpretation of observable nonlocal effects , the NHD model has drew much attention in recent years and become a popular model in the study of optical properties of metallic nanostructures. Meanwhile, numerical methods for solving Maxwell’s equations coupled to the NHD model have been extensively studied. In , the frequency-domain NHD model (frequency-domain Maxwell’s equations coupled to the NHD model) is solved for modeling nano-plasmonic structures with complex geometries by using the Nédélec elements based finite element method. In 
, a computational scheme based on the boundary integral equation and method of moments is developed for the frequency-domain NHD model to predict the interaction of light with metallic nanoparticles. The discontinuous Galerkin methods for the time-domain and frequency-domain NHD models have been considered in[11, 18, 24]. For other more numerical methods for this system of PDEs, we refer the reader to [6, 10, 19, 23] and references therein.
However, up to now, most existing studies on the NHD model focus on the development of numerical methods or the analysis of physical effects. The theoretical and numerical analysis of this model available in the literature is very limited. In , the well-posedness and the stability and convergence of numerical methods are proved for the modified time-domain NHD model. To the best of our knowledge, there seems to be no results on the mathematical and numerical analysis of the frequency-domain NHD model. In fact, just as it is for the frequency-domain and time-domain Maxwell’s equations, the proof of the well-posedness and numerical convergence for the frequency-domain NHD model is much difficult than its time-domain counterpart.
In this paper, we present a rigorous mathematical and numerical analysis of the frequency-domain NHD model for the first time. The existence and uniqueness of weak solutions to the equations are proved. A finite element method based on the Raviart–Thomas and Nédélec elements is developed for the equations and the convergence is proved. There are two aspects of this model that make the analysis challenging. First, the curl and div operators both have a large null space in the continuous and discrete levels, which must be removed from the function spaces by using the (discrete) Hemtholtz decompositions. To this end, an understanding of some properties of the continuous function spaces and the finite element spaces is required. Second, the bilinear forms in the weak formulation of the equations are not coercive, which brings difficulty to the analysis in both continuous and discrete levels. To overcome this problem, in the proof of the well-posedness, we first show the uniqueness of weak solutions and then apply the Fredholm alternative to prove the existence of weak solutions, while in the proof of the convergence of the finite element discretization, the theory of convergence of collectively compact operators developed in [13, 14] is used. Although our methods are somewhat similar to those used in  for the analysis of frequency-domain Maxwell’s equations, the proof presented in this paper is much more delicate due to the coupled system nature.
The rest of this paper is organized as follows. In section 2, we give a brief derivation of the NHD model and describe the problem considered in this paper. In section 3, we prove the existence and uniqueness of weak solutions to the equations. In section 4, we propose a finite element discretization for the equations and prove the convergence of the scheme. In section 5, we give some numerical examples to confirm our theoretical analysis.
2 Nonlocal hydrodynamic Drude model
In this section, we briefly introduce the NHD model and give the problem considered in this paper.
In the absence of external charge and current, macroscopic Maxwell’s equations for metals can be written as
The equations link four macroscopic fields (the electric field), (the magnetic field), (the dielectric displacement), and (the magnetic flux density) with the free charge and current densities and . Maxwell’s equations (1) are supplemented by the constitutive laws which link to and to via
Here is the magnetic permeability of metals and is the electric permittivity of metals that takes into account the polarization of bound electrons ( is the electric permittivity of vacuum).
We derive the NHD model starting from the hydrodynamic model within which the free electrons are modeled as a charged fluid with the Euler equations:
where is the electron charge, is the effective electron mass, is the electron density, is the hydrodynamic velocity, is the electron pressure, and is the damping constant. The term represents the Lorentz force. The polarization charge and current densities and of the free electrons are given as
The equations (1)-(4) form the self-consistent Euler–Maxwell coupled equations. In order to simplify the above equations, we linearize the equations (3) as in perturbation theory by expanding the physical fields in a non-oscillating term (e.g. the constant equilibrium electron density ) and a small first-order dynamic term. In this spirit, we can write the perturbation expansions for and
Similar expansions can be written for the electric and magnetic fields. Since in the absence of an external field , the nonlinear terms and vanish due to the linearization. By using the Thomas-Fermi model for the pressure term in (3), we can linearize it as
where is an important parameter representing the nonlocality related to the Fermi velocity . Using the assumptions above, we get the linearized hydrodynamic equation
and the linearized continuity equation
Replacing with in (10
) by Fourier transformation in the time domain, whereis the imaginary unit and is the angular frequency, and eliminating the magnetic field , we get Maxwell’s equations with the NHD model in frequency domain
By Fourier transformation in the space domain, we replace with in the second equation of (11), which gives
and then we obtain the spatially-dispersive (relative) permittivity for the metal
The parameter represents the level of nonlocality. As , we recover the classical Drude permittivity .
In this paper we consider the following equations
with the boundary conditions
Here and are the bounded, simply-connected, Lipschitz polyhedron domains in with . and are shown in Fig 2.1. The magnetic permeability and electric permittivity are two piecewise constant functions in the domain , namely
and are positive constants.
The hard-wall boundary conditions for the current density means that the electrons are confined within the metal and spill-out of electrons in free space is neglected. For Maxwell’s equations, we apply the first-order Silver–Müller boundary conditions 
3 Existence and uniqueness of the solutions
In this section, we study the well-posedness of the problem (14)-(15). To begin with, we introduce some notations. We denote as the conventional Sobolev spaces of complex-valued functions defined in and as the subspace of consisting of functions whose traces are zero on . Let and
be the Lebesgue spaces of complex-valued functions and vector-valued functions with 3 components, respectively.inner-products in and are denoted by without ambiguity. To avoid confusion, we use to denote the inner-products in and .
which are equipped with the norms
Functions in are equipped with the norm
For the sake of convenience, we denote by
hold for each , where and are given in (16) and with being the unit outward normal to . denotes the inner product in . and are defined as follows.
is the characteristic function of.
We now state the main result of this section. Let and be the bounded, simply-connected, Lipschitz polyhedron domains in with . The equations (22) exist a unique solution satisfying
where the constant might depend on , , , , , , and .
We first prove the uniqueness of solutions of (22). There is at most one solution of (22). Since (22) is a linear system, we only need to show that is the only solution of (22) with . To this end, we first choose in the first equation of (22) and take the imaginary part of the equation to obtain
Next by setting in the second equation of (22) and taking the imaginary part of the equation, we have
Thus we find that satisfies
Before proving the existence of solutions of (22), we give two useful lemmas. We have the following Helmholtz decompositions for and
The Helmholtz decomposition for was proved in Lemma 4.5 of . It is not difficult to show that is a closed subspace of . Therefore the Helmholtz decomposition for follows from the projection theorem. For every , we can write it as
where , . In particular, we can take to be divergence-free, i.e.,.
and are compactly embedded in and , respectively.
The compact embedding of was proved in Theorem 4.7 of  and the compactness property of can be proved by a similar trick.
Now that we know and , we can write any solution of (22) as
for some and . In addition, we assume that satisfies .
Now taking in (33), where and , we obtain
where we have used the fact that and . By introducing a Lagrangian multiplier , we can rewrite (35) as
Next we take in (33) to obtain
We now define the sesquilinear form by
for all . For convenience, we introduce the following notations.
It is not difficult to prove the following lemma. There exists a constant depending on , , , , , , and such that
To proceed further, we define a map : such that if then satisfies
where and are the solutions of
We have the following result. The operator is a bounded and compact map from into . Moreover,
This theorem can be proved by using the Lax–Milgram theorem. To begin with, we check the conditions of the Lax–Milgram theorem. It is not difficult to show that is bounded. That is, there exists a constant , such that
holds for all . Coercivity of has been given in Lemma 3. It remains to show that there exists a constant such that
Since satisfies (42), it is easy to see that
Combining (48) and (49), we have (47), which yields (46). Having verified the conditions of the Lax–Milgram theorem, we know is well defined and obtain (44). The compactness of can be proved by applying a similar argument in Theorem 4.11 of  and we omit the proof here.
Next we define a vector which satisfies
By using the Lax–Milgram theorem again, we see that is well defined and
By virtue of the operator , we find that the problem (38) is equivalent to finding such that
Note that (52) implies that , from which we deduce
Under the assumptions of Theorem 3, there exists a such that for all with , . If in addition, is a convex polyhedron, we have . Moreover, if is constant on the whole domain , then . By taking with in the second equation of (22) and using (24), we get
which implies that and thus . Applying Theorem 3.50 of , we know that there is a such that holds for all with . If is a convex polyhedron, then is continuously embedded into (Theorem 3.9 of ). Therefore we have .
4 Finite element approximation
In this section, we present the finite element approximation for the system (22) and prove the convergence of the scheme. Our proof relies on the theory of collective compact operators which has been used to prove the convergence of finite element approximations for Maxwell’s equations in Chapter 4 of .
Let be a quasiuniform triangulation of into tetrahedrons of maximal diameter which matches with the interface , i.e., both triangulations for and are combined into a standard triangulation of the whole domain . For convenience, we denote by the restriction of in the domain . Let be the spaces of polynomials of maximal total degree and be the spaces of homogeneous polynomials of total degree exactly . We define the finite element space of
and the Nédélec -conforming and Raviart–Thomas -conforming finite element spaces:
where is a subspace of homogeneous vector polynomials of degree
In addition, we use and to denote the degree- -conforming and