Disentangling multipole contributions to collective excitations in fullerenes
DDisentangling multipole contributions to collective excitations in fullerenes
M. Sch¨uler, ∗ J. Berakdar, and Y. Pavlyukh
Institut f¨ur Physik, Martin-Luther-Universit¨at Halle-Wittenberg, 06099 Halle, Germany (Dated: November 5, 2018)Angular resolved electron energy-loss spectroscopy (EELS) gives access to the momentum and the energydispersion of electronic excitations and allows to explore the transition from individual to collective excitations.Dimensionality and geometry play thereby a key role. As a prototypical example we analyze theoretically thecase of Buckminster fullerene C using ab initio calculations based on the time-dependent density-functionaltheory. Utilizing the non-negative matrix factorization method, multipole contributions to various collectivemodes are isolated, imaged in real space, and their energy and momentum dependencies are traced. A possibleexperiment is suggested to access the multipolar excitations selectively via EELS with electron vortex (twisted)beams. Furthermore, we construct an accurate analytical model for the response function. Both the model andthe ab initio cross sections are in excellent agreement with recent experimental data. PACS numbers: 79.20.Uv,31.15.A-,36.40.Gk
Plasmonics, a highly active field at the intersection ofnanophotonics, material science and nanophysics [1], has along history dating back to the original work of Gustav Mieon light scattering from spherical colloid particles [2, 3]. Forextended systems the plasmon response occurs at a frequencyset by the carrier density while in a finite system topology andfinite-size quantum e ff ects play a key role. E.g., for a nano-shell [4–6] in addition to the volume mode, two coupled ultra-violet surface plasmons arise having significant contributionsfrom higher multipoles, as demonstrated below. Such excita-tions can be accessed by optical means as well as by electronenergy-loss spectroscopy (EELS) [7, 8]. Particle-hole ( p - h )excitations and collective modes may “live” in overlappingmomentum-energy domains and couple in a size-dependentway that cannot be understood classically [9–11]. Giantplasmon resonances were measured in buckminster fullereneC [12–17] and explained, e.g., by assuming C to havea constant density of electrons confined to a shell with in-ner ( R ) and outer ( R ) radii (the spherical shell model) [18–20]. Refinements in terms of a semi-classical approximation(SCA) incorporate the quantum-mechanical density extendingout of the shell R < r < R (so-called spill-out density [21]).Time-dependent density functional theory (TDDFT) [21–23]was also employed in a number of calculations [24–26], how-ever, most of them use the jellium model, i.e., the ionic struc-ture is smeared out to a uniform positive background.We present here, to our knowledge, the first atomistic full-fledge TDDFT calculations for EELS from C at finite mo-mentum transfer. We demonstrate the necessity of the full ab-initio approach by unraveling the nature of the variouscontributing plasmonic modes and their multipolar character.This is achieved by analyzing and categorizing the ab ini-tio results by means of the non-negative matrix factorization method [27]. The results are in line with recent experimen-tal findings [28]. The analysis also allows for constructing anaccurate analytical model response function.In first Born approximation for the triply-di ff erential crosssection (TDCS) for detecting an electron with momentum p f , ∗ [email protected] i.e., measuring its solid scattering angle d Ω and energy (cid:15) p f isd σ d ω d Ω = γ q p f p i S ( q , ω ) . (1)Here, p i is the incidence momentum corresponding to an en-ergy (cid:15) p i , γ is the Lorenz factor, q = p f − p i is the momentumtransfer, and ω = (cid:15) p f − (cid:15) p i (atomic units are used through-out). S ( q , ω ) is the dynamical structure factor akin solelyto the target [29]. The fluctuation-dissipation [29] theoremlinks S ( q , ω ) with the non-local, retarded density-density lin-ear response function χ R ( r , r (cid:48) ; t − t (cid:48) ) [29–31] via S ( q , ω ) = − (1 /π )Im[ χ R ( q , − q ; ω )]. On the other hand, χ R ( r , r (cid:48) ; t − t (cid:48) ) de-scribes the change in the system density δ n ( r , t ) upon a small ℓ = ℓ = ℓ = ℓ = ℓ =
84 848484 84 | Im[ δ n ] | (a.u.) | Im[ δ n ] | (a.u.) | Im[ δ n ] | (a.u.) | Im[ δ n ] | (a.u.) | Im[ δ n ] | (a.u.) + + + + + (a) (b)(c) (d)(e) FIG. 1. (Color online) The (cid:96) -resolved constituents of the dynamicalstructure factor of C , | Im[ δ n (cid:96) ( q , ω )] | for (a) (cid:96) =
0, (b) (cid:96) =
1, (c) (cid:96) =
2, (d) (cid:96) =
3, and (e) (cid:96) =
4. The C molecule was treated in standardtruncated icosahedric geometry with bond lengths r C − C = .
445 Åand r C = C = .
390 Å. a r X i v : . [ phy s i c s . a t m - c l u s ] J u l perturbing potential δϕ ( r , t ), i.e. δ n ( r , t ) = − (cid:90) ∞−∞ d t (cid:48) (cid:90) d r (cid:48) χ R ( r , r (cid:48) ; t − t (cid:48) ) δϕ ( r (cid:48) , t (cid:48) ) . (2)The response function is determined by evaluating the den-sity variation with tunable perturbations, as accomplishedvia TDDFT which delivers δ n ( r , t ) upon solving the time-dependent Kohn-Sham (KS) equations [32].Along this line, we utilized the O ctopus package [33, 34],and propagated the KS equations. Kohn-Sham states are rep-resented on a uniform real space grid [35] (0.2 Å grid spac-ing) confined to a sphere with 10 Å radius. For the groundstate we checked the performance of di ff erent typical func-tionals and found that the local-density approximation (LDA)improved by self-interaction correction (SIC) yields fairlygood results. The HOMO (-9.2 eV) is located slightly toolow with respect to the experimental value (-7.6 eV) [15].The band width (which is typically underestimated in DFT)within the LDA + SIC scheme is the largest for the testedfunctionals [36]. LDA-type Troullier-Martins pseudopoten-tials are used to incorporate the influence of the two coreelectrons per C atom, such that only the 240 valence elec-trons accounted for. Gaussian smearing has been employedto deal with the degeneracy of the HOMO. In gas-phase themolecules are randomly oriented. Hence, we have to evalu-ate the spherically averaged structure factor S ( q , ω ). Techni-cally, this can be accomplished by choosing the perturbation δϕ ( r , t ) = I δ ( t ) j (cid:96) ( qr ) Y ∗ (cid:96) m ( Ω r ) [37] where j (cid:96) is the sphericalBessel function and Y (cid:96) m is the spherical harmonic. The pertur-bation strength lies with I = .
01 a.u. well within the regimeof linear response. The perturbed states are then propagatedby the AETRS propagator [38] up to T = (cid:126) / eV with a timestep of ∆ t = × − (cid:126) / eV, covering the range from 0 .
31 eV to3142 eV in frequency space. The large simulation box ensuredthe adequate representation of excited states. A mask wasmultiplied to the Kohn-Sham states at each time step in ordersmoothly absorb contributions above the ionization threshold.From the density variation δ n ( r , t ) = n ( r , t ) − n ( r , t = δ n (cid:96) m ( q , t ) = (cid:82) d r δ n ( r , t ) j (cid:96) ( qr ) Y (cid:96) m ( Ω r ) is then computed ineach time step and Fourier transformed to δ n (cid:96) m ( q , ω ) allow-ing to determine S ( q , ω ) as S ( q , ω ) = − I (cid:96) max (cid:88) (cid:96) = (cid:96) (cid:88) m = − (cid:96) Im (cid:2) δ n (cid:96) m ( q , ω ) (cid:3) . (3)The m -dependence is subsidiary. To a good approximationhenceforth m = ffi cient to consider | Im[ δ n (cid:96) ( q , ω )] | ≡ − Im[ δ n (cid:96), m = ( q , ω )] which stands for the (cid:96) -resolved dynamical structure factor depicted in Fig. 1. For q → molecule possessesa volume plasmon mode ( (cid:96) = (cid:96) ≥ (cid:96) ≥ ω V = (cid:112) / r s , ω , ,(cid:96) = ω (cid:2) ∓ (cid:96) + (cid:112) + (cid:96) ( (cid:96) + R / R ) (2 (cid:96) + (cid:3) . -0.002 0.0020 ω
40 50 60 7000.10.2 ℓ =1 ℓ =2 ℓ =4 ℓ =3 F S , ℓ ( / e V ) ω S1,1 = 21.59 eV ω S1,2 = 25.28 eV ω S1,3 = 26.03 eV ω S1,4 = 25.64 eV q (Å -1 ) G S , ℓ ρ S2, ℓ (a.u.) ℓ =1 ℓ =2 ℓ =3 ℓ =4 (eV) ω
40 50 60 70 (eV) ℓ =1 ℓ =2 ℓ =3 ℓ =4 (a) (c)(b) FIG. 2. (Color online) (a) Frequency-dependent part of the S1 modesobtained from the NMF (shaded curves) with fits (dashed lines). Forplasmon features we concentrate on the region ω >
18 eV. (b) q -dependent part of the S1 modes from the NMF (solid lines) alongwith fits using the model fluctuation density (symbols). (c) Modelfluctuation density ρ S1 ,(cid:96) in a plane cut through the center of themolecule for (cid:96) ranging from 1 (top) to 4 (bottom). Inspecting the (cid:96) = q ∼ . − , ω ∼
20 eV and q ∼ − , ω ∼
40 eV.As evident from Fig. 1, for higher q plasmonic modes(S1 , S2 , V) seem to merge and attain various multipoles con-tributions. This is a manifestation of electronic transitions be-tween the single-particle states with di ff erent angular momen-tum [39–41]. Thus, the question arises of how to disentanglethese modes and to unravel their multipolar nature.A suitable mathematical tool to tackle this task is thenon-negative matrix factorization (NMF), which is exten-sively used, e. g., for face recognition algorithms [27].Applied to our problem, the NMF delivers two functions F i ( ω ) ≥ G i ( q ) ≥ | δ n (cid:96) ( q , ω ) | = (cid:80) i F i ( ω ) G i ( q ) (see appendix A). This struc-ture follows namely from the Lehmann representation of χ R ( r , r (cid:48) ; ω ) as χ R ( r , r (cid:48) ; ω ) = (cid:88) α ξ α ( ω ) ρ α ( r ) ρ α ( r (cid:48) ) , ξ α ( ω ) = E α ( ω + i Γ α ) − E α (4)where ρ α is the real fluctuation density corresponding to atransition from the ground to an excited many-body state la-belled by α (with excitation energy E α ), and Γ α is the linewidth. Assuming spherical symmetry, excitations have angu-lar ( (cid:96) ) and radial ( ν ) components. Expanding ρ α ( r ) ρ α ( r (cid:48) ) = (cid:80) (cid:96) m R ν,(cid:96) ( r ) R ν,(cid:96) ( r (cid:48) ) Y (cid:96) m ( Ω r ) Y ∗ (cid:96) m ( Ω r (cid:48) ) Eq. (4) implies for the -0.014 0.0140 (b) ω S2,1 = 32.82 eV ω S2,3 = 32.15 eV ω S2,2 = 32.27 eV ω S2,4 = 31.89 eV q (Å -1 ) ρ S2, ℓ (a.u.) ℓ =1 ℓ =2 ℓ =3 ℓ =4 ℓ =1 ℓ =2 ℓ =3 ℓ =4 ω
40 50 60 70 (eV) ω
40 50 60 70 (eV) ℓ =1 ℓ =2 ℓ =3 ℓ =4 F S , ℓ ( / e V ) G S , ℓ ω S2,2 = 32.27 eV (a) (c) FIG. 3. (Color online) (a) Frequency-dependent part of the S2 modesfrom the NMF (shaded curves) and corresponding fits (dashed lines).For the latter, no constraint has been imposed on the frequency range.(b) q -dependent part of the S2 modes from the NMF (solid lines) andfits (symbols). (c) Model fluctuation density ρ S2 ,(cid:96) as in Fig. 2. structure factor S ( q , ω ) = (cid:88) ν(cid:96) (2 (cid:96) + F ν,(cid:96) ( ω ) G ν,(cid:96) ( q ) , F ν,(cid:96) ( ω ) = Im[ ξ ν,(cid:96) ( ω )] , G ν,(cid:96) ( q ) = (cid:32)(cid:90) ∞ d r r R ν,(cid:96) ( r ) j (cid:96) ( qr ) (cid:33) . In full generality the sum (4) contains infinite number of termscorresponding to the infinite number of excited states. For ho-mogeneous electron gas plasmons are strongly damped whentheir momentum enters the p - h continuum, where the non-interacting structure factor S (0) ( q , ω ) >
0. For electrons con-fined to a spherical shell the momentum can be representedby a magnitude q and an angular momentum (cid:96) . To mark thee ff ective region q max and (cid:96) max in which plasmon modes exist,we estimate the transverse momentum as 2 (cid:96)π/ R (with radius R ) and compare it to the critical momentum q crit = . k F [42] (the Fermi momentum is k F = (9 π/ / r − s ). We findso a critical (cid:96) ∼
3. Thus, any collective excitation beyond (cid:96) max = S (0) ( q , ω ) in SCA [43], for which the electrondensity enters as a central ingredient (we take the spherically-averaged DFT density n ( r )) [44]. This allows to estimate forwhich q the p - h pairs dominate the spectrum for each (cid:96) sepa-rately. For (cid:96) max = p - h domain at q (cid:38) . − .Note that due to geometrical confinement plasmons and p - h excitations intersect each other and couple so significantly.Now we separate the response into ν = S1 (Fig. 2) and ν = S2 (Fig. 3) for (cid:96) ≥
1, while the mode ν = V canbe found from (cid:96) = p - h excitations is also present in the spectra. Theplasmon frequencies ω ν,(cid:96) are identified from the maximum -0.006 0.0060 ω (eV) V2 V1 ω V2 = 24.17 eV ω V1 =42.69 eV
40 50 60 7000.02 G V F V ( / e V ) ρ V2 (a.u.) q (Å -1 ) ρ V1 (a.u.) (a)(b) (c) FIG. 4. (Color online) (a) Frequency-dependent part of the V1 (blueshaded curve) and V2 (purple shaded curve) mode from the NMFalong with corresponding fits (dashed lines). Fitting has been carriedout in the complete frequency range. (b) q -dependent part of the V1and V2 modes from the NMF (solid lines) and fits (symbols). (c)Model fluctuation densities in the same plane as in Fig. 2. of the ω -dependence spectra as obtained by the NMF in theform F fit ν,(cid:96) ( ω ) = Im[2 ω ν,(cid:96) / (( ω + i Γ ν,(cid:96) ) − ω ν,(cid:96) )]. InspectingFig. 2(a), we find the dipole plasmon at ω S1 , (cid:39) .
59 eV,this is a well established value. Increasing (cid:96) shifts the peakto larger energies (in line with the shell model); the sharppeak around 7.5 eV, which is known to consist of a seriesof p - h excitations [11], gains spectral weight until it dom-inates for (cid:96) =
4. Abundance of large angular momentumstates around HOMO-LUMO gap [41] increases the numberof channels for high-multipole electronic transitions and isresponsible for the peak’s enhancement. The plasmon fre-quency ω S1 , = .
64 eV on the other hand is smaller than ω S1 , = .
03 eV. This demonstrates the limitations of theSCA.The radial profile of the density oscillations R ν,(cid:96) ( r ) canbe inferred from G ν,(cid:96) ( q ) in that we assume R fitS1 ,(cid:96) ( r ) = A (cid:96) r exp[ − ( r − r (cid:96) ) / σ (cid:96) ] and extract the parameters ( A (cid:96) , r (cid:96) , σ (cid:96) )for which (cid:13)(cid:13)(cid:13) G S1 ,(cid:96) ( q ) − (cid:16)(cid:82) ∞ d r r R fitS1 ,(cid:96) ( r ) j (cid:96) ( qr ) (cid:17) (cid:13)(cid:13)(cid:13) is minimized.The e ff ective fluctuation densities are then given by ρ S1 ,(cid:96) ( r ) = R S1 ,(cid:96) ( r ) Y (cid:96) ( Ω r ), cf. figure 2(c).Analogous procedure for S2 modes (Fig. 3) reveals a de-crease of the plasmon energies in qualitative agreement withRef. [14]. However, the dispersion is less pronounced thanin the shell model. To characterize the fluctuation densi-ties, we use an Ansatz containing a node R fitS2 ,(cid:96) ( r ) = A (cid:96) r (1 − r / r (0) (cid:96) ) exp[ − ( r − r (cid:96) ) / σ (cid:96) ] and determine the parameters as tomatch G S2 ,(cid:96) ( q ) (Fig. 3(b)). The spatial structure of the plas-mon oscillation is shown in Fig. 3(c).A common and physically intuitive feature of the S1 andS2 modes is that the spatial extend of the fluctuation densityis growing with (cid:96) . This is a consequence of the increasing cen-trifugal force, ”pushing” the oscillation away from the center.Applying the NMF with two components to | Im[ δ n ( q , ω )] | totalexperiment
10 20 S1 S1+S2 θ =3° θ =4° θ =5° ω (eV) ω (eV) TDDFTmodel d σ / d ω d Ω ( a . u . ) d σ / d ω d Ω ( a . u . ) d σ / d ω d Ω ( a . u . ) (a) (d)(b)(c) (f)(e) FIG. 5. (Color online) TDCS for EELS of C at scattering angles θ = ◦ (a), θ = ◦ (b), and θ = ◦ (c). The energy loss ω is with respectto initial beam energy of (cid:15) = + S2and S1 + S2 + V. The thick curve shows experimental data [28]. (d–f):comparison of full TDDFT and model cross sections. shows (Fig. 4) that in addition to the expected volume plasmon(labelled by V1) around ω V1 = .
69 eV (which agrees wellwith density parameter r s ∼ ω V2 = .
17 eV appears. To clarify its origin we com-puted the response function from its non-interacting coun-terpart in the random-phase approximation and invoking theSCA (see appendix B). After obtaining | Im[ δ n ( q , ω )] | we ap-plied the NMF, as well. This procedure yields very similarspectra including the occurrence of V2. This feature is, how-ever, very sensitive to the details of the density distribution; itvanishes for a discontinuous step-like profile. Thus, it is theoscillations of the spill-out density taking place on the surfaceof the molecule that form V2. This is a pure quantum e ff ect.With the dynamical structure factor being fully character-ized, we proceed by computing the TDCS (Eq. (1)). Fig. 5compares calculated and measured [28] EELS spectra as afunction of the electron scattering angle θ which fixes themomentum transfer. The magnitudes of the measured spec-tra shown in Fig. 5 are determined up to an overfall factor.Thus, the theory-experiment comparison in Figs. 5 (b,c) is onan absolute scale. The classification of the plasmon modesaccomplished by the NMF analysis allows for plotting mode-resolved TDCS curves. As Figs. 5(a–c) demonstrate, the S1plasmons play the dominant role for small θ (which corre-sponds to the optical limit of small q ), while the S2 modesbecomes increasingly significant for larger θ (i.e., larger q ).The larger energy of the S2 with respect to the S1 plasmonsleads to the formation of a shoulder (clearly visible for θ = ◦ )and, thus, to the apparent shift of the maximum of the exper-imental EELS spectrum with growing θ . A similar e ff ect is also observed for the S1 modes due to their dispersion withrespect to (cid:96) .Furthermore, the extracted ω -dependencies and the modelfluctuation densities can be used to construct an approximatestructure factor S model ( q , ω ) = (cid:80) ν(cid:96) (2 (cid:96) + F fit ν,(cid:96) ( ω ) G fit ν,(cid:96) ( q ) thatreproduces the TDDFT results around the plasmon resonancesin a precise way by construction. Corresponding TDCSs arecompared in of Figs. 5(d–f).An important feature of the structure factor is the f -sumrule (cid:82) ∞ d ω ω S ( q , ω ) = Nq / N ).Checking for the (plasmon-dominated) S model ( q , ω ) shows thediscrepancy for larger q ; a critical value is reached when (cid:82) ∞ d ω ω S model ( q , ω ) decreases again after quadratic growth.We find q crit ∼ . − which is consistent with the estima-tion above. Hence, p - h excitations become more importantfor q > q crit and gradually diminish the plasmon contribution.In summary, we presented accurate TDDFT calculationsfor the dynamical structure factor and EELS spectra for C molecule underlining the role of higher multipole contribu-tions. Using NMF decomposition allowed to trace the evo-lution in q and ω of the symmetric and anti-symmetric sur-face and volume plasmons. In addition, we characterizedand modeled the fluctuation densities (i.e., the ingredients ofthe response function) and unveiled their multipolar character.These ingredients might, in principle, be accessed selectivelyby using electron beams carrying a definite angular momen-tum (electron vortex beams [45, 46]). By measuring the angu-lar momentum of the scattered beam the angular momentumtransfer ∆ (cid:96) becomes a control variable which the EELS spec-tra depends on. Particularly, provided the beam axis coincideswith the symmetry axis of spherical system, the plasmonicresponse upon scattering of such twisted electrons containsmultipole contributions for (cid:96) ≥ | ∆ (cid:96) | only [47]. Hence, specificmultipoles can be excluded or included by varying ∆ (cid:96) .Furthermore, we discussed the limitation of spherical-shellmodels in describing the quenching of the volume plasmonand identified the electronic density distribution as a key factordetermining its energy. We obtained excellent agreement withexperimental results and explained how the di ff erent plasmonmodes contribute to the spectra. Appendix A: Non-negative matrix factorization
As dictated by the fluctuation-dissipation theorem, theimaginary part of δ n (cid:96) ( q , ω ) for ω > | Im[ δ n (cid:96) ( q , ω )] | = − Im[ δ n (cid:96) ( q , ω )] to split (cid:12)(cid:12)(cid:12) Im (cid:2) δ n (cid:96) ( q , ω ) (cid:3)(cid:12)(cid:12)(cid:12) = N (cid:88) ν = F ν,(cid:96) ( ω ) G ν,(cid:96) ( q ) . (A1)Without imposing any restriction on the number of compo-nents ( N ) the expansion (A1) is exact and can be paralleledwith the singular value decomposition (SVD) of a general(complex or real) matrix M : M = U Σ V ∗ . The di ff erence is inthe additional requirements of positivity on the vectors form-ing U and V . The transition from continuous variables as ineq. (A1) to the matrix form is provided by discretizing the ω -and the q -points after smooth interpolation.We select N = | Im[ δ n (cid:96) ( q , ω )] | .The problem of non-negative matrix factorization can beformulated as a non-convex minimization problem for theresidue norm r = || A − WH || . Thus, the solution is notunique and may lead to local minima . Depending on the normused di ff erent algorithms can be formulated. A commonlyused method is the multiplicative update of D. Lee and S. Se-ung [27]: W ia ← W ia ( AH T ) ia ( WHH T ) ia , (A2a) H a j ← H a j ( W T A ) a j ( W T WH ) a j , (A2b)where i indexes the energy points and j numbers the timepoints. The method starts with some suitable guess for ma-trices W and H . Additionally, the vectors forming W are nor-malized each step: W ia ← W ia (cid:107) W a (cid:107) . Upon these prescriptions (A2) the Euclidean distance r monotonously decreases until the stationary point (local min-imum) has been reached. We initialized the vector W ( W )with cuts of | Im[ δ n (cid:96) m ( q , ω )] | along q direction at ω =
20 eV( ω =
40 eV), while H ( H ) is constructed by cuts at q = . − ( q = . − ). We found that typically 1000 itera-tions yield well converged results.The functions F ν,(cid:96) ( ω ) and G ν,(cid:96) ( q ) is then obtained from in-terpolating the data from H ν and W ν , respectively. We nor-malize the frequency spectra such that fitting by F fit ν,(cid:96) ( ω ) = Im[2 ω ν,(cid:96) / (( ω + i Γ ν,(cid:96) ) − ω ν,(cid:96) )] (as explained in the main text)can be performed without any additional prefactor. G ν,(cid:96) ( q ) isnormalized accordingly. This normalization procedure is con-sistent with the Lehmann representation. Appendix B: Semi-classical calculations
In order to eludicate the behavior of the volume plasmons,semi-classical calculations provide some insight. The starting point is the Dyson equation for the density-density responsefunction in random-phase approximation (RPA): χ ( r , r (cid:48) ; z ) = χ (0) ( r , r (cid:48) ; z ) + (cid:90) d r (cid:90) d r χ (0) ( r , r ; z ) × v ( r − r ) χ ( r , r (cid:48) ; z ) . (B1)We drop the superscript R and consider general complex argu-ment z here. In SCA, the non-interacting reference responsefunction χ (0) ( r , r (cid:48) ; z ) can be expressed in terms of ground-statedensity n ( r ) ( = n ( r ) as we assume spherical symmetry here)only. The subsequent derivations and the solution scheme foreq. (B1) are detailed in the Supplementary Material [48]. Theamount of spill-out density can be adjusted by varying thesmearing parameter ∆ r in the model density n ( r ) = N (cid:2) θ ∆ r ( r − R ) − θ ∆ r ( r − R ) (cid:3) ,θ ∆ r ( r ) = + exp[ − r / ∆ r ] , (B2)where R = R − ∆ R / R = R +∆ R / R = . N ensures thecorrect total valence charge. ∆ R is fixed to keep the mean den-sity constant. The scenario ∆ r → ∆ r = . ∆ r , the ( (cid:96) =
0) contributionto the structure factor, Im[ δ n ( q , z = ω + i Γ )] ( Γ = . ∆ r . While very pro-nounced for ∆ r = . ∆ r →
0. More details and graphs of volume plas-mon spectra can be found in the Supplementary Material.
ACKNOWLEDGMENTS
This work is supported by the DFG under Grants No.SFB762 and No. PA 1698 / [1] J. A. Schuller, E. S. Barnard, W. Cai, Y. C. Jun, J. S. White, andM. L. Brongersma, Nat. Mater. , 368 (2010).[2] G. Mie, Ann. Phys. , 377 (1908).[3] N. J. Halas, S. Lal, W.-S. Chang, S. Link, and P. Nordlander,Chem. Rev. , 3913 (2011).[4] L. M. Liz-Marzan, Langmuir , 32 (2006).[5] P. K. Jain and M. A. El-Sayed, Nano Lett. , 2854 (2007). [6] E. Prodan, C. Radlo ff , N. J. Halas, and P. Nordlander, Science , 419 (2003).[7] R. F. Egerton, Rep. Prog. Phys. , 016502 (2009).[8] F. J. Garc´ıa de Abajo, Rev. Mod. Phys. , 209 (2010).[9] C. Yannouleas and R. A. Broglia, Annals of Physics , 105(1992).[10] F. Alasia, R. A. Broglia, H. E. Roman, L. Serra, G. Colo, andJ. M. Pacheco, J. Phys. B , L643 (1994). [11] A. S. Moskalenko, Y. Pavlyukh, and J. Berakdar, Phys. Rev. A , 013202 (2012).[12] E. Sohmen, J. Fink, and W. Kr¨atschmer, Zeitschrift f¨ur PhysikB Condensed Matter , 87 (1992).[13] A. W. Burose, T. Dresch, and A. M. G. Ding, Zeitschrift f¨urPhysik D Atoms, Molecules and Clusters , 294 (1993).[14] P. Bolognesi, L. Avaldi, A. Ruocco, A. Verkhovtsev, A. V. Ko-rol, and A. V. Solovyov, Eur. Phys. J. D , 254 (2012).[15] I. V. Hertel, H. Steger, J. de Vries, B. Weisser, C. Menzel,B. Kamke, and W. Kamke, Phys. Rev. Lett. , 784 (1992).[16] S. W. J. Scully, E. D. Emmons, M. F. Gharaibeh, R. A. Phaneuf,A. L. D. Kilcoyne, A. S. Schlachter, S. Schippers, A. M¨uller,H. S. Chakraborty, M. E. Madjet, and J. M. Rost, Phys. Rev.Lett. , 065503 (2005).[17] A. Reink¨oster, S. Korica, G. Pr¨umper, J. Viefhaus, K. Gode-husen, O. Schwarzkopf, M. Mast, and U. Becker, J. Phys. B , 2135 (2004).[18] D. ¨Ostling, S. P. Apell, G. Mukhopadhyay, and A. Rosen, J.Phys. B , 5115 (1996).[19] B. Vasv´ari, Zeitschrift f¨ur Physik B Condensed Matter , 223(1996).[20] A. Verkhovtsev, A. V. Korol, and A. V. Solovyov, Eur. Phys. J.D , 253 (2012).[21] R. Esteban, A. G. Borisov, P. Nordlander, and J. Aizpurua,Nature Communications , 825 (2012).[22] E. Prodan and P. Nordlander, Nano Lett. , 543 (2003).[23] M. A. L. Marques and E. K. U. Gross, Annu. Rev. Phys. Chem. , 427 (2004).[24] R. Bauernschmitt, R. Ahlrichs, F. H. Hennrich, and M. M.Kappes, Journal of the American Chemical Society , 5052(1998).[25] M. E. Madjet, H. S. Chakraborty, J. M. Rost, and S. T. Manson,J. Phys. B , 105101 (2008).[26] E. Maurat, P.-A. Hervieux, and F. L´epine, J. Phys. B , 165105(2009).[27] D. D. Lee and H. S. Seung, Nature , 788 (1999).[28] A. V. Verkhovtsev, A. V. Korol, A. V. Solov’yov, P. Bolognesi,A. Ruocco, and L. Avaldi, J. Phys. B , 141002 (2012).[29] G. Giuliani and G. Vignale, Quantum Theory of the ElectronLiquid (Cambridge University Press, 2005).[30] G. Onida, L. Reining, and A. Rubio, Rev. Mod. Phys. , 601(2002).[31] M. A. L. Marques, N. T. Maitra, F. M. S. Nogueira, E. K. U.Gross, and A. Rubio, Fundamentals of Time-Dependent Den-sity Functional Theory (Springer, 2012).[32] The time-propagation method is more e ffi cient than the Casidalinear response scheme [31] for larger systems. The computa-tional cost scales cubically with the number of p - h pairs, whichis typically large for fullerenes [49].[33] M. A. L. Marques, A. Castro, G. F. Bertsch, and A. Rubio,Comp. Phys. Commun. , 60 (2003).[34] X. Andrade, J. Alberdi-Rodriguez, D. A. Strubbe, M. J. T.Oliveira, F. Nogueira, A. Castro, J. Muguerza, A. Arruabarrena,S. G. Louie, A. Aspuru-Guzik, A. Rubio, and M. A. L. Mar-ques, J. Phys. Condens. Matter , 233202 (2012). [35] X. Andrade, D. Strubbe, U. D. Giovannini, A. H. Larsen,M. J. T. Oliveira, J. Alberdi-Rodriguez, A. Varas, I. Theophilou,N. Helbig, M. J. Verstraete, L. Stella, F. Nogueira, A. Aspuru-Guzik, A. Castro, M. A. L. Marques, and A. Rubio, Phys.Chem. Chem. Phys. (2015).[36] Besides the LDA + SIC scheme, the ground-state calcula-tion was carried out using standard LDA, the generalized-gradient approximation functional LB94, and the hybrid func-tion B3LYP [50].[37] A. Sakko, A. Rubio, M. Hakala, and K. Hamalainen, J. Chem.Phys. , 174111 (2010).[38] A. Castro, M. A. L. Marques, and A. Rubio, J. Chem. Phys. , 3425 (2004).[39] M. Feng, J. Zhao, and H. Petek, Science , 359 (2008).[40] Y. Pavlyukh and J. Berakdar, Chem. Phys. Lett. , 313(2009).[41] Y. Pavlyukh and J. Berakdar, J. Chem. Phys. , 201103(2011).[42] G. Stefanucci and R. v. Leeuwen,
Nonequilibrium Many-BodyTheory of Quantum Systems: A Modern Introduction (Cam-bridge University Press, 2013).[43] Y. Pavlyukh, J. Berakdar, and K. K¨oksal, Phys. Rev. B ,195418 (2012).[44] For the non-interacting response χ R0 ( r , r (cid:48) ; ω ) function SCA, (cid:82) d r χ R0 ( r , r (cid:48) ; ω ) φ ( r ) = (1 /ω ) (cid:2) ∇ n ( r (cid:48) ) · ∇ φ ( r (cid:48) ) − n ( r (cid:48) ) ∇ φ ( r (cid:48) ) (cid:3) holds for any function φ ( r ). Setting φ ( r ) = exp[ − i q · ( r − r (cid:48) )] andassuming a spherically symmetric density n ( r ) yields for thestructure factor S (0) ( q , ω ) (cid:39) (4 /ω ) (cid:80) (cid:96) (2 (cid:96) + S (0) (cid:96) ( q , ω ), where S (0) (cid:96) ( q , ω ) = (cid:82) ∞ d r r j (cid:96) ( qr )[ q n ( r ) j (cid:96) ( qr ) + qn (cid:48) ( r ) j (cid:48) (cid:96) ( qr )].[45] M. Uchida and A. Tonomura, Nature , 737 (2010).[46] J. Verbeeck, H. Tian, and P. Schattschneider, Nature , 301(2010).[47] Employing the first Born approximation for an in-coming electron vortex beam with wave-function ψ k ⊥ , i (cid:96) i k z , i ( r ) = ζ k ⊥ , i (cid:96) i ( r , θ ) e i k z , i r cos θ e i (cid:96) i φ [51] scatteredfrom a spherical system with response function χ R ( r , r (cid:48) ; ω ) = (cid:80) (cid:96) m χ R (cid:96) ( r , r (cid:48) ; ω ) Y ∗ (cid:96) m ( Ω r ) Y (cid:96) m ( Ω r (cid:48) ) into the fi-nal state ψ k ⊥ , f (cid:96) f k z , f ( r ) yields the cross section proportional to (cid:80) (cid:96) ≥| ∆ (cid:96) | (cid:82) ∞ d r r (cid:82) ∞ d r (cid:48) r (cid:48) V (cid:96) ( r ) V ∗ (cid:96) ( r (cid:48) )Im[ χ R (cid:96) ( r , r (cid:48) ; ω )]. Explicitcalculations results in V (cid:96) ( r ) = π A (cid:96), ∆ (cid:96) (cid:90) ∞ d r (cid:48) (cid:90) π d θ (cid:48) sin θ (cid:48) ζ k ⊥ , i (cid:96) i ( r (cid:48) , θ (cid:48) ) ζ k ⊥ , f (cid:96) f ( r (cid:48) , θ (cid:48) ) × e i( k z , i − k z , f ) r (cid:48) cos θ (cid:48) ( r (cid:96)< / r (cid:96) + > ) P ∆ (cid:96)(cid:96) (cos θ (cid:48) ) . Here, Y (cid:96) m ( Ω r ) = A (cid:96), m P m (cid:96) (cos θ ) e i m φ , and ∆ (cid:96) = (cid:96) f − (cid:96) i .[48] Supplementary Material, available online.[49] G. Orlandi and F. Negri, Photochemical & Photobiological Sci-ences , 289 (2002).[50] S. F. Sousa, P. A. Fernandes, and M. J. Ramos, J. Phys. Chem.A , 10439 (2007).[51] J. Verbeeck, H. Tian, and A. B´ech´e, Ultramicroscopy113