Introduction

Following the discovery of graphene, the appearance of 2D transition metal dichalcogenides (TMD) significantly broadened the knowledge in the field of 2D materials, as well as opening potential optoelectronic applications. Besides this, vertical stacks of two TMD monolayers (ML) forming a bilayer demonstrate exciting optoelectronic properties, not present in individual MLs1. In such TMD bilayers, the two constituent MLs may possess different crystal directions creating in their overlapping region a moiré superlattice. The relative direction between the two MLs is called the twist-angle.

Photoluminescence (PL) studies on MoS2 TMD bilayers with different twist-angles reveal an interlayer electronic coupling, which corresponds to an indirect bandgap recombination which varies with twist-angle2. Recently, it was also discovered that the moiré periodic potential in twisted MoS2 bilayer can modify the properties of phonons in the respective ML constituents to generate Raman modes related to moiré phonons3. More recent studies revealed the presence of moiré excitons in twisted TMD homo- and hetero-bilayers4,5,6. The moiré pattern in the crystal symmetry of a twisted bilayer can be controlled through the rotation of the adjacent layers. In this context, the twist-angle is regarded as a new degree of freedom, enabling tuning of the physical properties of the TMD superlattices. Consequently by tuning the twist-angle in real-space, the change in the moiré pattern and consequently to the moiré periodic potential, one can control the interlayer coupling in order to obtain the desired superlattice properties. Obviously, the precise characterization of the twist-angle in moiré superlattices is essential for a global understanding and quality control in such 2D material systems, as well as precise tuning of the respective vdW devices’ performance.

Although TEM is the most commonly used technique to atomically reconstruct twisted TMD bilayers7, it requires tedious sample transfer on TEM grids, which is incompatible with most 2D materials fabrication and characterization techniques. Apart from being technically challenging, this procedure might eventually distort the relative lattice direction and alignment. Twisted TMD bilayers have also been imaged using atomic force microscopy (AFM)8, but this requires direct contact with the active area of the bilayer, thus entailing the risk of damaging the sample. While scanning electron microscopy (SEM) techniques do not generally suffer from these limitations, conventional SEM techniques used for crystallographic imaging rely on the detection of backscattered primary electrons, which is not enough to probe mono- or bilayer materials9.

As far as the optical techniques are concerned, up to date, the estimation of twist-angle in 2D TMD bilayers has been based either on simple optical microscopy observations, or on the production of SHG signals5,6. However, the approaches reported to date do not exhibit high accuracy and cannot image the twist-angle over extended bilayer areas. Having a tool that spatially resolves, with high-resolution and minimally invasively, the twist-angle in large area vdW heterostructures, would be therefore of great importance in the quality characterization of such structures. In this work, we demonstrate such a technique based on the areal imaging of polarization-resolved SHG (P-SHG) signals from TMD superlattices complemented with theoretical modeling that predicts the SHG signals interference from the respective overlapping areas of twisted-bilayers.

2D TMD materials like WS2 MLs belong to the D3h point symmetry group with broken inversion symmetry along the armchair direction. This lack of inversion symmetry in the TMD ML results in coherent SHG signals, when an intense field is incident on the 2D material10,11. The non-centrosymmetry which creates the SHG signals originates from the honeycomb lattice of the WS2 2D crystal, because of the alternating S and W atoms (top view in bilayer atoms configuration of Fig. 1a). The positions of the alternating S and W atoms define the broken symmetry axis of the crystal, which lies in the armchair direction. In contrast, the alternating S and S or W and W atoms define the zig–zag direction of the crystal.

Fig. 1: SHG signals originating from WS2/WS2 bilayers.
figure 1

a Schematic representation of the top view in the atomic configuration of a 2D WS2 bilayer, for the AA stacking sequence. Atoms of W–W or S–S are on top of each other and the SHG signals depend quadratically on the number N of the MLs (SHG~N2). b Coordinate system of the experimental configuration used for P-SHG imaging. θ1 and θ2 denote the armchair directions of ML-1 and ML-2, respectively. φ indicates the direction of the excitation linear polarization and ζ the axis of the analyzer. In our experiments ζ is fixed to 0(°) which greatly simplifies the experimentally retrieval of the P-SHG polar diagrams. c Atomic configuration for the AB stacking sequence. Alternating W–S atoms are on top of each other and the SHG signals cancel (SHG = 0), because centrosymmetry is restored. d Block diagram of the experimental setup used for P-SHG imaging. HWP: half-wave plate, GM: galvanometric mirrors, M: mirror, O: objective, SP: sample plane, C: condenser, SPF: short pass filter, BPF: bandpass filter, A: analyzer, PMT: photomultiplier tube.

Full size image

Τhe generated SHG from a WS2 ML (crystal class D3h) is described by its corresponding susceptibility tensor, χ(2). In our approach, we rotate the direction of the linear polarization \(\varphi\) of the excitation field and we detect the SHG component parallel to the X-axis (ζ = 0(°) in Fig. 1b). Then, the recorded SHG from a WS2 ML is given by (a four-leaved rose-like, polar diagram)12

$$I_{\left( {{\rm{ML}}} \right)}^{2\omega } = \left[ {A{\rm{cos}}\left( {3\theta _1 - 2\varphi } \right)} \right]^2.$$
(1)

here \({{{\mathrm{A}}}} = E_0^2\varepsilon _0{\upchi}_{{{{\mathrm{xxx}}}}}^{(2)}\), with ε0 being the dielectric constant, E0 the amplitude of the excitation field and θ1 ∈ [0(°)–60(°)], defines the armchair direction of the ML modulo 60(°). This means that the armchair directions with θ1(°) + k60(°) in Eq. (1) (where k is an integer) will provide the same P-SHG polar diagram.

When two TMDs, i.e. two WS2 MLs, are stacked to form a bilayer, the SHG signals from each ML interfere and the total SHG intensity from the WS2/WS2 superlattice is described by

$$I_{\left( {{\rm B}{\rm I}} \right)}^{2\omega } = A^2\left[ {{\rm{cos}}\left( {3\theta _1 - 2\varphi } \right) + {\rm{cos}}\left( {3\theta _2 - 2\varphi } \right)} \right]^2$$
(2)

Now, one can use the concept of effective armchair direction θeff in the overlapping region of the two WS2 monolayers and express the total SHG intensity produced by the 2 MLs as7

$$I_{\left( {{\rm{BI}}} \right)}^{2\omega } = \left[ {A_{{\rm{eff}}}{\rm{cos}}\left( {3\theta _{{\rm{eff}}} - 2\varphi } \right)} \right]^2$$
(3)

where

$$A_{{\rm{eff}}} = 2A{\rm{cos}}\frac{3}{2}\delta ,$$
(4)

where δ = θ1 − θ2 is the twist-angle, between the MLs and

$$\theta _{{\rm{eff}}} = \frac{{\theta _1 + \theta _2}}{2}.$$
(5)

As a result, the P-SHG modulation emerging from a bilayer region consisting of two WS2 MLs, at twist angle δ = θ1 − θ2, behaves as if it was the P-SHG modulation of a single ML with armchair direction θeff. The θeff can be extracted experimentally from the P-SHG polar obtained from the bilayer region. In the above, the excitation field propagates along the Z-axis and the x1 armchair direction of WS2 ML-1 is at angle θ1 with respect to the X-axis (Fig. 1b). While the second WS2 ML-2 has its armchair direction lying in the x2 direction, at angle θ2 with respect to the X-axis (Fig. 1b).

Note in Eq. (4), that the SHG intensity from the twisted bilayer depends on the twist-angle δ, being maximum for δ = 0(°) and zero for δ = 60(°). When the two WS2 MLs are perfectly aligned, e.g. for (δ = 0(°)), we have the AA stacking sequence (S or W atoms in one layer lie respectively above the S or W atoms of the second layer (S–S, W–W in Fig. 1a). In this case, the total SHG signal from the bilayer is the result of constructive interference, resulting in SHG intensity four times larger than that of the ML (for N number of MLs the produced SHG signal is analogous to N2). In the case of AB stacking sequence (δ = 60(°)) the S or W atoms in one layer lie, respectively, above the W or S atoms of the second layer (S–W, W–S in Fig. 1c). In this case, centrosymmetry is restored and the SHG signal from the bilayer vanishes (SHG = 0).

The above considerations, refer to the ideal cases of complete constructive or destructive interference in AA and AB bilayer stacking sequences, respectively. However, deviations from these ideal stacking sequences, can occur13. A direct consequence of such departures from the stacking sequences AA and AB is the incomplete constructive or destructive interference of the SHG from different MLs. In this case, the P-SHG signal modulation depends on the twist-angle between the MLs as shown in Fig. 2.

Fig. 2: P-SHG interference signatures in twisted WS2/WS2 bilayers.
figure 2

a Schematic showing the dependence of P-SHG signal on the twist-angle for fixed armchair direction ML-1, (θ1 = 0(°)) and varying armchair direction θ2 of ML-2, θ2 ∈ [0(°)–50(°)], with step 10(°). Top view atom configuration hexagonal lattice of the bilayer and simulated polar diagrams of the interference P-SHG originated from the WS2/WS2 bilayer (in red), as well as the simulated P-SHG signal modulation of each individual WS2 ML (green for ML-1 and blue for ML-2). b P-SHG signal modulation for fixed θ1 = 0(°) and varying θ2 ∈ [60(°)–110(°)], with step 10(°). By comparing a with b, we notice that the same green and blue P-SHG polar diagrams can provide two different interference P-SHG red polar diagrams. Therefore, we are able to define unequivocally the armchair direction θ2 that contributes to the SHG signal of the bilayer in the range of 0(°) ≤ θ2 ≤ 120(°).

Full size image

The graphical representation of Eq. (1) and the corresponding visualization in a polar diagram (presented in Fig. 2) demonstrates a fourfold symmetry of the P-SHG intensity modulation that rotates for different armchair directions θ. Thus, each armchair direction corresponds to a characteristic fourfold symmetric (four-leaved rose-like) polar diagram. Consequently, we can calculate θ by fitting Eq. (1) to experimentally retrieved P-SHG signal intensities with analyzer parallel to the X-axis and different linear excitation directions \(\varphi\). The experimental configuration of a fixed analyzer parallel to the X-axis used here (ζ = 0(°) in Fig. 1b), is much simpler to that of a rotating analyzer parallel to the rotating linear excitation polarization (ζ = φ in Fig. 1b), used previously6.

The green polar diagrams in Fig. 2 correspond to θ1 = 0(°) (created using Eq. (1)), while the blue polar diagrams correspond to 0(°) ≤ θ2 ≤ 110(°) with a step of 10(°) (created again using Eq. (1)). The red polar diagrams correspond to the product of interference P-SHG in the bilayer region (created using Eq. (2)).

By comparing Fig. 2a, with Fig. 2b, we note the effect of the modulo of 60(°) in the calculation of the armchair direction of individual MLs. That is, individual MLs with armchair directions θ(°)+k60(0), where k is an integer, will provide the same P-SHG polar diagrams. This implies that P-SHG measurements can calculate the armchair direction of an individual ML in the range 0(°)–60(°). Nevertheless, in the case of the bilayer, the SHG signals originate from constructive or destructive interference due to the atomic phase matching between the individual MLs (Fig. 2). Τhus, the produced P-SHG modulation from the bilayer (calculated using Eq. (2)) dictates a new angle range for the armchair direction θ2 of the second ML-2, i.e. 0(°)–120(°).

Additionally, as we note in Fig. 2, the SHG signals originating from the bilayers regions with twist-angles of 10(°) and 110(°), 20(°) and 100(°), 30(°) and 90(°), 40(°) and 80(°), are of equal SHG intensities. Nevertheless, the P-SHG interference polar diagrams from the bilayers’ regions (red lines in Fig. 2) for the twist-angles of 10(°) and 110(°), 20(°) and 100(°), 30(°) and 90(°), 40(°) and 80(°) are different, thus P-SHG is able to identify and discriminate twist-angles that produce similar SHG intensities from the bilayer regions.

Results and discussion

PSHG measurements

In order to create the WS2/WS2 bilayer, two WS2 monolayers were produced by mechanical exfoliation and were stacked with dry stamping on a Si3N4 support-grid (see Methods). This allows direct comparison between P-SHG and 4D STEM. The excitation source used for the SHG experiments is an fs oscillator, at 1030 nm and repetition rate in the order of MHz, which is adequate to excite non-linear signals like SHG (see Fig. 1b and Methods for the P-SHG microscope).

In Fig. 3a, a CCD image of two individual WS2 MLs (ML-1 and ML-2), forming a WS2/WS2 bilayer on the supporting TEM grid, is shown. One can identify two different types of regions, created by the above stacking, which produce different SHG signals (Fig. 3b, c). Specifically, we identify regions where only one of the two ML-1, and ML-2, WS2 is present, as well as the region where the two MLs spatially overlap and create the WS2/WS2 bilayer. We note that, in the absence of an analyzer in the detection path, the SHG intensity from the bilayer region is lower than the SHG from the MLs regions.

Fig. 3: P-SHG microscopy of WS2/WS2 bilayer on a TEM grid.
figure 3

a Wide field microscopy using the CCD camera of the microscope. We note two individual MLs of WS2 (ML-1 and ML-2), overlapping in a bilayer WS2/WS2 region. b SHG imaging of the same region seen in a., without the use of the analyzer in the detection path (see Fig. 1b, d). This ensures that all of the produced SHG signals will be acquired in the same image. Scale bar indicates 5 μm. c A zoom of the SHG image seen in (b). Three pixels-of-interest (POI) are chosen. One from the WS2 ML-1, one from the WS2 ML-2 and one from the WS2/WS2 bilayer region. Scale bar shows 1 μm. d P-SHG images of the region seen in c. for fixed analyzer (ζ = 0(°)), indicated by a blue double arrow and rotating linear excitation polarization (φ ∈ [10(°)–90(°)], with step 10(°)), indicated by a red double arrow. e, f Experimentally retrieved P-SHG polar diagrams for the POIs 1,2, respectively. The polarization rotates with φ ∈ [0(°)–360(°)], and step 1(°) and the analyzer is fixed at ζ = 0(°). The solid lines (green for ML-1 and blue for ML-2) are the fittings of Eq. (1) to the P-SHG data. This resulted to armchair directions θ1 = 44.69(°), R2 = 0.92 and θ2 = 29.17(°), R2 = 0.96 for ML-1 and ML-2, respectively, where R2 denotes the quality of the fitting. g Experimentally retrieved P-SHG data for POI 3 in the WS2/WS2 bilayer region using the same rotating polarization φ ∈ [0(°)–360(°)], with step 1(°) and fixed analyzer at ζ = 0(°) as in (e, f). The solid red line is the fitting of Eq. (3) to the P-SHG data and θeff = 6.88(°), R2 = 0.92 is the result of the fitting.

Full size image

In order to calculate the twist-angle of the bilayer, we utilize high-resolution P-SHG measurements with a step of φ = 1(°). Figure 3d presents the respective P-SHG images for fixed analyzer (ζ = 0(°)) and varying direction of the excitation linear polarization for φ ∈ [0(°)–360(°)]. The corresponding polar diagrams (Fig. 3e, f) obtained from two points of interest (POIs 1,2 in Fig. 3c), one for each ML, are fitted to Eq. (1) in order to calculate the armchair direction of the MLs. Then, by using the effective armchair θeff obtained from the bilayer (POI 3 in Fig. 3c), we can calculate the twist-angle in the superlattice (using the experimentally retrieved polar of Fig. 3g). We find that θ1 = 44.69(°), R2 = 0.92 in the ML-1 region and θeff = 6.88(°), R2 = 0.92 in the bilayer region.

In previous reports, calculation of twist-angle is performed by measuring the armchair direction of the individual MLs, outside the bilayer region and then deducing their twist angle in their overlapping area5,6. If we follow this strategy we obtain δ = θ1 − θ2 = 44.69(°) − 29.17(°) = 15.52(°). Nevertheless, as we see in Fig. 2a, a twist-angle of 15.52(°) should have resulted to constructive SHG interference, i.e. the SHG signal in the bilayer region should have been more intense than the SHG signals from the MLs. This is not the case in our experimental data, where the SHG signal in the bilayer region is less intense than the individual MLs (Fig. 3e–g). Consequently, following Fig. 2b we should use the modulo 60(°) in the calculation of θ1 and obtain, θ1 = 44.69(°) + 60(°) = 104.69(°). As a consequence, we unequivocally determine the twist-angle δ = θ1 − θ2 = 104.69(°) − 29.17(°) = 75.52(°). We could also use the θeff acquired from the bilayer region and Eq. (5) to retrieve the twist-angle as δ = θ1 − θ2 = 2(θ1 − θeff) = 2(44.69(°) − 6.88(°)) = 75.62(°).

In Fig. 4b we present pixel-wise mapping of the armchair direction for the region seen in Fig. 4a. We used 360 P-SHG images, each one acquired for 1(°) rotating excitation linear polarization for φ ∈ [0(°)–359(°)] and we have fitted them pixel-by-pixel to Eq. (1). We choose to keep only the pixels that exhibit a quality of fitting, R2 > 0.8. As we have explained in12 the P-SHG signal fails to fit in the borders among the regions of different armchair orientations in Fig. 4b, c. As a consequence, the pixels located at these borders are missing from the armchair map.

Fig. 4: Optical imaging of twist-angle in WS2/WS2 bilayer.
figure 4

a SHG image without the use of analyzer, showing the two ML-WS2 regions and the bilayer WS2/WS2 region. We note that the bilayer region has less SHG signal than the MLs regions. Scale bar indicates 1 μm. b Mapping of the armchair direction, obtained via pixel-by-pixel fitting of 360 P-SHG images, one for every φ ∈ [0(°)–359(°)], with step 1(°) and fixed analyzer (ζ = 0(°)) to Eq. (1). Only pixels that presented R2 > 0,8 are kept. c Pixel-by-pixel mapping of the armchair direction for the ROI seen in (b). Three ROIs are defined. ROIs-1,2 for the WS2 MLs and ROI-3 inside the bilayer WS2/WS2 region. d Image histogram showing the distribution of armchair directions seen in (c) and in the ROI seen in (b). The mean armchair directions and their standard deviation (σ) of the armchair values found in ROIs1-3 are presented.

Full size image

In Fig. 4c we present a magnification of the ROI marked in Fig. 4b. We have also marked two ROIs-1,2, one for each ML-1,2, respectively, and ROI-3 in the bilayer region. In Fig. 4d we show the image histogram of the armchair directions present in Fig. 4c and in the ROI seen in Fig. 4b. In the same Fig. 4d we present the mean and the standard deviation (σ) values of the armchair directions found in ROIs 1-3. We found <θ1 > = 44.64(°) with σ = 0.04(°), <θ2 > = 44.64(°), σ = 0.06(°) and <θeff > = 6.96(°), σ = 0.55(°). The iterative fitting procedure resulted in small standard deviation values for ROIs 1,3 (Fig. 4c). The very small size of ROIs, smaller than the optical resolution of our 1.3NA objective (~500 nm for the 1027 nm excitation), and the oversampling with many measurements (pixels) inside those small size ROIs are responsible for the resulted small standard deviations. The broadening of the angle distributions seen in Fig. 4d is attributed to erroneous values originating from the borders of the different regions with different crystal orientations. This is because during the laser raster-scanning the spot of approximately 500 nm size when passing from those borders, excites SHG from both regions of different crystal orientation. If we follow the strategy of deducing the twist-angle by measuring the armchair direction of the individual MLs we acquire δ = 44.64(°) − 29.22(°)) = 15.42(°). Nevertheless, as we see in Fig. 2 such twist-angle value should have resulted in constructive SHG interference, i.e. the SHG signal in the bilayer region should have been bigger than the SHG signal of the MLs (Fig. 3e–g). This is not the case for our experimental data, therefore we should follow Fig. 2 and use the modulo 60(°) in the calculation of θ1 and obtain, θ1 = 44.64(°) + 60(°) = 104.64(°). This enables us to unequivocally determine the twist-angle δ = θ1 − θ2 = 104.64(°) − 29.22(°) = 75.42(°). We could also use the mean θeff acquired from ROI-3 in the bilayer region and Eq. (5) to retrieve the twist-angle as, δ = θ1 − θ2 = 2(θ1 − θeff) = 2(44.64(°) − 6.96(°)) = 75.36 ± 0.55(°).

4D STEM measurements

4D STEM was employed to measure the crystal directions of the individual monolayers, as well to independently estimate the twist-angle in the WS2/WS2 bilayer region seen in Fig. 3a. The schematic representation of the application of the 4D STEM method for studying WS2 monolayers is shown in Fig. 5a. The microscope settings are presented in the Methods section. In particular, the bilayer region that has been analyzed by P-SHG is scanned with 256 × 256 probe positions with a step size of 25 nm. A virtual dark field image (VDF), calculated by summing the intensity of diffraction spots in each probe position over the selected virtual aperture, is shown in Fig. 5b. The inner and outer radii of the virtual aperture are chosen to select the second-order diffraction spots of the WS2 structure (Fig. 5c). From the VDF intensity distribution, we can identify four distinct regions corresponding to the silicon nitride support-gd (dark contrast) (mean diffraction pattern is shown in Fig. 5c), two monolayers (Fig. 5b, e) and an overlapped region with a bilayer (Fig. 5f). At the edges of the monolayers there are regions with increased intensity due to their folding (white arrows in Fig. 5b).

Fig. 5: 4D STEM measurements.
figure 5

a Schematic representation of the application of the 4D STEM method for studying WS2 monolayers. b Virtual dark field image where the intensity in each probe position is summed over selected virtual aperture indicated in (c). Scale bar 1 μm. c–f The mean diffraction from the four regions (P1-4). c The first WS2 monolayer (P1). d Silicon nitride supporting film (P2). e The bilayer region P3. f The second monolayer P4. Scale bars, 0.5 Å−1.

Full size image

To calculate a direction map of the region, a custom-made peak finding routine was used. Once peak positions are found, two reciprocal vectors can be fitted to describe all diffraction spots defining the relative crystal direction in every probe position. This algorithm works only for diffraction patterns from monolayers. For the bilayer region, four vectors are fitted to each acquired pattern when more than six peaks are detected inside the virtual aperture area. The relative direction map of the region with a bilayer area (white square in Fig. 5b) is shown in Fig. 6a. The tilt angle of every diffraction pattern for two sheets is presented in the histogram Fig. 6d. The twist-angle between the two monolayers was calculated with subpixel accuracy by fitting a Gaussian function for each peak on the histogram and was found to be 15.5 ± 0.3(°). These results are in good agreement with the experiments on the same area using the all-optical P-SHG imaging microscopy technique, presented above.

Fig. 6: Twist-angle calculation.
figure 6

a The relative color-coded direction map of a bilayer region selected in Fig. 6b. b Histogram of the relative direction of the monolayers and the bilayer indicating that the two layers are stacked upon each other since its direction in the bilayer correspond to both monolayers. c, d The diffraction pattern averaged over 400 probe position for the first and second monolayer. e The annular intensity profile of the first order spots where a modulation of intensity is observed. f, g The simulated diffraction patterns of the WS2 bilayer when the twist-angle is 15.5(°) and 75.5(°) (60(°)+15.5(°)). h Experimental diffraction of the bilayer region, averaged over 400 probe positions. i The annular intensity profile of the first order spots of experimental diffraction where the intensity modulation of the two layers is reversed indicating the twist angle of 75.5(°) instead of 15.5(°).

Full size image

However, for the direction determination only the positions of the diffraction spots were used without taking into account their intensity. By looking at the diffraction spot symmetry, one might think that WS2 crystal have six-fold symmetry along the [001] direction, however, due to dynamic electron scattering in combination with the non-centrosymmetry of the crystal, Friedel’s law is violated resulting in a symmetry reduction to three-fold symmetry and a nonequivalent angular range from 0(°) to 120(°) instead of 0(°) to 60(°). Figure 6c, d shows diffraction patterns of the first and second monolayer together with their azimuthal intensity profiles of the first order spots (Fig. 6e) clearly showing 3 fold symmetry. This allows us to identify whether the twist angle is 15.5(°) or 75.5(°) (15.5(°) + 60 (°)) by comparing to simulated diffraction patterns shown in Fig. 6f, g. The simulation is performed using the Multem software14,15 where the experimental conditions such as convergence angle and electron beam energy are chosen to closely resemble the experimental setup (see Methods section). A noticable difference between the simulated and experimental patterns (Fig. 6h) is caused by the presence of the Si3N4-support film (experimental), which gives rise to an isotropic background signal. As can be seen in the simulated patterns the intensity of the first order spots changes azimuthally from two low, two high for the twist angle 15.5(°), while it has alternating low-high character for the 75.5(°) angle. The annular intensity profile of the first order spots of the experimental pattern (Fig. 6i), shows a clear signature of a 75.5 ± 0.3(°) twist angle and rules out the 15.5 ± 0.3(°) hypothesis. This result is in excellent agreement with the P-SHG data presented above.

In terms of precision, the STEM data demonstrates 0.3(°) rotational precision compared to 0.55(°) for the P-SHG optical result with a significantly higher spatial resolution but with the downside of a technique which is far less attractive to employ in an inline production environment as compared to the all optical setup proposed in this paper.

Although P-SHG is fast, minimally invasive and can cover large areas, it cannot offer the atomic resolution of the STEM technique. In addition, SHG requires breaking of the symmetry and therefore vanishes in intrinsically centrosymmetric systems (such as graphene) or in systems where centrosymmetry is restored e.g. due to the number of layers (i.e., 2H-stacked TMDs with even number of layers) or e.g to homo-bilayers with twist-angle of 60(°). For the case of encapsulated TMDs (e.g. between hBN layers) the collected SHG signal from the overlapping region may include additional contributions (as in the case of hBN which is a polar material), and an additional measurement in an uncovered region is needed to unequivocally determine the relative crystal orientation of the individual layers comprising the bilayer.

Considering that the electronic and optical properties of 2D TMD bilayers can be tuned by changing their twist-angle, a robust and minimally invasive tool that can provide spatially resolved determination of the twist-angle, would be of great importance in research, production and large-scale characterization of 2D TMD bilayers. Here, we have used all-optical, laser raster-scanning, P-SHG imaging microscopy to precisely map the twist-angle in large areas of overlapping WS2/WS2 stacked monolayers and we benchmarked the results against 4D STEM electron microscopy. It is found that the twist-angle of WS2/WS2 bilayer obtained using P-SHG mapping is in excellent agreement, either in value and in precision, with that obtained using 4D STEM. It is additionally revealed that, given that the produced SHG signal from a bilayer is the vectorial addition of the SHG signals of the individual monolayers, the intensity modulation of the P-SHG signal can be used to deduce unequivocally the armchair direction. This is also in excellent agreement with 4D STEM microscopy analysis.

While STEM provides significantly higher spatial resolution, the P-SHG compensates for this with a wide range of advantages including: no need for vacuum, wide field of view, rapid data acquisition, significantly lower cost and instrument size, and most importantly its capability to work on TMD bilayers deposited on substrates without the need to transfer the films to a TEM grid. Our setup provides an accurate and robust all-optical twist-angle mapping of 2D TMD bilayers. Importantly, the technique is non-destructive, paving the way or directly correlating local twist-angle values with electronic properties, which is crucial for the development and scaling up of vdW bilayer devices with precisely controlled functionality.

Methods

WS2/WS2 bilayer fabrication on a TEM grid

Polydimethylsiloxane films (PDMS) were fabricated from 10:1 mixing ratio (SYLGARD 182 Silicone Elastomer Kit) with heat cure at 80 °C for two hours. High quality WS2 bulk crystals (HQ Graphene) were mechanically exfoliated directly on the aforementioned PDMS films. The films were placed on typical microscope glass slides using standard protocol16,17. Monolayers of these crystals were realized under an optical microscope. In order to produce the bilayer, at first, a glass slide with a WS2 monolayer was mounted on a XYZ micromechanical stage under a coaxially illuminated microscope and transferred on a silicon nitride (Si3N4) support-grid using viscoelastic stamping18. Finally, another WS2 monolayer was stamped on the previous one in a partial overlapping manner allowing for comparative P-SHG and 4D STEM imaging.

Custom-built P-SHG microscope

SHG imaging was performed in the forward-detection geometry using a custom-built laser raster-scanning microscope7,11,12. As shown schematically in Fig. 1d, a diode-pumped Yb:KGW fs oscillator (1027 nm, 90 fs, 76 MHz, Pharos-SP, Light Conversion, Lithuania), was inserted into a modified, inverted microscope (Zeiss Axio Observer Z1, Germany), after passing through a pair of silver-coated galvanometric mirrors (6215H, Cambridge Technology, UK). A motorized rotation stage (M-060.DG, Physik Instrumente, Karlsruhe, Germany), holding a zero order λ/2 wave plate (QWPO-1030-10-2, CVI-Laser, USA) was used to rotate the direction of the excitation linear polarization. Then, the beam was reflected on a silver mirror at 45(°) (PFR10-P01, ThorLabs, Germany), placed at the turret box of the microscope, just before the objective (Plan Apochromat 40 × 1.3NA, Zeiss, Germany). The SHG signals were collected from a high-numerical aperture (1.4NA) condenser lens (Zeiss) and guided into a photomultiplier tube (PMT) detector (H9305-04, Hamamatsu, Japan) using another silver mirror (CM1-P01, ThorLabs, Germany) at 45(°). The polarization extinction ratio was 28:1. In front of the PMT, a home-built mount was holding a bandpass filter (FF01-514/3-25, Semrock, USA) and a short pass filter (FF01-680/SP-25, Semrock, USA), appropriate for SHG imaging. After the filters, a film polarizer (LPVIS100- MP, ThorLabs) was inserted just in front of the PMT to measure the anisotropy of the SHG signals due to the rotation of the excitation linear polarization. Coordination of PMT recordings with the galvo-mirrors movements and with all the motors, as well as the image formation, was performed using LabView (National Instruments, USA). All of our SHG images were of 500 × 500 pixels and approximately 1.1 s was required for each image to be recorded. This resulted in pixel dwell time of approximately 4.4 μs.

4D STEM measurements

In order to determine local information on the single-layer direction and the bilayer twist angle of 2D materials, a scanning transmission electron microscope ThermoFisher Scientific (FEI) Titan X-Ant-EM was used. The electron microscope was operated in microprobe STEM mode at 300 kV at a convergence semi-angle α of 1 mrad resulting in a probe size of 1.2 nm in diameter. Local information was obtained by scanning the electron probe over the sample and acquiring an electron diffraction pattern at every probe position19,20,21,22,23,24 on a Medipix3 hybrid pixel direct electron detector (Quantum Detectors Merlin)25,26,27 with a camera length of 115 mm and exposure time of 5 ms. In comparison to the SHG pixel dwell time of 4.4 μs, the 5 ms are three orders of magnitude slower. This detector offers a high frame rate and high efficiency enabling the detection of individual electrons without dark or read-out noise and offering 24 bit dynamic range. The diffraction data were processed by custom-made scripts based on the open-source python library Pixstem28.