# Fileset

[Manuscript.pdf](https://mdr.nims.go.jp/filesets/c0e5bf60-299e-4e4c-8f95-a51b6c122415/download)

## Creator

Rina Mishima, [Takuro Nagai](https://orcid.org/0000-0001-5239-3334), [Hiroyo Segawa](https://orcid.org/0000-0002-7198-8410), Masahiro Ehara, Takashi Uchino

## Rights

This document is the Accepted Manuscript version of a Published Work that appeared in final form in JOURNAL OF PHYSICAL CHEMISTRY C, copyright © 2025 American Chemical Society after peer review and technical editing by the publisher. To access the final edited and published work see https://doi.org/10.1021/acs.jpcc.5c00804.[In Copyright](http://rightsstatements.org/vocab/InC/1.0/)

## Other metadata

[Atomic-Scale Observation of Moiré potential in Twisted Hexagonal Boron Nitride Layers by Electron Microscopy](https://mdr.nims.go.jp/datasets/97c3003f-9a7e-4d52-85a8-8d5c64717290)

## Fulltext

1  Atomic-Scale Observation of Moiré potential in Twisted Hexagonal Boron Nitride Layers by Electron Microscopy Rina Mishima,†Takuro Nagai,§ Hiroyo Segawa, Masahiro Ehara,‡ and Takashi Uchino†* †Department of Chemistry, Graduate School of Science, Kobe University, Nada, Kobe 657-8501, Japan §Research Network and Facility Services Division, National Institute for Materials Science, Tsukuba, Ibaraki 305-0044, Japan ǁResearch Center for Electronic and Optical Materials, National Institute for Materials Science, Tsukuba, Ibaraki 305-0044, Japan ‡Research Center for Computational Science, Institute for Molecular Science, Okazaki, Aichi, 444-8585, Japan  ABSTRACT: Moiré superlattices (MSLs) are an emerging class of two-dimensional functional materials whose electronic states can be tuned by the twist angle between two van der Waals layers and/or the relative placement of the layers. The intriguing properties of MSLs are closely correlated to the moiré potential, which is the electrostatic potential induced by interlayer coupling. Intensive efforts have been made to understand the nature and distribution of the moiré potential by using various experimental and theoretical techniques. However, the experimental observation of the moiré potential is still challenging because of the possible presence of the surface and/or interlayer contaminants. In this work, we develop a method to obtain hexagonal boron nitride (hBN) nanolayers (with or without twist) using a specially designed chemical exfoliation technique. The resulting hBN nanolayers are atomically clean and strain free, hence providing ideal MSLs for the investigation of their moiré potential. Aberration-corrected high resolution transmission electron microscopy measurements on the twisted hBN nanolayers allow us to observe moiré diffraction spots in Fourier space. Then, the moiré potential is reconstructed by the inverse fast Fourier transform of the moiré diffraction spots. It has been revealed that the local interlayer atomic overlap plays a decisive role in determining the periodicity and distribution of the moiré potential, as supported by density functional theory calculations. This work not only provides a general strategy to observe the moiré potential in MSLs, but it also expands the application of electron microscopy to the further study of MSLs with atomic resolution.  ■ INTRODUCTION   Vertically stacked twisted layers of two-dimensional (2D) materials provide a fruitful ground for the exploration of unique and novel quantum phenomena in van der Waals (vdW) material systems.114 Stacking of 2D vdW materials with a small twist angle θ (θ<~10°) and/or a slight lattice mismatch results in the reconstructed moiré superlattices (MSLs), accompanied by the enhanced interlayer cou-pling.4,1416 As for twisted graphene layers, these recon-structed MSLs especially with magic twist angles around θ≈1.1º can host flat electronic bands to form highly corre-lated phases, even showing superconductivity.1,8,9  In the small angle regime of MSLs, the commensuration cell lattice constant can be equal to the moiré periodicity D, which is given by15                           𝐷 =  𝑎2 sin𝜃2     (1),  where a is the lattice constant of the constituent layer. For large twist angles (θ>~15°), on the other hand, the commen-suration cell lattice constant is generally much greater than D.15 However, there exist a few special twist angles at, e.g., 13.2° and 21.8º for twisted hexagonal bilayers,15,17,18 where the “commensuration periodicityˮ is equal to D. Thus, from the view point of MSLs, any large twist angles except for the above special twist angles can be regarded as incommensu-rate twist angles. In the large angle regime, it is believed that atomic reconstruction hardly occurs, meaning that a rigid lat-tice model is a reasonable representation of the actual sys-tem.4 Although the interlayer coupling is supposed to be sup-pressed in these incommensurate MSLs,19 it has recently been revealed that this is not necessarily the case for 30°-twisted vdW bilayers, forming quasicrystalline patterns with a 12-fold rotational symmetry.2023 In these quasicrystalline layers, two twisted layers are coupled via Umklapp electron-2  electron scattering, which results in rich electronic structures, such as mirrored Dirac cones in graphene quasicrystals20 and mini-gaps near the valence band maximum in tungsten diselenide quasicrystals.23 Also, it has recently been demon-strated that boron nitride twisted bilayers exhibit flat bands in the quasicrystal limit.24 Thus, large twist-angle MSLs and related moiré potentials have attracted increasing attention in various MSLs,24,25,26,27 offering an interesting experi-mental platform to explore moiré physics beyond the small twist-angle regime. However, atomic-scale visualization of moiré potential is highly technically demanding. At present, scanning tunneling microscopy/spectroscopy (STM/STS) is one possible technique available for that purpose.23,28,29 Here we use an aberration-corrected high-resolution transmission electron microscopy (HRTEM) technique as a tool to explore the moiré potential of twisted hexagonal bo-ron nitride (hBN) nanolayers in atomic scale. It has been well recognized that hBN nanostructures are characterized by their excellent thermal and chemical stability and unique electronic and optical properties,30 providing an ideal build-ing block for the investigation of vdW potentials.31 Although HRTEM has been often employed to obtain atomic resolu-tion imaging of 2D materials and related MSLs,3237 we show that this technique is also applicable for imaging the moiré potentials of twisted hBN nanolayers by using their moiré diffraction spots in Fourier space. In previous works, the dif-fraction spots due to moiré potential have not been explicitly observed by HRTEM because of the general presence of hy-drocarbon contaminants which are often trapped during the assembly of 2D stacked materials.3840 These contaminants yield a severe background noise in Fourier space, obscuring weak moiré diffraction spots. Furthermore, any interlayer contaminants will modify the interlayer distance, leading to a fluctuation of moiré potential.41,42 In this work, we over-come the problem by developing a specially designed chem-ical exfoliation method to obtain atomically clean hBN nanolayers.      ■ EXPERIMENTAL SECTION Our exfoliation procedure is based on an intercalation-based method;43 however, unlike conventional methods, it does not require any sonication, sharing nor electrochemical process as a driving force, hence yielding mechanically intact and strain-free samples. Briefly, we start by creating hBN/H2SO4 intercalation compounds by heating hBN powders in an aqueous H2SO4 solution.44,45 The resulting solution is neu-tralized by a solution of NaHCO3 to form sodium sulphate salt in between the hBN layers. We found that this acid/base reaction leads to spontaneous exfoliation of high-quality hBN layers suitable for the characterization of moiré poten-tials by HRTEM (see the Supporting Information and Fig. S1 for details). HRTEM and scanning TEM (STEM) observations were carried out on a JEOL JEM-ARM300F instrument equipped with a cold field emission electron gun and a spherical aber-ration (Cs) corrector at an acceleration voltage of 80 kV, un-der 1 × 10−5 Pa in the specimen column. We analyzed strain of the exfoliated hBN layers based on the HRTEM images using the Peak Pairs Analysis (PPA) software package for Gatan Digital Micrograph.46 In order to calculate the interlayer electrostatic potential, we performed quantum chemical calculations using density functional theory (DFT) methods implemented in the Gauss-ian 16 suite of programs.47 In this work, hydrogen terminated clusters with various sizes and two different edges (zigzag and armchair) were used to model the local structure of twisted hBN bilayers with different twist angles. The morphology of the exfoliated layers was observed using a dynamic force mode atomic force microscope (AFM; SII SPA400). Further details about the characterization and calcula-tional procedures used in the work can be found in the Sup-porting Information.  ■ RESULTS AND DISCUSSION General characteristics of exfoliated hBN.  Figure 1 shows a typical AFM image of the exfoliated hBN sheets with a lateral size of a few m. Each sheet has a relatively smooth surface with a thickness t of 24 nm. Considering that the interlayer distance of hBN is ~0.33 nm, we can esti-mate that each nanosheet consists of ~515 BN layers. We show in Fig. 2a a low-magnification annular dark field scan-ning TEM (ADF-STEM) image of the exfoliated hBN, showing a random stacking of micrometer-sized hBN layers. A high-magnification ADF-STEM image of non-twisted hBN layers is given in Fig. 2b, illustrating 2D honeycomb structures of BN with asymmetric ADF scattering intensities. One may expect that three boron and three nitrogen atoms in each hexagonal ring are responsible for this asymmetric con-trast, which can be in principle distinguished by their inten-sity based on the Z-(atomic number) contrast principle.48 However, we should note that our exfoliated hBN nanosheets are likely to have more than 5 AA-stacked BN layers, which are formed by stacking many anti-aligned (B  Figure 1. (Left panel) A typical AFM image of the exfoliated hBN sheets. (Right panels) The corresponding thickness profiles along the lines denoted by (a), (b) and (c) in the left panel..  3  to N and N to B) layers. In such layers, the expected Z-con-trast will be averaged out to yield symmetric ADF scattering intensities.49 We hence assume that observed asymmetric ADF contrast results from a slight specimen tilt. To confirm the assumption, we performed simulations by changing x-tilt angle φ  (for details, see Fig. S2) using a comprehensive multipurpose crystallographic program Recipro.50 After sev-eral try-and-error attempts, we found that the specimen of t = ~8 nm and φ = ~1° yields the asymmetric contrast which is most comparable to that of the observed ADF-STEM im-age (see Fig. 2c). Although the thickness estimated from the simulation is somewhat greater than that obtained from an AFM image shown in Fig. 1, these estimated values are rea-sonably consistent with each other. Thus, we can conclude that the thickness of our exfoliated hBN sheets lies in the range from ~3 to ~8 nm, which will be thin enough to satisfy the weak phase object approximation (WPOA) especially for the samples consisting of light elements such as B and N. Under WPOA, the TEM image operated in the phase con-trast mode at optimum focus directly reveals the projected potential.51  HRTEM images shown in Fig. 3a,b demonstrate that the resulting layers have a clean and adsorbates-free area ex-tending several tens of nanometers (see also Fig. S3). The corresponding fast Fourier transform (FFT) image yields six spots of 0.21, 0.12 and 0.11 nm spacings attributed, respec-tively, to 101̅0, 112̅0 and 202̅0 Bragg peaks (see the inset in Fig. 3a). Figure 3c is the PPA strain mapping obtained for the HRTEM image shown in Fig. 3b. As shown in Fig. 3c, the strain tensor is estimated to be as low as 1 % over the entire region of interest (see Fig. S4 for more details of the PPA).  Figure 4a shows the edge area of exfoliated BN layers (see also the line profile in Fig. 4b). In general, the edge of 2D layers is folded or scrolled during the TEM specimen preparation process. As a result, c-planes become parallel to  Figure 2. (a) Low- and (b) high-magnification ADF-STEM images of exfoliated hBN layers. (c) A simulated ADF-STEM image of a 8 nm thick hBN, which is x-tilted by 1° (17.5 mrad).  Figure 3. (a) An HRTEM image showing a dark contrast for the B and N atoms. The corresponding FFT image is given in inset. (b) A magnified HRTEM image of the area marked by the white, dashed square in (a). (c) Shear strain component εxy of the lattice strain tensor obtained by applying peak-pairs analysis to the HRTEM image shown in (b). The map displays strain ranging from 20% to 20 %.   Figure 4. (a) An HRTEM image of the edge of an hBN multilayer. (b) Line profile as a function of distance along the yellow line in (a), showing a four-layer structure with a fringe separation of 0.325 nm. 4  the incident beam and appear as a series of lines in the TEM image. In Fig. 4a,b, one sees four lines with a separation of 0.325 nm, corresponding to the interlayer distance of hBN. Thus, the present exfoliated specimen consists of four BN layers. Note also that one cannot recognize any moiré pat-terns nor the related rotational stacking faults at the edge, but finds a well-ordered honeycomb structure inside the edge re-gion. This indicates that the edges are not folded back but rather exhibit the smallest possible scroll, as often observed in the edge structures of high-quality free-standing graphene layers.52  All the above results demonstrate that as far as the non-twisted hBN layers are concerned, atomically clean, flat and strain-free thin samples can be obtained by the present exfo-liation procedure. TEM observations on twisted hBN layers. In ad-dition to the non-twisted regions, we also found several con-taminant-free twisted regions at different twist angles. We show in Fig. 5a-d the HRTEM and the corresponding FFT images of the twisted hBN layers with four different twist angles of = 24.0º, 21.6º, 19.8º and 13.6 º. In this work, the twist angle was defined as the angle between the two clos-est 101̅0 spots, as often done in previous studies. The uncer-tainty of is within the range of ±0.4º. Note, however, that the alternative twist angle 𝜃̅ can be defined as the angle be-tween the two next-closest 101̅0 spots, i.e., 𝜃̅ = 60° − 𝜃 . Since we cannot identify which definition is appropriate from the diffraction pattern alone, the selection of is just for convenience. We will show that this selection is not cru-cial for the discussion of this work in a later subsection..  The observed HRTEM images illustrate a variety of moiré patterns typical to twisted bilayers with large twist an-gles. The two of the observed twist angles (13.6 ° and 21.6°) are comparable to two of the commensurate angles of twisted hexagonal bilayers (13.2° and 21.8°)17,18 within the experimental uncertainty. Thus, in what follows, we will re-fer to the samples with  =13.6 ° and 21.6° as nearly com-mensurate layers, whereas the rest of the samples as incom-mensurate ones. Note, however, that even in the incommen-surate cases, the partially eclipsed interlayer atomic overlap, i.e., N nearly on N (or B nearly on B) configurations, is ex-pected to be realized in the rigid lattice picture, as illustrated in Fig. S5. The FFT images given in Fig. 5 are especially worth not-ing. We see that in addition to the two sets of 101̅0 Bragg spots defining a twist angle, there exist several extra spots especially in the inner region of the 101̅0 spots, which are located symmetrically with respect to the center spot. These extra spots have not been explicitly reported in previous  Figure 5. HRTEM observations on twisted hBN layers with different twist angles. (a-d) (left panel) HRTEM image. Scale bar, 1 nm. (middle panel) FFT image, showing the twist angle determined from the angle between the two closer spots. Red and yellow circles show two sets of six 101̅0 diffraction spots belonging respectively to different layers. (right panel) A magnified FFT image of the area marked by the white dashed square shown in the middle panel.   Figure 6. (a) HRTEM image of twisted hBN layers with hydrocarbon adsorbates. Scale bar, 1 nm. (b) Correspond-ing FTT image, along with the twist angle determined from the angle between the two closest spots. Red and yellow circles show two sets of six 101̅0 diffraction spots belong-ing to different layers.  5  HRTEM observations on MSLs.3237 We also found that these extra spots are not observed in the TEM image of the twisted hBN layers with hydrocarbon adsorbates (Fig. 6). This allows us to expect that clean and contaminant-free samples are necessary for the observation of the extra dif-fraction spots.  What is the origin of these extra diffraction spots? One immediate answer is double diffraction, which often occurs in thick two-phase materials when a diffracted beam travel-ling through a crystal is rediffracted when it passes into a second crystal.53 Hence, it would be natural to expect that double diffraction induces several extra diffraction spots in MSLs. We, however, argue that the simple double diffrac-tion process cannot account for the extra diffraction spots shown in Fig. 5. Double-diffraction spots occur around each of the primary reflections, including  the direct beam spot, by keeping the symmetry of the constituent crystals.53 This general feature of double diffraction implies that as for twisted hBN layers, two types of primary reflections, each with 6-fold rotational symmetry, should be accompanied by a set of  double diffraction spots with the same 6-fold rota-tional symmetry, as indeed observed in 10° twist angle dou-ble-multilayer of ~100-nm-thick hBN sample.54 Contrary to the above expectation from double diffraction, the extra dif-fraction spots given in Fig. 5 does not show the 6-fold rota-tional symmetry except for the nearly commensurate sample  at  = 21.6°. Furthermore, for very thin specimens in which the WPOA holds, the double diffraction does not contribute to the entire diffraction process, as will be discussed below. For a better understanding of the extra diffraction spots, we analyzed the positions of these spots ke in reciprocal space. It has been revealed that ke can be found by connect-ing the two Bragg spots belonging to the different layers and translating the obtained vectors to the center of the diffrac-tion pattern, as demonstrated in Fig. S6. The relationship can be written as ke = G1 G2, where Gi (i = 1, 2) represents the reciprocal vectors from two different layers. This relation-ship implies that the interlayer interactions are responsible for the generation of these extra diffraction spots indeed, but the question is, why do particular sets of the reciprocal vec-tors contribute to the extra diffraction processes? Here, one should remind that under WPOA, the contrast in HRTEM images is proportional to the projected specimen potential, convoluted with the impulse response of the instrument.51 Note, however, that diffraction peaks due to moiré structures (or double diffraction) cannot be observed when acquired at typical TEM electron energies (60300 keV) under WPOA, provided that the projected specimen potential consists only of the potentials of individual atoms in the respective lay-ers55 (see 'Intensity of moiré peaks in diffraction pattern' sec-tion in the Supporting Information). We hence assume that the observed extra diffraction spots are originated from a well-defined moiré potential newly created by interlayer coupling and interactions. The moiré potential presumably has complicated spatial distributions, periodicity and ampli-tude modulation depending on the twist angle. It is most likely that specific sets of diffraction sports from different crystalline faces of the twisted hBN layers are required to  Figure 7. Atomic-scale moiré potentials of the twisted hBN layers. (a-d) (top panel) False-color IFFT images obtained by selecting the spots other than the Bragg and center spots shown in Fig. 5ad, yielding a triangular lattice-like pattern for the samples with twist angles of 24.0° (a), 21.6° (b), 19.8° (c), and 13.6° (d). Scale bar, 1 nm. (lower panels) Line profiles as a function of distance along the three yellow lines (1), (2) and (3) in the corresponding top panel.  6  represent such a spatially complicated moiré potential. When an electron beam is diffracted by the complicated moiré potential, the resulting diffraction pattern will not nec-essarily have 6-fold rotational symmetry except for the com-mensurate MSLs but should be represented by the combina-tion of specific set of ke vectors with particular frequencies and amplitudes. This accounts for the reason why the extra diffraction spots with particular frequencies derived from different layers are observed in the FFT images in Fig. 5. Construction of the Moiré Potential. If the assump-tion mentioned in the previous subsection is valid, the in-verse FFT (IFFT) of these extra spots will yield information on the spatial modulation of the moiré potential. It is hence interesting to reconstruct IFFT images from the FFT images by filtering in the frequency domain using the inclusive mask for the extra spots shown in Fig. 5, but excluding all the well-defined Bragg and center spots (for details, see Fig. S7). Figure 7a-d shows the thus obtained IFFT images using the FFT images shown in Fig. 5a-d. We found that regardless of the twisted angles, the IFFT images exhibit a triangular lattice-like pattern. Note also that a similar pattern can be obtained even if the contrast of the respective IFFT images is inverted (Fig. S8). In general, the atom contrast in HRTEM images can be inverted depending on various fac-tors, including thickness, orientation and scattering of the sample, objective lens focus and aberrations, electron beam coherence and convergence angle on the sample.51,53 Hence, the HRTEM imaging mode does not necessarily yields a structure image. If the IFFT images shown in Fig. 7 repre-sent the atom contrast of each layer, the image would look very different by inverting the contrast. As noted above, however, the apparent feature of the IFFT images are hardly affected by the contrast inversion. We thus suggest that the  Figure 8. DFT cluster calculations of electrostatic potential (ESP). (a-e) Optimized structure (left panel), ESP map (middle panel), and magnified ESP map (right panel) of the (B111N111N42)2 clusters with twist angles of (a) 0°, (b) 24.0º, (c) 21.6º, (d) 19.8°, and (e) 13.6°calculated at the B97XD/6-31G(d) level. In the left (right) panels, triangular lattice-like patterns created by the interlayer atomic overlap, as shown by the magnified views in the left panels, are indicated by red (white) triangles. The unit of the scale bar is e a.u.1. (f) Comparison between the ESP map (upper panels) and the corresponding IFFT image shown in Fig. 7 (lower panels) with the same twist angle. These images are created by cropping and rotating the original images for the ease of comparison.  7  resulting IFFT images are not associated with the atom con-trast of each layer but will represent interlayer 2D potential functions with comparable positive and negative amplitudes.  As for one nearly commensurate case at = 21.6° shown in Fig. 7b, a triangular lattice-like pattern is visible over the entire region observed. Note, however, that the detailed local structure varies as we go from top to bottom of the image. This is probably due to the fact the effect of incommensura-bility inherent in the original twisted layer is not completely negligible. Irrespective of the apparently inhomogeneous lattice structure, the corresponding line profiles show dis-tinct repeating patterns with positive and negative correla-tions along the three directions of the triangular lattice, which extends over the length scale of several to tens of na-nometers (see Fig.S9a for a wide-area IFFT image of = 21.6°). Thus, the periodicity of the MSLs is almost preserved in the IFFT image. Although such a lattice-like pattern is not clearly obvious in the other nearly commensurate case at =13.6° (see Fig. 7d), the line profiles indeed show a peri-odic pattern along the three sides of a triangle. The periodic-ities of the lattice for the IFFT images at = 21.6° and 13.6° are estimated to be ~0.6 nm and ~1.0 nm, respectively. These values are comparable to those deduced from Eq. (1) (D = 0.66 nm and 1.08 nm for = 21.6° and 13.6°, respectively). This allows us to confirm that the periodic patterns recog-nized in Fig. 7b,d represent the moiré periodicity or the moiré potential. For the samples with =24.0° (Fig. 7a) and 19.8° (Fig. 7c), the translational symmetry is restricted within the range of a few nm (see Fig. S9b for a wide-area IFFT image of = 24.0°). Accordingly, the corresponding line profile strongly varies from site to site, even showing a 1D stripe-like feature in some regions in Fig. 7c. The ob-served quasi-periodicity can be viewed as an incommensu-rate periodicity derived from the incommensurate layered structures of these MSLs. Calculations of electrostatic potential. The pre-sent IFFT images presumably represent the interlayer moiré potential in our incommensurate and nearly commensurate hBN layers. However, additional theoretical support would be necessary to confirm this argument. In theoretical studies of 2D materials, two models are usually employed, i.e., a pe-riodic model and a finite model. The periodic model is based on Blochʼs theorem and utilizes periodic boundary condi-tions to replicate the elementary cell.14 This approach is in principle applicable only to commensurate systems where the momentum is well defined. In the finite model approach, one performs, for example, quantum chemical calculations using a large enough finite cluster models, whose external dangling bonds are terminated usually by hydrogen at-oms.56,57 Although the results of the finite models may suffer from edge effects,57 this approach can be applied both to commensurate and incommensurate systems and is useful for calculating electrostatic potential (ESP) in real space. Thus, in this work, we carried out a series of quantum chem-ical density functional theory (DFT) calculations using suf-ficiently large clusters of atoms modeling the hBN bilayers with different twist angles as well as the one with the un-twisted AAstacking configuration (θ = 0°). Although the DFT calculations were performed on several model clusters with different sizes using various density functionals under the rigid lattice approximation, we only show the results on the (B111N111H42)2 cluster at the B97XD/6-31G(d) level, as explained in 'Computational details' section in the Support-ing Information. Figure 8 shows the resulting optimized structures of the H-terminated clusters with different twist angles, along with the 2D plot of the ESP in a plane bisecting the two hBN lay-ers of the respective clusters. One sees from Fig. 8a that as for the cluster with θ = 0°, the periodicity of the interlayer ESP matches that of the atomic configuration (for the cross-sectional map perpendicular to the atomic planes, see Fig. S10). The periodic pattern is created by the negative ESP on the N atoms and the positive ESP on the B atoms located exactly above (or below) the N atoms. As the twist angle in-creases from zero, the anti-aligned AA stacking (B on N or N on B) is altered to form various interlayer atomic align-ments, such as N on N and B on B configurations, especially for the nearly commensurate clusters (Fig. 8c,e). From the 2D plot of the ESP shown in Fig. 8c,e, one sees that such N/N (B/B) overlaps yield negative (positive) ESP. Conse-quently, the point where a fully eclipsed arrangement of N/N (B/B) shows the highest local negative (positive) ESP. This leads to the formation of triangular lattice-like patterns with pitches of 0.66 nm and 1.09 nm for θ = 21.6° and 13.6°, re-spectively, which are almost equal to those deduced from Eq. (1) (D = 0.66 nm and 1.08 nm for = 21.6° and 13.6°, re-spectively). Note also that even in the 2D ESP map of in-commensurate clusters with θ = 24.0° (Fig. 8b) and 19.8° (Fig. 8d), a similar periodic pattern with a pitch of 0.66 nm is created due to periodicity of N nearly on N (or B nearly on B) configurations. As for incommensurate MSLs, the pe-riodicity is not actually well defined in real systems. Hence, a similar periodicity obtained for the clusters with θ = 21.6°, 24.0° and 19.8° is valid only in the limitation of the present cluster calculations. However, the present results imply that the corresponding moiré potential is determined by the local interlayer atomic overlap for both the commensurate and in-commensurate hBN nanolayers. Further worth mentioning is that as shown in Fig. 8f, the triangular lattice-like patterns of the 2D ESP maps are very similar to those found in the IFFT images of the samples with the corresponding twist angles in terms both of the pitch and the pattern. A slight discrep-ancy between the observed and calculated pitch values may result from the rigid lattice approximation employed in the DFT calculations. This good correspondence confirms our assumption that a series of the IFFT images shown in Fig. 7 represent the interlayer moiré potential created by the local interlayer atomic overlap in the respective hBN bilayers. Thus far, the discussion was given based on the twist an-gle 𝜃 defined earlier. What happens when we employ an al-ternative definition of the twist angle, i.e., 𝜃̅ = 60° − 𝜃? To answer the question, we calculated the 2D ESP map of one of the nearly commensurate cluster with 𝜃̅ = 60° − 21.6° =8  38.4°. In the 38.4° cluster, the aligned N on N (or B on B) stacking points realized in the 21.6° cluster are replaced by the anti-aligned N on B (or B on N) stacking points, yielding locally nearly zero potential regions. Even in the 38.4° clus-ter, however, there exist a number of N nearly on N (or B nearly on B) configurations, whose periodicity is similar to that of the N on N (or B on B) stacking points in the 21.6° cluster. These N nearly on N (B nearly on B) configurations will yield locally negative (positive) EPS, although their am-plitudes are slightly smaller than those for the N exactly on N (B exactly on B) configurations. Consequently, as shown  in Fig. 9, the 2D ESP map at 38.4° show a triangular lattice-like pattern with nearly the same periodicity as that at 21.6°, implying that both the ESP maps at  21.6° and 38.4° fit the IFFT image shown in Fig. 7b. Although the moiré potentials at θ and 𝜃̅ are indistinguishable, their detailed atomic con-figurations are different. It is hence possible that the elec-tronic structures at θ and 𝜃̅ are not equivalent, as pointed out in commensurate twisted bilayer graphene in the gap feature (gapped or ungapped) and the related low-energy electronic spectra.17 However, the possible difference in the electronic structure of twisted hBN layers at θ and 𝜃̅  cannot be identified by this method. Thus, at present we cannot defi-nitely decide which angle, 𝜃 or 𝜃̅, is more consistent with the reconstructed IFFT image Finally, we discuss how the interlayer ESP varies as the interlayer distance dint increases from its optimal value dint(opt). Figure 10 shows a variation in the ESP with inter-layer distance obtained for the cluster with θ = 24.0°, which is characterized by dint(opt) = 0.329 nm. One sees from Fig. 10 that the positive ESP almost disappears for dint=0.375 nm, and the interlayer ESP eventually becomes almost structure-less for din=0.475 nm. This probably accounts for the reason why in previous HRTEM studies on MSLs, the extra FFT diffraction spots, such as those observed in Fig. 5, have not been reported to occur.3237 In the previous HRTEM meas-urements, it is likely that various adsorbates are trapped be-tween layers during their assembly and sample handling, which prevents the observation of interlayer moiré potential. Thus, as mentioned repeatedly before, the preparation of ad- Figure 10. Variation in the ESP with interlayer distance. (a) Enlarged representation of the BN bilayer model with a twist angle of 24.0º. (b-h) ESP for the interlayer distance dint of (b) 0.329 nm (optimized value) , (c) 0.350 nm, (d) 0.375 nm, (e) 0.400 nm, (f) 0.425 nm, (g) 0.450 nm, and (h) 0.475 nm. ESP is shown for a plane bisecting the two hBN layers of the respective clusters, and the unit of the scale bar is e a.u.1. A solid (dotted) triangle represents a portion of the triangular-like lattice consisting of  N nearly on N (B nearly on B) configurations.   Figure 9. ESP map of the (B111N111N42)2 clusters with twist angles of (a) 21.6º and (b) 38.4° calculated at the B97XD/6-31G(d) level. The unit of the scale bar is e a.u.1. 9  sorbate-free samples is a prerequisite to observe moiré po-tential as FFT diffraction spots. The chemical exfoliation procedure developed in this work is one useful method for obtaining clean hBN nanolayers and hence observing their moiré potential. We are now applying the present method to graphite and transition metal dichalcogenides. Currently, studies are in progress to confirm its applicability to other vdW systems.  ■ CONCLUSIONS We fabricated strain- and contaminant-free hBN layers with different twist angles using the intercalation-based ex-foliation technique and observed their HRTEM images. The FFT of the HRTEM images allowed us to recognize moiré diffraction spots, and their IFFT computations provided in-formation about the interlayer moiré potential, as corrobo-rated by the DFT calculations. Also, the DFT calculations highlight the importance of the local interlayer atomic over-lap in establishing the moiré potentials both in the commen-surate and incommensurate cases. The present results not only pave the way for observing the moiré potentials in twisted hBN layers, but they also expand the application of electron microscopy to visualizing the moiré potential in a variety of MSLs especially with large twist angles.  ASSOCIATED CONTENT  Supporting Information.   The supporting Information is available free of charge at http://pubs.acs.org. Experimental and computational details, scheme of  chemical exfoliation, additional HRTEM and IFFT images, PPA analysis, real and reciprocal representa-tions of  hBN bilayers, cross sectional ESP map, and addi-tional DFT results.  AUTHOR INFORMATION Corresponding Author *Takashi Uchino - Department of Chemistry, Kobe Univer-sity, Nada, Kobe 657-8501, Japan; orchid.org/0000-0002-4899-3078; Email: uchino@kobe-u.ac.jp. Author Contributions The manuscript was written through contributions of all au-thors. / All authors have given approval to the final version of the manuscript. /  ACKNOWLEDGMENT  This work was conducted in National Institute for Materials Science (NIMS) and Institute for Molecular Science, sup-ported by “Advanced Research Infrastructure for Materials and Nanotechnology in Japan (ARIM)” of the Ministry of Education, Culture, Sports, Science and Technology (Pro-posal Numbers JPMXP1223NM0091 and JPMXP1224 MS0004). This work was also supported by NIMS Joint Re-search Hub Program. The computations were performed in the Research Center for Computational Science, Okazaki, Japan (24-IMS-C194). We also acknowledge Research Fa-cility Center for Science and Technology, Kobe University, for providing access to the AFM facility.  REFERENCES (1) Park, J. M.; Cao, Y.; Xia, L.-Q.; Sun, S.; Watanabe, K.; Taniguchi, T.; Jarillo-Herrero, P. Robust Superconductivity in Magic-Angle Multilayer Graphene Family. Nat. Mater. 2022, 21 (7), 877–883. (2) Andrei, E. Y.; Efetov, D. K.; Jarillo-Herrero, P.; MacDonald, A. H.; Mak, K. F.; Senthil, T.; Tutuc, E.; Yazdani A.; Young, A. F. The Marvels of Moiré Materials. Nat. Rev. Mater. 2021, 6 (3), 201–206. (3)  Morales-Durán, N.; Shi, J.; MacDonald, A. H. Fractionalized Electrons in Moiré Materials. Nat. Rev. Phys. 2024, 6 (4), 349–351. (4) He, F.; Zhou, Y.; Ye, Z.; Cho, S.-H.; Jeong, J.; Meng, X.; Wang, Y. Moiré Patterns in 2D Materials: A Review. ACS Nano 2021, 15 (4), 5944-5958. (5) Mak K. F.; Shan, J. Semiconductor Moiré Materials. Nat. Nanotechnol. 2022, 17 (7), 686–695. (6)  Nuckolls K. P.; Yazdani, A. A microscopic Perspective on Moiré Materials. Nat. Rev. Mater. 2024, 9 (5), 460–480. (7) Cao, Y.; Fatemi, V.; Demir, A.; Fang, S.; Tomarken, S. L.; Luo, J. Y.; Sanchez-Yamagishi, J. D.; Watanabe, K.; Taniguchi, T.; Kaxiras, E. et al. Correlated Insulator Behaviour at Half-Filling in Magic-Angle Graphene Superlattices. Nature 2018, 556, 80–84. (8) Cao, Y.; Fatemi, V.; Fang, S.; Watanabe, K.; Taniguchi, T.; Kaxiras E.; Jarillo-Herrero, P. Unconventional Superconductivity in Magic-Angle Graphene Superlattices. Nature 2018, 556, 43–50. (9)  Saito, Y.; Ge, J.; Watanabe, K.; Taniguchi, T.; Young, A. F. Independent Superconductors and Correlated Insulators in Twisted Bilayer Graphene, Nat. Phys. 2020, 16 (6), 926–930. (10) Park, H.; Cai, J.; Anderson, E.; Zhang, Y.; Zhu, J.; Liu, X.; Wang, C.; Holtzmann, W.; Hu, C.; Liu, Z. et al. Observation of Fractionally Quantized Anomalous Hall Effect. Nature 2023, 622, 74–79. (11) Lu, Z.; Han, T.; Yao, Y.; Reddy, A. P.; Yang, J.; Seo, J.; Watanabe, K.; Taniguchi, T.; Fu, L.; Ju, L. Fractional Quantum Anomalous Hall Effect in Multilayer Graphene. Nature 2024, 626, 759–764. (12)  Dean, C. R.; Wang, L.; Maher, P.; Forsythe, C.; Ghahari, F.; Gao, Y.; Katoch, J.; Ishigami, M.; Moon, P.; Koshino, M. et al. Hofstadter’s Butterfly and the Fractal Quantum Hall Effect in Moiré Superlattices. Nature 2013, 497, 598–602. (13)  Ju, L.; MacDonald, A. H.; Mak, K. F.; Shan J.; Xu, X. The Fractional Quantum Anomalous Hall Effect. Nat. Rev. Mater. 2024, 9 (6), 455–459. (14) Carr, S.; Fang S.; Kaxiras, E. Electronic-Structure Methods for Twisted Moiré Layers. Nat. Rev. Mater. 2020, 5 (7), 748–763. (15) Shallcross, S.; Sharma, S.; Kandelaki, E.; Pankratov, O. A. Electronic Structure of Turbostratic Graphene. Phys. Rev. B 2010, 81 (16), 165105. (16)  Carr, S.; Fang, S.; Zhu, Z.; Kaxiras, E. Exact Continuum Model for Low-Energy Electronic States of Twisted Bilayer Gra-phene. Phys. Rev. Research 2019, 1 (1), 013001. (17) Mele, E. J. Commensuration and Interlayer Coherence in Twisted Bilayer Graphene. Phys. Rev. B 2010, 81 (16), 161405(R). (18) Conte, F.; Ninno, D.; Cantele, G. Electronic Properties and Interlayer Coupling of Twisted MoS2/NbSe2 Heterobilayers. Phys. Rev. B 2019, 99 (15), 155429. 10  (19)  Luican, A.; Li, G.; Reina, A.; Kong, J.; Nair, R. R.; Novose-lov, K. S.; Geim, A. K.; Andrei, E. Y. Single-Layer Behavior and Its Breakdown in Twisted Graphene Layers. Phys. Rev. Lett. 2011, 106 (12), 126802. (20)  Yao, W.; Wang, E.; Bao, C.; Zhang, Y.; Zhang, K.; Bao, K.; Chan, C. K.; Chen, C.; Avila, J.; Asensio, M. C. et al. Quasicrys-talline 30° Twisted Bilayer Graphene as an Incommensurate Super-lattice with Strong Interlayer Coupling. Proc. Natl. Acad, Sci. USA 2018, 115 (27), 6928-6933. (21) Ahn, S. J.; Moon, P.; Kim, T.-H.; Kim, H.-W.; Shin, H.-C.; Kim, E. H.; Cha, H. W.; Kahng, S.-J.; Kim, P.; Koshino, M. et al. Dirac Electrons in a Dodecagonal Graphene Quasicrystal. Science 2018, 361 (6404), 782-786. (22) Pezzini, S.; Mišeikis, V.; Piccinini, G.; Forti, S.; Pace, S.; Engelke, R.; Rossella, F.; Watanabe, K.; Taniguchi, T.; Kim, P. et al. 30°-Twisted Bilayer Graphene Quasicrystals from Chemical Vapor Deposition, Nano Lett. 2020, 20 (5), 3313-3319. (23)  Li, Y.; Zhang, F.; Ha, V.-A.; Lin, Y.-C.; Dong, C.; Gao, Q.; Liu, Z.; Liu, X.; Ryu, S. H.; Kim, H. et al. Tuning Commensurabil-ity in Twisted van der Waals Bilayers. Nature 2024, 625, 494–499. (24) Sponza, L.; Vu, V. B.; Richaud, E. S.; Amara, H.; Latil, S. Emergence of Flat Bands in the Quasicrystal Limit of Boron Ni-tride Twisted Bilayers. Phys. Rev. B 2024, 109 (16), L161403. (25)  Agarwal, H.; Nowakowski, K.; Forrer, A.; Principi, A.; Ber-tini, R.; Batlle-Porro, S.; Reserbat-Plantey, A.; Prasad, P.; Vistoli, L.; Watanabe, K. et al. Ultra-Broadband Photoconductivity in Twisted Graphene Heterostructures with Large Responsivity. Nat. Photon. 2023, 17 (9), 1047–1053. (26) Bucko, J.; Herman, F. Large Twisting Angles in Bilayer Gra-phene Moiré Quantum Dot Structures. Phys. Rev. B 2021, 103 (7), 075116. (27) Deng B.; Xia, F. Large-Angle Twist Effect. Nat. Photon. 2023, 17, 1021–1022. (28)  Woods, C. R.; Britnell, L.; Eckmann, A.; Ma, R. S.; Lu, J. C.; Guo, H. M.; Lin, X.; Yu, G. L.; Cao, Y.; Gorbachev, R. V. et al. Commensurate–Incommensurate Transition in Graphene on Hex-agonal Boron Nitride. Nat. Phys. 2014, 10 (4), 451–456. (29)  Deng, B.; Wang, B.; Li, N.; Li, R.; Wang, Y.; Tang, J.; Fu, Q.; Tian, Z.; Gao, P.; Xue, J. et al. Interlayer Decoupling in 30° Twisted Bilayer Graphene Quasicrystal. ACS Nano 2020, 14 (2), 1656-1664. (30)  Yin, J.; Li, J.; Hang, Y.; Yu, J.; Tai,G.; Li, X.; Zhang, Z.; Guo, W. Boron Nitride Nanostructures: Fabrication, Functionalization and Applications. Small 2016, 12 (22), 2942–2968. (31)  Kim, D. S.; Dominguez, R. C.; Mayorga-Luna, R.; Ye, D.; Embley, J.; Tan, T.; Ni, Y.; Liu, Z.; Ford, M.; Gao, F. Y. et al. Electrostatic Moiré Potential from Twisted Hexagonal Boron Ni-tride Layers. Nat. Mater. 2024, 23 (8), 65–70. (32) Robertson, A. W.; Bachmatiu, A.; Wu, Y. A.; Schäffel, F.; Büchner, B.; Rümmeli, M. H.; Warner, J. H. Structural Distortions in Few-Layer Graphene. ACS Nano 2011, 5 (12), 9984–9991. (33) Zan , R.; Bangert, U.; Ramasse, Q.; Novoselov, K. S. Imag-ing of Bernal Stacked and Misoriented Graphene and Boron Ni-tride: Experiment and Simulation. J. Microscopy 2011, 244 (2), 152-158. (34)  Yuan, S.; Linas, S.; Journet, C.; Steyer, P.; Garnier, V.; Bonnefont, G.; Brioude, A.; Toury, B. Pure & Crystallized 2D Bo-ron Nitride Sheets Synthesized via a Novel Process Coupling both PDCs and SPS methods. Sci. Rep. 2016, 6, 20388. (35)  Warner, J. H.; Rümmeli, M. H.; Gemming, T.; Büchner, B.; Andrew, G.; Briggs, D. Direct Imaging of Rotational Stacking Faults in Few Layer Graphene. Nano Lett. 2009, 9 (1), 102-106. (36)  Bachmatiuk, A.; Zhao, J.; Gorantla, S. M.; Martinez, I. G. G.; Wiedermann, J.; Lee, C.; Eckert, J.; Rummeli, M. H. Low Voltage Transmission Electron Microscopy of Graphene. Small 2015, 11(5), 515-542. (37)  Yamazaki, K.; Maehara, Y.; Gohara, K. Characterization of TEM Moiré Patterns Originating from Two Monolayer Graphenes Grown on the Front and Back Sides of a Copper Substrate by CVD Method. J. Phys. Soc. Jpn. 2018, 87 (6), 061011. (38) Haigh, S. J.; Gholinia, A.; Jalil, R.; Romani, S.; Britnell, L.; Elias, D. C.; Novoselov, K. S.; Ponomarenko, L. A.; Geim A. K.; Gorbachev, R. Cross-Sectional Imaging of Individual Layers and Buried Interfaces of Graphene-Based Heterostructures and Super-lattices. Nat. Mater. 2012, 11 (7), 764–767. (39) Gasparutti, I.; Song, S. H.; Neumann, M.; Wei, X.; Watanabe, K.; Taniguchi, T.; Lee, Y. H. How Clean Is Clean? Recipes for van der Waals Heterostructure Cleanliness Assessment. ACS Appl. Ma-ter. Interfaces 2020, 12 (6), 7701–7709. (40) Rosenberger, M. R.; Chuang, H.-J.; McCreary, K. M.; Han-bicki, A. T.; Sivaram, S. V.; Jonker, B. T. Nano-“Squeegee” for the Creation of Clean 2D Material Interfaces. ACS Appl. Mater. Inter-faces 2018, 10 (12), 10379–10387. (41)  Cai L.; Yu, G. Fabrication Strategies of Twisted Bilayer Gra-phenes and Their Unique Properties. Adv. Mater. 2021, 33 (13), 2004974. (42)  Liao, M.; Wei, Z.; Du, L.; Wang, Q.; Tang, J.; Yu, H.; Wu, F.; Zhao, J.; Xu, X.; Han, B. et al. Precise Control of the Interlayer Twist Angle in Large Scale MoS2 Homostructures. Nat. Commun. 2020, 11, 2153.  (43) Yang, R.; Fan, Y.; Mei, L.; Shin, H. S.; Voiry, D.; Lu, Q.; Li, J.; Zeng, Z. Synthesis of Atomically Thin Sheets by the Intercala-tion-Based Exfoliation of Layered Materials. Nat. Synth. 2023, 2 (2), 101–118. (44) Kovtyukhova, N. I.; Wang, Y.; Lv, R.; Terrones, M.; Crespi, V. H.; Mallouk, T. E. Reversible Intercalation of Hexagonal Boron Nitride with Brønsted Acids. J. Am. Chem. Soc. 2013, 135 (22), 8372. (45) Tsujimura, T.; Uchino, T. Oriented Crystallization of Ammo-nium Sulfate from Hexagonal Boron Nitride/Sulfuric Acid Interca-lation Compounds. ACS Omega 2021, 6 (9) , 6482-6487. (46) Galindo, P. L.; Kret, S.; Sanchez, A.M.; Laval, J.-Y.; Yáñez, A.; Pizarro, J.; Guerrero, E.; Ben, T.; Molina, S. I. The Peak Pairs Algorithm for Strain Mapping from HRTEM Images. Ultramicros-copy 2007, 107 (12), 1186-1193. (47) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Peters-son, G. A.; Nakatsuji, H. et al. Gaussian 16, Revision C.02; Gauss-ian, Inc.: Wallingford, CT, 2019. (48) Krivanek, O. L.; Chisholm, M. F.; Nicolosi, V.; Pennycook, T. J.; Corbin, G. J.; Dellby, N.; Murfitt, M. F.; Own, C. S.; Szilagyi, Z. S.; Oxley, M. P. et al. Atom-by-Atom Structural and Chemical Analysis by Annular Dark-Field Electron Microscopy. Nature 2010, 464, 571–574. (49) Odlyzko, M. L.;  Mkhoyan, K. A. Identifying Hexagonal Bo-ron Nitride Monolayers by Transmission Electron Microscopy. Mi-crosc. Microanal. 2012, 18, 558-567. (50)  Seto, Y.; Ohtsuka, M. ReciPro: Free and Open-Source Mul-tipurpose Crystallographic Software Integrating a Crystal Model Database and Viewer, Diffraction and Microscopy Simulators, and Diffraction Data Analysis Tools. J. Appl. Cryst. 2022, 55(2), 397410. (51) Spence, J. C. H. High-Resolution Electron Microscopy, 4th ed.; Oxford University Press: New York, NY, 2013. (52) Gass, M. H.; Bangert, U.; Bleloch, A. L.; Wang, P.; Nair, R. R.; Geim, A. K. Free-Standing Graphene at Atomic Resolution. Nat. Nanotechnol. 2008, 3 (9), 676–681. 11  (53) Williams, D. B; Carter, C. B. Transmission Electron Micros-copy: A Textbook for Materials Science, 2nd ed.; Springer: New York, NY, 2009.  (54)  Lee, H. Y.; Al Ezzi, M. M. ; Raghuvanshi, N. ; Chung, J. Y.; Watanabe, K.; Taniguchi, T.; Garaj, S.; Adam, S.; Gradečak, S. Tunable Optical Properties of Thin Films Controlled by the Inter-face Twist Angle, Nano Lett. 2021, 21(7), 2832–2839. (55)  Latychevskaia, T.; Escher, C.; Fink, H.-W. Moiré Structures in Twisted Bilayer Graphene Studied by Transmission Electron Microscopy. Ultramicroscopy 2019, 197, 46-52. (56) Mendelson, N.; Chugh, D.; Reimers, J. R.; Cheng, T. S.; Gottscholl, A.; Long, H.; Mellor, C. J.; Zettl, A.; Dyakonov, V.; Beton, P. H. et al. Identifying Carbon as the Source of Visible Sin-gle-Photon Emission from Hexagonal Boron Nitride. Nat. Mater. 2021, 20 (11), 321–328. (57) Karamanis, P.; Otero, N.; Pouchan, C. Electric Property Var-iations in Nanosized Hexagonal Boron Nitride/Graphene Hybrids. J. Phys. Chem. C 2015, 119 (21), 11872–11885.      Insert Table of Contents artwork here TOC Graphic                             1  Supporting Information: Atomic-Scale Observation of Moiré potential in Twisted Hexagonal  Boron Nitride Layers by Electron Microscopy  Rina Mishima,†Takuro Nagai,§ Hiroyo Segawa,ǁ Masahiro Ehara,‡ and Takashi Uchino†*  †Department of Chemistry, Graduate School of Science, Kobe University, Nada, Kobe 657-8501, Japan §Research Network and Facility Services Division, National Institute for Materials Science, Tsukuba, Ibaraki 305-0044, Japan ǁResearch Center for Electronic and Optical Materials, National Institute for Materials Science, Tsukuba, Ibaraki 305-0044, Japan ‡Research Center for Computational Science, Institute for Molecular Science, Okazaki, Aichi, 444-8585, Japan  *Corresponding author: uchino@kobe-u.ac.jp        Table of Contents  S1. Sample preparation..............................................................................................................  1 S2. TEM observations................................................................................................................  2  S3. Intensity of moiré peaks in diffraction pattern.................................................................  3 S4. Computational details..........................................................................................................  4 S5. Supporting figures................................................................................................................  7 S6. Supporting references.......................................................................................................... 20         2  S1. Sample preparation Boron Nitride was purchased from Kojundo Chemical Lab. Co., Ltd. (purity ∼99%, particle size ∼10 μm) without any additional treatment. We confirmed that the sample powders show Bragg X-ray diffraction peaks attributed only to the hBN phase. Sulfuric acid (95 %) was obtained from Wako Pure Chemical Industries, Ltd. and was used as purchased. We prepared h-BN/H2SO4 intercalation compounds according to the scheme shown schematically in Fig. S1. A 30 mg amount of h-BN was added to 0.5 mL of sulfuric acid in a 50 mL Teflon-lined stainless steel autoclave. The autoclave was heated at 200 °C for 24 h, and then naturally cooled to room temperature, producing a slightly brownish slurry. The slurry was transferred into a test tube and was neutralized with an aqueous solution of NaHCO3, resulting in an opaque solution. During neutralization, sodium sulfate salt (Na2SO4) is expected to be formed in between the hBN layers, leading to spontaneous exfoliation of hBN layers The solution was centrifuged at 15000 rpm for 45 min to remove insoluble materials, and then the supernatant was transferred to a new test tube. The thus obtained supernatant is an aqueous solution containing Na2SO4 and exfoliated hBN nanosheets. A solvent extraction method was utilized to isolate the exfoliated hBN nanosheets from the supernatant solution. That is, the solution was mixed with an equal volume of organic solvent in a separatory funnel, leading to migration of hBN nanosheets into the organic solvent. We found that 1-pentanol is the most suitable solvent for the present purpose. Although the extracted 1-pentanol solvent was transparent, a colloidal dispersion of the exfoliated hBN nanosheets was confirmed by the Tyndall effect (Fig. S1).  For the TEM observations, the extracts of a volume of ~20 mL were concentrated to a volume of about 0.7 mL, then ethanol of 20 mL was added. This solvent exchange is useful to minimize the amount of hydrocarbon adsorbates on the sample. Then, 5 μL of the suspension was collected with a pipette and drop casted onto a TEM microgrid. The TEM grid was dried in vacuum using a dry pumping station for ~2 h to remove solvent.   S2. TEM observations Atomic-resolution HRTEM and STEM observations were carried out on a JEOL JEM-ARM300F instrument equipped with a cold field emission electron gun and a spherical aberration (Cs) corrector at an acceleration voltage of 80 kV, under 1 × 10−5 Pa in the specimen column. Typically, we adjusted the Cs value to 1 μm. All HRTEM images were recorded with an exposure time of 4 s on a Complementary Metal-Oxide-Semiconductor (CMOS) camera (Gatan OneView, 4096 × 4096 pixels). In ADF STEM imaging, the probe size was 0.07 nm, the convergence angle was 32 mrad, and the inner and outer collection angles were 45 and 180 mrad, respectively. Although the operated voltage appears to be higher than the knock-on radiation damage threshold reported for h-BN (78 kV),58 we did not observe any radiation damage during the TEM observations. Image processing including the fast Fourier transformation (FFT) and inverse FFT (IFFT) was carried out by Gatan Digital Micrograph. We analysed strain of the exfoliated hBN layers using the Peak Pairs Analysis (PPA) software package (HREM research) for Gatan Digital Micrograph.46 For the analysis, we employed low-pass filtered HRTEM images of the non-twisted hBN layers. The peak positions were determined on the filtered image, and the relative 3  displacement fields (ux, uy) of the measured lattice with respect to a reference basis vector were calculated. The components of strain tensor 𝜀𝑥𝑥 , 𝜀𝑥𝑦 , 𝜀𝑦𝑦 , mean dilatation ∆𝑥𝑦 , rotation 𝜔𝑥𝑦  (in radians) are defined as follows:  𝜀𝑥𝑥 =𝜕𝑢𝑥𝜕𝑥, 𝜀𝑥𝑦 =12(𝜕𝑢𝑥𝜕𝑦+𝜕𝑢𝑦𝜕𝑥), 𝜀𝑦𝑦 =𝜕𝑢𝑦𝜕𝑦, ∆𝑥𝑦=12(𝜀𝑥𝑥 + 𝜀𝑦𝑦), 𝜔𝑥𝑦 =12(𝜕𝑢𝑦𝜕𝑥−𝜕𝑢𝑥𝜕𝑦).   S3. Intensity of moiré peaks in diffraction pattern We explain how the moiré peaks in the diffraction pattern of twisted bilayers are generated. First, we consider the electron transmission through a non-twisted 2D crystal. The total transmission function 𝑡(𝑥, 𝑦)  of the sample is described by:  𝑡(𝑥, 𝑦) = exp[−𝑖𝛷(𝑥, 𝑦)] = exp[−𝑖𝜎𝑉(𝑥, 𝑦)]    (S1)  where 𝛷(𝑥, 𝑦) is the sample phase shift, 𝜎 is the interaction parameter, 𝑉(𝑥, 𝑦) is the projected potential of the entire sample. The interaction parameter 𝜎 is given by σ =2𝜋𝑚𝑒𝜆ℎ2 , where m is the relativistic mass of the electron, e is the elementary charge,  is the wavelength of the electrons, and h is the Planckʼs constant. Eq. (S1) can further be approximated by the first order Taylor series as  𝑡(𝑥, 𝑦) ≈ 1 − 𝑖𝜎𝑉(𝑥, 𝑦) (S2)  under the weak-phase object approximation (WPOA),51 which is applicable to very thin crystals of light atoms, as in the case of hBN nanolayers.  Then we move on to the case of twisted bilayers. For the time being, we do not take into account the effect of the interlayer moiré potential in the electron transmission process. When neglecting the diffraction effects due to propagation between the two layers, the transmission function of the twisted bilayers can be written as under WPOA:  𝑡(𝑥, 𝑦) = [1 − 𝑖𝜎𝑉𝑡(𝑥, 𝑦)][1 − 𝑖𝜎𝑉𝑏(𝑥, 𝑦)]  = 1 − 𝑖𝜎𝑉𝑡(𝑥, 𝑦) − 𝑖𝜎𝑉𝑏(𝑥, 𝑦) − 𝜎2𝑉𝑡(𝑥, 𝑦)𝑉𝑏(𝑥, 𝑦)  (S3)  where 𝑉𝑡(𝑥, 𝑦) and 𝑉𝑏(𝑥, 𝑦) are the projected potentials of the top and the bottom layers, respectively. The last term describes the overlap of the top- and bottom-layer potentials, which can, in principle, contribute to the formation of moiré peaks in the diffraction pattern. In high-energy (20  300 keV) electron microscopy, however, 4  the related moiré peaks are not observed because of the following reasons.55 From Eq. (S3), one sees that the second and third terms, which describes the intensity of the main peaks, are proportional to 𝜎, whereas the last term to 𝜎2. Note, however, that for typical TEM electron energies (20  300 keV), the value of 𝜎 is relatively small, e.g., 𝜎 = 0.81 (keV ∙ Å)−1 at 100 keV. Consequently, the last term proportional to 𝜎2 becomes negligibly small, as compared to the second and third terms, resulting in the absence of moiré peaks in conventional TEM mode. When the electron energy is reduced below ~100 eV, moiré peaks should be visible as 𝜎 becomes relatively high,55 e.g., 𝜎 = 25.61 (keV ∙ Å)−1 for 100 eV. However, atomic-scale imaging is not possible for such low electron energies.   Finally, we take account of the interlayer moiré potential in the electron transmission process. We assume that the interlayer moiré potential is created in between the top and bottom layers as a result of the interlayer coupling and also that this potential explicitly contributes to the phase shift. The resulting transmission function of the twisted bilayers can be described as:  𝑡(𝑥, 𝑦) = [1 − 𝑖𝜎𝑉𝑡(𝑥, 𝑦)][1 − 𝑖𝜎𝑉𝑏(𝑥, 𝑦)][1 − 𝑖𝜎𝑉𝑖𝑛𝑡(𝑥, 𝑦)] = 1 − 𝑖𝜎𝑉𝑡(𝑥, 𝑦) − 𝑖𝜎𝑉𝑏(𝑥, 𝑦) − 𝑖𝜎𝑉𝑖𝑛𝑡(𝑥, 𝑦) …  (S4)  where 𝑉𝑖𝑛𝑡(𝑥, 𝑦) describes the spatial modulation of the moiré potential in the interlayer region. In this case, the last term in Eq. (S4) is proportional to 𝜎, and the intensity of the resulting peaks due to the moiré potential is expected to be comparable to that of the main peaks even for typical TEM electron energies. This argument provides a reasonable explanation for the occurrence of the extra diffraction spots shown in Fig. 5. Hence, we suggest that the absence of the moiré peaks in previous TEM studies on MSLs3035 should not be ascribed to the high TEM electron energies, but rather to the absence of interlayer moiré potential due to the strains and contaminants remained in the samples investigated. Considering that the effect of interlayer decoupling becomes more significant with increasing twist angle,59,60 we would assert that cleanliness of the sample is a matter of great importance in exploring the moiré potential especially for the samples with larger twist angles (see also Fig. 10).   S4. Computational details All quantum chemical calculations were performed by density functional theory (DFT) methods implemented in the Gaussian 16 suite of programs47 using the 6-31G(d) basis set. The suitable selection of functional, and particularly the treatment of long-range exchange interactions in the functional, is crucial for correctly describing vdW interaction.61 To provide comparison, two functionals were used: one is the B3LYP with the dispersion corrections for the non-bonding vdW interactions (B3LYP-D3),62,63 and other is a long-range corrected hybrid density functional (B97XD)64 with a parameterized classical dispersion term. In this work, hydrogen terminated clusters were used to model the local structure of twisted hBN bilayers with four different twist angles (𝜃 = 13.6°, 19.8°, 21.6° and 24.0°). For comparison, we also employed a cluster with zero twist angle, modeling the local structure of pristine hBN bilayers with the AA high-symmetry stacking. 5  To bear out the electronic structure characteristics of the systems considered, we carefully checked how the size and edge shape of the clusters affect or do not affect the calculated results. For this purpose, we employed three types of clusters. Two of them are (B75N75H30)2 and (B147N147H42)2 bearing zigzag edges, and the other one is (B111N111H42)2 featuring armchair edges (Fig. S11). Considering that the twist angles examined in this study are larger than 10° and circumventing possible edge effects of the clusters, we employed the rigid lattice assumption in the structural optimization process; that is, we assumed that each BN layer is completely flat and consists of ideal BN hexagons, in which all the BN distances are identical and all the BNB and NBN bond angles are 60°. As a result of the above structural constraints, the optimized structural parameters are the BN distance (dBN), the interlayer distance between the two BN layers (dint), and the terminal B-H and N-H bond distances. Since our main concern is how the optimized structure, especially the B-N bond distance dBN and the interlayer distance dint, are influenced by the size of the clusters and/or the type of edges (zigzag or armchair), we obtained the optimized values of dBN and dint for the above three types of clusters with twist angles of  = 0° and 24.0° at the B3LYP-D3/6-31G(d) level. Figure S11 summarizes the results of geometry optimization for the above three model clusters, showing that the optimized structural parameters depend hardly on the type of the clusters. We also found that irrespective of the types of the clusters, the optimized values of dBN and dint obtained for the cluster with  = 0° lie in the range 0.14451446 nm and 0.33320.3339 nm, respectively, which are in good agreement with the experimental values of bulk hBN (dBN=0.145 nm, dint= 0.333 nm),65 The good agreement between the calculated and observed bulk values of dint is likely to be fortuitous because the reported dint value for hBN multilayers is slightly shorter (dint= 0.325 nm)66 than the bulk one (see also Fig. 4b). When the value of  increases to 24.0°, the optimized values of dBN remain constant, whereas those of dint show a slight increase for all the clusters employed. This implies that an increase in leads to a decrease in the interlayer interaction For the calculation of the electrostatic potential (ESP), we employ the (B111N111H42)2 cluster with armchair edges because the (B147N147H42)2 cluster is too large to perform such calculations within a reasonable computational cost. The ESP at position r (ϕESP(𝒓)) is given as a sum of contributions from the nuclei and the electronic wave function Ψ given by,  ϕESP(𝒓) = ∑𝑍A|𝒓 − 𝑹A|nucleiA− ∫|Ψ(𝒓′)|2|𝒓 − 𝒓′|𝑑𝒓′,   (S5)  where A represents a nucleus, and 𝑍A and 𝑹A are the charge and position of the nucleus A, respectively. The ESP is generally insensitive to the level of sophistication, e.g., the size of the basis set and the level of electron correlation, since it depends directly on the electron density (ρ = |Ψ|2 ).67 2D ESP was rendered by Visual Molecular Dynamics (VMD)68 and/or Gaussview69 software packages. The effect of density functional on the optimized structure and ESP is also of interest. Figure S12 shows the optimized values of dBN and dint obtained for the (B111N111H42)2 cluster with  = 0°, 13.6°, 19.8°, 21,6° and 24.0° calculated at the B3LYP-D3/6-31G(d) and ωB97XD/6-31G(d) levels. The optimized values of 6  dBN and dint at the ωB97XD/6-31G(d) level yield smaller values than those at the B3LYP-D3/6-31G(d) level irrespective of the twist angle; that is, the value of dint (dBN) at the B97XD/6-31G(d) level is ~0.01 nm (~0.003 nm) smaller than that at the B3LYP-D3/6-31G(d) level. Consequently, as for the (B111N111H42)2 cluster with  = 0°, the ωB97XD/6-31G(d) level yields the dint value of 0.3241 nm, which is quite comparable to that of hBN multilayers (dint = 0.325 nm).66 This suggests that the long-range electron-electron exchange interactions are reasonably incorporated in the B97XD functional. The interlayer ESPs calculated at the B3LYP-D3/6-31G(d) and ωB97XD/6-31G(d) levels for the (B111N111H42)2 cluster with  = 21.6° and 13.6° are given in Fig. S13. We did not find any critical differences between the ESPs at the B3LYP-D3/6-31G(d) and ωB97XD/6-31G(d) levels, although these two levels of calculations result in somewhat different values of dint. This implies that the underlying features of the interlayer ESP are insensitive to the DFT functional. Hence, in the main text, we only show the ESP map for the (B111N111H42)2 cluster calculated at the ωB97XD/6-31G(d) level.                             7  S5. Supporting figures    Figure S1. Exfoliation of hBN layers from hBN powders. (a) Formation of hBN/H2SO4 intercalation compounds, exfoliation of hBN layers via neutralization with NaHCO3, and subsequent extraction of exfoliated hBN layers with 1-pentanol. (b) Microscopic scheme of the exfoliation process occurring spontaneously under the neutralization reaction.    8          Figure S2. (a) Schematic illustration of an x-tilt performed on a hBN 0001 plane. (b) Series of simulated ADF-STEM images for the x-tiled (x-tilt angle: 1°) hBN layers with a different thickness ranging from 4 to 14 nm. (c) An example of the observed ADF-STEM image of the exfoliated specimen. 9     Figure S3. HRTEM images of non-twisted hBN layer. (a) Raw HRTEM image, along with a FFT image (inset) and (b) a magnified HRTEM image of the area marked by the white, dashed square in (a). 10        Figure S4. PPA analysis. (a) HRTEM image used to perform the PPA analysis. (b) Rotation (in degree and anti-clockwise positive), (c) Mean dilatation, (d) Strain tensor. The PPA analysis demonstrates that rotation is within the range of ±1 degree, whereas mean dilatation and strain tensor within the range of ±1 %. 11    Figure S5. Real and reciprocal space representations of 2D hexagonal bilayers with different twist angles in the rigid lattice model. (a) hBN monolayer (left panel) and the corresponding diffraction patterns (right panel). G1 and G2 are reciprocal lattice vectors. (bd) Atomic configurations (left panels) of twisted hBN bilayers with twist angles of (b) 24.0º, (c) 21.6º, and (d) 19.8º, and the corresponding diffraction patterns (right panels). The Miller indices, hkl, are shown to indicate the set of lattice planes responsible for each diffraction peak. Even in the incommensurate cases in (b,d), the partially eclipsed interlayer atomic overlap can be recognized. 12     Figure S6. Identification of moiré peaks in FFT. (ad) FFT images of the hBN layers with different twist angles of (a) 24.0°, (b) 21.6°, (c) 19.8°, and (d) 13.6°. Red and yellow circles show sets of 101̅0, 112̅0 and 202̅0 diffraction spots belonging respectively to different layers. White or black vectors are created by connecting the two diffraction spots belonging to the different layers. These vectors provide the positions of the moiré peaks when translated to the center of the diffraction pattern. 13      Figure S7. Fourier filtering process to reconstruct the moiré potential in real space by inverse FFT (IFFT). Left panels show original FFT images for the sample with a twist angle of (a) 24.0°, (b) 21.6°, (c) 19.8°, and (d) 13.6°. White circles in right panels indicate the diffraction spots used for the inclusive mask in the IFFT process. Red and yellow circles in right panels show sets of 101̅0, 112̅0 and 202̅0 diffraction spots belonging respectively to different layers. These well-defined Bragg diffraction spots as well as the center spot were excluded in the IFFT process.  14       Figure S8. Comparison of the IFFT images between the normal and inverted contrast conditions. (ad) IFFT images of the hBN layers with different twist angles of (a) 24.0°, (b) 21.6°, (c) 19.8°, and (d) 13.6°. Left and right panels show the normal and inverted intensities, respectively. 15       Figure S9. Comparison of the wide-area IFFT images between the commensurate and incommensurate layers. (a,b) (upper panels) IFFT images of the hBN layers with different twist angles of (a) 21.6° and (b) 24.0°. (lower panels) Line profiles as a function of distance along the white broken line in the corresponding upper panel.  16                                          Figure S10. Cross-sectional ESP map perpendicular to the atomic planes. (a,b) (top panel) Optimized structure of the model clusters with twist angles of (a)  = 0.0º and (b) 21.6º. (bottom panel) The corresponding cross-sectional ESP contour map along the red line shown in the top panel. The blue and pink lines are positive and negative, respectively. Contour intervals are drawn at 0.002 e (a.u.)3. In the bottom panel in (b), negative potential regions are marked by red circles, which correspond to the N on N (or N nearly on N) regions shown in the top panel.  . 17                                          Figure S11. Comparison of the optimized structural parameters obtained for the three BN bilayer model clusters with different size and different twist angles ( = 0° and 24.0°) calculated at the B3LYP-D3/6-31G(d) level. (a) (B75N75H30)2 with zigzag edges. (b) (B111N111H42)2 with armchair edges. (c) (B147N147H42)2 with zigzag edges. The optimized parameters for the BN bond distance dBN and the interlayer distance dint are given. Schematic representations of the clusters illustrate the optimized geometry although the terminal H atoms are omitted for clarity. 18                                          Figure S12. Density functional dependence of the optimized structural parameters calculated for the (B111N111H42)2 clusters with different twist angles. (a) = 0°. (b)  = 13.6°, (c) = 19.8°, (d)  = 21.6°, and (e)  = 24.0°. 19            Figure S13. Density functional dependence of the electrostatic potential for the (B111N111H42)2 clusters with different twist angles. (a)  = 21.6°, (b)  = 13.6°. The 2D plot of the ESP is shown for a plane bisecting the two hBN layers of the respective clusters. The unit of the scale bar is e a.u.1.   20  S6. Supporting references  (58)  Zobelli, A.; Gloter, A.; Ewels, C. P. ; Seifert, G.; Colliex, C. Electron Knock-on Cross Section of Carbon and Boron Nitride Nanotubes. Phys. Rev. B 2007, 75 (24), 245402. (59)  Jeon, J.W.; Kim, H.; Kim, H.; Choi, S.; Kim, B. H. Experimental Evidence for Interlayer Decoupling Distance of Twisted Bilayer Graphene. AIP Advances 2018, 8 (7), 075228. (60)  Liao, M.; Wei, Z.; Du, L.; Wang, Q.; Tang, J.; Yu, H.; Wu, F.; Zhao, J.; Xu, X.; Han, B.; Liu, K.; Gao, P.; Polcar, T.; Sun, Z.; Shi, D.; Yang, R.; Zhang, G. Precise Control of the Interlayer Twist Angle in Large Scale MoS2 Homostructures. Nat. Commun. 2020, 11, 2153. (61)  Kamiya, M.; Tsuneda, T.; Hirao, K. A Density Functional Study of van der Waals Interactions. J. Chem. Phys. 2002, 117 (13), 6010–6015. (62)  Becke, A. D. Density‐Functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98 (7), 5648–5652. (63) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate ab initio Parameterization of Density Functional Dispersion Correction (DFT-D) for the 94 elements H-Pu. J. Chem. Phys. 2010, 132 (15), 154104. (64)  Chai, J.-D.; Head-Gordon, M. Long-Range Corrected Hybrid Density Functionals with Damped Atom–Atom Dispersion Corrections. Phys. Chem. Chem. Phys. 2008, 10 (44), 6615-6620. (65)  Pease, R. S. Crystal structure of boron nitride. Nature 1950, 165, 722–723. (66) Warner, J. H.; Rümmel, M. H.; Bachmatiuk, A.; Büchner, B. Atomic Resolution Imaging and Topography of Boron Nitride Sheets Produced by Chemical Exfoliation. Nano Lett. 2010, 4 (3), 1299-1305. (67)  Jensen, F. Introduction to Computational Chemistry, 3rd ed.; Wiley: Chichester, 2017; p 321.  (68) Humphrey, W.; Dalke, A;. Schulten, K. VMD: Visual Molecular Dynamics. J. Mol. Graph. 1996, 14 (1), 33–38. (69) Frisch, A.; Nielsen, A. B.; Holder, A. J. Gaussview Users Manual; Gaussian Inc.: Pittsburgh, PA, 2000.                 accepted_version_manuscript SupportingInformation_3