1. Introduction

In recent years, the presence of hyperbolic materials brought about a new and attractive field [1–4]. These materials have been paid great attention in optics and relevant technology fields due to their fascinating optical characteristics and applications, such as nano-imaging [5–7], sensing [8–9], guiding-waves [10–11], thermal conductivity and emission [12–14], and surface phonon polaritons [15–20], as well as guiding and confinement of electromagnetic waves [10–11,21]. Especially, the surface phonon polaritons supported by naturally hyperbolic crystals can propagate for a long distance since very small optical-phonon damping or optical absorption [15,22], unlike surface plasmon polaritons in metals or metal metamaterials [23–25]. Naturally hyperbolic crystals are very perfect examples for theoretical investigations, which have lower optical loss and whose physical parameters have been determined [2,21,26]. The most typical uniaxial hyperbolic crystal is the hexagonal boron nitride (hBN) [21,27–28]. It has the two separated reststrahlen frequency bands (RFBs), or RFB-I corresponding to the longitudinal principal-value of permittivity and RFB-II related to the transverse principal-values. The hBN is one hyperbolic material in the RFBs and is one ellipsoidal material outside the RFBs. The α-MoO₃ is a typical biaxial hyperbolic crystal that has three RFBs corresponding to the three principal values of its permittivity tensor, respectively [29–32]. In contract to the hBN, the three RFBs of α-MoO₃ are not separated so that there are partly overlaps.

Conventional surface polaritons (CSPs) can be put into two categories. CSPs in the first category propagate along the relevant material surface and monotonously attenuate along the surface normal, and meanwhile they are either TM surface waves whose magnetic-field normal to the propagation plane or TE surface waves whose electric-field vertical to the propagation plane [33–34]. For example, surface magnon polaritons supported by magnetic crystals in the Voigt geometry are TE surface waves [33,35–37], but surface plasmon polaritons at the metal surface and surface phonon polaritons supported by dielectric crystals are TM surface waves [34,38–39]. CSPs in the second category are of hybrid-polarization, or a mixture of TE and TM surface waves. Hybrid-polarization surface polaritons generally consist of two branch waves in the supporting materials. The Dyakonov-like surface polaritons [17,40–45] are typical hybrid-polarization surface polaritons, whose two branch waves are TM and TE waves in the supporting media or materials, respectively. For conventional surface polaritons, two very important propagation constants characterize them. One is the wavevector aligned along the surface and the other is the attenuation constant showing the exponential decay of surface-polariton fields with the distance from the surface, where the attenuation constant is a positive real quantity.

For a complex attenuation constant, it is not easy to find surface-wave solutions from the wave equation under the electromagnetic boundary conditions. However, surface polaritons with complex attenuation constant was predicted [46] in an ionic-crystal/dielectric metamaterial and was experimentally observed [47] at the surface of bulk calcite recently, where the optical axis of the material lies in the propagation plane and takes an angle with respect to the surface normal. The surface polariton was called the ghost surface polariton (GSP) and is a TM hyperbolic surface wave. The most difference between the CSPs and GSP is that the CSP electromagnetic fields are exponentially attenuate with the distance away from the surface, but the GSP electromagnetic fields not only attenuate but also oscillate with the distance. The anisotropy and oblique optical axis are the necessary conditions for the existence of GSP. More recently, another GSP was predicted at the surface of antiferromagnets and ferromagnets in an external magnetic field [44,48], which is a hybrid-polarization surface wave and whose two branch waves in the supporting magnetic crystals are coherent. This GSP results from the antiferromagnetic or ferromagnetic gyromagnetism, so the external magnetic field is its necessary condition.

The CSPs at the surface of naturally hyperbolic materials have been researched for many years [16,18,22,40–43] in the different geometries, where the principal axes are either normal to the surface or parallel to the surface. In this paper, our investigation will pay main attention to search ghost surface phonon polaritons supported by biaxial hyperbolic materials in the geometry where two principal axes are in the plane of propagation and the other lies in the surface-plane.

2. Theoretical mechanism

The geometry and coordinate system are shown in Fig. 1(a), where a biaxial hyperbolic material (BHM) occupies the space of $x > 0$ and the space of $x < 0$ is filled with air or vacuum. The two dielectric principal axes lie in the x-y plane (the plane of propagation) and the i-axis is at angle $\phi $ with respect to the surface normal. The l-axis is aligned along the z-axis. This geometry originates from the rotating transformation of the principal-axis framework around the l-axis (the z-axis), and effectively represents the three specific geometries (G-1, G-2 and G-3) as indicated by the right sketches in Fig. 1(a). The permittivity tensor of BHM is $\varepsilon = \textrm{diag}({\varepsilon _1},\; {\varepsilon _2},\; {\varepsilon _3}$) in the principal-axis framework, where ${\varepsilon _i}$ (i = 1, 2, 3) represents the three diagonal elements. Up to now, we have not seen these oblique surface geometries of BHMs in experiment.

 

Fig. 1. (a) The geometry and coordinate system, where three red bi-arrows represent the three principal-axis of BHM, and the i-axis is at an angle $\phi $ with respect to the surface normal and k is the wavenumber pointed along the z-axis (or the l-axis). This geometry represents the three specific geometries called G-i (i = 1,2 or 3), as shown by the right sketches. For the biaxial material α-MoO₃, one conventionally assumes that ${\varepsilon _1}$ is in [100], ${\varepsilon _2}$ is in [010] and ${\varepsilon _3}$ in [001] direction. (b) The frequency-dependence of principal values of the α-MoO₃ permittivity, where the black lines with bi-arrows indicate three reststrahlen frequency bands, respectively.

Download Full Size | PDF

Some of them may be easily realized, but some perhaps be uneasily achieved in experiment. Now, α-MoO₃ crystals of high-quality have obtained in experiment, whose size can reach to the order of cm [52]. It is one basis for one to realize different optical surfaces in experiment.

In the xyz coordinate system, the dielectric permittivity tensor of BHM is changed into a nondiagonal matrix

For the surface polariton, $\mathrm{\Gamma }^{\prime}$ and the real part of $\mathrm{\Gamma }$ must be positive. Therefore,

In the case of $D < 0$, ${\mathrm{\Gamma }_ \pm }$ are certainly two complex quantities. Letting ${\mathrm{\Gamma }_ \pm } = \alpha \pm i\beta $, we find

From the above discussion and analysis, we understand that the HSP is defined only in the case of $D > 0$ but the GSP can exist only in the case of $D < 0$, so $D = 0$ separates the HSP and GSP and determines the boundary between the HSP and GSP. If the dispersion curve of HSP is a finite segment, its either endpoint will correspond to ${\mathrm{\Gamma }_ - } = 0$ and ${\mathrm{\Gamma }_ + } = \sqrt { - a} $ with $a < 0$ or ${\mathrm{\Gamma }_ + } = {\mathrm{\Gamma }_ - } = \sqrt { - a/2} $ $({D = 0} )$ with $a < 0$. If the dispersion curve of GSP is a finite segment, its either endpoint will be related to $\beta = 0$ and $\alpha = \sqrt { - a/2} $ with $a < 0$ or $\alpha = 0$ and $\beta = \sqrt {a/2} $ with $a > 0$.

The three components of GSP or HSP electric field in the BHM couple together in the configuration in Fig. 1 and it is easy for us to find

The relevant magnetic amplitudes can be found to be

The existence of surface polariton naturally results in the dispersion relation since Eq. (11) must be satisfied. After eliminating the field amplitudes from Eq. (11), we find the dispersion relation to be

The dispersion Eq. (12) and field-amplitude relations ( Eq. 13) are available for both the GSP and HSP. However, Eq. (12) is a complex equation in the case of $D < 0$, so it will need to be further analyzed. The above results are obtained for $\phi \ne 0$ or $\pi /2$. For $\phi = 0$ or $\pi /2$, we can find only one TM surface-wave solution with

We see from the above analytical process and discussion that the electric field of the GSP or HSP can be determined only if the z-component ($E{^{\prime}_z}$) of the field-amplitude above the surface is fixed. Thus, the electric field in the BHM is expressed as

The relevant magnetic field can be determined with equation ${\boldsymbol h} ={-} i\delta \nabla \times {\boldsymbol e}/\omega $. It is obvious that it can be expressed in the BHM as,

As indicated above, Eq. (12) is a complex equation in the case of $D < 0$. However, it can be changed into a real equation since ${\mathrm{\Gamma }_ + }$ and ${\mathrm{\Gamma }_ - }$ are mutually conjugate. Due to the unique relation between ${\mathrm{\Gamma }_ + }$ and ${\mathrm{\Gamma }_ - }$, ${\lambda _ - } = \lambda _ + ^\ast $, ${\gamma _ - } = \gamma _ + ^\ast $ and $\mathrm{\Lambda }{\mathrm{\Lambda }^\ast } = 1$. Therefore, the two terms on the left-handed side of dispersion relation (12) are mutually conjugate, so Eq. (12) is reduced to (see Appendix 1)

The phase difference between ${e_z}$ and ${e_x}$ or ${e_y}$ is ${\pm} \pi /2\; $ and that between ${h_z}$ and ${h_x}$ or ${h_y}$ also is ${\pm} \pi /2$. The phase relations between the electric and magnetic fields also are very simple. According to the phase relations among the amplitude components, we conclude that the Poynting vector ${\boldsymbol S} = 0.5Re({{{\boldsymbol e}^{\boldsymbol \ast }} \times {\boldsymbol h}} )$ of the HSP or GSP is parallel to the surface and is aligned along the z-axis.

We see the main features of the GSP as follows. (i) It is a unique hybrid-polariton surface wave. Their two branch-waves are coherent in the BHM and the corresponding amplitude-components are identical in amplitude but are different in primary phase, a couple of corresponding amplitude-components has their primary phase difference not equal to that of another couple. (ii) Its electromagnetic fields and Poynting vector oscillate and attenuate with the distance from the surface. (iii) No GSP exists for $\phi = 0$ or $\pi /2$, so the hyperbolic anisotropy and oblique orientation of principal-axes to the surface are the necessary conditions for the existence of GSP.

3. Numerical results and discussions

We take an α-MoO₃ crystal as an example to numerically stimulate properties of GSP in biaxial hyperbolic materials. The physical parameters included in Eq. (2) are given in are given in The chirality in twisted bilayer α-MoO3 was discussed [49] and the near-field radiative modulator driven by anisotropic hyperbolic polaritons in this biaxial hyperbolic material recently [50,51].

Table 1 [30–32], where the damping constants of optical phonons are ignored. The three diagonal elements of the primary permittivity change with frequency, as illustrated in Fig. 1(b). The three reststrahlen frequency bands are RFB-I $({545\textrm{c}{\textrm{m}^{ - 1}} < \omega < 851.7\textrm{c}{\textrm{m}^{ - 1}}} )$, RFB-II $({820\textrm{c}{\textrm{m}^{ - 1}} < \omega < 972.4\textrm{c}{\textrm{m}^{ - 1}}} )$ and RFB-III $({958\textrm{c}{\textrm{m}^{ - 1}} < \omega < 1004.1\textrm{c}{\textrm{m}^{ - 1}}} )$, which correspond to ${\varepsilon _i}({i = 1,\textrm{}2,\textrm{}3} )$, respectively. There are overlaps among the three RFBs. The chirality in twisted bilayer α-MoO₃ was discussed [49] and the near-field radiative modulator driven by anisotropic hyperbolic polaritons in this biaxial hyperbolic material recently [50–51].

Table 1. Optical parameters of α-MoO₃

View Table

For theoretical completeness, we propose that the i-principal-axis is at angle $\phi $ with respect to the surface normal, so we have three geometries, i.e., the G-i (i = 1, 2, or 3), as illustrated in Fig. 1(a). We first discuss the phase diagram of surface polaritons propagating in the x-z plane for $\phi = 45^\circ $, as illustrated in Fig. 2. In the three Geometries, we clearly see three enclosed regions surrounded by the boundary lines $({D = 0} )$ in the $\omega - k$ space, respectively. We find one GSP and one HSP in each Geometry. The GSP is situated inside the region and the HSP is localized outside the region. The dispersion curves of GSP and HSP form a continuous curve, which means that the transition exists between the GSP and HSP across the boundary. The HSP presents on the right-side of the region in the G-1 or G-3, but we find it on both sides of the region in the G-2. Of course, the bulk polariton continuums should present outside the reststrahlen frequency bands.

 

Fig. 2. The surface polariton phase-diagram in the wavenumber-frequency space for $\phi = 45^\circ $. The three regions surrounded by the dashed lines are the GSP regions, and the solid curves represent the GSP and the dotted curves show the HSP, where G-i (i = 1, 2, or 3) indicates the three geometries. The two thin lines are the ATR scanning lines will be used in Fig. 7.

Download Full Size | PDF

The attenuation and oscillatory constants are important for the GSP. The attenuation constant determines the localization of GSP at the surface and the oscillatory constant controls the oscillatory behavior of GSP with x. The localization of the HSP is determined by its two attenuation constants. We apply Fig. 3 to illustrate the constants of GSP and HSP existing in Fig. 2. We find some differences among the three GSPs. In the G-1, the GSP starts at the point where $\alpha = 0$ and $\beta $ is very large. It terminates at the point related to $\beta = 0$ and $\alpha = {\mathrm{\Gamma }_ + } = {\mathrm{\Gamma }_ - }$. Its terminal is just the starting point of HSP. The terminal of HSP is at the point related to ${\mathrm{\Gamma }_ - } = 0$. In the G-2, the starting point and terminal point of GSP correspond to the two different points on the boundary $({D = 0} )$ where $\beta = 0$. The left HSP starts from the point with ${\mathrm{\Gamma }_ - } = 0$ and terminates at the point where $\beta = 0$, or at the starting point of GSP. The right HSP starts from the terminal of GSP and stretches to the right. In the G-3, the GSP starts from one point on the boundary and terminates at another point on the boundary. The attenuation constant is obviously larger than the oscillatory constant. The HSP appears on the right-side of the GSP region and stretches to the right.

 

Fig. 3. The attenuation and oscillatory constants of GSP or HSP in Fig. 2 for ф=45°. (a) in the G-1, (b) in the G-2 and (c) in the G-3, where the dotted curves represent the attenuations ${\mathrm{\Gamma }_ + }$ (the upper) and ${\mathrm{\Gamma }_ - }$ (the lower) of HSP, and the solid and dashed curves show the attenuation and oscillatory constants of GSP, respectively.

Download Full Size | PDF

The orientation of the i-principal-axis with respect to the surface normal is a key role for the existence of GSP. The GSP degenerates into a CSP with TM polarization when the orientation angle $(\phi )$ is decreased to $0^\circ $ or is increased to $90^\circ $. Figure 4 illustrates the $\phi $- dependence of GSP and HSP. Figure 4(a) indicates that the two surface polaritons in the G-1 are situated in the RFB-III for various values of $\phi $. The dispersion curves of GSP cross over mutually. For an smaller value of $\phi $, the GSP curve is sandwiched by the two HSP curve segments, and these curves form a continuous dispersion curve, as indicated by the black curves. The dispersion curve of GSP either lie between the two CSP curves obtained from Eqs. (13) and (14) or lie outside of the two CSP curves. Figure 4(b). illustrates the attenuation and oscillatory constants of GSP in the G-1. It demonstrates that the starting point and terminating point both correspond to $\mathrm{\beta } = 0$ for small values of $\phi $ or are on the region boundary of $\textrm{D} = 0$, but the two points correspond to $\mathrm{\alpha } = 0$ and $\mathrm{\beta } = 0$ for the two larger values of $\phi $, respectively. For a larger $\phi $, $\mathrm{\beta }$ is evidently larger than $\mathrm{\alpha }$, as indicated by the two blue curves. Figure 4(c) and (d) show the dispersion curves and attenuation and oscillatory constants of the GSP in the G-2. The GSP and HSP are situated in the RFB-I and below the CSPs. The HSP is localized on both sides of the GSP and their dispersion curves structure a smooth curve. All the dispersion curves of GSP start from a point on the boundary of $\textrm{D} = 0$ and terminate at another on the boundary. For the two larger values of $\phi $, $\mathrm{\beta }$ is much larger that $\mathrm{\alpha }$. Figure 4(e) and (f) show the dispersion curves, the attenuation and oscillatory constants of the GSP in the G-3. They demonstrate that the GSP and HSP appear in the RFB-II and are sandwiched between the two CSPs. Similar to the GSP in the G-1, the beginning point of the dispersion curve is related to $\mathrm{\alpha } = 0$ and the terminating point corresponds to $\mathrm{\beta } = 0$ for larger values of $\phi $. Unlike the GSP in the G-1 and G-2, the curves of attenuation and oscillatory constants are more complicated, which results from the overlap of the RFB-II and RFB-I. In general, the oscillatory constant can be much larger than the attenuation constant for a bigger orientation angle, so we can see the obviously oscillatory phenomenon of GSP in this situation.

 

Fig. 4. Dispersion curves of GSP and curves of $\mathrm{\alpha }$ and $\mathrm{\beta }$ for various values of $\phi $ in the three geometries. (a) and (b) dispersion curves and $\mathrm{\alpha }$- and $\mathrm{\beta }$-curves in the G-1. (c) and (d) those in the G-2. (e) and (f) those in the G-3. The cross-curves are the dispersion curves of CSPs.

Download Full Size | PDF

Subsequently, we numerically explore the polarization and energy-flux density distribution of the GSP at some specific points on the dispersion curves. The Poynting vector of a surface polariton reflects the surface-polariton energy propagation and intensity distribution in the space and the electric field indicates its polarization. We select some specific points on the GSP dispersion curves in Figs. 4 to exhibit behaviors of the electric field and energy-flux density distribution.

Figures 5(a) and (b) show the GSP electric field at the specific points on the dispersion curves and the Poynting vector on the blue solid curve in Fig. 4(a). We firstly find that the oscillatory phenomenon is too slight to be seen for $\phi = 15^\circ $ since Fig. 4(b) demonstrates that $\mathrm{\beta }$ is evidently smaller than $\alpha $, and meanwhile the energy-flux density mainly distributes outside the $\alpha $-MoO₃ crystal or above the surface. In this case, the electric field of GSP at the surface is approximatively lies in the plane of propagation. We secondly see the obviously oscillatory behavior of GSP in the $\alpha $-MoO₃ for $\phi = 75^\circ $. In the present case, the interferent fringes of two branch-waves in the $\alpha $-MoO₃ are clearly reflected in energy-flux density by Fig. 5(b), where the attenuation constant is much smaller than the oscillatory constant, as shown in Fig. 4(b). Therefore, we see an evidently oscillatory phenomenon in electric-field and energy-flux density for $\phi = 75^\circ $, but it is difficult to clearly identify this phenomenon for $\phi = 15^\circ $ in the G-1. In addition, the Poynting vector always is positive, or along the z-direction.

 

Fig. 5. The x-dependence of the GSP electric field and the energy-flux density distribution for $\phi = 75^\circ $ in the G-1, where x = 0 represents the position of the surface. (a) The x-dependence of the electric field at the specific point on the black and blue dispersion curves in Fig. 3(a). (b) that energy-flux density distribution on the blue dispersion curve in Fig. 3(a), where the energy-flux density is measured in ${|{E{^{\prime}_z}} |^2}$.

Download Full Size | PDF

Figure 6 illustrates the spatial distribution of GSP electric field and the relevant energy-flux density in the G-2. Figure 6(a) shows the distribution of electric field at the specific point on the two dispersion curves in the G-2. At the point $({\phi = 15^\circ } )$ on the black curve in Fig. 4(c), the GSP has the electric field approximately parallel to the plane of propagation above the surface, but its electric field is roughly pointed along the z-axis and chiefly exponentially decay with x beneath the surface. At the point $({\phi = 75^\circ } )$ on the cyan curve in Fig. 4(c), the oscillatory phenomenon is directly seen in electric field, where $\mathrm{\beta }$ is much larger than $\mathrm{\alpha }$. In order to elaborately observe the x-dependence of GSP energy-flux density, Fig. 6(b) illustrates the Poynting vector distribution on the cyan dispersion curve in Fig. 4(c), which is displayed in logarithm to base 10. The energy-flux density always is positive and further it exponentially attenuates and sinusoidally oscillates with x in the α-MoO₃. The interferent fringes in the α-MoO₃ are clearly seen, as shown in Fig. 6(b). In general, it is of attenuation-oscillation in the α-MoO₃.

 

Fig. 6. The electric field and energy-flux density as functions of x in the G-2. x = 0 corresponds to the position of the surface. (a) the distribution of electric field for the fixed wavenumber of $k = 1600\; \textrm{c}{\textrm{m}^{ - 1}}\; $ and two specific angles. (b) the x-dependence of energy-flux density for the fixed angle of $\phi = 45^\circ$ and two given wavenumbers, where the energy-flux density is measured in ${|{E{^{\prime}_z}} |^2}$.

Download Full Size | PDF

We have not offered the distribution of the GSP electric-field in the G-3 since its performance is qualitatively similar to that in the G-2. The electric field of HSP consists the two branch-electric-fields in the α-MoO₃, where either branch-field exponentially decays with the distance from the surface, but the attenuations of the two branch-fields are out of sync. In general, the electric field of HSP approximatively exhibits an exponential-attenuation feature. The attenuation-oscillation picture of the GSP Poynting vector results from the interference between the two branch-waves in the α-MoO₃ and is seen as interferent fringes in the plane vertical to the propagation direction, as shown in Figs.  5(b) and  6(b).

Finally, we numerically simulate the GSPs based on the attenuated total-reflection (ATR) technique to prove their existence and observability. The experimental set comprises a Si-prism (${\varepsilon _i} = 11.56)$, an air spacer and a substrate (α-MoO₃). An incident radiation is partly reflected on the interface between the prism and air-spacer (see Fig. A in Appendix 3). When the incident angle (θ) is larger than the critical angle (${\theta _c} = {17.1^o}$), the incident radiation is totally reflected and produces an evanescent wave in the spacer. However, the reflective radio can be obviously reduced when a surface polariton is excited in the substrate, and then the reflective spectrum exhibits a sharp dip corresponding to the excited surface polariton. In our numerical stimulation, the damping term must be considered since it produces optical absorption. We first use a transfer-matrix method to obtain the expression of reflective ratio (this process is put in Appendix 3) and then numerically calculate the reflective ratio (or the ATR spectrum). The TM-polarized incidence is used to excite the GSPs in the substrate and the two frequency-scanning lines stretch across the G-2 and G-3 enclosed regions, respectively, as shown in Fig. 2. We see from Fig. 7 that the two sharp dips accurately correspond to the two intersections of the scanning lines and GSP dispersion curves in Fig. 2. It proves that the GSPs can be excited and is observable based on the attenuated total reflection technique.

 

Fig. 7. The ATR spectra along the two scanning lines shown in Fig. 2, where the two lines correspond to incident angles $\theta = {45^o}\; and\; {25^o}$, respectively. $\tau $ is the damping constant and d is the thickness of the air spacer.

Download Full Size | PDF

4. Conclusion

We have found a new-type of surface phonon polaritons with hybrid-polarization in biaxially hyperbolic materials (BHMs), called ghost surface polariton (GSP). It consists of two coherent branch-waves in the BHMs and is a hybrid-polarization surface wave. Due to the biaxial anisotropy, the GSP appears in three enclosed regions in the wavenumber-frequency space, which corresponds to the three different geometries (G-i, i = 1, 2, or 3). The numerical results based on the α-MoO₃ crystal reflect that the three regions of GSP existence are situated in the three reststrahlen frequency bands of the α-MoO₃, respectively. An ordinary hybrid-polarization surface polariton (HSP) is seen outside each GSP-region. It is noteworthy that the GSP and HSP occupy different segments of one continuous dispersion curve. The orientation of the crystal-axis in each geometry relative to the α-MoO₃ surface is important for the existence of GSP. The transition between the GSP and HSP can be found as the wavenumber or orientation of the crystal-axis is changed. For the GSP, a series of interferent fringes in the plane vertical to the propagation direction can be seen in Poynting vector and electric field. This phenomenon is produced by the interference between the two branch-waves in the α-MoO₃. Despite the low-symmetric configuration in Fig. 1, the dispersion property of the GSP or HSP is reciprocal, i.e., $\omega (k )= \omega ({ - k} )$. The dependence of GSP on the orientation of crystal-axis is different in the different geometries.

The origination of the present GSP is different from that of the ghost surface magnon polariton found in antiferromagnets, i.e., the present GSP orientates from the hyperbolic anisotropy of BHMs and the used asymmetry geometries. But the GSP in the antiferromagnets results from the antiferromagnetic gyromagnetism or external magnetic field [48]. We think that such a kind surface polaritons may exist in other hyperbolic materials and it is worthy to further search them in different hyperbolic materials. This work opens new frequency-wavenumber windows of surface phonon polaritons and is helpful for surface-polariton applications.

Appendix 1: Derivation of dispersion equation of GSP

For a GSP, dispersion Eq. (12) in the text, or

(19)$${\gamma _ + }(\Gamma ^{\prime} + {\Gamma _ + })[({\Gamma _ - } - k{\lambda _ - }) - {\omega ^2}] - {\gamma _ - }(\Gamma ^{\prime} + {\Gamma _ - })[({\Gamma _ + } - k{\lambda _ + }) - {\omega ^2}] = 0,$$

is a complex relation. However, we can prove that this complex equation can be changed into a real equation. We first see that ${\mathrm{\Gamma }_ + }$ and ${\mathrm{\Gamma }_ - }$ are a pair of complex quantities with the same positive real part and are mutually conjugate (${\mathrm{\Gamma }_ - }$=$\mathrm{\Gamma }_ + ^\ast $). On this basis, we realize from Eq. (9) in the text that

(20)$${\lambda _ - } = \lambda _ + ^*,\,{\gamma _ - } = \gamma _ + ^*.$$

Appling the two relations and ${\mathrm{\Gamma }_ - } = \mathrm{\Gamma }_ + ^\ast $ to Eq. (19), we fortunately find that the second term is the conjugation of the first term on the left-handed side of Eq. (19). Therefore, dispersion relation (A1-1) is just reduced into

(21a)$${\mathop{\rm Im}\nolimits} \{ {\gamma _ + }(\Gamma ^{\prime} + {\Gamma _ + })[({\Gamma _ - } - k{\lambda _ - }) - {\omega ^2}]\} = 0,$$

or equally

(21b)$${\mathop{\rm Im}\nolimits} \{ {\gamma _ - }(\Gamma ^{\prime} + {\Gamma _ - })[({\Gamma _ + } - k{\lambda _ + }) - {\omega ^2}]\} = 0.$$

The two equations are equal for GSP. (21a) is the dispersion equation of GSPs in the text, or Eq. (17) in the text. According to the expression of Λ in the text, i.e.,

(22)$$\mathrm{\Lambda } ={-} {\gamma _ + }({\mathrm{\Gamma^{\prime}} + {\mathrm{\Gamma }_ + }} )/{\gamma _ - }({\mathrm{\Gamma^{\prime}} + {\mathrm{\Gamma }_ - }} ),$$

we easily find from (A1-2) and ${\mathrm{\Gamma }_ - } = \mathrm{\Gamma }_ + ^\ast $ that

(23)$${\mathrm{\Lambda }^\ast } ={-} {\gamma _ - }({\mathrm{\Gamma^{\prime}} + {\mathrm{\Gamma }_ - }} )/{\gamma _ + }({\mathrm{\Gamma^{\prime}} + {\mathrm{\Gamma }_ + }} )= {\mathrm{\Lambda }^{ - 1}},$$

Therefore,

(24)$$\mathrm{\Lambda }{\mathrm{\Lambda }^\ast } = 1,$$

which is useful for us to achieve the expressions of GSP electric field.

Appendix 2: Expressions of GSP electric field

We realize from the analytical process and discussion above Eq. (15) in the text that the electric field of the surface polaritons can be determined only if the z-component ($E{^{\prime}_z}$) of the field-amplitude above the surface is fixed. In the BHM, the electric field is expressed as [see Eqs. (15(a)) and (15(b)) in the text]

(25)$$\scalebox{0.9}{$\displaystyle{\mathbf e} = \frac{{{{E^{\prime}}_z}}}{{1 + \Lambda }}\left( {\begin{array}{{ccc}} {i{\lambda_ + }{e^{ - 2\pi {\Gamma _ + }x}} + i{\lambda_ - }\Lambda {e^{ - 2\pi {\Gamma _ - }x}},}&{i{\gamma_ + }{e^{ - 2\pi {\Gamma _ + }x}} + i{\gamma_ - }\Lambda {e^{ - 2\pi {\Gamma _ - }x}},}&{{e^{ - 2\pi {\Gamma _ + }x}} + \Lambda {e^{ - 2\pi {\Gamma _ - }x}}} \end{array}} \right){e^{2\pi i(kz - c\omega t)}}$},$$

and the electric field above the BHM is

(26)$${\mathbf e^{\prime}} = {E^{\prime}_z}{e^{2\pi \Gamma ^{\prime}x}}\left[ {\begin{array}{{ccc}} { - ik/\Gamma^{\prime},}&{i({\gamma_ + } + {\gamma_ - }\Lambda )/(1 + \Lambda ),}&1 \end{array}} \right]{e^{2\pi i(kz - c\omega t)}}.$$

Due to ${\mathrm{\Gamma }_ \pm } = \alpha \pm i\beta $ for a GSP, Eq. (25 ) is reduced into

(27)$$\scalebox{0.9}{$\displaystyle{\mathbf e} = \frac{{{{E^{\prime}}_z}{e^{ - 2\pi \alpha x}}}}{{1 + \Lambda }}\left( {\begin{array}{{ccc}} {i{\lambda_ + }{e^{ - 2\pi i\beta x}} + i{\lambda_ - }\Lambda {e^{2\pi i\beta x}},}&{i{\gamma_ + }{e^{ - 2\pi i\beta x}} + i{\gamma_ - }\Lambda {e^{2\pi i\beta x}},}&{{e^{ - 2\pi i\beta x}} + \Lambda {e^{2\pi i\beta x}}} \end{array}} \right){e^{2\pi i(kz - c\omega t)}}$}.$$

And further, assuming ${R_z} = real\; \{{1/({1 + \mathrm{\Lambda }} )} \}$ and ${I_z} = imag\; \{{1/({1 + \mathrm{\Lambda }} )} \}$ and using (A1-5), we find that

(28)$${e_z} = {E^{\prime}_z}{e^{ - 2\pi \alpha x}}\{ ({R_z} + i{I_z})[\cos (2\pi \beta x) - i\sin (2\pi \beta x)] + ({R_z} - i{I_z})[\cos (2\pi \beta x) + i\sin (2\pi \beta x)]\} {e^{2\pi i(kz - c\omega t)}},$$

which can be directly reduced into equation (18a) in the text, i.e.,

(29)$${e_z} = 2{E^{\prime}_z}{e^{ - 2\pi \alpha x}}\{ {R_z}\cos (2\pi \beta x) + {I_z}\sin (2\pi \beta x)\} {e^{2\pi i(kz - c\omega t)}}.$$

For the other two components ${e_x}$ and ${e_y}$, due to ${\lambda _ - } = \lambda _ + ^\ast $ and ${\gamma _ - } = \gamma _ + ^\ast $, as well as the relations between $E_{x,y}^ \pm $ and $E_z^ \pm $, we also can obtain the expressions of the two components. For example, assuming ${R_x} = real\; \{{{\lambda_ + }/({1 + \mathrm{\Lambda }} )} \}$ and ${I_x} = imag\; \{{{\lambda_ + }/({1 + \mathrm{\Lambda }} )} \}$, we have

(30)$${e_x} = i{E^{\prime}_z}{e^{ - 2\pi \alpha x}}\{ ({R_x} + i{I_x})[\cos (2\pi \beta x) - i\sin (2\pi \beta x)] + ({R_x} - i{I_x})[\cos (2\pi \beta x) + i\sin (2\pi \beta x)]\} {e^{2\pi i(kz - c\omega t)}},$$

which can be directly simplified to (18b) in the text, or

(31)$${e_x} = 2i{E^{\prime}_z}{e^{ - 2\pi \alpha x}}\{ {R_x}\cos (2\pi \beta x) + {I_x}\sin (2\pi \beta x)\} {e^{2\pi i(kz - c\omega t)}}.$$

Similarly, we can obtain the y-component to be

(32)$${e_y} = 2i{E^{\prime}_z}{e^{ - 2\pi \alpha x}}\{ {R_y}\cos (2\pi \beta x) + {I_y}\sin (2\pi \beta x)\} {e^{2\pi i(kz - c\omega t)}},$$

where ${R_y} = real\; \{{{\gamma_ + }/({1 + \mathrm{\Lambda }} )} \}$ and ${I_y} = imag\; \{{{\gamma_ + }/({1 + \mathrm{\Lambda }} )} \}$.

Appendix 3: Derivation of reflective ratio in the ATR method

In this appendix, we will present the expression of reflective ratio in the ATR technique in the Otto configuration. As shown in Fig. 8, a prism with a relative dielectric constant ${\varepsilon _i}$ and the a-MoO₃ substrate are separated by an air-spacer with thickness d. The incident radiation with incident angle $\theta $ illuminates on the interface between the prism and spacer and then is reflected. When the incident angle is larger than the critical angle, an evanescent wave exists in the spacer, which will couple with the relevant GSP.

 

Fig. 8. The Otto configuration for the calculation of attenuated total reflection, where I shows the incident radiation and R represents the reflected radiation and reflective ratio. The excited surface polariton propagates along the z-axis and d is the thickness of the air spacer.

Download Full Size | PDF

The electric field of incident wave is written as $({E_x^i,E_y^i,E_z^i} )$ and that of the corresponding reflected wave is expressed as $({E_x^r,E_y^r,E_z^r} )$ on the upper interface, where $E_x^i ={-} kE_z^i/{k_x} \textrm{and}E_x^r = kE_z^r/{k_x}$ with $k = \sqrt {{\varepsilon _i}} \omega sin\theta $ and ${k_x} = \sqrt {{\varepsilon _i}} \omega cos\theta $. The tangential components of the corresponding magnetic fields to be used are written as

(33)$$H_y^i = - \frac{{\omega \delta }}{{{k_x}}}E_z^i,\,H_y^r = \frac{{\omega \delta }}{{{k_x}}}E_z^r,\,H_z^i = \frac{{{k_x}\delta }}{\omega }E_y^i,\,H_z^r = - \frac{{{k_x}\delta }}{\omega }E_y^r.$$

Further, we assume the tangential components of the electric field in the air spacer to be

(34a)$${E_\textrm{y}} = A{e^{2\pi {\alpha _\textrm{0}}x}} + B{e^{ - 2\pi {\alpha _\textrm{0}}x}},$$

(34b)$${E_\textrm{z}} = C{e^{2\pi {\alpha _\textrm{0}}x}} + D{e^{ - 2\pi {\alpha _\textrm{0}}x}},$$

and then the tangential components in the air spacer are given by

(35a)$${H_\textrm{y}} = \frac{{\omega \delta }}{{i{\alpha _0}}}({C{e^{2\pi {\alpha_0}x}} - D{e^{ - 2\pi {\alpha_0}x}}} ),$$

(35b)$${H_\textrm{z}} = \frac{{{\alpha _0}\delta }}{{i\omega }}({A{e^{2\pi {\alpha_0}x}} - B{e^{ - 2\pi {\alpha_0}x}}} ),$$

where the common factor $\textrm{exp}[{2\pi i({kz - c\omega t} )} ]$ have been ignored, and ${\alpha _0} = \sqrt {{k^2} - {\omega ^2}} $. The electric and magnetic fields in the substrate α-MoO₃ have been expressed in the text. Subsequently, we use the electromagnetic boundary conditions that the tangential components of electric and magnetic fields are continuous at the two interfaces (x = 0 and x = d) to derive the reflected electric field from the incident electric field.

At the upper interface (x = 0), we have

(36a)$$E_y^i + E_y^r\textrm{ = }A + B,$$

(36b)$$E_z^i + E_z^r\textrm{ = }C + D,$$

(36c)$$i{k_x}({E_y^i - E_y^r} )\textrm{ = }{\alpha _0}({A + B} ),$$

(36d)$$\frac{{i{\varepsilon _i}}}{{{k_x}}}({E_z^i - E_z^r} )\textrm{ = }\frac{1}{{{\alpha _0}}}({C - D} ),$$

These equations lead to the following matrix equations i.e.,

(37a)$$\left( {\begin{array}{{c}} A\\ B \end{array}} \right) = \frac{1}{2}\left( {\begin{array}{{cc}} {1 + \frac{{i{k_x}}}{{{\alpha_0}}}}&{1 - \frac{{i{k_x}}}{{{\alpha_0}}}}\\ {1 - \frac{{i{k_x}}}{{{\alpha_0}}}}&{1 + \frac{{i{k_x}}}{{{\alpha_0}}}} \end{array}} \right)\left( {\begin{array}{{c}} {E_y^i}\\ {E_y^r} \end{array}} \right) = \frac{1}{2}{S_1}\left( {\begin{array}{{c}} {E_y^i}\\ {E_y^r} \end{array}} \right),$$

(37b)$$\left( {\begin{array}{{c}} C\\ D \end{array}} \right) = \frac{1}{2}\left( {\begin{array}{{cc}} {1 - \frac{{i{\varepsilon_i}{k_x}}}{{{k_x}}}}&{1 + \frac{{i{\varepsilon_i}{k_x}}}{{{k_x}}}}\\ {1 + \frac{{i{\varepsilon_i}{k_x}}}{{{k_x}}}}&{1 - \frac{{i{\varepsilon_i}{k_x}}}{{{k_x}}}} \end{array}} \right)\left( {\begin{array}{{c}} {E_z^i}\\ {E_z^r} \end{array}} \right) = \frac{1}{2}{S_\textrm{2}}\left( {\begin{array}{{c}} {E_z^i}\\ {E_z^r} \end{array}} \right),$$

At the lower interface (x = d), we find

(38a)$$A{e^{2\pi {\alpha _0}d}} + B{e^{ - 2\pi {\alpha _0}d}} = i{\gamma _ + }E_z^ +{+} i{\gamma _ - }E_z^ -,$$

(38b)$$C{e^{2\pi {\alpha _0}d}} + D{e^{ - 2\pi {\alpha _0}d}} = E_z^ +{+} E_z^ -,$$

(38c)$${\alpha _0}({A{e^{2\pi {\alpha_0}d}} - B{e^{ - 2\pi {\alpha_0}d}}} )={-} i({{\gamma_ + }{\Gamma _ + }E_z^ +{+} {\gamma_ - }{\Gamma _ - }E_z^ - } ),$$

(38d)$$\frac{{{\omega ^\textrm{2}}}}{{{\alpha _0}}}({C{e^{2\pi {\alpha_0}d}} - D{e^{ - 2\pi {\alpha_0}d}}} )= ({\Gamma _\textrm{ + }} - k{\lambda _ + })E_z^ + \textrm{ + }({\Gamma _ - } - k{\lambda _ - })E_z^ -,$$

Similarly, these equations can be transformed into tw matrix equations as follows

(39a)$$\left( {\begin{array}{{c}} A\\ B \end{array}} \right) = \frac{i}{2}\left( {\begin{array}{{cc}} {{\gamma_\textrm{ + }}(1 - \frac{{{\Gamma _ + }}}{{{\alpha_0}}}){e^{ - 2\pi {\alpha_0}d}}}&{{\gamma_ - }(1 - \frac{{{\Gamma _ - }}}{{{\alpha_0}}}){e^{ - 2\pi {\alpha_0}d}}}\\ {{\gamma_\textrm{ + }}(1 + \frac{{{\Gamma _ + }}}{{{\alpha_0}}}){e^{2\pi {\alpha_0}d}}}&{{\gamma_ - }(1 + \frac{{{\Gamma _ - }}}{{{\alpha_0}}}){e^{2\pi {\alpha_0}d}}} \end{array}} \right)\left( {\begin{array}{{c}} {E_z^ + }\\ {E_z^ - } \end{array}} \right) = \frac{1}{2}{T_1}\left( {\begin{array}{{c}} {E_z^ + }\\ {E_z^ - } \end{array}} \right),$$

(39b)$$\left( {\begin{array}{{@{}c@{}}} C\\ D \end{array}} \right) = \frac{1}{2}\left( {\begin{array}{{@{}cc@{}}} {\left[ {1 + ({\Gamma _ + } - k{\lambda_ + })\frac{{{\alpha_0}}}{{{\omega^2}}}} \right]{e^{ - 2\pi {\alpha_0}d}}}&{\left[ {1 + ({\Gamma _ - } - k{\lambda_ - })\frac{{{\alpha_0}}}{{{\omega^2}}}} \right]{e^{ - 2\pi {\alpha_0}d}}}\\ {\left[ {1 - ({\Gamma _ + } - k{\lambda_ + })\frac{{{\alpha_0}}}{{{\omega^2}}}} \right]{e^{2\pi {\alpha_0}d}}}&{\left[ {1 - ({\Gamma _ - } - k{\lambda_ - })\frac{{{\alpha_0}}}{{{\omega^2}}}} \right]{e^{2\pi {\alpha_0}d}}} \end{array}} \right)\left( {\begin{array}{{@{}c@{}}} {E_z^ + }\\ {E_z^ - } \end{array}} \right) = \frac{1}{2}{T_\textrm{2}}\left( {\begin{array}{{@{}c@{}}} {E_z^ + }\\ {E_z^ - } \end{array}} \right),$$

It is evident that the relation between the incident and reflected electric fields is achieved from Eqs. (A-5a, A-5b, A-7a and A-7b) to be

(40)$$\left( {\begin{array}{{c}} {E_y^i}\\ {E_y^r} \end{array}} \right) = S_1^{ - 1}T_1^{}T_2^{ - 1}S_2^{}\left( {\begin{array}{{c}} {E_z^i}\\ {E_z^r} \end{array}} \right) = \Lambda \left( {\begin{array}{{c}} {E_z^i}\\ {E_z^r} \end{array}} \right),$$

which leads to

(41a)$$E_y^r = \frac{{{\Lambda _{22}}}}{{{\Lambda _{12}}}}E_y^i + ({\Lambda _{21}} - \frac{{{\Lambda _{11}}{\Lambda _{22}}}}{{{\Lambda _{12}}}})E_z^i,$$

(41b)$$E_z^r = \frac{1}{{{\Lambda _{12}}}}(E_y^i - {\Lambda _{11}}E_x^i),$$

and

(41c)$$E_x^r = kE_z^r/{k_x}.$$

Therefore, the reflective ratio $R = ({{{|{E_x^r} |}^2} + {{|{E_y^r} |}^2} + {{|{E_z^r} |}^2}} )/({{{|{E_x^i} |}^2} + {{|{E_y^i} |}^2} + {{|{E_z^i} |}^2}} )$.

Funding

Natural Science Foundation of Heilongjiang Province (ZD2009103).

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. P. Shekhar, J. Atkinson, and Z. Jacob, “Hyperbolic metamaterials: fundamentals and applications,” Nano Converg. 1(1), 14 (2014). [CrossRef]  

2. P. Li, I. Dolado, F. J. Alfaro-Mozaz, et al., “Infrared hyperbolic metasurface based on nanostructured van der Waals materials,” Science 359(6378), 892–896 (2018). [CrossRef]  

3. K. Korzeb, M. Gajc, and D. A. Pawlak, “Compendium of natural hyperbolic materials,” Opt. Express 23(20), 25406–25424 (2015). [CrossRef]  

4. A. Poddubny, I. Iorsh, P. Belov, et al., “Hyperbolic metamaterials,” Nat. Photonics 7(12), 948–957 (2013). [CrossRef]  

5. P. Li, I. Dolado, F. J. Alfaro-Mozaz, et al., “Optical Nanoimaging of Hyperbolic Surface Polaritons at the Edges of van der Waals Materials,” Nano Lett. 17(1), 228–235 (2017). [CrossRef]  

6. T. G. Folland, A. Fali, S. T. White, et al., “Reconfigurable infrared hyperbolic metasurfaces using phase change materials,” Nat. Commun. 9(1), 4371 (2018). [CrossRef]  

7. P. Li, M. Lewin, A. V. Kretinin, et al., “Hyperbolic phonon-polaritons in boron nitride for near-field optical imaging and focusing,” Nat. Commun. 6(1), 7507 (2015). [CrossRef]  

8. N. Vasilantonakis, G. A. Wurtz, V. A. Podolskiy, et al., “Refractive index sensing with hyperbolic metamaterials: strategies for biosensing and nonlinearity enhancement,” Opt. Express 23(11), 14329–14343 (2015). [CrossRef]  

9. Y. Guo, W. Newman, C. L. Cortes, et al., “Applications of hyperbolic metamaterial substrates,” Adv. Optoelectron. 2012, 1–9 (2012). [CrossRef]  

10. S. Dai, Q. Ma, T. Andersen, et al., “Subdiffractional focusing and guiding of polaritonic rays in a natural hyperbolic material,” Nat. Commun. 6(1), 6963 (2015). [CrossRef]  

11. H. L. Xu, X. Wang, X. Jiang, et al., “Guiding characteristics of guided waves in slab waveguide with hexagonal boron nitride,” J. Appl. Phys. 122(3), 033103 (2017). [CrossRef]  

12. J. Liu and E. Narimanov, “Thermal hyperconductivity: Radiative energy transport in hyperbolic media,” Phys. Rev. B 91(4), 041403 (2015). [CrossRef]  

13. S. A. Biehs, M. Tschikin, and P. Ben-Abdallah, “Hyperbolic metamaterials as an Analog of a blackbody in the near field,” Phys. Rev. Lett. 109(10), 104301 (2012). [CrossRef]  

14. Y. Guo, C. L. Cortes, S. Molesky, et al., “Broadband super-Planckian thermal emission from hyperbolic metamaterials,” Appl. Phys. Lett. 101(13), 131106 (2012). [CrossRef]  

15. J. D. Caldwell, L. Lindsay, V. Giannini, et al., “Low-loss, infrared and terahertz nanophotonics using surface phonon polaritons,” Nanophotonics 4(1), 44–68 (2015). [CrossRef]  

16. Z. W. Zhao, H. W. Wu, and Y. Zhou, “Surface-confined edge phonon polaritons in hexagonal boron nitride thin films and nanoribbons,” Opt. Express 24(20), 22930–22942 (2016). [CrossRef]  

17. Q. Zhang, S. Zhou, S. F. Fu, et al., “Rich hybridized-polarization surface phonon polaritons in hyperbolic dielectric metamaterials,” AIP Adv. 7(10), 105211 (2017). [CrossRef]  

18. S. Dai, M. Tymchenko, Y. Yang, et al., “Manipulation and steering of hyperbolic surface polaritons in hexagonal boron nitride,” Adv. Mater. 30(16), e1706358 (2018). [CrossRef]  

19. Y. Zhang, X. Wang, D. Zhang, et al., “Unusual spin and angular momentum of Dyakonov waves at the hyperbolic-material surface,” Opt. Express 28(13), 19205 (2020). [CrossRef]  

20. Xiang-Guang Wang, Shu-Fang Yu-Qi Zhang, Sheng Zhou, Fu, et al., “Goos-Hanchen and Imbert Fedorov shifts on hyperbolic crystals,” Opt. Express 28(17), 25048–25059 (2020). [CrossRef]  

21. J. D. Caldwell, A. V. Kretinin, Y. Chen, et al., “Sub-diffractional volume-confined polaritons in the natural hyperbolic material hexagonal boron nitride,” Nat. Commun. 5(1), 5221 (2014). [CrossRef]  

22. W. Ma, P. Alonso-González, S. Li, et al., “In-plane anisotropic and ultra-low-loss polaritons in a natural van der Waals crystal,” Nature 562(7728), 557–562 (2018). [CrossRef]  

23. J. B. Khurgin and A. Boltasseva, “Reflecting upon the losses in plasmonics and metamaterials,” MRS Bull 37(8), 768–779 (2012). [CrossRef]  

24. J. B. Khurgin and G. Sun, “In search of the elusive lossless metal,” Appl. Phys. Lett. 96(18), 181102 (2010). [CrossRef]  

25. J. B. Khurgin and G. Sun, “Scaling of losses with size and wavelength in nanoplasmonics and metamaterials,” Appl. Phys. Lett. 99(21), 211106 (2011). [CrossRef]  

26. A. K. Geim and I. V. Grigorieva, “Van der Waals heterostructures,” Nature 499(7459), 419–425 (2013). [CrossRef]  

27. S. Dai, Z. Fei, Q. Ma, et al., “Tunable phonon polaritons in atomically thin van der Waals crystals of boron nitride,” Science 343(6175), 1125–1129 (2014). [CrossRef]  

28. A. J. Giles, S. Dai, O. J. Glembocki, et al., “Imaging of anomalous internal reflections of hyperbolic phonon-polaritons in hexagonal boron nitride,” Nano Lett. 16(6), 3858–3865 (2016). [CrossRef]  

29. Z. Zheng, N. Xu, S. L. Oscurato, et al., “A mid-infrared biaxial hyperbolic van der Waals crystal,” Sci Adv 5(5), eaav8690 (2019). [CrossRef]  

30. C. Wei, S. Abedini Dereshgi, X. Song, et al., “Polarization Reflector/Color Filter at Visible Frequencies via Anisotropic α-MoO3,” Adv. Opt. Mater. 8(11), 2000088 (2020). [CrossRef]  

31. G. Hu, Q. Ou, G. Si, et al., “Topological polaritons and photonic magic angles in twisted α-MoO3 bilayers,” Nature 582(7811), 209–213 (2020). [CrossRef]  

32. S. Abedini Dereshgi, T. G. Folland, A. A. Murthy, et al., “Lithography-free IR polarization converters via orthogonal in-plane phonons in α-MoO3 flakes,” Nat Commun 11(1), 5771 (2020). [CrossRef]  

33. M. G. Cottam and D. R. Tilley, Introduction to Surface and Superlattice Excitations (IOP Publishing Ltd, Philadelphia, 2005) Chaps. 4 and 6.

34. V. M. Agranovich and D. L. Mills, “Surface Polaritons,” (North-Holland Publishing Company, Amsterdam, 1982) Chapter.1.

35. R. L. Stamps and R. E. Camley, “Green’s function for antiferromagnetic polaritons, I Surface modes and resonances,” Phys. Rev. B 53(18), 12232–12241 (1996). [CrossRef]  

36. R. E. Camley and D. L. Mills, “Surface-Polaritons on Uniaxial Antiferromagnets,” Phys. Rev. B 26(3), 1280–1287 (1982). [CrossRef]  

37. M. R. F. Jensen, T. J. Parker, K. Abraha, et al., “Experimental observation of surface magnetic polaritons in FeF 2 by attenuated total reflection (ATR),” Phys. Rev. Lett. 75(20), 3756–3759 (1995). [CrossRef]  

38. P. V. Ratnikov, “Surface Plasmon Polaritons in Planar Graphene Superlattices,” Phys. Rev. B 101(12), 125301 (2020). [CrossRef]  

39. S. Baher and Z. Lorestaniweiss, “Propagation of surface plasmon polaritons in monolayer graphene surrounded by nonlinear dielectric media,” J. Appl. Phys. 124(7), 073103 (2018). [CrossRef]  

40. E. Cojocaru, “Comparative analysis of Dyakonov hybrid surface waves at dielectric–elliptic and dielectric–hyperbolic media interfaces,” J. Opt. Soc. Am. B 31(11), 2558–2564 (2014). [CrossRef]  

41. Q. Zhang, S. Zhou, S. F. Fu, et al., “Diversiform hybrid-polarization surface plasmon polaritons in a dielectric–metal metamaterial,” AIP Adv. 7(4), 045216 (2017). [CrossRef]  

42. E. E. Narimanov, “Dyakonov waves in biaxial anisotropic crystals,” Phys. Rev. A 98(1), 013818 (2018). [CrossRef]  

43. S.-F. Fu, S. Zhou, Q. Zhang, et al., “Completely confinement and extraordinary propagation of Dyakonov-like polaritons in hBN,” Optics and Laser Technology 125, 106012 (2020). [CrossRef]  

44. H. Peng and X.-Z. Wang, “Unique surface polaritons and their transitions in metamaterials,” Opt. Express 30(12), 20883 (2022). [CrossRef]  

45. J.-X. Ta, Y.-L. Song, and X.-Z. Wang, “Magneto-phonon polaritons of antiferromagnetic/ion-crystal superlattices,” J. Appl. Phys. 108(1), 013520 (2010). [CrossRef]  

46. S. Zhou, S. F. Fu, Q. Zhang, et al., “Ghost surface phononic polaritons in ionic-crystal metamaterial,” J. Opt. Soc. Am. B 35(11), 2764 (2018). [CrossRef]  

47. W. Ma, G. Hu, D. Hu, et al., “Ghost hyperbolic surface polaritons in bulk anisotropic crystals,” Nature 596(7872), 362–366 (2021). [CrossRef]  

48. H. Song, X. Wang, Y. Zhang, et al., “Ghost surface polariton in antiferromagnets,” Phys. Rev. B 106(2), 024425 (2022). [CrossRef]  

49. B.-Y. Wu, Z.-X. Shi, F. Wu, et al., “Strong chirality in twisted bilayer α-MoO₃,” Chin. Phys. B 31(4), 044101 (2022). [CrossRef]  

50. Y. Hu, B. Wu, H. Liu, et al., “Near-field radiative modulator driven by anisotropic hyperbolic polaritons in biaxial hyperbolic materials,” J. Quant. Spectr. & Radiat. Transfer 296, 108468 (2023). [CrossRef]  

51. L. Li, X. Wu, H. Liu, et al., “Near-field radiative modulator based on α-MoO₃ films,” Int. J. Heat Mass Transfer 216, 124603 (2023). [CrossRef]  

52. W. Shen, Y. Yu, Y. Huang, et al., “Origins and cavity-based regulation of optical anisotropy of α-MoO3 crystal,” 2D Mater. 10(1), 015024 (2023). [CrossRef]