Van der Waals thin films of WTe2 for natural hyperbolic plasmonic surfaces

A hyperbolic plasmonic surface supports highly directional propagating polaritons with extremely large density of states. Such plasmon polaritons have been realized in artificially structured metasurfaces. However, the upper bound of the achievable plasmon wave vector is limited by the structure size, which calls for a natural hyperbolic surface without any structuring. Here, we experimentally demonstrate a natural hyperbolic plasmonic surface based on thin films of WTe2 in the light wavelength range of 16 to 23 microns by far infrared absorption spectroscopy. The topological transition from the elliptic to the hyperbolic regime is further manifested by mapping the isofrequency contours of the plasmon. Moreover, the anisotropy character and plasmon frequency exhibit prominent temperature dependence. Our study demonstrates the first natural platform to host 2D hyperbolic plasmons, which opens exotic avenues for the manipulation of plasmon propagation, light-matter interaction and light emission in planar photonics.

H yperbolic plasmonic surfaces, whose isofrequency contour in the wave vector space is a hyperbola, have been realized at visible 1 , mid-infrared 2 , and microwave frequencies 3 in metasurfaces, created by artificial subwavelength structuring or self-assembling carbon nanotubes. This enables a series of potential applications for planar photonics 4,5 , including but not limited to nanoscale imaging 6,7 , negative refraction 1 , and enhancement of spontaneous emission 8 . Since the observation of highly confined and tunable plasmons in graphene [9][10][11] , polaritons, especially plasmon polaritons 12 and phonon polaritons [13][14][15] in two-dimensional (2D) materials, have attracted substantial attention. Recently, 2D materials with in-plane anisotropic interplay of intraband and interband responses are predicted to naturally sustain hyperbolic plasmon polaritons (HPPs) 16,17 . In contrast to artificial surfaces, whose wave vectors in the hyperbolic regime are limited by the inverse of the structure size, natural hyperbolic surfaces support higher electromagnetic confinement and more diverging photonic density of states 18,19 . More importantly, the hyperbolic regime can be extended to mid-infrared and terahertz (THz) frequencies in 2D materials 19,20 . This range corresponds to energies of molecular vibrations and thermal radiation, promising opportunities in chemical sensing and thermal management.
While some anisotropic 2D materials, such as black phosphorous, are predicted to represent a natural class of HPPs 16,17 , the experimental demonstration of any natural hyperbolic plasmonic surface has not been realized yet. In-plane hyperbolic phonon polaritons, on the other hand, have been shown in natural MoO 3 surfaces 13,21 and structured hexagonal boron nitride metasurfaces 22 . In comparison with phonon polaritons, plasmon polaritons exhibit inherently stronger coupling to light 23 , making it more versatile for the management of thermal and spontaneous emission. Moreover, the hyperbolicity of plasmons in 2D materials can be tuned through chemical doping and gating 16 , enabling active in-plane manipulation of polariton propagation.
Semimetal WTe 2 is one of the layered materials with in-plane anisotropic electrodynamics 24,25 . Moreover, WTe 2 thin films host a wide range of remarkable electronic properties, such as extremely high mobility 26 , tunable carrier density by electrostatic gating [27][28][29] and Mo doping 30 , and being a candidate of type-II Weyl semimetals 31 , which greatly facilitate the study of 2D plasmons. Previous reflection measurements on bulk WTe 2 crystals revealed temperature-dependent anisotropic bulk plasma frequencies with extremely low optical scattering rates down to 0.25 cm -1 at 6 K 25,32 . Interestingly, the fitted in-plane dielectric functions exhibit a signchanging regime in the far-IR range 25 , indicating the potential for realizing hyperbolic 2D plasmons in the thin films. However, the experimental observation of anisotropic plasmon modes in WTe 2 thin films and the defining evidence for the hyperbolic surface, the hyperbolic isofrequency contour, are still absent. Especially, the hyperbolic regime is located at the far-IR range indicated by the bulk reflection data, where, on one hand, there is no readily available setup for near-field experiments, and on the other, largearea single-crystal films are needed for far-field measurements, posing challenges in the observation of 2D plasmons.
Here, anisotropic 2D plasmons and hyperbolic plasmon dynamics are successfully observed in exfoliated WTe 2 thin films using Fourier transform-infrared spectroscopy (FTIR), demonstrating the anisotropic 2D material as an important platform to realize natural hyperbolic plasmonic surfaces.

Results
Sample characterization and the scheme for transmission spectra measurements. As shown in Fig. 1a, WTe 2 is crystallized in a T d -type layered structure with a tungsten chain running along the a axis. To maintain the lattice orientation, thin films were obtained through mechanical exfoliation from a single crystal (see "Methods"). Before investigation of the plasmon, we first explored the far-IR absorption spectra of an unpatterned thin film with thickness of about 60 nm on Si/SiO 2 substrate, as marked by dashed lines in Fig. 1c. The top edge is along the a axis, which has been confirmed by Raman measurement (see Supplementary Fig. 1). Figure 1b shows the scheme for measuring absorption spectra, which is characterized by the extinction 1 − T/T 0 , where T and T 0 are the transmission of the light through the film and the bare substrate, respectively. The amplitude of the extinction spectra is determined by the sheet optical conductivity σ(ω) 33 (see Methods), which is typically captured by the following expression 18 : The conductivity is generally constituted of intraband (Drude response) and interband (bound states) electronic transitions in the form of the first and second terms on the right-hand side of Eq. (1), where index j = a, b, D, and S denote the spectral weights, ω b is the frequency of interband resonance, and Γ and η are the scattering widths for intraband and interband transitions, respectively. For simplicity, we retain only one Lorentz-type bound state in Eq. (1), which is sufficient for our discussion in the far-IR region. Electrons here are treated as normal carriers with Drude weight D = πn eff e 2 / m eff 34 , since the Fermi energy in undoped WTe 2 is well below the Weyl points 31 . Although two types of carriers 24,35 are reported in WTe 2 , to simplify the discussion, we treat them as one type with averaged effective carrier density (n eff ) and mass (m eff ).
Anisotropic optical absorption in WTe 2 thin films. Figure 1d displays the measured extinction spectra and the corresponding fitting curves of the WTe 2 film in Fig. 1c. We can see that larger Drude response (dotted dashed lines) can be observed along the a axis, while the interband component (dashed lines) is more intense along the b axis. The different anisotropy of conductivity between the intraband and interband transitions, which tends to facilitate the formation of a relatively broad hyperbolic frequency range (see Supplementary note II for details), is further manifested by the polarization dependence of the fitted intraband and interband spectral weight (Fig. 1e, f). The mid-infrared absorption spectra exhibit consistent anisotropy, as shown in Fig. 1j (see also Supplementary Fig. 2 for spectra in a wider range). The fitted conductivity is summarized in the inset of Fig. 1d. The signs of the imaginary part of the optical conductivity along two principal axes determine the topology of the plasmonic surface 16,36 . A frequency interval where σ″ aa > 0 (metallic electron response) and σ″ bb < 0 (dielectric response) can be found in the shaded area, giving a hyperbolic regime from 427 to 623 cm -1 . Using the fitted Drude weight, we can extract the ratio of effective mass by m eff b / m eff a = D aa /D bb = 2.05, consistent with the value of 2.2 measured in the bulk 25 . The temperature dependence of Fig. 1d was also investigated (see Supplementary Fig. 3). The Drude scattering width becomes broader at higher temperature (Fig. 1g), which has been attributed to Fermi liquid properties 32 . When the temperature increases, the Drude weights along both axes are enhanced due to more thermal carriers, however with different enhancement ratios (Fig. 1h), leading to a less anisotropic effective mass at higher temperature (Fig. 1i).
Anisotropic plasmon resonance modes in WTe 2 disk arrays. When WTe 2 films are patterned into microstructures such as microdisk arrays (right inset of Fig. 2a), localized plasmons can be excited. Note that such patterning is only for plasmon detection in the far field, rather than creating metasurfaces, which need artificial structures much smaller than the plasmon wavelength to induce effective hyperbolicity. In the long-wavelength limit, the 2D plasmon dispersion of free carriers is given by 34 where ε 0 is the vacuum permittivity, ε env is the dielectric constant of the surrounding environment, m eff depends on the q direction, and q is the effective wave vector determined by the structure size.
For ribbons with width L and disks with diameter d, the effective q is π/L and 3π/4d, respectively 37 . As shown in the lower inset of ig. 1 Infrared absorption spectra in unpatterned exfoliated WTe 2 thin films. a Schematic illustration of the layered crystal structure of WTe 2 ; the tungsten chains along the a axis are shown as yellow zigzag segments. b Schematic of the setup for extinction spectra measurements. c Optical microscope image of the exfoliated WTe 2 thin film on Si/SiO 2 substrate. Scale bar is 100 μm. The optical measurements were performed in the uniform area marked by dashed lines, with thickness of about 60 nm. d Far-IR absorption spectra of the unpatterned exfoliated WTe 2 thin film in (c) with polarization along a and b axes at 20 K. Solid, dotted dashed, and dashed lines are the corresponding fitted curves of the total extinction spectra, Drude component, and interband transition component, respectively. Inset: conductivity along a and b axes calculated by the fitting results. The shaded area represents the hyperbolic frequency regime. e, f Polarization dependence of the fitted Drude weight and interband weight. Curves are the fitting results by cos 2 θ. g, h Temperature dependence of the fitted scattering width and weight of Drude response along a and b axes. i Temperature dependence of the effective mass ratio calculated by the fitted Drude weight. j Mid-infrared absorption spectra of an unpatterned WTe 2 thin film exfoliated on a polycrystalline diamond substrate with polarization along a and b axes. The film thickness is~100 nm.
according to Eq. (2) (carrier density n eff is independent of direction). The effective mass ratio can be calculated by the resonance frequencies as m eff b /m eff a = (ω P a /ω P b ) 2 = 1.91, in agreement with Drude response results for the unpatterned film. Similar to the Drude response, the plasmon modes are tunable by temperature as well, as shown in Fig. 2b. The fitted resonance frequency, width, and the inferred mass ratio exhibit prominent temperature dependence (Fig. 2c-e). Two-peak structure can be observed at high temperature due to the hybridization with the polar phonon on SiO 2 substrates. More discussions on the temperature dependence and possible indication of Lifshitz transition can be found in Supplementary note IV and Supplementary  Fig. 5.
The hyperbolic regime derived from the anisotropic plasmon dispersion. In order to investigate the plasmon dispersion at higher energy, where the plasmonic resonances are determined by both the intraband and interband transitions, we fabricated a set of rectangle arrays along the two optical axes of WTe 2 films on polycrystalline diamond substrates (see Methods). With such substrate, we avoided plasmon hybridization with substrate-polar phonons. The films investigated on diamond substrates all have similar thickness of 100 ± 20 nm to maintain similar sheet optical conductivities. We first studied rectangle arrays with a fixed aspect ratio L a = 2L b to compensate the effect of anisotropic effective mass. The results of the extinction spectra along a and b axes are shown in Fig. 3a. A typical scanning electron microscopy (SEM) image of the rectangle array is shown in Fig. 3b. Two trends can be identified from Fig. 3a. The first trend is the plasmon dispersion. For large rectangles (8 × 4 μm 2 ), the plasmon modes along the two axes have nearly the same frequency (about 217 cm -1 ). However, when the size of the rectangles is reduced, the plasmon frequencies along both axes blueshift but at a different pace, leading to a splitting of the plasmon frequency (see Supplementary Fig. 7 for details). The second trend is the spectral weight. As the plasmon moves to higher energies, the plasmon intensity along both axes gets reduced (with filling factors and film thickness taken into consideration), while plasmons along the b axis have a larger intensity reduction rate (see Supplementary Fig. 7 for details).
The two trends can be further demonstrated by the plasmon spectra of a rectangle array with a higher aspect ratio in Fig. 3c. Although plasmons along the two optical axes have nearly identical frequency, the aspect ratio increases from 2 to 4, compared with the rectangle array (8 × 4 μm 2 ) in Fig. 3a (the first pair of spectra from the top), further confirming the saturation tendency of plasmon frequency along the b axis. In addition, we can see a clear decrease in the intensity ratio between b and a axes in the two arrays (from 0.4 for 8 × 4 μm 2 to 0.13 for 1.8 × 0.45 μm 2 ). The peak frequency of the latter array is closer to the hyperbolic regime (fitted below), which indicates a significant reduction of plasmonic spectral weight along the b axis near the topological transition boundary between the elliptic and the hyperbolic regimes.
To study the plasmon dispersion in detail, all the measured plasmon frequencies of WTe 2 rectangle arrays on diamond substrates are summarized in Fig. 4a. Detailed plasmon absorption spectra of the red and blue spheres can be found in Supplementary Fig. 6. The dashed lines represent the standard 2D plasmon dispersion with ω / ffiffi ffi q p . As shown in Fig. 4a, at low energy, the measured plasmon peaks follow the ffiffi ffi q p scaling; thus, it is legitimate for us to determine the effective mass anisotropy of free carriers from plasmons in Fig. 2. However, as the wave vector increases, the dispersion along both in-plane axes softens and departs from the purely free carrier case due to the coupling to interband transitions, and approaches an energy limit around 400 cm -1 (600 cm -1 ) for plasmons along the b (a) axis at large wave vectors. Because of the lower limiting energy, the dispersion along the b axis has a larger softening rate than that along the a axis, leading to the frequency separation in Fig. 3a.
To fit the plasmon dispersion, the loss function ÀIm 1=ε ð Þ, defined from the imaginary part of the inverse of the dielectric function, is calculated by substituting Eq. (1) into the dielectric function with the following form in the 2D case 34 (see Supplementary note I): The fitting result is plotted as a pseudocolor map in Fig. 4a, which agrees well with the measured plasmon dispersion (see Methods), as indicated by the good match between the solid curves and data points. Figure 4b displays the imaginary parts of the conductivity extracted from the fitting of Fig. 4a. The hyperbolic regime is located at the shaded area in the range of 429-632 cm -1 (15.8-23.3 microns in wavelength), fully consistent with the regime obtained through the optical conductivity of the unpatterned film in Fig. 1d. The energies of σ″ jj = 0, where j = a, b, are marked by the white dashed lines in Fig. 4a, which can be demonstrated to determine the energy boundaries of the plasmon dispersion along a and b axes (see Supplementary note II for details), implying that plasmon resonance modes can be found only along the a axis in the hyperbolic regime. The hyperbolic regime can be further confirmed by plasmonintensity evolution in Fig. 4c, where plasmons along the two optical axes have intensity severely suppressed at the corresponding boundaries of the hyperbolic regime. This is consistent with the prediction by the loss function, as shown by the intensity in the pseudocolor plot in Fig. 4a. A phenomenological coupled oscillator model, in which the plasmon and the bound state are represented by two coupled oscillators, is applied to give a quantitative description of the plasmon-intensity evolution (see Supplementary note V for details). As shown in Fig. 4c, the agreement is good, manifesting that plasmons with nonzero intensity cannot be found along the b axis in the hyperbolic regime (the shaded area). This is consistent with the topology of hyperbolic plasmon isofrequency contour, illustrated in the inset of Fig. 4a. The decreased plasmon weight is supposed to go to the interband plasmons at higher energies, as predicted by the coupled oscillator model (see Supplementary note V and Supplementary Fig. 8 for details). In fact, the coupling between intraband and interband transitions makes the bound states to have plasmon-like features. Such hybrid modes are termed interband plasmons 38 . Supplementary Fig. 10 shows the mid-infrared extinction spectra along the b axis of rectangle arrays with different L b length at 10 K. The frequencies of the two bound states shift to higher energies as the wave vector increases, providing an evidence for interband plasmons.
Topological transition. After establishing the energy range for hyperbolic plasmons, the isofrequency contours of the plasmon at different energies are shown in Fig. 5b to exhibit the topological transition from the elliptic to the hyperbolic regime. The points in Fig. 5b with wave vector q away from the directions of the principal axes were measured in ribbon arrays with a skew angle of θ with respect to the a axis (Fig. 5a), in which q = π/L (L is the ribbon width), with the direction perpendicular to the ribbon. The plasmon dispersion at each skew angle can be found in Supplementary Fig. 9. All the points in Fig. 5b can be fitted well by the solid curves based on the optical conductivity in Fig. 4b (see Methods). As shown in Fig. 5b, at energy below 429 cm -1 , which is in the elliptic regime, the isofrequency contour is elliptic with the long axis running along the b axis due to the smaller conductivity. At higher energy, the ellipse becomes flatter along the b axis, consistent with the increasing anisotropy of conductivities along the two principal axes, as shown in Fig. 4b. At energy above 429 cm -1 , where the imaginary parts of the conductivity have opposite signs along the two principal axes, the isofrequency contours in Fig. 5b become hyperbolic, which provide a defining evidence for the existence of a hyperbolic plasmonic surface.

Discussion
The fitted plasmonic resonance width in Fig. 2a is about 50 cm -1 at 10 K along both axes, and increased to about 400 cm -1 at room temperature due to the Fermi liquid properties 32 . If we assume no inhomogeneous broadening, the lifetime is about 0.1 ps (0.013 ps) at 10 K (300 K). The low-temperature lifetime is comparable to that of graphene plasmon on Si/SiO 2 substrate (about 0.05-0.1 ps) 39 , but much lower than the values of plasmon polaritons in hBNencapsulated graphene (about 1 ps) 40 and phonon polaritons in hBN (about 2 ps) 41 and MoO 3 (about 8 ps) 13 films. At low temperature, the plasmon resonance width becomes broader in the hyperbolic regime, with an average value of 100 cm -1 , due to the coupling with interband transitions. For WTe 2 films with thickness of 100 nm on diamond substrates, the calculated propagation length based on the plasmon lifetime has the maximal value of about 0.5 μm in the hyperbolic regime, and decreases at smaller group velocities. In fact, such large resonance width in WTe 2 films is inconsistent with the low damping rate of the bulk plasma (down to 0.25 cm -1 ) 32 and the high mobility (about 10,000 cm 2 V -1 s -1 ) 26 of bulk WTe 2 . This is possibly due to the inhomogeneous broadening that causes an underestimation of the plasmon lifetime in our study. Moreover, the carrier scattering from the surface impurities and the charge inhomogeneity of the substrates 42 can cause plasmon broadening as well. However, these can be potentially alleviated by hBN encapsulation, as demonstrated in graphene 40 . It is noted that although the hyperbolicity is demonstrated in the WTe 2 film with larger thickness (about 100 nm) than that of the widely studied monolayer 2D materials (e.g., graphene, MoS 2 ), it can be treated as a real 2D plasmonic system since the plasmon wavelength is much larger than the film thickness. In fact, it is in principle possible to observe the same hyperbolic plasmons down to about 12 nm, above which the electronic structure is reported to remain the same as that of the thick film 43 . For few-layer WTe 2 films, the band structures are dramatically modified due to finitesize effects, and it becomes a 2D topological insulator in the monolayer limit 44 . It is interesting to check how the hyperbolic plasmons develop in the few-layer limit in future study.
Due to the relatively large film thickness, the light confinement factor (wavelength ratio between light and plasmon) is relatively low at low frequencies (about 6 at~200 cm -1 along the b axis in Fig. 3a), while this ratio is increased to about 40 at 512 cm -1 (ribbon width 240 nm at a skew angle of 33°) within our measurements, which is comparable to the values in graphene plasmons (50-60 in ref. 11 and 40 in ref. 10 ). It should be noted that such high value is realized in the 100-nm-thick film with the sheet optical conductivity of over one order of magnitude higher than that in graphene, manifesting the high electromagnetic confinement in the hyperbolic regime. The light confinement can be further increased by reducing the film thickness, and reaches a maximum value of about 2300 at film thickness of~10 nm. However, such factor will be reduced after considering the losses 36 and nonlocality 17 , which might close the otherwise hyperbolic isofrequency contour at large momentum. At high temperature, the light confinement capability will be weakened due to the increased carrier density. The tunability is an unique characteristic for the hyperbolic plasmonic surface in 2D materials 16 . In fact, there have been many works reporting tunable electrical properties by electrostatic gating for WTe 2 films, some of which have a film thickness of above 10 nm 28,29 . In this work, the plasmons are observed in the WTe 2 film down to about 50 nm, and thinner films are needed in future work to enable the electrostatic gating. Moreover, the recently developed solid ion-gating technique 45 provides a chance to tune the carrier density for relatively thick films. Besides, the chemical doping method is another possible way to tune the WTe 2 plasmons. For example, the carrier density and effective mass have been reported to be tuned in Mo x W 1-x Te 2 single crystals by different Mo-doping contents 30 .
In conclusion, we use FTIR to demonstrate the existence of a natural hyperbolic plasmonic surface in exfoliated WTe 2 thin films in the far-IR range. The same kind of phenomena, in principle, could be observed in near-field spectroscopy 13,21,22 . However, due to the restrictions of laser wavelength and cryogenic sample conditions, such experiment is challenging at this stage 46 . Our successful observation of 2D plasmons in films obtained by mechanical exfoliation will stimulate investigations of hyperbolic plasmons at low or even room temperature in other anisotropic 2D materials, such as black phosphorus, T d -MoTe 2 , and PtTe 2 , and topological plasmons in other layered topological semimetals (e.g., ZrSiS and ZrTe 5 ).

Methods
Sample preparation and fabrication. WTe 2 single-crystal films were prepared by a standard mechanical exfoliation technique from bulk WTe 2 crystals (HQ Graphene) onto substrates. The typical film dimension is about 300 μm, greater than the beam size under an infrared microscope. The film thickness was determined by a Bruker Dimension Edge AFM system (probe model RTESP-300, tapping mode).
Here, all the films we investigated have thickness of above 50 nm, in which the possible surface oxidation has a negligible effect on the spectrum measurement 47 . Actually, no degradation of plasmonic devices was found during the nanofabrication and measurement process. Two types of substrates were used in this work. One is Si/SiO 2 substrate with SiO 2 thickness of 300 nm. The other is polycrystalline diamond grown by chemical vapor deposition (CVD) method with thickness of 300 μm, in which no polar phonon exists in our range of interest. Using the latter substrate, one can eliminate the hybridization effect of surface-polar phonons with plasmon, and hence measure the intrinsic plasmon dispersion. Disk/rectangle/ ribbon arrays were fabricated from the exfoliated films by electron-beam lithography and subsequent reactive ion etching. Sulfur hexafluoride (SF 6 ) was used as the reaction gas. For the rectangle (ribbon) arrays, the gaps between adjacent rectangles along a and b axes (ribbons) are kept larger than half of L a and L b (ribbon width), so that the rectangles (ribbons) can be approximately regarded as isolated resonance cavities. Because of the inevitable lateral etching, the dimensions listed in Fig. 3a are nominal values, with real length of about 10% uncertainty and the corresponding aspect ratio of 2 ± 0.15 (see Supplementary Table 1). The film thickness d film is normalized to 100 nm with the wave vector q multiplied by d film /100 in Figs. 4a and 5b, given that the sheet conductivity σ is proportional to the thickness.
Far-IR optical spectroscopy. For the polarized far-IR measurements, we used a Bruker FTIR spectrometer (Vertex 70 v) integrated with a Hyperion 2000 microscope and a liquid-helium-cooled silicon bolometer as the detector. The incident light was focused on WTe 2 films with a 15× IR objective. A THz polarizer was used to control the light polarization. The low-temperature measurements were performed in a helium-flow cryostat (Janis Research ST-300) with a pressure of about 1 × 10 −6 mbar. Throughout the entire measurements, nitrogen gas was purged to an enclosed space housing the cryostat. This procedure minimized the absorption of infrared light by moisture in air and effectively increased the signal/noise ratio.
Fitting of the extinction spectra of the bare film and disk/rectangle/ribbon arrays. The extinction spectra are determined by the complex dynamic conductivity σ(ω) where Z 0 is the vacuum impedance, ω is the frequency of light, and n s is the refractive index of the substrate. The dynamic conductivity contributed by the plasmon mode in a disk/rectangle/ribbon array is given by where ω P and Γ P are the frequency and resonance width of the plasmon, S P is the spectral weight, and f is the filling factor (WTe 2 microstructure area over the total area). The extinction spectra of the bare film in Fig. 1d are fitted using Eqs. (1) and (4). Spectral weights D and S, and the resonance widths Γ and η, are determined along both axes. The frequency of the bound states ω b is fixed to 800 cm -1 , which is between the frequencies of the two bound states in Fig. 1j.
Plasmonic extinction spectra in the disk, rectangle, and ribbon arrays are fitted using Eqs. (4) and (5), with the contribution of possible residual Drude response and bound states taken into account. For plasmons, spectral weight S P , plasmonic resonance frequency ω P , and resonance width Γ P are fitting parameters. The fitted plasmonic resonance width of the disk array is 30% smaller than the Drude scattering width in the bare film. This is probably due to the fact that the measured spectra range at the low-energy side is limited by light intensity and available only above 100 cm -1 . Thus, only a portion of the Drude spectrum is measured, compromising the fitting accuracy. This also explains why the singularity in Fig. 2c is not clear enough in the temperature dependence of Drude weight in Fig. 1h.
Calculating the loss function and isofrequency contours. The pseudocolor plot in Fig. 4a shows the fitting result of loss function based on the dielectric function of Eq. (3) and the conductivity of Eq. (1) described in the main text. Spectral weights D and S, and the resonance widths Γ and η, are fitted along both principal axes. To simplify the calculation, only the interband transition with the lowest resonance energy of 710 cm -1 is considered. The dielectric constant of environment ε env is set to 3.3. The solid black curves in Fig. 4a are plotted to represent the peak position of the maximum value of the loss function at each given wave vector. The fitted Drude weights are 8.08 × 10 11 and 4.49 × 10 11 Ω -1 ·s -1 along a and b axes, respectively. The mass ratio calculated from the Drude weight is about 1.8, consistent with the ones obtained from the bare film and disk array. The fitted interband transition weights are 4.31 × 10 11 and 8.07 × 10 11 Ω -1 ·s -1 along a and b axes, respectively. The weight ratio is about 1.87, with a larger amplitude along the b axis, which is consistent with the results in the unpatterned film. The fitted scattering width of the Drude modes and bound states is nearly identical along both axes, with values of about 60 and 170 cm -1 , respectively. For a ribbon array with a skew angle of θ with respect to the a axis, in which case the wave vector q is not along the principal axes, the loss functions are calculated by substituting the conductivity σ = σ aa cos 2 θ + σ bb sin 2 θ into Eq. (3). See Supplementary Fig. 9 for the loss function calculation results and the measured plasmon dispersions in skew ribbon arrays. It should be noted that, in the hyperbolic regime, the uniaxial surface supports the hybrid transversemagnetic (TM)-transverse-electric (TE) polaritons. However, in the nonretarded regime (q ≫ ω/c), which is the case for our work, the hybrid polaritons are dominated by the TM field components, which are the plasmon modes studied in this work 48 .
At each skew angle, a line is plotted to represent the peak position of the maximum value of the loss function at each given wave vector, as shown in Supplementary Fig. 9. Thus, at a given frequency and skew angle, the wave vector length in the isofrequency contour (lines of Fig. 5b) is determined by selecting the q value in the loss function-fitting lines at the corresponding skew angle and frequency.

Data availability
The data that support the findings of this study are available from the corresponding author upon request.