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 