Plasmonic vortex is an optical field distribution with topological features formed by interfering surface plasmons (SPs), which enriches the class of vortex phenomena in nature^{1, 2}. Owing to their special orbital angular momentum (OAM) feature in the evanescent field region^{3, 4}, plasmonic vortices hold great promises for many cutting-edge applications, such as plasmonic tweezers for microparticle manipulations^5-12 and on-chip quantum information processing^{13, 14}. The generation methods and evolution dynamics of plasmonic vortices have thus elicited great research enthusiasm in the last decade. Inspired by the development of meta-optics^15-21, artificial structures have recently emerged as an ideal tool for the generation of plasmonic vortices. The most common method is constructing special couplers, such as Archimedean spiral-shaped long slits^{2, 22-25} and subwavelength slit resonators^26-30, to convert circularly polarized light carrying spin angular momentum into on-chip plasmonic vortex, utilizing the design degrees of freedoms of propagation phase and geometric phase. In parallel, state-of-the-art near-field microscopes including photoemission electron microscopy (PEEM)^31-36, nonlinear near-field optical microscopy^{37, 38}, and so forth^39-42 were successfully employed to image evanescent electromagnetic waves, which have also empowered the sub-optical-cycle study of plasmonic vortex evolution dynamics^{23, 31, 43}.

The above two branches of studies have provided many insights into the nature of plasmonic vortex and rapidly promoted the related applications forward^{7, 31}. Nevertheless, there are still some puzzles or shortages in this field. For plasmonic vortex generation, despite the sole or combined use of propagation phase and geometric phase can all achieve plasmonic vortex of target topological charge^{28, 29}, the actual differences of their spatiotemporal dynamics have remained unexplored. For characterization methods, the PEEM suffers from nonlinear spin-orbit coupling between the probing circularly polarized pulse and the propagating plasmonic vortex^{23, 31}, which can hardly obtain the exact evolution dynamics. On the other hand, the nonlinear near-field optical microscopy delivers evanescent field information to free-space far-field through nonlinear generation³⁸, whose operating wavelength range and spatiotemporal resolution are thus limited by the applied nonlinear process and the optical system itself.

Recently, with the developments of terahertz (THz) near-field microscopes^{42, 44-47}, THz frequencies have become a promising platform to explore fundamental properties and potential applications of SPs⁴⁸, where the THz SP fields can be imaged with very high spatial and temporal resolution. In this article, based on a near-field scanning terahertz microscopy (NSTM)⁴⁷, we experimentally demonstrate how to tailor the spatiotemporal dynamics of plasmonic vortex. More specially, the plasmonic vortices of the same topological charge can be endowed with distinct generation and evolution behaviors, as depicted in Fig. 1, where SP fields reach and decay simultaneously (Fig. 1(a) to 1(c)) or successively (Fig. 1(d) to 1(f)) at the same revolution orbit. As proof-of-concept experiments, we fabricated samples in THz frequencies and characterized them using NSTM. The full amplitude and phase information of SP fields and the exact evolution dynamics with ultrahigh temporal resolution were directly obtained. Our results introduce a new degree of freedom to tame OAM delivered by plasmonic vortices, i.e., the spatiotemporal dynamics, which is of fundamental importance for understanding the nature of plasmonic vortex and for designing practical applications in this field.

It has been demonstrated that, if a pair of slit resonators are separated by a distance of λ_SP/2 and arranged perpendicular to each other, under normal incidence with circular polarization (CP), the phase of the excited SPs can be tuned in the whole 2π range (λ_SP is the SP wavelength)^{49, 50}. Considering an Archimedean spiral-shape arranged plasmonic vortex coupler composed by m slit-pairs, as illustrated in Fig. 2(a). This coupler can be described by the geometric phase factor g and the propagation phase factor p. Firstly, g determines the orientation angles of the constituent slit resonators:

1θ1,m=gφm/2+π/2+α0,θ2,m=gφm/2+α0,

where, θ_{1, m} and θ_{2, m} represent the orientation angle of the inner and outer slit resonator of the mth slit-pair, respectively; φ_m is the azimuthal angle of the mth slit-pair with respect to the x-axis. α₀ is the initial orientation angle of the first outer slit resonator. Under the CP incidence, the overall phase shift of SPs propagating towards the center across a whole turn equals to 2σ(g − 1)π, where σ∈{+,−} stands for the spin-direction of the left- and right-handed CP (LCP and RCP)²⁸. Secondly, p determines the radii of the constituent slit-pairs:

2Rm=R0+pλSPφm/2π,

where, R_m is the radius of mth slit-pair and R₀ is the initial radius when φ_m= 0. As radius variation increases from 0 to pλ_SP, a phase gradient of 2pπ is generated along the azimuthal direction. Based on these designs, the phase of SP field at the center generated by mth slit-pair can be expressed as E∝ei(σ(g−1)+p)φm(see Supplemental information Section 1 for details). As thus, once g and p are designated, the plasmonic vortex coupler design is determined, where the topological charge of the generated plasmonic vortex can be expressed as l = σ(g − 1) + p. Clearly, under a given CP incidence, for generating plasmonic vortex of target topological charge l, there are infinite solutions of g and p combinations.

Each slit resonator functions as an individual SP source that can only be excited by the polarization component perpendicular to its long side and then produce SPs with in-plane dipolar feature^{28, 51}. Based on 2D Huygens-Fresnel principle⁵², the propagation behaviors of the excited SPs can be well predicted. Taking a single slit resonator as instance, which is orientated by an angle of θ with respect to the x-axis, as depicted in Fig. 2(b). In frequency-domain, under the CP normal incidence Ein=(2/2σi2/2), the SP field at an arbitrary point P can be calculated as:

3EP=η22exp(ikSP|r)|iλSP|r|(cosθ+σisinθ)cos(θ−ξ),

where, E_P is the complex amplitude of the SP field at point P; η is the conversion efficiency from incident light to SPs; k_SP = 2π/λ_SP is the SP wave number; r is the vector from the slit resonator position to the point P; ξ is the angle of r with respect to the x-axis. For investigating the frequency-domain performance of a given plasmonic vortex coupler, one can obtain the corresponding SP intensity and phase distributions by calculating the superposition of the SP excitations from all the constituent slit resonators.

In order to numerically reveal the spatiotemporal dynamics of plasmonic vortex formation, revolution, and decay processes, we generalize the 2D Huygens-Fresnel principle to time-domain. Suppose the excitation signal is an RCP pulse, which is composed by an x-polarized pulse and a y-polarized pulse with a relative time delay of Δ = T/4, as depicted in Fig. 2(c), and T = 1/f₀ represents the optical-cycle at f₀. Here, we choose f₀ = 0.75 THz to obtain RCP incidence at this frequency. In principle, the resonance of the slit resonator can be excited by the perpendicularly polarized component and then produce radiation forward, wherein part of the energy is coupled into SP excitations. When the slit resonator is orientated by an angle of 0 and π/2, as shown in Fig. 2(d), it can be only excited by the x- and y-component of the incident RCP pulse. Correspondingly, their radiation waveforms, S_x and S_y, are the same with a relative time delay of Δ. For a slit resonator with arbitrary orientation angle θ, its radiation waveform can be calculated as:

4Sθ=Sxcosθ+Sysinθ.

Next, based on the dipolar response of slit resonator and 2D Huygens-Fresnel principle, the SP amplitude at point P (as depicted in Fig. 2(b)) and moment t_s can be calculated as:

5AP(ts)={Sθ(ts−|r|/v)cos(θ−ξ)|r|,ts−|r|/v⩾00,ts−|r|/v<0,

where, v = c₀k₀/k_SP is the propagation speed of the SP waves with c₀ and k₀ representing the speed and wave vector of light in vacuum; |r|/v describes the required time that the SP waves propagate from the slit resonator to point P, clearly A_P = 0 whents−|r|/v<0;1/|r| describes the propagation attenuation of SP waves; cos(θ − ξ) describes the in-plane dipole source feature of slit resonator. Since the propagation distance is only several wavelengths, here we ignore the waveform variation caused by SP dispersion (see Supplemental information Section 2 for details). In this manner, one can numerically obtain the time-resolved SP field snapshots of a certain plasmonic vortex coupler, by calculating the amplitude superposition from all the constituent slit resonators over an area of interest.

Here, we numerically investigate the frequency- and time-domain performances of two plasmonic vortex couplers. Without loss of generality, they are designated of different p and g combinations to generate plasmonic vortices of the same topological charge of l = 4 under the RCP incidence, as shown in Fig. 3(a) and 3(d). Sample 1 is designed by introducing only the geometric phase through varying the orientation angle of slit resonators, corresponding to p = 0 and g = −3. Sample 2 employs the propagation phase by a varying radius of δR = 4λ_SP, corresponding to p = 4 andg = 1. The two samples both consist of 70 slit-pairs, where the orientation angle and spatial position of each constituent slit resonator can be calculated by using Eqs. (1) and (2), with the initial setting of α₀ = 0 and R₀ = 1 mm.

We first investigate their frequency-domain performances at 0.75 THz based on Eq. (3). The calculated SP intensity fields are shown in Fig. 3(b) and 3(e), which are similar to each other and exhibit concentric rings shaped intensity distributions. The intensity dependence along the radial direction can be approximately described by a Bessel function of order 4, that is, I(R)∝|J4(kSPR)|2^{53, 54}. The innermost bright intensity ring can be treated as the target orbit of 4th-order plasmonic vortex, whose radius is about 340 μm, corresponding to the first extremum of the Bessel function. As for the phase distributions (see Fig. 3(c) and 3(f)), they both increase clockwise and linearly from 0 to 8π along the target orbit ring, indicating the 4th-order vortex nature of the excited SP fields. These results together numerically confirm that different combinations of p and g can generate plasmonic vortices of the same topological charge, and the generated SP fields are similar in terms of both intensity and phase distributions. For plasmonic vortex related applications, in addition to frequency-domain performances, the actual dynamics in time domain including the formation of field patterns that carry OAM and their ultrafast revolution at target orbit are equally important.

Next, we investigate the spatiotemporal dynamics of Sample 1 and Sample 2 based on Eqs. (4) and (5). Figure 3(g) and 3(i) illustrate 11 time-resolved snapshots of the SP amplitude field distribution in the xy-plane, respectively, of Sample 1 and Sample 2. Here, t₀ represents the initial moment when SPs began to be excited from the slit resonator. These snapshots cover the formation, revolution, and decay stages in the lifetime of the generated plasmonic vortex³¹. For Sample 1, since the radial distances in all direction keeps consistent, the initial SP excitations synchronously form 8 converging spiraling wavefront threads (snapshots 1–3 in Fig. 3(g)). This is followed by the concentration of the SPs into the target orbit, where the inward and outward counter-propagation SPs interfere to form the radially standing but azimuthally rotating vortex field (snapshots 4–6 in Fig. 3(g)). At this stage, 8 rotating lobes (4 with positive maximum shown in red and 4 with negative maximum shown in blue) can be observed, corresponding to the 4 azimuthal wavefronts of the steady-state e^4iφ phase of the vortex. What follows is the decay stage (snapshots 7–11 in Fig. 3(g)), where the SP fields form outward-propagating spiraling wavefront and then dissolve. Notably, the handedness of the wavefront flips compared to the formation stage.

For Sample 2, the spatiotemporal dynamics are quite different. At the formation stage, as shown by snapshots 1–4 in Fig. 3(i), the converging spiraling wavefront threads successively appear and concentrate centripetally, since the SPs excited from Archimedean spiral-shape arranged slit resonators require different time to propagate to the target orbit. This also leads to the followed interference between inward and outward counter-propagation SPs only occurs at part of the orbit, forming uncomplete but azimuthally rotating vortex field (see snapshots 5–7 in Fig. 3(i)). At this stage, the 8 lobes cannot be seen at the same time. What remains is the handedness flip of the spiraling wavefront appearing in the decay stage (snapshots 8–11 in Fig. 3(i)).

In order to illustrate the spatiotemporal dynamics difference between Sample 1 and Sample 2 more clearly, we extract the absolute amplitude value of SPs on the target orbit. For Sample 1, as shown in Fig. 3(h), the SP fields of 8 lobes reach the target orbit synchronously with time, and gradually increase to their maximum then decrease to zero. During these processes, the SP fields collectively rotate along the orbit, forming relatively uniform OAM in both space and time. For Sample 2, as shown in Fig. 3(j), the SP fields first arise at the top-right side of the orbit, and the amplitude gradually increases with the SP lobes rotating counterclockwise. When the SP lobes rotate to the bottom-right side, the amplitude gradually diminishes and becomes zero ultimately. The SP fields vary successively during the whole evolution process, corresponding to a decomposition of the OAM in both space and time. The above calculation results are acquired using the MATLAB software and these results together verify that, with similar frequency-domain performances, the spatiotemporal dynamics of plasmonic vortices can be tailored very differently.

In order to experimentally verify the above numerical results, we fabricated the realistic samples and characterized them using the fiber-based NSTM system. The samples of plasmonic vortex couplers were fabricated by conventional photolithography, metallization, and spin-coating processing. First, thermal evaporation was employed to deposit a 200-nm-thick Aluminum film on the 2-mm-thick double-side polished high-resistivity silicon wafer. Then a layer of photoresist (AZ P4000) was spin-coated on the Aluminum film. The patterns of the plasmonic vortex couplers were exposed by standard lithography. The photoresist and Aluminum in the exposed area were removed by development and etching processes, respectively. Finally, the remaining photoresist in other area was washed away by acetone, leaving an Aluminum film perforated with structures of plasmonic vortex couplers on the silicon wafer. In order to enhance the confinement of the SPs on the metal surface, a 10-μm-thick polyimide layer was coated on the whole structure. The dimension of the slit resonators was set to be 70 μm × 10 μm, working at 0.75 THz, corresponding to λ_SP ≈ 400 μm.

The fiber-based NSTM was employed to characterize the generated plasmonic vortices. As schematically shown in Fig. 4, it is mainly composed of two parts, the transmitter modules and detector modules. In-fiber femtosecond pulse with ~50 fs pulse width and 1550 nm central wavelength is split into two beams. One of the femtosecond beams illuminates onto a commercial photoconductive antenna to generate broadband terahertz radiation. The other femtosecond beam passes through a frequency doubling crystal (BBO) and illuminates onto a commercial terahertz near-field probe⁵⁵ to detect the SP fields, in which the movement of excited carriers forms a current proportional to the detected terahertz SP field. The pump pulse at the emitter is delayed with regard to the one at the probe thanks to a fiber-based delay line (FDL) which is used to change the optical path. The basic idea of NSTM is to temporally read the waveform, step by step, by changing the optical delay line to sample the SP field at different moments and eventually obtain the whole time-domain signal⁵⁶. In experiment, the sampling step was chosen as 20 fs corresponding to ~1/66 of the optical-cycle at 0.75 THz, and the whole sample time of each signal is 60 ps. On the other hand, the terahertz probe is fixed onto a 2D electrically controlled translation stage to enable the scanning function of high spatial resolution. In experiment, all the measurements were obtained at the plane that about 75 μm above the sample surface, and the scanned range was set to be 2 mm × 2 mm with a 40 μm step.

Figure 5(a) and 5(d) show the microscopic photos of fabricated realistic Sample 1 and Sample 2, respectively. In experiment, as shown in Fig. 4, the generated terahertz radiation was first collimated by a terahertz lens and then passed through a terahertz polarizer (LP₁) to form a linearly polarized terahertz beam; through rotating another polarizer (LP₂) to be +45° and −45° with respect to LP₁, the experimental results under two orthogonal linear polarized incidences can be obtained, which can be used to composite the experimental results under CP incidence. More specifically, under +45° and −45° linearly polarized incidences, by controlling the optical delay line and the 2D electrically controlled translation stage, the SP field waveforms Sx,y+45∘(t) andSx,y−45∘(t) on the positions of a 2D lattice (x, y) can be obtained; the SP field waveforms under CP incidence can thus be composited asSx,yσ(t)=[Sx,y+45∘(t)+Sx,y−45∘(t+σΔ)]/2. Next, the broadband complex amplitudes in frequency-domain can be obtained by performing the Fourier transformEx,yσ(f)=F{Sx,yσ(t)}. Figure 5(b, c) and 5(e, f) illustrate the experimentally obtained SP intensity fields and phase distributions of Sample 1 and Sample 2, respectively, at 0.75 THz under the RCP incidence. Clearly, these experiment results agree well with the numerical results shown in Fig. 3(b), 3(c) and 3(e, f), forming similar plasmonic vortices from Sample 1 and Sample 2. To quantitatively analyze the quality of the generated plasmonic vortices, we extracted and integrated the complex amplitudes along the target orbit to get the OAM amplitudes (see Supplemental information Section 3 for details).Figure S3(a) and S3(b) respectively show the obtained OAM amplitudes of Sample 1 and Sample 2, where the target OAM component of order 4 is the strongest in both numerical and experimental results, and the other OAM components are quite weak. As for the generation efficiency of plasmonic vortices, it is mainly dependent on the coupling efficiency of adopted unit element, i.e., the slit resonator and cannot be directly measured in the experiments currently. In our previous work⁴⁹, the coupling efficiency of slit-pair constructed surface plasmon diffraction grating is estimated to be about 7.6% with 3.8% toward each side. Since the plasmonic vortex coupler can be seen as a bended surface plasmon diffraction grating, it is believed to have the similar coupling efficiency. To increase the coupling efficiency, one could apply the metal-dielectric-metal⁵⁷ or multiple-layered⁵⁸ resonators.

Figure 5(g, h) and 5(i, j) illustrate the experimentally obtained spatiotemporal dynamics of Sample 1 and Sample 2, respectively, including the snapshots of SP fields in xy-plane and the extracted absolute SP amplitudes on the target orbit. In addition to the fabrication and measurement induced errors, the deviations between numerical (Fig. 3(g, h) and 3(i, j)) and experimental (Fig. 5(g, h) and 5(i, j)) results can be mainly attributed to two aspects: the slight difference between the numerically generated SP waveform and the actual SP waveform generated from the slit resonator (see Supplemental information Section 4 for details); the near-field couplings between adjacent slit resonators also cause the slit resonator deviates from an ideal dipolar SP source. Nevertheless, the main features of experimental results complied well with that of numerical results, and experimentally verified that the plasmonic vortices of similar frequency-domain performances can be endowed with distinct spatiotemporal dynamics. It should be emphasized that since the waveforms Sx,y±45∘(t) were experimentally obtained with consecutive time steps of 20 fs, it is possible to reveal the plasmonic vortex evolution in deeper sub-optical-cycle tempol resolution. Supplemental Videos S1 and S2 (see Supplemental information Section 5) are provided to show the detailed revolution behaviors of the plasmonic vortices generated by Sample 1 and Sample 2, respectively, corresponding to the results of snapshots 3–6 in Fig. 5(g) and snapshots 3–7 in Fig. 5(i). These videos are set as 10 Hz and the interval between each frame corresponds to 20 fs in experimental results.

We have demonstrated that the spatiotemporal dynamics of plasmonic vortex can be manipulated by changing the combination of propagation phase and geometric phase in coupler design. In addition to the demonstrated samples, Supplemental information Section 6 illustrates the experimentally obtained frequency- and time-domain performances of another three samples. As demonstrated in Fig. S5, they are designed with p = 1, 2, 3 and g = –2, –1, 0, respectively. In frequency-domain (see Fig. S6), they can all form plasmonic vortex of the same topological charge l = 4. Whereas in time-domain, as shown in Fig. S7, their spatiotemporal dynamics are different from each other, corresponding to the intermediate transition cases between that of Sample 1 and Sample 2. On the other hand, the resonance response of slit resonator can also serve as another design degree of freedom in tailoring the spatiotemporal dynamics of plasmonic vortex, since it directly determines the waveform of the excited SPs. For reference, we carried out numerical studies in Supplemental information Section 7, where the coupler designs are the same as that of Sample 1 and Sample 2, but with two different excitation waveforms. Figures S9 and S10 illustrate the corresponding results including the snapshots of SP field in the xy-plane and the extracted absolute SP amplitudes on the target orbit. The excitation waveform mainly affects the revolution stage, that is, the more oscillations of the waveform the longer rotation period of the SP lobes along the target orbit. Specially, for designs like Sample 2, more oscillations of the waveform lead to more SP lobes that can be observed in the same snapshot.

For characterizing the spatiotemporal dynamics of surface plasmons, the main limitation of PEEM^{23, 31, 37, 43} is that it can only detect the in-plane electric components of surface plasmons, since it needs to be co-polarized and in phase with that of probe pulse to produce the two-photo photoemission process. However, for surface plasmons, the electric components perpendicular to the surface are of more interests. In particular, for charactering plasmonic vortices, when their in-plane electric components are radial, the PEEM has to use circular polarized pulse as probe light^{23, 31}. This will lead to a nonlinear spin-orbit coupling between the probing circularly polarized pulse and the propagating plasmonic vortex and it will be hard to obtain the exact evolution dynamics. More specifically, the number of detected spatiotemporally varying surface plasmon lobes will be one more or one less than its topological charge, depending on the actual nonlinear spin-orbit coupling process. In contrast, the NSTM we adopted can directly detect the electric components perpendicular to the surface, giving the exact evolution dynamics of surface plasmons.

Our research on the distinct generation and evolution behaviors of plasmonic vortices is of great importance for the practical applications related to time-varying characteristics. One of the most appropriate occasions is to be applied into the plasmonic tweezers⁷, which can overcome the diffraction limits and enable the trapping of nanoscale objects without contact, providing a way to manipulate particles over a wide length scale. In addition to the linear momentum that generates a push or pull force, plasmonic tweezers can also import angular momentum to nanoscale objects which tends to rotate the objects⁵. The electromagnetic force in the plasmonic tweezers can be further separated into gradient forces related to the field distribution and scattering forces related to the momentum transfer. Therefore, if the field distribution and plasmonic OAM evolution can be tailored in both the temporal and spatial domain, one can allocate the electromagnetic forces at will to manipulate particles with a new degree of freedom and higher precision.

In summary, we provide an experimental glance into the spatiotemporal dynamics of the plasmonic vortex and the corresponding evolution of OAM in the THz region. We verified the feasibility of tailoring plasmonic OAM in time-domain based on different plasmonic vortex coupler designs. We envision that our work will provide guidance for plasmonic vortex related applications, such as the torque design for plasmonic spinners/tweezers and the light-matter interaction design for on-chip information processing. In addition, the proposed plasmonic OAM manipulation strategy is general and can be directly applied to the infrared and visible regimes, providing a new degree of freedom to explore more intrinsic nature and potential applications of plasmonic vortices.