Tomographic Reconstruction of Quasistatic Surface Polariton Fields

We theoretically investigate the tomographic reconstruction of the three-dimensional photonic environment of nanoparticles. As input for our reconstruction we use electron energy loss spectroscopy (EELS) maps for different rotation angles. We perform the tomographic reconstruction of surface polariton fields for smooth and rough nanorods and compare the reconstructed and simulated photonic local density of states, which are shown to be in very good agreement. Using these results, we critically examine the potential of our tomography scheme and discuss limitations and directions for future developments.


■ INTRODUCTION
Nano optics deals with light confinement at the nanoscale. 1,2 This is achieved by binding light to surface resonances of nanoparticles, such as surface plasmon polaritons for metallic nanoparticles 3 or surface phonon polaritons for dielectric nanoparticles. 4,5 These resonances come along with strongly localized fields and allow squeezing light into extreme subwavelength volumes, which can be exploited for various applications. 6 Because of the diffraction limit of light, the strongly localized fields cannot be directly imaged in optical microscopy. In recent years, electron energy loss spectroscopy (EELS) has become a highly successful technique for imaging electromagnetic fields at the nanoscale and with high energy resolution. 7−10 In EELS, swift electrons pass by or through a nanoparticle and loose with a certain probability energy by exciting surface resonances. By raster-scanning the electron beam over the specimen and measuring the number of electrons that have lost a certain amount of energy, one obtains information about the electromagnetic fields at the nanoscale. 2,11 However, the technique does not provide direct information about the three-dimensional fields but only about the averaged interaction along the entire electron trajectory.
EELS tomography is a variant of electron tomography, 12 where the three-dimensional structure of a specimen is reconstructed from a collection of transmission electron micrographs for various tilt angles. In EELS, the reconstruction is complicated by the fact that the loss does not occur at a specific position of the specimen, but is a highly nonlocal process. 11 EELS tomography of surface plasmons was first suggested independently in refs 13 and 14, where the latter paper demonstrated experimentally the reconstruction of localized surface plasmon modes for a silver nanocube. While these seminal papers employed the quasistatic approximation, 2,11 successive work showed how to extend the scheme to full retardation 15 and demonstrated its applicability for single and coupled silver nanoparticles. 16,17 In a recent paper, 18 we have brought EELS tomography from the optical to the mid-infrared regime and have demonstrated experimentally the reconstruction of localized surface phonon polaritons for a MgO nanocube. Contrary to surface plasmon polaritons, the use of the quasistatic approximation is perfectly justified for surface phonon polaritons sustained by nanoparticles with dimensions of a few hundred nanometers. This considerably simplifies the methodology for the tomographic reconstruction. While going full circle from the quasistatic tomography of surface plasmon polaritons in our initial work 13 to quasistatic tomography of surface phonon polaritons, 18 we have gained quite some understanding of the critical elements in EELS tomography, and our approach has matured considerably. The time is ripe for a critical re-examination and reinterpretation of our tomography scheme.
In this paper we present a theoretical study of EELS tomography for prototypical dielectric nanoparticles. We submit a tilt series of simulated EELS maps to our tomography scheme in order to extract parameters characterizing the nanophotonic environment. For this parametrized photonic environment, we compute the photonic local density of states (LDOS), 1,2,19 which is compared with independent simulation results. From this comparison, we examine the strengths and weaknesses of our tomographic reconstruction scheme.
The photonic LDOS is a concept borrowed from solid state physics and accounts for the number of photonic modes per unit frequency and volume. In free space, the photonic LDOS is 1 (1) where ω is the angular frequency and c the speed of light. The photonic LDOS governs the power dissipated by an oscillating dipole through where p is the oscillator's dipole moment and ε 0 the free-space permittivity. Alternatively, we can relate via P 0 = ℏωγ 0 the power dissipation to the decay rate γ 0 of a quantum emitter. The concept of the photonic LDOS comes to full glory in nanophotonics, where the light−matter interaction becomes dramatically enhanced through surface excitations of nanoparticles, such as surface plasmon or phonon polaritons. The enhancement of the photonic LDOS ρ(ω) can be in the range of hundreds to thousands in comparison to its free-space value ρ 0 (ω). 20 Correspondingly, quantum emitters can transfer energy to the nanophotonic environment more efficiently, and their decay rate or power dissipation P is increased by the LDOS enhancement according to = P P : : Below we will compute the LDOS enhancement ρ:ρ 0 using the photonic environment reconstructed from EELS maps. It is obvious that electrons and oscillating dipoles couple quite differently to the nanophotonic environment. For this reason, the LDOS reconstruction from EELS data is quite delicate and provides a stringent testbed for our tomography approach.
We have organized our paper as follows. In the Theory section we present the theory and methodology of our tomographic reconstruction. We have tried to keep the presentation as compact and brief as possible and refer to the literature for the detailed derivations whenever possible. Some technical issues are transferred to Appendix, Details about Orthogonal Matrix. In the Results section we present the tomography results for smooth and rough nanorods and compare the reconstructed and the simulated photonic LDOS. Finally, in the Summary section we put our tomography into a broader context, examine critically the strengths and weaknesses of our approach, and identify lines for future research.

■ THEORY
For MgO nanoparticles the surface phonon polariton energies hν are of the order of 100 meV, corresponding to a free-space wavelength λ = c/ν ∼ 12 μm. For nanoparticle dimensions of approximately hundred nanometers we can thus safely introduce the quasistatic approximation, 2 where the electric field is expressed in terms of a quasistatic potential V(r) through E(r) = −∇V(r) and we keep the frequency dependence of the permittivity functions ε(ω).
Green's Functions. In the following we consider the problem depicted in Figure 1a, where a charge located at position r′ interacts with a dielectric nanoparticle situated in a background medium with dielectric constant ε 0 . Green's functions provide an elegant and efficient method for solving such problems. We first introduce the Green's function defined through 2,21 which gives the potential at position r for a unit charge located at position r′. In an unbounded medium, the Green's function would be given by the usual expression 5) and the potential associated with a charge distribution ρ(r) can be expressed as In presence of the nanoparticle, this incoming potential will induce a reflected potential associated with the particle response. To account for this, we split the total Green's function into two parts 0 refl (7) where the reflected part is a solution of Laplace's equation which is chosen such that Maxwell's boundary conditions are fulfilled at the nanoparticle boundary. Suppose for a moment that the reflected Green's function is at hand. It can then be shown that in EELS the loss probability is related to the reflected Green's function via 2,11 where R 0 = (x 0 , y 0 ) is the impact parameter of the electron beam propagating along the z direction (aloof geometry), ℏω is the loss energy, and ρ el (r) is the charge distribution of the swift electron. The term in brackets of eq 8 accounts for a selfinteraction process where the swift electron polarizes the nanoparticle and the polarization acts back on the electron. This nonlocal response is mediated by the reflected Green's function. Similarly, the power dissipated by a dipole oscillating with frequency ω becomes 2 In free space the Green's function G 0 (r, r′) gives the potential at position r for a unit charge located at position r′. In presence of a nanoparticle one must additionally add a reflected Green's function that accounts for the nanoparticle response. (b) The reflected Green's function can be expanded using a complete set of eigenpotentials V k (r). In our tomography scheme we can also start from the modes associated with a simpler reference boundary ∂Ω 0 rather than the actual nanoparticle boundary ∂Ω, and expand the reflected Green's function using the reference modes. For details, see text.
where P 0 is the free-space dissipation, p is the dipole moment, and r 0 is the position of the dipole. The ratio P:P 0 gives the enhancement of the photonic LDOS, see also eq 3. The expressions given in eqs 8 and 9 are two examples for the enhancement of light−matter interactions in the presence of nanoparticles and show that the nanophotonic environment is fully characterized upon knowledge of the reflected Green's function. Eigenmode Decomposition. A powerful and convenient representation of the reflected Green's function is in terms of geometric eigenmodes u k (s) and eigenvalues λ k , where s is a position located on the boundary of the nanoparticle. 2,22,23 These eigenmodes form a complete set of basis functions. To each eigenmode we can associate an eigenpotential which is a solution of Laplace's equation that fulfills Maxwell's boundary conditions at the nanoparticle boundary. We can then decompose the reflected Green's function outside the nanoparticle in terms of these eigenpotentials via 2,23 where Λ(ω) is an expression that solely depends on the permittivities of the nanoparticle and the embedding medium. Inserting eq 11 into the EELS loss probability of eq 8 leads us to with the line shape function Equation 12 is a particularly useful decomposition of the loss probability in terms of surface phonon polariton eigenmodes. Each eigenmode contributes with the line shape function L k (ω) and the oscillator strength given by the square modulus term, which is governed by the interaction energy between the charge distribution of the swift electron and the eigenpotential V k (r). Similarly, the power dissipated by an oscillating dipole of eq 9 can be decomposed into eigenmodes via with a corresponding interpretation in terms of line shape functions and oscillator strengths. From the dissipated power one can obtain the photonic LDOS using eqs 1 and 3, where one often additionally averages over all dipole orientations to account for the random orientation of quantum emitters in typical experiments. 1 Tomographic Reconstruction of Eigenmodes. It is apparent from eqs 12 and 14 that we can compute the EELS loss probability Γ(R 0 , ω) and the LDOS enhancement P:P 0 , or any other related response function, once the geometric eigenmodes u k (s) and the line shape function L k (ω) are at hand. Expressed differently, the nanophotonic environment is fully characterized upon knowledge of u k (s) and L k (ω). We can now formulate the goal of our tomography approach. Suppose that we are in possession of the EELS loss probabilities Γ(R 0 , ω), ideally for various impact parameters and electron propagation directions, but do not know the eigenmodes u k (s) and line shape functions L k (ω): can we obtain through solution of an inverse problem a viable approximation for u k (s) and L k (ω)? And if yes, how?
Optimization for Modes on the Nanoparticle Boundary. Consider first the situation that the nanoparticle boundary is known and that we are seeking for the linshape functions and eigenmodes L k , u k (s). This corresponds to the situation previously investigated in ref 18. Let s u ( ) 0 be a complete set of basis functions on the boundary. We shall refer to these modes as reference modes. As shown in the Appendix, Details about Orthogonal Matrix, the eigenpotentials of eq 10 can be expanded in terms of these modes via with being an orthogonal matrix. We can now formulate the tomographic reconstruction scheme for a given set of experimental EELS maps.
1. Find some reference modes s u ( ) 0 whose gross features are expected to be similar to those of the true eigenmodes u k (s). This point is irrelevant for a complete basis, but becomes crucial for actual reconstructions where the basis has to be truncated. 2. Start with some initial guess for the line shape function L k and orthogonal matrix , and compute the reprojected maps via eq 12. Use an optimization routine for L k , to obtain the best possible agreement between experiment and reprojection. Note that in principle L k (ω) depends on frequency, but for a fixed loss energy the line shape functions can be treated as mere numbers. 3. Use the optimized parameters to compute other quantities, such as the photonic LDOS. Optimization for Modes on Reference Boundary. The above scheme can be also generalized to cases where the true nanoparticle boundary ∂Ω is not known or is too complicated to be used in actual reconstructions. We start by introducing a reference boundary ∂Ω 0 that fully encapsulates the nanoparticle, see also Figure 1b. In our modified approach we are not aiming for a reconstruction of the eigenmodes u k (s) themselves, but of the eigenpotentials of eq 10 outside of the reference boundary. There they can be expressed as generic solutions of Laplace's equation 21 where σ k (s) specifies the normal derivative of the potential on ∂Ω 0 (von Neumann boundary condition). We can now use a complete set of basis functions s u ( ) 0 on ∂Ω 0 for the expansion of σ k (s) to arrive at where is a nonorthogonal matrix formed by the expansion coefficients. The tomographic reconstruction can now be performed in complete analogy to the scheme presented above, with the only exception that has to be replaced by a nonorthogonal matrix.

ACS Photonics
pubs.acs.org/journal/apchd5 Article Optimization Loop. In the following we discuss the optimization procedure in slightly more detail, see also Figures 2 and 3. We provide a unified description for the optimizations using modes defined on either the nanoparticle or the reference boundary. In the first case, is an orthogonal matrix. In our computational approach we have to truncate the basis and keep only n representative modes, where n is of the order of several tens to hundreds. Correspondingly, × n n is a matrix of size n × n, see also Appendix, Details about Orthogonal Matrix, for the parametrization of this matrix. In the case of a reference boundary, is a full matrix. In principle we can now use different truncation numbers m and n for the reconstructed eigenpotentials and basis functions, respectively, and × m n becomes a matrix of size m × n. In most cases it is sufficient to consider around 10 eigenpotentials, whereas the truncation number for the basis should be chosen considerably larger. Let be a set of impact parameters and tilt angles for a fixed loss energy, and Γ exp (x i ) is the corresponding experimental EELS maps. We only consider aloof electron trajectories that do not penetrate the nanoparticle. The interaction energy between the swift electron and a reference mode The experimental EELS maps Γ exp are obtained for a specific loss energy and for various rotation angles, we only keep aloof electron trajectories that do not penetrate the smooth rod. We start with some initial guess for the optimization parameters L k , and compute the reprojected EELS maps Γ using eq 20. These parameters are optimized until a local minimum is reached by the optimization algorithm. In the lowest row we show the relative error |Γ exp − Γ|:Γ exp between the experimental and optimized maps. The solid lines indicate the contours for an error of 0.1%. Once the parameters L k , are at hand, we can compute other quantities such as the photonics LDOS. the eigenmodes of the Laplace− Beltrami operator. Using the optimized parameters L k , we mix the modes to obtain the reconstructed modes shown on the right-hand side for the dipole and quadrupole resonances. For the smooth rod, is an orthogonal matrix of size n × n, where n is the truncation number of the basis. For the rough rod, is a full matrix of size m × n, where m is the number of eigenpotentials to be reconstructed. From the knowledge of we can compute the geometric eigenpotentials V k (r) outside the reference boundary. The bar plot on the right-hand side reports the reconstructed line shape parameters for the dipole (blue) and quadrupole (red) resonances, the modes are sorted in decreasing order of L k and the largest contributions are due to the modes shown in the insets.
where V el (r) is the potential associated with the charge distribution ρ el (r). When the nanoparticle boundary is known, the reference boundary in the above boundary integral is identical to ∂Ω. The loss probability of eq 12 can then be written in the compact form We can now define a cost function (21) that gives the "distance" between the experimental and reprojected EELS maps. This cost function is submitted to an optimization routine, such as a conjugate-gradient or quasi-Newton one, 24 which provides us with the optimized expressions for L k , . Some details about the parametrization of the orthogonal matrix, as well as the computation of the derivative of the cost function with respect to the optimization parameters are given in Appendix, Details about Orthogonal Matrix.

■ RESULTS
In ref 18 we have applied our tomography scheme to experimental EELS maps for a MgO nanocube. In this work we proceed differently and investigate the working principle of our tomography scheme using simulated data only.
1. We first compute for each loss energy EELS maps for a series of rotation angles, see also Figure 2. To be consistent with our previous notation, we denote these simulated EELS maps as Γ exp and will refer to them as experimental EELS maps. 2. These maps are submitted to our tomography scheme based on eq 21 in order to obtain the optimized parameters L k and that specify the nanophotonic environment. 3. Using eq 14 together with the optimized parameters, we compute the photonic LDOS and will refer to it as the reconstructed photonic LDOS. 4. Using eq 9, we compute the photonic LDOS directly, with a simulation approach to be discussed below, and will refer to it as the simulated photonic LDOS.
For ideal reconstruction, the simulated and reconstructed LDOS maps should be identical. Any deviation between the two maps can thus be attributed to deficiencies of our approach, caused for instance by the truncation of the reference basis s u ( ) 0 or a trapping of the optimization algorithm in a local minimum. We apply our tomography scheme to prototypical systems of a smooth and rough nanorod with a diameter to length ratio of approximately 1:2.5, see also Figure 4 and ref 25 for a detailed discussion of the rod modes. The rough rod has been generated by adding stochastic height variations to the smooth surface of an ideal nanoparticle following the prescription given in ref 26. We shall not be concerned whether such nanoparticles can indeed be fabricated with the material system under investigation. As we are working within the quasistatic regime, the actual size of the nanorods is irrelevant, and the results can be easily scaled to any size.
Computational Details. All our simulations are performed with the quasistatic classes of the NANOBEM toolbox, 27 which is based on a Galerkin scheme with linear shape elements. See, for example, ref 2 for a detailed discussion. The parametrization of the MgO dielectric function is the same as in refs 18 and 28. The nanorod boundaries are discretized using more than 3000 boundary elements of triangular shape. We checked that for such fine discretizations we obtained converged results. As for the EELS simulations, we consider the limit of large electron velocities v, where the potential for a swift electron with impact parameter R 0 takes the form with e being the elementary charge and R = (x, y). We have previously shown 29 that this simplified expression gives almost the same results as simulations based on the full Maxwell's equations.
As for the reference modes s u ( ) 0 , we did not choose the usual geometric eigenmodes 2,23 for two reasons. First, in order to demonstrate that our approach indeed works for any meaningful set of basis functions. Second, we observed that the geometric eigenmodes computed with the NANOBEM toolbox are often strongly localized around sharp corners or edges, such that a large number of such modes would be needed for a useful expansion. In this work, we choose for s u ( ) 0 the eigenmodes of the Laplace−Beltrami operator, which is a generalization of the Laplace operator for curved boundaries and is known to provide extremely smooth basis functions. 30 The modes were additionally orthogonalized using eq 25.
In our optimization approach we truncate the Laplace-Beltrami basis using the n modes of highest eigenvalue, where a value of n ≈ 100 turned out to be a good compromise between reasonably fast optimizations and sufficiently accurate results. The optimization was performed with the built in MATLAB function fminunc using a quasi-Newton algorithm together with a relatively small function and optimality tolerance of 10 −8 . In all our simulations, we typically needed about 2000 iterations to reach a local minimum.
Smooth Rod. We start by discussing the smooth rod shown in Figure 4. The loss spectra exhibit three peaks, which can be attributed to a dipolar mode (70 meV), a quadrupolar mode (80 . Loss spectra for smooth and rough nanorod, and for impact parameters located on the long (blue) and short (red) rod axis, see inset. We consider aloof electron trajectories with a propagation direction out of the image plane. One observes a dipole resonance around 70 meV, a quadrupole resonance around 80 meV, and a peak attributed to a multitude of modes around 90 meV. meV), and a peak that is composed of a multitude of modes (88 meV). For the smooth rod, the reference boundary ∂Ω 0 is identical to the true nanoparticle boundary ∂Ω. Note that the Laplace−Beltrami eigenmodes provide a (truncated) basis that does not coincide with the true geometric eigenmodes. Simulated and reprojected maps originating from our optimization algorithm are typically extremely similar, see lowest row in Figure 2 for the more difficult case of the rough nanorod. Figure 5 shows the simulated and reconstructed LDOS maps in the symmetry plane (left column), in planes away from the rod (other columns), and for the loss energies reported in the figure. We first consider the dipole mode shown in panel (a). The LDOS can be interpreted for an oscillating dipole as the enhancement of the dissipated power, see eq 9, throughout we average over all possible dipole orientations. Close to the rod, an oscillating dipole couples with comparable strength to all surface phonon polariton modes. This can be seen both in the symmetry plane of the rod (first column, logarithmic color scale), as well as in the plane closest to the rod (second column, linear color scale), where the photonic LDOS is large and unstructured close to the rod boundary. When moving away from the rod (other columns from left to right), the coupling strength between the oscillating dipole and the rod resonance modes have different distance dependencies, which are governed by the oscillator strengths given in eq 14. For the chosen loss energy the dipolar rod mode becomes strongest at larger distances, as can be inferred from the two lobes in the LDOS maps located at the rod caps. Figure 5b,c shows results for the reconstructed LDOS using (b) n = 100 and (c) 20 Laplace−Beltrami reference modes. Further away from the rod, the simulated and reconstructed results agree well for both truncation numbers n. For distances closer to the rod, the larger number of eigenmodes provides better agreement. This is in accordance to our previous reasoning that oscillating dipoles close to the rod couple to a larger number of eigenmodes, and thus a larger number of modes is needed for the reconstruction.
In Figure 6 we give a quantitative comparison between the simulated (full lines) and reconstructed (dashed lines) LDOS values for cuts along the long rod axis and for dipole positions outside the nanoparticle. The true LDOS enhancement would depend on the actual size of the nanorod; for simplicity, we give the results in arbitrary units. Also, the reconstruced LDOS cuts are scaled by a constant factor, where it is not obvious how this factor could be obtained in the absence of EELS loss probabilities given in absolute numbers. We here do not enter into the question of how to extract the absolute numbers of the reconstructed LDOS. Besides this unknown prefactor, the simulated and reconstructed LDOS values agree extremely well, with the possible exception of the smallest distances where a larger number of eigenmodes might be needed.
Finally, in the remaining panels of Figure 5, we compare the simulated and reconstructed LDOS for the (d, e) quadrupolar rod mode and the (f, g) multitude of modes. It can be seen that the reconstruction works well for the quadrupolar mode. Comparison with results for n = 20 (not shown) reveal that in this case a larger number of eigenmodes is strictly needed to obtain good agreement. For the multitude of modes shown in panels (f) and (g), the agreement between simulation and reconstruction is reasonable, but not overly good. In particular, for the smallest distances, the reconstructed maps show sharp or asymmetric features, which are absent in the simulated maps. From these results we conclude that the LDOS reconstruction works best for loss peaks that are governed by a few modes only.
In Figures S1 and S2 we also show results for a nanorod with reduced symmetry, which is obtained by squeezing the rod in all axes directions. Again, the simulated and reconstructed LDOS maps are in very good agreement. We also investigate in Figures  S3 and S4 the influence of the number of optimization iterations on the reconstructed LDOS maps. While the gross features of the maps are already reconstructed after a few tens to hundred iterations, it takes a few thousand iterations until reaching convergence.
Rough Rod. The case of the rough rod shown in Figure 4 is considerably more difficult. We keep considering the same reference modes as for the smooth rod, and select the reference boundary ∂Ω 0 such that it fully encapsulated the boundary ∂Ω of the rough rod. Note that this reference boundary is identical to the one of the smooth rod. Figure 2 shows for the dipolar mode the simulated ("experimental") EELS maps and the reprojected ones. The relative error between these maps is small throughout.
In Figure 7 we show the simulated and reconstructed LDOS maps for the rough nanorod. We compare different planes that are (a−f) parallel and (a*−f*) perpendicular to the electron beam direction. The main difference between these two configurations is that in the parallel case we reconstruct the LDOS throughout in regions through which swift electrons have traveled. In contrast, for the perpendicular case we reconstruct the LDOS also in planes above the nanoparticle through which no electron has traveled because of our restriction to aloof trajectories.
Let us consider the parallel case first. With the possible exception of the smallest distance, the agreement between simulated and reconstructed LDOS maps is extremely good, both for the dipolar and quadrupolar modes. Both asymmetries as well as hot spots, caused by localized fields in the vicinity of protrusions of the rough rod, are well reproduced by our tomography scheme. Things somewhat change for the perpendicular geometry shown in panels (b*) and (e*), where the comparison is reasonable but not overly good. We performed additional simulations where the tilt series for Γ exp is complemented by EELS maps where the nanorod is first rotated around the x-axis by 90°before being submitted to the same tilt series around y. As can be seen in panels (c*) and (f*), with this procedure we again obtain extremely good agreement between simulated and reconstructed LDOS maps. This shows that our tomography scheme works best for regions through which electrons have traveled.
We finally investigate in Figure 8 the impact of the cutoff parameters (m, n) on the reconstructed LDOS for the dipole mode. Recall that m is the number of eigenpotentials to be reconstructed and n is the cutoff parameter for the basis functions. Close to the particle (second column) a larger number n of basis states leads to a better agreement with the simulated LDOS, shown in the first row. When moving away from the nanoparticle, the agreement between simulated and reconstructed LDOS maps is very good for all chosen simulation parameters. This demonstrates that our reconstruction scheme is robust and that the optimization results do not depend decisively on the input parameters.

■ DISCUSSION
In the previous sections we have presented the methodology of our tomographic reconstruction scheme and have investigated the approach for prototypical nanophotonic structures. In this section we start by discussing our scheme within a broader context, and then address limitations, dos and don'ts, as well as extensions of our tomographic reconstructions.
Working Principle. The basic working principle of our tomographic reconstruction is shown in Figure 9 and consists of the triad formed by experiment, resonance modes, and reference modes. In short, the resonance modes are needed to formulate the theory, and the reference modes to provide a parametrization of the nanophotonic environment and to perform the actual reconstruction. The experimental data are the primary resource for the reconstruction. For this reason, the quality of the experimental data directly influences the quality of the tomographic reconstruction. Some further considerations about experiments will be given below.
Resonance Modes. The nanophotonic environment outside the nanoparticle is fully characterized in terms of the reflected Green's function of eq 11, which we repeat here in compact form (23) M k and L k are the real and imaginary parts of the term given in brackets of eq 11. The eigenpotentials V k (r) provide the preferred physical basis, only with this basis the reflected Green's function can be written in the diagonal form of eq 23. A similar decomposition of the Green's function can be also obtained in the retarded case when using quasinormal modes, 31−33 as will be discussed below. For this reason, from here on we use the more general expression of resonance modes rather than geometric eigenmodes, for which we have developed our theory so far.
With these modes, both the EELS loss probability of eq 12 as well as the power dissipation of an oscillating dipole, eq 14, can be written as the sum over individual loss channels. With any other basis one would obtain some kind of mixing between different modes. This particular form has the additional advantage that the line shape function L k is always positive, at least for lossy materials, which can be used in our optimization procedure as a constraint, see Appendix, Details about Orthogonal Matrix. Note that our tomography scheme only allows for the reconstruction of L k , which accounts for the loss properties of the nanophotonic environment, but not for the propagation properties described by M k . As EELS and LDOS account for energy losses of electrons and oscillating dipoles, respectively, this is not a problem here. However, additional experimental input or a reconstruction for various loss energies together with a Kramer−Kronig analysis would be needed for a reconstruction of M k .
To summarize this part, resonance modes are needed to formulate the abstract theory, without making contact to the actual shape or composition of the nanoparticle. Without resonance modes it would be unclear which properties of the nanophotonic environment govern EELS and LDOS, and which properties can be reconstructed using an inverse scheme. However, at no point of our approach we require explicit knowledge of the actual form of the resonance modes or line shape functions.
Reference Modes. The reference modes are the device needed for the actual reconstruction. They allow for a suitable parametrization of the nanophotonic environment, where the viable parameters can be extracted from the optimization loop using the experimental and reprojected EELS maps. In principle, Figure 7. Same as Figure 5, but for a rough nanorod and for the (a−c) dipolar and (d−f) quadrupolar rod resonances. In the reconstruction we consider n = 200 reference modes for a smooth nanrod (black contour shown on top) and m = 20 modes to be reconstructed. We compare LDOS values in planes (a−f) parallel and (a*−f*) perpendicular to the electron propagation direction. The lines and dots in the rod center indicate the tilt axis. In panels (b) and (e), we consider a tilt series where the nanoparticle is rotated around the y-axis only, whereas in panels (c) and (f), we additionally consider a rotation around the x-axis by 90°, followed by the same tilt series around y.

ACS Photonics pubs.acs.org/journal/apchd5
Article for a complete basis the choice of the reference modes is irrelevant. However, in all practical cases one has to truncate the basis, which should thus be based on an educated guess and should include the gross features of the expected resonance modes from the outset. Boundary Element Method Approach. Although all our calculations presented here and elsewhere 13,15−18 have been performed using a boundary element method (BEM) approach, it does not play an exceptional role in our tomography scheme. The reference modes are fixed by specifying their values on a properly chosen reference boundary, see Figure 9. Away from the boundary the modes propagate according to Laplace's equation, see eq 10, or the source-free Maxwell's equations in the retarded case. This propagation is reminiscent of Huygens' principle for the wavefront propagation in free space and can be well described within BEM, but otherwise, our tomography makes no particular use of it.
Frequently Asked Questions. All our cards are on the table now. Up to here, we have presented and examined our tomography approach in some detail, and have put it into a broader context. However, a number of open or not fully clear issues remains. In the following we address these issues in the form of frequently asked questions. As will become apparent, only some of these questions can be answered definitely while others remain open. In this sense, the following discussion is meant to summarize our present understanding of the field, to make aware where things can go wrong, and to identify directions for future research.
How much preknowledge is needed? Any tomography or inverse scheme requires some sort of preknowledge, less preknowledge usually makes an approach more general and powerful. In our case, we assume that the nanophotonic enviromnent can be expressed in terms of resonance modes, and that the potentials away from the boundary propagate as solutions of Laplace's equation.
How many reference modes (and which) are required? For any practical reconstruction one has to truncate the basis. The proper choice and truncation of the reference basis thus enters as an additional preknowledge. For spectrally isolated resonances often a few tens of modes suffice, while in other cases up to hundred modes might be needed. We are not aware of any general approach for determining the correct number of reference modes, so we advice potential users to vary the number and to obtain the best cutoff parameter on a case-to-case basis.
Are more reference modes always better? More modes slow down the optimization and require more iterations until reaching convergence. With the simulated EELS data used in

ACS Photonics
pubs.acs.org/journal/apchd5 Article this paper the quality of the reconstruction did not depend decisively on the truncation number.
What is the typical computational cost? Depending on the truncation number, typical optimization times range from one to several minutes on a normal computer. The code developments on top of a BEM solver, such as the NANOBEM one, 27 are moderate. In the future we consider publishing our code to make it accessible to interested users.
Are there dif ferences betweenEELSandLDOS? The question is somewhat odd, obviously EELS accounts for the energy loss of swift electrons and LDOS accounts for the enhancement of the decay rate of oscillating dipoles. On the other hand, the interaction potential of the swift electron to the nanoparticle has a log(r) spatial dependence while the dipole has a 1/r 2 dependence. For this reason, EELS maps are governed by the long-range features of the potential and LDOS maps by the short-range features. The prediction of LDOS maps from experimental EELS maps is thus a challenging and difficult task, and the good agreement between simulated and reprojected LDOS maps reported in this work should not be taken as granted.
Does the optimization always succeed? In all reconstructions considered in this work the optimization algorithm ended up in a minimum. However, there is no guarantee that this is a global minimum. Our results never depended decisively on the initial values for L k , , which we initially set equal to one. Note that zeros would be a bad choice because the derivatives of the cost function with respect to the optimization parameters would equate to zero then. We generally recommend using quasi-Newton optimizations rather than conjugate gradient ones, because they typically access larger portions of the parameter space.
How much experimental input is needed? We will not give too much advice on the experiments here, interested readers might consult our previous work 16−18 to see what worked for us. Depending on the electron microscope, contamination might play a role and might limit the amount of experimental data. As has been discussed before, our tomography seems to work best for regions through which swift electrons have traveled. The reconstruction of blind spots is possible, but the results should be handled with care.
What is the impact of noise? The main purpose of this paper has been to demonstrate that the photonic LDOS can be reconstructed in principle from EELS tomography. As a proof of principle, we have used simulated EELS maps without any noise. In Figure S5, we show EELS maps where noise with a Poissonian distribution has been added artificially. When submitting these noisy maps to our optimization procedure, see Figures S6 and S7, we observe that the LDOS can be only reconstructed reliably sufficiently far away from the nanorod. Close to the rod the added noise leads to artificial features that are not present in the simulated maps. However, we also observe that much better agreement could be achieved if the optimization was terminated after say 100 iterations, with increasing iterations the optimization ends up reconstructing the noisy features rather than the real physical ones. We expect that this deficiency could be overcome with more refined optimization procedures, but will not enter into this topic here.
How to address retardation? For surface phonon polaritons the quasistatic approximation works perfectly, but things might be more problematic for the reconstruction of other surface modes, such as surface plasmon polaritons. In the past we have developed a methodology for surface plasmon tomography, including retardation, and have applied the scheme to experimental EELS data. 15−17 There we used a biorthogonal basis, which shares many features with the resonance modes presented here, but provides no justification for a strictly positive line shape function. For this reason, we opted for a compressed sensing optimization that favors expansions with as few modes as possible, where luckily all of them contributed with a positive weight to the loss probability. In light of our present analysis, we suggest a slight modification of our previous approach. First, in the retarded case the preferred basis is given by quasinormal modes, 27,31−33 which have received considerable interest recently. With these modes we can decompose the dyadic Green's tensor for the full Maxwell's equations in a form similar to eq 23 tomographic reconstruction should be possible along the lines sketched in the present work. There remain a number of open issues, such as the proper choice of the reference modes or the consideration of complex mode functions, but we do not foresee any major roadblock. As a side remark, it is no surprise that the decomposition of the reflected Green's function in terms of resonance modes looks similar in the quasistatic and retarded case: such decompositions are in the spirit of generic singular value decompositions, where the special structure is due to the symmetry property of Green's function originating from the reciprocity theorem of optics.
How to address membranes and grids? Membranes or grids are needed in experiment to support the nanoparticle. One might wonder about the consequences of such a support in our tomographic reconstruction. First, modifications of the Figure 9. Working principle of our tomography scheme. The approach consists of the triad of experiment, resonance modes, and references modes. The experimental EELS maps for various tilt angles provide the basic resource for the reconstruction of the photonic environment. The resonance modes are used to formulate the theory underlying the reconstruction, the reference modes provide the parametrization of the photonic environment and are used for the actual reconstruction. The parameters L k , are obtained through an optimization procedure in order to minimize the difference between the measured and reprojected maps. As the potentials outside the nanoparticle are solutions of Laplace's equation, we can employ a boundary element method (BEM) scheme to express the potentials through their values on a boundary. resonance energies or surface charge distributions of the surface phonon polaritons can be already properly accounted for with the present approach, as has been demonstrated in ref 18. However, in principle also the free-space Green's function of eq 5 should be modified to account for the dielectric environment in the presence of a support. This modified Green's function should be used in eq 10 to propagate away the potentials from the boundary. It is the resulting modification of the electronnanoparticle interaction that has to be considered. Whether this modification has noticeable influence on the results has to be seen.

■ SUMMARY
To summarize, EELS tomography has become a successful scheme for reconstructing the three-dimensional photonic environment of nanoparticles with high spatial and energy resolution. In the past several case studies for plasmonic and photonic nanostructures have provided beautiful results, which would have been hard to achieve with other techniques. Yet, we feel that there is still enough room for improvements and further investigations. In this paper we have given an in-depth study of a prototypical nanophotonic system and have demonstrated that tomographic reconstructions work reliably and without major difficulties, at least for systems where the quasistatic approximation can be employed. We hope that this will motivate more research groups to enter the field, to investigate their systems with the tools presented here, and to continue developing EELS tomography with further improvements.

■ DETAILS ABOUT ORTHOGONAL MATRIX
In this appendix we first derive eq 15. We denote the nanoparticle boundary with ∂Ω and the geometric eigenmodes with u k (s). These eigenmodes form a complete set of basis functions and fulfill the orthogonality relation 2, 22 Inserting this expression into eq 24 and using the orthogonality relation of eq 25, we then immediately observe that is an orthogonal matrix. In our computational approach we employ Cayley's parametrization for orthogonal matrices where is a skew-symmetric matrix with x ij = −x ji . Equation 27 has the advantage that one can perform the derivative x / ij analytically. To show this, we start from Insertion into eq 28 gives We next use , which can be easily proven using + = 1 21 (1 ) and that both terms on the right-hand side commute with (1 ) 1 . We then arrive at our final expression In our optimization of the cost function we express the line shape function through L k = s k 2 , which guarantees that L k is always positive. The optimization algorithm can be significantly accelerated by providing in addition to the value of the cost function also the derivatives with respect to the optimization parameters. Using eq 20 together with eq 29, the derivative of the cost function (eq 21) with respect to the parameters s k and x k of the skew-symmetric matrix can be obtained analytically. Things are considerably easier for a nonorthogonal matrix where the derivatives with respect to the matrix elements can be performed straightforwardly.

Notes
The authors declare no competing financial interest.