The optical response of artificially twisted MoS$_2$ bilayers

Two-dimensional layered materials offer the possibility to create artificial vertically stacked structures possessing an additional degree of freedom - $the$ $interlayer$ $twist$. We present a comprehensive optical study of artificially stacked bilayers (BLs) MoS$_2$ encapsulated in hexagonal BN with interlayer twist angle ranging from 0 to 60 degrees using Raman scattering and photoluminescence spectroscopies. It is found that the strength of the interlayer coupling in the studied BLs can be estimated using the energy dependence of indirect emission versus the A$_\textrm{1g}$-E$_\textrm{2g}^1$ energy separation. Due to the hybridization of electronic states in the valence band, the emission line related to the interlayer exciton is apparent in both the natural (2H) and artificial (62$^\circ$) MoS$_2$ BLs, while it is absent in the structures with other twist angles. The interlayer coupling energy is estimated to be of about 50 meV. The effect of temperature on energies and intensities of the direct and indirect emission lines in MoS$_2$ bilayers is also quantified.

Two-dimensional layered materials offer the possibility to create artificial vertically stacked structures possessing an additional degree of freedomthe interlayer twist. We present a comprehensive optical study of artificially stacked bilayers (BLs) MoS 2 encapsulated in hexagonal BN with interlayer twist angle ranging from 0 to 60 degrees using Raman scattering and photoluminescence spectroscopies. It is found that the strength of the interlayer coupling in the studied BLs can be estimated using the energy dependence of indirect emission versus the A1g-E 1 2g energy separation. Due to the hybridization of electronic states in the valence band, the emission line related to the interlayer exciton is apparent in both the natural (2H) and artificial (62 • ) MoS 2 BLs, while it is absent in the structures with other twist angles. The interlayer coupling energy is estimated to be of about 50 meV. The effect of temperature on energies and intensities of the direct and indirect emission lines in MoS 2 bilayers is also quantified.

I. INTRODUCTION
Two-dimensional (2D) van der Waals (vdW) crystals have emerged as a new generation of materials with extraordinary properties. For instance, widely studied semiconducting transition metal dichalcogenides (S-TMDs) transform from indirect-to direct-band gap, optically-bright semiconductors when thinned down to a monolayer (ML), which results in unique electronic structures and consequent optical properties [1][2][3][4] . The family of 2D layered materials grows day by day, hugely expanding the scope of possible phenomena to be explored in two dimensions. Growing is also the number of possible vdW heterostructures that one can create. Such 2D materials currently cover a vast range of properties allowing potential applications in, i.a. spintronic devices 5 , optoelectronics 6,7 , tunnel field-effect transistors 8,9 , single-photon sources 10,11 , and quantum information processing 12,13 . Rapid advances in fabrication methods, like chemical vapor deposition (CVD) growth and mechanical exfoliation techniques, have contributed to increased interest in artificial stacking of different layered materials on top of each other. The simplistic approach of producing vertical vdW heterostructures without the constraints of crystal lattice mismatch enables integrating various 2D materials to create diverse systems with new electronic properties that are not present in pristine components. In addition to the selection of compounds in terms of their properties, a new degree of freedom has emerged: the twist angle between stacked layers, which gives rise to the group of the so-called twistronic materials 14,15 . The twist angle is responsible for the occurrence of moiré patterns, that leads to new and intriguing phenomena, like the forma-tion of secondary Dirac points in graphene on hexagonal boron nitride (hBN) 16,17 or hybridized (moiré) excitons in vdW heterostructures formed by stacked two S-TMD MLs [18][19][20][21][22][23] .
In this work, we investigate interlayer interactions in high-quality artificially stacked twisted MoS 2 bilayers (BLs) encapsulated in hBN flakes using the RS and PL spectroscopies. Our results indicate that the interlayer coupling can be determined by the comparison of the emission energies of the indirect transition versus the energy separation between two basic intralayer phonon modes (A 1g and E 1 2g ). The origin of the apparent emission due to the interlayer excitons in both the natural (2H) and artificial (62 • ) MoS 2 BLs and its absence for the BLs with other twist angles is associated with hybridization of electronic states in the valence band. We also investigate the evolution of energies and intensities of the emission lines due to direct and indirect transitions with temperature ranging from 5 K to 300 K. cal formula MX 2 where M=Mo or W and X=S, Se or Te, which the most common crystallographic structure is a hexagonal phase. In that case, X-Mo-X atoms in a monolayer (1 L) of MoX 2 are arranged in a trigonal prismatic structure, which does not exhibit inversion symmetry. A BL formed by the stacking of two MLs exhibits an additional degree of freedom -the twist angle between the layers. Due to the hexagonal symmetry of a MX 2 ML, the different arrangements are characterized by twist angles ranging from 0 • to 60 • . The schematic illustration of three patterns of MoS 2 BLs with the twist angles of 60 • , 15 • , and 0 • are presented in Fig. 1. The most stable structures, which exhibit the strongest coupling, are found for twist angles equal to 0 • and 60 • ascribed correspondingly to the 2H and 3R stackings 32,33 . We do not use the label 3R for the 0 • BL (see Fig. 1), because the 3R unit cell involves the atoms from three consecutive layers (2H unit cell is composed of two layers). While both the 2H and 3R polytypes are accessible in natural and CVD-grown MoS 2 multilayers 33 , the BLs with other twist angles can only be constructed by the artificial stacking of individual monolayers. It is important to point out that the inversion symmetry is restored in the 2H BL, while the zero-twist angle BL has neither inversion nor in-plane mirror.

B. Raman scattering spectroscopy
In order to investigate the interlayer coupling in the twisted BLs and therefore to distinguish their pristine properties, we compare their RS and PL spectra with those obtained for 1 L and 2 L exfoliated from 2H bulk MoS 2 . Note that the twist angle between MLs forming BLs was determined using the second-harmonic generation (SHG) technique, which permits to reveal the crystallographic orientation of individual MLs 34 . The measured RS spectra of the selected BLs with twist angles of 62 • , 18 • , and 6 • accompanied with the ones of 1 L and 2 L MoS 2 are presented in Fig. 2 , as a function of twist angle (θ) for all the studied twisted homobilayers. The corresponding ∆ω for 1 L and 2 Ls, represented by solid symbols, are also shown.
First, we focus on the low-energy range of the RS spectra, presented in Fig. 2. There are no modes related to rigid vibrations in 1 L affirming their interlayer nature. For 2 L, two modes are observed at ∼23 cm −1 and ∼44 cm −1 , which are associated with the in-plane shear (SM) and out-of-plane breathing (BM) vibrations of rigid layers, respectively 33,35,36 . In twisted BLs with angles of 62 • and 6 • , the SM and BM modes are also apparent. This indicates that the coupling between adjacent layers in these two BLs is relatively strong due to their resemblance to the natural 2H and 3R stackings. Note that the relatively smaller intensity of the SM peak in the 6 • BLs as compared with the 2H and 62 • ones can be attributed to the smaller interlayer bond polarizability for the 3R polytype 37 . This effect was applied to identify the stacking order of S-TMD layers 30,37 . On the other hand, only the BM-related peak can be observed for the sample with the 18 • angle as well as for other BLs with angles in the range from ∼10 • to ∼56 • (data not shown). The higher energy range of all RS spectra shown in Fig. 2(b) is dominated by two phonon modes ascribed to the in-plane (E 1 2g ) and out-of-plane (A 1g ) intralayer vibrations. In few layers 2H MoS 2 , the E 1 2g (A 1g ) phonon mode experiences a red (blue) shift with increasing the layer thickness. The energy separation between those modes, i.e. ∆ω = ω A1g − ω E 1 2g , is a useful tool to determine the number of MoS 2 layers 24 . As can be seen in the Figure, the ∆ω energy difference varies for different twist angles (Θ), which is summarized for all samples in Fig. 2  . The corresponding emission energies for the 2 L are also indicated with solid orange points.
to natural 2H and 3R BLs 38 . In the intermediate cases (6 • <Θ<60 • ), the ∆ω falls lower than for 2 L MoS 2 , but it is considerably larger than in 1 L. The ∆ω can be also used to characterize the effective interlayer mechanical coupling strength or, in other words, the distance between layers in twisted BLs 32 . In such case, we can conclude that the structures with the twist angle equal to 6 • and 62 • are characterized by the strongest coupling (the smallest interlayer distance), while the coupling strength is weaker (the interlayer distance is bigger) for other BLs with the twist angles substantially different from 0 and 60 degrees.
Note that our RS results are consistent with previous works devoted to the twisted MoS 2 BLs exfoliated on Si/SiO 2 substrates 32,39-42 .

C. Photoluminescence spectroscopy
Previous investigation of PL spectra of artificially twisted BLs MoS 2 were carried out at room temperature 32,39,42 . The PL spectroscopy at low temperature may provide however a more accurate analysis, which is due to the much lower emission linewidths 4,43 . The PL spectra of the selected BLs with twist angles of 62 • , 18 • , and 6 • accompanied with those of 1 L and 2 L MoS 2 mea-sured at low temperature (T =5 K) are shown in Fig. 3(a). Similarly to the aforementioned analysis of the RS, we begin with an examination of the PL spectra of the natural 1 L and 2 L MoS 2 . The 1 L spectrum consists of two narrow emission lines apparent in the vicinity of the optical band gap (so-called A exciton), which can be ascribed to the neutral and charged excitons in accordance with previous reports [44][45][46] . In contrast, the PL spectrum of the 2 L MoS 2 is composed of two distinct emission bands: (i) the transitions in the vicinity of the direct A (X A ) and B (X B ) excitons formed at the K ± points of the Brillouin zone (BZ), which are observed correspondingly at about 1.93 eV and 2.07 eV; (ii) the significantly much intense transitions, denoted as I and I', apparent at about 1.5 eV, which are ascribed correspondingly to an indirect recombination process between the Λ and K point in the conduction band (CB) and Γ points in the valence band (VB) of the BZ 32,43,47 . The PL spectra of the artificial twisted BLs also comprise two direct-and indirectrelated bands. While the energies of direct transitions are hardly affected by the twist angles, the energies of the indirect ones change significantly by about 150 meV. The low-intensity emission bands, denoted with * and apparent between the aforementioned I and X A lines, can be described as a remaining part of the defect-related emission, which is reported commonly for 1-L MoS 2 exfoliated on Si/SiO 2 substrates [44][45][46] .
To appreciate the effect of the twist angle on the PL spectra, the energy evolution of the X A , X B and I transitions for all the studied samples as a function of twist angle is presented in Fig. 3(b). As can be seen in the Figure, there is no effect of the twist angle on the X A and X B energies. This reflects a pure two-dimensional character of both the A and B excitons, i.e. these complexes are distributed spatially within a single layer even in a bulk form of MoS 2 48-50 . For the indirect transitions, we focus only on the I emission lines as they are observed for all studied samples. The emission reaches the lowest energies for border cases (6 • and 62 • ), which are similar to the value obtained for the natural 2 L MoS 2 . For BLs with other twist angles, the I energies are substantially higher by about 150 meV and they are almost independent of twist angle. This effect can be associated with changes in the interlayer distance between MLs forming the studied BLs as a function of twist angle, as it was reported previously in Refs. 32,39. As the energy separation between the A 1g and E 1 2g peaks (∆ω) can be considered as a probe of the interlayer distance in the studied artificial MoS 2 BLs, one can plot the energy dependence of the indirect transition I as a function of ∆ω in Fig. 3(c). The presented results can be divided into two groups: (i) for structures with twist angle close to 0 • and 60 • , i.e. ∆ω ∼23 cm −1 , the I emission energy equals approx. 1.5 eV; (ii) for samples with twist angle significantly different from 0 • and 60 • , i.e. ∆ω ∼22 cm −1 , the energy dispersion of the I line is of about 80 meV centered at ∼1.58 eV. These results are consistent with previously obtained for stacked triangle- shaped MoS 2 BLs grown on Si/SiO 2 substrate using CVD technique 32 . Summarizing, we can assume that the analysis of the indirect emission versus the ∆ω can be used as a useful tool to probe the interlayer distance in the twisted homo-BLs. This allows to distinguish between two cases of twist angles, which are in the vicinity of 0 • and 60 • or the twist angle is in-between.
While the X A and X B lines are observed for all studied samples, there is an additional emission line, labelled IL, apparent between the former ones for the 2 L and 62 • samples, which emission was not reported so far. According to recent results of reflectance contrast (RC) experiments carried out on natural 2H MoS 2 BLs at low temperature [51][52][53] , that line can be associated with the recombination of the so-called interlayer A exciton. The IL line originates from the hybridization of electronic states in the VB of the 2H-stacked BL due to the interlayer coupling (see Ref. 51 for details). The strength of the related VB coupling is described by the interlayer hopping term, t ⊥ . Based on a kp model of 2H BLs in the vicinity of K ± points 51 , the t ⊥ parameter is given by: where ∆ 2H A-B is the A-B energy difference for natural 2H BL and ∆ 2H v represents the spin-orbit VB splitting without interlayer coupling. Note that the detailed investigation of the band structure of the 2H and 0 • BLs in the vicinity of the K ± point of the BZ is performed in Appendices B and A. According to that analysis, we assume that: (i) the ∆ 2H v is an order of magnitude larger than its counterpart in the CB ∆ 2H c 54 ; (ii) the binding energies of the A and B excitons are comparable; (iii) the ∆ 2H v can be roughly approximated by the A-B energy difference for the BL with twist angle of 0 • . Consequently, the diagram of relevant subbands at the K + point of the BZ  in the MoS 2 BLs with twist angle of 60 • (2H) and 0 • is shown in Fig. 4. As can be appreciated from the Figure, the t ⊥ parameter can be evaluated using the energy separation between the A and B excitons measured in the 60 • (2H) and 0 • BLs. The extracted ∆ 2H A-B ∼180 meV (the corresponding ∆ 62 • A-B is about 10 meV smaller, which suggests the lower strength of interlayer coupling in twisted BL) and ∆ 6 • A-B ∼150 meV are in very good agreement with the reported ones for natural BLs 51-53 . Using these values, the extracted interlayer hopping parameter (t ⊥ ) is found to be on the order of 50 meV, which agrees very well with the value obtained using the RC experiment performed on the MoS 2 BLs grown using CVD technique and encapsulated in hBN flakes (49 meV) 53 .
The energy evolution of both the direct (X A , X B , and IL) and indirect (I and I') lines as a function of temperature measured on BL (2 L) and selected homobilayers with twist angles of 6 • , 18 • and 62 • is shown in Fig. 5. As can be appreciated in the Figure, all lines within the direct-and indirect-related transitions located at ∼2.0 eV and ∼1.5 eV are characterized by the same type of evolution upon increasing temperature, i.e. redshift and blueshift, respectively. The X A , X B , and IL peaks redshift when temperature is increased from 5 K to 300 K. This can be associated with the re-duction of the direct band gap resulting from the temperature expansion of those layers in lateral directions. Consequently, these evolutions can described by the relation proposed by O'Donnell et al. 55 , which expresses the temperature dependence of the band gap in terms of the average energy of acoustic phonons involved in the electron-phonon interaction ω . The relation reads , where E 0 stands for the band gap at absolute zero temperature, S is the coupling constant, and k B denotes the Boltzmann constant. We found that ω stayed on nearly the same level, ∼ 21 meV, for the X A , X B , and IL transitions in 2H BL 47 . As a consequence, we kept it fixed during fittings of all the experimental data shown in Fig. 5. It can be seen that the fitted curves correctly reproduce the energies of excitonic lines (see Fig. 5), which suggests that binding energies of investigated complexes do not depend on temperature. Interestingly, the overall redshift of X A lines are of about 60-70 meV, while the corresponding shift for the MoS 2 BL exfoliated on Si/SiO 2 substrate was found to be of almost 90 meV 43 . Simultaneously, the crystal expansion across the layers leads to a larger separation between the layers, which results in the blueshift of the indirect band gap. The I emission lines experience almost monotonic blueshifts in the range of about 20-40 meV, while the shift of ∼80 meV is reported in Ref. 43. Both these results indicate that the hBN encapsulation of the BLs induces much smaller expansion of the flake in both directions (in lateral directions and across the layers) with increasing temperature as compared with the BLs exfoliated on Si/SiO 2 substrate 43 .
The temperature evolution of the integrated intensities of the direct-and indirect-related transitions is presented in Fig. 6. It can be seen that the total intensities of all direct transitions mostly decrease monotonically with increasing temperature. The observed trends are in agreement with previous reports for direct emissions apparent in MoS 2 ML 47,56 and probably result from the competition between the efficiencies of the radiative and non-radiative recombination channels. Surprisingly, the temperature dependencies of the corresponding indirect transition intensities are clearly non-monotonic. First, quick reductions of the related intensities is observed from 5 K to about 100 K-150 K, which can be associated with increase of kinetic energies of excitons and interplay between radiative and non-radiative processes. Then it is followed by slower increase of the intensities up to room temperature, which can be explained in terms of the increased population of phonons at higher temperatures.

III. SUMMARY
A systematic investigation of optical properties of hBN-encapsulated artificially stacked MoS 2 BLs with varying interlayer twist angle has been conducted. It has been shown that RS spectroscopy is a sensible tool for determining the interlayer coupling and spacing in the BL. The strongest coupling has been observed in MoS 2 homobilayers with twist angles of 6 • and 62 • , which reproduce well the structure of 2H and 3R polytypes found in natural and CVD-grown MoS 2 multilayers. Increasing the twist angle of homobilayers from 6 • to ∼ 30 • leads to an increase in the interlayer spacing, and thus a decrease in the interlayer coupling. This effect can be observed as the absence of the low-energy interlayer phonon modes in the RS spectrum as well as the energy shift of ∼ 140 meV for the indirect transition between K-Γ points of the Brillouin zone. The PL of artificially assembled homobilayers, in comparison with the natural 2 L MoS 2 , features a single line associated with indirect emission, but its energy varies considerably depending on the interlayer twist. We conclude that gaining a new degree of freedom by inter-layer twisting of artificially assembled flakes permits the control over the energy of the indirect transition, This offers further insight in few-layer 2D systems and can be a useful tool for future device applications.

METHODS
The investigated MoS 2 BLs, MLs, and hBN flakes were fabricated by two-stage PDMS-based mechanical exfoliation of bulk crystals. An unoxidized Si wafer was used as a substrate. In order to ensure the best quality of the substrate surface, they were annealed at 200 • C and kept on hot-plate until the first non-deterministic transfer of h-BN flakes. Subsequent layers were transferred deterministically, to reduce inhomogeneities between each transfer sample were annealed. The complete structures were annealed at 160 • C for 1.5 hours to ensure the best layerto-layer and layer-to-substrate adhesion and to eliminate a substantial portion of air pockets on the interfaces between the constituent layers.
The SHG measurements performed to identify the relative twist angle of the exfoliated S-TMD MLs forming homobilayers were taken at T = 300 K using a homebuilt setup with a femtosecond Ti:Sapphire laser with excitation at 800 nm (1.55 eV). For each measurement, the laser light had a typical incident power of 500 µW, was linearly polarized, and focused to a spot size of 1 µm by a 50x objective lens. A set of a motorized half-wave plate and a fixed linear polarizer were used to analyse an SHG signal, which was detected by a Si avalanche photodiode.
The PL and RS measurements were performed using λ=515 nm (2.41 eV) radiation from continuous-wave Arion or diode lasers. For the PL and RS experiments at room temperature (T = 300 K), the excitation light was focused through a 100x long-working distance objective with a 0.55 numerical aperture (NA) producing a spot of about 1 µm diameter. The signal was collected via the same microscope objective, sent through a 1 m monochromator, and then detected by using a liquid nitrogen-cooled charge-coupled device (CCD) camera. To detect low-energy RS up to about ±10 cm −1 from the laser line, a set of Bragg filters was implemented in both excitation and detection paths. The temperaturedependent PL measurements were performed using an analogous setup with small modifications (a 50x longworking distance objective and a 0.5 m monochromator). Moreover, the studied samples were placed on a cold finger in a continuous flow cryostat mounted on x-y motorized positioners. The excitation power focused on the sample was kept at 200 µW during all measurements to avoid local heating. We model a S-TMD bilayer with 60 • -angle alignment as a pile of two monolayers (top and bottom), placed in parallel to xy plane. We define the positions of metal and chalcogen atoms of the bottom layer as respectively. Here we introduced the in-plane primitive lattice vectors of length a 0 the pair of integer numbers (m, n), the short notation for the (m, n)-th lattice vector R mn = a 1 m + a 2 n and unit vectors of Cartesian coordinate system e x , e y , e z . Vectors define in-plane positions of the metal and chalcogen atoms within a unit cell of S-TMD monolayer, respectively. Vectors ±ηe z with η > 0 define the out-of-plane positions of chalcogen atoms in the unit cell. We also introduce the primitive vectors of reciprocal lattice They satisfy the orthogonality property a j b k = 2πδ jk , where δ jk is the Cronecker delta. The top lattice of the bilayer can be obtained from the bottom one as a result of the shift along z direction on some distance l with subsequent rotation around Oz axis on 180 • degree. Then, the positions of the metal and chalcogen atoms of the top lattice become Note that all of the chalcogen atoms of bottom layer are aligned with the metal atoms of the top layer (and vice versa) along z-direction. Hence, the unit cell of the bilayer contains twice more atoms than in 1 L. The positions of metal and chalcogen atoms within the unit cell are defined by vectors {t M , t X + le z } and {t X ±ηe z , t M ±ηe z +le z }, respectively. This arrangement of atoms is called 2H-stacking and corresponds to thermodynamically stable form of all S-TMD crystals with any number of layers including bulk.
The bilayer possesses C 3 rotation symmetry (with Oz line as the rotational axis) and mirror symmetry P : x ↔ −x (the mirror's plane is yz-plane). Therefore, the crystal has the same hexagonal Brillouin zone as the Brillouin zone of the bottom layer. Hence, it is convenient to use the known Bloch states of the monolayer as the basis states. Namely we are interested in conduction (c) and valence (v ) band states at the K ± points. To this end, we introduce the vectors ±K = ±(b 1 + b 2 )/3, which define the position of the K ± points in the reciprocal space, respectively. In further we will use the notation ±K both for vectors and the positions of the edges of the Brillouin zone (K ± ), for clarity. In the vicinity of ±K points the Bloch states are predominantly made of the d-orbitals of metal atoms. The corresponding valence and conduction bands states of the bottom layer can be presented as Here N is the normalization factor and Y lm (r − R) is the value of the lm-th atomic orbital placed at the point R and calculated at the point r. The operator C 3 , which generates R 2π/3 rotation of the vectors in space, transforms the corresponding Bloch functions as The phases, which appear under the transformation of the basis states of the bottom layer correspond to the notation in Ref. 57.
In addition, the crystal has inversion symmetry I : r ↔ −r + 2R I , with the center of inversion in the point R I = le z /2. This symmetry together with the time-reversal symmetry induces the restriction on the band structure of the crystal. In accordance to the Kramers theorem, all the bands of the crystal become doubly degenerated by spin. Therefore, it is convenient to define the second pair of basis states, associated with the top layer, using the above-mentioned symmetry operations Ψ t ±K,n (r) = K 0 IΨ b ±K,n (r). Here K 0 and I are complex conjugation and inversion symmetry operators, respectively. Using the tight-binding representation of the basis states of the bottom layer we get The states satisfy the following transformation rules un- ±K,c (r), with phases which are opposite to the phases of the basis states of the bottom layer . It leads to the fact that bilayer crystal can absorb the light with both circular polarizations in K point as well as in −K one. This feature is a consequence of inversion symmetry of the crystal. Finally, the mirror symmetry operator P acts on the basis states as P Ψ α ±K,n (r) = [Ψ α ±K,n (r)] * = Ψ α ∓K,n (r), where n = c, v and α = b, t.
In further we focus on the states at the K point for brevity. We take into account the spin degree of freedom s =↑, ↓ of electron excitations and introduce the following set of 8 basis states According to the kp method, developed for S-TMD multilayers 4,51,58 the quasiparticles with the momentum k = k x e x + k y e y at the K point are described by the matrix elements Ψ α n , s| H|Ψ α n , s of the one-particle Hamiltonian Here m 0 is electron's mass, c -speed of light, -Planck's constant, σ = (σ x , σ y , σ z ) are Pauli matrices and p = −i ∇ is the momentum operator. The first term of the Hamiltonian defines the kinetic energy of an electron which propagates in the crystal field of the bottom U b (r) and top U t (r) layers of the bilayer. The next two terms describe the spin-orbital interaction in the system, induced by the potentials U b (r) and U t (r), respectively. The last kp term couples valence and conduction bands. This coupling is supposed to be small and we omit its effects for the current study. The detailed analysis of the impact of the kp term can be found in Refs. 51,58.
We consider first the matrix elements of the states of bottom layer. We present the Hamiltonian as where we introduced the Hamiltonian of the bottom and the term which affects the motion of quasiparticles of the bottom layer by the crystal field of the top layer The Hamiltonian H b 0 (r) has the diagonal matrix elements,which are nothing but the position of the conduction and valence bands in monolayer Here E v and E c are positions of the valence and conduction bands without spin splitting, ∆ v and ∆ c are their spin splittings, and σ s = +1(−1) for s =↑ (↓) states respectively. Note that in K point ∆ v is always positive, while ∆ c can be negative (bright type of S-TMD) and positive (darkish type of S-TMD).
The dominant contribution from of H t int (r) is the diagonal matrix elements in spin space We suppose that the corrections to the splitting are small |∆ c | |δ∆ c |, |∆ v | |δ∆ v |, and the type of S-TMD bilayer remains the same as the type of its constituents. H t int (r) term has also non-zero matrix element between valence and conduction bands of the same layer and opposite spins (see Refs. 59,60 for details). However, these matrix elements give the negligibly small contribution to the energies of the bands, proportional to ∼ 1/(E c − E v ). Therefore, we omit them from the study.
The matrix elements between the states of the top layer can be calculated in the same way. We present the total Hamiltonian as where the first term is the Hamiltonian of the top monolayer, while the second term affects the motion of quasiparticles of the top layer by the crystal field of the bottom one. With the help of Kramers theorem and the latter result one can immediately get the answer for the matrix elements of the above-mentioned terms Note that the sign before spin-splitting terms for the states of the top layer is opposite to the sign of the same terms of the bottom layer. This is the manifestation of the double degeneracy by spin of all the bands of the bilayer. Finally, we calculate the matrix elements of the Hamiltonian between the states of the different layers. Namely, we evaluate the following interlayer matrix elements Ψ t n , s|H(r)|Ψ b n , s , where n, n = c, v. The other matrix elements can be obtained by complex conjugation of the considered ones. We present the Hamiltonian in the following way and suppose the orthogonality of the states from the opposite layers Ψ t n , s|Ψ b n , s = 0. Then The p 2 operator is a spin singlet, hence the matrix elements are diagonal in spin space Ψ t n , s| p 2 |Ψ b n , s = δ ss Ψ t n | p 2 |Ψ b n , where |Ψ α n = Ψ α K,n (r). Using transformation properties of the basis states under C 3 rotation we get that only one matrix element is non-zero. The P symmetry of the crystal dictates that the parameter t ⊥ is a real number Im Here we introduced he absolute values of the binding energies of corresponding excitons E A , E B , E A , E B , and we the short notations E 2H The sketch of the bands position in MoS 2 bilayer with 2H-stacking is presented in Fig. 7.
For the particular case of MoS 2 the splitting in conduction band is supposed to be much smaller than the splitting in valence band |∆ 2H c | |∆ 2H v |, and the binding energies of A and B excitons are considered to be equal E A = E B . In this approximation we have the following result ∆ 2H Note that the intralayer and interlayer optical transitions in the same K (or −K) point are characterized by opposite circular polarizations and g-factors of corresponding excitons (see Refs. 51,58 for details). In order to compare the results of the measurements presented in the main text we describe the optical properties of the bilayer S-TMD with zero-angle alignment (in further "bilayer") in the similar way as it was done for 2H-stacked bilayer. Again, we consider the bilayer as a pile of two monolayers (top and bottom), placed in parallel to xy plane. We assume the positions of metal and chalcogen atoms of the bottom layer are the same as in the previous section The top lattice of the bilayer can be obtained from the bottom one as a result of two consequent shifts: along z direction on distance l (which is not equal to the distance l for the 2H-stacked bilayer) and then along y direction on distance a 0 / √ 3. Then, the position of the metal and chalcogen atoms of the top lattice can be presented as Note that the half of the chalcogen and half of metal atoms in this bilayer are aligned in z-direction. This type of stacking for hexagonal lattices is called Bernal or AB-stacking. The unit cell of the considering bilayer contains twice more atoms than in monolayer. The positions of metal and chalcogen atoms within the unit cell are defined by vectors {t M , a 1 +le z } and {t X ±ηe z , t M ± ηe z + le z }, respectively.
Note that the considering lattice has neither in-plane mirror symmetry (like AA-stacked case) nor inversion symmetry (like 2H-stacked bilayer). It possesses only C 3 rotation symmetry (with Oz line as a rotational axis) and mirror symmetry P : x ↔ −x (the mirror's plane is yz-plane). Again the crystal has the same hexagonal Brillouin zone as the Brillouin zone of the bottom layer. Hence, we choose the same basis states for valence and conduction bands in the ±K points of the bottom layer as we have in the previous section We define the Bloch states of the top lattice in the same They have the corresponding transformation rules ±K,c (r). After rotation the basis states of the top layer get the phases which deviates from the corresponding phases of the basis states of the bottom layer. This difference is the result of the shift in y direction of the atoms of the top lattice with respect to the bottom one. Despite this difference, the optical transition rules are not changed for each separate layer of bilayer. Namely the top and bottom layers in ±K points absorb the σ ± circular polarized light, respectively. These optical properties can be understood as a consequence of time-reversal symmetry which couples K and −K points of bilayer. This feature demonstrates the significant difference in optical properties of 0 • -and 60 • -aligned S-TMD bilayers. The mirror symmetry transformation P also couples ±K points P Ψ α ±K,n (r) = [Ψ α ±K,n (r)] * = Ψ α ∓K,n (r), where n = c, v and α = b, t.
Like in previous section we focus on the states at the K point, take into account the spin degrees of freedom and introduce the basis states |Ψ t c , s = Ψ t K,c (r)|s , |Ψ t v , s = Ψ t K,v (r)|s .
All the terms in this Hamiltonian has the same meaning as in previous section. We also omit the kp term, since we are focused on the optical transitions exactly in ±K points for clarity. Let us consider the states of bottom layer. We present the Hamiltonian in the following form where we introduced the Hamiltonian of the bottom monolayer H b 0 (r) = p 2 2m 0 + U b (r) + 4m 2 0 c 2 ∇U b (r), p σ, (B13) and the term which affects the motion of quasiparticles of the bottom layer by the crystal field of the top layer H t int (r) = U t (r) + 4m 2 0 c 2 ∇U t (r), p σ. (B14) The Hamiltonian H b 0 (r) has the diagonal matrix elements where all the parameters have the same meaning as in the previous section. The dominating contribution from of H t int (r) has also diagonal matrix elements Ψ b c , s|H t int (r)|Ψ b c , s = δE t c + σ s δ∆ t c /2, This term also has non-zero matrix element between valence and conduction bands of opposite spins. As in previous case they give a negligibly small contribution to the energy of the bands, proportional to ∼ 1/(E c − E v ). Therefore, we omit these matrix element from the current study.
The matrix elements between the states of the top layer can be calculated in the same way. Namely we present the Hamiltonian as with the Hamiltonian of the top monolayer H t 0 (r) = p 2 2m 0 + U t (r) + 4m 2 0 c 2 ∇U t (r), p σ, (B20) and the term which affects the motion of quasiparticles of the top layer by the crystal field of the bottom one Again, only diagonal matrix elements of the H t 0 (r) are nonzero Ψ t c , s|H t 0 (r)|Ψ t c , s = E c + σ s ∆ c /2, (B23) The dominating contribution from of H b int (r) has also diagonal matrix elements