\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}
