Nanoscale imaging of phonon dynamics by electron microscopy – Nature

Fabrication of SiGe QDs and preparation of TEM sample

The QD sample was grown at 600 °C in an ultra-high vacuum chemical vapour deposition (UHV/CVD) system. For Ge and Si depositions, pure GeH₄ and SiH₄ gases were used as precursors, respectively. The Si wafers were first etched in a diluted HF solution to create a hydrogen-passivated surface, before deposition. After deposition of a 50 nm thick buffer layer of Si, Ge layers were grown with 20 nm Si spacer layers in between each Ge QD layer for the formation of self-assembled conventional QDs. For investigating thermal stability and tuneability of the structural parameters, in situ post-deposition annealing was conducted at the growth temperature for 1 h. The final product was a thin-film-like material of 40-period multifold QD stacks with a thickness as high as approximately 1.2 µm. QD nanostructures spanned 70–90 nm across and were 6–8 nm thick (see Extended Data Fig. 2). The SiGe QDs have a diamond-like structure similar to that of Si and Ge, with a random and disordered arrangement of Si and Ge atoms. More details of the growth process can be found in ref. ³⁰. The cross-sectional TEM sample used in this study was prepared by focused-ion-beam milling, whereas the planar-view TEM sample was prepared by mechanical polishing. The QD interfaces were normal to the (001) crystallographic plane, whereas the zone axis was along the [110] crystallographic direction. Of all the QDs investigated, several were elongated by up to 10–15% (Extended Data Fig. 2g, h) in no particularly consistent direction. However, there were some quite symmetric SiGe QDs (Extended Data Fig. 2i).

Raman spectroscopy

Raman spectra were acquired using a Renishaw inVia confocal Raman microscope. Point spectra in the QD sample and Si wafer were acquired in 50 one second frames and then summed to achieve high signal-to-noise ratio.

STEM imaging

Images were taken with the Nion UltraSTEM 200 high-energy resolution monochromated EELS system (HERMES) operating at 60 kV and 200 kV acceleration voltages with 33 mrad and 34 mrad convergence semi-angles, producing 1.5 Å and 0.78 Å sized probes, respectively. The high-angle annular dark-field signals were collected using a high-angle annular detector with inner and outer collection angles of 103 mrad and 210 mrad, respectively. Beam current was approximately 100 pA in the imaging condition. The sizes of the QD interfaces were estimated using the 10–90% criterion³⁸.

Composition EELS mapping

2D planar and cross-sectional elemental maps were acquired using the double Cs-corrected 300 kV JEOL Grand ARM S/TEM system. A 0.5 eV per pixel dispersion was used for elemental mapping, with a convergence semi-angle of 25 mrad and a collection angle of 100 mrad. Spectra were acquired in the 500–2,355 eV range, including Ge L edge at 1,217 eV and Si K edge at 1,839 eV (Extended Data Fig. 2c). The core-loss spectrum in the interlayer Si shows only one peak, corresponding to the Si K edge, whereas the spectrum in the SiGe QD shows two peaks, corresponding to the Ge L edge and Si K edge. Single QD elemental maps were acquired with a 0.5 nm per pixel step size, whereas low-magnification spectra were acquired with a 2 nm step size and were used to determine the composition distribution in different QD layers. The core-loss signals of Si and Ge were background subtracted and integrated using Gatan’s Digital Micrograph software. The histogram data reported in Extended Data Fig. 2a were obtained by averaging 4 × 4 pixel centres of each QD.

HERMES experiments

Experimental data were acquired on a Nion UltraSTEM 200 microscope with the HERMES system operating at 60 kV, to achieve a balance between high spatial resolution and high inelastic scattering probabilities by electrons for vibrational spectroscopy. EEL spectra were produced by low-angle scattered electrons (Extended Data Fig. 1a). Combining the monochromator with a high-dispersion spectrometer, we achieved the best energy resolution of 5.7 meV with the probe in vacuum for a 10 ms acquisition at 60 keV. A CMOS camera collected a 2D spectrum image, which was then cropped and flattened in the undispersed direction to produce a 1D EELS spectrum.

Momentum-averaged 33 mrad condition

EELS acquisition was performed using a 33 mrad high-current probe with resolution and beam current as described in the section ‘STEM imaging’. A high probe current is necessary for attaining high signal-to-noise ratio spectra, at the cost of slightly decreasing the spatial resolution. Insertion of the monochromating slit reduces the probe current to about 3 pA, but gives an energy width of 8.3 meV for an acquisition time of 1 second with the probe placed on the sample. This enables accurate probing of the optical modes of germanium. An EELS entrance aperture of 1 mm subtends a 25 mrad EELS collection semi-angle. With the BZ of silicon having a Γ–X length of only 8.96 mrad, a collection semi-angle of 25 mrad collects inelastically scattered electrons from multiple neighbouring BZs as well.

Spectral data were acquired in the form of line scans and 2D maps with the monochromating slit inserted. For line scans, 100 spectral frames were collected at 1 second exposure, aligned by their respective zero-loss peak (ZLP) maximums and then summed for each probe position. The dark current spectrum was frequently obtained during experiments. Aligning individual frames by ZLP averages for any random noise present in each frame. Maximizing single-exposure time without saturating the EELS camera provided the best signal-to-noise ratio. Implementation of this acquisition scheme was achieved using Nion Swift software and custom Python scripts designed to directly control the necessary hardware parameters.

Phonon maps of an entire SiGe QD consist of nanometre-step 80 × 15 arrays of data points that were obtained with 1–1.5 second exposure and 5 frames per pixel. Due to the large size of the map, the number of integrated frames was optimized to limit sample drift and drops in emission current overtime. A typical acquisition time for maps was about 2–3 h and presented several challenges for acquiring high quality data: sample contamination, drift, gradual drop in beam current and deterioration of energy resolution due to low-order spectrometer aberrations. To mitigate these effects, the following precautions and procedures were undertaken:

1. 1.

   The instrument was tuned to optimal EELS acquisition parameters achieving an energy resolution of about 8.3 meV on the sample and was stabilized for at least 4 h. The experiment was then performed the morning after and, owing to the remarkable stability of the UltraSTEM, a sample drift of less than 5 Å per hour was present.

2. 2.

   Custom mapping code enabled the flashing of the tip every hour to ensure a high beam current, which was necessary for optimal signal-to-noise ratio as well as the automatic correction of EELS first-order aberrations every few points, to maintain a high-energy resolution throughout the experiment.

Momentum-resolved 3 mrad condition

With a 3 mrad convergence semi-angle and 2.5 mrad EELS collection semi-angle, the real-space probe size was estimated to be 2.6 nm and on-axis energy resolution to be 8.3 meV in vacuum and 9.1 meV on the sample. DPM mapping was carried out by collecting off-axis spectral intensity within the first Brillion Zone. Off-axis collection was achieved by using a combination of post-specimen lenses to shift the diffraction image while keeping the real-space probe stationary, so that the region of interest in reciprocal space lies on the EELS entrance aperture (Extended Data Fig. 9). Energy resolution in this condition was estimated to be 12–14 meV, making peak separation difficult. Acquired maps were 30 × 20 in size with a step size of 1 nm per pixel and 5 s exposure. The sample was stabilized to minimize drift so that data from both areas symmetric about the Γ could be acquired from the same QD with less than 1 nm of drift.

Spectral data processing

All acquired spectra were processed using custom Python codes. Spectra were first normalized with respect to the ZLP maximum so that each vibrational intensity represented a fractional scattering probability. Spectra were then background subtracted using a linear combination of a Lorentzian centred on the ZLP and an exponential polynomial e^p(x), where p(x) is an even fourth-order polynomial. This function was fitted to energy windows before and after the regions of interest. This combination produced a sharply decreasing left end and a slow and stable right tail, which were necessary for accurately extracting the low-energy optical modes of germanium. The fit was obtained using scipy.optimize, a Python library and fitting coefficients, and covariances were extracted. To determine the efficacy of the background fit, error values obtained from the square root of the diagonal terms of the covariance matrix, were examined and minimized. The background-subtracted signals were scrutinized for any negative, non-physical values.

By fitting the vibrational EELS signal with pseudo-Voight fits, which are linear combinations of Gaussians and Lorentzians, we can gather information about a mode’s excitation probability and its energy position, given by the fitting intensity and peak-position outputs. Individual inelastic probabilities were separated from background-subtracted spectra by performing pseudo-Voight fits using curve_fit() from scipy.optimize. Fitting parameters of individual peaks and their corresponding errors were extracted. We examined the value of the errors obtained from the covariance of the fit to validate our peak separation. Due to the nature of the convoluted signal, phonon dispersion curves and Raman scattering data were used to accurately deconvolute the overlapping peaks. Lower bounds of 0 and upper bounds of infinity and 40 meV were placed on the height and width of the pseudo-Voight fits, respectively. The energy position bounds were chosen to be ±5 meV from their initial reference position to avoid crossing of the separate pseudo-Voight fits.

Momentum-resolved spectra were first trimmed and processed by binned principal component analysis (PCA) before background subtraction due to the low phonon scattering cross-section under this beam geometry. After background subtraction, due to the poor energy resolution and low intensity, peak separation was not feasible, and a simple integration was carried out instead to obtain the Si optical mode intensity. The 2D maps were created by integrating the signal in the 55–65 meV energy region in each pixel.

Principal component analysis

Given the acquisition parameters of the line scans, the excellent signal-to-noise ratio was adequate to accurately separate the individual modes. The low acquisition parameters for mapping datasets required the use of PCA to improve the signal quality. Rather than using Fourier filtering or other smoothening techniques, PCA learns from the large map dataset and reconstructs the spectra, maintaining important features while improving their quality.

Map data were first background subtracted and trimmed to reduce feature size. The dataset was then arranged as N × D, where N represents background-subtracted spectra at different points and D the pixel intensity of each individual signal. Three principal eigenvectors described most of the data and including additional eigenvectors offered only a fraction of a percent increase in cumulative explained variance. PCA was used on map sizes of 80 × 15, producing a total of N = 1,200 data points and trimmed background-subtracted spectra of only D = 182 pixels, achieving an ideal condition for this machine-learning algorithm to produce excellent results (Extended Data Fig. 10).

Off-axis momentum-resolved data could not be smoothed the same way due to the smaller collection angle and modified phonon-scattering cross-section. In the standard implementation of PCA, not every sample is included and therefore it fails to smooth the data as the dominant feature is noise. Instead, we created a superset of the map data by binning the pixels and adding them to the original map dataset. Binning averages out the noise while enhancing the spectral features of interest. First, map data were binned by 2 then 3, etc., with each successively added to the original dataset to create a superset. With the addition of binned data back into the dataset, PCA performs as intended and accurately smooths the data. The same criteria as above were used in selecting the proper number of eigenvectors.

Data visualization

Contour plots of mapping data were created using matplotlib.pyplot, another Python library. Extracted parameters from the pseudo-Voigt fitting were then used to create a spatially resolved 2D map. Maps were generated using a Gaussian interpolation for better visualization. Phonon-mode intensity maps and peak-position maps were constructed for silicon optical modes in units of normalized intensity and meV, respectively.

Multimodal simulation approaches

We used three complementary computational methods (first-principles calculation, Green’s function and BTE) to simulate the non-equilibrium contributions from the beam and interface reflections. First, to understand the experimental vibrational EEL spectra from interlayer Si and SiGe QDs, first-principles calculations were used to carry out a precise simulation of phonon band structure and partial DOS (PDOS) of Si, Ge and disordered SiGe alloy (see the next section and Extended Data Figs. 1 and 4). First-principles calculations were also performed to consider the strain effect in interlayer Si (see Supplementary Section 3 and Supplementary Fig. 2). The weak tensile strain in the Si near the abrupt interface will slightly decrease the PDOS of Si OM, which cannot explain the enhancement of EELS intensity. Second, to further understand the enhancement of Si OM intensity in Fig. 3, we used Green’s function to investigate the effect of the QD geometry in a large supercell containing both gradual and abrupt interfaces (see Supplementary Section 2 and Supplementary Fig. 1), which cannot be handled by the first-principles calculations. However, the local phonon DOS results do not show a notable enhancement of PDOS near the interface, to interpret the change of vibrational signal of optical phonons in Fig. 3. The EELS cross-section is intimately related to the vibrational properties of atoms, including the phonon distribution function and phonon dispersion, whereas the non-equilibrium phonon part of the phonon population is closely related to the reflection coefficients (see equations (1) and (2) below). Finally, we narrowed down the source of the enhancement to the non-equilibrium distribution of phonons generated due to the energy transfer from fast electron to the sample. We used the BTE (see Supplementary Section 6, Fig. 3c, Supplementary Fig. 7, as well as Extended Data Fig. 6) to compute the non-equilibrium phonon population. In BTE simulations, the mode-resolved phonon reflection coefficients by the abrupt and gradual interfaces are obtained from atomistic Green’s function (see Supplementary Section 5 and Supplementary Figs. 5 and 6). The phonon Boltzmann transport describes the phonon dynamics, and we adopt the relaxation time approximation in solving the phonon BTE, which has been shown to be a good approximation in SiGe alloys³⁹.

Phonon dispersion and DOS simulations

Our first-principles calculations were carried out using the Vienna ab-initio simulation package with the projector augmented wave method. This method was adopted for the interaction between valence electrons and ionic cores^40,41, where the energy cut-off for the plane-wave basis expansion was set to 700 eV. The generalized-gradient approximation with the functional developed by Perdew–Burke–Ernzerhof was chosen for the exchange-correlation functional⁴². All atoms were fully relaxed using the conjugated gradient method for the energy minimization until the force on each atom became smaller than 0.01 eV Å⁻¹. The phonon spectrum and the corresponding phonon DOS were obtained using density functional perturbation theory⁴³. To compare with the experimental phonon signals, the phonon DOS was convolved with a Gaussian of width 7 meV to match the energy resolution of the ZLP.

DPM mapping

Well-established momentum-resolved approaches can measure phonon dispersion curves along certain reciprocal directions and reveal the local spectral variation induced by exotic phonon modes. However, such approaches cannot identify the direction of phonon propagation, which is essential for understanding the heat transport in real devices. To obtain the propagation direction of specific phonons, we developed a DMP method by subtracting phonon states at opposite reciprocal locations. Although state-of-the-art methods in momentum-resolved vibrational EELS have had great success in measuring phonon dispersion curves²⁶, and studying nanoscale defect modes²⁵, they have not yet taken advantage of the momentum polarization selection that becomes available in this configuration.

To recover directional and momentum information and elucidate phonon dynamics at the QD interfaces, a momentum-resolved beam-detector geometry was used (Fig. 4a) to obtain a DPM map. A convergence semi-angle of 3 mrad (momentum resolution of 0.5 Å⁻¹) was used here, which was about one-third the size of the Si BZ, as shown in Fig. 4b, enabling momentum resolution, whereas the 33 mrad area encompasses even neighbouring BZs (Extended Data Fig. 9a). Integrated spectral intensity in the 55–65 meV region (Si OM) in the red and blue regions represents Si OM phonons created by electrons deflected towards the 002 (Δ₊) and (00bar{2},) (Δ_–) regions in reciprocal space, respectively (Extended Data Fig. 9b, c).

A theoretical description of DPM imaging begins with the EELS scattering cross-section, which takes into account the non-equilibrium phonons by replacing the equilibrium phonon occupation ({n}_{{boldsymbol{q}},nu }) with the complete phonon population f_q,ν = n_q,ν + g_q,ν, where g_q,ν denotes only the non-equilibrium contribution²⁶:

Here, a₀ is the Bohr radius, ħ is reduced Planck’s constant, q is the momentum transfer, ({omega }_{{boldsymbol{q}},nu }) is the phonon frequency for momentum transfer q and phonon branch (nu ), ({M}_{I}) and ({{boldsymbol{tau }}}_{i}) are the atomic mass and position, ({F}_{I}({boldsymbol{q}})) is the component in the Fourier transform of the charge density associated with atom (I) and ({{boldsymbol{e}}}_{{boldsymbol{q}},nu }^{I}) represents the eigenvector for atom (I) with mass ({M}_{I}).

Then, we perform integration over the energy of the Si OM and over the physical aperture size and location in momentum space. We replace (nu ) with Si OM as well. We then obtain an expression for the EELS phonon intensity in the DMP spectra:

The key aspect of this expression is the differential across the diametrically opposed aperture locations and the non-equilibrium phonon population. The specularity of interfaces creates an anisotropy in the non-equilibrium phonon population, which provides the contrast in the DPM map.

When compared to the on-axis, momentum-averaged mapping in Fig. 3a, the intensity of the Si OM mode is higher in the QD for the Δ₊ region than the interlayer Si (Extended Data Fig. 9b), despite a lower Si concentration, suggesting that there is a stronger preference for phonons to have a downward momentum in the QD. The same intensity enhancement as seen in Fig. 3a is recovered when summing intensity from both Δ₊ and Δ_– regions (Extended Data Fig. 9d).

Inference of phonon MFP using reflection intensities

The reflection-induced non-equilibrium phonon population g′ decays as a function of distance from the interface, with the decay length being the MFP of the Si OM (Extended Data Fig. 6). A 20 nm × 15 nm mapping of the interlayer Si bounded by two QDs illustrates a gradual change in the Si OM intensity (Extended Data Fig. 7). The top of the region is bounded by an abrupt interface, whereas the bottom is bounded by a gradual interface. With the consideration that both interfaces have some degree of specularity, the data were fit with the sum of the two exponential decays arising from both the interfaces. With the fit, we obtained an MFP (Λ) value of 50.2 ± 19.4 nm for the Si OM, which is within an order of magnitude of accepted values, being in the range of tens of nanometres⁴⁴, demonstrating that the MFP of phonons can be measured at nanometre resolution.