## Appendix B Derivation of Surface Wave Modes Based on [PERSON]'s Method  
In [22], [PERSON] introduced an approximate method to describe how the guided waves propagate through rectangular dielectric optical waveguides. Based on his method, the surface waves are divided into \(E^{z}\) and \(E^{y}\) modes. For \(E^{z}\) mode, the field component \(E_{z}\) is strong, while \(H_{z}\) is assumed to be zero, which represents the quasi-TM mode in our discussion. Similarly, the \(E^{y}\) mode represents quasi-TE mode in our discussion, in which case the field \(E_{y}\) is strong and \(H_{y}\) is assumed to be zero. In addition, the fields in Regions V and VI as shown in Figure 2b are assumed to be negligible. Based on these assumptions, the quasi-TM and quasi-TE modes can be analyzed separately.  
For quasi-TM mode, \(H_{z}\) and \(H_{y}\) of Regions I, II, III, and IV can be expressed as:  
\[H_{z}^{\text{I}}=H_{z}^{\text{II}}=H_{z}^{\text{III}}=H_{z}^{\text{IV}}=0, \tag{7}\]  
\[H_{y}^{\text{I}}=A_{\text{I}}\Big{(}e^{-R_{y}y}+r_{1}e^{R_{y}y}\Big{)}e^{-j \beta_{x}x}\cos(k_{z}z), \tag{8}\]  
\[H_{y}^{\text{II}}=A_{\text{II}}\Big{(}e^{-R_{y}y}+r_{1}e^{R_{y}y}\Big{)}e^{- j\beta_{x}x}\text{e}^{-k_{0z}z}, \tag{9}\]  
\[H_{y}^{\text{III}}=A_{\text{III}}e^{-R_{0y}y}e^{-j\beta_{x}x}\cos(k_{z}z), \tag{10}\]  
\[H_{y}^{\text{IV}}=A_{\text{IV}}e^{R_{0y}y}e^{-j\beta_{x}x}\cos(k_{z}z). \tag{11}\]Based on the fields \(H_{z}\) and \(H_{y}\), all field distributions and dispersion relations in four regions can then be derived by solving [PERSON]'s equations.  
The continuity conditions of \(E_{x}\) and \(H_{y}\) at \(z=t\) give  
\[k_{z}\tan(k_{z}t)=\varepsilon_{x}k_{0z}. \tag{101}\]  
In addition, \(E_{z}\) and \(H_{x}\) are continuous at \(y=\pm a\), which leads to:  
\[R_{0y}\mu_{x}\Big{(}\varepsilon_{x}\mu_{y}k_{0}^{2}-k_{z}^{2}\Big{)}\frac{e^{- R_{y}a}+e^{R_{y}a}}{e^{-R_{y}a}-e^{R_{y}a}}=\varepsilon_{x}\mu_{y}R_{y} \Big{(}k_{0}^{2}-k_{z}^{2}\Big{)}. \tag{102}\]  
The dispersion relations in four regions are derived as:  
\[\beta_{x}^{2}-\frac{\mu_{y}}{\mu_{x}}R_{y}^{2}+\frac{\varepsilon_{z}}{ \varepsilon_{x}}k_{z}^{2}=\varepsilon_{z}\mu_{y}k_{0}^{2}, \tag{103}\]  
\[\beta_{x}^{2}-R_{y}^{2}-k_{0z}^{2}=k_{0}^{2}, \tag{104}\]  
\[\beta_{x}^{2}-R_{0y}^{2}+k_{z}^{2}=k_{0}^{2}. \tag{105}\]  
It is seen that \(\beta_{x}\), \(R_{y}\), \(R_{0y}\), \(k_{z}\), and \(k_{\alpha z}\) can be solved for given \(k_{0}\), \(t\), \(\stackrel{{=}}{{\varepsilon}}\), \(\stackrel{{=}}{{\mu}}\) based on Equations (101)-(105).  
For quasi-TE mode, \(H_{y}\) and \(H_{z}\) in Regions I, II, III, and IV can be expressed as:  
\[H_{y}^{\rm I}=H_{y}^{\rm III}=H_{y}^{\rm III}=H_{y}^{\rm IV}=0, \tag{106}\]  
\[H_{z}^{\rm I}=B_{\rm I}\Big{(}e^{-R_{y}y}+r_{2}e^{R_{y}y}\Big{)}e^{-j\beta_{x }x}\sin(k_{z}z), \tag{107}\]  
\[H_{z}^{\rm II}=B_{\rm II}\Big{(}e^{-R_{y}y}+r_{2}e^{R_{y}y}\Big{)}e^{-j\beta_ {x}x}e^{-k_{0z}z}, \tag{108}\]  
\[H_{z}^{\rm III}=B_{\rm III}e^{-R_{0y}y}e^{-j\beta_{x}x}\sin(k_{z}z), \tag{109}\]  
\[H_{z}^{\rm IV}=B_{\rm IV}e^{R_{0y}y}e^{-j\beta_{x}x}\sin(k_{z}z). \tag{110}\]  
Similarly, based on [PERSON]'s equations and boundary conditions, five equations with five unknowns \(\beta_{x}\), \(R_{y}\), \(R_{0y}\), \(k_{z}\), and \(k_{0z}\) can be derived, and they are summarized below:  
\[-k_{z}\cot(k_{z}t)=\mu_{x}k_{0z}, \tag{111}\]  
\[R_{y}\Big{(}k_{0}^{2}+R_{0}^{2}\Big{)}\frac{e^{-R_{y}a}-e^{R_{y}a}}{e^{-R_{y}a }+e^{R_{y}a}}=R_{0y}\varepsilon_{y}\Big{(}\varepsilon_{x}\mu_{x}k_{0}^{2}+R_{ y}^{2}\Big{)}, \tag{112}\]  
\[\beta_{x}^{2}-\frac{\varepsilon_{y}}{\varepsilon_{x}}R_{y}^{2}+\frac{\mu_{z}}{ \mu_{x}}k_{z}^{2}=\varepsilon_{y}\mu_{z}k_{0}^{2}, \tag{113}\]  
\[\beta_{x}^{2}-R_{y}^{2}-k_{0z}^{2}=k_{0}^{2}, \tag{114}\]  
\[\beta_{x}^{2}-R_{0y}^{2}+k_{z}^{2}=k_{0}^{2}. \tag{115}\]