# A hybridizable discontinuous Galerkin method for solving nonlocal optical response models

###### Abstract

We propose Hybridizable Discontinuous Galerkin (HDG) methods for solving the frequency-domain Maxwell’s equations coupled to the Nonlocal Hydrodynamic Drude (NHD) and Generalized Nonlocal Optical Response (GNOR) models, which are employed to describe the optical properties of nano-plasmonic scatterers and waveguides. Brief derivations for both the NHD model and the GNOR model are presented. The formulations of the HDG method are given, in which we introduce two hybrid variables living only on the skeleton of the mesh. The local field solutions are expressed in terms of the hybrid variables in each element. Two conservativity conditions are globally enforced to make the problem solvable and to guarantee the continuity of the tangential component of the electric field and the normal component of the current density. Numerical results show that the proposed HDG methods converge at optimal rate. We benchmark our implementation and demonstrate that the HDG method has the potential to solve complex nanophotonic problems.

###### keywords:

Maxwell’s equations; nonlocal hydrodynamic Drude model; general nonlocal optical response theory; hybridizable discontinuous Galerkin method## 1 Introduction

Nanophotonics is the active research field field concerned with the study of interactions between nanometer scale structures/media and light, including near-infrared, visible, and ultraviolet light. It bridges the micro and the macro worlds, and there are many connections between theoretical studies and feasible engineering. The many fascinating (potential) applications include invisibility cloaking, nano antennas, metamaterials, novel biological detection and treatment technologies, as well as new storage media Maier ().

All of the above applications of nanophotonics require elaborate control of the propagation of light waves. In order to do so, appropriate mathematical models are needed to predict the behavior of light-matter interactions. Metals are interesting for nanophotonics because they can both enhance and confine optical fields, making plasmonics of interest to emerging quantum technologies Tame (); Bozhevolnyi (); Fitzgerald (). This is enabled by the existence of Surface Plasmons (SPs). SPs are coherent oscillations that exist as evanescent waves at both sides of the interface between any two materials where the real part of the dielectric function changes sign across the interface. The typical example is a metal-dielectric interface, such as a metal sheet in air SaridChallener (). Maxwell’s equations can be employed to model the macroscale electromagnetic waves and armed with classical electrodynamics there are numerous approaches ranging from classical electrodynamics to ab initio treatments Gallinet2015 (); Varas2016 (). Ab initio techniques can be used to simulate the microscopic dynamics on the atomic scale, but with ab initio methods one can only deal with systems with up to about ten thousand atoms Varas2016 (), thus calling for semiclassical treatments MortensenGNOR (); ToscanoEtal2015 () or more effective inclusions of quantum phenomena into classical electrodynamics Luo2013 (); Yan2015 (); Christensen2016 (); Zhu2016 ().

If one models the interaction of light with metallic nanostructures classically or semiclassically, then this calls for appropriate modelling of the material response as described for example by the Drude model Drude1900 (); DresselScheffler2006 (), the Nonlocal Hydrodynamic Drude (NHD) model Bloch1933a (); RazaEtal2015JPCM (), or the Generalized Nonlocal Optical Response (GNOR) theory MortensenGNOR (), all in combination with and coupled to Maxwell’s equations. Except for some highly symmetric geometries, analytical solutions to the resulting systems of differential equations are not available. Thus, numerical treatment of these systems of PDEs is an important aspect of nanophotonics research. Numerical experiments help to find promising systems and geometries before real fabrication, to obtain optimized parameters, to visualize field distributions, to investigate the dominant contribution to a phenomenon, to explain experimental observations, and so on Busch2011 ().

Several numerical methods exist for computing the solution of Maxwell’s equations Gallinet2015 (). For time-dependent problems, the Finite-Difference Time-Domain (FDTD) algorithm is the most popular method Taflove2005 () among physicists and engineers. More recently, the Discontinuous Galerkin Time-Domain (DGTD) method has drawn a lot of attention because of several appealing features, for example, easy adaptation to complex geometries and material composition, high-order accuracy, and natural parallelism HesthavenWarburton2002 (). For time-harmonic problems, the Finite Element Method (FEM) is most widely used for the solution of Maxwell’s equations. In very recent years, the Hybridizable Discontinuous Galerkin (HDG) method appears as a promising numerical method for time-harmonic problems because it inherits nearly all the advantages of the DG methods while leading to a computational complexity similar to FEM CGL2009 (); NPC2011 (); LLP2014 (); HLLH2016 ().

Currently, FDTD (for time-dependent problems) and FEM (for time-harmonic problems) methods are still the methods most commonly adopted for the simulation of light-matter interactions. Most often, commercial simulation
software (such as Lumerical FDTD^{1}^{1}1https://www.lumerical.com/
and Comsol Multiphysics^{2}^{2}2http://www.comsol.com/) is used for that purpose. However, these methods and computer codes do not always offer the required capabilities for addressing accurately and efficiently the complexity of the physical phenomena underlying nanometer scale light-matter interactions. In the academic community, also the DGTD method has recently been considered in this context Busch2011 (); Viquerat2015 (); HuangYQ2016 (). In Ref. Schmittetal2016 (), some numerical results are presented for the NHD model using the DGTD method. In the present paper we are employing the HDG method to solve the frequency-domain NHD and GNOR models. The development of accurate and efficient numerical methods for computational nanophotonics is expected to be a long-lasting demand, both because new models are regularly proposed that require innovative numerical methods, and because there is demand for more accurate and faster simulation methods for existing models.

This paper introduces a HDG method for the solution the NHD and GNOR models. The rest of the paper is organized as follows. In section 2, we briefly introduce mathematical aspects both of the NHD model and of the GNOR model. HDG formulations are given in section 3. Numerical results are presented in section 4 to show the effectiveness of high-order HDG methods for solving problems in nanophotonics. We draw conclusions in section 5.

## 2 Physics problem: nonlocal optical response by nanoparticles

The problem considered is shown in Figure 1 where the nanometer-size metal is illuminated by an incident plane wave of light. The infinite scattering domain is truncated as a finite computational domain by employing an artificial absorbing boundary condition, which is designed to absorb outgoing waves.

### 2.1 Nonlocal hydrodynamic Drude model

There are a number of theories for the modeling of the light-matter interactions which are used under different settings. In this subsection, we briefly introduce the NHD model. The incoming light is described as a propagating electromagnetic wave that satisfies Maxwell’s equations. Without external charge and current, Maxwell’s equations of macroscopic electromagnetism for non-magnetic materials can be written as

(1) |

where and are respectively the magnetic and electric fields, is the permittivity constant, is the permeability constant, is introduced to account for the local response, and is the nonlocal hydrodynamic polarization current density which is due to the nonlocal material on the plasmonic scatterers HiremathEtal2012 (). In this paper, we will for simplicity set and , thereby focusing solely on the free-electron response to light. Equations (1) need to be completed to solve electromagnetic fields and because of the unknown polarization current density . The models that we will consider in this paper differ only in the assumed dynamics of the polarization current density, which we will now discuss in more detail.

The polarization current density due to the motion of the free-electron gas can be written as

(2) |

where is the charge of the electron, is the density of the electron gas (a scalar field), and is its hydrodynamic velocity (a vector field). Within the hydrodynamic model, the dynamics of the velocity field is given by Schmittetal2016 (); RazaEtal2015 ()

(3) |

where is the mass of an electron, is the Lorentz force with being the magnetic flux density, is a damping constant, is an energy functional of the fluid, and the term denotes the quantum pressure. Complementary to Eq. 3, the dynamics of the free-electron density is given by

(4) |

which is the well-known continuity relation that relates the velocity and the density .

The hydrodynamic dynamics described by Eq. (3) is obviously nonlinear in , but in the following we only consider the linear response of the electron gas on external fields. One can write a perturbation expansion and similarly for the electric and magnetic fields and for the density. Since in the absence of an external field , both the nonlinear term and the magnetic induction field disappear due to the linearization ToscanoEtal2015 (). If we furthermore assume the energy functional to be of the Thomas-Fermi form, then we obtain for the linearized quantum pressure

(5) |

where with being the Fermi velocity. The zero-order (i.e. equilibrium) density is constant within the plasmonic medium ToscanoEtal2015 (). Here in Eq. (5) and below, we write for the linearized density and similarly we will from now on simply write for the linearized velocity . As a result, we obtain the linearized hydrodynamic equation RazaEtal2015JPCM (); Schmittetal2016 ()

(6) |

as well as the linearized continuity relation

(7) |

Inserting Eqs. (2) (linearized as ) and (7) into (6), and taking the time-derivative , we obtain

(8) |

where is the plasma frequency with . By Fourier transformation we replace with , where is the imaginary unit and is the angular frequency, and obtain the frequency-domain relation between polarization current density and the electric field within the hydrodynamic model as

(9) |

This equation describes electron-field interaction within the plasmonic nanostructure . We will neglect spill-out of electrons outside the classical geometric surface of the structure, which for our purposes is a good assumption for noble metals such as silver and gold ToscanoEtal2015 (). Mathematically, this is arranged by imposing a hard-wall condition on the boundary , namely on Jewsbury:1981a (); Yan:2013a ().

### 2.2 General nonlocal optical response model

We also briefly present the mathematical derivation of the central equations of the GNOR model, based on Ref. MortensenGNOR (). In the GNOR model, also diffusion of the electron gas is taken into consideration. Let the density , where the last term is the induced density variation caused by a non-vanishing electric field , which we assume sufficiently small that holds. Instead of (7), we now consider the linearized convection-diffusion equation MortensenGNOR ()

(10) |

where is the diffusion constant for the charge-carrier diffusion. Then the current density is given by Fick’s law

(11) |

Multiplying (6) by the charge of the electron , the equilibrium density and taking the time-derivative we have

(12) |

Dividing (12) by and combining with Fick’s law (11) results in

(13) |

From the convection-diffusion equation (10), we have

(14) |

Like what we did for (8), transforming (14) to the frequency domain gives

(15) |

The physical predictions obtained by the GNOR and NHD models often differ substantially, as illustrated below. However, from a computational point of view the GNOR model only differs by the replacement in the frequency domain, whereby the nonlocal hydrodynamic parameter acquires an often non-negligible imaginary part. In the GNOR model we have the same additional boundary condition on as in the NHD model.

### 2.3 Specification to 2D TM mode

Now we can couple Maxwell’s equation (1) with (9) for the NHD model, or similarly with (15) for the GNOR model. We will compute light extinction by infinitely long nanowires. We take the wire axes along the -direction and consider TM-polarized incident light, i.e. polarized in the -plane. In this 2D setting, we can define to be a vector and a scalar function. Coupling the time-harmonic Maxwell’s equations and hydrodynamic Drude model (9), we have in 2D

(16) |

If the Silver-Müller boundary condition (first-order absorbing boundary condition) Stupfel1994 () is applied on the boundary of the computational domain, then we have the boundary conditions

(17) |

where and stand for the electromagnetic fields of the incoming light.

## 3 HDG formulations of nonlocal optical response models

### 3.1 The promise of hybridizable DG methods

In the Introduction some properties and advantages of DG and HDG methods were briefly mentioned, which we here explain in more detail. The classic DG method is seldomly employed for solving stationary problems, because it duplicates degrees of freedom (DOFs) on every internal edge. Thus the number of globally coupled DOFs is much greater than the number of DOFs required by conforming finite element methods for the same accuracy. Consequently, DG methods are expensive in terms of both CPU time and memory consumption. Hybridization of DG methods CGL2009 () is devoted to addressing this issue while at the same time keeping all the advantages of DG methods. HDG methods introduce additional hybrid variables on the edges of the elements. Then we define the numerical traces arising from partial integration in the DG formulations through the hybrid variables. We can thus define the local (element-wise) solutions by hybrid variables. Conservativity conditions are imposed on numerical traces to ensure the continuity of the tangential component of the electric field and the normal component of the current density and to make the problem solvable. As a result, HDG methods produce a linear system in terms of the DOFs of the additional hybrid variables only. In this way, the number of globally coupled DOFs is greatly reduced as compared to the classic DG method. In a recent study Yakovlev:2016 (), the authors showed that HDG methods outperform FEM in many cases.

### 3.2 Computational concepts and notations

In order to give a clear presentation of the HDG method, here we introduce some computational concepts and notations. We divide the computational domain into triangle elements. The union of all the triangles is denoted by . By we denote the union of all edges of . Furthermore, stands for the union of all the edges associated with the nanostructure. For an edge associated with two elements , let be the traces of on from the interior of , see Fig. 2, where we use the term trace to denote the restriction of a function on the boundaries of the elements Arnold2002 (). Note that from now on is used to describe a general vector function instead of velocity.

On every face, we define mean (average) values and jumps as

where denotes the outward unit norm vector to and denotes the unit tangent vectors to the boundaries such that and . For the boundary edges, either on or on , these expressions are modified as

Let denote the space of polynomial functions of degree at most on a domain . For any element , let be the space and the space . The discontinuous finite element spaces are then defined by

(18) | ||||

where is the space of square integrable functions on the domain . We also introduce a traced finite element space

Note that consists of functions which are continuous on an edge, but discontinuous at its ends. The restrictions of , and in are denoted by , and . For two vectorial functions and in , we introduce the inner product , where denotes the complex conjugation. Likewise for scalar functions and in , the inner product is defined as provided is a domain in . Finally we define the edge overlap , where is a specific edge. Accordingly, we can define the total edge overlap for the whole triangulation or for relevant subsets of edges. Important cases are

denoting, respectively, the total edge overlap on the computational domain, the cumulative edge overlap on the absorbing boundary of the computational domain, and finally the cumulative edge overlap on the nanostructure.

### 3.3 DG formulation of the coupled electrodynamical equations

We begin the construction of a DG implementation of the hydrodynamic Drude model by rewriting the coupled electrodynamical equations (16) into a system of first-order equations

(19) |

where we introduced the scalar function which coincides with a scaled charge density. In general, a DG method seeks an approximate solution in the space that for each element (in our case: for each discretization triangle) satisfies ElBouajaji2013 ()

(20) |

The application of appropriate Green’s formulas to this system of equations leads to terms on the element boundaries Arnold2002 (). These boundary terms are the keys to connect the elements, since the elements themselves are independent due to the nature of the discontinuous finite elements spaces of Eq. (18). In a DG method, one replaces the boundary terms by so-called numerical traces and CGL2009 (); LLP2013 (), which are also known as ‘numerical fluxes’ in the literature Busch2011 (). These numerical traces are defined as

(21) |

In these definitions there is still freedom to choose values for the parameters, and this corresponds to different DG schemes: by setting , one obtains the centered flux DG scheme. With , one obtains the upwind flux DG scheme ElBouajaji2013 (). For more validated DG schemes, we refer the interested readers to Ref. Arnold2002 (). Having defined the numerical traces, we finally form a global system of linear equations involving all the DOFs on all the elements

(22) |

which are coupled equations that are valid whatever DG scheme is adopted.

### 3.4 Hybridizable DG implementation of the electrodynamical equations

In Sec. 3.1 we mentioned that hybridized DG methods have advantages as compared to the classic DG schemes, and here we discuss the hybridized approach in more detail. Unlike in the above classic DG formulations where the numerical traces directly couple the values from the elements on both sides of the edges, in a HDG formulation the numerical traces are defined through hybrid variables. Introducing two hybrid variables and which live only on the boundaries of the elements, we define the numerical traces by

(23) |

where and are two stabilization parameters. Replacing the numerical traces in (22) with the expressions in (23) and applying Green’s formulas to the first and fourth equations in (22), we obtain the local formulation of the HDG method as

(24) |

One can solve the local fields element by element once the solutions for and are obtained. In order to make the problem solvable, we need to employ global conditions

(25) |

The first relation in (25) weakly enforces the continuity condition for the tangential component of the electric field across any edges, and also takes into account the Silver-Müller absorbing boundary condition. The other global condition in Eq. (25) weakly enforces the continuity condition for the normal component of the current density across any edges. The additional boundary condition on the surface of the nanostructure is implicitly contained in this relation.

Substituting and in (25) with the definitions in (23), we arrive at the global reduced system of equations

(26) |

Note that we used the fact that in (26). The two relations in Eq. (26) are not independent. They are coupled through the local solutions of , , and of the local equations (24).

Remark I. The proposed HDG formulation for the global system (26) is naturally consistent with the boundary conditions, both on the artificial boundary and on the medium boundary.

Remark II. Globally, we only need to solve Eq. (26), in which the fields , , and are replaced by the solutions in terms of and from the local problems (24). So the global DOFs are associated with in the whole computational domain, while they are associate with only within in the material medium. The discretization leads to a system of linear equations

(27) |

where and are vectors accounting for the degrees of freedom of the hybrid variables and respectively, and the coefficient matrix is large and sparse.

## 4 Numerical results

In this section we present numerical results to validate the proposed HDG formulations. All HDG methods have been implemented in Fortran 90. All our tests are performed on a Macbook with a 1.3 GHz Inter Core i5 CPU and 4 GB memory. We employ the multifrontal sparse direct solver MUMPS AmestoyEtal2000 () to solve the discretized systems of linear equations.

In HDG methods, we calculate the total fields and . The scattered fields are then calculated by subtracting the incident field from the total fields. We use HDG- to denote the HDG method with interpolation order . Here we choose fixed values for the stabilization parameters. Different choices are discussed in Ref. Gopalakrishnan2015 ().

### 4.1 Convergence study: Wave propagation in a cavity

While elsewhere in this article we focus on nanowire structures, here we first study the convergence of our method by considering wave propagation in a cavity. This cavity is assumed to be a square domain with the PEC boundary condition and hard-wall condition

This test case can be viewed as the frequency-domain version of the first test case in Schmittetal2016 (). The simplicity is achieved by introducing artificial current density and electric field, such that the analytical solutions coincide with Maxwell’s equations and with the hydrodynamic equation

(28) |

where is the wave number, with being the light speed. We make this modification to unify the scale of the electric and magnetic fields. The artificial terms and are also the analytical solution to this equation (28):

(29) | ||||

We only take the real part of and and the imaginary part of and into consideration. In order to have this analytical solution (29), one needs to set the length of the square and . The convergence history of the HDG method with interpolation order is given in Table LABEL:tbl:convRect and Figure 3. Mesh size is the edge length of elements associated to the boundary . The convergence orders are calculated by

where and denote a coarse and a refined mesh size, respectively. From Table LABEL:tbl:convRect and Figure 3, we observe that the proposed HDG method has an optimal convergence order which is for HDG-.

HDG- | HDG- | HDG- | ||||
---|---|---|---|---|---|---|

error | order | error | order | error | order | |

- | - | - | ||||

2.0 | 3.0 | 4.0 | ||||

2.0 | 3.0 | 4.0 | ||||

2.1 | 3.0 | 4.0 |

### 4.2 Benchmark problem: a cylindrical plasmonic nanowire

As our benchmark problem we consider the plasmonic behavior of a cylindrical nanowire. This has been used as a convenient benchmark problem for other numerical methods before ToscanoEtal2012 (); HiremathEtal2012 () because analytical solutions exist both for the local and for the NHD models, see the derivation in Ref. Ruppin2001 (). We make use of the fact that the analytical Mie solution of Ref. Ruppin2001 () allows making the nonlocal parameter complex-valued. This enables us to benchmark our HDG simulations against exact analytical results for the GNOR model as well. (For comparison, optical properties of a sphere in the GNOR model, based on exact Mie results, are discussed in Ref. RazaEtal2015JPCM ().)

For the NHD model, the configuration of the nanowire is taken to be the same as that in the first test in HiremathEtal2012 (): the radius of the cylinder is 2 nm, no interband transitions are considered, the plasma frequency , the damping constant , the Fermi velocity , and . For the GNOR model, we use the same parameters and furthermore we take MortensenGNOR (). An artificial absorbing boundary is set to be a concentric circle with a radius of nm.

As our benchmark observable we will calculate the Extinction Cross Section (ECS, ), which is given by the sum of the scattering cross section and the absorption cross section BergEtal2009 (),

More precisely, for the cylindrical nanowire we consider the extinction cross section per wire length, which actually has the units of a length. We scale this quantity by the diameter of the nanowire to obtain a dimensionless normalized extinction cross section that we denote by . It can be expressed as the sum of scaled scattering and absorption cross sections,

Here the integrations are performed along a closed path around the nanowire, and denotes the real part.

The simulations are performed on a mesh with 4,513 nodes, 8,896 elements and 13,280 edges of which 722 edges are located inside the nanostructure. The ECS is presented in Figure 4. Curvilinear treatment is employed for high-order accuracy, where the curved edges are geometrically approximated by second-order curves instead of straight lines LLP2013 (). From Figure 4 we can observe that the fourth-order HDG method produces an ECS curve that matches the analytical solution very well. By contrast, the first-order method is not accurate enough on this mesh. Contour plots of the electric field and the current density are presented in Figure 5. These results match well with corresponding results in Ref. HiremathEtal2012 () despite the lower resolution, probably because our simulation is performed on a coarser mesh. Comparing the two subfigures in Figure 4, we also find that the ECS curve for the GNOR model is smoother than for the NHD model. But this has a physical rather than a numerical origin. In particular the standing bulk plasmon resonances above the plasma frequency in the NHD model are essentially washed out by the introduced diffusion in the GNOR model. The ECS curves of HDG- and HDG- are are not presented in Figure 4, but we found that they they lie in between the displayed curves of HDG- and HDG-.

In our 2D simulations, we use a sparse direct solver MUMPS AmestoyEtal2000 () to solve the resulting systems of linear equations. We need to solve a linear system at each frequency. The computational performance mainly relies on the size of the coefficient matrices, i.e. the number of degrees of freedom (#DOF). The computational performance for one frequency is given in Table LABEL:tbl:perf, where denotes the CPU time for construction the matrices, denotes the CPU time used by MUMPS for the factorization of the coefficient matrix (27), and memory denotes the memory consumed by MUMPS. From Table LABEL:tbl:perf we can see that the HDG- is more expensive than HDG- in CPU time for both construction and factorization. However, high-order methods are preferable because they costs less for the same accuracy LLP2013 ().

#DOF | (second) | (second) | memory (MB) | |
---|---|---|---|---|

HDG- | 28,260 | 0.067 | 0.36 | 74 |

HDG- | 70,650 | 2.4 | 3.3 | 418 |

### 4.3 Dimer of cylindrical nanowires

Plasmonic dimer structures with small gaps are both experimentally interesting and computationally challenging because of high field enhancements in the gap region ToscanoEtal2012 (); Gallinet2015 (); RazaEtal2015JPCM (). Here we present our HDG simulations of a cylindrical gold dimer geometry as shown in Figure 6(a), and this particular configuration is from Ref. RazaEtal2015 (). A typical mesh is shown in Figure 6(b). On a mesh with 5,829 nodes, 11,520 triangles and 17,348 edges with 3,712 edges inside the nanostructure, we calculate the ECS curve by HDG-. The size of matrix for HDG- is , the matrix construction CPU time is 5.2 seconds, the factorization CPU time is 6.9 seconds for one frequency, and the memory cost is 717 MB.

For the material properties gold we use the same values as in Ref. ToscanoEtal2012 (): the plasma frequency , damping constant , the Fermi velocity , and the nonlocal parameter is determined by . The incoming plane wave of light is incident perpendicular to the line connecting the centers of the two circles, with a linear polarization parallel to this line (TM- or -polarization). A comparison of the ECS curves is presented in Figure 7. Overall, there are small but clear differences, illustrating that nonlocal response effects occur even for dimer structures for which the corresponding monomers ( nm nanowires) would show essentially no nonlocal effects ToscanoEtal2012 (). Both nonlocal models have blueshifted resonances as compared to the local model, and resonances in the GNOR model are less pronounced than in the local and NHD models. For smaller gap sizes, nonlocal blueshifts are larger and resonances are broadened more (the latter only in the GNOR model). Field distributions at the same particular frequency for the NHD and the GNOR models are compared in Figure 8. The figure illustrates the generic features that the GNOR model washes out some finer details of the field distributions, and also that minimal and maximal field values lie closer together in the GNOR model.

## 5 Conclusions

This paper introduces a HDG method to solve the nonlocal hydrodynamic Drude model and the GNOR model, both of which are often employed to describe light-matter interactions of nanostructures. The numerical fluxes are expressed in terms of two newly introduced hybrid terms. Only the hybrid unknowns are involved in the global problem. The local problems are solved element-by-element once the hybrid terms are obtained. The proposed HDG formulations naturally couple the hard-wall boundary condition. Numerical results indicate that the HDG method converges at the optimal rate. Our benchmark simulations for a cylindrical nanowire and our calculations for a dimer structure show that the HDG method is a promising method in nanophotonics. Building on these results, in the near future we plan to generalize our computations to 3D structures, and to introduce domain decomposition and model order reduction into nanophotonic computations.

## Acknowledgments

The first author was supported by the NSFC (11301057) and the Fundamental Research Funds for the Central Universities (ZYGX2014J082). N. A. M. and M. W. acknowledge support from the Danish Council for Independent Research (FNU 1323-00087). M. W. acknowledges support from the Villum Foundation via the VKR Centre of Excellence NATEC-II. The Center for Nanostructured Graphene is sponsored by the Danish National Research Foundation, Project DNRF103.

## References

## References

- (1) S.A. Maier, Plasmonics: Fundamentals and Applications, Springer, New York, 2007.
- (2) M.S. Tame, K.R. McEnery, S.K. Özdemir, J. Lee, S.A. Maier, M.S. Kim, Quantum Plasmonics, Nature Physics 9(6): 329-340, 2013.
- (3) S.I. Bozhevolnyi, N.A. Mortensen, Plasmonics for emerging quantum technologies, Nanophotonics 2016. doi: 10.1515/nanoph-2016-0179
- (4) J.M. Fitzgerald, P. Narang, R.V. Craster, S.A. Maier, V. Giannini, Quantum Plasmonics, Proceedings of the IEEE PP(99): 1-16, 2016. doi: 10.1109/JPROC.2016.2584860
- (5) D. Sarid, W.A. Challener, Modern Introduction to Surface Plasmons: Theory, Mathematica Modeling, and Applications, Cambridge University Press, 2010.
- (6) B. Gallinet, J. Butet, O. J. F. Martin, Numerical methods for nanophotonics: standard problems and future challenges Laser Photonics Rev. 9(6): 577–603, 2015.
- (7) A. Varas and P. García-González, J. Feist, F.J. García-Vidal, A. Rubio, Quantum plasmonics: from jellium models to ab initio calculations, Nanophotonics 5(3): 409-426, 2016.
- (8) N.A. Mortensen, S. Raza, M. Wubs, T. Søndergaard, S. I. Bozhevolnyi, A generalized non-local optical response theory for plasmonic nanostructures, Nature Communications, 5:3809, 2014.
- (9) G. Toscano, J. Straubel, A. Kwiatkowski, C. Rockstuhl, F. Evers, H. Xu, N.A. Mortensen, M. Wubs, Resonance shifts and spill-out effects in self-consistent hydrodynamic nanoplasmonics, Nature Communications, 6:7132, 2015.
- (10) Y. Luo, A.I. Fernandez-Dominguez, A. Wiener, S.A. Maier, J.B. Pendry, Surface plasmons and nonlocality: a simple model, Phys. Rev. Lett. 111(9): 093901, 2013.
- (11) W. Yan, M. Wubs, N.A. Mortensen, Projected dipole model for quantum plasmonics, Phys. Rev. Lett. 115(13): 137403. 2015.
- (12) T. Christensen, W. Yan, A.-P. Jauho, M. Soljačić, N.A. Mortensen, Quantum corrections in nanoplasmonics: shape, scale, and material, arXiv:1608.05421.
- (13) W. Zhu, R. Esteban, A.G. Borisov, J.J. Baumberg, P. Nordlander, H.J. Lezec, J. Aizpurua, K.B. Crozier, Quantum mechanical effects in plasmonic structures with subnanometre gaps, Nature Communications 7:11495, 2016.
- (14) P. Drude, Zur Elektronentheorie der Metalle, Annalen der Physik, 306(3):566-613, 1900.
- (15) M. Dressel, M. Scheffler, Verifying the Drude response. Ann. Phys., 15(7-8): 535-544, 2006.
- (16) F. Bloch, Bremsvermögen von Atomen mit mehreren Elektronen, Z. Physik 81:363, 1933.
- (17) S. Raza, S. I. Bozhevolnyi, M. Wubs, N.A. Mortensen, Nonlocal optical response in metallic nanostructures, J. Phys.: Condens. Matter, 27: 183204, 2015.
- (18) K. Busch, M. König, J. Niegemann, Discontinuous Galerkin methods in nanophotonics, Laser Photonics Rev., 5(6): 773-809, 2011.
- (19) A. Taflove, S.C. Hagness, Computational Electrodynamics - The Finite-Difference Time-Domain Method, Third Edition, Artech House Publishers, 2005.
- (20) J.S. Hesthaven, T. Warburton, Nodal high-order methods on unstructured grids: I. Time-domain solution of Maxwell’s equations, J. Comput. Phys., 181(1):186-221, 2002.
- (21) B. Cockburn, J. Gopalakrishnan, R. Lazarov, Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems, SIAM J. Numer. Anal., 47(2):1319-1365, 2009.
- (22) N.C. Nguyen, J. Peraire, B. Cockburn, Hybridizable discontinuous Galerkin methods for the time-harmonic Maxwell’s equations, J. Comput. Phys., 230(19): 7151-7175, 2011.
- (23) L. Li, S. Lanteri, R. Perrussel, A hybridizable discontinuous Galerkin method combined to a Schwarz algorithm for the solution of 3d time-harmonic Maxwell’s equation, J. Comput. Phys., 256(1): 563-581, 2014.
- (24) Y.X. He, L. Li, S. Lanteri, T.Z. Huang, Optimized Schwarz algorithms for solving time-harmonic Maxwell’s equations discretized by a hybridizable discontinuous Galerkin method, Comput. Phys. Commun., 200: 23-31, 2016.
- (25) J. Viquerat, Simulation of electromagnetic waves propagation in nano-optics with a high-order discontinuous Galerkin time-domain method, PhD thesis, University of Nice-Sophia Antipolis, December 2015.
- (26) Y. Q. Huang, J. C. Li, W. Yang, Theoretical and numerical analysis of a non-local dispersion model for light interaction with metallic nanostructures, Comput. Math. Appl., 72: 921-932, 2016.
- (27) N. Schmitt, C. Scheid, S. Lanteri, A. Moreau, 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(1): 396-415, 2016.
- (28) K.R. Hiremath, L. Zschiedrich, 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: 5890-5896, 2012.
- (29) S. Raza, M. Wubs, S. I. Bozhevolnyi, N. Asger Mortensen, Nonlocal study of ultimate plasmon hybridization, Opt. Lett., 40 (5): 839-842, 2015.
- (30) P. Jewsbury, Electrodynamic boundary conditions at metal interfaces, J. Phys. F: Met. Phys. 11(1): 195, 1981.
- (31) W. Yan, N. A. Mortensen, M. Wubs, Green’s function surface-integral method for nonlocal response of plasmonic nanowires in arbitrary dielectric environments, Phys. Rev. B 88: 155414, 2013.
- (32) B. Stupfel, Absorbing boundary conditions on arbitrary boundaries for the scalar and vector wave equations, IEEE T. Antenn. Propag., 42(6): 773-780, 1994.
- (33) S. Yakovlev, D. Moxey, R. M. Kirby, S. J. Sherwin, To CG or to HDG: A Comparative Study in 3D, J. Sci. Comput., 67(1): 192–220, 2016.
- (34) D.N. Arnold, B. Brezzi, B. Cockburn, L.D. Marini, Unified analysis of discontinuous Galerkin methods for elliptic problems, SIAM J. Numer. Anal., 39(5): 749-1779, 2002.
- (35) M. El Bouajaji and S. Lanteri, High order discontinuous Galerkin method for the solution of 2D time-harmonic Maxwell’s equations, Appl. Math. Comput., 219(13): 7241-7251, 2013.
- (36) L. Li, S. Lanteri, R. Perrussel, Numerical investigation of a high order hybridizable discontinuous Galerkin method for 2d time-harmonic Maxwell’s equations, COMPEL. 32: 1112-1138, 2013.
- (37) P. Amestoy, I. Duff, J. L’Excellent, Multifrontal parallel distributed symmetric and unsymmetric solvers, Comput. Methods Appl. Mech. Eng. 184: 501-520, 2000.
- (38) J. Gopalakrishnan, S. Lanteri, N. Olivares, and R. Perrussel, Stabilization in relation to wavenumber in HDG methods, Adv. Model. and Simul. in Eng. Sci. 2(13), 2015.
- (39) G. Toscano, S. Raza, A.P. Jauho, N.A. Mortensen, M. Wubs, Modified field enhancement and extinction by plasmonic nanowire dimers due to nonlocal response, Opt. Express, 20 (4): 4176-4188, 2012.
- (40) R. Ruppin, Extinction properties of thin metallic nanowires, Opt. Commun., 190: 205-209, 2001.
- (41) M.J. Berg, A. Chakrabarti, C.M. Sorensen, General derivation of the total electromagnetic cross sections for an arbitrary particle, J. Quant. Spectrosc. Radiat. Transfer, 10: 43-50, 2009.