## 3 Dispersion Relations, Group Velocities and Frequency Range of Backward Surface Waves  
In accordance with the transfer-matrix method, we depicted the dispersion relations of the TM-polarized surface states with a metallic layer filling factor \(f=0.7\) for the case of UDL, as shown in Figure 2a. In contrast, the dispersion curve of the SWs obtained through the EMA method, and one of the conventional SPPs at a metal-dielectric interface are also presented. Even in the condition, in which the unit-cell thickness is far smaller than the plasma wavelength (\(L=0.01\lambda_{p}\)), the dispersion curves acquired by the transfer-matrix method and the EMA method match each other only when the frequency is evidently below the pole frequency \(\Omega_{\text{E}}\), at which the propagation constant tends to infinity according to the EMA method. Otherwise, when the frequency is in the vicinity of \(\Omega_{\text{E}}\), the difference between dispersion curves obtained by the transfer-matrix method and the EMA method becomes significant. The propagation constant tends to infinity at \(\Omega_{\text{S}}\), which is the bulk plasmon frequency of the conventional SPPs and satisfies \(\varepsilon_{m}(\Omega_{\text{S}})+\varepsilon_{d}=0\), instead of \(\Omega_{\text{E}}\). There is a redshift of the dispersion curve calculated by the transfer-matrix method relative to that calculated by EMA method. More importantly, part of the former is anomalously dispersive, thus leading to negative group velocities, as shown in Figure 2b. Hence, the SWs in this condition are backward waves, of which the actual energy propagation direction is opposite to the wave vector. According to Figure 2a, for smaller propagation constant \(\beta\) close to the light line, the field variation within one unit cell can be regarded as quasi-static. Therefore, relevant dispersions keep consistent with the prediction of the EMA method. However, for larger \(\beta\), when the wavelength of the SWs is roughly ten times larger than the unit-cell thickness \(L\), the optical nonlocality effect, that is, the field variation within one unit cell becomes much more intensive, appears due to the excitation of surface plasmon modes; thus, the responding dispersions are beyond the EMA method. Under this circumstance, when \(L\) is rather small, not exceeding 10% of the plasma wavelength \(\lambda_{p}\), the coupling of the surface plasmon modes is relatively strong, leading to the anomalous dispersions. The upper threshold frequency \(\Omega_{\text{up}}\) of the anomalous dispersions can be located for the corresponding group velocity \(v_{g}=0\), as shown in Figure 2b, while the lower threshold frequency \(\Omega_{\text{low}}\) is equal to \(\Omega_{\text{S}}\). With \(L\) increased above 10% of \(\lambda_{p}\), the coupling of the surface plasmon modes weaken; thus, the corresponding dispersion curve turns back into the normal dispersion before \(\beta\) tends to infinity at \(\Omega_{\text{S}}\). In this case, \(\Omega_{\text{low}}\) is inferior to \(\Omega_{\text{S}}\). In both cases, the width of the anomalous dispersion frequency range always shrinks with the increasing \(L\), as shown in Figure 2c. When \(L\) is above 40% of \(\lambda_{p}\), the whole curve becomes normal dispersive, indicating that backward SWs are prohibited. In the limiting case with an extremely large \(L\), the plasmon modes become uncoupled, and hence, the dispersion curve eventually coincides with the one of the conventional SPPs excited at a single interface between a semi-infinite metal and a PbS dielectric medium. Moreover, the sensitivity of \(\Omega_{\text{up}}\) modulated by \(L\) is within \(0.28\omega_{p}/\lambda_{p}\sim 0.34\omega_{p}/\lambda_{p}\) (0.82 THz/nm\(\sim\)0.99 THz/nm for \(\lambda_{p}=320\) nm), while the one of \(\Omega_{\text{low}}\) is within \(0\sim 0.1\omega_{p}/\lambda_{p}\) (0-0.29 THz/nm for \(\lambda_{p}=320\) nm).