Spatio-spectral metrics in electron energy loss spectroscopy as a tool to resolve nearly degenerate plasmon modes in dimer plasmonic antennas

2.1 Plasmon modes of disc-shaped dimer

To facilitate the discussion of the spatial distribution of the loss probability, we introduce a cartesian coordinate system with the axes x and y parallel with the transverse and longitudinal directions of the disc-shaped dimer, respectively, and the origin located in the middle of the gap between both discs. Further, we utilize two polar coordinate systems originating in the centers of the discs, with the zero angle corresponding to the center of the dimer, and increasing angle toward the right side of the discs (clockwise for the bottom disc and counter-clockwise for the top disc). All coordinates are represented in Figure 1(a and b).

Figure 1(c–e) shows intuitive schemes of the LDB, TD, and LDA modes. The modes are represented by two oscillating dipoles, one in each disc, oriented along the longitudinal (y) or transverse (x) direction. Current oscillates along the arrows, and charge accumulates near the boundaries of discs in the areas marked with + and − signs. According to Gauss law, this charge is a source of the electric field, which acts at the electron beam and results in the energy loss of the probing electrons. Consequently, the loss probability is large for the electron beam passing through the areas of accumulated charge. The only exception is the gap region of the LDB mode, where the opposite charges accumulated at the opposite sides of the gap between the discs effectively cancel each other due to the long-range character of Coulomb interaction with the probing electron, and the loss probability exhibits low values there. The regions of the large loss probability are displayed by red patches in Figure 1(c–e). These regions are characteristic for each of the LDB, TD, and LDA modes (but not for the pair of TDB and TDA modes). The LDB mode features two loss-probability maxima along the longitudinal direction at the outer sides of the discs, the TD modes feature four maxima at the transverse sides of the discs, and the LDA mode exhibits three maxima along the longitudinal direction (two at the outer sides of the discs and one in the gap region between the inner sides).

Loss probability Γ(x, y, E), a three-dimensional function of the electron beam position x, y and the loss energy E, is typically visualized in the spectral domain as a spectrum for a specific beam position x ₀, y ₀,

or in the spatial domain as a map for a specific loss energy E ₀,

The loss probability spectra obtained for our disc-shaped dimer and three distinct beam positions are shown in Figure 2. To suppress the noise, the experimental spectra are integrated over multiple (several tens) pixels. The maps for several selected energies (later identified as the energies of LDB, TD, and LDA modes) are shown in Figure 3. To suppress the noise, the experimental maps are integrated over multiple energy slices in the range 0.1 eV around the central energy (corresponding to 10 energy slices with our energy step of 0.01 eV).

Figure 2:

Loss probability spectra of the disc-shaped dimer PA obtained from the numerical simulation (a) and the experiment (b). The insets show the electron beam positions for which the spectra were recorded. The color of the spot representing the electron beam position is identical to the color of the corresponding spectrum: pink for the outer longitudinal sides, orange for the transverse sides, and green for the gap.

Figure 3:

Loss probability maps of the disc-shaped dimer PA obtained from the numerical simulation (top) and the experiment (bottom) for several loss energies.

The theoretical loss spectra [Figure 2(a)] exhibit three well separate peaks at the energy of 1.0 eV, 1.2 eV, and 1.4 eV. The peaks can be assigned to specific modes by a simple inspection of spatial maps. The peak at 1.0 eV with two longitudinal maxima [Figure 3(c)] corresponds to the LDB mode, the peak at 1.2 eV with transverse maxima [Figure 3(e)] corresponds to the TD modes, and the peak at 1.4 eV with two longitudinal maxima and one central maximum [Figure 3(g)] corresponds to the LDA mode. The experimental spectrum [Figure 2(b)] exhibits three quite resolvable features at the energy of 0.8 eV, 1.0 eV, and 1.2 eV. Based on the examination of individual plasmonic discs (see Supplementary Material, Section S1), we can attribute the peaks to hybridized dipole modes since the quadrupole and higher modes are well outside this spectral region. However, the related spatial maps (Figure 3, bottom panels) do not clearly show features characteristic of considered modes and cannot be utilized for unambiguous peak-to-mode assignment. A natural question then arises whether any information on the modes is present and can be extracted from the experimental loss probability Γ(x, y, E). We will show that the answer is positive when inspecting the loss probability in both spatial and spectral domains simultaneously.

2.2 Spatio-spectral metrics

One of the characteristic features of the considered modes is the spatial spread of the loss probability in the x and y (transverse and longitudinal, respectively) directions. This feature can be quantified by the so-called absolute central moments, similar to standard statistical moments, defined as follows. We first introduce the loss image weight w(x, y, E), a quantity proportional to the loss probability but normalized to unity when integrated over a certain area of interest: ∫w(x, y, E) dx dy = 1. Clearly,

where the integration goes over the area of interest. For simulations, the area of interest is a full simulation area. For experimental data, though, we selected two rings concentric with the plasmonic discs with a width of 100 nm where the signal-to-noise ratio is high [see Figure 4(b and c)]. Next, we define the central coordinate of the loss image as

and the (energy-dependent) absolute central moments are defined as

Figure 4:

Area of interest for experimental values of the spatio-spectral metrics. (a) ADF STEM image of the disc-shaped dimer PA, with white/dark color representing gold/substrate. The red circles determine the boundaries of the discs utilized in the processing of the experimental loss probability. (b) A map of the loss probability at the energy of 1.0 eV. The black circles determine the boundaries of the discs and the gray circles mark the area with the high signal-to-noise ratio utilized for the calculation of the spatio-spectral metrics. (c) Same as panel b showing only the area with a high signal-to-noise ratio.

Before inspecting the actual values of the moments, we shall qualitatively discuss their expected properties. The longitudinal spread of the loss function shall be larger for the LDB mode than for the LDA mode due to the presence of a strong central maximum for the latter with near zero longitudinal (i.e., y) coordinate. Therefore, the longitudinal absolute central moment M _y (E) shall be larger at the energy of the LDB mode compared to that at the energy of the LDA mode. Similarly, the longitudinal spread and the moment M _y (E) of the TD modes shall be smaller than that of the LDB mode since the transverse maxima of the former are related with a lower longitudinal coordinate than the longitudinal maxima of the latter. We note that such a qualitative comparison is not possible for the TD and LDA modes. In the transverse direction, the TD mode shall feature a rather large transverse spread, accompanied by the large values of the transverse central moment M _x (E), presumably larger than for the LDB and LDA modes.

The spectral dependence of moments M _x (E), M _y (E) (normalized so that the maximum value in the inspected spectral region equals to unity) is shown in Figure 5. The theoretical spectra are shown in the left panel. We will first discuss the longitudinal moment M _y (E). As expected, its large value at the energy 1.0 eV corresponding to the LDB mode gradually decreases over the energy 1.2 eV corresponding to the TD modes to the energy 1.4 eV corresponding to the LDA mode. Strikingly, the maximum and minimum of the longitudinal moment corresponds exactly to the energies of the LDB and LDA modes, respectively. This suggests that the energies of the maximum and minimum can be used as the operational definition for the LDB and LDA modes energy if no more straightforward measure is available (such as in the case of the experiment). Similarly, the maximum of the transverse central moment at the energy of 1.2 eV agrees with the energy of the TD modes determined from the loss spectra (cmp. Figure 2 and related discussion) and can be used as the operational definition for the TD mode energy.

Figure 5:

Spectral dependence of the longitudinal absolute central moment M _y (E) (blue) and the transverse absolute central moment M _x (E) (red). The insets show the schemes of the LDB, TD, and LDA modes, and their position on the energy scale qualitatively corresponds to the actual mode energies. The colors of the patches representing the loss probability maxima correspond to their significance for the longitudinal (blue) and transverse (red) spectral features of the moments.

The experimental spectra of the moments correspond very well to the theoretical spectra, only with the reduced variance (3–5 times) over the examined spectral range, which is attributed to the experimental noise and/or background. This seemingly trivial statement represents an important finding. In contrast to theoretical maps of the loss probability, the experimental maps do not exhibit the apparent formation of the modes-identifying features such as the longitudinal, transverse, and central maxima. Still, the information allowing to identify the modes and to determine their energy can be extracted with judiciously defined metrics. The experimental energies of the LDB, TD, and LDA modes read 0.9 eV, 1.0 eV, and 1.2 eV, respectively (with a precision of 0.05 eV set by the half-width of the energy window used for the spectral integration of the experimental data), and correspond reasonably well to the positions of the peaks in Figure 2(b) at 0.8 eV, 1.1 eV, and 1.3 eV. We suppose that the energies based on the absolute central moments are more relevant, since they take into account the spatial distribution of the loss probability, while the energies of the peaks in Figure 2 are more prone to the experimental noise and possible incomplete background subtraction.

The spatial distribution of the loss probability can be inspected also in polar coordinates. A nontrivial question arises how to set the origin of the polar coordinate system. The examined modes of the dimer are formed through the hybridization of the modes of individual discs. Therefore, it is meaningful to inspect the mode of a specific disc in its local coordinate system related to that disc. We use the following approach. We utilize the same area of interest as in the case of absolute central moments [i.e., two rings concentric with the plasmonic discs with a width of 100 nm; see Figure 4(b and c)]. We also obtain the loss image weight as the loss intensity normalized to the unit integral over the area of interest. The loss image weight in each ring of the area of interest is processed separately in a local coordinate system of the relevant disc shown in Figure 1(b) (the pixels from the intersection of the rings are processed twice). The loss image weight w(x, y, E) is transformed to the (local) polar coordinates r, ϕ,

and then integrated over the radial coordinate, leaving only angular dependence,

The indices U and L refer to the rings around the upper and lower disc, respectively, and the related local coordinate system. Finally, we define the zeros and positive directions of the angular coordinates in a way that respects the symmetry of the dimer. In this way, the angles 0, π/2, π, and 3π/2 correspond for both local coordinate systems to the center, right-hand (transverse) side, outer (longitudinal) side, and left-hand (transverse) side [see Figure 1(b) and a related discussion]. This also allows us to inspect the total angular distributions of the loss image weight

instead of individual components. For brevity, we will refer to the quantity W(ϕ, E) as angular loss weight in the following. Now, the W(ϕ, E) shall feature maxima at the angle π for the LDB mode, the angles π/2 and 3π/2 for the TD modes, and the angles 0, π, and 2π for the LDA mode.

The angular loss weight calculated and experimentally retrieved for the dimer under study is shown in Figure 6(a and b). The calculated angular loss weight evolves with the energy as expected, featuring a single peak at the angle of π corresponding to the LDB mode at 0.9–1.0 eV, two peaks at the angles of π/2 and 3π/2 corresponding to the TD modes at 1.2 eV, and two peaks at the angles 0 (and also 2π) and π corresponding to the LDA mode at 1.4–1.5 eV. The energies correspond well to the energies determined from the spectra (Figure 2) and absolute central moments (Figure 5). The experimental weight angular distribution is qualitatively similar, with the signatures of the LDB, TD, and LDA modes at the energies of 0.8 eV, 1.0–1.1 eV, and 1.2–1.3 eV, respectively. However, the distribution profiles are distorted. Most prominently, they feature a strong peak around 4 rad for all energies, originating in the irregular shape of the dimer [a protrusion at the left upper corner of the dimer with a strong loss intensity observable in Figure 2(a–h)]. To compensate for such irregularities, we proposed a renormalized total angular distribution of the loss weight.

i.e., we divide the raw weight angular distribution by the sum of the weight angular distributions over a range of energies. In this way, the renormalized weight W ̃ ( ϕ , E ) has a sum over the energy independent of the angle, and all angular variations are related purely to the properties of a specific mode. The renormalized total angular distribution of the loss weight is shown in Figure 6(c and d). For calculated data, we observe only insignificant differences between raw and renormalized weights, supporting the feasibility of the renormalization. For experimental data, the renormalized distributions are more symmetric and have peaks closer to the expected angles. The energies of the modes estimated from the renormalized weights are identical to those estimated from raw weights.

Figure 6:

Angular loss weights. (a and b) Angular loss weight obtained from the calculated (a) and experimental (b) data, taken at the energy indicated by the numbers (in eV units). The solid lines correspond to the energies of 0.8 eV, 1.0 eV, 1.2 eV, and 1.4 eV. The dashed lines correspond to the energies of 0.9 eV, 1.1 eV, 1.3 eV, and 1.5 eV. (c and d) Renormalized angular loss weight obtained from the calculated (c) and experimental (d) data.

The spatio-spectral metrics introduced in the manuscript predict the mode energies in agreement with each other as well as in agreement with EEL spectra, as summarized in Table 1. The metrics offer several benefits over the traditional EEL spectra. First, they are tailored to specific modes and describe their features in more detail than the spectra. Second, they utilize a larger part of the EELS data cube. In consequence, spatio-spectral metrics allow reliable identification of the nearly degenerate modes and a determination of their energy. With EEL spectra, this is possible using calculated data but not using noisy experimental data. Both metrics, the absolute central moment and the angular weight distribution, have their advantages and disadvantages. The absolute central moment is well suited for systems with naturally defined longitudinal and transverse directions, such as dimers or single rods. However, its predictive power will be limited for systems lacking such directions, e.g., triangular prisms or single spheres. On the other hand, the angular weight distribution is more general (as long as the natural central point of the system exists).

It is worth comparing the performance of the spatio-spectral metrics with the methods of statistical data analysis (e.g., principal component analysis or non-negative matrix factorization) that are frequently used for processing EELS data including resolving LSP modes [22, 23]. The metrics benefit from the a-priori knowledge of the modes, they can be tailored to a specific problem, and in this respect, they can overperform the unbiased statistical data analysis. This is indeed the case of our data, where the non-negative matrix factorization rather successfully resolves the three modes but fails to identify them, as discussed in Supplementary Material, Section S2. On the other hand, statistical data analysis might perform better when the modes are a-priori unknown or not compatible with the metrics. For illustration, let us consider plasmonic aperture heptamers discussed in Ref. [24]. By inspecting the calculated near fields of the modes supported by heptamers, we infer that neither of our metrics (ACS and ALW) is characteristic for the modes and cannot provide their resolution and identification. On the other hand, it might be possible to use the concept of the spatio-spectral metrics and, considering the expected mode properties, design the metrics tailored to the heptamer.

We foresee several classes of plasmonic structures for which the spatio-spectral metrics might be beneficial: (i) in plasmonic dimers with a weak coupling, they can help to resolve the bonding and antibonding modes. The antibonding modes differ from the bonding ones by the presence of the central maximum in the loss probability. In addition to disc-shaped dimers, such structures include, e.g., bowtie antennas [11], dimers of squares [25] and nanocrosses [26], and heterodimers [27]. (ii) The same task of the mode resolving and identification can be performed in highly symmetrical plasmonic monomers as well as in the case of accidental mode quasi-degeneracy, as encountered, e.g., for the longitudinal quadrupole and transverse dipole modes of a plasmonic nanorod studied in [28]. (iii) In highly symmetric plasmonic structures, e.g., discs and squares, the metrics can be employed to detect minor asymmetry (e.g., elongation along one of the axes).

In addition to LSP modes, spatio-spectral metrics are applicable to other local excitations with nontrivial spectral and spatial dependence, such as the optical modes in dielectric nanoparticles [29] and localized vibrational modes supported by nanoparticles [30, 31]. With some adaptation, the metrics can be applied also to data with additional dimensions, as in the case of spectral, spatial, and angular dependence of localized surface plasmon modes characterized by cathodoluminescence [32, 33].

Finally, we comment on the impossibility to resolve the TDB and TDA modes. First, the energy separation of these modes is very small due to weak interaction between the charges induced at the sides of the discs [Figure 1(d)], much smaller than the spectral resolution of EELS for our system. On the contrary, the energy separation of the LDB and LDA modes is dominated by considerably stronger interaction of the charges across the gap between the discs [Figure 1(c–e)] and is, therefore, much larger. Second, the near fields of the TDB and TDA modes differ only in phase (or sign). Common EELS is not sensitive to the phase of the near field and does not allow for resolving TDB and TDA modes. This issue might be overcome with further development of EELS and introduction of shaped electron beams [34]. We also remark that the far-field optical microscopy is sensitive to the TDA mode only, while the TDB mode with a zero net electric dipole moment does not interact with a plane wave.

3.1 Fabrication of plasmonic antenna

A gold disc-shaped dimer PA was fabricated on the substrate consisting of a silicon nitride membrane with the thickness of 30 nm and the lateral dimensions of 250 × 250 μm². First, a gold layer with the thickness of 30 nm has been deposited by magnetron sputtering. The employed growth protocol produces a large-grain polycrystalline gold layer with optical properties comparable to monocrystalline gold in terms of Q factors of localized plasmon resonances [28]. The dimer was fabricated by focused ion beam milling (using Ga⁺ ions at 30 keV) of the gold film in a dual beam system FEI Helios. The diameter of individual discs has been set to 275 nm and the edge-to-edge distance of the discs in the dimer has been set to 30 nm. The distance of the dimers from the boundary of the metal-free square has been at least 500 nm, which is a sufficient separation to prevent the interaction between the dimers or between the dimer and the surrounding metallic frame [13]. Annular dark-field STEM image of the dimer is shown in Figure 4(a).

Prior to the EELS measurements, the sample was cleaned in an argon-oxygen plasma cleaner for 20 s to prevent carbon contamination [35].

3.2 Electron energy loss spectroscopy

EELS measurements were carried out in a scanning transmission electron microscope FEI Titan equipped with a GIF Quantum spectrometer. The microscope was operated at 120 kV in the scanning monochromated mode with the convergence semi-angle set to 10 mrad and the collection semi-angle set to 11.4 mrad. These settings are optimized for EELS characterization of PAs [36]. The dispersion of the spectrometer was set to 0.01 eV per channel and the full-width at half-maximum of the zero-loss peak was found in the range from 0.10 to 0.12 eV. The probe current was adjusted to around 100 pA. The acquisition time of every spectrum was set to 0.5 ms to use the full intensity range of the CCD camera in the spectrometer and avoid its overexposure. The spatial resolution of the EELS data cube is determined by the pixel size, which was set to 5 nm. Such settings allowed the acquisition of one EELS data cube with a stable electron beam in a reasonable time.

The raw data containing electron counts recorded by the CCD shall be divided by the total number of electrons impinging the sample to obtain the loss probability, taking into account also the sensitivity of the detector. Since this approach is impractical, we utilize instead the electron counts of the zero-loss peak (with the energy integration window from −1 eV to +1 eV) to represent the total number of electrons. The division is then performed pixel-wise, i.e., independently for each position of the electron beam. In this way, we obtain the quantity proportional to the loss probability per channel, which is further divided by the energy interval of the channel (0.01 eV) to obtain the usual loss probability per electronvolt.

Further processing depends on the utilization of the EELS data cube. For EEL maps [Figure 3] and metrics [Figures 5 and 6], we integrate the data over the energy window of 0.1 eV around the target energy to suppress the noise. The zero-loss peak and background are not subtracted. They are both assumed to be constant over the mapped area, the zero-loss peak due to the normalization, and the background due to the homogeneity of the membrane. The EEL spectra [Figure 2] are integrated over tens of pixels to reduce the noise. Next, experimentally determined background (including zero-loss peak) is subtracted to obtain a pure contribution of LSP. The background is measured for the electron beam impinging a bare membrane far from any plasmonic antennas.

3.3 Electromagnetic simulations

Electron energy loss spectra have been calculated with the boundary element method using a software package MNPBEM [37, 38]. In all simulations, the dimers are represented by two gold discs of the height of 30 nm on top of a 30-nm-thick silicon nitride membrane. The diameters of the discs (275 nm) and their edge-to-edge distance (30 nm) were matching the values set in the fabrication process. The dielectric function of gold was taken from Ref. [39], and the dielectric constant of the silicon nitride membrane was set equal to 4, which is a standard approximation in the considered spectral region [12].