Twist Angle mapping in layered WS2 by Polarization-Resolved Second Harmonic Generation

Stacked atomically thin transition metal dichalcogenides (TMDs) exhibit fundamentally new physical properties compared to those of the individual layers. The twist angle between the layers plays a crucial role in tuning these properties. Having a tool that provides high-resolution, large area mapping of the twist angle, would be of great importance in the characterization of such 2D structures. Here we use polarization-resolved second harmonic generation (P-SHG) imaging microscopy to rapidly map the twist angle in large areas of overlapping WS2 stacked layers. The robustness of our methodology lies in the combination of both intensity and polarization measurements of SHG in the overlapping region. This allows the accurate measurement and consequent pixel-by-pixel mapping of the twist angle in this area. For the specific case of 30° twist angle, P-SHG enables imaging of individual layers.

The graphene-related atomically thin 2D TMDs show great promise for high-tech optoelectronic applications [1][2][3][4] . In particular, they exhibit unique nonlinear optical properties owing to their reduced dimensionality and lack of centrosymmetry, that give rise to pronounced SHG [5][6][7][8][9] . Analysis of the emitted SHG signal provides information on the crystal orientation [5][6][7] and homogeneity 8,9 as well as the thickness 10 and stacking sequence 11,12 of TMD structures. At the same time, there has been an increasing scientific interest in twisted TMD structures that can either occur during chemical vapour deposition (CVD) growth or be prepared artificially. In the latter case, one can tailor the interlayer coupling, that is strongly twist angle-dependent, and thus reveal new physical phenomena [13][14][15][16] . Therefore, the twist angle can be regarded as a new degree of freedom, enabling tuning of the physical properties of stacked 2D materials. A prominent example was recently reported in a ground-breaking work, showing that bilayer graphene exhibits unconventional superconductivity for a small value of the twist angle between the layers 17 . It has also been demonstrated that the twist angle allows control of the valley and band alignment of stacked 2D TMDs 11 and it results in ultraflat bands and shear solitons in twisted bilayer MoS 2 18 , thus enabling ultrafast charge transfer between the 2D layers 19 . The strain induced SHG 20 , the twist angle-dependent moiré-templated strain patterning 21 , the interlayer valley excitons in TMD heterobilayers 22,23 , the twist angle-dependent conductivities across MoS 2 /graphene heterojunctions 24 , and the moiré excitons in heterobilayers [25][26][27][28] are just a few more studies that have also been reported recently. These observations indicate the strong potential to harness and tune the physical properties of layered 2D materials, via the adjustment of the twist angle of the stacked layers 29 . In this context, the development of an optical technique capable to map the twist angle with high precision over large areas would be an invaluable tool for the construction and characterization of such new materials.
Earlier studies on the SHG interference from artificially stacked TMD bilayers have shown that the differences in SHG intensity can be attributed to differences in the armchair orientation between the two twisted TMD layers 12 . In addition, phase-resolved SHG techniques have also been used for the determination of the relative orientation between monolayers 30 . In all these cases however, solely variations in SHG intensity were used to identify the armchair angle difference, a criterion that cannot unambiguously exclude other phenomena as a source of these changes, such as structural transformations and inhomogeneities. The novelty of our method lies in the combination of intensity and polarization-resolved SHG measurements in the overlapping region of stacked 2D materials. SHG intensity-only measurements are insufficient for the determination of the twist angle since variations in intensity may also be caused by changes in the stacking sequence of the layers (e.g. from 2 H to 3 R). On the other hand P-SHG modulation may be due to imperfections of the crystal quality that result in local changes of the main crystallographic axis 9 . Our work aims to resolve this issue by offering, for the first time, a combined SHG intensity and P-SHG study of stacked 2D structures, that allows the pixel-by-pixel mapping of the twist angle.
As a proof of concept, we demonstrate the advantages of our all optical nonlinear imaging method, by application both to a pair of overlapping CVD-grown WS 2 layers with different armchair orientation, and an artificially dry stacked WS 2 twisted bilayer that has been produced with mechanical exfoliation. Interestingly, we show that there is a specific twist angle of 30°, for which the SHG signal originating from the stacked single-armchair layers (SLs) can be selectively switched-on and -off. This enables the effect of optical discrimination of atomically thin layers and therefore provides a form of axial super-resolution SHG imaging of each individual layer of stacked 2D TMDs.

Results and Discussion
Within our approach, the SHG emission from a layered TMD is considered as the result of interference between SLs. As a SL is considered a 2D layer of 2 H or 3 R stacking sequence 29,31 . Each one of those SLs is acting as a surface phase array antenna inside the SHG active volume 32 . In order to examine the effect of the twist angle of the stacked SLs on their combined SHG pattern, we have performed polarization-in, polarization-out SHG measurements on a raster-scanned area of a CVD-grown WS 2 sample and an artificially prepared WS 2 /WS 2 bilayer.
In our recent work on P-SHG of TMDs, we have demonstrated that the SHG signal modulates as the angle of linear polarization of the excitation field, ϕ, rotates and that the modulation depends on the armchair orientation. Using this P-SHG modulation one can map with high precision (~0.15°) the armchair angle, θ, at every pixel of the image 9 . In the case of two SLs that partially overlap, the P-SHG from each individual layer modulates with a phase depending on the armchair orientation of the layer, while in the overlapping area, the P-SHG modulation follows the SHG interference of the two SLs. Based on this effect, it can be shown that upon using appropriate linear polarization for the excitation of layered WS 2 with two different armchair angles overlapping in a region, we are able to decompose the interfered SHG signal originating from each individual layer. This is shown in Fig. 1 and Video S1, presenting the P-SHG modulation measured from a multilayered CVD-grown WS 2 triangular flake that exhibits a central region of enhanced SHG intensity. In particular, by rotating the linear polarization of the excitation field, incident to the overlapping area of two SLs, we are able to switch-on the SHG signal from one SL, while the SHG signal from the other is switched-off and vice versa. As shown in detail below, such SHG switchon/off effect occurs when the armchair angle difference, i.e. the twist angle, between the layers equals to 30°.
At the stacking level, the produced SHG originates from the interference between at least two stacked SLs. It is therefore determined by the relative orientation of the two SLs, i.e., the twist angle. In a CVD-grown TMD, as in the case of our WS 2 sample, the crystal could be a mixture of 3 R and 2 H phases 33 . A direct consequence of such deviation from the ideal stacking sequence is the incomplete constructive or destructive interference of the SHG fields from different layers. Furthermore, the surface dipoles of the SLs are misaligned and the total SHG signal, as well as its polarization, depend on the twist angle. Regardless such crystal structure deviations, here we show that by using P-SHG measurements one can map for each pixel the armchair orientation of each SL constituting the multilayered structure. In order to accomplish that, the measurement utilizes the rotation of the linear excitation polarization with respect to the X-lab axis and a polarization analyzer parallel to X-lab axis, prior to SHG signal detection (Fig. 2). The armchair orientation of the first SL (x 1 y 1 z 1 coordinate system) is at an angle θ 1 with respect to the X-axis, whereas the armchair orientation of the second SL (x 2 y 2 z 2 coordinate system) is at an angle θ 2 with www.nature.com/scientificreports www.nature.com/scientificreports/ respect to the same axis ( Fig. 2). In the polarization-in, polarization-out measurements the propagation of the laser beam is along Z = z 1 = z 2 .
Based on our previous work 9 , it is straightforward to obtain an expression for the detected SHG signal in the case of two overlapping SLs (Fig. 3). The total SHG field in the overlapping area is given by the vector sum of the SHG signal from each layer. For the case of the analyzer orientation parallel to the lab X-axis and rotating linear polarization ϕ the recorded SHG intensity is described by: Here, θ 1 , θ 2 correspond to the armchair orientations of the SL1 and SL2 respectively, and , with C i being constants that depend on the local fields and the number of monolayers comprising the SL 33 .
, in the overlapping region of two monolayers we obtain SHG equal to the SHG from each individual monolayer. This results in a form of axial super-resolution imaging of individual 2D layer and occurs at the 'magic'-SHG twist angle: δ θ θ = − = ±°30 1 2 . The total detected SHG intensity from both SLs, is also given by: where θ = θ θ + eff 2 1 2 is the effective armchair orientation in the overlapping region. Since θ θ , Note that the SHG intensity depends on the armchair angle difference δ between the SLs, being maximum for δ = 0°and zero for δ = 60°, while for δ ≤ 30° and δ ≥ 30° we have partially constructive and partially destructive SHG interference, respectively. The SHG of N number of such SLs for the case of the analyzer orientation parallel to X-axis and rotating linear polarization ϕ is described by: photomultiplier tube; XYZ: Lab coordinate system; x 1 y 1 Z and x 2 y 2 Z: coordinate systems of SL1 and SL2, respectively; θ 1 , θ 2 : angles between XYZ and x 1 y 1 Z, x 2 y 2 Z, respectively; ϕ: angle between X-axis and excitation linear polarization E(ϕ). We note three distinct regions: clear SL1, clear SL2 and their overlapping region producing SHG 1 (θ 1 , ϕ), SHG 2 (θ 2 , ϕ) and SHG 1+2 (θ 1 , θ 2 , ϕ), respectively.
We can therefore employ Eq. (3) with N = 2, using the total SHG from a pair of layers, to deduce the actual armchair orientation of the second SL by considering the first SL as a reference. This calculation can be performed provided that there is non-zero SHG signal and that the SLs are built as shown in Fig. 3, where part of the first SL at the bottom (SL1) is not overlapping with the second (SL2). In particular, the P-SHG signal detected from the uncovered region of SL1 can be used to calculate the armchair θ 1 for SL1 (Eq. (3) with N = 1). Subsequently, the armchair orientation θ 2 of the second SL2 (red area in Fig. 3), can be derived upon using the known θ 1 and Eq.
(2). In this case, the resolution of θ 2 , calculated as error propagation of θ 1 resolution 0.15°, is ~0.33°. In the same manner, by knowing θ 1 , θ 2 , the armchair of a third SL3 can be computed upon using again Eq. (3) for N = 3. This procedure can be repeated for an arbitrary number N of SLs, provided that there is always a non-overlapping region among the stacked layers (Fig. 3).
For example, in Fig. 4, we present the theoretical P-SHG modulation of two stacked SLs, as well as the modulation of their SHG interference signal in the overlapping region, for several twist angles of interest, assuming A 1 = A 2 (i.e., layers of equal SHG intensity). This assumption is realistic since layers of the same symmetry and similar composition should possess similar χ (2) . The blue polar diagrams shown in Fig. 4 are fixed and correspond to the SHG from an individual SL (Eq. (3) for N = 1, A 1 = 1 and θ 1 = 0°), whereas the green polar diagrams show the SHG from a second SL (Eq. (3) for N = 1, A 2 = 1, and varying armchair orientation θ 2 . Finally, the red polar diagrams in Fig. 4 correspond to the SHG interference from the two SLs in their overlapping area.
Note that for δ = ±°40 , in Fig. 4, all three polar diagrams, corresponding to the SHG signal from SL1, SL2 and their overlapping area, reach the same maximum intensity. Upon rotating the excitation linear polarization E(ϕ), the SHG from SL1 (blue polar) reaches its maximum first, while the SHG signals from both SL2 (green polar) and the overlapping region (red polar) are zero. Subsequently, the SHG from both SL1 and SL2 goes to zero while that from the overlapping region maximizes. Finally, both the SHG from SL1 and the overlapping region go to zero as the SHG from SL2 becomes maximum. This means that in the case of two overlapping layers at δ = ±°40 , one could selectively: (i) switch on the SHG from SL1, while that of the overlapping region and SL2 are switched-off, or (ii) switch-on the SHG from the overlapping region, while that of SL1 and SL2 are both switched-off, or (iii) switch-on the SHG from SL2, while that of the overlapping region and SL1 are both switched-off. Additionally, we note in Fig. 4 that when δ = ±°30 the polar diagram produced by the interference of the two SLs (red line), passes from the maximum of the two polar diagrams produced by each individual layer SL1, SL2 (blue and green lines, respectively). Note also that when the SHG from the overlapping region is equal to the SHG from SL1 (SL2), the SHG from SL2 (SL1), is zero. This results in the complete switching-on of the SHG from one layer, while the SHG from the other layer is completely switched-off. This is also confirmed experimentally below, as for δ = ±°30 one can selectively switch-on the SHG from the SL of preference.
Based on the above analysis, the amplitude of the P-SHG modulation (see Eq. (2)) originating from layered regions depends on both the number of SLs and their relative armchair orientation, i.e., the twist angle δ. Consequently, a change in the SHG amplitude in a SL is an indication of either the presence of a second TMD SL, or a change in the stacking order of the same SL (e.g. from 2 H to 3 R stacking 34,35 ). Figure 5a for ϕ = 64 ο shows the SHG signal originating from one SL for the region of interest (ROI) shown in Fig. 1a, whereas Fig. 5b for ϕ = 17 ο presents the SHG from an assumed second SL (see Eq. (2)). Finally, Fig. 5c is the summation of the two SHG images. In order to interpret the SHG signal variations and identify the existence of a second SL we perform a SHG-intensity profile analysis along the line-of-interest (LOI) in Fig. 5a-c. In particular, in Fig. 5d-f we plot the intensity profile for each pixel along this line. In Fig. 5a we mark with a white dashed line the region where the intensity drops close to zero. We set that lack of SHG signal as reference corresponding to the substrate. In Fig. 5b we mark with a light blue dashed line the region of the assumed second SL2. Finally, in Fig. 5c we note the overlapping region, contained within the yellow dashed line, as well as the empty and SL2 areas enclosed within the white and light blue dashed lines, respectively. We finally note the yellowish change in the colour in the overlapping region caused by the addition of green and red colours of Fig. 5(a,b), respectively. In Fig. 5d the SHG intensity everywhere along the LOI corresponds to that of one SL (~150 a.u.), except from the dark region (ranging from 200 to 350 pixels) where the material is absent and the SHG drops to zero. In Fig. 5e the SHG intensity again corresponds to one SL (~150 a.u. for pixels ~200-~500). The summation of the profile intensities shown in Fig. 5d,e, results in the profile intensity shown in Fig. 5f, where it is obvious that for pixels 0-~350 and ~500-~600 the SHG intensity corresponds again to one SL (~150 a.u.). Interestingly, the SHG intensity in the region where it is assumed that the two SLs overlap (pixels ~350-~500), is double (~300 a.u.) with respect to the one SL region, indicating the presence of the second SL, or a change in the stacking order. A change from 3 R to 2 H stacking is excluded, considering that it would give rise to a reduction of the SHG signal intensity 32 . It is therefore concluded that the increase of the SHG signal observed in the overlapping region, is a signature of an additional SL or a change to 3 R stacking order (presence of second SL with equal armchair orientation) within the same layer.
In order to rule out the possibility of a change in the stacking order, intensity-only measurements are not sufficient and one has to perform a P-SHG analysis. Figure 6a-d present the corresponding P-SHG polar plots for the same ROI as above and the armchair angles for points of interest (POI) lying in the different regions (SL1, SL2 and overlapping). Indeed, from the P-SHG signal obtained from POI1, represented with the blue polar diagram in Fig. 6b, we can calculate, using Eq. (3) with N = 1, the armchair angle of the bottom SL1 layer to be θ 1 = 39.5°. Using Eq. (2) with θ 1 = 39.5°, we obtain for POI2 of Fig. 6a the green polar shown in Fig. 6c, which corresponds to armchair angle of θ 2 = 12.8° for the SL2. Notably, the P-SHG signal from POI3 corresponds to the data points and www.nature.com/scientificreports www.nature.com/scientificreports/ the fitted red polar in Fig. 6d which leads, using Eq. (2) to θ eff = 25.1°. This signal is the result of the interference between the two SLs in POIs 1 and 2. By using the measured θ eff = 25.1° and by fixing to the measured θ 1 = 39.5°, we can calculate, using θ eff = (θ 1 + θ 2 )/2, that θ 2 = 2θ eff − θ 1 = 10.7°, which corresponds to the armchair angle of the second SL that produces the SHG interference (see Video S2). Consequently the twist angle can be calculated using two measurements, one in the overlapping and another in the monolayer region, as δ = 2(θ 1 − θ eff ) = 28.8°.
From the above observations we can exclude the possibility of a change in the stacking order within the same SL in the overlapping region, since in that case θ eff should be equal to θ 1 (3 R stacking occurs in lattices of the same armchair orientation). Indeed, using P-SHG measurements we calculate a θ eff that indicates the presence of a second SL with armchair orientation that fits very well with that calculated for the overhanging SL2 region. These www.nature.com/scientificreports www.nature.com/scientificreports/ results denote a continuous region of the same armchair orientation and thus the presence of a second SL. This region can also be optically isolated in the SHG image of Fig. 5b. It should be noted here that the four-leave rose patterns shown in Fig. 6b-d correspond to the P-SHG signatures of WS 2 9 . Therefore, any detected SHG signal that does not comply with this pattern modulation, is excluded during the data fitting.
It should also be emphasized that the experimental twist angle calculated above in the overlapping region, δ = 28.8°, is very close to the theoretical one δ = 30° required for complete suppression of the SHG signal from one SL. Confirming our theoretical prediction (see Fig. 4), the experimentally inferred δ = 28.8°, results in almost complete suppression of the SHG from one layer, while the SHG of the other is maximum, as it is experimentally demonstrated in Fig. 1 and Video S1. This is also shown in Fig. 6d and Video S2, where the experimental polar diagram produced by the interference of the two SLs (red curve) passes very close to the maximum of the two polar diagrams produced by each individual layer (blue and green curves). Note that when the SHG interference intensity is equal to the maximum of the SHG intensity from one SL, the SHG intensity from the other SL is almost zero. This results in the suppression of the SHG from the first layer, while at the same time the SHG from the second layer is completely switched-on and vice-versa. It is stressed out that the second SL was not intentionally placed in the ideal 30° twist angle but it was naturally grown during the CVD, measured at δ = 28.8° twist angle.
In this context, we can interpret in detail the P-SHG map obtained from the multilayered CVD-grown WS 2 triangular flake shown in Fig. 1 as follows. Figure 6e presents the map of armchair angles θ for the ROI in Fig. 1a. This map is obtained by performing pixel-wise fitting using Eq. (3) for N = 1 (assuming only one layer). Using this map one can create the armchair histogram of specific areas, for example the one shown in Fig. 6g. Notably, the armchair mapping gives values of θ ~ 40 ο for the outer region, but only θ ~ 25 ο for the central one (yellowish). In addition in Fig. 6e, there is a part of the central region that exhibits θ ~ 12 ο (blue). The above results confirm the presence of more than one SL. If we now assume two different SLs, fix the measured θ 1 = 39.5° in the monolayer and use the measured θ eff in the overlapping region, we obtain the map of the twist angle shown in Fig. 6f and its corresponding armchair histogram shown in Fig. 6h. In this new histogram <δ> = 30.12° and σ = 3.44°.
We applied the above analysis in the simplest case of an artificially prepared WS 2 homobilayer (see Methods and Video S3), presented in Fig. 7. Our technique readily mapped the twist angle in the overlapping region of the two SLs (see Fig. 7h). Figure 7a shows the total SHG intensity from a large area of two WS 2 SLs overlapping in a region. Focusing on a smaller ROI, that contains both overlapping and monolayer regions, we obtain in Fig. 7b the corresponding SHG within the region defined by the white rectangle. Choosing three points of interest (POIs) within the ROI, we present in Fig. 7c-f the polar diagrams of the SHG intensity as function of the polarization angle of the fundamental field. Repeating the same procedure for all the pixels inside the ROI we obtain in Fig. 7g the armchair angle mapping and in Fig. 7h the corresponding armchair angle distribution within the ROI. Finally, using the procedure described above we deduce the twist angle mapping and its corresponding distribution, shown in Figs. 7i,j, respectively. we obtain the green polar diagram that corresponds to armchair θ 2 = 12.8°; (d) P-SHG from POI3 (red dots) and fitted polar diagram (red curve) corresponding to armchair θ eff = 25.1°. This is the result of the interference between the two SLs shown in the POIs of (b) and (c). By using the measured θ eff = 25.1° and by fixing θ 1 = 39.5°, we calculate θ 2 = 10.7° (using θ 2 = 2θ eff − θ 1 ), which corresponds to the armchair of a second SL that produces the SHG interference (see Video S2). www.nature.com/scientificreports www.nature.com/scientificreports/ In particular, the measured armchair angle of SL1 in the non-overlapping region (POI1 in Fig. 7b) is θ 1 = 37.8° (Fig. 7c), whereas the effective armchair angle in the overlapping region (POI3 in Fig. 7b) is θ eff = 30.4° (Fig. 7e). The overlapping region in the ROI shown in Fig. 7a, provided the twist angle mapping shown in Fig. 7h. This www.nature.com/scientificreports www.nature.com/scientificreports/ results in θ 2 = 2θ eff − θ 1 = 23.0° (Fig. 7f). This value is indeed very close to the experimentally measured one in the non-overlapping region (POI2 in Fig. 7b) θ 2 = 23.1 o (Fig. 7d). We can finally extract the mean of the distribution of the twist angle values shown in Fig. 7h, using <δ> = 2(θ 1 − θ eff ) = 14.9 o , σ = 1.3° (Fig. 7j).

conclusions
In conclusion, we have used P-SHG imaging microscopy to map the twist angle in stacked TMD layers. In particular, the effect of the SHG interference in the overlapping region of two individual WS 2 SLs was described in terms of the newly introduced concept of the effective orientation θ eff = (θ 1 + θ 2 )/2 which dictates the P-SHG modulation in an overlapping region of two stacked twisted 2D layers. This novel concept of the effective orientation allowed the determination of the armchair orientation of a second layer that contributed to the total SHG signal detected in the overlapping region. Consequently, by determining the crystal orientation of the second layer that resulted to the measured SHG in the overlapping region we were able to calculate the twist angle between the two layers and for the first time create its pixel-by-pixel mapping both in CVD-grown and in artificially stacked WS 2 bilayers. Thus, we have demonstrated an all-optical technique that identifies the presence of stacked SLs, calculates and maps their twist angle. In addition, we have shown experimentally and interpreted theoretically, that when the twist angle between the two SLs is 30° one can selectively suppress the SHG from one layer, while at the same time the SHG from the other is switched-on and vice-versa. This SHG "magic" twist angle enables axial super-resolution imaging and consequently provides a quality characterization of the layered 2D-structure. We envisage that our methodology will provide a new and easy characterization tool of twisted 2D crystals towards their numerous optoelectronic applications.

custom-Built polarization-in, polarization-out SHG microscope. Our experimental apparatus
is based on a diode-pumped Yb:KGW fs oscillator (1030 nm, 70-90 fs, 76 MHz, Pharos-SP, Light Conversion, Vilnius, Lithuania) inserted in a modified Axio Observer Z1 (Carl Zeiss, Jena, Germany) inverted microscope (Fig. 2). The laser beam is passing through a zero-order half-wave retardation plate (QWPO-1030-10-2, CVI Laser), placed at a motorized rotation stage (M-060.DG, Physik Instrumente, Karlsruhe, Germany) that rotates with high accuracy (1°) the orientation of the excitation linear polarization. Raster-scanning of the beam at the sample plane is performed using a pair of silver-coated galvanometric mirrors (6215 H, Cambridge Technology, Bedford, MA, USA). The beam is reflected on a silver-coated mirror, at 45° (PFR10-P01, ThorLabs, Newton, NJ, USA), placed at the motorized turret box of the microscope, just below the objective (Plan-Apochromat 40×/1.3NA, Carl Zeiss). The choice of the silver coating of all the mirrors (PF 10-03-P01, ThorLabs), including the galvanometric mirrors, makes our setup insensitive to the laser beam polarization and its angle of incidence. The mean polarization extinction ratio of the different linear polarization orientations, calculated using crossed polarization measurements at the sample plane, was 28:1. In the forward direction, the SHG is collected using a high numerical aperture (NA) condenser lens (achromatic-aplanatic, 1.4NA, Carl Zeiss). The SHG is separated from the laser using a short-pass filter (FF01-720/SP, Semrock, Rochester, NY, USA) and from any unwanted signal using a bandpass filter (FF01-514/3, Semrock). A rotating film polarizer (LPVIS100-MP, ThorLabs) is placed just in front of the PMT (H9305-04, Hamamatsu, Hamamatsu City, Japan) to measure the anisotropy due to the polarization of the SHG signals. The P-SHG imaging (i.e. in the exfoliated WS 2 /WS 2 stacked structure) is performed in the epi-detection using a dichroic mirror (DMSP805R, ThorLabs), with the forward P-SHG detection module described above, placed in an epi-detection port of the microscope. Coordination of PMT recordings with the galvanometric mirrors for the image formation, as well as the movements of all the motors, is carried out using LabView (National Instruments, Austin TX, USA) software.
Samples. The WS 2 samples were grown by the low-pressure chemical vapor deposition method (LP-CVD) on a c-cut (0001) sapphire substrate (2D Semiconductors). Note that in the case of CVD-grown samples the stacking of layers is not artificial like in 12 but occurs naturally during the growth 35 . Nevertheless, we expect a similar behavior like in 35 from CVD-grown layered samples. This effect has been observed previously and is attributed to the nucleation that occurs during the CVD-growth and commences from the center 36 . WS 2 bulk crystals were exfoliated by micromechanical cleavage on a polydimethylsiloxane (PDMS) stamp, placed on top of a glass slide for optical inspection. The first monolayer was transferred on a Si/SiO 2 (285 nm) substrate by an all-dry viscoelastic stamping and then it was mounted on a XYZ micromechanical stage 37 . The stage was placed under a coaxially illuminated microscope and following the same procedure, a second WS 2 monolayer on a different PDMS was carefully aligned and then stamped slowly on top of the first monolayer. The final step included a controlled release of the PDMS stamp for the fabrication of the WS 2 bilayers.

Data availability
The data that support this study are available from the corresponding author upon reasonable request.