Fabrication of SiGe QDs and preparation of TEM sample
The QD sample was grown at 600 °C in an ultrahigh vacuum chemical vapour deposition (UHV/CVD) system. For Ge and Si depositions, pure GeH_{4} and SiH_{4} gases were used as precursors, respectively. The Si wafers were first etched in a diluted HF solution to create a hydrogenpassivated 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 selfassembled conventional QDs. For investigating thermal stability and tuneability of the structural parameters, in situ postdeposition annealing was conducted at the growth temperature for 1 h. The final product was a thinfilmlike material of 40period 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 diamondlike 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. ^{30}. The crosssectional TEM sample used in this study was prepared by focusedionbeam milling, whereas the planarview 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 signaltonoise ratio.
STEM imaging
Images were taken with the Nion UltraSTEM 200 highenergy resolution monochromated EELS system (HERMES) operating at 60 kV and 200 kV acceleration voltages with 33 mrad and 34 mrad convergence semiangles, producing 1.5 Å and 0.78 Å sized probes, respectively. The highangle annular darkfield signals were collected using a highangle 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^{38}.
Composition EELS mapping
2D planar and crosssectional elemental maps were acquired using the double Cscorrected 300 kV JEOL Grand ARM S/TEM system. A 0.5 eV per pixel dispersion was used for elemental mapping, with a convergence semiangle 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 coreloss 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 lowmagnification spectra were acquired with a 2 nm step size and were used to determine the composition distribution in different QD layers. The coreloss 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 lowangle scattered electrons (Extended Data Fig. 1a). Combining the monochromator with a highdispersion 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.
Momentumaveraged 33 mrad condition
EELS acquisition was performed using a 33 mrad highcurrent probe with resolution and beam current as described in the section ‘STEM imaging’. A high probe current is necessary for attaining high signaltonoise 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 semiangle. With the BZ of silicon having a Γ–X length of only 8.96 mrad, a collection semiangle 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 zeroloss 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 singleexposure time without saturating the EELS camera provided the best signaltonoise 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 nanometrestep 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 loworder spectrometer aberrations. To mitigate these effects, the following precautions and procedures were undertaken:

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.
Custom mapping code enabled the flashing of the tip every hour to ensure a high beam current, which was necessary for optimal signaltonoise ratio as well as the automatic correction of EELS firstorder aberrations every few points, to maintain a highenergy resolution throughout the experiment.
Momentumresolved 3 mrad condition
With a 3 mrad convergence semiangle and 2.5 mrad EELS collection semiangle, the realspace probe size was estimated to be 2.6 nm and onaxis energy resolution to be 8.3 meV in vacuum and 9.1 meV on the sample. DPM mapping was carried out by collecting offaxis spectral intensity within the first Brillion Zone. Offaxis collection was achieved by using a combination of postspecimen lenses to shift the diffraction image while keeping the realspace 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 fourthorder 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 lowenergy 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 backgroundsubtracted signals were scrutinized for any negative, nonphysical values.
By fitting the vibrational EELS signal with pseudoVoight 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 peakposition outputs. Individual inelastic probabilities were separated from backgroundsubtracted spectra by performing pseudoVoight 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 pseudoVoight fits, respectively. The energy position bounds were chosen to be ±5 meV from their initial reference position to avoid crossing of the separate pseudoVoight fits.
Momentumresolved spectra were first trimmed and processed by binned principal component analysis (PCA) before background subtraction due to the low phonon scattering crosssection 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 signaltonoise 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 backgroundsubtracted 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 backgroundsubtracted spectra of only D = 182 pixels, achieving an ideal condition for this machinelearning algorithm to produce excellent results (Extended Data Fig. 10).
Offaxis momentumresolved data could not be smoothed the same way due to the smaller collection angle and modified phononscattering crosssection. 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 pseudoVoigt fitting were then used to create a spatially resolved 2D map. Maps were generated using a Gaussian interpolation for better visualization. Phononmode intensity maps and peakposition maps were constructed for silicon optical modes in units of normalized intensity and meV, respectively.
Multimodal simulation approaches
We used three complementary computational methods (firstprinciples calculation, Green’s function and BTE) to simulate the nonequilibrium contributions from the beam and interface reflections. First, to understand the experimental vibrational EEL spectra from interlayer Si and SiGe QDs, firstprinciples 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). Firstprinciples 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 firstprinciples 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 crosssection is intimately related to the vibrational properties of atoms, including the phonon distribution function and phonon dispersion, whereas the nonequilibrium 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 nonequilibrium 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 nonequilibrium phonon population. In BTE simulations, the moderesolved 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^{39}.
Phonon dispersion and DOS simulations
Our firstprinciples calculations were carried out using the Vienna abinitio 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 cutoff for the planewave basis expansion was set to 700 eV. The generalizedgradient approximation with the functional developed by Perdew–Burke–Ernzerhof was chosen for the exchangecorrelation functional^{42}. 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 Å^{−1}. The phonon spectrum and the corresponding phonon DOS were obtained using density functional perturbation theory^{43}. 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
Wellestablished momentumresolved 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 stateoftheart methods in momentumresolved vibrational EELS have had great success in measuring phonon dispersion curves^{26}, and studying nanoscale defect modes^{25}, 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 momentumresolved beamdetector geometry was used (Fig. 4a) to obtain a DPM map. A convergence semiangle of 3 mrad (momentum resolution of 0.5 Å^{−1}) was used here, which was about onethird 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 crosssection, which takes into account the nonequilibrium 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 nonequilibrium contribution^{26}:
$$frac{{{rm{d}}}^{2}sigma ({boldsymbol{q}},omega )}{{rm{d}}Omega {rm{d}}omega }=frac{4}{{{a}_{0}}^{2}}frac{hslash }{{q}^{4}}{sum }_{v}frac{{f}_{{boldsymbol{q}},nu }+1}{{omega }_{{boldsymbol{q}},nu }}{{sum }_{I}frac{1}{sqrt{{M}_{I}}}{F}_{I}({boldsymbol{q}}){boldsymbol{q}}cdot {{boldsymbol{e}}}_{{boldsymbol{q}},nu }^{I}{e}^{i{boldsymbol{q}}cdot {{boldsymbol{tau }}}_{I}}}^{2}delta (omega {omega }_{{boldsymbol{q}},nu }).$$
(1)
Here, a_{0} 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:
$${I}_{{rm{D}}{rm{M}}{rm{P}}}=frac{4hslash }{{{a}_{0}}^{2}}{int }_{{omega }_{{rm{S}}{rm{i}}{rm{O}}{rm{M}}}}{rm{d}}omega ({int }_{{Delta }_{+}}{int }_{{Delta }_{}})frac{{rm{d}}Omega }{{q}^{4}}frac{{f}_{{boldsymbol{q}},{rm{S}}{rm{i}}{rm{O}}{rm{M}}}+1}{{omega }_{{boldsymbol{q}},{rm{S}}{rm{i}}{rm{O}}{rm{M}}}}{sum _{I}frac{1}{sqrt{{M}_{I}}}{F}_{I}({boldsymbol{q}}){boldsymbol{q}}cdot {{boldsymbol{e}}}_{{boldsymbol{q}},{rm{S}}{rm{i}}{rm{O}}{rm{M}}}^{I}{e}^{i{boldsymbol{q}}cdot {{boldsymbol{tau }}}_{I}}}^{2}delta (omega {omega }_{{boldsymbol{q}},{rm{S}}{rm{i}}{rm{O}}{rm{M}}}).$$
(2)
The key aspect of this expression is the differential across the diametrically opposed aperture locations and the nonequilibrium phonon population. The specularity of interfaces creates an anisotropy in the nonequilibrium phonon population, which provides the contrast in the DPM map.
When compared to the onaxis, momentumaveraged 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 reflectioninduced nonequilibrium 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^{44}, demonstrating that the MFP of phonons can be measured at nanometre resolution.