Advances in modelling electron energy loss spectra from first principles

Electron energy loss (EEL) spectroscopy carried out within a (scanning) transmission electron microscope can provide chemical and bonding information with atomic resolution. The information that lies within the spectrum can be difficult to extract, and often reference spectra are used to identify atomic bonding environments. First principles simulations are able to relate features in spectra to atomistic models and are particularly important in the interpretation of spectra where there are no appropriate bulk references, such as those from nanomaterials or interfaces. This paper reviews the recent developments in first principles simulations of EEL spectra and highlights the potential for advances in our understanding of materials.


Introduction
When a beam of electrons passes through the sample in a (scanning) transmission electron microscope ((S)TEM), it can excite, amongst other things, phonons, plasmons or core-electrons within the material. Once the beam has interacted with the material, it can be dispersed to produce an electron energy loss, or EEL, spectrum. Different amounts of energy are required to excite phonons (up to hundreds of meV), plasmons (up to about 50 eV) and core-electrons (10-1000 s of eV), and so peaks corresponding to the different inelastic events appear in different parts of the EEL spectrum.
The approach taken to modelling an EEL spectrum depends on several factors: the part of the spectrum being studied (the phonon, plasmon or core-loss region), the material being studied, including if it is amorphous, molecular or crystalline (and the elements contained within it), and the specific edge being investigated. In an EEL experiment, we are probing the loss function [1], L(q, ω), which is related to the dielectric function, ϵ(q, ω), by L(q, ω) = Im[−1/ϵ(q, ω)] 1 . For modelling the plasmon region, the full dielectric function needs to be modelled, whilst when modelling the core-loss region it is possible to approximate the loss function by the imaginary part of the dielectric function [1]. Historically, rather than a full calculation of the EEL spectrum, core-loss edges have been interpreted using the local density of states, LDOS. This is still a common method for understanding the electronic properties of nano-materials. In addition, even in cases where the full spectrum is used to interpret experiment, LDOS can be an important tool for relating peaks in the spectrum to the structural model and electronic structure.
Several codes exist which can calculate EEL spectra from first principles including WIEN2k [7,8], FEFF [9,10], CASTEP [11,12], Elk [13], Yambo [14,15], ONETEP [16,17] and MADNESS [18]. Different codes use different approximations and are suitable in different situations. An excellent description of the mechanics of modelling an EEL spectrum is given in Chapter 14 of Transmission Electron microscopy: Diffraction, Imaging, and Spectrometry [19] and readers are referred there (and to references within it) for more details about how to model EEL spectra.
No matter the spectrum being simulated, modelling EEL spectra relies on constructing a model and then running a computer code based on a theoretical treatment of the excitation process. There are three important ingredients: (a) The model of the system; (b) The theoretical treatment of the excitation process; (c) The computational power available.
The review aims to cover: (a) Improvements / developments in structure prediction techniques, including why it is important and what the effect on our ability to model spectra is. (b) New theoretical developments: The need for theory to interpret the phonon region of the EELS spectrum and the development of new functionals for density functional theory (DFT) which have the potential to allow us to model a wider range of materials. (c) Ways of making larger systems more tractable.

High-resolution EEL spectroscopy
Whilst this review is mainly concerned with EEL spectroscopy carried out within the (S)TEM, it is useful to mention the very closely related technique, reflection, or high-resolution, EEL spectroscopy (HREELS). In HREELS, the beam of electrons (of energy 2-20 eV [20] rather than the 10-100 s of keV energies used for (S)TEM) is reflected from the sample. The technique is used to investigate surfaces, often with molecules absorbed onto them [21], and has been used for looking at phonons [21,22], plasmons [22] and magnons [23]. Details of the technique can be found in [21,22]. Recent advances in modelling HREELS data from first principles includes work showing the importance of including both the elastic and inelastic signals [24] and using a quantum-mechanical treatment of the probe electron [25]. In addition, Nazarov has discussed the correct form of the energy loss function for both reflection and transmission EEL spectroscopy for quasi-2D materials and developed an ab initio approach for modelling these materials using periodic codes [26].

Comparison of experimental and simulated data
The comparison of experiment and theory is not necessarily straightforward. It is usually done in a qualitative manor [27] and there is a hierarchy of comparison: Are there the correct number of peaks? Are the peaks in the right places? Are the intensities correct? Some authors have tried a more quantitative method for comparing different approximations and experiment, such as the χ 2 method employed by McCulloch et al [27]. Often, however, theory does not reproduce experiment exactly. There are a large variety of reasons for this, for example the level of approximation in the theory being used may not be adequate or the experimental data may contain more than one bonding environment. Often what is more important than reproducing exactly the experimental data is understanding the origin of a difference between spectra or relating a peak to a structural feature. Understanding something about the material or sample is generally the aim of the modelling.

Structure prediction
The first step in simulating an EEL spectrum is to construct a model of the sample you are interested in. If the material is well known and well characterised, there are databases of crystal structures, such as the Inorganic Crystallographic Structural Database, ICSD [28,29], which can provide atomic co-ordinates. Much characterisation, however, is done on materials about which less is known; often that is the motivation behind the research. It then becomes necessary to find an alternative approach to constructing a model. Sometime it is possible to use our chemical knowledge about similar materials to construct viable structures. These can then be compared using DFT to see which structure is most energetically favourable.
In the past decade computational methods have been developed to discover new materials [30,31]. It is now possible to systematically use existing knowledge and databases in data mining approaches, which rely less on imagination and intuition [32,33]. In addition, structure prediction from first principles which needs little or no pre-existing knowledge is now also possible [30]. There are several different approaches to first principles structure prediction, including random sampling [34][35][36], evolutionary algorithms (see references in [30,37]) and minima hopping [38,39]. Experimental data can help structure prediction techniques by providing constraints or data to compare with. Structure prediction is not confined to bulk structures, but can be used for interface [40,41] and defect structure prediction [42][43][44]. Structure prediction is becoming more utilised in materials science [45], and with initiatives such as the Materials Project [46,47], the trend is likely to continue. The Materials Project itself includes an x-ray absorption module which uses FEFF9 [48] to simulate x-ray absorption spectra, which are very similar to EEL spectra. TEM and EEL spectroscopy can provide constraints as well as information to test structure prediction against. Constraints can include, for example, composition or symmetry information. Imaging and spectroscopy can be used to compare model and experimental data as a way of validating a model structure which is important in the case where the structure of a material is unknown [49].
Work combining EEL spectroscopy and an ab initio random structure searching (AIRSS) algorithm developed by Pickard and Needs [34,50] on ZrO show the power of combining structure prediction and electron microscopy [51]. Zirconium based alloys are used for fuel cladding and structural components in water-cooled nuclear reactors. Much effort has gone into optimising these materials to withstand the harsh environment and the lifetime of the zirconium cladding is limited by oxidation caused by water corrosion [52]. The oxidation process had been linked to a sub-oxide phase which forms at the interface between Zr metal and ZrO 2 [53] and, though the composition of the phase had been determined using EEL spectroscopy, the structure was unknown. Low-loss EEL spectra taken across the interface showed distinct signals from the metal, sub-oxide and oxide [51]. To use modelling to understand the structure of the ZrO sub-oxide it was first necessary to produce candidate structures and this was done using AIRSS. For a series of randomly generated structures with the correct stoichiometry, geometry optimisations were carried out using DFT to find a local minima in the energy landscape. All optimised structures were then ranked according to energy to produce a convex hull, as shown in figure 1. One of the most interesting aspects of this work was that the AIRSS produced two structures which could not be distinguished energetically using DFT as the difference in their formation energies was 0.002 eV per formula unit. EEL spectra and diffraction patterns were simulated from the two structures and compared with experiment to determine the structure of the ZrO. Only one of the structures predicted from first principles was able to explain all the experimental data. Both structure prediction and experimental EEL spectra were necessary to determine the structure which is now routinely identified (for example by EBSD [54]) at the grain boundary.
The first important ingredient in modelling EEL spectra is the model of the system. Being able to predict structures from first principles gives us a much wider variety of structures to consider and releases us from the constraints of our own imagination.

Developments in theoretical descriptions of the excitation process
EEL spectra span a large energy range and, as mentioned in the introduction, includes different types of energy loss events including phonons, plasmons and core-electron excitations. Although the loss function contains information about all of these different excitations, the different energy ranges and dominating physics often mean that different approaches are taken to model spectra from the different types of excitation.

Phonon spectroscopy in the transmission electron microscope
Phonons are a materials specific property and phonon dispersions have been probed experimentally using neutrons since the 1950s [55]. The theory of lattice dynamics is well developed [56] and first principles calculations can be used to calculate phonon dispersions [57,58]. The challenge in interpreting different types of experimental data is knowing which modes contribute to the spectrum, and with what intensity.
Electrons have been used to probe phonons in materials since the 1960s using HREELS [21] (see section 1.1). Recent developments in hardware have allowed phonon vibrations to be seen in the TEM using EEL spectroscopy [59]. Theoretical treatments developed for HREELS are not applicable to transmission EEL spectroscopy due to the geometrical constraints present in reflection EEL spectroscopy. Advances in theory of the excitation of phonons in transmission EEL spectroscopy have allowed first-principles calculations to aid interpretation of the experimental data.
For phonon scattering in transmission EEL spectroscopy, two different regimes have been identified by Dwyer et al [60]: dipole scattering and localised vibrational scattering. Dipole scattering involves small momentum transfer and can be carried out in aloof mode, i.e. with the beam next to, but not traversing, the sample. Several approaches to modelling spectra in the dipole regime have been developed [60][61][62][63][64]. In one of the earliest works, Rez [61] looked at the theoretical response of a diatomic molecule and predicted that the signal associated with phonons would be of a similar strength to that from core-loss excitations. Saavedra and García de Abajo [62] have used a Hamiltonian formalism to describe the interaction between the fast (beam) electron and an atomic cluster. They use a classical external potential to describe the electron and a second-quantisation formalism for the phonons. They applied the theory to graphene clusters of varying size and made the link between the spectrum and density of vibrational states.
Radtke et al [63] have used a classical approach to relate the EELS spectrum obtained in an aloof geometry to the low-frequency dielectric tensor. By using DFT to calculate the dielectric tensor they were able to interpret data from a guanine crystal, see figure 2. Radtke et al [65] have also shown how the symmetry of the phonon modes can be retrieved from orientation dependent aloof beam experiments on anisotropic B 12 P 2 crystals. A semi-classical approach was also used by Hohenester et al to calculate the q = 0 spectra obtained from MgO nanocubes [64,66]. A classical response theory has been used by Govyadinov et al [67] to show the existence of hyperbolic phonon polaritons in hexagonal boron nitride, which is important for the interpretation of aloof spectra for strongly anisotropic materials.
In terms of the localised vibrational scattering regime, Forbes and Allen [69] have used a transition potential approach and looked at the spatial effect of the beam geometry. Nicholls et al have extended the scattering theory for neutrons and x-rays to include electrons used it to interpret momentum resolved measurements from cubic and hexagonal boron nitride [70,71], see figure 3. A similar approach has been used by Senga et al to compare insulating hexagonal boron nitride and its semi-metal carbon counterpart graphite [72]. Hage et al [73] have used the phonon density of states to interpret the vibrational signal associated with a single impurity atom in graphene. They show the importance of constructing a large supercell to reduce artificial localisation for defects.
Transmission phonon EEL spectroscopy is an exciting new development in the exploration of lattice dynamics and, to make the most of this new experimental technique, interpretation of spectra using the newly developed theoretical frameworks will be key.

Advances in DFT functionals
DFT provides the basis of many EEL spectra calculations and the quality of the predicted spectrum is a reflection of the ability of DFT to model a particular material. One of the main approximations in DFT is the treatment of the exchange-correlation term and this is the main limitation in the prediction of ground state structures by DFT [74]. The earliest and simplest treatment of this term was to use the local density approximation (LDA) when the potential is a function of the local density only [75]. The generalised gradient approximation (GGA) also includes information about the gradient of the charge density, but, unlike LDA, there is no single universal form [75]. The Perdew-Burke-Ernzerhof [76] functional is most widely used form of GGA functional and is currently the 'standard' for a DFT calculation.
It is possible to go beyond this and include the second derivative of the electron density in what is known as a meta-GGA functional and many types of meta-GGA have been developed [77]. In addition to the semi-local functionals, there are also non-local and hybrid functionals (which include an amount of Hartree-Fock exchange energy) and the variety can be daunting. Many of the meta-GGAs have been developed empirically (based on the energies, rather than electron densities) for a specific class of materials. As a result, the functionals are not generally good for a wide variety of materials. Medvedev et al [78] have highlighted the dangers of empirical fitting for functionals. SCAN [79,80] is a recently developed meta-GGA which includes orbital kinetic energy density terms and satisfies all known constraints for a semi-local functional. Because it has not been fitted to a particular system it can predict diverse types of bonding, including intermediate van der Waals interactions [81]. This is a general functional and so far has been shown to describe materials at least as well as GGA [82,83]. It has been used for a variety of materials including metal surfaces [84], defects in semiconductors [81], ice [81] and ferroelectrics [85]. It is the first functional to be able to describe cuperate superconductors without any free parameters [86]. This functional has the potential to allow the simulation of a larger variety of systems without the requirement of using system-specific functionals, making them more accessible to modelling of EEL spectra.

Approaches to making larger calculations more tractable
EEL spectra simulations can get larger computationally for two main reasons: one is that the number of atoms, i.e. the size of the chunk of material being modelled, increases. This could be due to decreasing core-hole-core-hole interactions, looking at a defect, or simply looking at a material with a large unit cell. The second way calculations can get more computationally intensive is to increase the complexity of the calculation. A standard DFT calculation calculates the ground state of a system. The very nature of doing EEL spectroscopy involves looking energy lost to excitations within the system-the system is no longer in the ground state. Standard DFT calculations of core-loss spectra carried out with a core hole (an electron missing from the core state of the excited atom) are often useful for interpreting spectra [19]. For some materials, it is necessary to use using methods that go beyond standard DFT [6] such as the Bethe-Salpheter Whilst it is possible to conceive and construct ever larger and more complex models and simulations, what can actually be calculated is restricted by the amount of computational power available. Improvements in hardware and software mean that calculations are continually able to become faster and more complex. The increase in the availability of computing power is one of the driving forces behind the development of structure prediction methods discussed above.
As more atoms are included in a DFT calculation, it gets larger in terms of the computational resources it requires. Different approaches have been used to try and decrease the computational expense of calculations. Several authors have developed approaches to combine all electron and pseudopotential methods. In all-electron methods, all the electrons in the atom are included in the calculation, whilst in pseudopotential methods the core electrons (which do not take part in bonding) are replaced with a more slowly-varying pseudopotential. As such, pseudopotentials reduce the computational cost of a simulation. Donval et al [87] have developed a hybrid approach for light elements where the corehole effect can be neglected by combining WIEN2k and VASP [88]. They have applied this method to study the silicon L 2,3 edge. MADNESS is a molecular DFT code [18] which is able to assign atoms as all-electron or pseudopotential within one calculation. It is able to calculate core-edge spectra, including a core-hole [89]. In addition, MADNESS also uses a wavelet basis set, which is computationally efficient [18] and makes it ideal for modelling non-periodic systems.
In general the size of a calculation scales with the cube of the number of atoms. Linear scaling DFT methods have been developed which scale linearly with the number of atoms involved. One such code, ONETEP [16,17], can calculate up to 3-4000 atoms, and is also able to calculate EEL spectra [90]. Potential applications of such large structures include large supercells to minimise core-hole-core-hole interactions and calculating the spectra of defect atoms, surface atoms or structures such as interfaces and grain boundaries. ONETEP has already been used to investigate anatase surfaces and defects [90] (note: in these examples wavefunction truncation was applied, but kernal truncation was not as the systems were not large enough for it to be advantageous).
Computationally, when using DFT, one of the difficulties with running core-loss calculations is that you are normally interested in the energy range 20-30 eV above the edge onset. The number of states needed to correctly describe this energy range increases as you increase the size of the unit cell leading to a possibly worse than cubic scaling with number of atoms. Core-loss EELS calculations can still be computationally demanding, i.e. at the limit of what is reasonably possible, and whilst that means that currently some calculations are too big to run, we can look forward to improvements in the future.

Conclusions
This review paper has highlighted the potential for impact for several recent advances in first principles simulations of EEL spectra. Advances in structure prediction can open up the field of candidate structures for materials where the structure is not known. Combining structure prediction with EEL simulations can be very powerful and it is likely to be seen more often in the future. Advances in microscope hardware which have allowed phonons to be detected within the TEM have been the catalyst for advances in the theory to describe the spectra obtained. Already we are seeing the importance of having theory to explain the experimental spectra. In addition, the meta-GGA functional SCAN has the potential to allow a wider variety of materials to be studied using the combination of EEL spectroscopy and modelling.
It is still the case that EEL calculations can be computationally demanding. Going beyond DFT is important for a variety of materials, and there are still areas in which there is a lot of opportunity for advancing and improving simulations. Each advance in modelling opens up the variety of problems that can be tackled and the materials and phenomena that can be understood, and that is the driving force behind advancing simulation.