Atomic-scale double-slit interferometry with a focused electron probe

Since Thomas Young’s original demonstration with light in the early 1800s¹, double-slit experiments have evolved far beyond a proof of wave–particle duality into a precision framework for extracting quantitative information from interference fringe visibility and phase. Over the past century, this model has extended from photons to matter waves, progressing to increasingly massive particles—electrons², neutrons³, atoms⁴ and even complex molecules^5,6—underpinning a diverse array of transformative measurement platforms. Today, electron holography^7,8, atom interferometry^9,10 and gravitational-wave observatories¹¹ exemplify how fringe visibility and phase serve as precision readouts for electromagnetic fields, inertial forces and spacetime curvature, respectively. Pushing interferometry towards atomic scales, recent experiments have used crystal lattices¹² and diatomic molecules^13,14 as effective slits. The unifying power of interferometry lies in its ability to encode physical quantities into interference patterns; shrinking the interferometer to atomic dimensions thus offers the prospect of investigating matter at the single-bond level. This atomic confinement is expected to bring intrinsic spatial localization, potentially granting access to symmetry-broken regions such as defects and interfaces.

Crystalline solids present an ideal architecture for this goal: atoms locked at lattice sites could act as double slits with well-defined separations, offering direct access to local structure and dynamics. However, extracting such local information from real crystals requires circumventing the spatial averaging inherent to the periodic lattice by strictly confining the interaction to a focused atomic-scale volume. This challenge is addressed by recent advances in STEM¹⁵. Aberration-corrected electron probes allow for the targeting of individual atomic pairs, whereas segmented and pixelated electron detectors capture detailed convergent-beam electron diffraction (CBED) patterns at each probe position (for example, differential phase contrast STEM^16,17, optimum bright-field STEM¹⁸ and four-dimensional STEM (4D-STEM)¹⁹) with sufficient sensitivity and speed. Here we demonstrate that crystalline materials can function as true atomic-scale double-slit interferometers. We show that the interference fringes encode hidden signatures of correlated thermal vibrations between the selected pair of atomic columns—synchronous atomic motions that preserve coherence of the electron wave—unlocking a new window into local lattice dynamics at the atomic scale.

We realize atomic-scale interferometry by a focused electron probe transmitting through adjacent silicon atomic columns using STEM equipped with a pixelated detector. From a 4D-STEM dataset acquired with a 1.1 Å full width at half maximum (FWHM) probe, we extracted CBED patterns corresponding to probe positions at the midpoint between Si [110] atomic column pairs separated by 1.36 Å—downscaling Young’s experiment by seven orders of magnitude (Fig. 1a,b). As the electrons propagate through the crystal, the attractive potential of the atomic nuclei acts as a waveguide, confining the beam intensity onto the two columns (Fig. 1c,d). This channelling effect^20,21,22 effectively creates two separated point sources on the exit surface at the atomic scale (Fig. 1e), which generate an interference pattern on the detector. Although this localization is rooted in the same physics as the atomic-focuser concept²², and these fringes share the dynamical-diffraction origin of conventional CBED²³, our sub-unit-cell, position-selected geometry turns the crystal into the interferometric source itself, making the fringe visibility sensitive to the relative motion of a single atomic-column pair.

Figure 2a shows the experimental CBED pattern acquired under this double-slit geometry. Distinct interference fringes appear perpendicular to the intercolumnar direction, with a periodicity of 0.736 Å⁻¹ corresponding to the projected Si-column distance of 1.36 Å. Averaging 356 crystallographically equivalent patterns across the same thickness region (Fig. 2b; Methods) extends fringe visibility up to the third bright maximum. The geometric agreement between fringe spacing and slit separation demonstrates that the single pair of atomic columns acts as a pair of coherent electron sources, effectively functioning as an atomic-scale double slit. By contrast, these periodic fringes vanish when the probe is placed on a single Si column (Extended Data Fig. 1), confirming that the intercolumn position constitutes a specialized interferometric condition unique to the closely spaced atomic pair.

To confirm that the observed fringes genuinely arise from atomic-scale double-slit interference, we first performed a benchmark simulation using a simplified model. We defined two circular apertures with diameters determined from the point-source intensity distribution (Fig. 1e), separated by 1.36 Å to match the Si [110] dumbbell geometry. We then calculated the far-field interference pattern. The resulting pattern (Fig. 2c) exhibits interference fringes with spacing identical to the experimental observations, reproducing both the fringe periodicity from the 1.36 Å slit separation and the intensity envelope arising from finite aperture size. Furthermore, multislice simulations²⁴ using a static Si structural model, in which atoms are frozen at their equilibrium positions, also produced identical interference fringes (Fig. 2d). These correspondences confirm the atomic-scale double-slit origin of the interference.

To calculate realistic interference patterns, we performed multislice frozen-phonon simulations²⁵ (Methods). The simplest assumption—independent isotropic vibrations (Einstein model)—produces a pattern in which all fringes higher than the first order vanish (Fig. 2e), in clear disagreement with experiment. Including correlated displacements derived from the phonon modes of silicon recovers the experimental fringe pattern, including the high-order maxima (Fig. 2f). Quantitative line profiles extracted from inner and outer detector windows (Fig. 2g,h) further confirm this picture: the correlated phonon model accurately reproduces the experimental fringe profiles in both windows, whereas the Einstein model deviates substantially. This discrepancy indicates that correlated thermal vibrations are essential for preserving atomic-scale interference.

Quantum excitation of phonons (QEP²⁶) simulations, which separately calculate elastic and thermal diffuse scattering (inelastic) components, further explain this behaviour (Extended Data Fig. 2). At finite temperatures, atomic displacements owing to phonon excitations generate thermal diffuse scattering, which reduces the elastic scattering intensity, particularly at higher angles, compared with the static case²¹. Notably, when correlated displacements are incorporated, coherent interference fringes persist within the thermal diffuse scattering itself, indicating that the experimentally observed visibility is largely carried by the inelastic scattering, rather than only by the residual elastic signal.

Crucially, this interference is robust against thermal agitation. We extended our measurements to increased temperatures up to 900 K (Extended Data Fig. 3a–c). Even under these high-temperature conditions, in which thermal displacements are much larger, interference fringes persist in the experimental data: whereas the third-order fringe fades, the second-order maximum remains visible at 900 K (Extended Data Fig. 3c). This persistence is structurally impossible to reproduce with the Einstein model, which yields featureless diffuse scattering (Extended Data Fig. 3g–i) but is well reproduced by the correlated model (Extended Data Fig. 3d–f).

To characterize these correlations quantitatively, we define x along the projected intercolumnar axis and y perpendicular to it, in the [110] projection in which the two probed columns form a zigzag chain of nearest-neighbour Si–Si bonds (Fig. 3a,b). Sampling the bivariate distribution of nearest-neighbour displacements from the phonon-based model (Fig. 3d) yields correlation coefficients (see Methods for formal definition) of ρ_x = 0.39 and ρ_y = 0.24, with negligible cross-correlation (ρ_xy ≈ 0.0014). Furthermore, systematic sensitivity tests examining how strongly the atomic-scale interference pattern is affected by spatially localized correlated vibrations (Supplementary Note 1) confirm that only intercolumn correlation between the two investigated Si columns preserves interference. We therefore capture this correlation characterized with just two independent parameters ρ_x and ρ_y, which map directly onto the bond-stiffness ratio κ_α = k_α/K_α between effective interatomic and on-site (lattice) stiffnesses, with the explicit relation κ_α = ρ_α/(1 − ρ_α)². The diffraction patterns calculated with this bivariate model accurately reproduce both the full-phonon result (ρ_x = 0.39 and ρ_y = 0.24; Fig. 3f) and the Einstein model (ρ_x = ρ_y = 0.0; Fig. 3e). This result confirms that the simple two-parameter model fully captures the coherence-preserving correlations. Note that, when combined with the local mean squared displacement, this model provides the absolute local stiffnesses k_α within the harmonic regime, revealing a possible local alternative to macroscopic elastic measurements (Supplementary Note 2).

Systematically varying ρ_x and ρ_y (Fig. 3g,h) reveals that the two parameters modulate the interference pattern along the perpendicular directions. Increasing ρ_x extends the visibility of high-order interference fringes along the intercolumnar axis, whereas increasing ρ_y extends the fringes perpendicular to that axis. These two effects appear in the line profiles of two different angular windows (compare Fig. 2g,h). In the inner window (Fig. 3i), the profile is dominated by changes in ρ_x, whereas in the outer window (Fig. 3j) it is governed by ρ_y and is essentially insensitive to ρ_x. This orthogonality enables independent extraction of both coefficients from experimental data by matching the intensity distributions.

The underlying mechanisms are illustrated in Extended Data Fig. 4. Positive ρ_x corresponds to in-phase motion along the bond direction (Extended Data Fig. 4a): the instantaneous slit separation of 1.36 Å is preserved, so a consistent fringe periodicity survives ensemble averaging. Weak or negative ρ_x implies out-of-phase motion that causes the effective slit separation to fluctuate, washing out high-order fringes. Positive ρ_y maintains the alignment of the atomic pair relative to the beam axis (Extended Data Fig. 4b); a lack of y correlation introduces relative displacements that tilt the double-slit axis configuration by configuration and averaging over these rotated patterns produces the angular blurring observed at high scattering angles.

With ρ_x and ρ_y encoded in spatially orthogonal regions of the diffraction pattern, the two coefficients can be extracted by joint mean squared error (MSE) minimization (Methods). Figure 4a presents the resulting MSE landscape for the 300 K dataset. A single, well-defined minimum appears near ρ_x ≈ 0.4 and ρ_y ≈ 0.2, within an elliptical low-MSE basin whose principal axes are aligned with the parameter grid. This shape provides an experimental confirmation that ρ_x and ρ_y are independently extractable from the data, validating our geometric interpretation of the previous section. Continuous optimization with the spline-interpolated landscape yields the correlation coefficients at 300 K to be ρ_x = 0.362 ± 0.005 and ρ_y = 0.184 ± 0.005 (95% confidence interval (CI)), consistent with the theoretical predictions from lattice dynamics (ρ_x = 0.39, ρ_y = 0.24). Because the extraction is dominated by the high-angle thermal diffuse scattering background, it is robust against residual aberrations, mistilt and surface amorphization (Supplementary Note 3).

We also applied this quantitative extraction workflow to the datasets at 500 and 900 K. Figure 4b,c shows the MSE landscapes for each temperature. In all cases, the landscapes exhibit well-defined basins. The coefficients were determined to be (0.402 ± 0.012, 0.263 ± 0.004) at 500 K and (0.418 ± 0.009, 0.288 ± 0.010) at 900 K. To verify the accuracy of these estimates, we performed CBED simulations using the experimentally extracted (ρ_x, ρ_y) values for each temperature (Extended Data Fig. 5a–c). These simulations show excellent visual agreement with the experimental data (Extended Data Fig. 3). For a rigorous quantitative check, we calculated ratio maps between the experimental and simulated patterns (Extended Data Fig. 5d–f). These maps show a uniform distribution near unity (white) without any discernible residual interference patterns or systematic deviation, in strong contrast to the pronounced periodic residuals observed when assuming uncorrelated atomic motion (Extended Data Fig. 5g–i). This confirms that our nearest-neighbour correlation model, parameterized by the extracted coefficients, successfully captures the quantitative details of the atomic-scale interference across the entire temperature range.

Figure 4d summarizes the temperature dependence of the extracted parameters (circles). Comparison with theoretical predictions derived from the phonon-based correlation model (solid lines connecting calculated points at the discrete temperatures, marked by ×; Extended Data Fig. 6) reveals two key features. First, the experimental results confirm the predicted anisotropy in which the bond-parallel correlation ρ_x consistently exceeds the perpendicular correlation ρ_y. Second, despite the mean squared displacements increasing substantially with temperature, the correlation coefficients show only a moderate temperature dependence, exhibiting a relatively flat trend comparable with the theoretical prediction. This stability validates that these coefficients are intrinsic signatures of the interatomic potential.

The origin of this stability lies in the way different phonon modes contribute to the projected mean squared relative displacement (MSRD), defined as \(\langle {({u}_{A,\alpha }-{u}_{B,\alpha })}^{2}\rangle \): the quantity that controls fringe attenuation, rather than the overall vibrational amplitude (Supplementary Note 4). For a given phonon mode, its contribution to the MSRD is determined by two factors: (1) the thermal population (that is, the mode amplitude set by temperature and frequency) and (2) the extent to which the mode generates relative motion between the two columns, encoded in the mode eigenvector. Using these two factors, we quantify the mode-resolved contributions to the MSRD and thereby identify the specific phonon modes that predominantly control the interference visibility.

Figure 4e,h shows the mode-resolved contributions to the MSRD (x and y components, respectively; see also Methods and Extended Data Fig. 7). These spectra demonstrate that the lattice does not degrade coherence uniformly across all vibrations: only a narrow subset of phonon modes contributes appreciably to the MSRD and therefore to the fringe-visibility loss. Long-wavelength acoustic modes near the Γ point are strongly thermally excited but they drive the two columns predominantly in phase, producing little relative displacement and hence a small MSRD contribution. Conversely, optical modes tend to generate relative motion between neighbouring atoms but their high frequencies make their thermal population small at the experimental temperatures, again limiting their contribution to the MSRD. As a result, the dominant MSRD contribution is concentrated in short-wavelength acoustic modes, specifically those with wavevectors near the Brillouin zone boundary. Representative examples of these strongly contributing modes are shown in Fig. 4f,g,i,j, in which the relative-motion character is visually evident.

In the complementary real-space picture, this zone-boundary selectivity provides a direct readout of local bond mechanics. In our geometry, relative displacements measured along the x-axis contain a bond-parallel (stretching/compressing) component, whereas the y-direction predominantly examines transverse (shear) motion of the same atom pair. Under the harmonic approximation, the atomic displacement covariances are determined by the (pseudoinverse of the) force constant matrix through the equipartition relation²⁷, establishing a direct correspondence between measured vibrational correlations and bond stiffness. The observed anisotropy (ρ_x > ρ_y) therefore indicates stronger constraints (that is, higher stiffness) for bond-parallel distortions than for transverse distortions. Consistent with this picture, the mode-resolved maps (Fig. 4e,h) show that high-frequency longitudinal acoustic (LA) modes involving bond stretching (Fig. 4f) are stiff and thus thermally suppressed, resulting in minimal visibility loss along the x-axis, which is partly aligned with the bond axis. Instead, the visibility loss is dominated by softer transverse acoustic (TA) modes—such as the L-point phonons (Fig. 4g,j) and the X-point phonon (Fig. 4i)—which induce shear deformations at lower energy costs. Our interferometer thereby functions as a phonon-mode-selective probe, particularly sensitive to these low-energy acoustic modes that affect local thermal transport. Notably, whereas the phonon-mode picture relies on crystal periodicity, the real-space bond-stiffness interpretation is local and therefore extends naturally to aperiodic environments such as defects, interfaces and disordered regions. From another point of view, this mode selectivity has a natural counterpart in the ‘which-path’ reading of phonon-induced electron decoherence²⁸ (Supplementary Note 5).

By recasting STEM as an atomic-scale interferometer, this work enables extraction of the directional correlation coefficients (ρ_x, ρ_y) between specific atomic pairs—a local, anisotropic, direction-resolved observable that can be interpreted as an effective bond stiffness for that pair. Such a quantity is not accessible by established techniques based on the pair-correlation function^29,30. Unlike inelastic neutron or X-ray scattering methods, which measure phonon dispersion in momentum space by averaging over macroscopic volumes, our approach visualizes the corresponding correlations in real space. It also complements a growing set of atomic-scale electron methods for accessing lattice dynamics: multislice electron ptychography provides sub-atomic-resolution imaging that is increasingly sensitive to subtle signatures of atomic vibrations³¹; quantitative electron diffraction with an atomic-scale probe yields vibrational amplitudes and anisotropies³²; and monochromated electron energy-loss spectroscopy resolves local phonon states at the atomic scale³³. Whereas these methods examine single atomic columns or energy-resolved spectra, our interferometry directly provides complementary information on a selected atomic pair—its vibrational correlation and the associated bond stiffness—and, because fringe visibility is governed by the relative displacement between the atoms, is intrinsically weighted towards the large-amplitude, low-energy acoustic modes (Supplementary Note 6).

Because the present approach relies on channelling of the electron probe along the atomic columns, it can be applied to a relatively broad range of specimen thicknesses—from thin (7.4 nm) to thick (30.3 nm) samples—when appropriate conditions are selected (Supplementary Note 7). In principle, beam-splitting interferometric STEM (for example, STEM holography³⁴) could further generalize this concept from a fixed nearest-neighbour geometry to user-defined atomic pairs by forming two mutually coherent probes with a controllable separation. This ability to target correlations at the level of a chosen atomic pair is particularly powerful where the local structure deviates from bulk periodicity, including chemically inequivalent column pairs as in GaAs (Supplementary Note 8). Combined with established atomic-resolution STEM structural analysis^35,36, this method opens a route to extracting atomic-scale vibrational correlations at grain boundaries, interfaces and other defects, for which the assumptions of bulk periodicity no longer apply. Establishing STEM as an atomic-scale interferometer thus transforms the diffraction pattern itself into a direct probe of correlated lattice dynamics, enabling bond-resolved studies at sites at which local atomic fluctuations govern materials functionality.

Sample preparation and 4D-STEM

A high-purity (99.999%), undoped, commercially available Si single crystal (Crystal Base Co., Ltd.) was mechanically crushed in a mortar and the resulting fragments dispersed onto a molybdenum TEM grid with a carbon support film and the grid was immediately transferred into the microscope to avoid 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 and subsequently cooled to room temperature before observation. STEM observations were performed using a JEM ARM300CF (JEOL, Ltd.) equipped with a cold field emission gun and a JEOL DELTA corrector, operated at an accelerating voltage of 300 kV. The probe-forming aperture semi-angle was set to 9.1 mrad and the probe current was estimated to be approximately 9.1 pA. This convergence semi-angle was chosen so that the probe size matches the Si [110] dumbbell spacing; for substantially larger or smaller probes, the double-slit condition is not satisfied and the characteristic interference fringes do not appear (Supplementary Note 9). 4D-STEM datasets were acquired using a pixelated detector (ARINA³⁷, DECTRIS Ltd.) with 192 × 192 pixels to record the full CBED pattern at each probe position. The scan sampling interval was set to 12 pm for all 4D-STEM measurements.

To realize the atomic-scale double-slit geometry, the crystal was precisely oriented to 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 determined a posteriori with 12-pm precision from the acquired 4D-STEM dataset using a two-step numerical procedure. We note that the high-angle interference fringes used for correlation extraction are inherently robust against small probe positioning offsets, as the fringe visibility is governed by the intercolumn optical-path difference rather than the absolute probe position, and any residual positional uncertainty is already subsumed within the finite effective source size accounted for in the simulations (Supplementary Note 3). First, initial atomic column positions were estimated from the reconstructed annular dark-field images using 2D Gaussian peak fitting (Extended Data Fig. 8a). Then, to overcome precision limits imposed by scan noise and sample drift, we refined the dumbbell centre positions by exploiting the crystallographic symmetry of the CBED patterns³⁸. Using the twofold rotational symmetry (C₂) of the Si [110] projection, we calculated a symmetry score S(r) at each probe position r, defined as 1 minus the normalized MSE between the CBED intensity and its 180°-rotated counterpart:

in which S(r) = 1 (S ∈ [−1, 1]) corresponds to exact twofold rotational symmetry. As shown in Extended Data Fig. 8b, this score exhibits sharp local maxima at high-symmetry points. We identified the refined dumbbell centres by locating these maxima in the vicinity of the coarse estimates. To achieve high signal-to-noise ratios while minimizing the impact of sample drift at each temperature, we acquired several 4D-STEM datasets (typically 10–15 scans within the same tens of nanometres field of view) under identical experimental conditions. For each temperature, equivalent CBED patterns corresponding to the refined dumbbell centres identified across these datasets were extracted and averaged to produce the high signal-to-noise ratio experimental patterns used for quantitative visibility analysis.

In situ heating experiments

Temperature-dependent observations were conducted using the same in situ TEM heating holder. The sample was sequentially heated to 300 K (room temperature), 500 K and 900 K, with sufficient time allowed at each temperature set point to ensure thermal equilibrium before image acquisition. STEM images and 4D-STEM datasets were recorded at each temperature under identical optical conditions to enable quantitative comparison of temperature-dependent changes.

Scattering simulations

CBED patterns were simulated using the multislice method implemented in the abTEM³⁹ code. The microscope parameters were set to match the experimental conditions: an accelerating voltage of 300 kV and a probe-forming aperture semi-angle of 9.1 mrad. The electron probe was positioned at the centre of the Si [110] dumbbell. To account for the finite source size and effective probe instability, we mix adjacent diffraction patterns using Gaussian weights with a FWHM of 0.8 Å as a function of the distance between the probe positions at which the patterns were simulated. No aberrations (defocus, spherical or chromatic) were applied, as the effect of typical residual aberrations was found to be negligible (Supplementary Note 3). The sample thickness for each temperature dataset was determined by maximizing the cross-correlation coefficient between experimental and simulated PACBED patterns⁴⁰ (Extended Data Fig. 9). We confirmed by simulation that whether the left or the right atomic column terminates last at the exit surface makes no notable difference to the resulting patterns (Supplementary Note 10). The estimated thicknesses were 12.7 nm at 300 K, 10.8 nm at 500 K and 10.4 nm at 900 K. For all simulations, thermal diffuse scattering was calculated by averaging over 1,000 frozen-phonon configurations. In the full correlated model, phonon vibrations in all directions are included. Note that atomic displacements parallel to the beam have a negligible first-order effect on the projected potential and hence on the CBED intensities, so the lateral displacements are the dominant factors governing fringe visibility (Supplementary Note 11).

Phonon-displacement models

We used three distinct approaches to model atomic displacements for the frozen-phonon calculations, ranging from independent vibrations to full ab initio-derived correlations.

Independent displacements (Einstein) model: as a baseline for uncorrelated motion, atomic displacements were sampled from isotropic Gaussian distributions with zero interatomic correlation (ρ = 0). The vibrational amplitudes were determined from the full phonon calculations described below to provide a consistent reference for comparison with the correlated model. We used nominal root mean square displacements of \(\sqrt{\langle {u}^{2}\rangle }=0.08\,\mathring{{\rm{A}}}\) at 300 K, 0.10 Å at 500 K and 0.13 Å at 900 K, which are in good agreement with the experimental Si Debye–Waller values tabulated by Peng et al.⁴¹, derived from experimentally determined phonon densities of states.

Full phonon-based correlation model: to capture realistic vibrational correlations of the periodic crystal, we performed phonon calculations. We used a machine-learned Gaussian approximation potential⁴² trained on density functional theory simulations for silicon⁴³. Interatomic force evaluations were performed in QUIP⁴⁴ through its Python interface, quippy⁴⁵. Second-order harmonic force constants were then extracted by means of finite-displacement calculations in a 4 × 4 × 4 supercell using the hiPhive⁴⁶ package, including two-body terms with a cut-off of 4 Å (a higher cut-off or the inclusion of three-body terms did not appreciably alter the results). The calculated phonon dispersion relation derived from these force constants using phonopy^47,48 (Extended Data Fig. 10) shows excellent agreement with established theoretical and experimental data, confirming that the force constants extracted from the Gaussian approximation potential appropriately reproduce the vibrational properties of silicon.

Thermal atomic displacements with correlated phonons were generated by superimposing harmonic normal modes with amplitudes and phases sampled to satisfy canonical ensemble statistics. Correlated displacement snapshots were generated in a 10 × 14 × N_z supercell constructed by repeating the Si [110] conventional unit cell, in which N_z is the supercell dimension along the beam-propagation direction and was set to correspond to the specimen thickness used in the experiment. To account for nuclear quantum zero-point motion, which is non-negligible especially at lower temperatures, we used the QM_statistics option in hiPhive⁴⁶, in which the classical phonon amplitudes are replaced by quantum-statistical harmonic-oscillator amplitudes⁴⁹. Equivalently, this can be written in terms of a mode-dependent effective temperature T_eff (in which ħ is the reduced Planck constant, k_B the Boltzmann constant and ω the mode frequency):

Nearest-neighbour chain model: for parametric studies, we constructed a simplified model that captures the essential physics of nearest-neighbour coupling within the two adjacent Si columns. This model enables us to generate atomic displacements for frozen-phonon simulation at a given temperature using specific correlation coefficients. Although the full phonon calculation involves the entire crystal, we modelled the lattice vibrations along the beam direction effectively as a harmonic chain with nearest-neighbour interactions. In this picture, the cumulative interaction with the surrounding bulk crystal—atoms other than those in the two columns—is renormalized into an effective on-site potential. We separately treated the intercolumn (x) and perpendicular (y) correlated displacements as independent 1D harmonic chains. The effective Hamiltonian is defined by an on-site spring constant K (representing the mean-field stiffness provided by the surrounding lattice) and an interatomic spring constant k (coupling adjacent atoms within the chain):

in which u_α,i represents the displacement of the ith atom along the α ∈ {x, y} axis.

In the language of statistical mechanics, we consider the collective displacement vector u_α = [u_α,1, u_α,2,…, u_α,N]^T for a chain of N atoms. The potential energy term can be written in matrix form as \({{V}}_{\alpha }=\frac{{\rm{1}}}{{\rm{2}}}{{{\bf{u}}}_{\alpha }}^{{\rm{T}}}{{\Phi }}_{\alpha }{{\bf{u}}}_{\alpha }\). Here Φ_α is the tridiagonal force-constant matrix along the α axis, which explicitly takes the form:

Here we assume periodic boundary conditions (or focus on the bulk limit N → ∞), so that each atom has two nearest neighbours along the chain. The diagonal elements (K_α + 2k_α) represent the total stiffness felt by an atom owing to the on-site potential and bonds to both neighbours, whereas the off-diagonal elements (−k_α) represent the coupling. In the classical canonical ensemble, harmonic vibrations at temperature T follow the Boltzmann distribution P(u_α) ∝ exp(−V_α/k_BT), which is mathematically equivalent to a multivariate normal distribution with the precision matrix \({{\varSigma }_{\alpha }}^{-1}={{\Phi }}_{\alpha }/{k}_{{\rm{B}}}T\). The correlation coefficient ρ_α between nearest-neighbour atoms A and B is formally defined as the normalized covariance of their displacements:

in which ⟨…⟩ denotes the ensemble average. Analytical solution of this linear-chain system yields a direct relationship between the correlation coefficient ρ_α and the stiffness ratio κ_α ≡ k_α/K_α:

with the exact inverse relation

This relationship allows us to bridge the statistical and physical pictures: extracting the correlation coefficient from the interference pattern is equivalent to determining the local stiffness ratio κ_α of the atomic bonds. For the frozen-phonon simulations, we generated the displacements of atoms within two Si columns from the multivariate normal distribution defined by the precision matrix \({{\varSigma }_{\alpha }}^{-1}\) derived from ρ_α.

Correlation extraction workflow

To extract the correlation coefficients from experimental data (Fig. 4), we performed a systematic grid search. For each temperature, we simulated CBED patterns on a 41 × 41 grid of (ρ_x, ρ_y) values ranging from 0.0 to 0.8 in steps of 0.02. To specifically isolate the atomic-scale double-slit interference fringes emerging on the thermal diffuse scattering background, we masked out the bright-field disc and the low-angle Bragg diffraction regions (<27.3 mrad). This masking strategy effectively excludes low-angle intensities that are highly sensitive to experimental imperfections (such as residual aberrations and sample mistilt). The agreement between experiment and simulation was evaluated using the MSE of the normalized intensity distributions within this unmasked high-angle region. The optimal parameters were determined by fitting the discrete MSE landscape with a bicubic spline function and finding the global minimum of the function using the L-BFGS-B (ref. ⁵⁰) optimization algorithm implemented in the SciPy (ref. ⁵¹) package. To estimate the statistical uncertainties of the extracted coefficients, we randomly partitioned the total averaged CBED patterns into five independent subsets. The entire extraction procedure was repeated for each subset and the 95% CI for both ρ_x and ρ_y was calculated from the resulting variance. Note that this precision was achieved from a total dose of about 3 × 10⁷ electrons. Because the uncertainty is limited by shot noise and scales with the inverse square root of the dose, and averaging over equivalent pairs simply accumulates dose, tuning the total dose to the precision required for a given problem could make single-position acquisition at an individual column pair feasible.

Spectral analysis of vibrational correlations

To identify which phonon modes contribute to the interference visibility, we decomposed the projected MSRD (see Supplementary Note 4 for derivation) into individual mode contributions and visualized them on the phonon dispersion relation, presented in Extended Data Fig. 7.

In STEM, the electron beam interacts with the atomic potential integrated along the column. To account for this projective geometry, we define the column-averaged projected-mass-normalized eigenvector \({\mathop{{\bf{e}}}\limits^{ \sim }}_{j,{\bf{q}}\nu }\) for atom j (j = A, B for the two columns) in mode (q, ν) as

in which e_j,qν is the phonon eigenvector from phonon calculations, m_j is the atomic mass and R_l denotes lattice translation vectors along the beam direction. The summation runs 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 (Fig. 4e,h), we extract the Cartesian component α ∈ {x, y}:

in which \({\hat{{\bf{n}}}}_{x}\) is aligned with the Si–Si bond axis and \({\hat{{\bf{n}}}}_{y}\) is perpendicular to it within the imaging plane.

For visualization of the phonon dispersion, we group modes into spectral bins \({{\mathcal{D}}}_{{\bf{q}},\omega }\), defined as the set of modes at wavevector q with frequency within a small window centred at ω. The contribution of each bin to the projected MSRD along direction α is given by:

in which

is the quantum harmonic oscillator amplitude factor, arising from the 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 structure, inversion symmetry about the bond midpoint ensures \({|{\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 factorizes exactly into a thermal population term \({B}_{{\mathcal{D}}}^{(\alpha )}\) (i) and a spectral correlation coefficient term \({c}_{{\mathcal{D}}}^{(\alpha )}\) (ii):

1. (i)

   Thermal population \({B}_{{\mathcal{D}}}^{(\alpha )}\): the thermally excited mean squared displacement is quantified by

   $${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 quantity represents the vibrational power available from each mode, combining the zero-point amplitude (∝ 1/ω) and thermal occupation. In the high-temperature (classical) limit, the thermal part scales as \({B}_{{\mathcal{D}}}^{(\alpha )}\propto 1/{\omega }^{2}\).

2. (ii)

   Spectral correlation coefficient: the extent to which a mode generates relative displacement between the two columns is quantified by the correlation coefficient:

This coefficient ranges from +1 (perfectly in-phase motion, generating no relative displacement) to −1 (perfectly out-of-phase motion, generating maximum relative displacement). The factor \((1-{c}_{{\mathcal{D}}}^{(\alpha )})\) thus represents the capacity of the mode to induce relative motion between the two columns along direction α.

This factorization confirms the spectral filtering mechanism: marked contribution to the MSRD (and hence to visibility loss) requires both substantial thermal population (\({B}_{{\mathcal{D}}}\) large) and substantial relative-motion capacity (\((1-{c}_{{\mathcal{D}}})\) large).

The mode-resolved analysis is shown in Fig. 4e,h and Extended Data Fig. 7. In Extended Data Fig. 7a,c, the thermal population factor is encoded in the line width of the dispersion curves and the spectral correlation coefficient is shown by the line colour (red: in-phase, \({c}_{{\mathcal{D}}}\simeq 1\); grey: no correlation, \({c}_{{\mathcal{D}}}\simeq 0\); blue: out-of-phase, \({c}_{{\mathcal{D}}}\simeq -1\)). The resulting MSRD contribution W_MSRD is visualized by the colour intensity in Fig. 4e,h and Extended Data Fig. 7b,d, highlighting the modes that dominate visibility loss.

The simulation code used in this study has been deposited in Zenodo (https://doi.org/10.5281/zenodo.18189080)⁵².

1. Young, T. II. The Bakerian Lecture. On the theory of light and colours. Philos. Trans. R. Soc. Lond. 92, 12–48 (1802).

   ADS  Google Scholar 

2. Tonomura, A., Endo, J., Matsuda, T., Kawasaki, T. & Ezawa, H. Demonstration of single-electron buildup of an interference pattern. Am. J. Phys. 57, 117–120 (1989).

   Article  ADS  Google Scholar 

3. Zeilinger, A., Gähler, R., Shull, C. G., Treimer, W. & Mampe, W. Single- and double-slit diffraction of neutrons. Rev. Mod. Phys. 60, 1067–1073 (1988).

   Article  ADS  CAS  Google Scholar 

4. Keith, D. W., Ekstrom, C. R., Turchette, Q. A. & Pritchard, D. E. An interferometer for atoms. Phys. Rev. Lett. 66, 2693–2696 (1991).

   Article  ADS  CAS  PubMed  Google Scholar 

5. Arndt, M. et al. Wave–particle duality of C₆₀ molecules. Nature 401, 680–682 (1999).

   Article  ADS  CAS  PubMed  Google Scholar 

6. Fein, Y. Y. et al. Quantum superposition of molecules beyond 25 kDa. Nat. Phys. 15, 1242–1245 (2019).

   Article  CAS  Google Scholar 

7. Tonomura, A. (ed) in Electron Holography 29–49 (Springer, 1999).

8. Lichte, H. & Lehmann, M. Electron holography—basics and applications. Rep. Prog. Phys. 71, 016102 (2007).

   Article  ADS  Google Scholar 

9. Peters, A., Chung, K. Y. & Chu, S. Measurement of gravitational acceleration by dropping atoms. Nature 400, 849–852 (1999).

   Article  ADS  CAS  Google Scholar 

10. Cronin, A. D., Schmiedmayer, J. & Pritchard, D. E. Optics and interferometry with atoms and molecules. Rev. Mod. Phys. 81, 1051–1129 (2009).

    Article  ADS  CAS  Google Scholar 

11. LIGO Scientific Collaboration and Virgo Collaboration. Observation of gravitational waves from a binary black hole merger. Phys. Rev. Lett. 116, 061102 (2016).

    Article  ADS  MathSciNet  Google Scholar 

12. Kanitz, C. et al. Diffraction of helium and hydrogen atoms through single-layer graphene. Science 389, 724–726 (2025).

    Article  ADS  CAS  PubMed  Google Scholar 

13. Zhou, H., Perreault, W. E., Mukherjee, N. & Zare, R. N. Quantum mechanical double slit for molecular scattering. Science 374, 960–964 (2021).

    Article  ADS  CAS  PubMed  Google Scholar 

14. Akoury, D. et al. The simplest double slit: interference and entanglement in double photoionization of H₂. Science 318, 949–952 (2007).

    Article  ADS  CAS  PubMed  Google Scholar 

15. Pennycook, S. J. & Nellist, P. D. Scanning Transmission Electron Microscopy: Imaging and Analysis (Springer, 2011).

16. Shibata, N. et al. Differential phase-contrast microscopy at atomic resolution. Nat. Phys. 8, 611–615 (2012).

    Article  CAS  Google Scholar 

17. Kohno, Y., Seki, T., Findlay, S. D., Ikuhara, Y. & Shibata, N. Real-space visualization of intrinsic magnetic fields of an antiferromagnet. Nature 602, 234–239 (2022).

    Article  ADS  CAS  PubMed  Google Scholar 

18. Ooe, K. et al. Direct imaging of local atomic structures in zeolite using optimum bright-field scanning transmission electron microscopy. Sci. Adv. 9, eadf6865 (2023).

    Article  PubMed  PubMed Central  Google Scholar 

19. Ophus, C. Four-dimensional scanning transmission electron microscopy (4D-STEM): from scanning nanodiffraction to ptychography and beyond. Microsc. Microanal. 25, 563–582 (2019).

    Article  ADS  CAS  PubMed  Google Scholar 

20. Uggerhøj, E. & Frandsen, F. Fast-electron channeling investigated by means of Rutherford scattering. Phys. Rev. B 2, 582–590 (1970).

    Article  ADS  Google Scholar 

21. Pennycook, S. J. & Jesson, D. E. High-resolution Z-contrast imaging of crystals. Ultramicroscopy 37, 14–38 (1991).

    Article  ADS  Google Scholar 

22. Cowley, J. M. Electron holography with atomic focusers. Phys. Rev. Lett. 84, 3618–3621 (2000).

    Article  ADS  CAS  PubMed  Google Scholar 

23. Cowley, J. M. Coherent interference in convergent-beam electron diffraction and shadow imaging. Ultramicroscopy 4, 435–449 (1979).

    Article  CAS  Google Scholar 

24. Cowley, J. M. & Moodie, A. F. The scattering of electrons by atoms and crystals. I. A new theoretical approach. Acta Crystallogr. 10, 609–619 (1957).

    Article  CAS  Google Scholar 

25. Loane, R. F., Xu, P. & Silcox, J. Thermal vibrations in convergent-beam electron diffraction. Acta Crystallogr. A 47, 267–278 (1991).

    Article  ADS  Google Scholar 

26. Forbes, B. D., Martin, A. V., Findlay, S. D., D’Alfonso, A. J. & Allen, L. J. Quantum mechanical model for phonon excitation in electron diffraction and imaging using a Born-Oppenheimer approximation. Phys. Rev. B 82, 104103 (2010).

    Article  ADS  Google Scholar 

27. Dove, M. T. Introduction to Lattice Dynamics (Cambridge Univ. Press, 1993).

28. Van Dyck, D. Is the frozen phonon model adequate to describe inelastic phonon scattering? Ultramicroscopy 109, 677–682 (2009).

    Article  PubMed  Google Scholar 

29. Muller, D. A., Edwards, B., Kirkland, E. J. & Silcox, J. Simulation of thermal diffuse scattering including a detailed phonon dispersion curve. Ultramicroscopy 86, 371–380 (2001).

    Article  CAS  PubMed  Google Scholar 

30. Vila, F. D., Rehr, J. J., Rossner, H. H. & Krappe, H. J. Theoretical x-ray absorption Debye-Waller factors. Phys. Rev. B 76, 014301 (2007).

    Article  ADS  Google Scholar 

31. Zhang, Y. et al. Atom-by-atom imaging of moiré phasons with electron ptychography. Science 389, 423–428 (2025).

    Article  ADS  CAS  PubMed  Google Scholar 

32. Tabata, K. et al. Direct imaging of atomic rattling motion in a clathrate compound. Small Sci. 4, 2300254 (2024).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

33. Hage, F. S., Radtke, G., Kepaptsoglou, D. M., Lazzeri, M. & Ramasse, Q. M. Single-atom vibrational spectroscopy in the scanning transmission electron microscope. Science 367, 1124–1127 (2020).

    Article  ADS  CAS  PubMed  Google Scholar 

34. Leuthner, T. H., Lichte, H. & Herrmann, K.-H. STEM-holography using the electron biprism. Phys. Status Solidi A 116, 113–121 (1989).

    Article  ADS  Google Scholar 

35. Buban, J. P. et al. Grain boundary strengthening in alumina by rare earth impurities. Science 311, 212–215 (2006).

    Article  ADS  CAS  PubMed  Google Scholar 

36. Futazuka, T., Ishikawa, R., Shibata, N. & Ikuhara, Y. Grain boundary structural transformation induced by co-segregation of aliovalent dopants. Nat. Commun. 13, 5299 (2022).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

37. Zambon, P. et al. High-frame rate and high-count rate hybrid pixel detector for 4D STEM applications. Front. Phys. 11, 1308321 (2023).

    Article  Google Scholar 

38. Krajnak, M. & Etheridge, J. A symmetry-derived mechanism for atomic resolution imaging. Proc. Natl Acad. Sci. USA 117, 27805–27810 (2020).

    Article  ADS  CAS  PubMed  PubMed Central  Google Scholar 

39. Madsen, J. & Susi, T. The abTEM code: transmission electron microscopy from first principles. Open Res. Eur. 1, 24 (2021).

    Article  PubMed  PubMed Central  Google Scholar 

40. LeBeau, J. M., Findlay, S. D., Allen, L. J. & Stemmer, S. Position averaged convergent beam electron diffraction: theory and applications. Ultramicroscopy 110, 118–125 (2010).

    Article  ADS  CAS  PubMed  Google Scholar 

41. Peng, L.-M., Ren, G., Dudarev, S. L. & Whelan, M. J. Debye–Waller factors and absorptive scattering factors of elemental crystals. Acta Crystallogr. A 52, 456–470 (1996).

    Article  ADS  Google Scholar 

42. Bartók, A. P., Payne, M. C., Kondor, R. & Csányi, G. Gaussian approximation potentials: the accuracy of quantum mechanics, without the electrons. Phys. Rev. Lett. 104, 136403 (2010).

    Article  ADS  PubMed  Google Scholar 

43. Bartók, A. P., Kermode, J., Bernstein, N. & Csányi, G. Machine learning a general-purpose interatomic potential for silicon. Phys. Rev. X 8, 041048 (2018).

    Google Scholar 

44. Csányi, G. et al. Expressive programming for computational physics in Fortran 95+. Institute of Physics, Computational Physics Group. Newsletter, Spring 2007, 1–24 (2007).

45. Kermode, J. R. f90wrap: an automated tool for constructing deep Python interfaces to modern Fortran codes. J. Phys. Condens. Matter 32, 305901 (2020).

    Article  ADS  CAS  PubMed  Google Scholar 

46. Eriksson, F., Fransson, E. & Erhart, P. The Hiphive package for the extraction of high-order force constants by machine learning. Adv. Theory Simul. 2, 1800184 (2019).

    Article  Google Scholar 

47. Togo, A., Chaput, L., Tadano, T. & Tanaka, I. Implementation strategies in phonopy and phono3py. J. Phys. Condens. Matter 35, 353001 (2023).

    Article  CAS  Google Scholar 

48. Togo, A. First-principles phonon calculations with phonopy and phono3py. J. Phys. Soc. Jpn. 92, 012001 (2023).

    Article  ADS  Google Scholar 

49. West, D. & Estreicher, S. K. First-principles calculations of vibrational lifetimes and decay channels: hydrogen-related modes in Si. Phys. Rev. Lett. 96, 115504 (2006).

    Article  ADS  CAS  PubMed  Google Scholar 

50. Zhu, C., Byrd, R. H., Lu, P. & Nocedal, J. Algorithm 778: L-BFGS-B: Fortran subroutines for large-scale bound-constrained optimization. ACM Trans. Math. Softw. 23, 550–560 (1997).

    Article  MathSciNet  Google Scholar 

51. Virtanen, P. et al. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nat. Methods 17, 261–272 (2020).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

52. Tabata, K. Atomic-scale double-slit interferometry with a focused electron probe. Zenodo https://doi.org/10.5281/zenodo.18189080 (2026).

Download references

Authors and Affiliations

1. Institute of Engineering Innovation, School of Engineering, The University of Tokyo, Tokyo, Japan

   Koudai Tabata, Takehito Seki, Ryo Ishikawa & Naoya Shibata

2. PRESTO, Japan Science and Technology Agency, Kawaguchi, Japan

   Takehito Seki

3. Faculty of Physics, University of Vienna, Vienna, Austria

   Toma Susi

4. Nanostructures Research Laboratory, Japan Fine Ceramics Center, Nagoya, Japan

   Naoya Shibata

Contributions

K.T., T. Seki and N.S. designed the study and wrote the paper. K.T. performed the STEM experiments, scattering simulations and analysis. T. Susi established the full phonon-based correlation model for use in the scattering simulations. T. Susi and R.I. contributed to the discussion and comments. T. Seki and N.S. directed the entire study.

Corresponding authors

Correspondence to Takehito Seki or Naoya Shibata.

Competing interests

The authors declare no competing interests.

Peer review information

Nature thanks Xingxu Yan and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions