Atomic-Scale Vibrational Mapping and Isotope Identification with Electron Beams
AAtomic-Scale Vibrational Mapping and Isotope Identification with Electron Beams
Andrea Konečná, Fadil Iyikanat, and F. Javier García de Abajo
1, 2, ∗ ICFO-Institut de Ciencies Fotoniques, The Barcelona Institute ofScience and Technology, 08860 Castelldefels (Barcelona), Spain ICREA-Institució Catalana de Recerca i Estudis Avançats,Passeig Lluís Companys 23, 08010 Barcelona, Spain
Transmission electron microscopy and spectroscopy currently enable the acquisition of spatiallyresolved spectral information from a specimen by focusing electron beams down to a sub-Ångstromspot and then analyzing the energy of the inelastically scattered electrons with few-meV energyresolution. This technique has recently been used to experimentally resolve vibrational modes in 2Dmaterials emerging at mid-infrared frequencies. Here, based on first-principles theory, we demon-strate the possibility of identifying single isotope atom impurities in a nanostructure through thetrace that they leave in the spectral and spatial characteristics of the vibrational modes. Specifi-cally, we examine a hexagonal boron nitride molecule as an example of application, in which thepresence of a single isotope impurity is revealed through dramatic changes in the electron spectra, aswell as in the space-, energy-, and momentum-resolved inelastic electron signal. We compare theseresults with conventional far-field spectroscopy, showing that electron beams offer superior spatialresolution combined with the ability to probe the complete set of vibrational modes, including thosethat are optically dark. Our study is relevant for the atomic-scale characterization of vibrationalmodes in novel materials, including a detailed mapping of isotope distributions.
CONTENTS
I. Introduction 1II. Theoretical Description of Spatially ResolvedVibrational EELS 2III. Results and discussion 4A. Valence-Electron Charge Gradients 4B. Vibrational EELS vs Infrared Spectroscopy 5C. Vibrational Mapping at the Atomic Scale 6IV. Conclusions 7APPENDIX 8EELS Probability 8Optical Response Associated with AtomicVibrations 8Optical Extinction Cross Section 9DFT Calculations 9Acknowledgments 10References 10Supplementary figures 12
I. INTRODUCTION
The ability of exciting vibrational modes in crystalsand molecules with localized probes has attracted muchattention over the last decade because of the possibility of ∗ [email protected] investigating chemical composition and atomic bondingwith high spatial resolution. Numerous theoretical andexperimental works have demonstrated that vibrationalspectroscopy is feasible with nanoscale or even atomic-scale resolution using tip-based spectroscopic techniques,such as scanning near-field optical microscopy [1–3] andtip-enhanced Raman spectroscopy [4–7]. These ap-proaches rely on the electromagnetic optical-field en-hancement produced at the probed sample area by intro-ducing sharp metallic tips, such as those that are com-monly used in atomic force and scanning tunneling mi-croscopies. However, despite the substantial efforts madein engineering the geometric properties of the probing tipto improve resolution and sensitivity [8–10], tip-based mi-croscopies are still unable to map atomic vibrations, canonly examine a fixed orientation of the specimen, andgenerally require complex analyses to subtract the unde-sired effects associated with tip-sample coupling. In ad-dition, these techniques are only sensitive to vibrationalmodes that are optically or Raman active, while dark ex-citations without a net dipole moment remain difficult todetect.Enabled by recent advances in instrumentation [11,12], electron energy-loss spectroscopy (EELS) performedin scanning transmission electron microscopes (STEMs)has emerged as an alternative, versatile technique capa-ble of mapping vibrational modes with an atomic levelof detail. In STEM-EELS, the sample is probed by abeam of fast electrons (typically with energies of 30 –300 keV) focused below 1 Å, thus allowing for the identi-fication of individual atoms. Early theoretical predictions[13–17] were followed by experimental studies demon-strating the spectral and spatial characterization of low-energy excitations (10s – 100s meV), such as phonons ingeneral, phonon polaritons in nanostructured polar crys-tals [12, 18–26], molecular vibrations [27–30], and hybrid a r X i v : . [ c ond - m a t . m t r l - s c i ] F e b modes resulting from the coupling between vibrationaland plasmonic excitations [31]. Recent achievements in-clude the detection of vibrations at truly atomic resolu-tion in hexagonal boron nitride (h-BN) [32], silicon [33],and graphene [34], as well as the visualization of a single-atom impurity in the latter.We depict common STEM-EELS experimental ar-rangements and types of probed infrared (IR) excitationsin Figure 1. Dispersion relations of vibrational modescan be retrieved under broad electron beam (e-beam) ir-radiation of extended crystals combined with angle- andenergy-resolved electron detection [22, 25] (left panel). Inaddition, localized vibrational modes (upper panel) canbe probed with atomic detail using focused e-beams [32–34], while the spatial and spectral distribution of phononpolaritons –hybrids of atomic vibrations and photons–can also be mapped in structured samples [16, 19, 20](lower panel). Interestingly, by changing the collectionconditions or the orientation of the specimen, one can fil-ter different contributions to the inelastic signal [14, 35–37], and similarly, selected excitations of specific symme-try may be addressed by shaping the electron wave func-tion (right panel) in coordination with energy and mo-mentum post-selection, as shown in a proof-of-principleexperiment at higher excitation energies for the triggereddetection of plasmons with dipolar and quadrupolar char-acter [38].Here, we present first-principles calculations demon-strating the potential of e-beams for atomic-scale vi-brational mapping, including the identification of sin-gle isotope-atom impurities. Specifically, we introducea general computational methodology based on density-functional theory that can be applied to model spatiallyresolved EELS spectra at the atomic level. Such ap-proach goes beyond the macroscopic dielectric formal-ism, which is successfully used to simulate low-loss EELSwithout atomic detail [39, 40], but fails to model the mi-croscopic characteristics of the EELS signal. We showthat e-beams are capable of exciting both bright anddark vibrational modes, while the latter are missed byfar-field IR spectroscopy, as we illustrate by comparingour results to simulated optical absorption spectra. Ad-ditionally, we discuss polarization selectivity for differente-beam orientations with respect to the sample.We note that atomically resolved vibrational mappinghas so far been demonstrated only in samples that areresistant to e-beam damage, which is not the case of or-ganic molecules, for which aloof EELS needs to be consid-ered [27–30]. However, we foresee that less energetic elec- trons together with expected improvements in the sensi-tivity of electron analyzers may eventually allow us toprobe molecular vibrations with high spatial resolution(see upper-mid panel in Figure 1). In preparation forthese advances, we analyze here a rather stable h-BN-likemolecule, which is purposely chosen as a model systemthat can naturally incorporate different boron isotopes(in particular B and B) [41], so it serves to studythe effect of a single isotope defect on the resulting high-resolution EELS maps. In particular, we demonstratethat such impurity can significantly affect the vibrationalmodes and produce clearly discernible variations in theenergy-filtered maps compared to those obtained for anisotopically pure sample. We thus conclude the feasibilityof isotope identification at the single-atom level. In ad-dition, we show that post-selection of scattered electronsdepending on the acceptance angle of the spectrometercan improve these capabilities for isotope site recognition.
II. THEORETICAL DESCRIPTION OFSPATIALLY RESOLVED VIBRATIONAL EELS
For swift electrons of kinetic energy >
30 keV interact-ing with thin specimens (<10s nm) or under aloof con-ditions ( i.e. , without actually traversing any material),coupling to each sampled mode is sufficiently weak asto be described within first-order perturbation theory,so that the EELS spectrum is contributed by electronsthat have experienced single inelastic scattering events.As a safe assumption for e-beams, the initial ( i ) andfinal ( f ) electron wave functions can be separated as ψ i | f ( r ) = (1 / √ L ) e i p i | f,z z ψ i | f ⊥ ( R ), where R = ( x, y )denotes the transverse coordinates, L is the quantiza-tion length along the e-beam direction z , the longitudi-nal wave functions are plane waves of wave vectors p i | f ,z ,and ψ i | f ⊥ ( R ) are the transverse wave functions. In ad-dition, we focus on scattering events that produce negli-gible changes in the electron momentum relative to theinitial value, such that the electron velocity vector v canbe considered to remain constant (nonrecoil approxima-tion). Finally, the atomic vibrations under study takeplace over small spatial extensions compared to the wave-length of light with the same frequency, so we adopt thequasistatic limit to describe the electron-sample interac-tion. Under these approximations, the EELS probabilityis given by [39] (see Appendix)Γ EELS ( ω ) = e π (cid:126) v X f Z d R Z d R ψ i ⊥ ( R ) ψ ∗ f ⊥ ( R ) ψ ∗ i ⊥ ( R ) ψ f ⊥ ( R ) (1) × Z ∞−∞ dz Z ∞−∞ dz e i ω ( z − z ) /v Im {− W ( r , r , ω ) } , FIG. 1. Experimental approaches to probing vibrational excitations by EELS. Left: phonon dispersion in bulk materials can bestudied through momentum-resolved EELS using an extended electron beam (e-beam) with a focal size of tens of nm [22, 25].Middle top: localized vibrational modes in finite or defective structures can be probed at atomic resolution with tightly focusede-beams [32–34]. Middle bottom: phonon polaritons –electromagnetic waves coupled to the optical phonons of ionic crystals–can efficiently be excited by e-beams to probe their spatial (nm-to-micron scale) and spectral (10s-to-100s meV) properties,which strongly depend on sample and probing geometry [16, 19, 20]. Right: pre-shaping and post-selection of the electrontransverse wave function enables the identification of specific vibrational-mode shapes and symmetries. Images are illustrationstaken from our calculations. where we sum over final transverse states f and the spec-imen enters through the screened interaction W ( r , r , ω ),defined as the potential created at r by a unit chargeplaced at r and oscillating with frequency ω . For atomicvibrations, the screened interaction reduces to a sumover the contributions of different vibrational modes, asshown in eqs A10 and A11 in the Appendix. Insertingthese expressions into eq 1 and using the identity [42] R ∞−∞ dz e i qz /r = 2 K ( | q | R ), where K is a modified Bessel function, we findΓ EELS ( ω ) = 4 e π (cid:126) v X nf Im ( | N nfi ( ω/v ) | ω n − ω ( ω + i γ ) ) , (2)where N nfi ( q ) = Z d R ψ i ⊥ ( R ) ψ ∗ f ⊥ ( R ) I n ( R , q ) ,I n ( R , q ) = X l √ M l Z d r K ( | q || R − R | ) e i qz × [ e nl · ~ρ l ( r )] , (3)the n sum runs over the sampled vibrational modes, ω n and e nl are the associated frequencies and normalizedatomic displacement vectors, respectively, the l sum ex-tends over the atoms in the structure, M l is the massof atom l , the gradient of the charge distribution withrespect to displacements of that atom is denoted ~ρ l ( r ),and we have incorporated a phenomenological dampingrate γ . As described in the Appendix, we use density-functional theory to obtain ~ρ l ( r ) and the dynamical ma- trix. The latter is then used to find the natural vibrationmode frequencies ω n and eigenvectors e nl by solving thecorresponding secular equation of motion. Because coreelectrons in the specimen are hardly affected by the swiftelectron, we assimilate them together with the nuclei topoint charges eZ l located at the equilibrium atomic po-sitions r l , so that they contribute to ~ρ l ( r ) with a term ~ρ nucl l ( r ) = eZ l ∇ r l δ ( r − r l ), which upon insertion into eq3 leads to I n ( R , q ) = I val n ( R , q ) + X l eZ l √ M l e i qz l e nl · (cid:20) | q | ( R − R l ) | R − R l | K ( | q || R − R l | ) + i qK ( | q || R − R l | ) ˆ z (cid:21) . Here, I val n ( R , q ) is computed from eq 3 by substituting ~ρ l ( r ) = ~ρ nucl l ( r ) + ~ρ val l ( r ) by the gradient of the valence-electron charge density ~ρ val l ( r ) and carrying out the inte-gral over a fine spatial grid. Incidentally, close encounterswith the atoms in the structure produce an unphysicaldivergent contribution to the loss probability, which isavoided when accounting for the maximum possible mo-mentum transfer [39] and also when averaging over thefinite width of the e-beam. For simplicity, we regularisethis divergence in this work by making the substitution | R − R | → p | R − R | + ∆ with ∆ = 0 .
16 Å in the argument of the Bessel functions of the above equations.Measurement of the inelastic electron signal in theFourier plane corresponds to a selection of transmittedelectron wave functions ψ f ⊥ ( R ) = e i Q f · R / √ A (normal-ized by using the transverse quantization area A ), hav-ing a well-defined final transverse wave vector Q f ⊥ ˆ z .Inserting this expression into eq 2 and making the sub-stitution P f → [ A/ (2 π ) ] R d Q f , we find Γ EELS ( ω ) = R d Q f [ d Γ EELS ( ω ) /d Q f ], where d Γ EELS ( ω ) d Q f = e π (cid:126) v X n (cid:12)(cid:12)(cid:12)(cid:12)Z d R ψ i ⊥ ( R )e − i Q f · R I n ( R , ω/v ) (cid:12)(cid:12)(cid:12)(cid:12) Im (cid:26) ω n − ω ( ω + i γ ) (cid:27) (4)is the momentum-resolved EELS probability. A depen-dence on the initial electron wave function is observed,which is however known to be lost when performing theintegral over the entire Q f space [39, 43], leading toΓ EELS ( ω ) = R d R | ψ i ⊥ ( R ) | Γ EELS ( R , ω ), whereΓ EELS ( R , ω ) = 4 e π (cid:126) v X n Im ( | I n ( R , ω/v ) | ω n − ω ( ω + i γ ) ) (5)( i.e. , the loss probability reduces to that experienced bya classical point electron averaged over the transverse e-beam density profile | ψ i ⊥ ( R ) | ).The above results are derived at zero temperature.When the vibrational modes are in thermal equilib-rium at a finite temperate T , the loss probabili-ties given by eqs 2, 4, and 5 need to be correctedas Γ T EELS ( ω ) = Γ EELS ( | ω | ) [ n T ( ω ) + 1] sign( ω ), where n T ( ω ) = 1 / (e (cid:126) ω/k B T −
1) is the Bose-Einsten distri-bution function and ω > ω <
0) describes elec-tron energy losses (gains). This expression, which has been used to determine phononic [21, 24] and plasmonic[44] temperatures with nanoscale precision using EELS,can be derived from first principles for bosonic modes[45], whose excitation and de-excitation probabilities areproportional to the quantum-harmonic-oscillator factors n T ( ω n ) + 1 and n T ( ω n ), respectively. In what follows,we ignore thermal corrections as a good approximationfor energy losses (cid:126) ω (cid:38)
50 meV at room temperature( k B T ∼ III. RESULTS AND DISCUSSIONA. Valence-Electron Charge Gradients
We focus on the model system sketched in Figure 2a,which consists of 7 boron (green) and 6 nitrogen (blue)atoms arranged in a h-BN-like configuration, and in-cludes 9 hydrogen edge atoms (gray) to passivate theedges and give stability to the structure [46] (see Ap-
FIG. 2. Charge density gradients in a h-BN-like molecule.(a) We show the unperturbed valence-electron charge densityintegrated over the direction z normal to the atomic planefor the molecular structure overlaid above the density plot.(b) Gradient of the z -integrated valence charge density withrespect to in-plane atomic displacements along x and y forselected B, N, and H individual atoms (identified by the rapidvariations of the gradient in these plots). The gradient ismultiplied by a different factor in each plot (see labels) tomaintain a common color scale. pendix and supplementary Figure 11). Molecules like thisone are likely produced during chemical vapor depositionbefore a continuous h-BN film is formed [46], and theycan also emerge when destroying continuous h-BN layersusing physical methods [47]. In what follows, we comparethe EELS signal from the isotopically pure molecule (tak-ing all boron atoms as B) with that obtained when oneof the boron atoms is replaced by the isotope B (darkgreen ball in Figure 2a). The isotopic composition af-fects the vibrational modes, but not the valence-electroncharge density, which is represented for the unperturbedmolecule in the underlying color plot of Figure 2a (in-tegrated over the z direction, normal to the x - y atomicplane). We observe electron accumulation around nitro-gen atoms that reflects the ionic nature of the N-B bonds,similar to what is observed in extended h-BN monolayers.Coupling to the electron probe is mediated by the chargegradients ~ρ l ( r ) entering eq 3, the valence contribution ofwhich is shown in Figure 2b after integration over the out-of-plane direction ( i.e. , R dz ~ρ val l ( r )) for the three types of atoms under consideration. Each atomic displacementleads to a dipole-like pattern centered around the dis-placed atom, and is qualitatively similar for other atomsof the same kind. As valence electrons pile up aroundnitrogen atoms, their displacements lead to the highestvalues of ~ρ val l ( r ). B. Vibrational EELS vs Infrared Spectroscopy
The gradients of the charge distribution in the nanos-tructure together with a vibrational eigenmode analy-sis provide the elements needed to calculate EELS andoptical-extinction spectra in the IR range. In Figure 3a,we show EELS spectra obtained by using eq 5 for 60 keVelectrons focused at four different positions within the x - y atomic plane. We compare results for the isotopicallypure molecule (solid curves) and the same structure withan isotopic impurity (dashed curves). The presence ofthe impurity leads to slight energy shifts of the mode en-ergies ( ∼
10s meV, see details in supplementary Figures 6and 7), as well as substantial changes in the correspond-ing spectral weights. In addition, we observe that mostfeatures are enhanced when the e-beam is moved closer tothe atomic positions, while several of them persist evenwhen the beam passes outside the molecule ( cf. blueand green spectra). Such persisting features are associ-ated with modes that exhibit a net dipole moment, sothey can also be revealed through far-field optical spec-troscopy, as shown in the extinction cross sections plottedin Figure 3c (see Appendix for details of the calculation).Dipole-active modes couple to the e-beam over long dis-tances, and consequently, they can be detected in thealoof configuration, which has been recently employedin EELS experiments to study beam-sensitive molecules[27–30]. Indeed, we corroborate a strong resemblanceof the optical extinction (Figure 3c) and the aloof EELS(Figure 3a, green curves) spectra. Incidentally, the three-fold symmetry of the molecules renders the optical crosssection independent of the orientation of the polarizationvector within the plane of the atoms. The most intenseEELS peaks are observed at around 130 meV and 170-180 meV, corresponding to modes involving B-N bondvibrations (see supplementary Figures 6 and 7 for de-tails of the atomic displacement vectors associated withdifferent vibrational modes). Weaker features at around320 meV (see supplementary Figure 10) arise from N-Hbond stretches. Interestingly, the presence of the isotopeimpurity only has a small effect on dipole-active modes( i.e. , a small reduction of intensity and a weak energyshift). In contrast, some of the dark modes probed byEELS are severely affected.By rotating the nanoflake with respect to the e-beam direction (or analogously, the light polarization foroptical measuremensts), a different set of modes con-tributes to the spectral features, dominated by out-of-plane atomic motion under the conditions of Figure 3b,d.In EELS, the strengths of these features strongly depend
FIG. 3. Effect of an isotope impurity in EELS and infrared (IR) extinction. We compare calculated EELS (a,b) and opticalextinction (c,d) spectra for an isotopically pure molecule (solid curves, only B atoms) and a molecule containing a singleisotope atom impurity (dashed curves, B impurity marked by black arrows in the insets). For EELS (a,b), we take 60 keVelectrons and consider different e-beam positions (see color-coordinated vertically offset curves) and orientations (normal tothe plane of the atoms in (a); parallel to and 1 Å away from that plane in (b)). For optical IR spectroscopy (c,d), we considerin-plane (c) and out-of-plane (d) light polarization (see red arrows in the insets). All vibrational mode energies obtained forthe isotopically pure structure are indicated by vertical gray lines. We incorporate a spectral broadening of 1 meV to accountfor intrinsic losses and instrument resolution. on e-beam position and orientation. In particular, thespectra shown in Figure 3b (with the electrons passingparallel to and 1 Å away from the plane of the atoms) re-veal low-energy peaks that are absent under normal inci-dence (Figure 3a), associated with optical modes that arealso missed in the optical extinction for out-of-plane po-larization (Figure 3d), which confirms their dark nature,in contrast to the excitations observed around ∼
100 meVenergy.
C. Vibrational Mapping at the Atomic Scale
To visualize the complete spatial dependence of thevibrational EELS signal at specific energies correspond-ing to selected modes, we calculate energy-filtered mapsby scanning the beam position over an area covering thestudied structure (placed in the x - y plane). In Figure4, we show maps calculated for a beam of 0.6 Å fwhmfocal size and selected modes in isotopically pure and de-fective molecules (see supplementary Figures 8 and 9 for additional maps). By inspecting the results for ∼
44 meV(Figure 4a) and ∼
179 meV (Figure 4a) energy losses, to-gether with the corresponding atomic displacement vec-tors (two degenerate modes at each of these energies), weobserve a strong correlation in symmetry and strengthbetween the mode displacements and the EELS maps,thus corroborating that the latter provide a solid basisto reconstruct the contribution of each atom to the vi-brational modes. In addition, by introducing a boronisotope impurity (indicated by black arrows), the three-fold symmetry of the mode displacements and the result-ing EELS maps are severely distorted. The impurity canproduce either depletion or enhancement of the EELS in-tensity around its position, but in general, it influencesthe inelastic electron signal over the entire area of in-terest. However, for the B-H stretching modes around320 meV, which are nearly decoupled from vibrations indifferent parts of the molecule, only the area close to theimpurity is affected.In Figure 5, we explore the effect of post-section onthe resulting EELS intensity, as calculated for transverse-
FIG. 4. Isotope-impurity determination through energy-filtered EELS imaging. We plot the calculated EELS probabilityintegrated over different energy-loss ranges (see labels) as a function of e-beam position under the conditions of Figure 3a forthe isotopically pure molecule (left plots in each panel) and a molecule with an isotope atom defect (right plots, B atomicimpurity indicated by black arrows). In (a,b), we show the atomic displacements associated with the involved vibrationalmodes. The probability is multiplied by a different factor in each density plot (see lower-right labels) to maintain a commoncolor scale. Frame colors are coordinated with the vertical bands indicated in Figure 3a. We spatially average the probabilityover a transverse Gaussian e-beam profile of 0.6 Å fwhm. momentum-resolved scattered electrons using eq 4. Thee-beam is taken to have a Gaussian density profile | ψ i ( R ) | of 0.6 Å fwhm. We present calculations fortwo of the most prominent features in EELS at ener-gies around ∼
44 meV and ∼
179 meV, and compareresults from isotopically pure (upper row) and defec-tive (lower row) samples. The obtained momentum-dependent energy-filtered maps clearly demonstrate thatthe directly transmitted electrons ( Q f = 0) carry long-range information associated with the polarization ofvalence-electron charges (see supplementary Figures 8and 9, showing the separate contributions of valence-electron and nuclear charges). Such long-range signal isparticularly prominent for the ∼
179 meV modes, whichare dipole-active (see IR spectra in Figure 3b). In con-trast, the modes around ∼
44 meV are dark, and thus, themap filtered at low momentum transfers reveals a weakersignal, except when symmetry is broken by the presenceof the isotope impurity, which introduces a nonzero netdipole moment. In general, when collecting electronsexperiencing larger transverse momentum transfers, thespatial localization of the signal is increased, enabling abetter determination of the atomic positions and their contributions to the observed modes with a resolutionlimited by the finite size of the e-beam spot. Inciden-tally, we can link the energy-integrated low- and large-angle signal to bright- and dark-field images, respectively,as collected with different apertures in STEMs [32, 36].
IV. CONCLUSIONS
STEM-EELS has evolved into a leading technique en-abling both atomic-resolution imaging and spectroscopyover a broad frequency range extending down to the mid-IR. In particular, atomic-scale mapping of IR vibrationalmodes, which has recently become accessible [32–34], isimportant for understanding and manipulating the op-tical response in such spectral range, as well as for de-termining the effect of phonons and phonon-polaritonsin the electrical and thermal conductivities of nanostruc-tured materials. We have shown that STEM-EELS canprobe the complete set of vibrational modes, includingoptically bright and dark modes, the polarization char-acteristics of which can be addressed by tilting the sam-ple relative to the e-beam or by resolving the inelastic
FIG. 5. Isotope-impurity determination through momentum-resolved EELS. Under the conditions of Figure 3a, we plot thecalculated EELS probability as a function of e-beam position for the isotopically pure molecule (upper plots) and a moleculewith an isotope atom defect (lower plots, B atomic impurity indicated by black arrows) after integration over both the specifiedenergy-loss range (around ∼
44 meV and ∼
179 meV, respectively) and the transverse wave vector of the transmitted electronswithin a narrow annular aperture centered around the e-beam direction. Different radii of the latter ( Q f ) are considered ineach column, as indicated by the lower labels. The probability is multiplied by a different factor in each color plot (see labels)to maintain a common color scale. The incident e-beam has a transverse Gaussian profile of 0.6 Å fwhm. electron signal in scattering angle. The latter can benefitfrom new advances in hybrid pixel detection technology[48]. As an application of these methods, we have demon-strated through first-principles calculations that an indi-vidual isotope impurity in a h-BN-like molecule producesradical changes in the spectral and spatial characteristicsof the electron signal associated with the excitation of itsvibrational modes. Our study supports the use of STEM-EELS to map isotope distributions with atomic precision,which can be important for understanding phononic life-times, as well as thermal and electrical transport at thenanoscale. APPENDIXEELS Probability
Under the assumptions discussed in the main text, theloss rate reduces to [39]Γ
EELS ( ω ) dt = 2 e (cid:126) X f Z d r Z d r ψ i ( r ) ψ ∗ f ( r ) ψ ∗ i ( r ) ψ f ( r ) × Im {− W ( r , r , ω ) } δ ( ε f − ε i + ω ) , where ψ i ( r ) and ψ f ( r ) denote initial and final electronwave functions of energies (cid:126) ε i and (cid:126) ε f , respectively.We now separate longitudinal and transverse compo-nents as specified in the main text ( i.e. , ψ i | f ( r )) =(1 / √ L ) e i p i | f,z z ψ i | f ⊥ ( R )), use the nonrecoil approxima-tion to write ε f − ε i ≈ ( q f,z − q i,z ) v , transform thesum over final states through the prescription P f → ( L/ π ) R dq f,z P f ( i.e. , the remaining sum over f nowrefers to transverse degrees of freedom), carry out the q f,z integral using the δ function, and multiply the result by the interaction time L/v to convert the loss rate intoa probability. Following these steps, we readily find eq 1.
Optical Response Associated with AtomicVibrations
The quasistatic limit is adopted here under the as-sumption that the studied structures are small com-pared with the light wavelength at the involved oscil-lation frequencies. We consider a perturbation poten-tial φ ext ( r , t ) due to externally incident light or a swiftelectron, in response to which the atoms in the struc-ture (labeled by l = 1 , · · · , N ) oscillate around theirequilibrium positions r l with time-dependent displace-ments u l ( t ). The charge density ρ ( { u } , r ), which ob-viously depends on { u } ≡ { u , · · · , u N } , is calculatedfrom first principles as discussed below. Following astandard procedure to describe atomic vibrations [49],we Taylor-expand the configuration energy (also calcu-lated from first principles) for small displacements andonly retain the lowest-order u -dependent contribution(1 / P ll u l ( t ) · D ll · u l ( t ), where D ll is the so-calleddynamical matrix. We now write the Lagrangian of thesystem as L = (1 / P l M l | ˙ u l ( t ) | − (1 / P ll u l ( t ) ·D ll · u l ( t ) − R d r ρ ( { u } , r ) φ ext ( r , t ), where M l is the mass ofatom l and the rightmost term accounts for the potentialenergy in the presence of the perturbing electric poten-tial φ ext ( r , t ). The equation of motion then follows from ∂ t ∇ ˙ u l L = ∇ u l L , which leads to M l ¨ u l ( t ) = − X l D ll · u l ( t ) − Z d r ~ρ l ( r ) φ ext ( r , t ) , (A6)where we have used the property D T ll = D l l , while thevector field ~ρ l ( r ) = ∇ u l ρ ( { u } , r ) represents the gradientof the electric charge density with respect to displace-ments of atom l . Linear response is assumed by using thedynamical matrix. In addition, to be consistent with thisapproximation, we evaluate ~ρ l ( r ) at the equilibrium posi-tion u = 0. It is then convenient to treat each frequencycomponent separately to deal with the corresponding ex-ternal potential φ ext ( r , ω ) = R ∞−∞ dt e i ωt φ ext ( r , t ), so thateq A6 becomes ω ( ω + i γ ) M l u l ( ω ) = X l D ll · u l ( ω )+ Z d r ~ρ l ( r ) φ ext ( r , ω ) , (A7)where we have introduced a phenomenological dampingrate γ . The solution to eq A7 can be found by first con-sidering the symmetric eigenvalue problem for the freeoscillations, X l √ M l M l D ll · e nl = ω n e nl , (A8)where n labels the resulting vibration modes of fre-quencies ω n . The eigenvectors e nl form a complete( P n e ∗ nl ⊗ e nl = δ ll I , where I is the 3 × P l e ∗ nl · e n l = δ nn ) basis set, whichwe use to write the atomic displacements as u l ( ω ) = 1 √ M l X n c n ( ω ) e nl with expansion coefficients c n ( ω ) = 1 ω ( ω + i γ ) − ω n × X l √ M l e ∗ nl · Z d r ~ρ l ( r ) φ ext ( r , ω ) . Incidentally, because the dynamical matrix is real andsymmetric, the eigenvectors e nl can be chosen to bereal, but we consider a more general formulation usingcomplex eigenvectors that are convenient to describe ex-tended crystals, where the mode index n may naturallyincorporate a well-defined Bloch momentum. Now, intro-ducing the above expression for u l in the induced chargedensity ρ ind ( r , ω ) = P l u l · ~ρ l ( r ), we can obtain the sus-ceptibility χ ( r , r , ω ) = X nll √ M l M l [ e nl · ~ρ l ( r )] [ e ∗ nl · ~ρ l ( r )] ω ( ω + i γ ) − ω n , (A9)which is implicitly defined by the relation ρ ind ( r , ω ) = R d r χ ( r , r , ω ) φ ext ( r , ω ). Finally,the screened interaction, defined by W ( r , r , ω ) = R d r R d r χ ( r , r , ω ) | r − r | − | r − r | − , reducesto W ( r , r , ω ) = X n S n ( r ) S ∗ n ( r ) ω ( ω + i γ ) − ω n , (A10) where S n ( r ) = X l √ M l Z d r e nl · ~ρ l ( r ) | r − r | . (A11)Equations A10 and A11 are used to evaluate eq 1 in themain text and produce eqs 2, 4, and 5. Optical Extinction Cross Section
For a small structure such as that considered in Fig-ure 2, we describe the far-field response in terms of thepolarizability tensor, which we in turn calculate from thesusceptibility by writing the dipole induced by a unit elec-tric field, that is, ¯¯ α ( ω ) = − R d r R d r χ ( r , r , ω ) r ⊗ r .Using eq A9, we find¯¯ α ( ω ) = 2 (cid:126) X n ω n d n ⊗ d ∗ n ω n − ω ( ω + i γ ) , where d n = X l r (cid:126) ω n M l Z d r [ e nl · ~ρ l ( r )] r (A12)is the transition dipole associated with mode n . Wethen calculate the optical extinction cross section as [50] σ ext ( ω ) = (4 πω/c )Im { α ii ( ω ) } for light polarization alongeither an in-plane ( i = x ) or an out-of-plane ( i = z )symmetry direction. In practice, as described in themain text, we separate the charge gradient ~ρ l ( r ) = ~ρ nucl l ( r ) + ~ρ val l ( r ) into the contributions of valence elec-trons ( ~ρ val l ( r )) and the rest of the system ( ~ρ nucl l ( r ), thatis, nuclei and core electrons), so the transition dipoles re-duce to d n = d val n + e P l Z l p (cid:126) / (2 ω n M l ) e nl , where d val n is calculated from eq A12 by replacing ~ρ l ( r ) by ~ρ val l ( r ). DFT Calculations
We perform density functional theory (DFT) cal-culations using the projector-augmented-wave (PAW)method [51] as implemented in the Vienna ab initio sim-ulation package (VASP) [52–54] within the generalizedgradient approximation in the Perdew-Burke-Ernzerhof(PBE) form to describe electron exchange and correla-tion [55]. The cut-off energy for the plane waves is setto 500 eV. The atomic positions of the molecular struc-tures under consideration with and without an isotopicimpurity are determined by minimizing total energies andatomic forces by means of the conjugate gradient method.Atomic positions are allowed to relax until the atomicforces are less than 0.02 eV/Å and the total energy dif-ference between sequential steps in the iteration is be-low 10 − eV. The edges of the molecule are passivatedwith hydrogen terminations to maintain structural sta-bility, leading to a nearly hexagonal structure (see supple-mentary Figure 11). Additionally, hydrogen passivation0eliminates a strong effect associated with the danglingbonds observed in the vibrational spectra. We calculatethe dynamical matrix D ll using the small displacementmethod: each atom in the unit cell is initially displacedby 0.01 Å and the resulting interatomic force constantsare then determined to fill the corresponding matrix ele-ments. Vibrational eigenmodes and eigenfrequencies areobtained by diagonalizing the dynamical matrix, as dis-cussed above. We assimilate nuclear and core-electroncharges in each atom to a point charge. However, thedistribution of the valence-electron charge density is in-corporated using a dense grid in the unit cell to tabu-late the charge gradients ~ρ val l ( r ) from the change in thevalence electron density produced for each small atomicdisplacement. ACKNOWLEDGMENTS
This work has been supported in part by the EuropeanResearch Council (Advanced Grant 789104-eNANO), theEuropean Commission (Horizon 2020 Grants 101017720FET-Proactive EBEAM and 964591-SMART-electron),the Spanish MINECO (MAT2017-88492-R and SeveroOchoa CEX2019-000910-S), the Catalan CERCA Pro-gram, and Fundaciós Cellex and Mir-Puig.
REFERENCES [1] R. Hillenbrand, T. Taubner, and F. Keilmann, Nature , 159 (2002).[2] F. Huth, A. Govyadinov, S. Amarie, W. Nuansing,F. Keilmann, and R. Hillenbrand, Nano Lett. , 3973(2012).[3] I. Amenabar, S. Poly, W. Nuansing, E. H. Hubrich,A. A. Govyadinov, F. Huth, R. Krutokhvostov, L. Zhang,M. Knez, J. Heberle, et al., Nat. Commun. , 1 (2013).[4] K. Kneipp, Y. Wang, H. Kneipp, L. T. Perelman,I. Itzkan, R. R. Dasari, and M. S. Feld, Phys. Rev. Lett. , 1667 (1997).[5] R. Zhang, Y. Zhang, Z. Dong, S. Jiang, C. Zhang,L. Chen, L. Zhang, Y. Liao, J. Aizpurua, Y. Luo, et al.,Nature , 82 (2013).[6] J. Lee, K. T. Crampton, N. Tallarida, and V. A. Apkar-ian, Nature , 78 (2019).[7] R. B. Jaculbia, H. Imada, K. Miwa, T. Iwasa, M. Take-naka, B. Yang, E. Kazuma, N. Hayazawa, T. Taketsugu,and Y. Kim, Nat. Nanotech. , 105 (2020).[8] B.-S. Yeo, W. Zhang, C. Vannier, and R. Zenobi, Appl.Spectrosc. , 1142 (2006).[9] M. Asghari-Khiavi, B. R. Wood, P. Hojati-Talemi,A. Downes, D. McNaughton, and A. Mechler, J. RamanSpectrosc. , 173 (2012).[10] S. Mastel, A. A. Govyadinov, C. Maissen, A. Chuvilin,A. Berger, and R. Hillenbrand, ACS Photonics , 3372(2018).[11] O. L. Krivanek, J. P. Ursin, N. J. Bacon, G. J. Corbin,N. Dellby, P. Hrncirik, M. F. Murfitt, C. S. Own, andZ. S. Szilagyi, Philos. Trans. Royal Soc. A , 3683(2009). [12] O. L. Krivanek, T. C. Lovejoy, N. Dellby, T. Aoki, R. W.Carpenter, P. Rez, E. Soignard, J. Zhu, P. E. Batson,M. J. Lagos, et al., Nature , 209 (2014).[13] P. Rez, Micros. Microanal. , 671 (2014).[14] C. Dwyer, Phys. Rev. B , 054103 (2014).[15] J. R. M. Saavedra and F. J. G. de Abajo, Phys. Rev. B , 115449 (2015).[16] H. Lourenço-Martins and M. Kociak, Phys. Rev. X ,041059 (2017).[17] A. Konečná, T. Neuman, J. Aizpurua, and R. Hillen-brand, ACS Nano , 4775 (2018).[18] C. Dwyer, T. Aoki, P. Rez, S. L. Y. Chang, T. C. Lovejoy,and O. L. Krivanek, Phys. Rev. Lett. , 256101 (2016).[19] M. J. Lagos, A. Trügler, U. Hohenester, and P. E. Batson,Nature , 529 (2017).[20] A. A. Govyadinov, A. Konečná, A. Chuvilin, S. Vélez,I. Dolado, A. Y. Nikitin, S. Lopatin, F. Casanova, L. E.Hueso, J. Aizpurua, et al., Nat. Commun. , 1 (2017).[21] J. C. Idrobo, A. R. Lupini, T. Feng, R. R. Unocic, F. S.Walden, D. S. Gardiner, T. C. Lovejoy, N. Dellby, S. T.Pantelides, and O. L. Krivanek, Phys. Rev. Lett. (2018).[22] F. S. Hage, R. J. Nicholls, J. R. Yates, D. G. McCulloch,T. C. Lovejoy, N. Dellby, O. L. Krivanek, K. Refson, andQ. M. Ramasse, Sci. Adv. , eaar7495 (2018).[23] A. Konečná, K. Venkatraman, K. March, P. A. Crozier,R. Hillenbrand, P. Rez, and J. Aizpurua, Phys. Rev. B , 205409 (2018).[24] M. J. Lagos and P. E. Batson, Nano Lett. , 4556(2018).[25] R. Senga, K. Suenaga, P. Barone, S. Morishita, F. Mauri,and T. Pichler, Nature pp. 247–250 (2019).[26] R. Qi, R. Wang, Y. Li, Y. Sun, S. Chen, B. Han, N. Li,Q. Zhang, X. Liu, D. Yu, et al., Nano Lett. , 5070(2019).[27] P. Rez, T. Aoki, K. March, D. Gur, O. L. Krivanek,N. Dellby, T. C. Lovejoy, S. G. Wolf, and H. Cohen, Nat.Commun. , 10945 (2016).[28] D. M. Haiber and P. A. Crozier, ACS Nano , 5463(2018).[29] J. R. Jokisaari, J. A. Hachtel, X. Hu, A. Mukherjee,C. Wang, A. Konecna, T. C. Lovejoy, N. Dellby, J. Aizpu-rua, O. L. Krivanek, et al., Adv. Mater. , 1802702(2018).[30] J. A. Hachtel, J. Huang, I. Popovs, S. Jansone-Popova,J. K. Keum, J. Jakowski, T. C. Lovejoy, N. Dellby, O. L.Krivanek, and J. C. Idrobo, Science , 525 (2019).[31] L. H. G. Tizei, V. Mkhitaryan, H. Lourenço-Martins,L. Scarabelli, K. Watanabe, T. Taniguchi, M. Tencé,J. D. Blazit, X. Li, A. Gloter, et al., Nano Lett. ,2973 (2020).[32] F. S. Hage, D. M. Kepaptsoglou, Q. M. Ramasse, andL. J. Allen, Phys. Rev. Lett. , 016103 (2019).[33] K. Venkatraman, B. D. Levin, K. March, P. Rez, andP. A. Crozier, Nat. Phys. , 1237 (2019).[34] F. S. Hage, G. Radtke, D. M. Kepaptsoglou, M. Lazzeri,and Q. M. Ramasse, Science , 1124 (2020).[35] A. Asenjo-Garcia and F. J. García de Abajo, Phys. Rev.Lett. , 066102 (2014).[36] P. M. Zeiger and J. Rusz, Phys. Rev. Lett. , 025501(2020).[37] P. Rez and A. Singh, Ultramicroscopy , 113162(2021). [38] G. Guzzinati, A. Beche, H. Lourenco-Martins, J. Martin,M. Kociak, and J. Verbeeck, Nat. Commun. , 14999(2017).[39] F. J. García de Abajo, Rev. Mod. Phys. , 209 (2010).[40] G. Radtke, D. Taverna, M. Lazzeri, and E. Balan, Phys.Rev. Lett. , 027402 (2017).[41] A. J. Giles, S. Dai, I. Vurgaftman, T. Hoffman, S. Liu,L. Lindsay, C. T. Ellis, N. Assefa, I. Chatzakis, T. L.Reinecke, et al., Nat. Mater. , 134 (2018).[42] I. S. Gradshteyn and I. M. Ryzhik, Table of Integrals,Series, and Products (Academic Press, London, 2007).[43] R. H. Ritchie and A. Howie, Philos. Mag. A , 753(1988).[44] V. Mkhitaryan, K. March, E. Tseng, X. Li, L. Scara-belli, L. M. Liz-Marzán, S.-Y. Chen, L. H. G. Tizei,O. Stéphan, J.-M. Song, et al., p. arXiv:2011.13410(2020).[45] F. J. García de Abajo and V. Di Giulio, p.arXiv:2010.13510 (2020).[46] M. Maruyama and S. Okada, Sci. Rep. , 16657 (2018).[47] P. Valerius, C. Herbig, M. Will, M. A. Arman, J. Knud-sen, V. Caciuc, N. Atodiresei, and T. Michely, Phys. Rev.B , 235410 (2017).[48] B. Plotkin-Swing, G. J. Corbin, S. De Carlo, N. Dellby,C. Hoermann, M. V. Hoffman, T. C. Lovejoy, C. E.Meyer, A. Mittelberger, R. Pantelic, et al., Ultrami-croscopy , 113067 (2020).[49] N. W. Ashcroft and N. D. Mermin, Solid State Physics (Harcourt College Publishers, Philadelphia, 1976).[50] H. C. van de Hulst,
Light Scattering by Small Particles (Dover, New York, 1981).[51] P. E. Blöchl, Phys. Rev. B , 17953 (1994).[52] G. Kresse and J. Hafner, Phys. Rev. B , 558 (1993).[53] G. Kresse and J. Furthmüller, Phys. Rev. B , 11169(1996).[54] G. Kresse and J. Furthmüller, Comput. Mater. Sci. , 15(1996).[55] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev.Lett. , 3865 (1996). SUPPLEMENTARY FIGURES
FIG. 6. Plots showing all of the eigenmodes and eigenvectors for the isotopically pure h-BN-like molecule in Figure 2a. Opencircles denote the positions of B (green), N (blue), and H (gray) atoms, whereas the arrows show the magnitude and directionof the atomic discplacement vectors projected on the plane of the atoms. Mode energies are indicated above the plots, with(quasi-)degenerate ones sharing the same label color. We have excluded the 9 lowest energy modes, 6 of which emerge fromrigid translations or rotations of the molecule. FIG. 7. Same as Fig. 6, but for the molecule containing a single B isotope impurity, represented by a dark green circles (seealso black arrow in the upper-left panel and Figure 2a). FIG. 8. Energy-filtered maps calculated as a function of the position of the focused electron beam by including differentcontributions to the polarization of the molecule associated with its vibrational modes (i.e., using the notation of the mainpaper, we consider partial contributions to ~ρ l ( r ) = ~ρ nucl l ( r )+ ~ρ val l ( r ), where "nucl" refers to the sum of nuclear and core-electroncharges, while "val" indicates the contribution of valence electrons). Different columns show calculations performed includingvalence-electron charges (valence), nuclear and core-electron charges (nuclear), and the sum of the two of them (total), asindicated by the upper labels. The three columns on the left (right) correspond to the isotopically pure (defective) molecule(see Figures 6 and 7). Each map is obtained by integrating the EELS probability over the indicated energy range. FIG. 9. Continuation of Figure 8. FIG. 10. Same as Figure 3a,c in the main text, but plotted over a wider energy range. (a)(b)(a)(b)