## Appendix A Derivation of Numerical-Dispersion Relation for FO FDTD  
For clarity of presentation, the one-dimensional dispersion relation and stability condition are derived first, with subsequent generalisations to the three-dimensional case in Appendix C.  
Let us consider one-dimensional discretized FO [PERSON]'s equations  
\[\frac{E_{x}^{n}(k+1)-E_{x}^{n}(k)}{\Delta z}=-\mu_{a}\left(\frac{1}{\Delta t} \right)^{\alpha}\sum_{l=0}^{n}w_{l}^{\alpha}H_{y}^{n-l+\frac{1}{2}}(k+\frac{1 }{2}) \tag{11}\]  
\[-\frac{H_{y}^{n+\frac{1}{2}}(k+\frac{1}{2})-H_{y}^{n+\frac{1}{2}}(k-\frac{1}{ 2})}{\Delta z}=\varepsilon_{x}\left(\frac{1}{\Delta t}\right)^{\alpha}\sum_{l =0}^{n}w_{l}^{x}E_{x}^{n-l+1}(k). \tag{12}\]  
Let us assume the existence of harmonic solutions with the field distributions  
\[H_{y}^{n}(k)=H_{y0}e^{j(\tilde{\omega}n\Delta t-\tilde{k}_{z}k\Delta z)} \tag{13}\]  
\[E_{x}^{n}(k)=E_{x0}e^{j(\tilde{\omega}n\Delta t-\tilde{k}_{z}k\Delta z)} \tag{14}\]  
where \(\tilde{\omega}\in\mathbb{C}\), \(\tilde{k}_{z}\in\mathbb{R}\). Substituting (14) into (11) gives  
\[H_{y0}\frac{\sum_{l=0}^{N}[w_{l}^{\alpha}e^{j(\tilde{\omega}n-l+\frac{1}{2}) \Delta t}e^{-j\tilde{k}_{z}(k+\frac{1}{2})\Delta z}]}{(\Delta t)^{\alpha}}=- \frac{1}{\mu_{\alpha}}E_{x0}\frac{e^{j(\tilde{\omega}n\Delta t}e^{-j\tilde{k}_ {z}(k+1)\Delta z}-e^{j\tilde{k}_{z}k\Delta z}}{\Delta z} \tag{15}\]  
\[H_{y0} =-\frac{(\Delta t)^{\alpha}E_{x0}}{\mu_{\alpha}\Delta z}\frac{e^{j( \tilde{\omega}n\Delta t}e^{-j\tilde{k}_{z}(k+\frac{1}{2})\Delta z}\left(e^{-j \tilde{k}_{z}(\frac{1}{2})\Delta z}-e^{-j\tilde{k}_{z}(-\frac{1}{2})\Delta z }\right)}{e^{-j\tilde{k}_{z}(k+\frac{1}{2})\Delta z}\sum_{l=0}^{N}w_{l}^{ \alpha}e^{j\tilde{\omega}(n-l+\frac{1}{2})\Delta t}} \tag{16}\] \[=\frac{(\Delta t)^{\alpha}E_{x0}}{\mu_{\alpha}\Delta z}\frac{e^{ j(\tilde{\omega}n\Delta t}\left(e^{j\tilde{k}_{z}\frac{1}{2}\Delta z}-e^{-j \tilde{k}_{z}\frac{1}{2}\Delta z}\right)}{e^{j\tilde{\omega}n\Delta t}\sum_{l =0}^{N}w_{l}^{\alpha}e^{j\tilde{\omega}(-l+\frac{1}{2})\Delta t}}\] \[=\frac{(\Delta t)^{\alpha}E_{x0}}{\mu_{\alpha}\Delta z}\frac{2j \frac{j(\tilde{k}_{z}\frac{1}{2}\Delta z}{2}-e^{-j\tilde{k}_{z}\frac{1}{2} \Delta z}}{2j}\] \[=E_{x0}\frac{(\Delta t)^{\alpha}}{\mu_{\alpha}\Delta z}\frac{2j \sin\left(\frac{\tilde{k}_{z}\Delta z}{2}\right)}{\sum_{l=0}^{N}w_{l}^{\alpha }e^{j\tilde{\omega}(-l+\frac{1}{2})\Delta t}}.\]  
After substituting (14) into (12), one obtains  
\[E_{x0}\frac{\sum_{l=0}^{N}[w_{l}^{\alpha}e^{j\tilde{\omega}(n-l+1)\Delta t}e^ {-j\tilde{k}_{z}k\Delta z}]}{(\Delta t)^{\alpha}}=-\frac{1}{\varepsilon_{x}}H _{y0}\frac{e^{j\tilde{\omega}(n+\frac{1}{2})\Delta t}e^{-j\tilde{k}_{z}(k+ \frac{1}{2})\Delta z}-e^{j\tilde{\omega}(n+\frac{1}{2})\Delta t}e^{-j\tilde{k} _{z}(k-\frac{1}{2})\Delta z}}{\Delta z} \tag{17}\]  
After substituting \(\hat{\omega}=\hat{\omega}_{r}+j\hat{\omega}_{i}\), one obtains  
\[-\frac{H_{y}^{n+\frac{1}{2}}(k+\frac{1}{2})-H_{y}^{n+\frac{1}{2}}(k-\frac{1}{ 2})}{\Delta z}=\varepsilon_{x}\left(\frac{1}{\Delta t}\right)^{\alpha}\sum_{l =0}^{n}w_{l}^{x}E_{x}^{n-l+1}(k). \tag{12}\]0  
Then, \(1-e^{\hat{\omega}_{i}\Delta t}e^{-j\hat{\omega}_{r}\Delta t}\) can be expressed in the exponential form as follows:  
\[ \begin{split} 1-e^{\hat{\omega}_{i}\Delta t}e^{-j\hat{\omega}_{r} \Delta t}&=1-e^{\hat{\omega}_{i}\Delta t}[\cos(-\hat{\omega}_{r} \Delta t)+j\sin(-\hat{\omega}_{r}\Delta t)]\\ &=1-e^{\hat{\omega}_{i}\Delta t}\cos(\hat{\omega}_{r}\Delta t)+je ^{\hat{\omega}_{i}\Delta t}\sin(\hat{\omega}_{r}\Delta t)\\ &=\sqrt{[1-e^{\hat{\omega}_{i}\Delta t}\cos(\hat{\omega}_{r} \Delta t)]^{2}+[e^{\hat{\omega}_{i}\Delta t}\sin(\hat{\omega}_{r}\Delta t)]^{2 }}e^{j\text{arctg}2\big{[}e^{\hat{\omega}_{i}\Delta t}\sin(\hat{\omega}_{r} \Delta t),1-e^{\hat{\omega}_{i}\Delta t}\cos(\hat{\omega}_{r}\Delta t)\big{]} },\end{split}  \tag{145}\]  
Hence, one obtains  
\[