Capturing 3D atomic defects and phonon localization at the 2D heterostructure interface.

The three-dimensional (3D) local atomic structures and crystal defects at the interfaces of heterostructures control their electronic, magnetic, optical, catalytic, and topological quantum properties but have thus far eluded any direct experimental determination. Here, we use atomic electron tomography to determine the 3D local atomic positions at the interface of a MoS2-WSe2 heterojunction with picometer precision and correlate 3D atomic defects with localized vibrational properties at the epitaxial interface. We observe point defects, bond distortion, and atomic-scale ripples and measure the full 3D strain tensor at the heterointerface. By using the experimental 3D atomic coordinates as direct input to first-principles calculations, we reveal new phonon modes localized at the interface, which are corroborated by spatially resolved electron energy-loss spectroscopy. We expect that this work will pave the way for correlating structure-property relationships of a wide range of heterostructure interfaces at the single-atom level.

INTRODUCTION

A major challenge in materials design and engineering is to tailor the 3D atomic structures at the interface to achieve the desired properties. Although 2D lateral and vertical heterostructures have been actively studied for fundamental interest and practical applications (–), our current understanding of the atomic structure at the heterointerface has primarily relied on aberration-corrected electron microscopy and scanning probe microscopy (–), which provide either 2D projection images or surface structure. On the computational side, density functional theory (DFT) can be used to predict the properties of heterostructures (–) but requires knowledge of the 3D local atomic coordinates. Because of the difficulty in directly measuring these 3D coordinates, such studies often use perfect crystal lattices (–), statistically incorporate crystal defects into the interface (), and relax the atomic configurations to the minimum energy states. However, real heterointerfaces neither have perfect crystal lattices nor are in the minimum energy states. Here, using a MoS2-WSe2 lateral heterojunction as a model, we applied atomic electron tomography (AET) (–) to determine the 3D atomic coordinates and crystal defects at the heterointerface with picometer precision. We observed various crystal defects—including vacancies, substitutional defects, bond distortion, and atomic-scale ripples—and quantitatively characterized the 3D atomic displacements and full strain tensor across the heterointerface. The experimentally measured 3D atomic coordinates, representing a metastable state of the heterojunction, were used as direct input to first-principles calculations to reveal new phonon modes localized at the heterointerface, which were corroborated by the measurements of spatially resolved electron energy-loss spectroscopy (EELS). In contrast, the phonon dispersion derived from the minimum energy state of the heterojunction is absent of the local interface phonon modes, indicating the importance of using experimental 3D atomic coordinates as direct input to better predict the properties of heterointerfaces.

RESULTS

The experiment was conducted with an aberration-corrected scanning transmission electron microscope (STEM), operated at 60 kV in annular dark-field (ADF) mode. A tilt series of 12 images was acquired from an interface region of an epitaxial MoS2-WSe2 lateral heterojunction (see Materials and Methods, Fig. 1A, fig. S1, and tables S1 and S2). Using scanning AET (sAET) (), we determined the 3D atomic coordinates and atomic species of an interface region, containing 488 Mo, 991 S, 150 W, and 257 Se atoms and 16 S/Se vacancies (Fig. 1, B and C; fig. S2; and table S1). From the experimental 3D atomic model, we combined multislice simulations with sAET to estimate the 3D precision of the method to be 4, 15, 6, and 15 pm for the Mo, S, W, and Se atoms, respectively (see Materials and Methods). We found that there are 38 S atoms and 8 Se atoms at the interface, indicating that during the two-step epitaxial growth, the first grown WSe2 layer had W-terminals and favored W─S bonds. Most of the S and Se vacancies (9 of 16) are located along the interface, suggesting that the 2D interface lowers the vacancy formation energy. We also observed S/Se substitutional defects near the interface (Fig. 1C).

Fig. 1.

3D atomic coordinates of a monolayer MoS2-WSe2 heterostructure.

(A) An optical microscope image of the as-grown sample on a sapphire substrate. Scale bar, 2 μm. (B) Top view of the 3D AET reconstruction of an interface region in the heterostructure. The 3D atomic coordinates and species in the octagon were determined with picometer precision. Scale bar, 1 nm. (C) Experimental 3D atomic model of the square region in (B). (D) Atomic-scale ripples in the interface region, where the dots represent the Mo/W atoms. The 3D box indicates the interface region with large ripples. (E) Side view of the 3D box in (D) along the interface (i.e., the red rectangle), showing a sinusoidal oscillation at the interface. (F) Distribution of the roughness, i.e., the standard deviation of the z coordinates, within each of the nine strips (inset) as a function of the distance to the interface.

Since the lattice constant of WSe2 is 5.53% larger than that of MoS2 (), the lattice mismatch is expected to create crystal defects at the interface. Figure 1D shows atomic-scale ripples at the heterostructure interface. A side view perpendicular to the interface shows a sinusoidal oscillation of the ripples with an amplitude of ~1.5 Å in the z axis and a wavelength of ~5 nm along the interface (Fig. 1E). We divided the reconstructed region into nine strips along the interface and calculated the standard deviation (roughness) of the z coordinates within each strip. The roughness is largest on the WSe2 side of the interface, propagating into the MoS2 side, and decaying away from the interface (Fig. 1F). The atomic-scale ripples induced bond distortion at the interface, which we quantified by measuring the angles of the S─Mo─S and Se─W─Se bonds. As the interface in the reconstructed region mainly falls in the valley of the ripple (Figs. 1D and 2A), the average bond angle difference between the bottom and top layers is 8.1° for MoS2 and 6.4° for WSe2 (Fig. 2, B and C).

Fig. 2.

3D atomic displacements and strain tensor at the 2D heterostructure interface.

(A) Side view of the bond distortion along the MoS2-WSe2 interface. Distribution of the in-plane angles of the S─Mo─S (B) and Se─W─Se bonds (C). The average bond angle difference between the bottom (θ2) and top layer (θ1) is 8.1° for MoS2 and 6.4° for WSe2. (D) 3D atomic displacements of Mo/W atoms along the x, y, and z axes. (E) Six components of the strain tensor in the Mo (weak dots)/W (dark dots) layer. The strain tensor of the two S/Se layers is shown in fig. S3. Scale bar, 2 nm.

The lattice mismatch also creates atomic displacements from the perfect crystal lattice (Fig. 2D). From the displacements of all atoms, we determined the full 3D strain tensor at the heterointerface. Figure 2E and fig. S3 show the six components of the strain tensor in the three atomic layers. The in-plane components ε and ε are positive on the MoS2 side and negative on the WSe2 side, while the shear component ε is negative along the interface. This is consistent with the fact that the smaller lattice constant of MoS2 than of WSe2 induces an expansion in MoS2 and a compression in WSe2 to achieve an epitaxial interface. We also observed that ε is the largest among all six components, which is due to the lack of constraint along the z axis, and consistent with observations of the out-of-plane ripple. Although strain engineering in 2D heterointerfaces has been previously reported (, ), our results represent the first experimental measurements of the 3D strain tensor in a 2D heterostructure with near-atomic resolution, showing that the out-of-plane component ε is at least as important as the other components.

Next, we investigated the vibrational properties at the heterointerface by using the experimental 3D atomic coordinates as direct input to first-principles calculations without relaxation. Two atomic strips perpendicular to interface were chosen and stitched together to create an interface supercell, consisting of 144 atoms with periodic boundary conditions (see Materials and Methods and Fig. 3A). The phonon dispersion was calculated for the interface supercell, which was unfolded onto an effective band structure with the symmetry of the primitive hexagonal cell (Fig. 3B and fig. S4) (). The color scale represents the supercell spectral function with higher intensity arising from bulk MoS2 and WSe2 (). The presence of lower-intensity, dispersion-less modes such as those labeled by the polygon markers suggests the emergence of localized, interface modes. There are also shadow bands in the spectral function, introduced primarily by the 3D atomic defects at the heterointerface (). Prompted by the emergence of these lower-intensity, dispersion-less modes, we further investigated the phonon eigen displacements of those modes at the Γ, M, and K points, respectively, indeed revealing localized interface modes (Fig. 3, C and D). These are not present in the first-principles calculations of a perfect crystallographic heterojunction, i.e., one without point defects, and with all atoms at minimum energy configurations (fig. S4D).

Fig. 3.

Phonon dispersion and localized interface modes.

(A) Top view of the experimental 3D atomic model, where two atomic strips (red and blue) perpendicular to interface were stitched together to create an interface supercell for DFT calculations. (B) Phonon dispersion calculated from the interface supercell and unfolded onto a hexagonal primitive cell. (C) Eigen displacements of three interface modes at high symmetry points with increasing frequencies, labeled with a triangle, square, and pentagon in (B). Red and blue horizontal rectangles represent the stitched interface supercell, where the top row shows three supercells. (D) Magnified views of the eigen-displacements of three interface modes, corresponding to the vertical rectangles in (C). The displacement arrows are magnified by 10× for visual clarity. (E) Orange and blue curves represent the PDOS obtained from (B) (top) and by the least-squares fit of a linear combination of the bulk MoS2 and WSe2 dispersion (bottom), respectively. (F) Inverse PRs for the interface supercell (top) and weighted by the distance to the interface (bottom). In the latter, the emergence of high interface participation contributions at the three energy windows (labeled with a triangle, square, and pentagon) becomes apparent with intensities larger than the average (1/144, where 144 is the number of atoms in the stitched supercell.)

To quantitatively characterize the phonon dispersion, we used the least-squares fit of a linear combination of the bulk MoS2 and WSe2 phonon density of states (PDOS) (Fig. 3E, bottom). Two distinctive peaks were observed in the PDOS of the interface dispersion, namely, 17.5 and 28.5 meV, but are absent in the least-squares fit PDOS (Fig. 3E, top). We note that the absence of the localized K-mode around 41.5 meV in the PDOS is due to the fact that PDOS is, by definition, momentum-averaged, and subtle features can be obfuscated by bulk modes with different momentum transfers. To further analyze the phonon localization, we computed the participation ratios (PRs), a measure of the deviation away from plane-wave like phonon vibrations expected in perfect crystal lattices (). Figure 3F shows the PR as a function of phonon frequency. The phonon modes with PR smaller than unity illustrate the disordered nature introduced by the 3D atomic defects at the interface (top). We correlated this localization with the interface by computing an interface-weighted PR (see Materials and Methods). Increased density of phonon modes with interface-weighted PR above the average (dashed line), labeled with a triangle, square, and pentagon in Fig. 3F (bottom), further confirm the existence of interface-localized phonon modes (Fig. 3, B to D).

To experimentally corroborate the interface phonon modes, we measured the vibrational spectra of the MoS2-WSe2 heterostructure using spatially resolved EELS with an energy resolution of 5.7 meV at 60 keV (see Materials and Methods and fig. S5) (, ). By scanning continuously from a MoS2 to a WSe2 region (Fig. 4, A and B), we collected a series of vibrational spectra across the heterointerface (Fig. 4C), which reflect the PDOS throughout the entire Brillouin zone (, ). The sample remains the same after the EELS data collection (fig. S6). On the basis of the different features at the MoS2, interface, and WSe2 regions, the vibrational spectra were divided into three groups via k-means clustering (fig. S8). To enhance the signal-to-noise ratio, the spectra in each group were summed up to produce an average spectrum of which the peak positions were fit with multiple Gaussians (Fig. 4D, fig. S9, and table S3). The peak positions of the average spectra in the MoS2 and WSe2 regions agree with those of previously reported PDOS curves for bulk MoS2 and WSe2 (table S3). However, the average interface spectrum shows two peaks at 27.9 and 41.1 meV, which are absent in the PDOS curves of bulk MoS2 and WSe2. Both peaks are consistent with the second and third interface modes of 28.5 and 41.5 meV, obtained by DFT calculations with the experimental 3D atomic coordinates (Fig. 3). The first interface mode of 17.5 meV in the DFT calculations could not be detected in the vibrational spectrum due to the influence of the zero-loss peak (ZLP) tail. To verify the generality of our experimental measurements and analysis, we acquired the vibrational spectrum from another MoS2-WSe2 interface, exhibiting similar results (fig. S10 and table S3).

Fig. 4.

Vibrational spectra of the MoS2-WSe2 heterointerface.

(A) Low-magnification ADF-STEM image of the heterointerface. Scale bar, 10 nm. (B) Magnified ADF-STEM image of the orange box in (A), showing the local atomic structure at the interface. Scale bar, 2 nm. (C) Vibrational spectra measured from scanning the MoS2 region (green) through the interface (orange) to the WSe2 region (brown) with a step size of 1.6 nm. The background was subtracted by fitting a power function to each spectrum (fig. S7A). K-means clustering was used to separate the spectra in the three regions (see Materials and Methods). Our experimental results show that the residual polymers do not affect the phonon spectra or the peak assignments (fig. S7B). (D) Average vibrational spectra of the WSe2, interface, and MoS2 region, overlapped with fitted Gaussian peaks and peak positions (fig. S9). The shaded Gaussian peaks in WSe2 (brown) and MoS2 (green) spectra represent the bulk phonon modes based on the comparison with corresponding DFT results. The interface phonon modes are shaded in the interface spectrum (orange). All nonshaded peaks are caused by the delocalization effect of low energy loss excitation (see more details in table S3).

DISCUSSION

The observation of localized interface phonon modes has implications on the thermal transport mechanism across the interface. When phonons propagate across an interface of two materials with different vibrational properties, phonons undergo additional scattering, which impedes heat transport across the interface (–). In the case of the MoS2-WSe2 heterojunction, the difference in PDOS indicates that several optical phonon modes in WSe2 cannot directly propagate into MoS2 because of the requirement of conservation of energy and momentum. As our experimental-based first-principles calculations suggest, the emergence of localized modes at the interface leads to increased scattering and thus an increased thermal interface resistance across the heterostructure. Furthermore, our results show that the interface phonon modes cannot be obtained by a linear combination of the bulk phonon modes of two different materials forming the interface but rather depend on the 3D local atomic structure and crystal defects at the interface, including point defects, bond distortions, atomic-scale ripples, and strain. Thus, atomic-scale design and engineering of interface phonon modes could potentially improve the thermal transport in 2D and 3D heterostructures, mitigate the heat dissipation, and increase the lifetime of modern microelectronic devices.

Using a monolayer MoS2-WSe2 heterojunction as a model, we demonstrated a correlative experimental and first-principles method to determine the 3D atomic coordinates and crystal defects and to capture the localized vibrational properties at the epitaxial interface. We observed the local atomic structure reconstruction of the interface by means of out-of-plane ripples and bond distortion to accommodate the lattice mismatch. We measured the full 3D strain tensor at the interface, showing that ε is the largest among the six components. The experimentally measured 3D atomic coordinates, representing a metastable state of the heterojunction, were directly used in first-principles calculations of vibrational properties. The unfolded phonon dispersion relation and the PDOS show distinctive phonon modes emerging at the interface, which were verified by spatially resolved EELS. Although, in this study, we revealed 3D atomic defects and local interface phonon modes in a MoS2-WSe2 heterojunction, we expect that the correlative experimental and first-principles method can be used to probe the structure property relationship of a broad class of heterostructure interfaces at the single-atom level.

MATERIALS AND METHODS

Sample preparation

Single-crystal WSe2 monolayer was first grown by the chemical vapor deposition method. WO3 power of 0.6 g was placed in a quartz boat located in the heating zone center of the furnace. A sapphire substrate was placed at the downstream side, just next to the quartz boat. The Se powders were placed in a separate quartz boat at the upper stream side of the furnace, and the local temperature was maintained at 260°C during the reaction. The gas flow was brought by a mixture of Ar [90 standard cubic centimeter per minute (sccm)] and H2 (6 sccm) with a chamber pressure of 20 torr. After reaching the required growth temperature of 925°C, the heating was kept for 15 min, and the furnace was then naturally cooled down to room temperature. After optical characterizations for the as-grown WSe2, the sample was then put into a separate furnace for the second step of MoS2 growth. The setup for MoS2 growth is similar to that of WSe2, by switching the source to MoO3 power (0.6 g) and S powers. The Ar gas flow was set at 70 sccm, and the pressure was controlled at 40 torr. The sapphire substrate with WSe2 sample was placed at the downstream side of the MoO3 boat, and the distance between the sample and quartz boat was 9 cm to achieve the best Mo and S source ratio to form the WSe2/MoS2 heterojunction. The center zone and S source region were heated to 755° and 190°C, respectively, and held for 15 min for synthesis and then naturally cooled down to room temperature. The as-grown WSe2-MoS2 heterojunction was then transferred onto a QUANTIFOIL holey carbon film TEM grid by a poly(methyl methacrylate) (PMMA) (950 PMMA A4, MicroChem)–assisted transfer method. PMMA thin film was spin-coated on top of the sample, and then the PMMA/sample/sapphire was dipped in a 6 M hydrogen fluoride solution to etch the sapphire substrate. PMMA/sample was lifted from the etching solution, diluted and cleaned in deionized water, and then transferred onto the TEM grid. The PMMA layer was rinsed out by acetone and isopropanol. The as-prepared TEM grid was vacuum annealed at 160°C for more than 8 hours to remove most of the residual polymer before loading into the electron microscope.

Tomographic data acquisition of a MoS2-WSe2 heterojunction

The AET experiment was performed on a MoS2-WSe2 heterojunction using a Nion UltraSTEM 200 equipped with C3/C5 corrector and high-energy resolution monochromated EELS system. ADF-STEM images were acquired at 60 kV with an averaged beam current of 17 pA and a probe size of 1.3 Å. The convergence semiangle of electron probe was 38 mrad, and the inner and outer collection semiangles of ADF detector were 50 and 210 mrad, respectively, where a small detector inner angle was chosen to reduce the electron damage to the sample. Table S1 summarizes the experimental parameters. A tilt series of 12 angles were acquired from the heterojunction (table S2). To reduce the total electron dose and eliminate the sample drift issue, 10 images per tilt angle were collected with a dwell time of 4 μs per pixel and 512 × 512 pixels for each image. Because of the use of a low-voltage electron beam and a low-dose acquisition scheme, there was minimal structural change of the heterojunction during the data acquisition (fig. S1). The total electron dose for the whole tilt series was estimated to be 3.8 × 105 e/Å2.

Drift correction, denoising, and background subtraction

Drift correction was applied to each image as described elsewhere (–). At each tilt angle, 10 images were registered, and the remaining 9 images were aligned to the first image with cross correlation to subpixel precision. The drift between two neighboring images was typically less than one pixel, and the drift vector followed a linear trajectory. As the drift vector was uniformly distributed along the slow scan direction during data acquisition, the scan distortion was implemented by interpolating the noncorrected image onto the drift-corrected pixel positions. After the drift correction, the 10 images were summed to obtain an average image for the tilt angle. As each image contains mixed Poisson and Gaussian noise, a block-matching and 3D filtering algorithm was applied to denoise the average image of each tilt angle (). From the denoised image, the background was estimated on the basis of the fact that the center of the hexagonal atomic lattice should be empty. The local minimum of each hexagonal lattice was chosen to form a sparse map, and the background was obtained by interpolating the local minimum map with a radial basis function (). The estimated background was subtracted from the denoised image.

Angle calibration, vibration correction, and image alignment

The nominal tilt angles were recorded during data acquisition (table S2) but must be calibrated because of the hysteresis of the stage. For each image, the x and y coordinates of the W and Mo atoms were determined by fitting a 2D Gaussian. To calibrate the tilt angles, a least-square method was used to minimize E (), defined aswhere P and P are the functions of projecting the 3D coordinates of the W and Mo atoms to the x and y coordinates in the image, respectively; is the 3D coordinates of the ith W or Mo atom; α and β are the tilt angles of the jth projection; and and are the measured x and y coordinates in the ith atom in the jth projection, respectively. By minimizing E, we calibrated the tilt angles of the 12 images (table S2).

Because of the free-standing nature of the monolayer 2D heterojunction during data acquisition, we observed that the 2D heterojunction vibrated in the vertical direction. As a result, the atoms in the ADF-STEM images are elongated and blurred at relatively high tilt angle. The vibration blurring is a type of motion blurry in image processing and can be corrected by deconvolution if the vibration kernel is known. In this experiment, we estimated the vibration kernel based on the fact that all atoms should be spherical. First, each experimental image was interpolated by a factor of 2 with linear interpolation, from which a high-quality MoS2 region of 100 × 100 pixels was cropped. Second, a MoS2 atomic model was used to generate reference images at 12 tilt angles. Third, a harmonic oscillator kernel was constructed of which the vibration direction, amplitude, and width were optimized by using the Lucy-Richardson algorithm (, ) to match the cropped experimental images with the reference ones. Last, the optimized kernel was applied to the whole experimental image, and the image was then binned to its original size. The procedure was repeated for all the 12 experimental images. After vibration correction, the 12 images were aligned with the center of the mass methods (, ). Using the x and y coordinates of the W and Mo atoms in the region of interest determined during angle calibration, the center of mass of all W and Mo atoms in each image was aligned to the center of mass of the image with subpixel precision.

3D reconstruction with sAET

We implemented sAET reconstruction as described elsewhere (). Briefly, a 3D window of 60 × 60 × 60 voxels was chosen and scanned along the x and y axes with a step size of 30 voxels. At each step, the corresponding regions in all 12 images were cropped and grouped into an image stack. Each image stack consists of 12 images with varied shapes, corresponding to the projection of the 3D window along different tilt angles. After a full 2D scan was completed, all the images were partitioned into image stacks. The image stacks were aligned and reconstructed in parallel by the generalized Fourier iterative reconstruction (GENFIRE) algorithm (, ). Each GENFIRE reconstruction used a 33-voxel support along the z axis and ran 1000 iterations. To remove these artifacts, we stitched together only the central 30 × 30 × 33 voxels of the reconstructed windows to produce a full 3D reconstruction.

Determination of 3D atomic coordinates and chemical species

The 3D atomic coordinates and chemical species of the 2D heterojunction were traced from the 3D reconstructions with a polynomial fitting method (, , ). The reconstruction was first interpolated onto a finer grid with 3× oversampling using the spline method. All local maxima in the reconstruction were identified, and the positions of potential atoms were extracted with a 3D volume (5 × 5 × 5 voxels). For every potential atom, a minimum distance of 1.6 Å to its neighboring atoms was enforced. These positions were manually checked to correct for unidentified or misidentified atoms due to fitting failure or artifacts. We then assigned the initial chemical species based on the MoS2 and WSe2 region and the intensity of the atoms.

The traced 3D atomic coordinates and chemical species were refined by a gradient descent method as described elsewhere (–). The 3D atomic coordinates in each model were refined by minimizing the error between the experimental images and those calculated from the model. From the refined 3D atomic coordinates, we applied an atom pair flipping method to identify the S/Se vacancies (, ), which consists of the following four steps. First, we randomly chose a pair of S/Se atoms between the top and bottom atomic layers. For each selected S/Se pair, projection images were calculated for all the 12 tilt angles by flipping the pair among nine cases: (i) both S atoms, (ii) both vacancies, (iii) both Se atoms, (iv) top S and bottom Se atom, (v) top Se and bottom S atom, (vi) top S atom and bottom vacancy, (vii) top vacancy and bottom S atom, (viii) top Se atom and bottom vacancy, and (ix) top vacancy and bottom Se atom. Nine atomic models were generated for the nine cases. Second, nine sets of 12 images with the experimental tilt angles were calculated from the atomic models. An R1 factor was computed between measured and calculated images (–). By comparing the R1 factor among the nine cases, the one with the smallest error R1 was chosen and updated in the atomic model. Third, we repeated steps one and two for all the S/Se atom pairs and obtained an updated 3D atomic model. Last, we iterated steps one to three for all the S/Se atom pairs until there was no further change in the chemical species.

Because of the use of a low-voltage electron beam (60 keV) and a small detector inner angle, we observed the dynamic scattering effect of the W atoms, which influenced some nearby S atoms. We identified 19 outlier S atoms that violated a minimum distance of 1.6 Å (during refinement) or deviated from the S atomic layer by more than 1 Å. For each outlier S atom, the coordinates of the S atom were scanned with a range of ±40 pm along the x and y axes and ±120 pm along the z axis. The step size was 20, 20, and 30 pm along the x, y, and z axes, respectively. At each scanning step, we used multislice simulations to calculate 12 images of 60 × 60 pixels in size at different tilt angles and computed the R1 factor between the calculated and measured images (–). The S atom in the atomic model was updated to the position corresponding to the minimum R1 factor. This dynamic refinement procedure was repeated for all outlier S atoms.

Multislice simulations for 3D precision estimation

We performed multislice simulations to estimate the 3D precision of the atomic coordinates. Twelve multislice images were computed from the refined atomic model using the same experimental parameters (table S1). Each multislice image was convolved with a Gaussian kernel (σ = 1.1) to account for the electron probe size, thermal vibrations, and other incoherent effects. The mixture of Gaussian and Poisson noise determined from the experimental images was added to the multislice images. From the 12 multislice images, we used the same imaging processing, sAET reconstruction, atom tracing, and refinement procedures to obtain a new 3D atomic model. By comparing it with the experimental 3D atomic model, we found that all the chemical species were correctly identified, including all S vacancies. The root mean square deviation of the S, Mo, Se, and W atoms was estimated to be 15, 4, 15, and 6 pm, respectively.

Measurements of the 3D strain tensor in the MoS2-WSe2 heterojunction

The strain tensor was determined from the experimental 3D atomic coordinates using a procedure described elsewhere (). The experimental atomic model was separated to MoS2 and WSe2 regions, and each region was aligned to an ideal MoS2-WSe2 atomic model with a least-squares fit. The displacement vectors were calculated as the difference in the atomic positions between the experimental and ideal model and interpolated to a cubic grid using a radial basis function (). A 3D Gaussian kernel with σ = 3.16 Å was convolved with the displacement field to increase the precision of the strain tensor measurement.

First-principles calculations using experimental 3D atomic coordinates without relaxation

To investigate the vibrational properties of the lateral heterojunction using first-principles calculations, a periodic supercell across the interface is required. The reconstructed region consists of two crystallographically distinct interfaces with Pmc2 (armchair) and Pmm2 (zigzag) space groups, respectively. We constructed a periodic supercell by stitching together two regions across a zigzag interface. To this end, imagine starting in the MoS2 region and moving perpendicularly toward the middle zigzag interface and into the WSe2 region, stopping at a heavy (W) atom (Fig. 3A, red box). Having moved sufficiently far into the WSe2 region, the effects of the interface on the atomic coordinates should be minimal, and thus, we assumed that we could equivalently use the displacement of a different W atom, also sufficiently far from the interface. We chose a different W atom and moved perpendicularly toward the bottom zigzag interface into the MoS2 region (Fig. 3A, blue box).

We avoided relaxing the experimental coordinates to extract vibrational properties by sampling the Born-Oppenheimer energy surface, around the experimental coordinate local minimum, using ab initio molecular dynamics (AIMD). In particular, we used the Vienna ab initio simulation package (, ) for AIMD and extract the harmonic force constants using the ALM code (), as implemented in phonopy (). We used a cutoff energy of 500 eV on a Monkhorst-Pack k-points grid of 8 × 8 × 1.

The frequency-dependent PRs and interface-weighted PRs were evaluated usingwhere eα(q) and Mα are the eigenvectors and atomic masses of species α, respectively (). The interface weight wα is a normally distributed weight of the distance of an atom in the supercell from the interface, normalized to sum up to one interface, thus enabling the 1/144 comparison in the work.

To capture the influence of the 3D atomic defects on lattice vibrations, the experimental atomic coordinates of the stitched supercell were not relaxed. Figure S4 (C and D) compares the unfolded phonon dispersions against bulk phonon dispersions and the unfolded phonon dispersion of perfect crystallographic zigzag interface, i.e., one where both sides have been independently relaxed to their respective global minimum. While the perfect crystallographic PDOS is more consistent with a linear combination of bulk MoS2 and WSe2, the PDOS directly calculated from the experimental atomic coordinates shows localized interface modes, corroborated by spatially resolved EELS (see section below). This study highlights the importance of using the experimental coordinates as direct input to first-principles calculations. What is more, the increased density of soft (imaginary-frequency) modes suggests that the out-of-plane ripple and structural reconstructions present in the stitched supercell stabilize the structure’s lattice vibrations (Fig. 3B).

Spatially resolved vibrational EELS

The vibrational spectra were measured by the same Nion microscope at 60 kV with a convergence semiangle of 33 mrad and a probe current of 100 pA, which provides a spatial resolution of 1.5 Å. This higher probe current is beneficial to the acquisition of the vibrational spectra. ADF-STEM images of Fig. 4 (A and B) were acquired with an ADF collection semiangle ranging from 70 to 210 mrad. With the implementation of an alpha-type monochromator and newly designed EELS spectrometer, the best energy resolution of the EELS is 5.7 meV with a short exposure time of 3 ms. To obtain a sufficient vibrational signal from the 2D heterojunction and avoid the saturation of ZLP, an exposure time of 300 ms per frame was used, which corresponds to an energy resolution of 8.3 meV (fig. S5). The EELS collection semiangle was set to be 20 mrad, which is larger than the size of Brillouin zone of either MoS2 or WSe2 along the [001] direction. The EELS dispersion was 0.44 meV per channel. The high-resolution vibrational spectra (Fig. 4) were collected by running our custom-developed Python script on top of Nion’s Swift software (, ). A 40-point line-scan experiment was performed with a step size of 1.6 nm by aligning and summing 200 frames of 300-ms exposures per point. EELS acquisition was conducted under low-magnification condition, where the ADF-STEM images acquired before and after EELS acquisition show consistent material morphology (fig. S6). Background was subtracted from the vibrational spectra by fitting a power function to two fitting windows of 16 to 20 meV and 60 to 70 meV (fig. S7). A k-means clustering algorithm () was used to divide the vibrational spectra into three groups (#1 to 16 for MoS2, #17 to 26 for interface, and #27 to 40 for WSe2), where the definition of “distance” (the minimization parameter) was one minus the correlation between the background-subtracted spectra and the centroid. Each centroid was the component-wise mean of the points in that cluster after centering and normalizing all spectra to zero mean and unit standard deviation. The process of minimization increased the correlation between each spectrum and its corresponding centroid. For the peak separation of the averaged vibrational spectra, the Levenberg-Marquardt iteration algorithm () was applied to fit spectra by five Gaussian-type peaks without any bounds on all fitting parameters.