## 3 Dispersion Relations, Group Velocities and Frequency Range of Backward Surface Waves
Hence, the optical nonlocality in this case is mainly caused by the surface plasmon modes excited at the boundary between the dielectric and metallic layers. However, for the UML case, optical nonlocality is mostly due to the surface plasmon modes at the boundary between the upper metallic layer and the partner medium according to the location of the field confinement shown in Figure 5e. In addition, from Figure 5b,d,f, it should be noted that \(S_{t}\) remains positive in the dielectric layers and the partner dielectric medium but negative in the metallic  
Figure 4: Dispersion relations of the TM-polarized surface states supported by a homogenous dielectric medium and an MMDM in the UML case. The values of material parameters are the same as those in Figure 2. **(a)** Dispersion relations with different unit-cell thicknesses when the filling factor is set as \(f=0.7\). **(b)** Dispersion relations with different filling factors when the unit-cell thickness is set as \(L=0.3\lambda_{p}\). layers, indicating that the energy flow in the dielectric layers and the partner dielectric medium is in the same direction with the wave vector of the SW, while the energy flow in the metallic layers is opposite to the direction of the wave vector of the SW. For the anomalous dispersive SW shown in Figure 5b, the energy flow in metallic layers is dominant; therefore, the SW unfolds backward wave characteristics. On the contrary, for the normal dispersive SWs shown in Figure 5d,f, the energy flow in the partner dielectric medium predominates; thus, these SWs unfold forward wave characteristics. Numerical Simulations of Backward Surface Waves Excited by a TM-Polarized Gaussian Beam Based on the Prism-Coupled Configuration  
The prism-coupled configuration has been widely used in the excitation of various types of surface modes, such as SPPs [9], Dyakonov SWs [11] and Dyakonov-Tamm waves [20]. Our aim in this section is to verify the backward SWs of anomalous dispersion relations based on a prism-coupled configuration through the finite-element method. As depicted in Figure 6, a prism with a large refraction index \(n_{p}\) is placed onto the partner medium with a confined thickness \(h\). Both the permittivities of the partner medium and the dielectric layers remain unchanged. Additionally, we extended the simulation to a more practical case when the permittivity of the metallic layers becomes complex in the Drude form as \(\varepsilon_{m}=1-\omega_{p}^{2}/[\omega(\omega+i\Gamma)]\)[39], considering the metallic losses. Figure 5: The profiles of the tangential components of the magnetic field and the time-averaged Poynting vector of TM-polarized surface states with **(a,b)**\(f=0.7\), \(L=0.1\lambda_{p}\), \(\Omega=0.299\) and \(\beta=2.392\times 10^{7}\)m\({}^{-1}\) in the UDL case, (**c,d**)\(f=0.7\), \(L=0.1\lambda_{p}\), \(\Omega=0.2452\) and \(\beta=0.5613\times 10^{7}\)m\({}^{-1}\) in UDL case, (**e,f**) \(f=0.7\), \(L=0.1\lambda_{p}\), \(\Omega=0.2942\) and \(\beta=0.6628\times 10^{7}\)m\({}^{-1}\) in the UML case. The yellow regions represent the metallic layers while the green regions represent the dielectric layers. The white regions on the left side of figures denote the partner dielectric medium. Here, \(\Gamma\) is the damping coefficient, which is assumed to be \(1\times 10^{11}\) s\({}^{-1}\) for low loss metallic material. In this condition, the imaginary part of the metallic permittivity only leads to a limited propagation length of the SWs, while corresponding dispersion relations calculated in the last section remain unchanged as the phase transition in the metallic layers is still insignificant [39]. A Gaussian beam with a beam waist \(G_{w}\) radiates obliquely from the prism to the partner medium, and the SWs can be excited in the vicinity of the interface separating the MMDM and the partner medium (at the \(x\)-\(z\) plane in Figure 6) by properly modulating \(h\).