## Appendix A Stability  
In this appendix, we study the stability of the algorithm by means of the combination of the von Neumann method with the Routh-Hurwitz criterion [40]. The von Neumann method employs a Fourier series expansion of the error at the mesh nodes at a given time instant \(t=s\Delta_{t}\). Due to linearity, only a single term of this expansion needs to be considered, i.e.,  
\[\vec{E}^{s}|_{i,j,k} =\widetilde{\vec{E}}_{0}y^{\delta}e^{j\left(i\Delta_{x}\widetilde {k}_{x}+j\Delta_{y}\widetilde{k}_{y}+k\Delta_{x}\widetilde{k}_{z}\right)} \tag{11}\] \[\vec{D}^{s}|_{i,j,k} =\widetilde{\vec{D}}_{0}y^{\delta}e^{j\left(i\Delta_{x}\widetilde {k}_{x}+j\Delta_{y}\widetilde{k}_{y}+k\Delta_{x}\widetilde{k}_{z}\right)} \tag{12}\]  
where \(\widetilde{\vec{E}}_{0}\) and \(\widetilde{\vec{D}}_{0}\) are complex amplitudes, indexes \(i,j,k\) denote the position of the nodes in the mesh, \(\Delta_{x}\), \(\Delta_{y}\), \(\Delta_{z}\) are the sizes of the discretization cell, and \(\widetilde{k}_{x}\), \(\widetilde{k}_{y}\) and \(\widetilde{k}_{z}\) are the wavenumbers of the discrete modes. The factor \(y\) in (11) and (12) is a complex variable, often called the amplification factor, which gives the growth of the error in a time updating. To ensure that a finite-difference algorithm is stable, the error must not grow as the time increases, thus, the condition \(|y|\leq 1\) must be satisfied. The stability condition for a particular FDTD scheme leads to a characteristic polynomial in \(y\):  
\[\Phi(y)=\sum_{\ell=0}^{N_{p}}a_{\ell}y^{N_{p}-\ell} \tag{13}\]  
Assuming that \(y_{\ell}\) are the roots of \(\Phi(y)\), the condition for stability implicates that the roots \(y_{\ell}\) must lie inside or on the unit circle in the Y-plane. In this stability purpose, the Routh-Hurwitz criterion provides a method that ensures the roots \(y_{\ell}\) on the unit circle in the Y-plane. It establishes that the polynomial \(\Phi(r)=\sum_{\ell=0}^{N_{p}}b_{\ell}y^{N_{p}-\ell}\) with constant real coefficients \(b_{\ell}\) has no roots in the right-half of the r-plane if all the entries of the first column of the Routh table are non-negative quantities. In [40], there is a error that makes it impossible to build the Routh table, so we explain here how to build up the Routh table. The first and second rows of the coefficients correspond to even and odd powers of r, respectively. For simplicity in the notation, and without loss of generality, it is usual assuming \(N_{p}\) to be an even number, therefore \(c_{1,\ell}=b_{2\ell}\), with \(\ell\)= 0, 1, 2,...,\(\frac{N_{p}}{2}\) and \(c_{2,\ell}=b_{2\ell+1}\), with \(\ell\)= 0, 1, 2,...,\(\frac{N_{p}}{2}-1\). The remaining entries of the table are obtained by using the following expression \(c_{m,\ell}=\frac{c_{m,-2\ell+1}c_{m,-1,0}-c_{m,-2\ell}c_{m,-1,\ell+1}}{c_{m,-1,0}}\). The Routh-Hurwitz criterion can be used to determine if any root of the stability polynomial \(\Phi(y)\) lies outside the unit circle in the y-plane. To this end, the bilinear transformation \(y=\frac{r+1}{r-1}\) is carried out in \(\Phi(y)\). This transformation maps the exterior of the unit circle in the y-plane onto the right-half of the r-plane. To satisfy the von Neumann stability condition for a different scheme, as a function of the parameters of interest, all the entries of the first column of the Routh table are forced to be non-negative quantities. Nevertheless, before beginning the search for stability conditions, some definitions are necessary.