## Appendix A Derivation of Numerical-Dispersion Relation for FO FDTD
\begin{split}\Big{|}1-e^{\hat{\omega}_{i}\Delta t}e^{-j\hat{ \omega}_{r}\Delta t}\Big{|}&=\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}}\\ &=\sqrt{1-2e^{\hat{\omega}_{i}\Delta t}\cos(\hat{\omega}_{r} \Delta t)+e^{2\hat{\omega}_{i}\Delta t}\cos^{2}(\hat{\omega}_{r}\Delta t)+e^{ 2\hat{\omega}_{i}\Delta t}\sin^{2}(\hat{\omega}_{r}\Delta t)}\\ &=\sqrt{1-2e^{\hat{\omega}_{i}\Delta t}\cos(\hat{\omega}_{r} \Delta t)+e^{2\hat{\omega}_{i}\Delta t}}\\ \end{split}  \tag{146}\]  
\[\text{Arg}\left(1-e^{\hat{\omega}_{i}\Delta t}e^{-j\hat{\omega}_{r}\Delta t }\right)=\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{]}. \tag{147}\]  
Finally, one obtains  
\[\Big{|}e^{\frac{-\hat{\omega}_{i}\Delta t}{2\hat{\omega}}}e^{j\frac{\hat{ \omega}_{i}\Delta t}{2\hat{\omega}}}\Big{(}1-e^{\hat{\omega}_{i}\Delta t}e^{- j\hat{\omega}_{r}\Delta t}\Big{)}\Big{|}=e^{-\frac{\hat{\omega}_{i}\Delta t}{2 \hat{\omega}}}\sqrt{1-2e^{\hat{\omega}_{i}\Delta t}\cos(\hat{\omega}_{r}\Delta t )+e^{2\hat{\omega}_{i}\Delta t}} \tag{148}\]  
\[\text{Arg}\left(e^{\frac{-\hat{\omega}_{i}\Delta t}{2\hat{\omega}}}e^{j\frac{ \hat{\omega}_{i}\Delta t}{2\hat{\omega}}}\Big{[}1-e^{\hat{\omega}_{i}\Delta t} e^{-j\hat{\omega}_{r}\Delta t}\Big{]}\right)=\frac{\hat{\omega}_{r}\Delta t}{2 \alpha}+\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{]}. \tag{149}\]  
Comparing modules and arguments of the left- and right-hand sides of (144), one obtains  
\[e^{-\frac{\hat{\omega}_{i}\Delta t}{2\hat{\omega}}}\sqrt{1-2e^{\hat{\omega}_{ i}\Delta t}\cos(\hat{\omega}_{r}\Delta t)+e^{2\hat{\omega}_{i}\Delta t}}=2^{ \frac{1}{2}}\xi_{\alpha} \tag{150}\]  
\[\frac{\hat{\omega}_{r}\Delta t}{2\alpha}+\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{]}=\frac{\pi}{2\alpha}. \tag{151}\]  
Analyzing the argument relation (151), one obtains  
\[\frac{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)}=\text{tg}\left(\frac{\pi- \hat{\omega}_{r}\Delta t}{2\alpha}\right). \tag{152}\]  
The sampling time \(\Delta t\) must respect the Nyquist theorem; therefore, the only proper solution to (152) is \(\hat{\omega}_{r}=\frac{\pi}{\Delta t}\). Let us substitute \(\cos(\hat{\omega}_{r}\Delta t)=\cos(\frac{\pi}{\Delta t}\Delta t)=\cos(\pi)=-1\) into the left-hand side of (150), which gives  
\[ \begin{split} e^{-\frac{\hat{\omega}_{i}\Delta t}{2\alpha}}\sqrt{1 +2e^{\hat{\omega}_{i}\Delta t}+e^{2\hat{\omega}_{i}\Delta t}}&=e^{- \frac{\hat{\omega}_{i}\Delta t}{2\alpha}}\sqrt{(e^{\hat{\omega}_{i}\Delta t}+1) ^{2}}\\ &=e^{-\frac{\hat{\omega}_{i}\Delta t}{2\alpha}}(e^{\hat{\omega}_{i} \Delta t}+1).\end{split}  \tag{153}\]  
Let us introduce auxiliary variables \(e^{\hat{\omega}_{i}\Delta t}=u\) and \(\frac{1}{2\alpha}=\beta\). Using (150) and (153), one then obtains  
\[u^{-\beta}(1+u)=2^{\frac{1}{2}}\xi_{\alpha} \tag{154}\]  
where \(\frac{1}{2}<\alpha<1\), \(\frac{1}{2}<\beta<1\). Equation (154) can be written as  
\[u+1=2^{\frac{1}{2}}\xi_{\alpha}u^{\beta}. \tag{155}\]A simulation is stable when \(\hat{\omega}_{i}\geq 0\), which is equivalent to \(u=e^{\hat{\omega}_{i}\Delta t}\geq 1\). Equation (A25) is solved graphically and visualised in Figure A1. For \(\xi_{\alpha}\) larger than the boundary value equal 2, (A25) has two solutions, i.e., one smaller and one larger than 1. In this case a simulation is considered unstable. A smaller value of \(\xi_{\alpha}\) means that either two solutions larger than 1, one such solution, or even no such solutions exist. All these cases are considered to be stable, as they lack any unstable term smaller than 1.