## Appendix C Generalization of Stability Condition Towards Three-Dimensional FO FDTD
\begin{bmatrix}\mathbf{A}&\mathbf{B}\\ \mathbf{C}&\mathbf{D}\end{bmatrix}  \tag{102}\]  
where \(\mathbf{A}\), \(\mathbf{B}\), \(\mathbf{C}\), \(\mathbf{D}\) are \(3\times 3\) matrices. They commute because \(\mathbf{B}=-\mathbf{C}\) and \(\frac{1}{\mu_{\alpha}}\mathbf{A}=\frac{1}{\varepsilon_{\alpha}}\mathbf{D}=- \mathbf{I}\), where \(\mathbf{I}\) is the identity matrix. In such a case, the determinant of matrix \(\mathbf{U}\) can be calculated using the following property [28]:  
\[\det(\mathbf{U})=\det(\mathbf{A}\mathbf{D}-\mathbf{B}\mathbf{C}). \tag{103}\]  
Hence, one obtains  
\[\det(\mathbf{U}) = \det\left( \begin{bmatrix}\mu_{\alpha}\varepsilon_{\alpha}-{ \delta_{y}}^{2}-{\delta_{z}}^{2}&\delta_{x}\delta_{y}&\delta_{x}\delta_{z}\\ \delta_{x}\delta_{y}&\mu_{\alpha}\varepsilon_{\alpha}-{\delta_{x}}^{2}-{ \delta_{z}}^{2}&\delta_{y}\delta_{z}\\ \delta_{x}\delta_{z}&\delta_{y}\delta_{z}&\mu_{\alpha}\varepsilon_{\alpha}-{ \delta_{x}}^{2}-{\delta_{y}}^{2}\end{bmatrix} \right) \tag{104}\] \[= \mu_{\alpha}\varepsilon_{\alpha}({\delta_{x}}^{2}+{\delta_{y}}^{2 }+{\delta_{z}}^{2}-\mu_{\alpha}\varepsilon_{\alpha})^{2}.\]  
Equating the determinant to 0, one obtains  
\[\mu_{\alpha}\varepsilon_{\alpha}={\delta_{x}}^{2}+{\delta_{y}}^{2}+{\delta_{z }}^{2}. \tag{105}\]  
Using (102) one obtains  
\[\left(\frac{e^{i\omega\Delta t\frac{1}{2}}\big{(}1-e^{-j\omega \Delta t}\big{)}^{\alpha}}{2j}\right)^{2} = \left((\Delta t)^{\alpha}\frac{e^{-j[\frac{\gamma}{y}\frac{1}{2} \Delta x]}-e^{j[\frac{\gamma}{y}\frac{1}{2}\Delta x]}}{2j\sqrt{\mu_{\alpha} \varepsilon_{\alpha}}\Delta x}\right)^{2}+ \tag{106}\] \[\left((\Delta t)^{\alpha}\frac{e^{-j[\frac{\gamma}{y}\frac{1}{2} \Delta y]}-e^{j[\frac{\gamma}{y}\frac{1}{2}\Delta y]}}{2j\sqrt{\mu_{\alpha} \varepsilon_{\alpha}}\Delta y}\right)^{2}+\left((\Delta t)^{\alpha}\frac{e^{- j[\frac{\gamma}{y}\frac{1}{2}\Delta z]}-e^{j[\frac{\gamma}{y}\frac{1}{2}\Delta z]}}{2j \sqrt{\mu_{\alpha}\varepsilon_{\alpha}}\Delta z}\right)^{2}.\]  
From (106), one directly obtains the numerical-dispersion relation for the three-dimensional FO FDTD algorithm  
\[\left(\frac{e^{i\omega\Delta t\frac{1}{2}}\big{(}1-e^{-j\omega \Delta t}\big{)}^{\alpha}}{2j}\right)^{2} = \frac{(\Delta t)^{2\alpha}}{\mu_{\alpha}\varepsilon_{\alpha}( \Delta x)^{2}}\bigg{[}\sin\!\left(\frac{\gamma}{x}\Delta x}{2}\right)\bigg{]}^ {2}+ \tag{107}\] \[\frac{(\Delta t)^{2\alpha}}{\mu_{\alpha}\varepsilon_{\alpha}( \Delta y)^{2}}\bigg{[}\sin\!\left(\frac{\tilde{\gamma}_{y}\Delta y}{2}\right) \bigg{]}^{2}+\frac{(\Delta t)^{2\alpha}}{\mu_{\alpha}\varepsilon_{\alpha}( \Delta z)^{2}}\bigg{[}\sin\!\left(\frac{\tilde{\gamma}_{z}\Delta z}{2}\right) \bigg{]}^{2}.\]  
If one omits the \(x\)- and \(y\)-related terms in (107), it simplifies that formula to the one-dimensional case (103). Substituting  
\[\xi_{\alpha} = \left(\frac{(\Delta t)^{2\alpha}}{\mu_{\alpha}\varepsilon_{ \alpha}(\Delta x)^{2}}\bigg{[}\sin\!\left(\frac{\tilde{\gamma}_{x}\Delta x}{2} \right)\bigg{]}^{2}+\frac{(\Delta t)^{2\alpha}}{\mu_{\alpha}\varepsilon_{ \alpha}(\Delta y)^{2}}\bigg{[}\sin\!\left(\frac{\tilde{\gamma}_{y}\Delta y}{2} \right)\bigg{]}^{2}+\right. \tag{108}\] \[\left.\frac{(\Delta t)^{2\alpha}}{\mu_{\alpha}\varepsilon_{ \alpha}(\Delta z)^{2}}\bigg{[}\sin\!\left(\frac{\tilde{\gamma}_{x}\Delta z}{2} \right)\bigg{]}^{2}\right)^{\frac{1}{2\alpha}}\]into (A54), and continuing as in the one-dimensional case, one obtains  
\[ \begin{split} 2^{\frac{1}{\lambda}}\left(\frac{(\Delta t)^{2\alpha}}{ \mu_{\alpha}\varepsilon_{\alpha}(\Delta x)^{2}}\bigg{[}\sin\left(\frac{\tilde{ \imath}_{x}\Delta x}{2}\right)\bigg{]}^{2}+\frac{(\Delta t)^{2\alpha}}{\mu_{ \alpha}\varepsilon_{\alpha}(\Delta y)^{2}}\bigg{[}\sin\left(\frac{\tilde{ \jmath}_{y}\Delta y}{2}\right)\bigg{]}^{2}+\right.\\ \left.\frac{(\Delta t)^{2\alpha}}{\mu_{\alpha}\varepsilon_{\alpha} (\Delta z)^{2}}\bigg{[}\sin\left(\frac{\tilde{\kappa}_{z}\Delta z}{2}\right) \bigg{]}^{2}\right)^{\frac{1}{2\alpha}}\leq 2.\end{split}