---
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.'
address:
- 'School of Mathematical Sciences, University of Electronic Science and Technology of China, 611731, Chengdu, P.R. China'
- 'INRIA, 2004 Route des Lucioles, BP 93 06902 Sophia Antipolis Cedex, France'
- 'Department of Photonics Engineering, Technical University of Denmark, [Ø]{}rsteds Plads 343, DK-2800 Kgs. Lyngby, Denmark.'
- 'Center for Nanostructured Graphene, Technical University of Denmark, [Ø]{}rsteds Plads 343, DK-2800 Kgs. Lyngby, Denmark'
author:
- Liang Li
- Stéphane Lanteri
- 'N. Asger Mortensen'
- Martijn Wubs
title: A hybridizable discontinuous Galerkin method for solving nonlocal optical response models
---

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

Introduction {#sec:intro}
============

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] and Comsol Multiphysics[^2]) 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 \[sec:models\], we briefly introduce mathematical aspects both of the NHD model and of the GNOR model. HDG formulations are given in section \[sec:hdg\]. Numerical results are presented in section \[sec:tests\] to show the effectiveness of high-order HDG methods for solving problems in nanophotonics. We draw conclusions in section \[sec:con\].

Physics problem: nonlocal optical response by nanoparticles {#sec:models}
===========================================================

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

![Sketch of the incident electromagnetic wave illuminating the scatterer $\Omega_S$ that has a subwavelength size and is surrounded by free space. $\Omega_S$ is usually filled with metals, such as gold, silver or sodium. An artificial absorbing boundary $\partial\Omega$ is introduced to make a computational domain $\Omega$.[]{data-label="fig:struct"}](./struct.pdf)

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 $$\label{eq:max}
\left\{
\begin{array}{l}
\nabla\times {\mathbf{H}}=   \varepsilon_0\varepsilon_{\text{loc}}\partial_t{\mathbf{E}}+{\mathbf{J}},\\
\nabla\times{\mathbf{E}}= -\mu_0\partial_t\mathbf{H},
\end{array}
\right.$$ where $\mathbf{H}$ and ${\mathbf{E}}$ are respectively the magnetic and electric fields, $ \varepsilon_0$ is the permittivity constant, $\mu_0$ is the permeability constant, $\varepsilon_{\text{loc}}=\varepsilon_{\infty}+\varepsilon_{\text{inter}}$ is introduced to account for the local response, and ${\mathbf{J}}$ 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 $\varepsilon_{\text{inter}}=0$ and $\varepsilon_{\infty}=1$, thereby focusing solely on the free-electron response to light. Equations need to be completed to solve electromagnetic fields ${\mathbf{E}}$ and $\mathbf{H}$ because of the unknown polarization current density ${\mathbf{J}}$. 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 ${\mathbf{J}}$ due to the motion of the free-electron gas can be written as $$\label{eq:jHD}
{\mathbf{J}}= -en{\mathbf{v}},$$ where $e$ is the charge of the electron, $n$ is the density of the electron gas (a scalar field), and ${\mathbf{v}}$ is its hydrodynamic velocity (a vector field). Within the hydrodynamic model, the dynamics of the velocity field is given by [@Schmittetal2016; @RazaEtal2015] $$\label{eq:hd}
m_e(\partial_t+{\mathbf{v}}\cdot\nabla){\mathbf{v}}=-e({\mathbf{E}}+{\mathbf{v}}\times\mathbf{B}) - m_e\gamma{\mathbf{v}}-\nabla\left(\frac{\delta g[n]}{\delta n}\right),$$ where $m_e$ is the mass of an electron, $-e({\mathbf{E}}+{\mathbf{v}}\times\mathbf{B})$ is the Lorentz force with $\mathbf{B}$ being the magnetic flux density, $\gamma$ is a damping constant, $g[n]$ is an energy functional of the fluid, and the term $\nabla\left(\frac{\delta g[n]}{\delta n}\right)$ denotes the quantum pressure. Complementary to Eq. \[eq:hd\], the dynamics of the free-electron density is given by $$\partial_t n+\nabla\cdot(n{\mathbf{v}})=0,$$ which is the well-known continuity relation that relates the velocity ${\mathbf{v}}$ and the density $n$.

The hydrodynamic dynamics described by Eq. (\[eq:hd\]) is obviously nonlinear in ${\bf v}$, but in the following we only consider the linear response of the electron gas on external fields. One can write a perturbation expansion ${\bf v} \simeq {\bf v}_{0} + {\bf v}_{1}$ and similarly for the electric and magnetic fields and for the density. Since in the absence of an external field ${\bf v} = {\bf v}_{0} = {\bf 0}$, both the nonlinear term ${\mathbf{v}}\cdot\nabla{\mathbf{v}}$ and the magnetic induction field $\mathbf{B}$ 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 $$\label{eq:qpressure}
-\nabla\left(\frac{\delta g[n]}{\delta n}\right)= -m_e\beta^2\frac{1}{n_0}\nabla n,$$ where $\beta^2=\frac{3}{5}v_F^2$ with $v_F$ being the Fermi velocity. The zero-order (i.e. equilibrium) density $n_0$ is constant within the plasmonic medium [@ToscanoEtal2015]. Here in Eq.  and below, we write $n$ for the linearized density ${n}_{1}$ and similarly we will from now on simply write ${\bf v}$ for the linearized velocity ${\bf v}_{1}$. As a result, we obtain the linearized hydrodynamic equation [@RazaEtal2015JPCM; @Schmittetal2016] $$\label{eq:hds} 
m_e\partial_t{\mathbf{v}}=-e{\mathbf{E}}-m_e\gamma{\mathbf{v}}-m_e\beta^2\frac{1}{n_0}\nabla n,$$ as well as the linearized continuity relation $$\label{eq:nv}
\partial_t n = -n_0 \nabla\cdot {\mathbf{v}}.$$ Inserting Eqs.  (linearized as ${\mathbf{J}}=-en_0{\mathbf{v}}$) and into , and taking the time-derivative $\partial_t$, we obtain $$\label{eq:hdT} 
\partial_{tt} {\mathbf{J}}+\gamma\partial_t{\mathbf{J}}-\beta^2\nabla(\nabla\cdot{\mathbf{J}})-\omega_p^2\varepsilon_0\partial_t{\mathbf{E}}=0,$$ where $\omega_p$ is the plasma frequency with $\omega_p^2=n_0e^2/(m_e\varepsilon_0)$. By Fourier transformation we replace $\partial_t$ with $-{\mathrm{i}}\omega$, where ${\mathrm{i}}$ is the imaginary unit and $\omega$ is the angular frequency, and obtain the frequency-domain relation between polarization current density and the electric field within the hydrodynamic model as $$\label{eq:hdFre}
\omega(\omega+{\mathrm{i}}\gamma){\mathbf{J}}+\beta^2\nabla(\nabla\cdot{\mathbf{J}})={\mathrm{i}}\omega\omega^2_p\varepsilon_0{\mathbf{E}}.$$ This equation describes electron-field interaction within the plasmonic nanostructure $\Omega_S$. 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 $\partial\Omega_S$, namely ${\mathbf{n}}\cdot{\mathbf{J}}=0$ on $\partial\Omega_S$ [@Jewsbury:1981a; @Yan:2013a].

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 $n({\mathbf{r}},t)=n_0+n_1({\mathbf{r}},t)$, where the last term is the induced density variation caused by a non-vanishing electric field ${\mathbf{E}}$, which we assume sufficiently small that $n_1\ll n_0$ holds. Instead of , we now consider the linearized convection-diffusion equation [@MortensenGNOR] $$\label{eq:cd}
\partial_t (-e_1) n_1 = D\nabla^2(-e)n_1-\nabla\cdot[(-e)n_0{\mathbf{v}}]=-\nabla\cdot{\mathbf{J}},$$ where $D$ is the diffusion constant for the charge-carrier diffusion. Then the current density is given by Fick’s law $$\label{eq:current}
{\mathbf{J}}= (-e)n_0{\mathbf{v}}-D\nabla(-e)n_1.$$ Multiplying by the charge of the electron $-e$, the equilibrium density $n_0$ and taking the time-derivative we have $$\label{eq:hddv}
m_e(\partial_t+\gamma)\partial_t[(-e)n_0{\mathbf{v}}]=n_0e^2\partial_t{\mathbf{E}}-m_e\beta^2\nabla[\partial_t(-e)n_1].$$ Dividing by $m_e$ and combining with Fick’s law  results in $$\label{eq:hddj}
(\partial_t+\gamma)\{\partial_t{\mathbf{J}}+D\nabla[\partial_t(-e)n_1]\}=\frac{n_0e^2}{m_e}\partial_t{\mathbf{E}}-\beta^2\nabla[\partial_t(-e)n_1].$$ From the convection-diffusion equation , we have $$\label{eq:hdd}
(\partial_t+\gamma)[\partial_t{\mathbf{J}}+D\nabla(\nabla\cdot{\mathbf{J}})]=\frac{n_0e^2}{m_e}\partial_t{\mathbf{E}}-\beta^2\nabla(\nabla\cdot{\mathbf{J}}).$$ Like what we did for , transforming to the frequency domain gives $$\label{eq:gnor}
\omega(\omega+{\mathrm{i}}\gamma){\mathbf{J}}+[\beta^2+D(\gamma-{\mathrm{i}}\omega)]\nabla(\nabla\cdot{\mathbf{J}})={\mathrm{i}}\omega\omega^2_p\varepsilon_0{\mathbf{E}}.$$ 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 $\beta^2\rightarrow\beta^2+D(\gamma-{\mathrm{i}}\omega)$ 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 ${\mathbf{n}}\cdot{\mathbf{J}}=0$ on $\partial\Omega_S$ as in the NHD model.

Specification to 2D TM mode
---------------------------

Now we can couple Maxwell’s equation with for the NHD model, or similarly with for the GNOR model. We will compute light extinction by infinitely long nanowires. We take the wire axes along the $z$-direction and consider TM-polarized incident light, i.e. polarized in the $(x,y)$-plane. In this 2D setting, we can define ${\mathbf{E}}= (E_x, E_y)^T$ to be a vector and $H = H_z$ a scalar function. Coupling the time-harmonic Maxwell’s equations and hydrodynamic Drude model , we have in 2D $$\label{eq:hyd}
\left\{
\begin{aligned}
\nabla\times H & =  -{\mathrm{i}}\omega\varepsilon_0{\mathbf{E}}+{\mathbf{J}}, \text{ in } \Omega, \\
\nabla\times{\mathbf{E}}& = {\mathrm{i}}\omega\mu_0 H, \text{ in } \Omega,\\
\nabla(\nabla\cdot{\mathbf{J}}) + \frac{\omega(\omega+{\mathrm{i}}\gamma)}{\beta^2}{\mathbf{J}}& =\frac{{\mathrm{i}}\omega\omega^2_p\varepsilon_0}{\beta^2}{\mathbf{E}}, \text{ in } \Omega_S.
\end{aligned}
\right.$$ If the Silver-M[ü]{}ller boundary condition (first-order absorbing boundary condition) [@Stupfel1994] is applied on the boundary $\partial\Omega$ of the computational domain, then we have the boundary conditions $$\label{eq:bc}
\left\{
\begin{aligned}
{\mathbf{n}}\times{\mathbf{E}}-H & = {\mathbf{n}}\times{\mathbf{E}}^{\text{inc}}-H^{\text{inc}} = g^{\text{inc}}, \text{ on } \partial\Omega, \\
{\mathbf{n}}\cdot{\mathbf{J}}& =0, \text{ on } \partial\Omega_S,
\end{aligned}
\right.$$ where ${\mathbf{E}}^{\text{inc}}$ and $H^{\text{inc}}$ stand for the electromagnetic fields of the incoming light.

HDG formulations of nonlocal optical response models  {#sec:hdg}
=====================================================

The promise of hybridizable DG methods {#Sec:promiseHDG}
--------------------------------------

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.

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 $\Omega$ into triangle elements. The union of all the triangles is denoted by ${\mathcal{T}}_h$. By ${\mathcal{F}}_h$ we denote the union of all edges of ${\mathcal{T}}_h$. Furthermore, ${\mathcal{F}}^I_h$ stands for the union of all the edges associated with the nanostructure. For an edge associated with two elements $F =\overline{K^+}\cap\overline{K^-}\in{\mathcal{F}}_h$, let $({\mathbf{v}}^\pm, v^\pm)$ be the *traces* of $({\mathbf{v}}, v)$ on $F$ from the interior of $K^{\pm}$, see Fig. \[fig:interEdge\], 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 ${\mathbf{v}}$ is used to describe a general vector function instead of velocity.

![Two neighboring discretization elements (here: triangles) within the computational domain. An edge $F$ is shared by two elements $K^+$ and $K^-$. The outward normal vectors $n^+$ and $n^-$ point in opposite directions. A characteristic property of the DG method is that computed functions are allowed to be discontinuous across $F$ (hence the “D” in DG). For example, for a function $v$, be it a scalar or a vector, its value on $F$ from $K^+$ is $v^+$, while its value on $F$ from $K^-$ is $v^-$, and these $v^+$ and $v^-$ are not necessarily equal. By contrast, the hybrid variables in the HDG method [*are*]{} single-valued on $F$.[]{data-label="fig:interEdge"}](./traceDem.pdf)

On every face, we define *mean* (*average*) *values* ${\left\{ \cdot \right\}}$ and *jumps* ${\llbracket \cdot \rrbracket}$ as $$\left\{
    \begin{aligned}
      {\left\{ {\mathbf{v}}\right\}}_F  & = \dfrac{1}{2} ({\mathbf{v}}^+ + {\mathbf{v}}^-),
      \\
      {\left\{ v \right\}}_F &  = \dfrac{1}{2}(v^+ + v^-),
      \\
      {\llbracket {\mathbf{n}}\times{\mathbf{v}}\rrbracket}_F &  = {\mathbf{n}}^+\times{\mathbf{v}}^+ +
      {\mathbf{n}}^-\times{\mathbf{v}}^-, \\
       {\llbracket {\mathbf{n}}\cdot{\mathbf{v}}\rrbracket}_F &  = {\mathbf{n}}^+\cdot{\mathbf{v}}^+ +
      {\mathbf{n}}^-\cdot{\mathbf{v}}^-,
      \\
      {\llbracket v{\mathbf{t}}\rrbracket}_F &  = v^+ {\mathbf{t}}^+ + v^- {\mathbf{t}}^-,
    \end{aligned}
  \right.$$ where ${\mathbf{n}}^{\pm}$ denotes the outward unit norm vector to $K^\pm$ and ${\mathbf{t}}^{\pm}$ denotes the unit tangent vectors to the boundaries $\partial    K^{\pm}$ such that ${\mathbf{t}}^+\times{\mathbf{n}}^+=1$ and ${\mathbf{t}}^-\times{\mathbf{n}}^-=1$. For the boundary edges, either on $\partial\Omega$ or on $\partial\Omega_S$, these expressions are modified as $$\left\{
    \begin{aligned}
      {\left\{ {\mathbf{v}}\right\}}_F &  = {\mathbf{v}}^+,
      \\
      {\left\{ v \right\}}_F &  = v^+,
      \\
      {\llbracket {\mathbf{n}}\times{\mathbf{v}}\rrbracket}_F &  = {\mathbf{n}}^+ \times{\mathbf{v}}^+, \\
       {\llbracket {\mathbf{n}}\cdot{\mathbf{v}}\rrbracket}_F & = {\mathbf{n}}^+ \cdot{\mathbf{v}}^+,
      \\
      {\llbracket v{\mathbf{t}}\rrbracket}_F & = v^+ {\mathbf{t}}^+.
    \end{aligned}
  \right .$$ Let $\mathbb{P}_p(D)$ denote the space of polynomial functions of degree at most $p$ on a domain $D$. For any element $K\in{\mathcal{T}}_h$, let $V^p(K)$ be the space $\mathbb{P}_p(K)$ and ${\mathbf{V}}^p(K)$ the space $(\mathbb{P}_p(K))^2$. The discontinuous finite element spaces are then defined by $$\label{eq:dgSpace}
  \begin{aligned}
    V^p_h & = \left\{ v \in L^{2}(\Omega) \ \lvert \ v|_K\in
    V^p(K), \ \forall K\in{\mathcal{T}}_h\right\},
    \\
    {\mathbf{V}}^p_h & = \left\{ {\mathbf{v}}\in (L^2(\Omega))^2 \ \lvert \
    {\mathbf{v}}|_K\in {\mathbf{V}}^p(K), \ \forall K\in{\mathcal{T}}_h \right\},
  \end{aligned}$$ where $L^2(\Omega)$ is the space of square integrable functions on the domain $\Omega$. We also introduce a traced finite element space $$\label{SpaceM}
M^p_h = \left\{ \eta \in L^2({\mathcal{F}}_h) \ \lvert
  \ \eta|_F \in\mathbb{P}_p(F),
  \ \forall F\in{\mathcal{F}}_h  \right\}.$$ Note that $M^p_h$ consists of functions which are continuous on an edge, but discontinuous at its ends. The restrictions of $V^p_h$, ${\mathbf{V}}^p_h$ and $M^p_h$ in $\Omega_S$ are denoted by $\widetilde{V}^p_h$, $\widetilde{{\mathbf{V}}}^p_h$ and $\widetilde{M}^p_h$. For two vectorial functions $\mathbf{u}$ and ${\mathbf{v}}$ in $(L^2(D))^2$, we introduce the inner product $(\mathbf{u},{\mathbf{v}})_{D}   =   {\displaystyle}\int_D\mathbf{u}\cdot\overline{{\mathbf{v}}}  \,  {\rm d}x$, where $\overline{\cdot}$ denotes the complex conjugation. Likewise for scalar functions $u$ and $v$ in $L^2(D)$, the inner product is defined as $(u,v)_D = {\displaystyle}\int_D u  \overline v \,{\rm d}x$ provided $D$ is a domain in $\mathbb{R}^2$. Finally we define the edge overlap $\langle  u,v\rangle_F ={\displaystyle}\int_F u \overline v \, {\rm d}s$, where $F$ 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 $$\langle\cdot,\cdot\rangle_{{\mathcal{F}}_h}  =
    \sum_{F\in{\mathcal{F}}_h}\langle\cdot,\cdot\rangle_F, \qquad
     \langle\cdot,\cdot\rangle_{\partial\Omega}  =
    \sum_{F\in{\mathcal{F}}_h\cap\partial\Omega}\langle\cdot,\cdot\rangle_F,
     \qquad
    \langle\cdot,\cdot\rangle_{{\mathcal{F}}^I_h}  =
    \sum_{F\in{\mathcal{F}}^I_h}\langle\cdot,\cdot\rangle_F.$$ 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.

DG formulation of the coupled electrodynamical equations  {#subsec:dg}
---------------------------------------------------------

We begin the construction of a DG implementation of the hydrodynamic Drude model by rewriting the coupled electrodynamical equations  into a system of first-order equations $$\label{eq:MaxHyd}
\left\{
\begin{aligned}
{\mathrm{i}}\omega \varepsilon_0{\mathbf{E}}+ \nabla\times H -{\mathbf{J}}&= 0 \qquad \text{ in } \Omega, \\
{\mathrm{i}}\omega\mu_0 H - \nabla\times{\mathbf{E}}& =0 \qquad \text{ in } \Omega, \\
\nabla q +\frac{ \gamma-{\mathrm{i}}\omega}{\beta^2}{\mathbf{J}}- \frac{\omega_p^2\varepsilon_0}{\beta^2}{\mathbf{E}}& = 0 \qquad \text{ in } \Omega_S, \\
{\mathrm{i}}\omega q - \nabla\cdot{\mathbf{J}}& = 0 \qquad \text{ in } \Omega_S,
\end{aligned}
\right.$$ where we introduced the scalar function $q = ({\mathrm{i}}\omega)^{-1}\nabla\cdot{\mathbf{J}}$ which coincides with a scaled charge density. In general, a DG method seeks an approximate solution $({\mathbf{E}}_h, H_h, {\mathbf{J}}_h, q_h)$ in the space ${\mathbf{V}}^{p}_{h}\times V^{p}_{h}\times\widetilde{{\mathbf{V}}}^{p}_{h}\times \widetilde{V}^{p}_{h}$ that for each element $K$ (in our case: for each discretization triangle) satisfies [@ElBouajaji2013] $$\label{eq:MaxHydInt}
\left\{
\begin{aligned}
({\mathrm{i}}\omega \varepsilon_0{\mathbf{E}}_h, {\mathbf{v}})_{K} + (\nabla\times H_h, {\mathbf{v}})_{K} - ({\mathbf{J}}_h, {\mathbf{v}})_{K} & = 0 \qquad \forall {\mathbf{v}}\in {\mathbf{V}}^{p}(K), \\
({\mathrm{i}}\omega\mu_0 H_h, v)_{K} - (\nabla\times{\mathbf{E}}_h, v)_{K} & =0 \qquad \forall v\in V^{p}(K), \\
(\nabla q_h, {\mathbf{v}})_{K} + \Big(\frac{\gamma-{\mathrm{i}}\omega}{\beta^2}{\mathbf{J}}_h, {\mathbf{v}}\Big)_{K} - \Big(\frac{\omega_p^2\varepsilon_0}{\beta^2}{\mathbf{E}}_h, {\mathbf{v}}\Big)_{K} & = 0 \qquad \forall {\mathbf{v}}\in \widetilde{{\mathbf{V}}}^{p}(K), \\
({\mathrm{i}}\omega q_h, v)_{K} - (\nabla\cdot{\mathbf{J}}_h, v)_{K} & = 0 \qquad \forall v\in \widetilde{V}^{p}(K).
\end{aligned}
\right.$$ 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. . In a DG method, one replaces the boundary terms by so-called *numerical traces* $\hat{{\mathbf{E}}}_{h}, \hat{H}_h, \hat{{\mathbf{J}}}_{h}$ and $\hat{q}_{h}$ [@CGL2009; @LLP2013], which are also known as ‘numerical fluxes’ in the literature [@Busch2011]. These numerical traces are defined as $$\label{eq:dgNumericalTrace}
\left\{
\begin{aligned}
\hat{H}_h & ={\left\{ H_h \right\}} + \alpha_{E}{\llbracket {\mathbf{n}}\times{\mathbf{E}}_h \rrbracket}, \\
{\mathbf{n}}\times\hat{{\mathbf{E}}}_h & = {\left\{ {\mathbf{n}}\times{\mathbf{E}}_h \right\}}+\alpha_{H}{\llbracket H_h \rrbracket},\\
\hat{q}_h & = {\left\{ q_h \right\}}+\alpha_J{\llbracket {\mathbf{n}}\cdot{\mathbf{J}}_h \rrbracket},\\
{\mathbf{n}}\cdot\hat{{\mathbf{J}}}_h & ={\left\{ {\mathbf{n}}\cdot{\mathbf{J}}\right\}}+\alpha_q{\llbracket q_h \rrbracket}.
\end{aligned}
\right.$$ In these definitions there is still freedom to choose values for the $\alpha$ parameters, and this corresponds to different DG schemes: by setting $\alpha_{E}=\alpha_{H}=\alpha_{J}=\alpha_{q}=0$, one obtains the *centered flux* DG scheme. With $\alpha_{E}=\alpha_{H}=\alpha_{J}=\alpha_{q}=1$, 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 $$\label{eq:MaxHydIntGreen}
\left\{
\begin{aligned}
({\mathrm{i}}\omega \varepsilon_0{\mathbf{E}}_h, {\mathbf{v}})_{K} + (H_h, \nabla\times {\mathbf{v}})_{K} -  \langle\hat{H}_{h},{\mathbf{n}}\times{\mathbf{v}}\rangle_{\partial K} - ({\mathbf{J}}_h, {\mathbf{v}})_{K} & = 0 \quad  \forall {\mathbf{v}}\in {\mathbf{V}}^{p}(K), \\
({\mathrm{i}}\omega\mu_0 H_h, v)_{K} - ({\mathbf{E}}_h, \nabla\times v)_{K} - \langle{\mathbf{n}}\times\hat{{\mathbf{E}}}_{h}, v\rangle_{\partial K} & =0 \quad \forall v\in V^{p}(K), \\
-(q_h, \nabla\cdot{\mathbf{v}})_{K} + \langle \hat{q}_{h}, {\mathbf{n}}\cdot{\mathbf{v}}\rangle_{\partial K} +  \Big(\frac{\gamma-{\mathrm{i}}\omega}{\beta^2}{\mathbf{J}}_h, {\mathbf{v}}\Big)_{K} - \Big(\frac{\omega_p^2\varepsilon_0}{\beta^2}{\mathbf{E}}_h, {\mathbf{v}}\Big)_{K} & = 0 \quad \forall {\mathbf{v}}\in \widetilde{\mathbf{V}}^{p}(K), \\
({\mathrm{i}}\omega q_h, v)_{K} + ({\mathbf{J}}_h, \nabla v)_{K} -\langle{\mathbf{n}}\cdot\hat{{\mathbf{J}}}_{h}, v\rangle_{\partial K} & = 0 \quad \forall v\in\widetilde V^{p}(K),
\end{aligned}
\right.$$ which are coupled equations that are valid whatever DG scheme is adopted.

Hybridizable DG implementation of the electrodynamical equations  {#subsec:hdg}
-----------------------------------------------------------------

In Sec. \[Sec:promiseHDG\] 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 $\lambda_h$ and $\eta_h$ which live only on the boundaries of the elements, we define the numerical traces by $$\label{eq:numTrac}
\left\{
\begin{aligned}
\hat{H}_h & = \lambda_h, \\
\hat{{\mathbf{E}}}_h & = {\mathbf{E}}_h + \tau_{\lambda}(\lambda_h-H_h){\mathbf{t}}, \\
\hat{q}_h & = \eta_h, \\
\hat{{\mathbf{J}}}_h & = {\mathbf{J}}+ \tau_{\eta}(q_h-\eta_h){\mathbf{n}},
\end{aligned}
\right.$$ where $\tau_{\lambda}$ and $\tau_{\eta}$ are two stabilization parameters. Replacing the numerical traces in with the expressions in and applying Green’s formulas to the first and fourth equations in , we obtain the local formulation of the HDG method as $$\label{eq:MaxHydIntGreen2}
\left\{
\begin{aligned}
({\mathrm{i}}\omega \varepsilon_0{\mathbf{E}}_h, {\mathbf{v}})_{K} + (H_h, \nabla\times {\mathbf{v}})_{K} - \langle\lambda_{h},{\mathbf{n}}\times{\mathbf{v}}\rangle_{\partial K} - ({\mathbf{J}}_h, {\mathbf{v}})_{K} & = 0, \ \forall {\mathbf{v}}\in {\mathbf{V}}^{p}(K), \\
({\mathrm{i}}\omega\mu_0 H_h, v)_{K} - (\nabla\times{\mathbf{E}}_h, v)_{K} + \langle\tau_{\lambda}(H_h-\lambda_h), v\rangle_{\partial K} & =0, \ \forall v\in V^{p}(K), \\
-(q_h, \nabla\cdot{\mathbf{v}})_{K} + \langle\eta_{h}, {\mathbf{n}}\cdot{\mathbf{v}}\rangle_{\partial K} +  \Big(\frac{\gamma-{\mathrm{i}}\omega}{\beta^2}{\mathbf{J}}_h, {\mathbf{v}}\Big)_{K} - \Big(\frac{\omega_p^2\varepsilon_0}{\beta^2}{\mathbf{E}}_h, {\mathbf{v}}\Big)_{K} & = 0, \ \forall {\mathbf{v}}\in\widetilde {\mathbf{V}}^{p}(K), \\
({\mathrm{i}}\omega q_h, v)_{K} - (\nabla\cdot{\mathbf{J}}_h, v)_{K} -\langle\tau_{\eta}(q_h-\eta_h), v\rangle_{\partial K} & = 0, \ \forall v\in\widetilde V^{p}(K)  .
\end{aligned}
\right.$$ One can solve the local fields element by element once the solutions for $\lambda_h$ and $\eta_h$ are obtained. In order to make the problem solvable, we need to employ global conditions $$\label{eq:globalRaw}
\left\{
\begin{aligned}
\langle{\llbracket {\mathbf{n}}\times\hat{{\mathbf{E}}}_h \rrbracket}, v\rangle_{{\mathcal{F}}_h} - \langle\lambda_h, v\rangle_{\partial\Omega} & = \langle g^{\text{inc},v}\rangle_{\partial\Omega}, \ \forall v\in M^p_h, \\
\langle{\llbracket {\mathbf{n}}\cdot\hat{{\mathbf{J}}}_h \rrbracket}, v\rangle_{{\mathcal{F}}_h^I} & = 0, \ \forall v\in\widetilde M^p_h.
\end{aligned}
\right.$$ The first relation in  weakly enforces the continuity condition for the tangential component of the electric field across any edges, and also takes into account the Silver-M[" u]{}ller absorbing boundary condition. The other global condition in Eq.  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 $\hat{{\mathbf{E}}}_h$ and $\hat{{\mathbf{J}}}_h$ in with the definitions in , we arrive at the global reduced system of equations $$\label{eq:global}
\left\{
\begin{aligned}
\langle{\mathbf{n}}\times{\mathbf{E}}_h-\tau_{\lambda}(\lambda_h-H_h), v\rangle_{{\mathcal{F}}_h} - \langle\lambda_h, v\rangle_{\partial\Omega} & =  \langle g^{\text{inc}},v\rangle_{\partial\Omega}, \ \forall v\in M^p_h, \\
\langle{\mathbf{n}}\cdot{\mathbf{J}}_h+\tau_\eta(q_h-\eta_h), v\rangle_{{\mathcal{F}}_h^I} & = 0, \ \forall v\in\widetilde M^p_h.
\end{aligned}
\right.$$ Note that we used the fact that ${\mathbf{n}}\times{\mathbf{t}}=-1$ in . The two relations in Eq.  are not independent. They are coupled through the local solutions of ${\mathbf{E}}$, $H$, ${\mathbf{J}}$ and $q$ of the local equations .

**Remark I.** The proposed HDG formulation for the global system 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. , in which the fields ${\mathbf{E}}_h$, $H_h$, ${\mathbf{J}}_h$ and $q_h$ are replaced by the solutions in terms of $\lambda_h$ and $\eta_h$ from the local problems . So the global DOFs are associated with $\lambda_h$ in the whole computational domain, while they are associate with $\eta_h$ only within in the material medium. The discretization leads to a system of linear equations $$\label{eq:linSys}
A\begin{bmatrix}
\underline{\lambda}_h \\[5pt]
\underline{\eta}_h
\end{bmatrix}=\begin{bmatrix}
\underline{g}^\text{inc}_h \\[5pt]
0
\end{bmatrix},$$ where $\underline{\lambda}_h$ and $\underline{\eta}_h$ are vectors accounting for the degrees of freedom of the hybrid variables $\lambda$ and $\eta$ respectively, and the coefficient matrix $A$ is large and sparse.

Numerical results {#sec:tests}
=================

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 ${\mathbf{E}}^\text{tot}$ and $H^\text{tot}$. The scattered fields are then calculated by subtracting the incident field from the total fields. We use HDG-$\mathbb{P}_p$ to denote the HDG method with interpolation order $p$. Here we choose fixed values $\tau_{\lambda}= \tau_{\eta}=1$ for the stabilization parameters. Different choices are discussed in Ref. [@Gopalakrishnan2015].

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 $\Omega_{\Box}=\{(x,y)\in[0,L]\times[0,L]\}$ with the PEC boundary condition and hard-wall condition $${\mathbf{n}}\times {\mathbf{E}}= 0, \text{ and } {\mathbf{n}}\cdot{\mathbf{J}}=0, \text { on } \partial\Omega_{\Box}.$$ 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 $$\label{eq:MaxHydArt}
\left\{
\begin{aligned}
{\mathrm{i}}k{\mathbf{E}}+ \nabla\times H & = {\mathbf{J}}-{\mathbf{J}}^{a}, \\
{\mathrm{i}}k H - \nabla\times{\mathbf{E}}& =0, \\
\nabla q -\frac{ {\mathrm{i}}\omega}{\beta^2}{\mathbf{J}}&  = -\frac{\gamma}{\beta^2}{\mathbf{J}}+\frac{\gamma}{\beta^2}{\mathbf{J}}^{a} + \frac{\omega_p^2\varepsilon_0}{\beta^2}{\mathbf{E}}-\frac{\omega_p^2\varepsilon_0}{\beta^2}{\mathbf{E}}^{a}, \\
{\mathrm{i}}\omega q - \nabla\cdot{\mathbf{J}}&  = 0,
\end{aligned}
\right.$$ where $k=\frac{\omega}{c}$ is the wave number, with $c$ being the light speed. We make this modification to unify the scale of the electric and magnetic fields. The artificial terms ${\mathbf{J}}^a$ and ${\mathbf{E}}^a$ are also the analytical solution to this equation : $$\label{eq:anaSol}
\begin{aligned}
{\mathbf{E}}^a & = \frac{\sqrt{2}}{2}{\mathrm{i}}\begin{bmatrix}
-\cos(\frac{\sqrt{2}}{2}kx)\sin(\frac{\sqrt{2}}{2}ky)\\[5pt]
\sin(\frac{\sqrt{2}}{2}kx)\cos(\frac{\sqrt{2}}{2}ky)
\end{bmatrix},\\{\mathbf{J}}^a & = -\frac{\sqrt{2}\mu_0 k\beta^2}{2\omega}\begin{bmatrix}
\sin(\frac{\sqrt{2}}{2}kx)\cos(\frac{\sqrt{2}}{2}ky)\\[5pt]
\cos(\frac{\sqrt{2}}{2}kx)\sin(\frac{\sqrt{2}}{2}ky)
\end{bmatrix} .
\end{aligned}$$ We only take the real part of $H$ and ${\mathbf{J}}$ and the imaginary part of ${\mathbf{E}}$ and $q$ into consideration. In order to have this analytical solution , one needs to set the length of the square $L=\frac{\sqrt{2}\pi}{k}$ and $\beta^2=\frac{\omega^2}{k^2}$. The convergence history of the HDG method with interpolation order $\mathbb{P}_p \ (p=1,2,3)$ is given in Table \[tbl:convRect\] and Figure \[fig:convRect\]. Mesh size $h$ is the edge length of elements associated to the boundary $\partial\Omega_{\Box}$. The convergence orders are calculated by $$\frac{\log(\|{\mathbf{E}}^a-{\mathbf{E}}^h\|_{\Omega_{\Box}}^{h_2}/\|{\mathbf{E}}^a-{\mathbf{E}}^h\|_{\Omega_{\Box}}^{h_1})}{\log(h_2/h_1)},$$ where $h_1$ and $h_2$ denote a coarse and a refined mesh size, respectively. From Table \[tbl:convRect\] and Figure \[fig:convRect\], we observe that the proposed HDG method has an optimal convergence order which is $p+1$ for HDG-$\mathbb{P}_p$.

    \[tbl:convRect\]                                                                                        
  -------------------- ---------------------- ------- ---------------------- ------- ---------------------- -------
                                                                                                            
          $h$                  error           order          error           order          error           order
    $5\time10^{-8}$     $1.67\times10^{-9}$      -     $4.52\times10^{-10}$     -     $2.04\times10^{-11}$     -
   $2.5\time10^{-8}$    $4.10\times10^{-10}$    2.0    $5.61\times10^{-11}$    3.0    $1.28\times10^{-12}$    4.0
   $1.25\time10^{-8}$   $9.98\times10^{-11}$    2.0    $7.52\times10^{-12}$    3.0    $7.78\times10^{-14}$    4.0
   $6.25\time10^{-9}$   $2.40\times10^{-12}$    2.1    $9.11\times10^{-13}$    3.0    $5.03\times10^{-15}$    4.0

  : Convergence results for the cavity problem.

![Convergence history of the proposed HDG method for the cavity problem.[]{data-label="fig:convRect"}](./convRect.pdf)

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 $\beta$ 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 2nm, no interband transitions are considered, the plasma frequency $\omega_p=8.65\times10^{15}$, the damping constant $\gamma=0.01\omega_p$, the Fermi velocity $v_F=1.07\times 10^{6}$, and $\beta^2=\frac{3}{5}v_F^2$. For the GNOR model, we use the same parameters and furthermore we take $D=2.04\times10^{-4}$  [@MortensenGNOR]. An artificial absorbing boundary is set to be a concentric circle with a radius of $100$nm.

As our benchmark observable we will calculate the Extinction Cross Section (ECS, $\sigma_{\text{ext}}$), which is given by the sum of the scattering cross section $\sigma_\text{sca}$ and the absorption cross section $\sigma_\text{abs}$ [@BergEtal2009], $$\sigma_\text{ext} = \sigma_\text{sca} + \sigma_\text{abs}.$$ 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 $2 r$ of the nanowire to obtain a dimensionless normalized extinction cross section that we denote by $\sigma_\text{ext}$. It can be expressed as the sum of scaled scattering and absorption cross sections, $$\sigma_\text{sca} = \frac{1}{2r}Re \oint_S({\mathbf{E}}^\text{sca}\times \overline{H}^\text{sca})\cdot{\mathbf{n}}\, dS, \ \text{ and } \ \sigma_\text{abs} = -\frac{1}{2r}Re \oint_S({\mathbf{E}}^\text{tot}\times \overline{H}^\text{tot})\cdot{\mathbf{n}}\, dS.$$ Here the integrations are performed along a closed path around the nanowire, and $Re$ 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 \[fig:Cext\_nanowire\]. 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 \[fig:Cext\_nanowire\] 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 \[fig:Fields\_nanowire\]. 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 \[fig:Cext\_nanowire\], 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-$\mathbb{P}_2$ and HDG-$\mathbb{P}_3$ are are not presented in Figure \[fig:Cext\_nanowire\], but we found that they they lie in between the displayed curves of HDG-$\mathbb{P}_1$ and HDG-$\mathbb{P}_4$.

![The electric-field and current-density distributions of the light-matter interaction of a Na nanowire. In the upper row, we show the distributions on the fourth-order nonlocal resonance at $\omega/\omega_p=1.227$ for the NHD model. For comparison, in the bottom row we show the corresponding distributions for the GNOR model.[]{data-label="fig:Fields_nanowire"}](./FieldsNanoWire.pdf){width="95.00000%"}

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 \[tbl:perf\], where $t_\text{construction}$ denotes the CPU time for construction the matrices, $t_\text{factorization}$ denotes the CPU time used by MUMPS for the factorization of the coefficient matrix $A$ , and memory denotes the memory consumed by MUMPS. From Table \[tbl:perf\] we can see that the HDG-$\mathbb{P}_4$ is more expensive than HDG-$\mathbb{P}_1$ in CPU time for both construction and factorization. However, high-order methods are preferable because they costs less for the same accuracy [@LLP2013].

      \[tbl:perf\]                                                                                     
  -------------------- -------- ---------------------------------- ----------------------------------- -------------
                        \#DOF    $t_\text{construction}$ (second)   $t_\text{factorization}$ (second)   memory (MB)
   HDG-$\mathbb{P}_1$   28,260                0.067                               0.36                      74
   HDG-$\mathbb{P}_4$   70,650                 2.4                                 3.3                      418

  : Computational performance of the nanowire problem.

Dimer of cylindrical nanowires {#sec:dimer}
------------------------------

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 \[fig:dimerConfig\](a), and this particular configuration is from Ref. [@RazaEtal2015]. A typical mesh is shown in Figure \[fig:dimerConfig\](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-$\mathbb{P}_4$. The size of matrix for HDG-$P1$ is $105,300\times105,300$, 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 $\omega_p=1.34\times10^{16}$, damping constant $\gamma=1.14\times 10^{14}$, the Fermi velocity $v_F=1.39\times 10^{6}$, and the nonlocal parameter $\beta$ is determined by $\beta^2=\frac{3}{5}v_F^2$. 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 $p$-polarization). A comparison of the ECS curves is presented in Figure \[fig:dimECS\]. Overall, there are small but clear differences, illustrating that nonlocal response effects occur even for dimer structures for which the corresponding monomers ($r= 30$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 \[fig:dimerFields\]. 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.

![Comparison of extinction cross sections of a gold dimer as calculated with the local Drude model, the NHD model and the GNOR model. The configuration is shown in Figure \[fig:dimerConfig\] and the material parameters are given in subsection \[sec:dimer\]. []{data-label="fig:dimECS"}](./dimerECS.pdf){width="45.00000%"}

![Various field distributions in the gold dimer when illuminated by a plane wave of light. On the top line we show the distributions at the third SPR of the NHD model, at $\omega/\omega_p=0.66$. On the second line we show the corresponding distributions in the GNOR model at the same frequency.[]{data-label="fig:dimerFields"}](./FieldsDimer.pdf){width="95.00000%"}

Conclusions {#sec:con}
===========

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 {#acknowledgments .unnumbered}
===============

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 .unnumbered}
==========

[10]{} S.A. Maier, *Plasmonics: Fundamentals and Applications*, Springer, New York, 2007.

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.

S.I. Bozhevolnyi, N.A. Mortensen, Plasmonics for emerging quantum technologies, *Nanophotonics* 2016. doi: 10.1515/nanoph-2016-0179

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

D. Sarid, W.A. Challener, *Modern Introduction to Surface Plasmons: Theory, Mathematica Modeling, and Applications*, Cambridge University Press, 2010.

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.

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.

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.

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.

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.

W. Yan, M. Wubs, N.A. Mortensen, Projected dipole model for quantum plasmonics, *Phys. Rev. Lett.* 115(13): 137403. 2015.

T. Christensen, W. Yan, A.-P. Jauho, M. Solja[č]{}i[ć]{}, N.A. Mortensen, Quantum corrections in nanoplasmonics: shape, scale, and material, arXiv:1608.05421.

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.

P. Drude, Zur Elektronentheorie der Metalle, *Annalen der Physik*, 306(3):566-613, 1900.

M. Dressel, M. Scheffler, Verifying the Drude response. *Ann. Phys.*, 15(7-8): 535-544, 2006.

F. Bloch, Bremsverm[" o]{}gen von Atomen mit mehreren Elektronen, Z. Physik 81:363, 1933.

S. Raza, S. I. Bozhevolnyi, M. Wubs, N.A. Mortensen, Nonlocal optical response in metallic nanostructures, *J. Phys.: Condens. Matter*, 27: 183204, 2015.

K. Busch, M. König, J. Niegemann, Discontinuous Galerkin methods in nanophotonics, *Laser Photonics Rev.*, 5(6): 773-809, 2011.

A. Taflove, S.C. Hagness, *Computational Electrodynamics - The Finite-Difference Time-Domain Method*, Third Edition, Artech House Publishers, 2005.

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.

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.

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.

L. Li, S. Lanteri, R. Perrussel, A hybridizable discontinuous [G]{}alerkin method combined to a Schwarz algorithm for the solution of 3d time-harmonic Maxwell’s equation, *J. Comput. Phys.*, 256(1): 563-581, 2014.

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.

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.

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.

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.

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.

S. Raza, M. Wubs, S. I. Bozhevolnyi, N. Asger Mortensen, Nonlocal study of ultimate plasmon hybridization, *Opt. Lett.*, 40 (5): 839-842, 2015.

P. Jewsbury, Electrodynamic boundary conditions at metal interfaces, *J. Phys. F: Met. Phys.* 11(1): 195, 1981.

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.

B. Stupfel, Absorbing boundary conditions on arbitrary boundaries for the scalar and vector wave equations, *IEEE T. Antenn. Propag.*, 42(6): 773-780, 1994.

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.

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.

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.

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.

P. Amestoy, I. Duff, J. L’Excellent, Multifrontal parallel distributed symmetric and unsymmetric solvers, *Comput. Methods Appl. Mech. Eng.* 184: 501-520, 2000.

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.

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.

R. Ruppin, Extinction properties of thin metallic nanowires, *Opt. Commun.*, 190: 205-209, 2001.

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.

[^1]: https://www.lumerical.com/

[^2]: http://www.comsol.com/
