Near-field transmission matrix microscopy for mapping high-order eigenmodes of subwavelength nanostructures

As nanoscale photonic devices are densely integrated, multiple near-field optical eigenmodes take part in their functionalization. Inevitably, these eigenmodes are highly multiplexed in their spectra and superposed in their spatial distributions, making it extremely difficult for conventional near-field scanning optical microscopy (NSOM) to address individual eigenmodes. Here, we develop a near-field transmission matrix microscopy for mapping the high-order eigenmodes of nanostructures, which are invisible with conventional NSOM. At an excitation wavelength where multiple modes are superposed, we measure the near-field amplitude and phase maps for various far-field illumination angles, from which we construct a fully phase-referenced far- to near-field transmission matrix. By performing the singular value decomposition, we extract orthogonal near-field eigenmodes such as anti-symmetric mode and quadruple mode of multiple nano-slits whose gap size (50 nm) is smaller than the probe aperture (150 nm). Analytic model and numerical mode analysis validated the experimentally observed modes.

V arious types of resonant optical interactions with nanostructures have been exploited in nanophotonics to engineer the performance of nanoscale devices for applications ranging from sensing 1 to optoelectronics and metamaterials [2][3][4] . By tailoring the collective optical modes associated with the resonant interactions 5 , important figure of merits such as the coupling efficiency of the far-field illumination to the near-field, local field enhancement, spatial field localization, and control of nanoscale light flow 6 have been significantly improved. In these studies, NSOM has served as a general and robust tool to map the local near-field for unveiling the underlying process of nanoscopic physical phenomena and optimizing the performance of nanophotonic devices 7 . It has been especially useful for visualizing the formation of these optical modes and their coupling with the other spatially and spectrally overlapped modes [8][9][10][11][12] .
In NSOM, subwavelength probes are used to convert the nearfield evanescent waves carrying spatial frequencies exceeding the far-field limit set by k 0 = 2π/λ, where λ is the wavelength of the light source 13 , to the detectable far-field waves or vice versa. To map the optical eigenmodes in nanostructures, the probe is placed near the device to either collect local near-field electromagnetic wave or launch a spatially localized electromagnetic wave 8 . Importantly, the excitation wavelength is tuned to the resonance of the particular optical mode for selectively mapping the target modes. This approach has been widely used to analyze the plasmonic and photonic modes of various systems 7,14,15 .
With advancements in fabrication technology, there have been growing interest in designing smaller nanostructures and integrating them at a greater density for merging multiple functionalities within a smaller chip 7,11 . This often makes it necessary to engineer the hybridized modes having distinct resonances with minimal cross-talks by properly engineering the coupling among the modes of the constituent nanostructures 12,16,17 . However, mode hybridization is more complicated because the modes of the basic building blocks such as nanoparticles, nanorods and nano-slits are spectrally broadened and spatially overlapped with the reduction of their size 18 . In this respect, most of the existing NSOM modalities are not well suited for elucidating the detailed mode formation mechanism. They tend to map either the dominant modes or integration of all the modes at the given excitation/detection wavelength. To better understand the formation of the hybridized modes in the small-scale nanophotonic devices, it will be important to have experimental tools that can separately map the multiple superposed modes.
Here, we proposed a near-field imaging method termed a nearfield transmission matrix microscopy for simultaneously mapping multiple hybridized eigenmodes that are formed by the coupling between the modes of the constituent nanostructures. We constructed an interferometric NSOM system based on the selfinterference geometry and measured both the amplitude and phase maps of the near-field wave at the surface of nanostructures. Various phase imaging NSOM techniques have been developed in the past 14,[19][20][21] , but here we made an important addition. We scanned the angle of the far-field illumination and recorded the near-field complex-field maps for various illumination angles. These measurements allow us to construct a fully phase-referenced far-to near-field transmission matrix (FNTM), which describes the far-field input to near-field output response of the given nanostructures. By performing the singular value decomposition (SVD) of the measured matrix, we could extract individual orthogonal near-field eigenmodes contributing to the input-output response of the device. With this approach, multiple high-order modes of the double nano-slits whose spacing (50 nm) is significantly smaller than the diameter of the NSOM probe aperture (150 nm) were made visible. This is an important addition of information to the conventional NSOM, in which only the integrated signal from all the contributing modes at the excitation wavelength is detected. We provided a theoretical model to support the working principle of the proposed method and conducted numerical mode analysis to validate the experimentally observed modes. In addition, we found that the mapping of the antisymmetric mode allows us to better resolve the existence of multiple nano-slits whose gap size is smaller than the probe aperture.

Results
Experimental measurement of a FNTM. A transmission matrix describes the coherent linear interaction between light and arbitrary device including disordered media 22,23 . It describes the complex-field map, i.e., phase and amplitude maps, at the output plane of the device depending on the excitation of free modes at the input plane. It has been widely used in the far-field regime in the past and offered unique opportunities. For example, the inversion of the transmission matrix led to image delivery through a scattering layer 24 , and the eigenmodes of the transmission matrix were used for efficient light energy delivery and mapping target objects within the scattering medium [24][25][26][27][28] . The transmission matrix approach also allows for the control and image delivery of near-field waves [29][30][31][32] . In the present study, we developed an experimental method to record a fully phasereferenced FNTM, t(x, y; k in ), which describes the near-field complex-field maps at the upper surface (x, y, z = 0) of the nanostructures for the far-field illumination at the bottom surface (z = −z 0 ) with various transverse wavevectors, k in ¼ k x in ; k y in À Á , set by the illumination angle. Figure 1 shows a simplified experimental scheme where an NSOM aperture probe picked up the near-field wave at the upper surface (see also Supplementary Note 1 for the detailed experimental setup). The amplitude and phase of the near-field wave were recorded by using the selfinterference phase-shifting interferometry. To achieve this, we installed a spatial light modulator (SLM) at the conjugate plane to the sample plane in the illumination beam path and wrote a phase pattern on the SLM to generate two plane waves at the bottom of the nano-slits, one with a normal incidence angle (i.e., k in = 0) and the other with nonzero k in Here, the normally incident plane wave, whose wavefronts are indicated by the red lines in Fig. 1a, serves as a reference wave. Its relative phase Δϕ R with respect to the other sample wave, indicated by the blue lines in Fig. 1a, was controlled by the phase pattern on the SLM. The transmitted wave on the upper surface of the nanostructures can be expressed as where E R (x, y) and E S (x, y; k in ) are the complex-field amplitudes of the reference and sample waves, respectively, at the upper surface. The NSOM probe recorded the interference intensity I out ¼ E out x; y; 0; k in ; Δϕ R À Á 2 , which is expressed as a sinusoidal function of Δϕ R (see the example in Supplementary Note 1). By measuring the intensities at four incremental steps of Δϕ R , i.e., Δϕ R ¼ 0; π 2 ; π and 3π 2 , we could demodulate the amplitude and phase of E S (x, y; k in ) 33 . Here, we accounted for the amplitude of E R (x, y) by separately measuring its intensity map taken only for the normal illumination. The phase of E R (x, y) was assumed to be flat in space, which is typical for the symmetrically driven subwavelength nanostructures (see Supplementary Note 3). In the experiment, we chose an illumination wavevector k in and sequentially displayed four phase patterns on the SLM setting Δϕ R to multiples of π 2 . The NSOM probe was scanned across the lateral plane at the upper surface as shown in the left image of Fig. 1a to obtain four corresponding near-field intensity maps. By applying the phase-shifting interferometry algorithm, we obtained the near-field complex-field map for the corresponding k in , i.e., E S (x, y; k in ), and used it to construct the FNTM.
To validate the proposed method, we prepared a pair of nanoslits whose characteristic physical dimensions are significantly smaller than the diameter of the NSOM probe aperture (100 or 150 nm). Nano-slits are one of the fundamental building blocks of complex devices 12,34 . The formation of hybridized eigenmodes by the mode coupling between the nearby nano-slits was well understood both theoretically and experimentally 34,35 . Therefore, devices composed of nano-slits can serve as ideal platforms for validating the proposed method. We fabricated double nano-slits on a 100-nm-thick gold film by the focused ion beam milling, as shown in the upper image in Fig. 1b (see Supplementary Note 4 for the sample fabrication details). The width W and length L of each slit were approximately 20 and 160 nm, respectively. The gap D between the two slits was approximately 50 nm. The long axis of the nano-slits was rotated by θ = 34°clockwise with respect to the vertical line, and the polarization direction of the far-field illumination was set to x-direction (red arrow in Fig. 1b). In this arrangement, both the transverse and longitudinal modes with respect to the long axis of the nano-slits were excited. The aperture diameter α of the NSOM probe was 150 nm, which was much smaller than λ exc = 637 nm, but three times larger than the slit gap. The NSOM probe was made of an Au-Cr coated tapered waveguide, and its aperture dimension is shown in Fig. 1b for the direct comparison with the dimension of the nano-slits. The bending direction of the tapered fiber probe used in the experiment determines the polarization of near-field collection, which is indicated by the green arrow in Fig. 1b 36 . Similar to the excitation geometry, this allows the detection of both the transverse and longitudinal modes. For each k in , we scanned the NSOM fiber probe around the center of the nano-slits over an area of 800 × 400 nm 2 with a scanning step of 25 nm. The distance of the NSOM probe to the sample surface was maintained at approximately 20 nm throughout the set of measurements. We repeated the same measurements for 100 different k in , which uniformly covered the full numerical aperture (NA = 0.6) of the bottom objective lens (see Supplementary Note 1 for the detailed coverage of k in ). The measurements of the entire angular set of complex-field maps lasted around 27 min. Figure 1c shows the representative near-field complex-field maps for each k in in unit of k 0 . These individual NSOM images, each of which corresponds to a conventional NSOM image, mainly visualized the fundamental symmetric mode. Multiple higher order transverse and longitudinal modes were spectrally overlapped to such a degree that individual modes cannot be excited by the selected excitation wavelength 37 . Therefore, relatively weak higher order modes were obscured by the fundamental mode.
Extraction of high-order near-field eigenmodes. The near-field complex-field map (Fig. 1c), E S (x, y; k in ), resulted from the linear superposition of the orthogonal near-field eigenmodes. With the increase of k in j j, the phase difference of the excitation wave at the two nano-slits is increased such that the way the eigenmodes are superposed varies with k in . In this section, the extraction of individual eigenmodes contributing to E S (x, y; k in ) is described. First, we constructed a FNTM, t(x, y; k in ), by assigning 100 complex-field maps taken by different values of k in as constituent columns. It is worth noting that the maintenance of the phase stability is crucial as the measurement was performed on a pointby-point basis. As all of the phase measurements were carried out with respect to the normal illumination whose relative phase was well controlled by the SLM, our measurements were sufficiently robust to link multiple measurements together in their phases. This ensured that the measured FNTM was fully phasereferenced. Figure 2a shows the phase part of the FNTM. To identify near-field modes, we performed the SVD of the matrix, i.e., t x; y; k in ð Þ¼ UΣV y , where V and U are the unitary matrices whose columns are the input and output eigenchannels, respectively, and ∑ is a diagonal matrix whose diagonal elements are nonnegative real numbers referred to as singular values 23 . V y indicates a conjugate transpose of V. The squares of the singular values are the eigenvalues of t y t, which are shown in Fig. 2b after sorting them in the descending order with respect to the eigenchannel index. Essentially, SVD could identify the orthogonal basis at the input (V) that is mapped onto the orthogonal output basis (U) for the linear transformation operator, t. Therefore, the columns of U contain near-field eigenmodes at the upper plane of the nano-slits (see Supplementary Note 2 for further explanation). In addition, the physical meaning of the eigenvalue is the coupling efficiency of each mode from the far-field input to the near-field output. We found that the first six largest eigenvalues were above the noise level for the sample with θ = 34°; their eigenchannels were the near-field eigenmodes of the double nano-slits.
In Fig. 2d, we visualized six near-field eigenmodes obtained from the first six columns of U. The first column of U, associated with the largest eigenvalue, corresponds to the symmetric mode (TE 00 mode). The mode indices n and m in TE nm indicate the orders of transverse and longitudinal modes, respectively, with respect to the long axis of the nano-slits, and the associated eigenvalue is represented as τ nm . The fourth column of U corresponds to the antisymmetric transverse electric mode (TE 10 mode), according to the spatial phase distribution and sharp dark line in the middle. Notably, the eigenvalue of TE 10 mode was 13.4 times smaller than that of the symmetric mode (TE 00 mode) as shown in Fig. 2b. This explains why it was not visible in the individual near-field maps in Fig. 1c. Based on the shape of this antisymmetric mode, we could estimate the position and the rotation angle of the double nano-slits. The long axis of the slit was rotated by 34°with respect to the polarization of the incident wave, which agrees well with our experimental preparation. The two white rectangular boxes in the TE 00 mode map indicate the positions and orientations of the slits identified by the antisymmetric mode. Furthermore, we identified other longitudinal highorder modes, the TE 01 , TE 02 , and TE 03 modes from the second, third, and fifth largest eigenvalues, respectively. These eigenmodes were excited by the polarization component of the illumination along the long axis of the nano-slits. We could even identify the quadrupole mode (TE 11 mode), whose near-field map shows that the directions of the dipole moments at the two slits were opposite. This mode was rarely observed in conventional near-field imaging as the coupling efficiency of the far-field energy to this mode was extremely low, which was reflected by the extremely small eigenvalue (28.4 times smaller than the eigenvalue of TE 00 mode). All these high-order modes revealed the local subaperture-scale near-field waves as well as the structural details of the nanostructures that are invisible in the symmetric mode.
We performed a separate measurement for the sample rotated by an angle θ = 14°, smaller than the previous one. The eigenvalue distribution for this sample geometry is shown in Fig. 2c and its eigenmodes are shown in Supplementary Note 1. We observed that the eigenmodes were rotated according to the rotation of the sample. There were two noteworthy observations that support the validity of our measurements. The eigenvalue ratio τ 01 =τ 00 was reduced from 0.4 to 0.1 as θ was reduced from a t (x,y ; k in )  34°to 14°. This is because the excitation and collection of the longitudinal TE 01 mode were less efficient with the reduction of θ.
On the contrary, the eigenvalue ratio τ 10 /τ 00 was reduced from 0.07 to 0.05, which remained similar with the rotation of the sample as the TE 00 and TE 10 modes are both transverse mode. It is noteworthy that the scanning step of the NSOM probe (25 nm) was much finer than the resolving power set by the probe aperture diameter (150 nm) 38 to ensure high mode reconstruction fidelity and minimize the pixelation artifact. From the separate analysis, we confirmed that the scanning step of 50 nm was good enough to map all the observed modes. The underlying principle of extracting high-order eigenmodes by the SVD of the measured FNTM can be understood by a simple double-slit model, where we considered only the symmetric and antisymmetric modes (see "Methods" and Supplementary Note 2 for the detailed theoretical model). As the slit separation is too small for far-field illumination, i.e., the slit gap D is much smaller than the wavelength, far-field illumination is mainly coupled to the symmetric mode, and the antisymmetric mode is barely excited. For a given far-field incident wavevector k in , the phase difference of the incident wave between the two slits is given by Δφ k in ð Þ ¼ k in Á D j j, where D is a vector connecting the centers of the two slits. As k in j j≤ k 0 , Δφ ≤ 2πD λ $ 0:5, much smaller than π. Therefore, the incident wave is coupled mostly to the symmetric mode even at the maximum incidence angle. This explains the visibility of only the symmetric mode in conventional NSOM imaging (Fig. 1c) and the eigenvalues of the higher-order modes were tens of times smaller than the symmetric mode (Fig. 2b). In the experiment, the linear combination of these orthogonal modes was measured, and SVD served as the means to identify individual orthogonal modes from the superposed measurements. In the simple double-slit model, we constructed a FNTM after accounting for the far-field excitation limit and analytically demonstrated that SVD can identify the symmetric and antisymmetric modes separately. The eigenvalue ratio τ 10 /τ 00 between the antisymmetric and symmetric modes is estimated as 1 16 Δφ 2 in the weak-coupling regime. This model explains the experimentally measured ratio, which was on the order of 10 −1 -10 −2 . A more general model incorporating the spatial shapes of the near-field modes was also developed (see Supplementary Note 2) to confirm the capability of our FNTM approach for extracting high-order near-field eigenmodes.
As the eigenvalues of higher-order modes were extremely small in the far-field excitation, it was crucial to increase the sensitivity of the measurements. This was especially the case given the noisy nature of the near-field recording due to weak signal strength. A large number of k in measurements were important in this instance. Although it is sufficient for the number of k in 's to be equal to the number of orthogonal modes in noise-free measurements, we used 100 different k in values for the matrix measurement, far larger than the required number of measurements. This increase in the number of independent measurements increased the fidelity of mode mapping, thereby enhancing the signal-to-noise ratio, particularly for the mapping of higherorder modes. To confirm this, we constructed multiple FNTMs by varying the number N in of k in measurements used for matrix construction. By performing the SVD of each matrix with N in columns, we identified the observable near-field eigenmodes. As shown in Fig. 3, the FNTM with N in = 2 showed only the TE 00 mode because the eigenvalue of TE 01 was smaller than the experimental noise level. In contrast, the FNTM with N in = 10 columns revealed the TE 01 mode. When we increased N in , various higher-order near-field modes were observed as the eigenvalues of the corresponding eigenmodes were increased beyond the noise level. Observing the TE 02 mode extracted in the N in range of 10-100, we could recognize that the image quality was improved with the increase in N in . This clearly supports that a large number of k in measurements played a crucial role in identifying the nearfield eigenmodes. Conventional NSOM corresponds to N in = 1 case for the plane wave excitation or coherent integral of all the measurements for the focused illumination, which explains why it is unable to extract individual modes.
Numerical validation of the experimentally observed modes. We validated the experimentally observed near-field eigenmodes by numerical mode analysis based on the finite-difference timedomain (FDTD) method. As shown in Fig. 4a, the geometric parameters of the sample were taken from the fabricated sample (L = 160 nm, W = 20 nm, and D = 50 nm). Here, the double nano-slits were prepared on a 100-nm-thick gold film on a SiO 2 substrate. The real and imaginary dielectric constants of gold were taken from Johnson and Christy's experimental study 39 and fitted with the Drude model to conduct FDTD simulations. The refractive index of glass was set to 1.45. To ensure the precision of simulation, we defined the fine meshes around the nano-slits (Δx = 1 nm, Δy = 1 nm and Δz = 5 nm), and the far-field illumination was set normal to the surface of the sample. To extract multiple hidden eigenmodes under normal-incident pumping conditions, we intentionally restricted the field symmetry conditions in FDTD simulation along the xzand yz-planes at the center of the nano-slits (Fig. 4a) and controlled the polarization direction of the excitation source to xor y-direction. Figure 4b shows the maximum |E| 2 spectra of the double nanoslits for three representative symmetry conditions with two different polarization directions. Here, the notation of "(odd, even), x-pol", for example, indicates that the simulation had an odd and even mirror symmetry of electric fields with respect to the yzand xz-plane, respectively, and the excitation source was set to x-polarization. The maximum |E| 2 was detected above the upper surface of the sample over an area of 10 × 10 μm 2 around the nano-slits. Notably, the six prominent individual eigenmodes (TE 00 , TE 01 , TE 02 , TE 10 , TE 03 , and TE 11 ) were observed within a spectral window of 450-800 nm, and their spectra overlapped with respect to one another. This is mainly due to the small physical dimension of the nano-slits, where the Q factor of each mode was extremely low (Q < 10) and the bandwidth of each mode was around 100 nm in wavelength. Consequently, multiple modes could contribute simultaneously to the measured nearfield signals at a given excitation wavelength 37 . It is worth noting that the far-field excitation is favorably coupled to lower-order modes as their spatial field is better overlapped with the far-field excitation than the higher-order modes at the input plane. For this reason, the near-field maps obtained in experiment (Fig. 1c) showed mainly the fundamental symmetric mode with minimal signatures of high-order eigenmodes.
For each eigenmode, we obtained its near-field complex field map, which is shown in Fig. 4c in descending order according to the resonance wavelength for each polarization. Here, the field profiles were taken 20 nm above the upper surface of the sample. The polarization state of each mode is indicated as "x-pol" or "y-pol" in each figure panel. The TE 00 and TE 10 modes show an odd and even mirror symmetry with respect to the yz-plane, respectively, where the x-components of the electric field are dominant. Therefore, they have an even and odd phase symmetry with respect to the yz-plane, respectively. In addition, the TE 01 , TE 02 , and TE 03 correspond to high-order eigenmodes with polarization parallel to the long axis of the nano-slits. The TE 11 mode shows a quadrupole mode profile with polarization orthogonal to the long axis of the nano-slits. Plasmonic modes are often subtle to explain due to their short lifetime and low Q factor 40 . In our study, these observed eigenmodes are related to the antenna modes of the double nano-slits whose charge distributions are indicated in Fig. 4c. We found that these mode profiles are in excellent agreement with those obtained in the experiment shown in Fig. 2c. This result confirmed the validity of our experimental eigenmode mapping method. (odd, even), x-pol (even, even), y -pol (even, odd), y -pol (even, even), x-pol (even, odd), x-pol = 541 nm x -pol (even, odd) = 717 nm y-pol (even, even) = 646 nm y-pol (even, odd) = 590 nm y-pol (even, even) Fig. 4 Numerical mode analysis of the double nano-slits. a Schematic sample geometry along with far-field excitation geometry. The polarization direction of the excitation source is selected to either xor y-direction. Geometric parameters were taken from the experimentally fabricated double nano-slits (L = 160 nm, W = 20 nm, and D = 50 nm). b Maximum |E| 2 spectra measured above the upper surface of the sample under mirror symmetry conditions of electric field along the xzand yz-plane at the center of the nano-slits. Plots were taken for three representative symmetry conditions with two different polarization directions. Here, the notation of "(odd, even), x-pol", for example, indicates that the simulation had odd and even mirror symmetry of electric field with respect to the yzand xz-plane, respectively, and the excitation source was set to x-polarization. c Spatial electric field (x or y) maps of various near-field eigenmodes TE nm , obtained 20 nm above the gold surface, identified from Fig. 4b. The maps are displayed in descending order according to the resonance wavelength for each polarization. Plus and minus signs indicate the charge distributions of the antenna mode.

Discussion
In summary, we developed a near-field transmission matrix microscopy that can extract the eigenmodes of subwavelength nanostructures whose spectra are too overlapped to selectively excite individual modes by the choice of an excitation wavelength. Under this condition, conventional NSOM can visualize only the dominant modes as it detects the integral of signals from all the contributing modes. In our study, we exploited a new degree of freedom, which is the angle of far-field excitation, to disentangle multiple superposed modes. We experimentally measured the complex-field maps of near-field waves for various angles of farfield excitation and constructed a fully phase-referenced FNTM. This is in contrast with the conventional NSOM, where only a single spatial mode of far-field excitation (normal or a focused illumination) is used. By performing the SVD of the measured matrix, we extracted individual orthogonal eigenmodes. We could visualize the antisymmetric mode, quadruple mode, and other high-order modes of the double nano-slits, which were completely hidden under the conventional NSOM images, and quantify their relative coupling efficiency from their eigenvalues. An analytic model for explaining the working mechanism of eigenmode extraction was presented, and numerical mode analysis confirmed the experimentally observed modes.
The proposed approach is an important advancement in the context of near-field imaging. It provides a new dimension of information that is inaccessible to conventional NSOM, which is the decomposition of superposed modes. In fact, apertureless NSOM can also map the near-field modes of nanostructures 14,41,42 . Cathodoluminescence 43,44 , photoemission electron microscopy 45,46 , and electron-energy-loss spectroscopy 47,48 that are based on electron microscopy are also useful tools to map local optical modes with high spatial resolution. These methods could often visualize dark modes having no net dipole moments 41,49 as well as bright modes. However, these other methods usually map the dominant modes at the given excitation/detection wavelength. The reason these methods could map the dark modes is because the spectral bandwidths of the dark modes are well separate from those of bright modes. On the contrary, our method can separate multiple spectrally and spatially superposed hybridized bright/dark modes by means of the SVD. This capability can be especially useful in interrogating small-scale nanostructures where the modes of basic building blocks are spectrally so broad that the hybridized modes can have spectral overlaps. The double nano-slits used in our study is one of the good examples.
Our method can potentially be combined with other existing NSOM modalities 14,19,[50][51][52][53] , particularly those based on interferometric detection, and help extract near-field eigenmodes information of various quantities such as electric near-field vector components, magnetic near-field, and time/frequency-resolved measurements 19,54 . The ability to extract the spectrally and spatially overlapped eigenmodes can help designing the nanophotonic devices containing multiple hybridized modes 55,56 and providing new insights into the development of novel nanoscale photonic devices 11,57 . In the present study, the mode decomposition experiment was conducted only for single wavelength excitation. Repeating the same measurements with the scanning of the excitation wavelength could allow the full mapping of the spectra of individual resonance modes constituting the total nearfield scattering spectrum. Another noteworthy observation in this study was that higher-order eigenmodes exhibited multiple subaperture nodes due to the destructive interference of local nearfield waves. This opened a new possibility to resolve the fine structural details of nano-slits whose gap is smaller than the physical size of the probe aperture 58,59 (see Supplementary Note 1). Although this type of resolution enhancement is applicable in limited cases where there are distinct resonance modes, it is still valuable considering that it can retrieve deep subwavelength structural information hidden with conventional NSOM. Considering that the steep reduction in near-field collection efficiency accompanied by the use of the smaller probe sets the practical limit of the spatial resolving power, the capability of imaging with the same resolution by the use of a larger probe aperture will push the ultimate limit of the resolving power in the near-field imaging.

Methods
Experimental setup. To perform the experimental mapping of near-field modes, we integrated a far-field phase modulation system into an NSOM system (Nanonics MV2000) (see Supplementary Note 1 for the detailed experimental setup). The output beam of a laser diode (Thorlabs Inc., LP637-SF70) was enlarged and collimated to uniformly illuminate the SLM (Hamamatsu LCOS-SLM X10468). To measure both the amplitude and phase of the near-field waves on the surface of the nano-slits with their phases fully referenced, we employed a self-interference measurement system without an additional setup for the reference beam line. For this purpose, the SLM generated both the sample and reference beams whose relative phase was controlled in multiples of π/2. They were demagnified and delivered to the bottom of the NSOM sample stage using an objective lens (Nikon ELWD 40×, NA: 0.6). The overall magnification from the SLM plane to the sample stage was 1/1000×. Near-field waves generated on the upper surface of the nanostructure were measured by scanning the NSOM fiber probe. The aperture diameter of the NSOM probe was 150 or 100 nm. The near-field light captured by the NSOM probe was delivered to a photomultiplier tube (Hamamatsu, H8259-01).
Simple double-slit model. The identification of the antisymmetric mode by SVD of the measured FNTM can be understood by a simple double-slit model. The transmission of the incident field through a simple double-slit system can be described by where, the vectors E 1 in ; E 2 in À Á and E 1 out ; E 2 out À Á represent the electric fields at the input and output planes, respectively (the subscripts 1 and 2 indicate the left-and righthand slits, respectively), J is the coupling constant to the output of the same slit as the input, and K is the coupling to the other slit's output. In the experiment, we cannot directly measure J and K by far-field excitation because the slit separation is too small for far-field illumination to couple light to individual slits. Likewise, the symmetric and antisymmetric modes cannot be individually addressed by the farfield excitation. Instead, the combination of orthogonal modes was measured in the experiment, and SVD was performed to identify the orthogonal modes from the superposed measurements. For the given far-field incident wavevector k in , the phase difference of the incident wave between the two slits is Δϕ(k in ), as defined in the main text. The incident electric field at the two slits can be expressed as E 1 in ; E 2 in À Á ¼ E 0 ; E 0 e ÀiΔϕðk in Þ À Á for any given k in . We sent 100 different incident wavevectors in the experiment; however, for simplicity, we consider two representative incident wavevectors, one with normal illumination and the other with k in . The incident electric fields at the two slits are respectively written as (E 0 , E 0 ) and E 0 ; E 0 e ÀiΔϕ k in ð Þ À Á . Their corresponding output electric fields, which correspond to the quantities measured by the interferometric NSOM in the experiment, can be obtained by inserting these vectors into Eq. (3). Using these two output fields, we can construct a FNTM and analytically calculate the eigenvalues and eigenmodes of t y t (see Supplementary Note 2 for details). We could verify that the two output eigenvectors are v S % 1 ffiffi  Fabrication of nano-slits. To fabricate the nano-slits with sub-50-nm spacing, we used proximal milling techniques in Ga + -based focused ion-beam processes on a sputtered gold film prepared on a silica coverslip (see Supplementary Note 4 for the fabrication details). We intentionally off-designed the milling patterns from the original nano-slit design to make use of the proximity effect of focused ion beam milling. By optimally controlling the distance between rectangular milling patterns and the total milling time, we could realize sub-50-nm spacing between the nanoslits.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.