## 2 Methods  
### Formulation  
We define a partial differential operator \(\mathbf{\mathcal{L}}\) as  
\[\mathbf{\mathcal{L}}=\mathbf{\
abla}\times(\vec{\mu}_{r}^{-1}\mathbf{\
abla}\times)-k_{0 }^{2}\vec{e}_{r}, \tag{1}\]  
where \(\vec{\mu}_{r}\) is the relative magnetic permeability tensor, \(\vec{e}_{r}\) is the relative electric permittivity tensor, and \(k_{0}\) is the wavenumber in vacuum. If we apply \(\mathbf{\mathcal{L}}\) to the electric field \(\mathbf{E}\), we obtain  
\[\mathbf{\mathcal{L}}\mathbf{E}=\mathbf{f}, \tag{2}\]  
where \(\mathbf{f}\) is the excitation of our problem, i.e.,  
\[\mathbf{f}=-jk_{0}\eta_{0}\mathbf{J}-\mathbf{\
abla}\times(\vec{\mu}_{r}^{-1}M), \tag{3}\]  
being \(\eta_{0}\) the vacuum impedance, whereas \(\mathbf{J}\) and \(\mathbf{M}\) denote the electric and magnetic currents, respectively. We use bold characters for vector magnitudes. Thus, we use as the mathematical model for the finite element analysis the well-known double-curl wave equation in a given domain \(\Omega\)[1],  
\[\mathbf{\
abla}\times(\vec{\mu}_{r}^{-1}\mathbf{\
abla}\times\mathbf{E})-k_{0}^{2}\vec{e} _{r}\mathbf{E}=-jk_{0}\eta_{0}\mathbf{J}-\mathbf{\
abla}\times(\vec{\mu}_{r}^{-1}\mathbf{M}). \tag{4}\]  
To have a unique solution, we use as boundary conditions  
\[\mathbf{\hat{n}}\times\mathbf{E} =0,\text{ on }\Gamma_{\text{D}}, \tag{5a}\] \[\mathbf{\hat{n}}\times(\vec{\mu}_{r}^{-1}\mathbf{\
abla}\times\mathbf{E}) =0,\text{ on }\Gamma_{\text{N}},\] (5b) \[\mathbf{\hat{n}}\times(\vec{\mu}_{r}^{-1}\mathbf{\
abla}\times\mathbf{E})+ \gamma\mathbf{\hat{n}}\times(\mathbf{\hat{n}}\times\mathbf{E}) =\mathbf{\Phi},\text{ on }\Gamma_{\text{C}}, \tag{5c}\]  
where \(\mathbf{\hat{n}}\) is the outward normal vector, \(\gamma\) and \(\mathbf{\Phi}\) are the propagation constant and the excitation term that depend on the boundary condition (e.g., \(\gamma=\omega\sqrt{\mu\epsilon},\mathbf{\Phi}=0\) for absorbing boundary conditions), and \(\Gamma_{\text{D}}\), \(\Gamma_{\text{N}}\), and \(\Gamma_{\text{C}}\) stand for Dirichlet, Neumann, and Cauchy boundary conditions, respectively. The vector electric field \(\mathbf{E}\) belongs to the curl-conforming space [24],\[\mathbf{H}(\text{curl})=\{\mathbf{N}\in(L^{2}(\Omega))^{3}:\mathbf{\
abla}\times\mathbf{N}\in(L^{2 }(\Omega))^{3}\}, \tag{6}\]  
where \(L^{2}(\Omega)\) is the space of square-integrable functions and, specifically, here we use the space  
\[\mathbf{H}_{0}(\text{curl},\Omega)=\{\mathbf{N}\in\mathbf{H}(\text{curl},\Omega):\mathbf{\hat{n }}\times\mathbf{N}=0\in\Gamma_{\text{D}}\} \tag{7}\]  
with homogeneous Dirichlet boundary conditions on \(\Gamma_{D}\). Now, we define the residual \(\mathbf{R}=\mathbf{\mathcal{L}}\mathbf{E}-\mathbf{f}\) that we test with the weight function \(\mathbf{F}\) to obtain the variational formulation, i.e.,  
\[<\mathbf{F},\mathbf{\mathcal{L}}\mathbf{E}-\mathbf{f}>=0, \tag{8}\]  
so we need to find \(\mathbf{E}\in\mathbf{H}_{0}(\text{curl},\Omega)\) such that  
\[ \begin{split}&\int_{\Omega}(\mathbf{\
abla}\times\mathbf{F})\cdot(\bar{ \mu}^{-1}\mathbf{\
abla}\times\mathbf{E})d\Omega+\gamma\int_{\Gamma_{c}}(\mathbf{\hat{n}} \times\mathbf{F})\cdot(\mathbf{\hat{n}}\times\mathbf{E})d\Gamma_{c}\\ &-k_{0}^{2}\int_{\Omega}\mathbf{F}\cdot(\bar{\mathbf{\epsilon}}\mathbf{E})d \Omega-jk_{0}\eta_{0}\int_{\Omega}\mathbf{F}\cdot\mathbf{f}d\Omega\\ &-\int_{\Omega}\mathbf{F}\cdot\mathbf{\
abla}\times(\bar{\mu}_{r}^{-1} \mathbf{M})d\Omega+\int_{\Gamma_{c}}\mathbf{F}\cdot\mathbf{\Phi}d\Gamma_{c}=0,\end{split}  \tag{9}\]  
for all \(\mathbf{F}\in\mathbf{H}_{0}(\text{curl},\Omega)\). This formulation is discretized with a collection of finite elements (e.g., tetrahedral mesh) that we assemble for each element \(e\) obtaining  
\[\sum_{\forall e}<\mathbf{F}^{e},\mathbf{\mathcal{L}}\mathbf{E}^{e}-\mathbf{F}^{e}>=0, \tag{10}\]  
where we have used the superscript \(e\) to denote magnitudes relative to each element.