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 