
==== Front
Proc Natl Acad Sci U S A
Proc Natl Acad Sci U S A
PNAS
Proceedings of the National Academy of Sciences of the United States of America
0027-8424
1091-6490
National Academy of Sciences

38466854
202319465
10.1073/pnas.2319465121
research-articleResearch ArticlephysPhysics426
Physical Sciences
Physics
Spatiotemporal beating and vortices of van der Waals hyperbolic polaritons
Zhang Tianning a b 1
Yan Qizhi a b 1 https://orcid.org/0000-0001-5440-5888

Yang Xiaosheng a b 1 https://orcid.org/0000-0002-7632-0401

Ma Weiliang a b
Chen Runkun a b
Zhang Xin a b
Janzen Eli c
Edgar James H. c https://orcid.org/0000-0003-0918-5964

Qiu Cheng-Wei chengwei.qiu@nus.edu.sg
d 2 https://orcid.org/0000-0002-6605-500X

Zhang Xinliang xlzhang@mail.hust.edu.cn
a b e 2
Li Peining lipn@hust.edu.cn
a b 2
aWuhan National Laboratory for Optoelectronics and School of Optical and Electronic Information, Huazhong University of Science and Technology, Wuhan 430074, China
bOptics Valley Laboratory, Wuhan 430074, China
cTim Taylor Department of Chemical Engineering, Kansas State University, Manhattan, KS 66506
dDepartment of Electrical and Computer Engineering, National University of Singapore, Singapore 117583, Singapore
eOffice of the President, Xidian University, Xi’an 710126, China
2To whom correspondence may be addressed. Email: chengwei.qiu@nus.edu.sg, xlzhang@mail.hust.edu.cn, or lipn@hust.edu.cn.
Edited by David Weitz, Harvard University, Cambridge, MA; received November 7, 2023; accepted February 19, 2024

1T.Z., Q.Y., and X.Y. contributed equally to this work.

11 3 2024
19 3 2024
11 9 2024
121 12 e231946512107 11 2023
19 2 2024
Copyright © 2024 the Author(s). Published by PNAS.
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This article is distributed under Creative Commons Attribution-NonCommercial-NoDerivatives License 4.0 (CC BY-NC-ND).

Significance

Two-dimensional van der Waals (vdW) materials can squeeze light to subwavelength scale and support so-called hyperbolic polaritons (HPs), a peculiar light–matter interaction mode where the light is steered along conical rays. Here, we present a time-domain near-field interferometry approach to characterize the spatiotemporal beating of HPs as they undergo total internal reflections while propagating in the vdW material. These reflections produce spatiotemporal optical vortices occurring on nanometer and femtosecond scales along the trajectory. We further demonstrate the potential of utilizing highly confined HP rays to transfer spectroscopic information of molecules with a vdW material.

In conventional thin materials, the diffraction limit of light constrains the number of waveguide modes that can exist at a given frequency. However, layered van der Waals (vdW) materials, such as hexagonal boron nitride (hBN), can surpass this limitation due to their dielectric anisotropy, exhibiting positive permittivity along one optic axis and negativity along the other. This enables the propagation of hyperbolic rays within the material bulk and an unlimited number of subdiffractional modes characterized by hyperbolic dispersion. By employing time-domain near-field interferometry to analyze ultrafast hyperbolic ray pulses in thin hBN, we showed that their zigzag reflection trajectories bound within the hBN layer create an illusion of backward-moving and leaping behavior of pulse fringes. These rays result from the coherent beating of hyperbolic waveguide modes but could be mistakenly interpreted as negative group velocities and backward energy flow. Moreover, the zigzag reflections produce nanoscale (60 nm) and ultrafast (40 fs) spatiotemporal optical vortices along the trajectory, presenting opportunities to chiral spatiotemporal control of light–matter interactions. Supported by experimental evidence, our simulations highlight the potential of hyperbolic ray reflections for molecular vibrational absorption nanospectroscopy. The results pave the way for miniaturized, on-chip optical spectrometers, and ultrafast optical manipulation.

hyperbolic polariton
boron nitride
spatiotemporal optical vortex
==== Body
pmcPolaritons—hybrid light–matter excitations—in layered van der Waals (vdW) materials provide opportunities for nanoscale light manipulation (1–5). Specifically, hyperbolic polaritons (HPs) in these materials showcase distinct hyperbolic dispersions in the momentum space, facilitating unusual ray-like propagation accompanied by enhanced density of optical states (6–9). Within thin slabs, HP rays follow a zigzag trajectory, dictated by total internal reflection at boundaries. This gives rise to an unlimited number of subdiffractional waveguide modes at a single frequency for lossless cases (8, 10–12). The unique abilities of HP rays and their composite subdiffractional waveguide modes have potential applications in hyperlensing (8, 13, 14), high-sensitivity sensing (15, 16), negative refraction (17–19), twisted polaritonic metasurfaces (20–24), and photo-induced programmable waveguiding (25).

While parallels may be drawn between HP rays’ directional propagation behavior and the light beam in conventional attenuated total reflectance (ATR, see SI Appendix, Note S1 and Fig. S1) (26, 27), there are important differences that fundamentally shape their potential uses. Unlike ATR where different frequency components align in a single light beam, deep-subwavelength confined HP rays typically manifest strong frequency dispersion (28–30). This frequency dispersion is demonstrated in hexagonal boron nitride (hBN), an exemplary polar vdW material, see Fig. 1A. Extraordinary HP rays exist in the type-II Reststrahlen band (1,395 to 1,650 cm−1) of hBN with opposing in- and out-of-plane permittivities (εx = εy < 0, εz > 0) (7, 31). For these rays, as the operation frequency ω of the dipole source increases, the angle of propagation with respect to the optic axis (z axis) decreases (8, 9). This frequency dependence raises doubts about utilizing highly confined HP rays to transfer broadband information to a certain position.

Fig. 1. Transient nanoimaging of HP rays in a thin hBN slab. (A) Propagation of HP rays excited by a continuous-wave dipole source, where the angle of HP ray with respect to the optic axis (z axis) decreases with increasing frequency ω. (B) Propagation of HP rays excited by a pulsed laser. (C) Schematic of interferometric near-field nanoimaging of HP rays launched by a gold antenna on hBN. (D) Images of second order near-field amplitude |S2| with different time delay t around the upper end of the gold antenna. The area shaded with yellow color corresponds to the gold antenna, where the experimental signals are omitted for better illustration. Green dashed lines show the wavefronts of fringe patterns around the antenna.

Here, we explore this challenge by using a pulsed laser to generate HP rays in an hBN slab (depicted in Fig. 1B). The resulting HP ray pulses encompass a broad spectrum of frequency components (32). Their trajectory along the reflection pathway within the hBN slab shares similarities to the conventional ATR setup, but on a significantly smaller nanoscale. We show in the following that this approach holds potential for transmitting spectroscopic data through vdW systems using pulsed HP rays.

Results

To experimentally validate our proposal, we utilized the measurement setup depicted in Fig. 1C. The sample consisted of an 80-nm-thick hBN slab with a 2.3-μm-long gold rod functioning as a resonant infrared antenna (33) (see monochromatic nanoimaging experiments, cf. SI Appendix, Note S2 and Fig. S2). We sent infrared laser pulses, covering the type-II Reststrahlen band of hBN (Materials and Methods), into an asymmetric Michelson interferometer. Upon illuminating the sample, the resonant gold antenna efficiently triggered the excitation of HP ray pulses, which propagated along a zigzag trajectory within the hBN slab. These were detected by a scattering-type scanning near-field optical microscope (s-SNOM) equipped with a sharp metalized tip. We captured images of the near fields Ep of antenna-launched polariton pulses alongside the background illumination Ein with varying time delay t (adjusted by moving the reference mirror, see Materials and Methods). This interferometric technique allowed us to record time-resolved, complex-valued near-field optical images with unprecedented spatiotemporal resolutions. While recent advances in time-resolved transmission electron microscope have enabled a direct observation of polariton wavepackets (34), the beating between multiple dispersion branches has yet been thoroughly studied. Our approach provides deep insights into the physics of pulsed HP rays and dynamics arising from the coherent beating of different order hyperbolic waveguide modes.

The Fig. 1D displays near-field amplitude images measured at varying time delays, revealing distinct bright and dark circular fringes of HP pulses in close proximity to one extremity of the gold antenna. We precisely tracked the time-evolving positions of the fringes (e.g., dark fringe marked by the dashed green line), showing that the pulse fringes gradually shift inward (toward the antenna) at time intervals t = 320, 333, and 347 fs. However, at t = 360 fs, a surprising development occurs where the fringes jump abruptly outward (away from the antenna) compared to the earlier moment at t = 347 fs. Intriguingly, this alternating inward–outward movement repeats as the time delay is increased, challenging the common sense of energy flow as the polaritonic rays are emitted from the antenna. These behaviors stem potentially from the interference of polariton pulse fields (Ep) and background illumination (Ein), or from the intrinsic exotic dynamics of HP ray pulses. The latter explanation presents a more intriguing avenue for in-depth investigation.

To understand such anomalous fringe movements, we performed detailed line-scan imaging of the complete space-time evolution of the antenna-launched HP pulse fields, accounting for the in-plane isotropic (circular) pulse propagation. By applying Fourier transform data filtering processing (see SI Appendix, Note S3 and Figs. S3 and S4 for details), we obtained background-subtracted HP pulse fields, which are shown as the amplitude and real part of the signal in Fig. 2 A and B, respectively. This thus allowed us to determine the intrinsic propagation dynamics of HP ray pulses. There are two primary branches in the space-time data (Fig. 2 A and B), showing intensity maxima at the bottom-left corner of the image, near the extremity of the gold antenna (located at the origin). The first, broader branch showed a diagonal pattern and decays gradually toward the top-right corner, signifying that the pulse fringe moved away from the antenna over time. In contrast, the second, more tightly confined branch had a steeper slope (lower velocity) than the first but showed alternating oscillations of high and low intensity (empty circles or squares marking the local minima), as indicated by dashed gray lines in Fig. 2A. These oscillation patches have negative slopes (∂x/∂t < 0, negative velocities), confirming the existence of backward-moving fringes observed in Fig. 1D. The presence of multiple alternating oscillations with negative slopes also explains the apparent forward-jumping of fringe patterns when oscillations with longer time delays appear, which is consistent with the experimental results shown in Fig. 1D.

Fig. 2. Capturing the spatiotemporal beating of HPs in hBN. (A and B) Amplitude (A) and real part (B) of the filtered experimental near-field data S(x, t) in the space-time domain. Dashed lines of negative slope indicate backward-moving fringes. Empty circles and squares correspond to the local amplitude minima. (C and D) Line profiles (blue curves) extracted at fixed spatial positions (C) or time delays (D). Corresponding envelopes of the fringe patterns are shown as black curves. ve, envelope velocity; vf, fringe velocity. (E) Separated M0 and M1 modes in the space-time domain through spatial filtering. (F) Simulated electric field distributions in air and hBN in the xz plane at two different time delays. Corresponding local minima of amplitude in air are marked as in (A and B). Purple dashed lines indicate the extent of the zigzag propagating HP ray in hBN. Insets are sketches showing a train of HP ray pulse, whose leading or trailing edge undergoes the reflection on the interface, producing local minima of amplitude above hBN. (G) Short-time Fourier transform signals as a function of time delay, extracted at the position shown as vertical dashed line in (B). (H) Line profiles extracted from the spectrogram in (G) at different time delays. a.u., arbitrary units.

To further support our observations, we extracted line profiles from Fig. 2B at fixed values of x or t, respectively. The curves are shown in Fig. 2 C and D, together with their envelopes (black curves). By comparing the local maxima of the envelopes at different x values (Fig. 2C, x = 0.56 μm and 0.71 μm), we detected the presence of both positive and negative envelope velocities (ve, black arrows) associated with forward and backward propagating fringes, respectively. Conversely, line profiles at three different t values (Fig. 2D, t = 504, 514, and 524 fs) all showed positive fringe velocities (vf, blue arrows), which are in line with previous reports of positive phase velocities for type-II HPs in hBN (32). However, the presence of both positive and negative envelope velocities ve for type-II HPs suggests the co-existence of positive and negative group velocities vg, which could imply the simultaneous forward and backward energy flow. This finding contradicts the conventional understanding of HPs propagating in vdW materials.

To gain more insights into this unexpected phenomenon, we isolated the two branches of pulse fringes by applying a filtering approach on the line-scan imaging data, as described in SI Appendix, Note S3. Analysis of the filtered data in the space-time domain (Fig. 2E) revealed the two branches of pulse fringes correspond to two distinct polariton modes. As proved by the filtered data shown in the momentum-frequency domain (SI Appendix, Fig. S4), the dispersions of these two modes in our experiments are in excellent agreement with the theoretical dispersion of the first two orders of waveguide modes supported in the hBN slab, M0 and M1. Therefore, the two modes are assigned to the M0 and M1 modes, respectively.

After isolating the M0 and M1 modes, we observed both modes moving away from the gold antenna with time (Fig. 2E), which is consistent with causality. However, their superposition during propagation, as shown in Fig. 2 A and B, indicates multiple alternating oscillations with negative slopes. This can be attributed to the phenomenon that when two waves with slightly different frequencies interfere, the so-called beating effect between the two waves produces constructive and destructive interference and therefore perceived oscillation in intensity. In fact, the HP rays in the slab are essentially the coherent beating of different waveguide modes (8, 35). Here, we show that the beating between the M0 and M1 modes can effectively capture the majority of the characteristics of HP rays. This gains further validation through experimental and simulated near-field distribution, as illustrated in SI Appendix, Fig. S2. The beating patterns of the M0 and M1 modes can indeed construct the periodic zigzag trace of HP rays. Building on the understanding of the beating of M0 and M1 modes already effectively constructing the HP rays, we estimated the decay length (~0.5 μm), lifetime (~0.4 ps) and ultraslow group velocity (~4 × 10−3 c) of HP rays via quantifying the overlaps of M0 and M1 modes in the space-time domain (SI Appendix, Note S4 and Fig. S5).

To comprehend the measured motion of HP ray pulses, analytical calculations of the time-dependent near-field distribution both inside and above the hBN slab are needed (SI Appendix, Note S5 and Fig. S6). In this analysis, we take into account the contributions of the M0 and M1 modes in constructing the HP rays. The two-dimensional simulation was performed for an air/hBN interface in the (x, z) plane. Fig. 2F shows the simulated field distributions at two different time delays, which correspond to two local minima, marked as empty circle (at 530 fs) and square (at 640 fs), respectively, of the experimental near-field amplitude in Fig. 2A. Since our near-field nanoimaging experiment detects signal primarily from the region close to the hBN surface, these two local minima are well reproduced in the simulated field distribution above hBN. Also, when increasing the time delay from 530 fs to 640 fs, the local minimum in simulated field distribution moves backward in the −x direction, in good agreement with the observation of backward-moving fringe patterns in experiment. Inside hBN, however, the simulated electric field maximum shifts in the +x direction with increasing time, which implies that the HP ray pulse moves forward in the hBN slab. Here, we have used purple dashed lines as guides to the eyes to indicate the broad zigzag path of HP ray pulses. These findings demonstrated that an apparent backward movement of fringe patterns and the forward propagation of HP ray pulse indeed occur simultaneously. Interestingly, we found that the local minima in the near-field amplitude (Fig. 2A) appear in pairs and they locate at the moments when the leading or trailing edge of HP ray pulse gets reflected on the interface, as schematically shown in the Inset of Fig. 2F. Later, we will show that these local minima marked as empty circles or squares are related to opposite phases in the space-time domain.

According to the analysis above, we have determined the spatial positions where the HP ray pulse undergoes a total internal reflection, such as the near-field local minimum marked with empty circle in Fig. 2F. This prompts us to explore the frequency contents carried by HP ray pulses at these positions for potential spectroscopic investigation. We turn to the first reflection point (x = 0.45 μm, vertical dashed line in Fig. 2B) where the first local minimum is located in the space-time data. We performed short-time Fourier transform (STFT, see Materials and Methods) over time delay from 0 to 2 ps for the signal at this spatial position. The resulted spectrogram is shown in Fig. 2G and three line-scans extracted from the spectrogram are shown in Fig. 2H. In the line-scans at t = 300 fs and 740 fs, we can see broadband frequency components ranging from 1,400 cm−1 to 1,600 cm−1. Their magnitudes vary with the time delay and the frequency content also exhibits a redshift, as indicated by the change of peak from ~1,500 cm−1 to ~1,460 cm−1, which is possibly due to the faster decay of higher frequency components. In the line-scan at t = 560 fs, however, the frequency components are strikingly different, showing a clear minimum at approximately 1,480 cm−1. The varying frequency contents of HP rays, together with the time dependence, may bring about spectroscopic studies with selected spatiotemporal and spectral resolutions.

Inspecting the spatiotemporal characteristics of the reflection points in close-up, we uncovered another intriguing feature: the presence of spatiotemporal optical vortices (STOVs), a unique variety of optical vortices in the spatiotemporal domain that has garnered significant attention lately (36–40). We extracted the amplitude, the real part and the corresponding phase around two local minima (empty circle at 530 fs and empty square at 640 fs) from Fig. 2 A and B and plotted the enlarged views in Fig. 3. With both minima in amplitude (Fig. 3 A and D), we found different fork-shaped dislocations in the real-part distributions (Fig. 3 B and E) in a relatively small spatial (<60 nm) and temporal (<40 fs) window, according to the full width at half maximum in the extracted line profiles at these positions. Since our time-domain interferometric measurements simultaneously acquire near-field amplitude and phase, we also show the corresponding phase images in Fig. 3 C and F, where spiral phase of opposite directions can be identified. The presence of null intensity surrounded by phase circulation is one of the major features of optical vortices, which are commonly detected in the real space (41–43). More specifically, as the data are in the spatiotemporal domain, the measured optical vortices are STOVs. Orbital angular momentum (OAM) is also conserved, as is characteristic of STOVs (44–47). As shown in Figs. 2 and 3 (red arrows showing clockwise and anticlockwise phase circulations), the STOVs of topological charge l = −1 (empty square) and l = +1 (empty circle) always appear in pairs in experiment. The experimental generation of STOVs with transverse OAM defined on the spatiotemporal (x, t) plane is considered more complex than the generation of conventional optical vortices with longitudinal OAM on the (x, y) plane. Here, we have demonstrated the generation of STOVs as a result of the beating between propagating HP modes of different phase velocities. These occur on the much smaller nanometer and femtosecond scales in a vdW material, compared to the common implementation in free space where STOVs locate on millimeter and picosecond scales (39, 40). Note that our simulations reproduce our experimental findings on STOVs (SI Appendix, Note S6 and Fig. S7), verifying the generic character of STOVs present along the HP ray pulse’s trajectory.

Fig. 3. Spatiotemporal optical vortices of HP rays measured on the reflection point. Filtered experimental space-time near-field data at local amplitude minima (empty circles and squares, the same as in Fig. 2) with spatiotemporal optical vortices of topological charge l = −1 (A–C) and l = +1 (D–F). Columns from left to right are amplitude S(x, t) (A and D), real part Re[S(x, t)] (B and E) and phase ϕ(x, t) (C and F) signals. Color bars are rescaled for better demonstration. Plots with black and gray lines next to (A and D) are intensity profiles extracted from the amplitude data along the horizontal black line and the vertical gray line, respectively.

The next question that arises is whether such nanoconfined pulses can probe and carry spectroscopic information of the environment when undergoing reflection. We performed time-domain simulations (Materials and Methods and SI Appendix, Note S7) to investigate nanoscale transient ATR spectroscopy of molecular vibrational absorption. We consider molecules with a C–H vibration mode at 1,507 cm−1, indicated by the peak in the imaginary part dielectric permittivity of CBP molecule (16) (Fig. 4A). As illustrated in Fig. 4B, in the simulations, we place the molecules at the reflection position underneath an 80-nm-thick hBN slab and use a dipole source on the top surface to excite pulses of HP rays (complete rays, not only M0 and M1 modes, see SI Appendix, Fig. S8). The excited rays propagate to the bottom surface and interact with molecules via evanescent field coupling, before reflecting back to the top surface for detection. Similar to Fig. 2G, we obtain at the first reflection point STFT signals for cases with and without molecules (Fig. 4C). With molecules, an amplitude dip appears in the frequency of 1,507 cm−1 and at the time t = 3.0 ps. The dip at 1,507 cm−1 implies the coupling between the HP ray fields and the molecular vibrations. The duration for 3.0 ps agrees well with the time HP ray needed to travel from the origin to the first reflection point (SI Appendix, Figs. S5 and S8). At the same time, the detection sensitivity is dependent on the lateral position of molecules (SI Appendix, Fig. S8E). The transient spectrum (Fig. 4D, blue) shows a large frequency split (Δω) and dip depth (Δh, defined in Fig. 4D), with both Δω and Δh peaking at t = 3.0 ps (Fig. 4 E and F), indicating great visibility (and thus sensitivity) in the ω-t domain. The results suggest the potential of using HP rays for high-sensitivity transient nano-ATR spectroscopy.

Fig. 4. Simulations for demonstrating transient nanospectroscopy with HP rays. (A) Imaginary part of the dielectric permittivity εmol. of simulated molecules. (B) Schematic of realizing remote sensing with HP rays for molecules under a vdW slab. (C) Short-time Fourier transform signals as a function of time delay with and without molecules, detected at the position where the first internal reflection takes place as shown in (A). On the right side are line profiles extracted at 1,507 cm−1, where the spectrum with molecules present shows a dip at t = 3.0 ps. (D) Transient infrared spectrum at t = 3.0 ps (blue curve), in comparison with the spectrum when no molecule is included in simulation (gray curve). Δω, frequency range of the peak splitting; Δh, depth of the dip in spectra. (E and F) Change of the peak splitting Δω (E) and the dip’s depth Δh (F) as a function of time. a.u., arbitrary units.

Discussion

In summary, we have experimentally explored the rich spatiotemporal dynamics of HP rays in thin vdW materials, including beating, negative and positive velocities of fringe motions, and vortex phenomena. Compared to conventional surface polaritons, hyperbolic polaritons have the advantage of sub-diffractional localization and frequency-dependent directionality within the hyperbolic medium. This enables more potential applications where a precise nanoscale guiding of light is of interest. This work can be extended to other in-plane anisotropic vdW materials (48–50) (e.g., α-MoO3, V2O5, and others) and artificial nanostructures, where more fundamental physics and applications could be uncovered. Experiments and numerical simulations verified the transient spectral analyzing capabilities of HP rays for nanoscale ATR spectroscopy. This could expand conventional nanospectroscopy into the time domain, enabling ultrafast detection of dynamic molecular processes. Furthermore, the STOVs in a narrow spatiotemporal parameter space (on nanometer-femtosecond scales) could open opportunities in ultrafast optical manipulation, chirality probing, and quantum optics by harnessing polaritonic angular momentum on the deep-subwavelength scale.

Materials and Methods

Materials and Fabrication.

For experiments, we used isotopically (10B) enriched hBN following the growth recipe described in refs. 51 and 52. After the mechanical exfoliation, the hBN thin films were dry-transferred onto a silica substrate. Gold rod antennas were fabricated on the hBN surface via standard electron beam lithography. The antenna patterns were written on the spin-coated photoresist (PMMA: 495/A4, thickness ~100 nm). Standard lift-off procedure was conducted after depositing Ti (3 nm)/Au (40 nm) onto the developed photoresist.

Time-Domain Interferometric Measurements.

We employed a commercial nano-FTIR system (Neaspec GmbH, Germany) based on an s-SNOM comprising an asymmetric interferometer as shown in Fig. 1C. A broadband laser pulse of 100 fs duration ranging from 1,000 to 2,000 cm−1 illuminates the atomic force microscope (AFM) tip (Arrow NCPt, NanoWorld, oscillation amplitude 74 nm, frequency 268 kHz) and the tip-scattered light is recombined with the reference beam at the detector. The two-dimensional nanoimaging maps as a function of time delay t are measured by fixing the reference mirror at different positions. At each reference mirror position, the interferometric detector signal was demodulated at the second-order harmonic, yielding near-field amplitude S2 and phase ϕ2 images.

Numerical Simulations.

For simulating the near-field distributions, we used two-dimensional finite-difference time-domain software (Lumerical FDTD) using the same cross-section structure in the (x, z) plane as schematically shown in Fig. 4B. We place molecules (height 50 nm, length 0.1 μm) whose absorption peak is at 1,507 cm−1 below the hBN layer. The sample is illuminated by an electric dipole located 55 nm above the hBN top surface. The spectrum of the electric dipole is Gaussian shaped and central at 1,507 cm−1. The simulation range in the (x, z) plane is 55 μm × 35 μm with perfectly matched layers set as boundary conditions. Note that the simulation cell is large enough to prevent interference that could arise from the boundary, and no AFM tip is included in simulations. We monitor the time-dependent near-field distributions (Ez) and perform STFT to obtain the data in frequency–time (ω, t) domain, see SI Appendix, Note S7 for details. The dielectric permittivity tensor of hBN is modeled using the parameters from ref. 12. The parameters for the dielectric permittivity of SiO2 are taken from ref. 53.

STFT Calculations.

We used the Matlab function stft to perform STFT calculations on the signal along the time axis, obtaining the data in the frequency-time (ω, t) domain. For Fig. 2G, The signal within 2 ps was linearly interpolated to 10,000 points and broken down into 2,000 sample segments, with 1,900 samples of overlap between adjoining segments and a 20,000-point FFT (fast Fourier transform). For Fig. 4C, the signal within 6 ps was linearly interpolated to 10,000 points and broken down into 2,500 sample segments, with 2,495 samples of overlap between adjoining segments and a 90,000-point FFT. For both cases, Hamming windows were applied to each segment. The appropriate number of sample segments were chosen such that the obtained spectral width after Fourier transform is equivalent to the original spectral width.

Supplementary Material

Appendix 01 (PDF)

We acknowledge the support from National Natural Science Foundation of China (62075070), National Key Research and Development Program of China (2021YFA1201500), Hubei Provincial Natural Science Foundation of China (2022CFA053), and Innovation Fund of WNLO. C.-W.Q. is supported by the Competitive Research Program Award (NRF-CRP22-2019-0006) from the NRF, Prime Minister’s Office, Singapore, and by a grant (A-0005947-16-00) from Advanced Research and Technology Innovation Centre (ARTIC), National University of Singapore. The support for the hexagonal boron nitride crystal growth was provided by the Office of Naval Research, Award No. N00014-20-1-2474. We thank the Optoelectronic Micro&Nano Fabrication and Characterizing Facility, Wuhan National Laboratory for Optoelectronics of Huazhong University of Science and Technology for the support in sample fabrication.

Author contributions

C.-W.Q., Xinliang Zhang, and P.L. designed research; T.Z., W.M., Xin Zhang, E.J., and J.H.E. performed research; Q.Y. and R.C. contributed new reagents/analytic tools; T.Z., Q.Y., X.Y., Xin Zhang, and P.L. analyzed data; and T.Z., X.Y., C.-W.Q., and P.L. wrote the paper.

Competing interests

The authors declare no competing interest.

Data, Materials, and Software Availability

All study data are included in the article and/or SI Appendix.

Supporting Information

This article is a PNAS Direct Submission.
==== Refs
1 D. N. Basov, M. M. Fogler, F. J. García de Abajo, Polaritons in van der Waals materials. Science 354 , aag1992 (2016).27738142
2 W. Ma , Anisotropic polaritons in van der Waals materials. InfoMat 2 , 777–790 (2020).
3 Q. Zhang , Interface nano-optics with van der Waals polaritons. Nature 597 , 187–195 (2021).34497390
4 Y. Wu , Manipulating polaritons at the extreme scale in van der Waals materials. Nat. Rev. Phys. 4 , 578–594 (2022).
5 X. Guo , Polaritons in Van der Waals Heterostructures. Adv. Mater. 35 , 2201856 (2023).
6 D. N. Basov, A. Asenjo-Garcia, P. J. Schuck, X. Zhu, A. Rubio, Polariton panorama. Nanophotonics 10 , 549–577 (2020).
7 J. D. Caldwell , Sub-diffractional volume-confined polaritons in the natural hyperbolic material hexagonal boron nitride. Nat. Commun. 5 , 5221 (2014).25323633
8 S. Dai , Subdiffractional focusing and guiding of polaritonic rays in a natural hyperbolic material. Nat. Commun. 6 , 6963 (2015).25902364
9 P. Li , Hyperbolic phonon-polaritons in boron nitride for near-field optical imaging and focusing. Nat. Commun. 6 , 7507 (2015).26112474
10 A. Poddubny, I. Iorsh, P. Belov, Y. Kivshar, Hyperbolic metamaterials. Nat. Photonics 7 , 948–957 (2013).
11 S. Dai , Tunable phonon polaritons in atomically thin van der Waals crystals of boron nitride. Science 343 , 1125–1129 (2014).24604197
12 A. J. Giles , Ultralow-loss polaritons in isotopically pure boron nitride. Nat. Mater. 17 , 134–139 (2018).29251721
13 Z. Liu, H. Lee, Y. Xiong, C. Sun, X. Zhang, Far-field optical hyperlens magnifying sub-diffraction-limited objects. Science 315 , 1686–1686 (2007).17379801
14 P. Li , Collective near-field coupling and nonlocal phenomena in infrared-phononic metasurfaces for nano-light canalization. Nat. Commun. 11 , 3663 (2020).32694591
15 M. Autore , Boron nitride nanoresonators for phonon-enhanced molecular vibrational spectroscopy at the strong coupling limit. Light Sci. Appl. 7 , 17172 (2018).30839544
16 A. Bylinkin , Real-space observation of vibrational strong coupling between propagating phonon polaritons and organic molecules. Nat. Photonics 15 , 197–202 (2021).
17 X. Lin , All-angle negative refraction of highly squeezed plasmon and phonon polaritons in graphene–boron nitride heterostructures. Proc. Natl. Acad. Sci. U.S.A. 114 , 6717–6721 (2017).28611222
18 H. Hu , Gate-tunable negative refraction of mid-infrared polaritons. Science 379 , 558–561 (2023).36758071
19 A. J. Sternbach , Negative refraction in hyperbolic hetero-bicrystals. Science 379 , 555–557 (2023).36758086
20 G. Hu , Topological polaritons and photonic magic angles in twisted α-MoO3 bilayers. Nature 582 , 209–213 (2020).32528096
21 G. Hu, M. Wang, Y. Mazor, C.-W. Qiu, A. Alù, Tailoring light with layered and Moiré metasurfaces. Trends Chem. 3 , 342–358 (2021).
22 Z. Zheng , Phonon polaritons in twisted double-layers of hyperbolic van der Waals crystals. Nano Lett. 20 , 5301–5308 (2020).32574060
23 J. Duan , Twisted nano-optics: Manipulating light at the nanoscale with twisted phonon polaritonic slabs. Nano Lett. 20 , 5323–5329 (2020).32530634
24 M. Chen , Configurable phonon polaritons in twisted α-MoO3. Nat. Mater. 19 , 1307–1311 (2020).32661384
25 A. J. Sternbach , Programmable hyperbolic polaritons in van der Waals semiconductors. Science 371 , 617–620 (2021).33542134
26 P. Larkin, Infrared and Raman Spectroscopy (Elsevier, 2011).
27 R. Lu , High-sensitivity infrared attenuated total reflectance sensors for in situ multicomponent detection of volatile organic compounds in water. Nat. Protoc. 11 , 377–386 (2016).26820794
28 G. Hu, J. Shen, C.-W. Qiu, A. Alù, S. Dai, Phonon polaritons and hyperbolic response in van der Waals materials. Adv. Opt. Mater. 8 , 1901393 (2020).
29 G. Hu , Real-space nanoimaging of hyperbolic shear polaritons in a monoclinic crystal. Nat. Nanotechnol. 18 , 64–70 (2023).36509927
30 X. Ni , Observation of directional leaky polaritons at anisotropic crystal interfaces. Nat. Commun. 14 , 2845 (2023).37202412
31 A. J. Giles , Imaging of anomalous internal reflections of hyperbolic phonon-polaritons in hexagonal boron nitride. Nano Lett. 16 , 3858–3865 (2016).27159255
32 E. Yoxall , Direct observation of ultraslow hyperbolic polariton propagation with negative phase velocity. Nat. Photonics 9 , 674–678 (2015).
33 P. Pons-Valencia , Launching of hyperbolic phonon-polaritons in h-BN slabs by resonant metal plasmonic antennas. Nat. Commun. 10 , 3242 (2019).31324759
34 Y. Kurman , Spatiotemporal imaging of 2D polariton wave packet dynamics using free electrons. Science 372 , 1181–1186 (2021).34112689
35 S. Dai , Efficiency of launching highly confined polaritons by infrared light incident on a hyperbolic material. Nano Lett. 17 , 5285–5290 (2017).28805397
36 A. Chong, C. Wan, J. Chen, Q. Zhan, Generation of spatiotemporal optical vortices with controllable transverse orbital angular momentum. Nat. Photonics 14 , 350–354 (2020).
37 J. Ni , Multidimensional phase singularities in nanophotonics. Science 374 , eabj0039 (2021).34672745
38 C. Wan, Q. Cao, J. Chen, A. Chong, Q. Zhan, Toroidal vortices of light. Nat. Photonics 16 , 519–522 (2022).
39 C. Wan, A. Chong, Q. Zhan, Optical spatiotemporal vortices. eLight 3 , 11 (2023).
40 W. Chen, Y. Liu, Y.-Q. Lu, Spatiotemporal optical vortices: Toward tailoring orbital angular momentum of light in full space-time. ACS Photonics 10 , 2011–2019 (2023).
41 G. Spektor , Revealing the subfemtosecond dynamics of orbital angular momentum in nanoplasmonic vortices. Science 355 , 1187–1191 (2017).28302854
42 W.-Y. Tsai , Twisted surface plasmons with spin-controlled gold surfaces. Adv. Opt. Mater. 7 , 1801060 (2019).
43 M. Wang , Spin-orbit-locked hyperbolic polariton vortices carrying reconfigurable topological charges. eLight 2 , 12 (2022).
44 S. W. Hancock, S. Zahedpour, A. Goffin, H. M. Milchberg, Free-space propagation of spatiotemporal optical vortices. Optica 6 , 1547–1553 (2019).
45 S. Huang, P. Wang, X. Shen, J. Liu, Properties of the generation and propagation of spatiotemporal optical vortices. Opt. Express 29 , 26995–27003 (2021).34615122
46 G. Gui, N. J. Brooks, H. C. Kapteyn, M. M. Murnane, C.-T. Liao, Second-harmonic generation and the conservation of spatiotemporal orbital angular momentum of light. Nat. Photonics 15 , 608–613 (2021).
47 S. W. Hancock, S. Zahedpour, H. M. Milchberg, Second-harmonic generation of spatiotemporal optical vortices and conservation of orbital angular momentum. Optica 8 , 594–597 (2021).
48 W. Ma , In-plane anisotropic and ultra-low-loss polaritons in a natural van der Waals crystal. Nature 562 , 557–562 (2018).30356185
49 J. Taboada-Gutiérrez , Broad spectral tuning of ultra-low-loss polaritons in a van der Waals crystal by intercalation. Nat. Mater. 19 , 964–968 (2020).32284598
50 N. C. Passler , Hyperbolic shear polaritons in low-symmetry crystals. Nature 602 , 595–600 (2022).35197618
51 T. B. Hoffman, B. Clubine, Y. Zhang, K. Snow, J. H. Edgar, Optimization of Ni–Cr flux growth for hexagonal boron nitride single crystals. J. Cryst. Growth 393 , 114–118 (2014).
52 S. Liu , Single crystal growth of millimeter-sized monoisotopic hexagonal boron nitride. Chem. Mater. 30 , 6222–6225 (2018).
53 E. D. Palik, G. Ghosh, Handbook of Optical Constants of Solids (Academic Press, 1998).
