## 2 Theoretical Modal and Transfer-Matrix Method for TM-polarized Surface States  
The geometry of a system formed by a periodic MMDM in contact with a homogenous isotropic dielectric medium is depicted in Figure 1. The partner dielectric medium fills the space of \(y>0\) while the MMDM fills the space of \(y<0\). The permittivities and thicknesses of the two alternating layers are \(\varepsilon_{p}\), \(L_{p}\) (blue region) and \(\varepsilon_{q}\), \(L_{q}\) (yellow region) respectively. The thickness of one unit cell is represented by \(L=L_{p}+L_{q}\). For the case of the upper dielectric layer (UDL) being adjacent to the partner medium, \(p=d\), \(q=m\), while for the case of upper metallic layer (UML) being adjacent to the partner medium, \(p=m\), \(q=d\). Here, the subscripts \(m\) and \(d\) represent the metallic and dielectric layers, respectively. We defined \(f=L_{m}/L\) as the filling factor of the metallic layer, and the permittivity of the metallic layer described by \(\varepsilon_{m}=1-1/\Omega^{2}\) was acquired by the Drude model, neglecting the effect of losses with relative frequency \(\Omega=\omega/\omega_{p}\), where \(\omega_{p}\) and \(\omega\) are the plasma frequency and the frequency of wave in vacuum, respectively. The plasma wavelength can be written as \(\lambda_{p}=2\pi c/\omega_{p}\), where \(c\) is the light velocity in vacuum. The non-dispersive dielectric layers are PbS, which has a relatively large permittivity (\(\varepsilon_{d}=18.8\)) and has been used in the similar multilayer metamaterials according to a previous report [29]. Without loss of generality, the non-dispersive partner medium was set as air with a relative permittivity \(\varepsilon_{i}=1\). The wave vector of the SWs \(\mathbf{k}_{SW}=(\beta,0,0)\) is along the \(x\) axis, where \(\beta\) is the propagation constant of the SWs. Next, we derived the dispersion equation of TM-polarized SWs supported by an MMDM and a homogenous dielectric medium based on the transfer-matrix method. The electromagnetic field of plane waves in single layer \(j\) (\(j=d,m\)) of the MMDM can be written as:  
\[\mathbf{E}_{j\pm}=\mathbf{E}_{0j\pm}\exp\big{[}ik_{0}\big{(}\pm\alpha_{j}y+ \beta_{j}x\big{)}-i\omega t\big{]}, \tag{1}\]  
\[\mathbf{H}_{j\pm}=\mathbf{k}_{j\pm}\times\mathbf{E}_{j\pm}, \tag{2}\]  
where the wave vector \(\mathbf{k}_{j\pm}=k_{0}\big{(}\beta_{j},\pm\alpha_{j},0\big{)}\), and \(\alpha_{j}\) and \(\beta_{j}\) are the projections of \(\mathbf{k}_{j\pm}\) on the \(y\) and \(x\) axes in the unit of \(k_{0}\), respectively. The \(\pm\) sign represents the positive and negative running waves along the \(y\) axis and \(\mathbf{E}_{0j\pm}\) is the polarization vector, \(\mathbf{E}_{j\pm}\) is in the \(x\)-\(y\) plane while \(\mathbf{H}_{j\pm}\) is along the \(z\) axis. Owing to the momentum conservation principle, the tangential component of the wave vector along the \(x\) axis in both dielectric and metallic layers should remain the same as the propagation constant of the SWs, that is, \(\beta_{d}=\beta_{m}=\beta\). The normal component along \(y\) axis can be written as:  
\[\alpha_{j}=\sqrt{\varepsilon_{j}-\beta^{2}}. \tag{3}\]  
The transfer matrix of single layer \(j\) for TM polarization can be given by [38]:  
\[\mathbf{\hat{M}}_{j}=\left( \begin{array}{cc}\cos\Phi_{j}&\frac{-i}{\gamma_{ j}}\sin\Phi_{j}\\ -i\gamma_{j}\sin\Phi_{j}&\cos\Phi_{j}\end{array} \right), \tag{4}\]  
where \(\Phi_{j}=k_{0}\alpha_{j}L_{j}\) and \(\gamma_{j}=\alpha_{j}z_{0}/\varepsilon_{j}\). \(z_{0}=\sqrt{\mu_{0}/\varepsilon_{0}}\) is the impedance of free space. From Equation (4), we obtained the transfer matrix for a period of the MMDM as:  
\[\mathbf{\hat{M}}_{p}=\mathbf{\hat{M}}_{d}\mathbf{\hat{M}}_{m}, \tag{5}\]  
in the UDL case, or  
\[\mathbf{\hat{M}}_{p}=\mathbf{\hat{M}}_{m}\mathbf{\hat{M}}_{d}, \tag{6}\]  
Figure 1: The geometry of a MMDM in contact with a homogenous dielectric medium. MMDM: the multilayered metal–dielectric metamaterial. in the UML case.