Full Three-Dimensonal Reconstruction of the Dyadic Green Tensor from Electron Energy Loss Spectroscopy of Plasmonic Nanoparticles.

Electron energy loss spectroscopy (EELS) has emerged as a powerful tool for the investigation of plasmonic nanoparticles, but the interpretation of EELS results in terms of optical quantities, such as the photonic local density of states, remains challenging. Recent work has demonstrated that, under restrictive assumptions, including the applicability of the quasistatic approximation and a plasmonic response governed by a single mode, one can rephrase EELS as a tomography scheme for the reconstruction of plasmonic eigenmodes. In this paper we lift these restrictions by formulating EELS as an inverse problem and show that the complete dyadic Green tensor can be reconstructed for plasmonic particles of arbitrary shape. The key steps underlying our approach are a generic singular value decomposition of the dyadic Green tensor and a compressed sensing optimization for the determination of the expansion coefficients. We demonstrate the applicability of our scheme for prototypical nanorod, bowtie, and cube geometries.

Electron energy loss spectroscopy (EELS) is a powerful tool for the investigation of plasmonic nanoparticles.[1,2] EELS is a technique based on electron microscopy and measures the probability of a swift electron to lose part of its kinetic energy through plasmon excitation as a function of electron beam position. Following first proof of principle experiments,[3,4] in the last couple of years EELS has been exhaustively used for the investigation of plasmon modes in single and coupled nanoparticles.

Despite its success, the interpretation of EELS data in terms of optical quantities, such as the photonic local density of states[5] (LDOS), remains challenging.[6,7] To overcome this problem, in ref (8) we formulated EELS as a tomography scheme[9] and showed that under certain assumptions a collection of EELS maps can be used to reconstruct the three-dimensional mode profile of plasmonic nanoparticles. A similar approach was presented independently by Nicoletti and co-workers,[10] who demonstrated the applicability of the scheme for a silver nanocube. Extracting three-dimensional information through sample tilting was also shown for a split-ring resonator[11] and a nanocrescent using cathodoluminescence imaging.[12]

The problem with EELS tomography is that the measurement signal (the loss probability) is not simply the integral of local losses along the electron trajectory but involves a two-step process where the swift electron first excites a particle plasmon and then performs work against the induced particle plasmon field. This leads to a nonlocal response function, which allows for a tomographic reconstruction only under restrictive assumptions, such as the applicability of the quasistatic approximation or a plasmonic response governed by a single mode. In this paper we use additional preknowledge, namely, that the particle plasmon fields are solutions of Maxwell’s equations and that the dyadic Green tensor[5] can be decomposed into modes, in order to rephrase EELS in terms of an inverse problem. We develop a rather generic model for the EELS probabilities, which depends on a few parameters, and determine the parameters such that the model data match as closely as possible the measured data. Within this approach we are able to obtain the most accurate reconstructions of the dyadic Green tensor, which, in turn, allows us to extract the three-dimensional photonic LDOS from a collection of tilted EELS maps. We demonstrate the applicability of our scheme for prototypical nanorod, bowtie, and cube geometries.

Theory

We start by analyzing EELS within a semiclassical framework,[1] where a swift electron propagating with velocity loses a tiny part of its kinetic energy by performing work against the electric field E[re(t)] produced by itself. For sufficiently large velocities, we can ignore velocity changes in the electron trajectory re(t) ≈ R0 + vt, with R0 being the impact parameter. It is convenient to split E = Ebulk + Esurf into a bulk contribution[13]Ebulk, corresponding to the electric field within an unbounded homogeneous medium, and a surface contribution Esurf, corresponding to field modifications (including surface plasmons) from the interfaces between different materials. Bulk losses are due to Cherenkov radiation and electronic excitations,[1] and the loss probability is obtained by simply multiplying the loss probability per unit length γbulk(ω), inside material j and for loss energy ℏω, with the path length of the electron inside material j,Bulk losses can be interpreted in terms of local scatterings where the electron emits a photon or excites electrons in the dielectric material and loses part of its kinetic energies. To compute the surface loss probability, we integrate the work dW = eEsurf·vdt performed by the electron over the entire trajectory and decompose it into the different loss energies ℏω according toThus, the energy loss probability becomes[1]where we have explicitly indicated the dependence on the electron propagation direction and the impact parameter through R = (v̂,R0). To understand the physical process underlying eq , it is convenient to introduce the current distribution J(r,t) = −evδ(r – re(t)) of the swift electron and the dyadic Green tensor[5]G(r,r′,ω) that relates for a given frequency ω a current source at position r′ to an electric field at position r via E(r,ω) = iωμ0G(r,r′,ω)·J(r′,ω). The loss probability of eq can then be rewritten in the formwhere dr denotes integration over the spatial variable r. Contrary to eq , the above expression describes a genuinely nonlocal self-interaction process where the electron first induces a field (through excitation of a surface plasmon) and then performs work against the induced field.

In ref (6), the authors tried to interpret eq in terms of the photonic local density of states[5] (LDOS)which is of paramount importance in the field of nanooptics and describes how the decay rate of a quantum emitter located at position r and with dipole moment oriented along n̂ becomes modified in the presence of a structured dielectric environment. While such interpretation can be formally established for nanostructures with translational symmetry along one spatial dimension, it becomes problematic for nanoparticles with generic shape.[7]

A different interpretation of eq in terms of a tomography scheme was formulated independently in refs (8 and 10). As a preliminary step, let us consider the bulk losses of eq for a given R value. Then, each point r inside a medium j contributes with γbulk to the total loss rate. Within the field of tomography[9] it is well-known that the three-dimensional profile of γbulk(r) can be uniquely reconstructed from a sinogram, where bulk losses are recorded for all possible propagation directions v̂, using the inverse Radon transform. Such tomography reconstruction is significantly more complicated for the surface losses of eq since Γsurf is not the sum of local losses (as in the bulk case) but governed by the self-interaction process of excitation and back-action. Only for certain, rather restrictive simplifications, a viable tomography scheme can be formulated:[8,10] the nanoparticles must be small enough such that the quasistatic approximation can be employed; the plasmonic response must be governed by a single plasmonic eigenmode; the sinogram must only consist of electron trajectories that do not penetrate the particle; the sign of the eigenmode potentials must be unique. Although it has been demonstrated that reconstruction is possible in certain cases,[8,10] it is obvious that the above restrictions provide a serious bottleneck for general plasmon field tomography.

In this paper we formulate a significantly more general scheme, which approaches the reconstruction as an inverse problem rather than a tomography scheme. We first describe our approach and discuss possible problems and generalizations at the end. First, we decompose the dyadic Green tensor into a number of modes E(r,ω)where C controls how much the different modes contribute to the decomposition. In the following we only consider positions r and r′ outside the plasmonic nanoparticle and assume that E(r,ω) is a solution of Maxwell’s equations. The expansion of eq is generally possible because G is a symmetric matrix that can be submitted to a singular value decomposition, with C being the singular values and E being the orthogonal matrices. In this respect, eq is similar to a wave function expansion in quantum mechanics into a complete set of basis functions.

To be useful as a reconstruction scheme the modes E(r,ω) should be sufficiently well adapted to the problem such that a limited number n suffices for a suitable representation of G(r,r′,ω). Possible modes are quasi normal modes of the plasmonic nanoparticles,[14−17] which have recently received considerable interest, or natural oscillation modes of our boundary element method approach (see Methods). With these modes, the surface losses of eq becomewhere A±(R,ω) = ∫–∞∞e±v̂·E(R0 + v̂z,ω)dz is the averaged mode profile along the electron propagation direction. We can now formulate our inverse problem as follows. Suppose that one has measured EELS spectra Γexp for a given loss energy and for various impact parameters and electron propagation directions. We then determine the coefficients C such that the entity of measurement data differs as little as possible from the model data of eq ,resulting in a least-squares optimization (we adopt the norm definitions ∥x∥2 = ∑|x|2 and ∥x∥ = ∑|x|). Alternatively, in this work we will use a compressed sensing optimization[18,19]which attempts to minimize the moduli of the expansion coefficients, therefore the scheme is often referred to as a L1-optimization, and μ is a parameter that allows to switch between genuine compressed sensing and least-squares optimizations.[19] For a sufficiently small number of expansion modes E, the determination of the expansion coefficients C is a highly overdetermined problem since the measured loss data can be assembled for many propagation directions and impact parameters R. The only preknowledge entering our optimization is the self-interaction-type scattering process of the electron loss, eq , and the assumption that the dynamics of the electric fields outside the plasmonic nanoparticles is governed by Maxwell’s equations. Importantly, once the coefficients C are determined, we have (approximately) reconstructed the dyadic Green tensor of eq , which allows us to compute all electrodynamic properties including the photonic LDOS.

Results

To prove the applicability of our reconstruction scheme, we generate the “experimental” EELS data Γexp using the simulation toolbox MNPBEM for plasmonic nanoparticles.[20,21] We first consider a silver nanorod with dimensions 200 × 65 × 30 nm3 and compute the loss spectra for the three selected impact parameters indicated in Figure a. The two prominent loss peaks at low energies can be attributed to the dipole and quadrupole plasmon modes. Corresponding EELS maps at the resonance frequencies are shown for a few selected electron propagation directions (rotation angles) in Figure c. The mode profiles are reminiscent of the dipole and quadrupole surface charge distributions.[8] For the decomposition of eq into modes E(r,ω), we use the information about the nanoparticle shape, which in experiment can be obtained from additional high-angle annular dark-field (HAADF) data[22,23] and compute the 50 natural oscillation modes of lowest energy (see Methods). Figure b shows the modulus of coefficients C obtained from either a compressed sensing or least-squares optimization. Although the two approaches give quite different C distributions, the back-projected EELS maps, obtained by assembling the dyadic Green tensor using eq and computing Γ̃surf from eq , both are in almost perfect agreement with the original Γexp maps.

Figure 1

EELS spectra and maps for a silver nanorod. (a) EELS spectra recorded at the positions indicated in the inset. The peaks at approximately 1.5 and 2.7 eV are attributed to the dipole and quadrupole plasmon mode. (b) Mode decomposition of the dipole and quadrupole mode from the collection of rotated EELS maps, using either the least-squares minimization of eq or the compressed sensing optimization of eq . For each mode, the coefficients C are normalized to unity. (c) Selected EELS maps for dipole (upper part) and quadrupole (lower part) mode and for different electron propagation directions (rotation angles), as computed with the MNPBEM toolbox.[20,21] (d) Back projected EELS maps for the C distribution obtained from the compressed sensing optimization, using eq for the Green function decomposition and eq for the calculation of the loss probabilities. (e) Same as panel (d) but for C distribution obtained from the least-squares optimization.

Having obtained the C values from the optimizations of eqs and 9, we can use eq to approximately reconstruct the dyadic Green tensor, which allows us to compute any electrodynamic response function for the plasmonic nanorod. In the following we consider the projected photonic LDOS of eq . Figure shows the true and reconstructed LDOS maps and compares the quality of compressed sensing and least-squares optimizations. In particular, the inspection of panels (b) and (c), which report the LDOS in a plane 20 nm above the nanorod, reveals that the compressed sensing results are in very good agreement with the true LDOS values, whereas the least-squares optimization completely fails to provide even qualitative agreement. This finding seems at first sight surprising since both optimization approaches were previously capable of reconstructing the experimental EELS data almost perfectly, as shown Figure c–e. We attribute the least-squares shortcoming to the fact that the EELS loss of eq is governed by the long-range tails of the particle plasmon field distributions, with which the passing electron predominantly interacts, whereas the LDOS of eq is governed by the short-range evanescent field components. Thus, when the optimization has no strong bias on the C determination, it comes up with the proper long-range components, resulting in high-quality EELS maps shown in Figure e, but fails for the short-range components, which contribute little to the minimization function of eq . In contrast, the compressed sensing optimization of eq seeks for a C distribution with as few nonzero components as possible. For suitable basis functions E, this bias helps to properly select those modes that contribute little but still noticeably to the loss probability of eq . We emphasize that such a bias for selecting a sparse expansion distribution is by no means unique to the problem of our present concern, but has been previously highlighted in various studies, for example, in the context of plasmon tomography[10] or single-pixel cameras,[24] and lies at the heart of the compressed sensing algorithm.

Figure 2

Photonic LDOS of eq and reconstructed LDOS. (a) Three-dimensional LDOS distribution, as computed with the MNPBEM toolbox (LDOS),[20] and the distributions reconstructed from the compressed sensing (CS) and least-squares (LSQ) optimizations. The projected LDOS ρ(r,ω) is shown for different projection directions n̂ = x̂,ŷ,ẑ. (b) LDOS density map in a plane 20 nm above the nanoparticle, as reconstructed from the compressed sensing optimization. The lower (upper) part of each panel shows the dipole (quadrupole) mode, the left (right) part shows the true (reconstructed) LDOS. (c) Same as panel (b) but for least-squares optimization. The reconstructed least-squares LDOS has also negative contributions, which are set to zero for clarity.

An advantage of compressed sensing is that the reconstruction can, in general, be performed, even with a very limited amount of measurement data, and the quality of the reconstructed data is usually not strongly affected by noise.[18] In Figure we show reconstructed EELS and LDOS maps for the small number of impact parameters and rotation angles shown in the first row of measurement data. As can be seen, the quality of the reconstructed data is extremely good, despite the limited amount of measurement data. This might be beneficial for EELS experiments that typically suffer from a limited amount of rotation angles (missing wedge problem) and where the number of measurement points is often kept low to avoid sample contamination.

Figure 3

Compressed sensing reconstruction for a strongly reduced number of measurement points. The first row shows the measurement data for a few rotation angles. In the second row we compare the EELS data for a finer sampling mesh (upper part of panel) with the reconstructed signal (lower part), finding almost perfect agreement. The last row reports the true (upper part of panel) and reconstructed (lower part) LDOS maps in a plane 20 nm above the nanorod.

Finally, in Figure we compare LDOS maps with reconstructed maps for (a,b) a bowtie nanoparticle and (c,d) a cube. For the bowtie geometry, we show the LDOS for the two plasmon modes of lowest energy, which can be labeled as bonding and antibonding according to the parallel and antiparallel orientation of the dipole moments of the individual nanotriangles.[25] The agreement between the true and reconstructed LDOS maps is very good; in particular, one can clearly observe the strongly increased LDOS enhancement in the gap region. For the cube, we show the dipole and corner modes of lowest energy,[10] finding fair agreement between the true and the reconstructed LDOS maps. We attribute the small differences to problems of our algorithm when dealing with degenerate modes of symmetric particles, which might be improved by explicitly accounting for mode symmetries.[26]

Figure 4

(a) True (upper row) and reconstructed (lower row) LDOS for a bowtie geometry (total size 215 × 85 × 30 nm3 and 10 nm gap) and for the bonding and antibonding modes of lowest energy. Color code is identical to Figure . (b) Density map of LDOS in a plane 20 nm above the bowtie structure. (c) True (left) and reconstructed (right) LDOS for a cube with 150 nm side length and for the dipole and corner modes of lowest energy.[10] (d) Density map of LDOS in a plane 30 nm above the cube.

Summary and Discussion

To summarize, we have shown how to extract the dyadic Green tensor of Maxwell’s theory from a collection of EELS maps recorded for different electron propagation directions (rotation angles). Our reconstruction scheme is based on a singular-value decomposition of the Green tensor and a compressed-sensing optimization for the expansion coefficients. We have demonstrated the applicability of our approach for various elementary nanoparticle shapes. We foresee several improvements for plasmon tomography based on EELS. On the experimental side, electron holography[22] can provide additional information and could allow to disentangle the excitation and measurement channels of plasmonic EELS. On the theoretical side, the presented reconstruction scheme works surprisingly well for most nanoparticle geometries, but further work is needed to clarify the role of various ingredients.

First, there are several possibilities for choosing the basis functions for the decomposition of the dyadic Green tensor, eq . In this work we have chosen biorthogonal “constant flux states”[27] that are the eigenstates of the Green function evaluated for real frequencies (see Methods). They have the advantage that they can be computed rather straightforwardly, even in the case of degenerate or near-degenerate modes; on the other hand, they have to be computed for each loss energy separately, and several of these modes can govern the plasmonic response. Another possibility for a basis are the quasi normal modes evaluated at the poles of the Green function in complex frequency space.[14−17] The computation of these modes requires an iterative solution scheme,[17] however, once they are computed, they can be used for a large frequency range, and in general, the plasmonic response is only governed by very few of these modes.

In this work we have considered the situation where the basis is already computed for the true nanoparticle shape and have shown that even in this case the EELS tomography scheme can be quite tricky. However, our approach is less restrictive than it may appear: in principle, for electron beams not penetrating the nanoparticle, any basis with modes being solutions of the free-space Maxwell’s equations can be employed. Thus, even if a slightly different particle shape or dielectric material is considered in the computation of the basis, this will not necessarily degrade the quality of the reconstruction. In this case, it might be beneficial to adapt our approach such that (i) the modes for the Green function decomposition are expanded in a given nonideal basis and (ii) the compressed sensing algorithm seeks for a minimum number of decomposition modes. Here it might be advantageous to use quasi normal modes, because the same few modes could be optimized for a whole range of loss energies, thus, imposing stronger restrictions in comparison to an independent optimization at individual loss energies.

Although further work is needed to establish EELS tomography of plasmonic nanoparticles as a robust and out-of-the-box scheme, we believe that our present work provides an important step forward for reconstructing electrodynamic quantities from EELS measurements and makes significant progress with respect to the recently developed tomography schemes that were bound to quasistatic approximation and other restrictive assumptions.

Methods

Simulations

In our simulation approach, we compute the LDOS and EELS spectra using the MNPBEM toolbox[20,21] and a silver dielectric function extracted from optical experiments.[28]

Mode Decomposition

For the mode decomposition of eq , we follow the prescription of García de Abajo et al.[29] and compute the natural oscillation modes through diagonalization of the Σ matrix, see eq 21 of ref (29) for details, keeping for the solution of the inverse problem the 50 modes of lowest energy. A higher number of modes did not show a significant improvement in the reconstruction results. For our mode decomposition it turns out to be convenient to use a biorthogonal basis, similarly to the quasistatic case.[30] Our approach closely follows recent related work,[17] and we introduce the right and left eigenmodes E(r,ω) and Ẽ(r′,ω) associated with the Σ matrix, respectively. Instead of the decomposition of eq , we then useand, accordingly, also modify eq . The biorthogonal expansion turns out to be advantageous in particular for nanoparticles with degenerate modes, as it automatically guarantees proper mode orthonormalization.

Compressed Sensing

The least-squares optimization is performed with the built-in Matlab functions, for the compressed sensing optimization we use the YALL1 software freely available at http://yall1.blogs.rice.edu/. We set the mixing parameter μ = 5 × 10–2, and the stopping tolerance has a value of 10–4. We take 12 rotated EEL-maps for each structure with equidistant angles between 0 and 180°, each map consisting of 31 × 51 points. To speed up the optimization process, we take only 2000 random measurement points of the generated maps. Further, only measurement points with distance more than 15 nm away from the particle surface are used for optimization. For the volume visualization of the LDOS, we use the MatVTK software freely available at http://hdl.handle.net/10380/3076.