 # Mathematical and numerical analysis of a nonlocal Drude model in nanoplasmonics

In this paper, we consider the frequency-domain Maxwell's equations coupled to a nonlocal Drude model which describes the nonlocal optical response in metallic nanostructures. We prove the existence and uniqueness of weak solutions to the coupled equations. A Galerkin finite element method based on the Raviart--Thomas and Nédélec elements is proposed to solve the equations and the associated error estimates are given. This is the first work on the mathematical and numerical analysis of this model. Numerical examples are presented to verify our theoretical analysis.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

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

 ∇×E=−∂tB,∇⋅B=0,∇×H=∂tD+J,∇⋅D=ρ. (1)

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

 B=μH,D=ϵ0ϵ∞E. (2)

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:

 ⎧⎨⎩∂tn+∇⋅(nv)=0,me(∂t+v⋅∇+γ)v=−e(E+v×B)−∇pn, (3)

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

 ρ=−en,J=−env. (4)

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

 n(x,t)≈n0+n1(x,t),v(x,t)≈v0+v1(x,t). (5)

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

 ∇pn≈meβ2∇nn0, (6)

where is an important parameter representing the nonlocality related to the Fermi velocity . Using the assumptions above, we get the linearized hydrodynamic equation

 ∂tv=−emeE−γv−β2∇nn0 (7)

and the linearized continuity equation

 ∂tn+n0∇⋅v=0. (8)

Differentiating (7) with respect to time , inserting the linearized current density and using (8), we obtain

 ∂ttJ+γ∂tJ−β2∇(∇⋅J)−ω2pε0∂tE=0, (9)

where is the plasma frequency. Combining (1), (2), and (9), we have Maxwell’s equations with the NHD model for metals

 (10)

Replacing with in (10

) by Fourier transformation in the time domain, where

is 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

 {∇×(μ−1∇×E)−ε0ε∞ω2E=iωJ,ω(ω+iγ)J+β2∇(∇⋅J)=iωω2pε0E. (11)
###### Remark 2.1

By Fourier transformation in the space domain, we replace with in the second equation of (11), which gives

 ω(ω+iγ)J−β2k2J=iωω2pε0E, (12)

and then we obtain the spatially-dispersive (relative) permittivity for the metal

 ϵ(ω,k)=ϵ∞−ω2pω(ω+iγ)−β2k2. (13)

The parameter represents the level of nonlocality. As , we recover the classical Drude permittivity .

In this paper we consider the following equations

 {∇×(μ−1∇×E)−εω2E=iωJ,inΩ,ω(ω+iγ)J+β2∇(∇⋅J)=iωω2pε0E,inΩs (14)

with the boundary conditions

 (15)

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

 μ={μ1,in Ωs,μ2,in Ω/Ωs,ε={ε1,in Ωs,ε2,in Ω/Ωs, (16)

and are positive constants.

###### Remark 2.2

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 

 n×E−H=n×Einc−Hinc,on∂Ω, (17)

where and represent the electromagnetic fields of the incoming light. By substituting into (17) and denoting by , we get the boundary conditions (15) for the electric field .

## 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 .

We define

 H(curl;Ω)={u∈L2(Ω)|∇×u∈L2(Ω)},H(div;Ω)={u∈L2(Ω)|∇⋅u∈L2(Ω)}, (18)

which are equipped with the norms

 ∥u∥H(curl;Ω)=∥u∥L2(Ω)+∥∇×u∥L2(Ω),∥u∥H(div;Ω)=∥u∥L2(Ω)+∥∇⋅u∥L2(Ω).

 HT(curl;Ω)={u∈H(curl;Ω)|uT=(n×u)×n∈L2(∂Ω)on∂Ω},H0(curl;Ω)={u∈H(curl;Ω)|u×n=0on∂Ω},H0(div;Ω)={u∈H(div;Ω)|u⋅n=0on∂Ω}. (19)

Functions in are equipped with the norm

 ∥u∥HT(curl;Ω)=∥u∥L2(Ω)+∥∇×u∥L2(Ω)+∥uT∥L2(∂Ω). (20)

For the sake of convenience, we denote by

 X(Ω)=HT(curl;Ω),∥u∥X(Ω)=∥u∥H(curl;Ω),Y(Ωs)=H0(div;Ωs),∥v∥Y(Ωs)=∥v∥H(div;Ωs). (21)

We now give the weak formulation for the problem (14)-(15). Given , find , such that the equations

 (22)

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.

 ¯¯¯J={J,in Ωs% ,0,in Ω/Ωs,E|Ωs=1ΩsE, (23)

where

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

 ∥E∥X(Ω)+∥J∥Y(Ωs)≤C∥g∥L2(∂Ω), (24)

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

 −∥ET∥2L2(∂Ω)=Re(¯¯¯J,E). (25)

Next by setting in the second equation of (22) and taking the imaginary part of the equation, we have

 γω2pε0∥J∥2L2(Ωs)=Re(E|Ωs,J)s. (26)

Since , from (25) and (26), we deduce that

 J=0,ET=0. (27)

Thus we find that satisfies

 (μ−1∇×E,∇×u)−ω2(εE,u)−iω⟨ET,uT⟩=0,forallu∈X(Ω). (28)

It was proved in Theorem 4.12 of  that the homogeneous problem (28) exists the only solution . Consequently, is the only solution of (23) with .

Before proving the existence of solutions of (22), we give two useful lemmas. We have the following Helmholtz decompositions for and

 X(Ω)=X0(Ω)⊕∇H10(Ω),Y(Ωs)=Y0(Ωs)⊕∇×H0(curl;Ωs), (29)

where

 X0(Ω)={u∈X(Ω)|(εu,∇ξ)=0forallξ∈H10(Ω)},Y0(Ωs)={v∈Y(Ωs)|(v,∇×w)s=0forallw∈H0(curl;Ωs)}. (30)

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

 v=v0+∇×A, (31)

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

 E=E0+∇φ,J=J0+∇×A (32)

for some and . In addition, we assume that satisfies .

Substituting (32) into (22), we find that

 (33)

for all .

Now taking in (33), where and , we obtain

 −ω2(ε∇φ,∇ξ)=iω(¯¯¯J0,∇ξ),forallξ∈H10(Ω), (34)

and

 (35)

where we have used the fact that and . By introducing a Lagrangian multiplier , we can rewrite (35) as

 (36)

Next we take in (33) to obtain

 {(μ−1∇×E0,∇×u)−ω2(εE0,u)−iω⟨E0,T,uT⟩=⟨g,uT⟩+iω(¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯J0+∇×A,u),−ω(ω+iγ)(J0,v)s+β2(∇⋅J0,∇⋅v)s=−iωω2pε0((E0+∇φ)|Ωs,v)s. (37)

Combining (33)-(37), we now reformulate the problem (22) as follows. Find , such that

 ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩(μ−1∇×E0,∇×u)−ω2(εE0,u)−iω⟨E0,T,uT⟩+β2(∇⋅J0,∇⋅v)s−ω(ω+iγ)(J0,v)s+iωω2pε0(E0|Ωs,v)s−iω(¯¯¯J0,u)=iω(¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∇×A,u)−iωω2pε0((∇φ)|Ωs,v)s+⟨g,uT⟩ (38)

for all , where and satisfy (34) and (36), respectively.

We now define the sesquilinear form by

 (39)

for all . For convenience, we introduce the following notations.

 ∥^u∥X(Ω)×Y(Ωs):=∥u1∥X(Ω)+∥u2∥Y(Ωs),∥^u∥L2(Ω)×L2(Ωs):=∥u1∥L2(Ω)+∥u2∥L2(Ωs).

It is not difficult to prove the following lemma. There exists a constant depending on , , , , , , and such that

 |a+(^u,^u)|≥α∥^u∥X(Ω)×Y(Ωs),forall^u=(u1,u2)∈X0(Ω)×Y0(Ωs). (40)

To proceed further, we define a map : such that if then satisfies

 a+(K^u,^v)=−2ω2(εu1,v1)−2ω2(u2,v2)s−iω(¯¯¯¯u2,v1)+iωω2pε0(u1|Ωs,v2)s−iω(¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∇×A,v1)+iωω2pε0((∇φ)|Ωs,v2)s,forall^v=(v1,v2)∈X0(Ω)×Y0(Ωs), (41)

where and are the solutions of

 −ω2(ε∇φ,∇ξ)=iω(¯¯¯¯u2,∇ξ),forallξ∈H10(Ω), (42)

and

 (43)

We have the following result. The operator is a bounded and compact map from into . Moreover,

 ∥K^u∥X(Ω)×Y(Ωs)≤C∥^u∥L2(Ω)×L2(Ωs). (44)

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

 |a+(^u,^v)|≤C∥^u∥X(Ω)×Y(Ωs)∥^v∥X(Ω)×Y(Ωs) (45)

holds for all . Coercivity of has been given in Lemma 3. It remains to show that there exists a constant such that

 |−2ω2(εu1,v1)−2ω2(u2,v2)s−iω(¯¯¯¯u2,v1)+iωω2pε0(u1|Ωs,v2)s−iω(¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∇×A,v1)+iωω2pε0((∇φ)|Ωs,v2)s|≤C∥^u∥L2(Ω)×L2(Ωs)∥^v∥L2(Ω)×L2(Ωs), (46)

where and are the solutions of (42) and (43), respectively. To prove (46), it suffices to deduce that

 |((∇φ)|Ωs,v2)s|+|(¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯∇×A,v1)|≤C∥^u∥L2(Ω)×L2(Ωs)∥^v∥L2(Ω)×L2(Ωs). (47)

Since satisfies (42), it is easy to see that

 ∥∇φ∥L2(Ω)≤C∥u2∥L2(Ωs). (48)

By using the classical theory of variational problems , we see that the mixed problem (43) exists a unique solution . In addition,

 ∥∇×A∥L2(Ωs)≤C∥u1∥L2(Ω). (49)

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

 a+(^F,^v)=⟨g,v1,T⟩,forall^v=(v1,v2)∈X0(Ω)×Y0(Ωs). (50)

By using the Lax–Milgram theorem again, we see that is well defined and

 ∥^F∥X(Ω)×Y(Ωs)≤C∥g∥L2(∂Ω). (51)

By virtue of the operator , we find that the problem (38) is equivalent to finding such that

 (I+K)^Z=^F. (52)

Since is compact, by applying the Fredholm alternative theorem and Lemma 3, we see that (52) exists a unique solution with the following estimate

 ∥^Z∥L2(Ω)×L2(Ωs)≤C∥^F∥L2(Ω)×L2(Ωs). (53)

Note that (52) implies that , from which we deduce

 ∥^Z∥X(Ω)×Y(Ωs)≤C(∥^F∥X(Ω)×Y(Ωs)+∥K^Z∥X(Ω)×Y(Ωs))≤C(∥^F∥X(Ω)×Y(Ωs)+∥^Z∥L2(Ω)×L2(Ωs)), (54)

where we have used (44). Substituting (53) into (54) and applying (51), we arrive at

 ∥E0∥X(Ω)+∥J0∥Y(Ωs)=∥^Z∥X(Ω)×Y(Ωs)≤C∥g∥L2(∂Ω). (55)

Since and satisfy (34) and (36), respectively, it follow from (55) that

 ∥∇φ∥L2(Ω)+∥∇×A∥L2(Ωs)≤∥E0∥L2(Ω)+∥J0∥L2(Ωs)≤C∥g∥L2(∂Ω). (56)

Combining (55)-(56) and recalling the Helmholtz decompositions for and in (32), we have (24) and complete the proof of Theorem 3.

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

 |(J,∇×w)s|≤C|(∇×E|Ωs,w)s|≤C∥∇×E∥L2(Ω)∥w∥L2(Ωs)≤C∥w∥L2(Ωs), (57)

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 .

If is constant on the whole domain , selecting with in the first equation of (22) and applying a similar argument, we find

 |(E,∇ξ)|≤C∥ξ∥L2(Ω). (58)

Consequently, we have , which yields that (Theorem 3.47 of ).

## 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

 Sh={uh∈H10(Ω):uh|K∈Pr,∀K∈Th}, (59)

and the Nédélec -conforming and Raviart–Thomas -conforming finite element spaces:

 Xh={uh∈HT(curl;Ω):uh|K∈(Pr−1)3⊕Rr,∀K∈Th},Yh,s={uh∈H0(div;Ωs):uh|K∈(Pr−1)3⊕˜Pr−1x,∀K∈Th,s}, (60)

where is a subspace of homogeneous vector polynomials of degree

 Rr={u∈(˜Pr)3:x⋅u=0}.

In addition, we use and to denote the degree- -conforming and