Atomic-Scale Double-Slit Interferometry Using a Focused Electron Probe

Sample Preparation and 4D-STEM Data Acquisition

A high-purity (99.999%), undoped, commercially available silicon (Si) single crystal (Crystal Base Co., Ltd.) was mechanically crushed in a mortar. The resulting fragments were dispersed onto a molybdenum transmission electron microscopy (TEM) grid with a carbon support film. The grid was transferred immediately into the microscope to minimize surface oxidation. To remove hydrocarbon contamination, the TEM grid was annealed overnight at 300 °C under high-vacuum conditions using an in situ TEM heating holder (JEOL, Ltd.) inside the microscope. The sample was then cooled to room temperature before imaging.

Scanning transmission electron microscopy (STEM) measurements were performed using a JEM ARM300CF (JEOL, Ltd.) equipped with a cold field emission gun and a JEOL DELTA corrector. The microscope operated at an accelerating voltage of 300 kV. The probe-forming aperture semi-angle was 9.1 mrad, and the probe current was approximately 9.1 pA. This convergence semi-angle was selected so that the electron probe size matched the Si [110] dumbbell spacing. When the probe is substantially larger or smaller, the double-slit condition is not satisfied and the characteristic interference fringes disappear (Supplementary Note 9).

Four-dimensional scanning transmission electron microscopy (4D-STEM) datasets were collected with a pixelated ARINA³⁷ detector (DECTRIS Ltd.). The detector provided 192 × 192 pixels, enabling the full convergent-beam electron diffraction (CBED) pattern to be recorded at every probe position. The scan sampling interval was 12 pm for all 4D-STEM measurements.

To establish the atomic-scale double-slit geometry, the silicon crystal was accurately aligned with the [110] zone axis using Kikuchi lines and position-averaged convergent beam electron diffraction (PACBED). The electron-probe position at the centre of the Si [110] dumbbell was subsequently determined with 12-pm precision from the 4D-STEM data using a two-step numerical procedure.

The high-angle interference fringes used for correlation analysis are intrinsically resistant to small probe-positioning errors. Fringe visibility is primarily determined by the optical-path difference between atomic columns rather than by the absolute probe position. Any remaining positional uncertainty is incorporated into the finite effective source size used in the simulations (Supplementary Note 3).

First, initial atomic-column positions were obtained from reconstructed annular dark-field images using two-dimensional Gaussian peak fitting (Extended Data Fig. 8a). The dumbbell centre positions were then refined using the crystallographic symmetry of the CBED patterns³⁸, overcoming the precision limitations caused by scan noise and sample drift. Specifically, the twofold rotational symmetry (C₂) of the Si [110] projection was used to calculate a symmetry score S(r) at each probe position r. This score was defined as one minus the normalized mean squared error (MSE) between the CBED intensity and its 180°-rotated counterpart:

Here, S(r) = 1, with S ∈ [−1, 1], indicates exact twofold rotational symmetry. As shown in Extended Data Fig. 8b, the symmetry score produces sharp local maxima at high-symmetry positions. Refined dumbbell centres were identified by locating these maxima near the initial estimates.

To improve the signal-to-noise ratio while limiting the effects of sample drift at each temperature, several 4D-STEM datasets were acquired under identical experimental conditions. Typically, 10–15 scans were collected within the same field of view spanning several tens of nanometres. At each temperature, equivalent CBED patterns from the refined dumbbell centres identified across all datasets were extracted and averaged. These averaged patterns were used for the quantitative analysis of interference-fringe visibility.

In Situ TEM Heating Experiments

Temperature-dependent STEM observations were performed using the same in situ TEM heating holder. The silicon sample was sequentially heated to 300 K (room temperature), 500 K and 900 K. Sufficient time was allowed at each temperature set point to establish thermal equilibrium before data acquisition. STEM images and 4D-STEM datasets were recorded at every temperature using identical microscope and probe conditions, enabling a quantitative comparison of temperature-dependent changes in atomic-scale interference and vibrational correlations.

CBED Scattering Simulations

Convergent-beam electron diffraction (CBED) patterns were simulated using the multislice method implemented in the abTEM³⁹ code. The simulation parameters matched the experimental conditions, including an accelerating voltage of 300 kV and a probe-forming aperture semi-angle of 9.1 mrad. The simulated electron probe was positioned at the centre of the Si [110] dumbbell.

To account for finite source size and effective probe instability, adjacent diffraction patterns were mixed using Gaussian weights. The weighting function depended on the distance between the probe positions used to generate the patterns and had a full width at half maximum (FWHM) of 0.8 Å. No defocus, spherical aberration or chromatic aberration was applied because typical residual aberrations were found to have a negligible effect on the results (Supplementary Note 3).

The specimen thickness for each temperature dataset was determined by maximizing the cross-correlation coefficient between the experimental and simulated PACBED patterns⁴⁰ (Extended Data Fig. 9). Simulations confirmed that whether the left or right atomic column terminated last at the exit surface did not significantly affect the resulting diffraction patterns (Supplementary Note 10). The estimated sample thicknesses were 12.7 nm at 300 K, 10.8 nm at 500 K and 10.4 nm at 900 K.

Thermal diffuse scattering was calculated for all simulations by averaging 1,000 frozen-phonon configurations. The full correlated model included phonon vibrations in all spatial directions. Atomic displacements parallel to the electron beam have a negligible first-order effect on the projected potential and CBED intensities. Consequently, lateral atomic displacements are the dominant factors controlling interference-fringe visibility (Supplementary Note 11).

Atomic Phonon-Displacement and Correlation Models

Three approaches were used to model atomic displacements in the frozen-phonon simulations. These models ranged from independent atomic vibrations to fully correlated displacements derived from first-principles-informed phonon calculations.

Independent-displacement (Einstein) model: As a baseline for uncorrelated atomic motion, displacements were sampled from isotropic Gaussian distributions with zero interatomic correlation (ρ = 0). The vibrational amplitudes were calculated from the full phonon model described below, providing a consistent reference for comparison. The nominal root-mean-square displacements were \(\sqrt{\langle {u}^{2}\rangle }=0.08\,\mathring{{\rm{A}}}\) at 300 K, 0.10 Å at 500 K and 0.13 Å at 900 K. These values agree well with the experimental silicon Debye–Waller factors reported by Peng et al.⁴¹, which were derived from experimentally measured phonon densities of states.

Full phonon-based correlation model: Realistic vibrational correlations in the periodic silicon crystal were modelled using phonon calculations. A machine-learned Gaussian approximation potential⁴², trained on density functional theory simulations of silicon⁴³, was used to evaluate interatomic forces. Force evaluations were performed in QUIP⁴⁴ through its Python interface, quippy⁴⁵.

Second-order harmonic force constants were obtained from finite-displacement calculations in a 4 × 4 × 4 supercell using the hiPhive⁴⁶ package. Two-body interactions with a 4 Å cut-off were included. Increasing the cut-off or adding three-body terms did not appreciably change the results. The phonon dispersion relation calculated from these force constants using phonopy^47,48 (Extended Data Fig. 10) showed excellent agreement with established theoretical and experimental data. This agreement confirms that the Gaussian approximation potential reproduces the vibrational properties of silicon accurately.

Thermal displacement configurations incorporating correlated phonons were generated by superimposing harmonic normal modes with amplitudes and phases sampled according to canonical ensemble statistics. The correlated snapshots were generated in a 10 × 14 × N_z supercell constructed by repeating the Si [110] conventional unit cell. The dimension N_z along the beam-propagation direction was selected to match the experimental specimen thickness.

Nuclear quantum zero-point motion, which is particularly important at low temperatures, was included using the QM_statistics option in hiPhive⁴⁶. This option replaces classical phonon amplitudes with quantum-statistical harmonic-oscillator amplitudes⁴⁹. Equivalently, the quantum correction can be expressed using a mode-dependent effective temperature T_eff, where ħ is the reduced Planck constant, k_B is the Boltzmann constant and ω is the mode frequency:

Nearest-neighbour chain model: A simplified nearest-neighbour model was developed for parametric studies of the two adjacent silicon columns. This model captures the essential physics of coupled atomic motion and generates frozen-phonon displacement configurations at a specified temperature and correlation coefficient.

Although the full phonon calculation includes the entire crystal, the lattice vibrations along the beam direction were approximated as a harmonic chain with nearest-neighbour interactions. The cumulative influence of the surrounding bulk crystal, including atoms outside the two selected columns, was represented by an effective on-site potential. Intercolumn (x) and perpendicular (y) correlated displacements were treated as independent one-dimensional harmonic chains. The effective Hamiltonian was defined using an on-site spring constant K, representing the mean-field stiffness of the surrounding lattice, and an interatomic spring constant k, representing coupling between adjacent atoms:

Here, u_α,i is the displacement of the ith atom along the α ∈ {x, y} direction.

For a chain containing N atoms, the collective displacement vector is defined as u_α = [u_α,1, u_α,2,…, u_α,N]^T. The potential energy can then be written in matrix form as \({{V}}_{\alpha }=\frac{{\rm{1}}}{{\rm{2}}}{{{\bf{u}}}_{\alpha }}^{{\rm{T}}}{{\Phi }}_{\alpha }{{\bf{u}}}_{\alpha }\), where Φ_α is the tridiagonal force-constant matrix:

Periodic boundary conditions were assumed, corresponding to the bulk limit N → ∞, in which each atom has two nearest neighbours along the chain. The diagonal elements, K_α + 2k_α, represent the combined stiffness from the on-site potential and bonds to both neighbours. The off-diagonal elements, −k_α, describe the coupling between neighbouring atoms.

In the classical canonical ensemble, harmonic vibrations at temperature T follow the Boltzmann distribution P(u_α) ∝ exp(−V_α/k_BT). This distribution is equivalent to a multivariate normal distribution with precision matrix \({{\varSigma }_{\alpha }}^{-1}={{\Phi }}_{\alpha }/{k}_{{\rm{B}}}T\). The nearest-neighbour correlation coefficient between atoms A and B is defined as:

where ⟨…⟩ denotes the ensemble average. Analytical solution of the linear-chain model gives the following relationship between the correlation coefficient ρ_α and the stiffness ratio κ_α ≡ k_α/K_α:

The exact inverse relationship is:

This relationship connects the statistical and mechanical descriptions of atomic vibrations. In other words, extracting the correlation coefficient from the double-slit interference pattern is equivalent to determining the local bond-stiffness ratio κ_α. For the frozen-phonon simulations, atomic displacements within the two silicon columns were generated from the multivariate normal distribution defined by the precision matrix \({{\varSigma }_{\alpha }}^{-1}\) derived from ρ_α.

4D-STEM Correlation-Extraction Workflow

Correlation coefficients were extracted from the experimental 4D-STEM data shown in Fig. 4 using a systematic grid-search procedure. At each temperature, CBED patterns were simulated on a 41 × 41 grid of (ρ_x, ρ_y) values ranging from 0.0 to 0.8 in increments of 0.02.

To isolate the atomic-scale double-slit interference fringes from the thermal diffuse scattering background, the bright-field disc and low-angle Bragg diffraction regions below 27.3 mrad were masked. This approach removes low-angle intensities that are particularly sensitive to experimental imperfections, including residual aberrations and sample mistilt.

Experimental and simulated patterns were compared using the MSE of their normalized intensity distributions within the unmasked high-angle region. The optimal correlation parameters were obtained by fitting the discrete MSE landscape with a bicubic spline and locating its global minimum using the L-BFGS-B algorithm (ref. ⁵⁰) implemented in SciPy (ref. ⁵¹).

Statistical uncertainties were estimated by randomly dividing the total set of averaged CBED patterns into five independent subsets. The complete correlation-extraction procedure was repeated for each subset, and the 95% confidence intervals (CIs) for ρ_x and ρ_y were calculated from the resulting variance.

This precision was achieved using a total electron dose of approximately 3 × 10⁷ electrons. Because the uncertainty is limited primarily by shot noise and scales with the inverse square root of the dose, averaging equivalent atomic-column pairs effectively accumulates dose. Adjusting the total dose to meet the precision requirements of a particular experiment could therefore enable single-position measurements at an individual column pair.

Phonon Spectral Analysis of Vibrational Correlations

To determine which phonon modes contribute to the visibility of the atomic-scale interference fringes, the projected mean-square relative displacement (MSRD; see Supplementary Note 4 for the derivation) was decomposed into contributions from individual modes. The resulting contributions were mapped onto the phonon dispersion relation shown in Extended Data Fig. 7.

In STEM, the electron beam interacts with the atomic potential projected along the beam direction. To account for this projection geometry, the column-averaged, projected-mass-normalized eigenvector \({\mathop{{\bf{e}}}\limits^{ \sim }}_{j,{\bf{q}}\nu }\) for atom j (j = A, B) in phonon mode (q, ν) was defined as:

Here, e_j,qν is the phonon eigenvector obtained from the phonon calculations, m_j is the atomic mass and R_l is the lattice translation vector along the electron-beam direction. The summation extends over N_z unit cells corresponding to the specimen thickness. The Bloch phase factor exp(iq · R_l) determines whether displacements in successive unit cells interfere constructively or destructively in the projected signal.

For directionally resolved analysis, as shown in Fig. 4e,h, the Cartesian component α ∈ {x, y} was calculated as:

In this expression, \({\hat{{\bf{n}}}}_{x}\) is aligned with the Si–Si bond axis, while \({\hat{{\bf{n}}}}_{y}\) is perpendicular to the bond within the imaging plane.

For phonon-dispersion visualization, modes were grouped into spectral bins \({{\mathcal{D}}}_{{\bf{q}},\omega }\). Each bin contains modes at wavevector q with frequencies within a narrow window centred at ω. The contribution of each bin to the projected MSRD along direction α is:

where

is the quantum harmonic-oscillator amplitude factor. It arises from the displacement variance \(\langle {|{\bf{u}}|}^{2}\rangle =\frac{\hbar }{2m\omega }\coth \left(\frac{\hbar \omega }{2{k}_{{\rm{B}}}T}\right)\). For the Si [110] dumbbell, inversion symmetry about the bond midpoint ensures that \({|{\mathop{e}\limits^{ \sim }}_{A,{\bf{q}}\nu }^{(\alpha )}|}^{2}={|{\mathop{e}\limits^{ \sim }}_{B,{\bf{q}}\nu }^{(\alpha )}|}^{2}\).

Under this symmetry, the MSRD contribution can be factorized exactly into two terms: the thermal population factor \({B}_{{\mathcal{D}}}^{(\alpha )}\) and the spectral correlation coefficient \({c}_{{\mathcal{D}}}^{(\alpha )}\):

1. (i)

   Thermal population factor \({B}_{{\mathcal{D}}}^{(\alpha )}\): The thermally excited mean-square displacement associated with each spectral bin is:

   $${B}_{{\mathcal{D}}}^{(\alpha )}({\bf{q}},\omega )=\left(\sum _{\nu \in {\mathcal{D}}}{|{\mathop{e}\limits^{ \sim }}_{A,{\bf{q}}\nu }^{(\alpha )}|}^{2}+\sum _{\nu \in {\mathcal{D}}}{|{\mathop{e}\limits^{ \sim }}_{B,{\bf{q}}\nu }^{(\alpha )}|}^{2}\right)\,A(\omega ).$$

   This factor describes the vibrational power available in each mode, combining its zero-point amplitude, which is proportional to 1/ω, with its thermal occupation. In the high-temperature classical limit, the thermal contribution scales as \({B}_{{\mathcal{D}}}^{(\alpha )}\propto 1/{\omega }^{2}\).

2. (ii)

   Spectral correlation coefficient: The ability of a phonon mode to generate relative displacement between the two atomic columns is described by:

The spectral correlation coefficient ranges from +1, representing perfectly in-phase motion, to −1, representing perfectly out-of-phase motion. In-phase motion produces little or no relative displacement, whereas out-of-phase motion produces the maximum relative displacement. Therefore, the factor \((1-{c}_{{\mathcal{D}}}^{(\alpha )})\) describes the capacity of a phonon mode to drive relative motion between the two columns along direction α.

This factorization demonstrates the spectral filtering mechanism governing interference visibility. A phonon mode makes a substantial contribution to the MSRD, and therefore to fringe-visibility loss, only when it has both a high thermal population, represented by a large \({B}_{{\mathcal{D}}}\), and a strong capacity to produce relative atomic motion, represented by a large \((1-{c}_{{\mathcal{D}}})\).

The mode-resolved results are presented in Fig. 4e,h and Extended Data Fig. 7. In Extended Data Fig. 7a,c, the thermal population factor is represented by the line width of the phonon-dispersion curves. The spectral correlation coefficient is shown by line colour: red indicates in-phase motion (\({c}_{{\mathcal{D}}}\simeq 1\)), grey indicates negligible correlation (\({c}_{{\mathcal{D}}}\simeq 0\)) and blue indicates out-of-phase motion (\({c}_{{\mathcal{D}}}\simeq -1\)). The resulting MSRD contribution, W_MSRD, is represented by colour intensity in Fig. 4e,h and Extended Data Fig. 7b,d, highlighting the phonon modes that dominate temperature-dependent interference-fringe loss.