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}
