### Experimental setup for THz near-field imaging

Figure S1 (A and B) schematically shows the experimental setup. We used a mode-locked Ti:Sapphire regenerative amplifier (Solstice, Spectra-Physics) that delivers 100-fs optical pulses (center wavelength, 780 nm) with 3.2-mJ pulse energy at 1-kHz repetition rate. Linearly polarized coherent THz pulses with Gaussian beam profile were generated by optical rectification in a LiNbO_{3} (LN) crystal using tilted-pulse-front excitation scheme (fig. S1C) (*28*, *29*). Two lenses (lens 1: Tsurupica lens, *f* = 50 mm; lens 2: 3-mm-thick, 2-mm-radius silicon bullet lens) were used to match the spot size of the excitation beam to the size of the metallic structure for efficient excitation. We used 1-μm-thick x-cut LN mounted on a glass substrate as an electro-optic (EO) crystal to detect an electric field of THz pulses (*30*). 2D EO imaging was performed using a probe beam with a large spot size (*19*). The image of the probe beam at the EO crystal was relayed to a 16-bit complementary metal-oxide semiconductor camera from PCO (model pco.edge) with a polarization analyzer unit. To perform EO sampling in the reflection geometry, the top and bottom surfaces of this substrate are coated with high-reflection (HR) and anti-reflection (AR) films (3 μm thick) for the probe pulse at 780 nm, respectively. The short pass (SP) and long pass (LP) filters were used to enhance the detection sensitivity using the probe spectrum filtering technique (*20*). An x-cut LN with a thickness of 2 μm is used to compensate for the birefringence of the EO crystal. When performing vortex beam experiments, we put an SPP in the collimated part of the THz beam propagating in the *z* direction (*31*–*33*). We made two SPPs to generate vortex beams with the *z* component of the OAM of ±*ħ* or ±2*ħ* centered at around 0.5 THz (SPP-1 and SPP-2). We used ZEONEX (cyclo-olefin polymer; refractive index, 1.52) as a material. SPP-1 is nominally designed to convert the Gaussian beam into the vortex beam (±*ħ*) at 0.45 THz. The step height is 1.29 mm. SPP-2 has a step of 2.18 mm, which generates the vortex beam (±2*ħ*) at 0.53 THz. The total step heights shown above are discretized to 16 small steps. The sign of the OAM can be inverted by reversing the SPP in the beam path.

### Sample fabrication

Corrugated gold disks are fabricated on the top surface (high-reflection coating material) of the THz detector crystal by a standard photolithographic technique. Chromium (10 nm thick) is deposited under the gold as an adhesion layer.

### Azimuthal angle dependence of the local electric field around multipolar LSPs

As schematically shown in Fig. 2 (B, D, and F), the surface charge distributions of LSPs along the circumference (radius *r* = *r*_{0}) can be expressed as follows

(1)

Here, *l* = 1, 2, … corresponds to the dipole, quadrupole, and higher-order LSPs, respectively. We ignore the time dependence (rotation of the charge distribution) at this moment and calculate the static electric field, **E**_{l}(*r*, φ). The scalar potential Φ(*r*, φ) can be obtained by solving the Laplace equation (at *r* ≠ *r*_{0}) in cylindrical coordinates (*r*, φ, *z*). Because of the cosine-type charge distribution, the scalar potential is an even function about φ. We first consider the 2D case. When the scalar potential is independent of *z*, the general solution is

(2)where *A*_{0}, *B*_{0}, *A _{l}*, and

*A′*are coefficients determined by boundary conditions. To avoid the divergence of the electric field in the limit of

_{l}*r*→ 0 and

*r*→ ∞, the electric field distributions have to take the following forms

(3)

From the boundary condition, [**E**(*r* ≥ *r*_{0}, φ) − **E**(*r* ≤ *r*_{0}, φ)]**·e**_{r} = ρ* _{l}*(φ) at

*r*=

*r*

_{0}, we obtain the following result for the static electric field around ρ

*(φ),*

_{l}(4)

By taking the projection onto the detection polarization axis **e**_{0} = **e**_{r}cos(φ) − **e**_{φ}sin(φ), the azimuthal angle dependence of the electric field detected in the experiment is

(5)

When ρ* _{l}*(φ) rotates as ρ

*(φ,*

_{l}*t*) = cos(

*l*φ ∓ 2π

*ft*) = cos(±

*l*φ − 2π

*ft*), additional electric field components such as radiation field and induced field are generated. However, in the near-field regime, the following quasi-static field is the dominant one

(6)

By replacing ±(*l* + 1) by *m* = ±2 (dipole), ±3 (quadrupole), and so on, we have

(7)

This shows that different LSP modes correspond one-to-one to cosine functions that are orthogonal to each other. Similar calculation for *l* = 0, i.e., a uniform surface charge ρ_{0}(φ) ∝ const., leads to the following form

(8)

These functions *E _{m}*(φ,

*t*) (

*m*= ±1, ±2, ±3, …) plus

*E*

_{0}(φ,

*t*) = const. form a complete set for symmetric field distributions about φ = 0 with the period of 2π. In the actual 3D case (100-nm-thick disk), the electric field distribution starts to deviate from the simple cosine functions (Eq. 7), as the detection plane gets away from the center of the sample along the

*z*direction. In our experiment, this effect is minimized by measuring the electric field distribution in the plane as close to the sample as possible. Therefore, we assumed that the electric field distribution can be approximately expressed by Eq. 7 and performed the mode expansion analysis. The effectiveness of this approximation is supported by the excellent agreement between the frequency spectrum obtained from the experimental data (based on the 2D model) and the numerical simulation (full 3D) shown in the Supplementary Materials (fig. S7).

### Dispersion relation of spoof LSPs

The dispersion relation of the spoof surface plasmons that propagate on the periodically corrugated metal surfaces is given as follows (*34*, *35*)

(9)

Here, *k* is the wave vector along the surface and *c* is the speed of light. The asymptotic (cutoff) frequency in the limit of *k* → ∞ is *c*/4*n*_{g}(*R* − *r*). By folding the surface to make a corrugated metallic disk, *k* becomes the wave vector along the circumference, which is quantized because of the periodic boundary condition. We used Eq. 9 as a fitting function in Fig. 3D. The fitting parameters determined by the method of least squares are *n*_{g} = 4.8 and *n*_{out} = 2.1. These numbers represent the effective refractive indices for the electromagnetic waves inside the groove and outside the disk. These should be expressed by combinations of the refractive indices of air (*n* = 1.00), silicon (*n* = 3.41), glass (*n* = 1.96), high-reflection coating materials (*n* ~ 2), and LN (*n* = 5.11) with appropriate weightings depending on the spatial extent of the electric field (see fig. S1B for the structure around the metallic disk). The high effective index inside the groove (*n*_{g} = 4.8) indicates that the electromagnetic mode is strongly localized around the LN (*n* = 5.11). The effective index outside the disk (*n*_{out} = 2.1), which is close to the index of the glass substrate (*n* = 1.96), suggests that the localization is weaker.