Ultrafast anisotropic dynamics of hyperbolic nanolight pulse propagation

Polariton pulses—transient light-matter hybrid excitations—traveling through anisotropic media can lead to unusual optical phenomena in space and time. However, studying these pulses presents challenges with their anisotropic, ultrafast, and nanoscale field variations. Here, we demonstrate the creation, observation, and control of polariton pulses, with in-plane hyperbolic dispersion, on anisotropic crystal surfaces by using a time-resolved nanoimaging technique and our developed high-dimensional data processing. We capture and analyze movies of distinctive pulse spatiotemporal dynamics, including curved ultraslow energy flow trajectories, anisotropic dissipation, and dynamical misalignment between phase and group velocities. Our approach enables analysis of polariton pulses in the wave vector time domain, demonstrating a time-domain polaritonic topological transition from lenticular to hyperbolic dispersion contours and the ability to study the polariton-induced time-varying optical forces. Our findings promise to facilitate the study of diverse space-time phenomena at extreme scales and drive advances in ultrafast nanoimaging.

Here, we demonstrate the experimental visualization of the inplane anisotropic propagation of HP pulses on calcite surfaces, using a phase-sensitive time-resolved infrared nanoimaging technique.We select calcite crystals for their ability to exhibit highly anisotropic HP propagation patterns at sizes ranging up to tens of micrometers (12).Scattering-type scanning near-field optical microscopy (s-SNOM) provides high spatial resolution reaching down to ~10 nm in broad spectral ranges from visible to terahertz frequencies (32,33), benefiting the material characterization and nanophotonic structure diagnosis (34).Our technique combines time-domain interferometry and s-SNOM, providing unparalleled imaging resolution for measuring the complex-valued HP pulse fields.To address the challenges associated with the high-dimensional space-time [i.e., two-dimensional (2D) space and time] imaging nature of our approach, including time cost, instability, and background noise, we develop a comprehensive strategy encompassing data acquisition, processing, and interpretation.Our nanoimaging approach and high-dimensional data processing enable the simultaneous assessment of anisotropic velocities v p and v g and dynamic visualization of their vectorial misalignment, providing unprecedented insights into the anisotropic spatiotemporal dynamics of HP pulses.

RESULTS
The experimental setup used in this study is illustrated in Fig. 1A.Mid-infrared pulses of 100-fs duration are directed into an asymmetric Michelson interferometer, with the pulse spectral range covering the Reststrahlen band of calcite (see Materials and Methods).The object beam illuminates a gold disk that is fabricated on the (100) surface of a calcite crystal, acting as an optical antenna and converting free-space light pulses to highly anisotropic HP pulses that propagate along the surface.A metallic SNOM tip is then brought close to the sample surface to scatter the excited HP pulse fields, which are superimposed with background illumination.The tip-scattered fields are cross-correlated with reference pulses and recorded by an infrared detector, while the interferometer time delay τ is adjusted by moving the reference mirror.To address instability caused by collecting the complete interferogram for each tip position on the sample surface, we develop a more robust scanning strategy where the tip is scanned over the surface at each fixed time delay τ, resulting in a 2D near-field snapshot of the HP pulse fields with nanoscale spatial resolution (see Materials and Methods).These snapshots are repeated when varying the time delay τ, creating a nano-movie-a high-dimensional space-time imaging dataset-that displays the anisotropic in-plane propagation of HP pulses, including their electric-field variations, anisotropic energy flow, and dissipation.
Figure 1B showcases experimental near-field amplitude snapshots that vividly capture bright and dark butterfly-like fringes of HP pulses emitted from the gold disk.These raw transient data (fig.S1) illustrate the dynamic evolution of the fringes in space and time, from the initial emergence through the build-up to the eventual annihilation.To gain a complete understanding of the evolution, a video comprising 317 frames, reconstructed via data interpolation (note S1), is available in movie S1.At short time delays (τ < 0.4 ps), a single fringe emerges near the disk, and with τ increasing, HP pulses propagate away from the disk, creating multiple fringes of larger concave curvatures.The fringes disappear at longer time delays (τ > 1.3 ps), indicating their relatively long lifetime.The validity and high quality of our spatiotemporal experiments are further confirmed by time-domain simulations that reproduce all observed HP pulse features (see fig.S2).
To accurately quantify the propagation of HP pulses, we use high-dimensional Fourier transform (FT) filtering to subtract the background fields and the field components of frequencies outside the Reststrahlen band (fig.S5 and detailed descriptions in note S2).The resulting filtered near-field data (real part) exhibit HP pulse fringes with intrinsic phase variations, which we captured in a sequence of snapshots shown in movie S2 and in Fig. 2 (A and B).We determine the fringe velocity v f by analyzing the gradient variations in the time-domain image sequence (see Materials and Methods), with v f parallel to the phase velocity v p and the momentum k.In Fig. 2 (A and B), we plot the time-dependent local distributions of v f as arrows (see zoom-ins in Fig. 2, C and D), each aligned with its corresponding v f .At short time delays (Fig. 2C), we observe diverse orientations of v f due to the appearance of small-k HPs.However, as time progressed, the arrows point uniformly in the same direction (Fig. 2D), indicating the appearance and dominance of highly directional large-k HPs in pulse-field distribution.
Figure 2E provides background-subtracted near-field amplitude snapshots (movie S3) that illustrate the anisotropic energy flow of HP pulses.This image sequence features butterfly-shaped polariton beams with four nearly symmetric branches, providing a clear visualization of how HP energies are emitted from the disk (τ = 0 ps), spread out along the sample surface (τ ≥ 0.2 ps), and dissipate over time.To analyze these dynamic processes, we track the centroid trajectory of each beam branch in the high-dimensional space-time (see Materials and Methods), represented by four yellow lines in Fig. 2E.From these trajectories, we extract the centroid velocities v c (indicated by pale gold arrows in Fig. 2), which correspond to the averaged v g that gauges the overall energy flow (see Materials and Methods).We find that v c is not parallel to v f (i.e., v g being not parallel to v p ; see Fig. 2, C and D), confirming the most prominent feature of HPs-the vectorial misalignment between the velocities v p and v g .
In Fig. 3A, we show in detail the trajectory of a specific centroid in the x-y space, revealing a curved path that goes against the ray-like trajectories usually observed in monochromatic imaging experiments (12).The observed pulse beams and their curved centroid paths can be demonstrated by our time-domain simulations, which can be explained by the coherent superposition of different frequency components of HP pulses (fig.S8).We further analyze this anomalous curved energy flowing by examining the centroid velocity v c .The orientation angle θ c of v c continuously increases with time during the pulse propagation (Fig. 3B), while the fringe velocity v f exhibits a nearly constant angle (|θ f | ⁓ 45°) taken from positions of the centroid.By comparing the θ f value with the theoretical HP dispersion, we infer that it corresponds to the frequency component ω = 1470 cm −1 (the center of the Reststrahlen band) of the pulse fields.The misalignment angle θ between v c and v f increases with time, up to more than 90°after τ ⁓ 0.8 ps, which corroborates the extreme HP anisotropy.By examining the timevarying magnitude of v c (the color gradient in Fig. 3A), we demonstrate the ultraslow energy flow of HP pulses with the velocity |v c | as slow as 0.005c (c being the speed of light in free space).Moreover, we unveil the dynamic phenomena of the energy-flowing acceleration (τ < 0.6 ps) and deceleration (τ > 0.6 ps), which can be explained by the frequency dispersion of HP pulses.The observations of the curved energy-flowing trajectories, along with the acceleration and deceleration, unambiguously demonstrate the potential of our nanoimaging approach in investigating the peculiar spatiotemporal physics of anisotropic polaritons.
Figure 3C delves into the anisotropic dissipations of HP pulses by showing two cross-section slices from the high-dimensional dataset (the one depicted in Fig. 2E).These slices reveal the timevarying near-field amplitudes of a single HP pulse beam in two azimuthal directions (15°and 30°, with respect to the x axis).In both directions, we observe the pulse packet (shown in purple) moving away from the disk (located at the coordinate origin) with increasing time.The maximum amplitudes (marked by white circles) do not appear at the disk, but rather at a distance, indicating the energy accumulation during propagation.To quantify the anisotropic decay, we extract the 1/e value of the maximum of the overall beam in the high-dimensional space-time (dashed lines in Fig. 3C; see also fig.S9) and derive the maximum values of the profiles along the position axis and the time axis as the maximum decay lengths L m and maximum decay times τ m (see symbols marked in Fig. 3C), respectively.It is remarkable that L m and τ m do not coincide (see symbols marked in Fig. 3C), indicating the large group velocity dispersion of HP pulses.
Figure 3D displays the extracted τ m for various directions (top panel), which promptly reveals the in-plane anisotropic decay of HP pulses, in good agreement with numerical simulations (bottom panel; see note S3 for details).Note that we did not consider the 2D spreading decay of disk-excited polaritons (6,8,9).Therefore, the experimentally obtained quantities are smaller than the intrinsic values.These spatiotemporal observations directly unveil the extreme space-time anisotropy and the strong group velocity dispersion of HP pulses, which hold potential applications in anisotropic waveguiding (21,22), nanoscale heat transfer (35), and peculiar light-matter interactions (36).In Fig. 3E, we demonstrate the capability of regulating the anisotropic propagation of HP pulses by altering their excitation direction with respect to the optic axis (OA) of the calcite crystal.We present two experimental near-field amplitude snapshots, acquired at different incident angles relative to the crystal's OA (i.e., the y direction).Our observations show that when HP pulses are excited parallel to the y axis, they exhibit propagation fringe patterns with mirror symmetry relative to the y direction.Conversely, varying the incidence angle to 45°results in asymmetry fringe patterns of HP pulses.Our simulation results validate the effectiveness of controlling the propagation symmetry and dynamics of HP pulses, showing excellent agreement with the experimental data.This regulating effect is attributed to the disk antenna acting as an in-plane dipole emitter that is able to selectively excite the momentum components of HPs [see details in (37)].These findings offer opportunities for shaping the propagation behavior of HP pulses, which could have implications for applications in topological photonics (38), ultrafast information transfer, and processing (19).
The spatiotemporal characteristics of polariton pulses are determined by their extremely anisotropic wave vectors k.In our experiments, we are able to monitor polariton wave vectors in the (k, t) domain, allowing us to demonstrate a time-dependent topological transition for isofrequency contours (IFCs) in k-space.We achieve this by using the 2D FT to convert near-field data into domains of in-plane momenta (k x and k y ) and time delay τ (as detailed in note S4, Fig. 4, A to H, and movie S4).In this analysis, we include the data with laser frequencies outside the Reststrahlen band.The processed dataset reveals a lenticular-shaped IFC at short time delays (τ < 0.25 ps), which undergoes a topological transition at around τ = 0.25 ps.During this transition, the lenticular IFC flattens out, and the hyperbolic IFC emerges at larger k values.Subsequently, the lenticular IFC disappears, and the remaining hyperbolic IFCs exhibit a decreasing opening angle with increasing time.We also analyze the polariton field confinement by obtaining the largest wave vector of the IFC as a function of τ, as shown in fig.S12.We show that this time-domain topological transition can be tuned by adjusting the initial phases of the illumination pulses, as discussed in fig.S13.Topological transitions arising from the modifications of the polariton dispersion in anisotropic materials enable the enhancement of the local density of photonic states for light-matter interactions, typically reported by static experiments (12,27).Our results represent the time-domain topological transition, which holds potential for various applications on ultrafast time scales, such as sensing, waveguiding, and imaging.
After characterizing the spatiotemporal dynamics of HP pulses and observing the time-domain topological transition, we demonstrate the ability to explore the optical pulling forces of HP pulses.Recent theoretical studies predict that variations in the wave vector k, specifically the transition topological of their IFCs, can induce optical pulling forces when light scatters off an object within a designed photonic environment (39).By directly accessing the timevarying wave vectors and light momentum P (with P = ℏk, where ℏ is the reduced Planck constant), we are able to experimentally evaluate the optical forces of polariton pulses induced on the disk antenna (note S5).In Fig. 4I, we show the distribution of the time-dependent variations of light momentum |∆P| at τ = 0.2 ps.We observe that the patterns of |∆P| consist of two parts, with the low-k lenticular parts experiencing a larger increment than the high-k hyperbolic ones.The patterns of |∆P| are asymmetric, which we attribute to the oblique incidence and imperfect alignment of the laser pulses.We estimate the resulting force induced on the disk by summing the changes in momentum for all wave vector components using the equation of where ∆N is the change for each wave vector component k (see note S5 for details).The direction of the force F is opposite to that of the momentum of the incident light (illustrated in Fig. 4J), representing a pulling force.Figure 4J shows that the value of | P ∆P| (blue line, indicating the magnitude of F) and the direction angle of the force F (red line) vary over time.Both values show maximum at time delays close to the topological transition (indicated by a dashed gray line), which is in line with the theoretical prediction that the topological transition can enhance the optical forces (39).The value of | P ∆P| reaches a maximum at around τ = 0.21 ps and subsequently decreases.The angle ϕ of the F direction shows its maximum at around τ = 0.31 ps.Taking into account the parameters of the pulse laser used in the experiment and assuming an efficiency of 1% for the antenna converting free-space pulses to polariton pulses (40), we estimate the maximum transient force F to be on the order of approximately 100 pN (see note S5 for details).

DISCUSSION
In summary, we have demonstrated the spatiotemporal characteristics of extremely anisotropic nanolight pulses using the timedomain interferometric nanoimaging technique.In combination with diverse polarized light excitations, such as circularly polarized light (41,42), we envision that our approach can be applied to reveal various exotic polariton dynamics.It also holds potential for multiphysics space-time analysis, as demonstrated by the evaluation of momentum topology-induced optical pulling forces exerted by polariton pulses on an object situated on a naturally occurring hyperbolic surface.The time-dependent evolution of polariton wave vectors and associated topological transitions observed in our study offers exciting opportunities for dynamically tailoring the local density of photonic states in various applications (19) and exploring high-dimensional time-varying optics (43,44).We anticipate that the use of artificial intelligence techniques (45) and machine learning algorithms (46) will become increasingly important to accelerate the acquisition of high-dimensional spatiotemporal data and improve the imaging quality and accuracy.The spatiotemporal dynamics of HP pulses demonstrated here are essential for developing anisotropy-driven ultrafast nanophotonic and thermal applications.

Materials and sample fabrications
We used 1 cm-by-1 cm-sized calcite substrates with characteristic (100) plane, which were commercially available by mechanical cleavage from bulk calcite single crystals.Gold disk antennas were fabricated on the calcite surface via standard electron beam lithography.The disk patterns were written on the resist (polymethyl methacrylate: 495/A4, thickness ≈ 100 nm) spin-coated on the calcite substrate.A conductive polymer (AR-PC 5090) was spincoated on the resist to prevent charge accumulation because of the weak conductivity of calcite.The standard lift-off procedure was carried out after the e-beam evaporation of Ti (3 nm)/Au (60 nm) onto the developed resist.

Time-domain interferometric measurements
For time-domain interferometric nanoimaging experiments, we used a commercial nano-FTIR system from Neaspec GmbH based on an atomic force microscope (AFM).The Pt-coated AFM tip oscillates vertically with an amplitude of about 70 nm at a frequency Ω ≈ 250 kHz (Arrow NCPt, NanoWorld).It is illuminated by a broadband laser pulse of 100-fs duration (range from 1000 to 2000 cm −1 ; see fig.S3).A piezo-controlled linearly moving mirror is used to precisely control the time delay τ between the tip-scattered field and the reference field.The 2D mappings were carried out at different time delays with a temporal resolution of Δτ = 16.7 fs (corresponding to the reference-mirror moving interval of 2.5 μm) satisfying the band-pass sampling (see note S1 for details).The interferometric detector signal was demodulated at a higher harmonic nΩ (n ≥ 2), yielding near-field amplitude s n and phase ϕ n images as a function of both the 2D (x, y) coordinates and time delay τ.In this work, we recorded the signal of the second order for producing a raw space-time imaging dataset (see Fig. 1B and fig.S1), e.g., the real part of near fields, E(x,y,τ) = s 2 cos(ϕ 2 ).We display the mappings of s 1 in Fig. 3E since the higher-order signal is weak for an incidence angle of 45°.

Polariton velocity vectors reconstruction
We used the Horn-Schunck (HS) algorithm to extract the pulse fringe velocities v f .The HS algorithm is well-known for computing the optical flow from image sequences (47), which was recently used to analyze the plasmonic flow in the photoemission electron microscopy experiment (48).It estimates a velocity vector field (v x , v y ) by minimizing the equation E HS ≈ ∫∫(I x v x + I y v x + I t ) 2 dxdy, where I x,y and I t are the spatial and temporal image intensity derivatives in each pixel, respectively.Specifically, we applied the HS algorithm to analyze the image sequences of the real part of filtered spatiotemporal data after normalizing them to the range of [0,1] for each frame, thus extracting the pixel-wise local distributions of time-dependent fringe velocities.In Fig. 3 (A and B), the fringe velocities v f at the centroid location were obtained by averaging fringe velocities in the range of 0.5 μm by 0.5 μm to make the trend smoother.
Taking the antenna center as the coordinate origin, we divided the 2D real space into four quadrants and processed each quadrant (i.e., each polariton beam), respectively.At each time delay, we removed the data whose absolute value is less than the 1/e value of the maxima of this quadrant so that the processed polariton beam is more distinct and less noisy.Then, we calculated the centroid position by the equation x c ¼ P I i x i P I i and y c ¼ , where I i is the absolute value of processed near-field data in each pixel, x i and y i are the corresponding coordinates.We extracted the centroid velocities v c by calculating the centroid movement between time steps, i.e., v c ðτÞ ¼ pðτÞÀ pðτÀ dτÞ dτ , where p(τ) is the position vector (x c (τ), y c (τ)) at each time delay τ.

Numerical simulations
We used a finite-difference time-domain simulation software (Lumerical FDTD) to simulate the dynamics of near-field distributions (E z ) of polariton pulses propagating along the calcite's surface.
To validate the measured space-time features of HP pulses, we simulated disk-launched polariton distributions when a gold disk (height, 50 nm; diameter, 1 μm) was placed on top of a calcite substrate (with an angle of 23.3°between the OA and the surface).A ppolarized broadband plane-wave pulse was directed obliquely onto the sample, illuminating a gold disk and launching HP pulses.The incident angle is 30°relative to the surface, and the OA of calcite is in the plane of incidence.The spectrum of the broadband pulse is Gaussian-shaped and central at 1500 cm −1 , approximating the infrared pulse used in the experiment.The horizontal simulation range covers 50 μm by 50 μm, and Bloch boundary conditions are applied.Note that the simulation range is sufficiently large to avoid interference that could arise from the boundary.We recorded the time-varying near-field distributions (E z ) at 100 nm above the surface.Simulated results shown in Fig. 3E were obtained by changing the illuminating orientation.
To study the decay time of HPs in different lateral directions shown in Fig. 3D, we monitored the time-dependent near-field distributions at 50 nm above the calcite substrate with a dipole (polarized perpendicular to the surface) exciting at 500 nm above the sample.Perfectly matched layers are used as boundary conditions.

Supplementary Materials
This PDF file includes: Notes S1 to S5 Figs.S1 to S13 Legends for movies S1 to S4 References Other Supplementary Material for this manuscript includes the following: Movies S1 to S4

Fig. 1 .
Fig. 1.Time-domain interferometry nanoimaging of HP pulse propagation.(A) Experimental setup.IR, infrared; BS, beam splitter.Note that the incident laser pulses are perpendicular to the propagating direction of polaritons (see the illustration) for suppressing the asymmetry of polariton fields excited on either side of the gold disk.(B) Experimental near-field snapshots (amplitude signal s, raw data).The entire raw datasets are up to 80 frames (fig.S1).

Fig. 2 .
Fig. 2. Space-time mapping of HP pulse velocity vectors.(A and B) Two background-subtracted experimental near-field snapshots showing the velocity pseudovector textures of HP pulses (see the complete movie shown in movie S2).Small colored arrows show local distributions of fringe velocity v f , with the colors showing the real part of the near fields of HPs.The large pale gold arrows indicate the centroids of four HP beams, respectively.Their directions are aligned with the corresponding centroid velocity v c .(C and D) Zoom-in images taken from the dashed boxes in (A) and (B), respectively.(E) Time-dependent background-subtracted experimental nearfield amplitude snapshots revealing the space-time evolution of butterfly-shaped polariton beams (the complete movie provided in movie S3), giving rise to centroid trajectories of each beam (yellow lines).Arrows indicate the directions of v c in each moment.

Fig. 3 .
Fig. 3. Anisotropic spatiotemporal dynamics and incidence-dependent control of HP pulse propagation.(A) Centroid trajectory of the HP pulse beam (marked by a dashed box in Fig. 2E).The color gradient of the trajectory curve indicates the magnitude of the centroid velocity v c that is obtained from time intervals.Blue and yellow arrows indicate the direction of the fringe velocity v f and the centroid velocity v c obtained at the centroid, respectively.Their orientation angles θ c and θ f are defined with respect to the x axis.θ is the angle between v c and v f .(B) Time-dependent variations of θ c , θ f , and θ during the HP pulse propagation.(C) Space-time evolution slices of near-field amplitudes of the HP pulse beam for two different propagating directions (the angle α defined according to the inset).Black dashed lines show the profiles where the pulse amplitudes decay to 1/e of the maximum of the whole beam.Maximum decay lengths L m and maximum decay times τ m are marked by blue squares and red rhombs, respectively.White circles indicate the maximum amplitude in each slice.(D) Experimental anisotropic decay time of HP pulses (top), showing good agreement with numerical simulation (bottom).The slight difference between the left and right in the experimental result is attributed to the imperfect alignment of the laser pulses with respect to the y direction.(E) Near-field distributions of HP pulses launched in different excitation directions (illustrated by arrows).Left: Experimental nearfield amplitude snapshots.Right: Simulated near-field amplitude snapshots.

Fig. 4 .
Fig. 4. Time-domain topological transition and time-varying optical pulling forces.(A to H) Time-dependent wave vector distribution of HP pulses revealing the IFCs from lenticular to hyperbolic.The time-domain topological transition occurs at around τ = 0.25 ps when the hyperbolic IFCs gradually dominate.Note that the antenna excites the polariton pulses composed of broad frequency components, resulting in the simultaneous presence of both hyperbolic and lenticular polariton components.(I) Momentum variation at τ = 0.2 ps.The blue dashed arrow shows the orientation of summing momentum variation and the red arrow indicates the direction of the corresponding optical force.(J) Summing momentum variation (blue line) and the orientation of optical forces (red line, the angle ϕ defined according to the inset) as a function of time.The gray dashed line marks the time when the topological transition occurs.a.u.arbitrary units.