Hyperbolic Plasmons Propagate through a Nodal Metal

Metals are canonical plasmonic media at infrared and optical wavelengths allowing one to guide and manipulate light at sub-diffractional length scales. A special form of optical waveguiding is offered by highly anisotropic crystals revealing different signs of the dielectric function along orthogonal directions. These latter types of media are classified as hyperbolic and many crystalline insulators, semiconductors and artificial metal-based metamaterials belong to that class. Layered anisotropic metals are also anticipated to support hyperbolic waveguiding. Yet this behavior remains elusive primarily because interband processes introduce extreme losses and arrest light propagation. Here, we report on the observation of propagating hyperbolic waves in a prototypical layered nodal-line semimetal ZrSiSe. The unique electronic structure with touching energy bands at nodal points/lines suppresses losses and enables a hyperbolic regime at the telecommunications frequencies. The observed waveguiding in metallic ZrSiSe is a product of polaritonic hybridization between near-infrared light and long-lived nodal-line plasmons. By mapping the energy-momentum dispersion of the nodal-line hyperbolic modes in ZrSiSe we inquired into the role of additional screening associated with the surface states. Main text: The dielectric permittivity tensor (ε = ε# + iε&) in layered materials exhibits strong anisotropy such that the real part ε# may display opposite signs along different crystal axes. In this so-called hyperbolic regime, the isofrequency surface of the electromagnetic wave is an open hyperbola instead of the usual closed sphere or ellipse. In the hyperbolic regime, the interaction of light with collective modes of the matter results in hyperbolic polaritons exhibiting exotic properties including ray-like waveguiding inside the crystals at deeply subdiffractional wavelength. Propagating hyperbolic polaritons were previously explored in polar insulators, including: hBN, MoO3, V2O5 , calcite and in semiconducting WSe2. Layered anisotropic metals have been predicted to manifest hyperbolicity across a broad spectrum. The prerequisite anisotropy in a layered metal stems from suppressed interlayer hopping leading to a stark directional-dependence of the (screened) plasma frequency ω(, which marks the zero-crossing of the dielectric function. Layered metals, in principle, could reveal a vast frequency range of hyperbolicity between ω() < ω < ω(+, , where ω(+,and ω() are in-plane and out-of-plane (screened) plasma frequencies. Yet, insofar propagating hyperbolic plasmon polaritons (HPP) in layered metals have eluded faithful detection by direct imaging. Even though preconditions for hyperbolic plasmons are satisfied in a wide variety of anisotropic conductors, the significant electronic losses in common metals have precluded waveguiding. Here, we use scanning near-field optical microscopy (SNOM) to visualize both the principal and the higher-order HPP modes. We demonstrate that the nodal-line semimetal ZrSiSe hosts propagating hyperbolic plasmon polaritons across the entire telecommunication wavelength range (λ ≈ 1.3 – 1.7 μm, or ω ≈ 7700 – 5900 cm3#). The unique nodal-line band structure leads to a strong suppression of interband optical absorption, enabling near-infrared (IR) waveguiding in the interior of a layered metallic slab. To access hyperbolic modes in ZrSiSe we performed two types of near-field optical experiments (Fig. 1a). The first one involves placing thin crystals of ZrSiSe on patterned gold antennas, which serve as launchers of hyperbolic modes into the interior of the sample. The second approach utilizes the sharp sample edge to launch the HPPs and reveals characteristic higherorder hyperbolic modes. The two complementary experiments produce consistent results.


Main text:
The dielectric permittivity tensor ( = # + & ) in layered materials exhibits strong anisotropy such that the real part # may display opposite signs along different crystal axes. In this so-called hyperbolic regime, the isofrequency surface of the electromagnetic wave is an open hyperbola instead of the usual closed sphere or ellipse [1][2][3][4][5][6] . In the hyperbolic regime, the interaction of light with collective modes of the matter results in hyperbolic polaritons exhibiting exotic properties including ray-like waveguiding inside the crystals at deeply subdiffractional wavelength. Propagating hyperbolic polaritons were previously explored in polar insulators, including: hBN 1,7 , MoO3 8,9 , V2O5 10 , calcite 11 and in semiconducting WSe2 12 . Layered anisotropic metals have been predicted to manifest hyperbolicity across a broad spectrum 4,[13][14][15] . The prerequisite anisotropy in a layered metal stems from suppressed interlayer hopping leading to a stark directional-dependence of the (screened) plasma frequency ( , which marks the zero-crossing of the dielectric function. Layered metals, in principle, could reveal a vast frequency range of hyperbolicity between ( ) < < ( +, , where ( +, and ( ) are in-plane and out-of-plane (screened) plasma frequencies 15 . Yet, insofar propagating hyperbolic plasmon polaritons (HPP) in layered metals have eluded faithful detection by direct imaging. Even though preconditions for hyperbolic plasmons are satisfied in a wide variety of anisotropic conductors 4,[16][17][18][19] , the significant electronic losses in common metals have precluded waveguiding.
Here, we use scanning near-field optical microscopy (SNOM) to visualize both the principal and the higher-order HPP modes. We demonstrate that the nodal-line semimetal ZrSiSe hosts propagating hyperbolic plasmon polaritons across the entire telecommunication wavelength range (λ ≈ 1.3 -1.7 µm, or ≈ 7700 -5900 cm 3# ). The unique nodal-line band structure leads to a strong suppression of interband optical absorption, enabling near-infrared (IR) waveguiding in the interior of a layered metallic slab.
To access hyperbolic modes in ZrSiSe we performed two types of near-field optical experiments (Fig. 1a). The first one involves placing thin crystals of ZrSiSe on patterned gold antennas, which serve as launchers of hyperbolic modes into the interior of the sample 2,3,12 . The second approach utilizes the sharp sample edge to launch the HPPs and reveals characteristic higherorder hyperbolic modes 2,7 . The two complementary experiments produce consistent results. We first focus on experiments involving an Au disk launcher underneath the crystal. As illustrated in Fig. 1a, the HPPs propagate as conical rays and emerge on the surface of the sample as "hot rings" with enhanced nano-optical contrast surrounding the edge of the Au antennas. The propagation angle θ (with respect to surface normal) is controlled by the anisotropic permittivities of the sample [1][2][3]12 : where δ is the separation of the rings on the top surface and d is the sample thickness. We examined the topography (Fig. 1b) and the corresponding near-field optical amplitude ( 7 (ω), Fig. 1c) image of a thin ZrSiSe crystal partially covering the gold disk. At ω = 6667 cm 3# , a clear double-ring like pattern (red dashed lines) emerges along the Au antenna boundary. This double-ring pattern is distinct from a featureless response in the interior of the sample (away from the Au antenna) and also distinct from the internal resonances of the Au antenna at much larger length scales (Extended Data Fig.3). In Fig. 1d we plot the simulated z-component of the electric field for a 26 nm thick hyperbolic medium fully covering a 25 nm thick gold antenna (see Methods). A clear double-ring feature appears in the simulation, in agreement with the experiment (Fig. 1c). The separation of the hot rings on the sample surface is anticipated to systematically change with varying frequency of the illuminating beam in a fashion predicted by Eq. (1). The varying ring-separation is confirmed by the experimental near-field amplitude images taken as a function of frequency with selected images shown in Fig. 1e. The blue and red dashed lines mark the positions of the hot rings at ω = 7143 cm 3# and ω = 5556 cm 3# , with gradually changing ring-separations for frequencies in between. In contrast, the double-ring feature is completely absent at ω = 2222 cm 3# (outside of the hyperbolic range), showing instead a homogeneous near-field response (see also Extended Data Fig.2).
We now inquire into quantitative details of propagating HPPs in ZrSiSe. We averaged the radial line profiles within the sectors depicted in Fig Fig. 2a). With the ab-plane permittivity known from previous measurements on bulk crystals 20 , the experimental ring-separation δ together with the sample thickness d allows for the extraction of the c-axis permittivity of ZrSiSe from Eq. (1). These latter data are displayed in Fig. 2b along with the experimental ab-plane permittivity (squares). The hyperbolic regime in ZrSiSe extends between ≈ 2837 − 9091 cm 3# (see Supporting Information Sec. S1, S2). Finally, we observe yet another hallmark of the hyperbolic ray propagation which is the linear scaling of inter-peak separation δ with the increasing sample thickness (Fig. 2c). Broadband hyperbolic electrodynamics in layered nodal metal ZrSiSe is therefore firmly established. The natural edges of thin hyperbolic materials can also launch polaritons (Fig. 1a), albeit with weaker launching efficiencies 21 . In order to explore HPPs near the edges, we focused on the phase contrast images. The enhanced sensitivity of the phase contrast identified the weak higherorder HPP modes: yet another electrodynamics signature of a hyperbolic medium 2,7 . In Figures  3a-3c, we present near-field phase contrast images obtained for a triangular-shaped 20 nm thin ZrSiSe crystal on Si/SiO2 substrate for a few representative laser frequencies.
At ω = 8333 cm 3# (Fig. 3a), the phase contrast reveals a strong fringe localized very close to the right edge, which gradually extends into the interior of the sample as the laser frequency decreases. Multiple fringes of different periodicity are particularly apparent at ω = 6250 cm 3# (Fig. 3c), analogous to the principle and higher-order polariton modes on the edge of hBN 2,7 . In order to accurately extract the HPP wavelength, we utilized a previously developed electromagnetic solver 22 (See Methods) to simulate the phase contrast with the complex polariton momentum Here, λ Y is the polariton wavelength and γ accounts for damping of the polariton wave. The experimental and simulation results for hyperbolic modes probed at ω = 6250 cm 3# are displayed in Fig. 3d, where the blue dotted line is the phase contrast line profile along the line indicated in Fig. 3c. A weaker peak ( # ) is visible around 70 nm from the edge of the crystal (grey shaded area), followed by a longer wavelength oscillation ( c ). To better resolve the shorter wavelength oscillations ( # ), we inspected the derivative of the phase line profile (orange dotted line). Both weaker and stronger peaks in the derivative line profiles are well described by the simulation (black dashed line) involving two polariton modes. As we show below, the two modes correspond to the principal and higher-order HPPs in ZrSiSe. Importantly, the location of higher-order polariton modes with shorter wavelengths is uniquely determined by the location of the weaker peak in the phase derivative line profile.
To map the dispersion of HPP modes, we simulated the phase derivative profiles for multiple frequencies. In Fig. 3e, we show the experimental and simulated phase derivative traces for frequencies between ω = 8333 cm 3# , and 5000 cm 3# . The phase derivative line profile obtained at ω = 8333 cm 3# can be adequately reproduced with a single damped polariton of wavelength λ Yc . However, for laser frequencies smaller than ω = 6944 cm 3# an additional mode with much shorter wavelength λ Y# is needed to fully account for the experimental data. The polariton momentum can then be extracted from the simulated polariton wavelength using Re q Y = &Z [ \ , allowing a direct comparison with predictions based on full experimental dielectric tensors in Fig. 2b.  Fig. 3. The data points are superimposed on the calculated imaginary part of the reflection coefficient, Im(rp), using experimental dielectric functions (Fig. 2b). The grey dashed line represents the free-space light cone. Green and black dashed lines indicate analytical solutions for the divergence of Im(rp) for the higher-order HPP branches, with and without surface state metallicity, respectively. Red dashed line is a guide for the dispersion of the principal branch. b, The band structure of ZrSiSe: bulk (blue) and 5-layer slab (red). Compared to the bulk bands, additional interband transitions (green arrow) associated with the surface states contribute to the plasmonic response. Inset shows the nodal-line structure (red lines) in the Brillouin zone of ZrSiSe. c, Schematic band structure vs +, , ) away from the high-symmetry plane (grey-shaded plane in the top inset of b). In partially gapped nodal metals, interband optical transitions at energies above the van Hove energy (Δ h ) are suppressed. d, Schematic of the real optical conductivity (blue) and the joint density-of-states (gray) for bulk ZrSiSe. Optical losses quantified by the conductivity are likewise rapidly suppressed at energies above Δ h .
The extracted HPP momenta are organized in a frequency-momentum ( , ) dispersion plot in Figure 4. It is customary to identify HPPs via the divergences of the reflection coefficient ( ( , ) 1,23 . A colormap of Im(rp) provides an instructive way to visualize both the dispersion and the damping of the HPP modes. The colormap in Fig. 4a is calculated for a 20 nm thick crystal of ZrSiSe residing on SiO2/Si substrate using experimental dielectric functions in Fig. 2b. As expected, multiple dispersive branches develop in the hyperbolic frequency range, corresponding to the principal and higher-order modes that are propagating within the bulk. The calculated hyperbolic dispersions agree with the experimental momenta (colored circles and triangles), unequivocally corroborating the notion of hyperbolic plasmon polaritons in ZrSiSe.
We remark that the dispersion of the higher-order branch (triangles) in Fig.4 slightly deviates from the predicted maximum of Im(rp) (black dashed lines). These latter deviations are more apparent in the ( )-dispersion at > 6400 cm 3# (≈ 0.8 eV), marked by a green dot in Fig.  4a. We recall that ZrSiSe shows a significant degree of electronic correlations 20,24,25 revealed through the reduction of the Drude spectral weight. Such interaction effects are expected to renormalize the plasmon dispersions at finite momenta and to slower the group velocity 26,27 . However, the observed dispersion anomalies Fig.4a reveal an opposite trend, implying an enhanced group velocity.
Searching for the root cause of the dispersion anomaly in Fig.4a, we focus on the possible role of surface states. ZrSiX (X=Se, S, Te) family along with other topological semimetals are known to host various forms of conducting surface states [28][29][30][31] . In ZrSiSe, these surface states produce enhanced metallicity attested by band structure calculations for a thin slab (Fig.4b) and therefore could modify the polaritonic dispersion. Indeed, the HPP mode dispersion admits a modified dispersion prompted by the surface state metallicity (see Methods). By adding surface state layers with a (screened) plasma frequency (,k = 1 eV 30 and an interband transition term at c = 0.78 eV, the agreement between the data (triangles) and prediction (green dashed lines) is much improved (Fig. 4a).
We now comment on the preeminent role that the nodal-line structure of the ZrSiSe family of materials 31,33 plays in both broadband hyperbolicity and reduced interband loss. Dirac nodal-line semimetals possess band-touching points extended through lines or loops in the Brillouin zone [34][35][36] . In ZrSiSe, the unique cage-like nodal-lines (Fig. 4b inset) imply that the bands are essentially "flat" along the c-axis while remaining highly dispersive in the ab-plane 20,31,37 , thus giving rise to ( +, ≫ ( ) . More importantly, the combination of quasi-2D joint density-of-states (JDOS) and the abundance of van Hove singularities suppress the interband loss in ZrSiSe. Specifically, the optical conductivity of ZrSiSe follows # ( ) ∼ . across a broad range of frequencies 20,36,38,39 ; this latter behavior is a product of linear-in-frequency JDOS that resembles graphene. Van Hove singularities exist between the nodal-lines 20,37 and the saddle-point structure (Fig. 4c) leads to a decrease in JDOS above the van Hove energy (Δ h ). As a result, the corresponding # ( ) is suppressed above 0.4 eV (Fig. 4d, Fig. S1). The dip in the dissipative part of the optical conductivity # for ZrSiSe leads to the regime of high & / # ≈ 3 (Extended Data Fig.1) ultimately enabling propagation of only moderately damped plasmon polaritons in near-IR (Extended Data Fig.5). In contrast, the & / # is less than 1 for other hyperbolic materials candidates in the near-IR range 4,16 . The net effect is that nodal metals represent a unique type of hyperbolic conductors with losses low enough to support propagating HPPs at near-IR frequencies, as we have demonstrated by direct nano-imaging for ZrSiSe.
The low-loss hyperbolic plasmons in the near-IR range benefit from the combination of the high density of Dirac fermions and the unique nodal band structure in Fig.4. Remarkably, the suppression of interband optical transition offers a concrete example to the hypothesis of the "band structure engineering" approach to low-loss hyperbolic plasmonics 40,41 . The propagating hyperbolic plasmon polaritons reported here demonstrate yet untapped potential of nodal-line semimetals as low-loss natural hyperbolic media. We anticipate similar effects in the closely related ZrSiS 42 and other nodal semi-metals 43 . We envision future band structure engineering based on these nodal-line semimetals could lead to further increase of the polariton quality factor and propagation length in the technologically important near-IR wavelength. Beyond applications, the multiple hyperbolic branches reported here constitute a crucial first step in exploring the intertwined bulk and surface states effects in topological semimetals with polaritonic observables.

Methods
Single crystal growth and device fabrication. The ZrSiSe single crystals were synthesized using a chemical vapor transport method as described previously 20,33 . For Au antenna patterned devices, Au/Cr (25 nm/1 nm) disks were e-beam deposited on SiO2/Si substrates following standard e-beam lithography processes using a lift-off resist. ZrSiSe flakes were then directly exfoliated on Au/Cr disks in a glovebox filled with inert gas (O2 < 1 ppm, H2O < 0.1 ppm). Prior to exfoliation, the substrates were annealed in glovebox at 250 °C for 1 hr to remove any residual moisture on the surface.

Near-infrared nano-optical measurements.
We used a scattering-type scanning near-field optical microscope (s-SNOM, Neaspec) based on an atomic force microscope (AFM) operating in tapping mode. The tapping frequency of the AFM tip is around 70 kHz and near-field data are collected at higher harmonic (n = 3 or 4) of the tapping frequency to suppress the far-field background. For the gold antenna launcher experiment the difference frequency generation outputs of a pulsed laser source (Pharos, Light Conversion) were used. For the edge-launching experiments, we used a continuous-wave tunable laser from M Squared. Tunable outputs between 1140 nm -2200 nm are generated by frequency mixing of a high-power 532 nm diode laser (EQUINOX) and a Ti:Sapphire laser tunable between 700 nm -1000 nm (SolsTiS).
Geometrical correction considering finite antenna thickness. Due to the finite thickness of the underlying Au antenna, the sample exhibits a downward slope outside the antenna boundary ( Fig.1b). This small slope (tan ) leads to a finite asymmetry ( with the operation * denoting convolution, Φ( ) being the full potential, and ( ) being the in-plane Coulomb potential. The external field is taken to be constant with an incidence angle = 60 ∘ from the surface normal. Using the translation symmetry of the problem in the lateral direction, Eq. (S2) is reduced to a one-dimensional integral equation. Replacing the derivatives by finite differences further simplifies Eq. (S2) to a matrix inversion. The SNOM signal is then found by computing the dipole moment induced on the probe by the distribution ˆ( ). The complex signal (S) is calculated for a range of tip positions from each quantity ˆ( ) and then demodulated to the 4 th harmonic in order to compare with the experimental line profiles shown in and in the lossless limit (Re &™ = 0) we obtain: š ( ) = I ) I| +, | + 2 | +, | + 8 Im( &™ )/ ( 5)