Fine structure of $\mathrm{K}$-excitons in multilayers of transition metal dichalcogenides

Reflectance and magneto-reflectance experiments together with theoretical modelling based on the $\mathbf{k\cdot p}$ approach have been employed to study the evolution of direct bandgap excitons in MoS$_2$ layers with a thickness ranging from mono- to trilayer. The extra excitonic resonances observed in MoS$_2$ multilayers emerge as a result of the hybridization of Bloch states of each sub-layer due to the interlayer coupling. The properties of such excitons in bi- and trilayers are classified by the symmetry of corresponding crystals. The inter- and intralayer character of the reported excitonic resonances is fingerprinted with the magneto-optical measurements: the excitonic $g$-factors of opposite sign and of different amplitude are revealed for these two types of resonances. The parameters describing the strength of the spin-orbit interaction are estimated for bi- and trilayer MoS$_2$.

Reflectance and magneto-reflectance experiments together with theoretical modelling based on the k · p approach have been employed to study the evolution of direct bandgap excitons in MoS2 layers with a thickness ranging from mono-to trilayer. The extra excitonic resonances observed in MoS2 multilayers emerge as a result of the hybridization of Bloch states of each sub-layer due to the interlayer coupling. The properties of such excitons in bi-and trilayers are classified by the symmetry of corresponding crystals. The inter-and intralayer character of the reported excitonic resonances is fingerprinted with the magneto-optical measurements: the excitonic g-factors of opposite sign and of different amplitude are revealed for these two types of resonances. The parameters describing the strength of the spin-orbit interaction are estimated for bi-and trilayer MoS2.

I. INTRODUCTION
Scientific curiosity to uncover the properties of new materials and to demonstrate their possible novel functionalities drive the research efforts focused on atomically-thin matter, and, in particular, on thin layers of semiconducting transition metal dichalcogenides (S-TMD) [1][2][3][4]. Intense works have been devoted to studies of S-TMD monolayers which appeared to be the efficient light emitters, the two-dimensional semiconductors with a direct bandgap positioned at the K ± points of their 1-st hexagonal Brillouin zone (BZ) [5][6][7]. New and rich possibilities of tuning the band structure, the strength of Coulomb interaction, and thus the optical properties, are opened when stacking the S-TMD monolayers into a form of multilayers and/or hetero-layers [8][9][10][11][12][13][14][15][16]. The properties of the archetypes of S-TMD stacks which are the thermodynamically stable 2H-stacked multilayers are to be well understood first.
In 2H stacks of N monolayers (N ML), the electronic bands are known to be effectively modified, with N , in the range outside the K ± points of the BZ [17][18][19][20][21][22][23]. This, in particular, implies the indirect bandgap in N MLs when N > 1, what strongly affects the emission spectra of these multilayers [24][25][26][27][28][29]. Instead, more subtle effects of the hybridization of electronic states around the direct bandgap which appears at K ± points of the BZ in any N ML are relevant for the absorption-type processes [30][31][32]. Understanding the absorption response of S-TMD multilayers might be of special importance for their potential applications in photo-sensing or photo-voltaic devices [33][34][35][36]. * artur.slobodeniuk@lncmi.cnrs.fr † marek.potemski@lncmi.cnrs.fr In this paper we present the theoretical outline (based on k · p approach [17,27,30,37]) and the experimental data (results of reflectance and magnetoreflectance measurements), which, in a consistent manner, unveil the nature of direct bandgap (K ± points) excitonic transitions in 2H-stacked S-TMDs multilayers. The geometry of 2H stacking, wherein each subsequent layer is 180 • rotated around previous one, induces interaction and hybridization of the Bloch states in K + /K − /K + . . . (K − /K + /K − . . .) points of subsequent monolayers, which form the K + (K − ) valleys of the multilayer. Such non-trivial interlayer coupling delocalizes the electron states in the out-of-plane direction, which leads to the formation of the new type of excitons, associated with K ± valleys of multilayers as well as bulk crystals [27,38,39]. We classify such excitons, associate them with the symmetry of the crystal and describe their properties in terms of the simple theoretical model. Notably, the intra-or interlayer character of the excitonic transitions is fingerprinted with, correspondingly, positive or negative sign of the g-factor associated to these transitions.
Due to the different symmetries of the conduction and valence band orbitals, the hybridization of K ± -electronic states in N ML occurs predominantly in the valence band and is more effective when the spin-orbit splitting ∆ v in the valence band is small. MoS 2 crystals, with the smallest ∆ v among all other S-TMDs, have been therefore chosen for the investigations. In the experiments, we have largely profited of the significantly improved optical quality of MoS 2 layers when they are encapsulated in between the hexagonal boron nitride (hBN) [11,[40][41][42][43][44]. We consider the bi-and trilayer systems as the simplest multilayer representatives with different spatial symmetry. Comparing the experimental data with the theoretical expectation on a more quantitative level, we discuss the characteristic parameters which reflect the effects of spin-orbit interaction in our bi-and tri-layer MoS 2 . Intriguingly, we estimate/speculate that in these multilayers the spin-orbit splitting in the conduction band is quite large (∼ 50 meV) whereas the spin-orbit coupling parameter for the valence band (∼ 120...140meV) is somewhat small, as compared to the expected values for the MoS 2 monolayer.
The paper is organized as follows. The section II introduces the theoretical description of the interlayer coupling and optically active transitions in 2H stacked multilayers. In the section III we present the experimental data for excitonic resonances in bi-and trilayer of MoS 2 . In the section IV we outline the properties of multilayers and the possible applications of such materials. In the Appendix A the samples, instrumentation and experimental details are presented. The Appendices B and C contain the derivation and discussion of the exciton g-factors in bi-and trilayer of MoS 2 respectively.

II. THEORETICAL DESCRIPTION
We consider the optically active transitions at the K ± points in 2H-stacked multilayer S-TMD crystals encapsulated in hBN. Our investigation is based on k · p approximation. We focus on the optical properties of such crystals at K + point for brevity [27]. The results for K − point can be obtained by analogy.
We briefly remind the features of S-TMD monolayer which will be used in the subsequent description of the multilayer structures. Namely, a single layer crystal is a direct band gap semiconductor. The maximum of valence (VB) and minimum of conduction (CB) bands are located at the K ± points of BZ. Due to strong spin-orbit interaction both bands are spin-split (from hundreds of meV in VB, up to tens meV in CB). Hence, Bloch states in the corresponding points can be presented as a tensor product of spin and spinless states. The spinless valence |Ψ v and conduction |Ψ c band states at the K + point are made predominantly from d x 2 −y 2 + id xy and d z 2 orbitals of transition metal atoms respectively [45,46]. Due to time-reversal symmetry (TRS) the analogous states at the K − point are complex conjugated to the previous ones, i.e. they are made from d x 2 −y 2 − id xy and d z 2 orbitals accordingly. The TRS also dictates that Bloch states with the same band index (c or v) but with opposite spins in different valleys have equal energies. The crystal spatial symmetry together with TRS define the optical properties of monolayers -only σ ± polarized light can be absorbed or emitted at the K ± points respectively. Since the VB and CB are split, there are two possible optical transitions (which conserve spin) at the K ± points -the lower-energy T A and the higher-energy T B ones. All the described features are depicted on Fig. 1 for the case of MoS 2 .
A bilayer crystal can be presented as two monolayers separated by a distance l ∼ 6Å [47], and in which the top layer is rotated with respect to the bottom one by |∆c| and ∆v define the absolute values of spin-orbit splitting in the CB and VB, respectively. Eg is the single-particle band gap. 180 • . We arrange them in z = l/2 and z = −l/2 planes, respectively (see Fig. 2). In this presentation, the crystal possesses the centrosymmetry I with an inversion center in the z = 0 plane. The bilayer embedded in between two flakes of the same material (in our case, hBN), thickness and size conserves this symmetry. Therefore, one can extend the k · p model proposed in [30] also for this case.
Let us examine the Bloch states of the valence and conduction bands at the K + point of bilayer. We construct them from the Bloch states of top and bottom layers considering them separately. Namely, we introduce the states |Ψ  (1) c are made predominantly from d x 2 −y 2 +id xy and d z 2 orbitals of transition metal atoms respectively [45,46]. They coincide with monolayer spinless states |Ψ v and |Ψ c mentioned before. The top (second) layer states |Ψ The initial basis states in a given layer are affected by crystal fields of another layer and surrounding hBN medium. Such fields being considered as a perturbation in k · p model produce intra-and interlayer corrections to the bilayer Hamiltonian. Namely, the intralayer ones renormalize the band gap E g and spin-splittings ∆ c , ∆ v for the electron excitations in the considered layer. Due to the symmetry of the system, the parameters of the other layer get the same modifications. The interlayer corrections link the states from different layers. The symmetry analysis of such terms demonstrates i) the strong coupling between VB basis states with the same spins; ii) the quasi-momentum dependent coupling between CB basis states with the same spins; iii) the coupling between VB and CB basis states of the opposite spins, which appears due to spin-orbit interaction. The latter term is supposed to be small and is omitted from our study. In this case, the spin-up and spin-down states of bilayer are decoupled and can be considered separately.
As a result, the VB Hamiltonian written in the basis {|Ψ |∆c| and √ ∆ 2 v + 4t 2 denote the splitting in the CB and VB, respectively. Eg is the single particle band gap.

The bilayer VB have the energies E
The corresponding upper-energy eigenstates are where we introduced cos(2θ) = ∆ v / ∆ 2 v + 4t 2 . The lowenergy eigenstates |Φ − v↑ and |Φ − v↓ can be derived from the first ones by replacing cos θ → − sin θ, sin θ → cos θ. The new VB are doubly degenerated by spin.
The CB states do not interact with each other in K + valley (i.e. in k x , k y → 0 limit). Hence, in leading order they are not mixed and form doubly degenerated bands with energies E g ± ∆ c /2. In further we focus on MoS 2 bilayer with ∆ c < 0. For this case, the upper and lower energy CB states are {|Ψ All new energy states of bilayer system are depicted graphically in Fig. 3. We divide them on spin-up (left) and spin-down (right) subsets for clarity. The single and doubled bumps in the figure represent the new CB and VB states. The size of the bump describes the probability to observe the new composite quasiparticle in the bottom (green part) or the top (red part) layers. Note that the valence band states can be found in both layers. That makes these states optically active in both polarizations. Namely, the transitions from any valence band state into |Ψ There are two intralayer exciton transitions in bilayer All of them have the same intensity I = I 0 cos 2 θ, where I 0 is the intensity of exciton line in monolayer.
There are two interlayer exciton transitions in bilayer They have the intensities I ′ = I 0 sin 2 θ. Finally, since the intra-and inter-layer exciton transitions are active in opposite circular polarization at a given K point, they should have the opposite signs of g-factors (see Appendix B). For a trilayer case, a crystal can be presented as three monolayers each separated by distance l, with the middle layer is 180 • rotated around the external ones (see Fig. 4). Similarly to the bilayer case, we focus on properties of quasiparticles at the K + point. The basis states of a trilayer are {|Ψ They belong to z = −l, z = 0 and z = l layers, respectively. In this case, the crystal has mirror symmetry σ h , with the mirror plane z = 0. It determines the symmetry relations |Ψ n . The properties of |Ψ (1) n according to the crystal symmetry. Hence, in K + point, the 1-st and 3-d layers absorb only σ + polarised light, while the 2-nd one is active in σ − polarization. Note that the trilayer encapsulated in between to equal hBN flakes conserves the mirror symmetry.
We derive the effective k · p Hamiltonian for considering system within several approximations: i) the interlayer spin-orbit coupling between CB and VB states with opposite spins from different layers is neglected, like in bilayer; ii) the intralayer crystal field corrections are equal for each layer, i.e. the bands of each layer are characterised by the same band-gap E g and spinsplitting ∆ v , ∆ c parameters, the magnitudes of which, however, can deviate from their bilayer analogs; iii) we neglect the coupling between the states of the 1-st and 3-d layers. In this approach, like in bilayer case, the spin-up and spin-down states of trilayer are decoupled. The Hamiltonian for VB states, written in the basis {|Ψ The Hamiltonian for conduction bands, written in the basis {|Ψ The system possesses a mirror symmetry. Hence, it is convenient to introduce the new basis states which are even |Υ The 2 × 2 blocks have the structure of bilayer Hamiltonians. Therefore, all the properties of the trilayer even states can be obtained from bilayer result by replacing t → √ 2t, u → √ 2u. Consequently, there are doubly degenerate valence and conduction bands, which define the intensive T A and T B and weak T ′ A and T ′ B groups excitonic transitions, which are active in opposite polarizations of light (see Fig. 5(a)). The 1 × 1 block has the structure of monolayer Hamiltonian. Hence, there are only T o A and T o B transitions between odd states of trilayer. They are active in the same polarization as T A and T B transitions of monolayer (see Fig. 5(b)). Summarising, the trilayer possesses two types of "even" intralayer excitons -T A and T B , two types of "odd"intralayer ones -T o A and T o B . All of them are active in σ + polarization at the K + point and have the g-factors of the same sign. The interlayer "even" excitons -T ′ A and T ′ B , conversely are active in σ − polarization at the K + point and have an opposite sign to intralayer excitons g-factors. The detailed description of trilayer's g-factors is presented in Appendix C.
Note that the current theoretical description predicts the number of exciton resonances and their polarization properties. However, the relative position of T A and T o A exciton remains an open question. Indeed, both lines are active in the same polarization and have similar g-factors, that makes impossible to identify each of them from the experiment. The analysis of the interference effects together with the consideration of the Coulomb interaction between new quasiparticles answers this question, which is, however, beyond the scope of the current paper. In our case, we took the position of exciton resonances following the position of the valence and conduction bands, derived in single-particle k · p approximation, for definiteness.
Our single-particle considerations are summarized in Fig. 6. This picture exhibits the position and number of optically active transitions in mono-, bi-and trilayer S-TMD crystals.

III. ABSORPTION RESONANCES IN MOS2 MULTILAYERS
In the experiments, we have investigated a set of van der Waals heterostructures with optically active parts consisting of MoS 2 mono-(1ML), bi-(2ML) and trilayer (3ML). The structures were prepared on ultraflat unoxidized silicon substrates by combination of exfoliation and dry transfer techniques [48]. In order to achieve the highest possible quality of our samples, as well as to preserve the characteristic inverse and mirror symmetries of, respectively, 2H-stacked bi-and trilayer, the microcleaved MoS 2 flakes were sandwiched between two hBN flakes exfoliated from high-quality bulk crystals grown under high-pressure conditions [49]. Further information on the preparation of samples and their characterization is provided in Appendix A.
Experiments consisted of measurements of the re-flectance contrast spectra, which are defined as where R(E) is the reflectance spectrum when the light is focused on the MoS 2 flake and R 0 (E) is the reflectance from the region outside the flake. Two different microoptical setups have been used for measurements: one setup for measurements in the absence of magnetic fields and the second one for measurements in magnetic fields supported by either a 14 T superconducting magnet or a 29 T resistive magnet. A spatial resolution of about 1 µm (light spot diameter) is characteristic of our freebeam-optics setup used for experiments at zero magnetic field. Instead, the fiber-optics-based arrangement applied for experiments in magnetic fields provides a poorer spatial resolution, of about 10 µm. A spectral resolution of 0.32 nm has been assured for both setups. A combination of a quarter wave plate and a polarizer are used to analyse the circular polarization of signals. The measurements are performed with a fixed circular polarization, whereas reversing the direction of magnetic field yields the information corresponding to the other polarization component due to time-reversal symmetry. More on experimental details can be found in Appendix A. The collection of low temperature (∼ 5K) RC spectra using our B = 0 setup, measured in the spectral range of the optical gap (onset of strong absorption) of the MoS 2 layers is presented in Fig. 7. The observed, more or less pronounced resonances, show the dispersive spectral shape, correspond to typically, strongly bound direct bandgap excitons in S-TMD layers. Although our theoretical considerations neglect the effects of Coulomb binding and account only for interband transitions, a distinct resemblance between the measured spectral evolution (Fig. 7) and that theoretically concluded in Fig. 6 can be recognized.
In the following we assume that each interband transition T X is associated with the corresponding excitonic resonance X and we proceed with the assignment of the resonances observed in the experiment. We suppose that the resonating excitons are shifted in energy, by their binding energies, with respect to the associated interband transitions: E(X) = E(T X ) − E b (X). Binding energies might be, however, different for different excitons and this fact must be taken into account when comparing the energy position of X-resonances and T X -transitions. Important for the assignment of the observed resonances are the results of the polarization resolved measurements carried out in magnetic fields (see Figs. 8 and 9). As discussed above, we expect that our multilayers host two different types of intra-and interlayer excitons, each of them being distinguishable by their different polarization properties (g-factors of opposite sign and various magnitudes). At this point we must admit that all subtle spectral features which are well visible in the RC spectra measured at B = 0 with high spatial resolution are somewhat less pronounced in the magneto-optical measure- ments which imply worse spatial resolution. Certain degree of sample inhomogeneity is an obvious cause of this drawback. In consequence, in the case of weak and/or broad resonances the information about their g-factors is not easily extractable from the raw magneto-RC data. This information becomes more apparent if we inspect the RC-polarization spectra. These spectra have been constructed as a difference between the RC-spectra measured at the same strength but for two opposite direction of the magnetic field (that mimic, due to the time reversal symmetry, the spectra corresponding to σ + and σ − circular polarization of the reflected light [39,50]). As can be deduced from the data shown in Figs. 8 and 9, the apparent dips in our RC polarization spectra correspond to excitonic resonances with negative g-factors whereas the characteristic upswings mark the resonances with positive g-factors.
Very first classification of the observed resonance takes into account a relatively large spin-orbit splitting in the MoS 2 valence band, rising two groups of excitons associated with well-separated upper (A-excitons) and lower (B-excitons) valence band subbands. The 1ML spectrum shown in Fig. 7, resembles that previously reported for a high quality 1ML MoS 2 encapsulated in hBN [44,51]. It depicts three resonances: one well-separated resonance due A-exciton (ground state), X A , and two other superimposed resonances, due to B-exciton (X B ) and the first excited (2s) state of A-exciton (X 2s A ). X A and X B resonances are associated with T A and T B transitions sketched in Fig. 1; the excited excitonic states are beyond the frame of our single particle theoretical approach. As expected, the X A , X B and X 2s A excitonic resonances display similar polarization properties in the magneto-optical experiments (see Fig. 8). We estimate g XA ≈ −4 in agreement with previous reports [52][53][54]; the accurate estimation of the amplitudes of g XB and g X 2s A is more cumbersome though both these values are also negative.
The RC spectrum of the 2ML (Fig. 7) shows 4 resonances [38] which, in reference to our theoretical expectations (see Figs. 3 and 6), are assigned to the pair of intra-(X A ) and interlayer (X ′ A ) A-excitons, and to the analogous pair of intra-(X B ) and interlayer (X ′ B ) Bexcitons. We use the signs and magnitude of the exciton g-factors as fingerprints the intra-and interlayer nature of the X A -and X ′ A -excitons respectively [39]. Namely, we conclude that g XA ≈ −4, while g X ′ A ≈ 8. Again, the exact values for the g-factor amplitudes of B-excitons are hard to be precisely estimated. Nonetheless it is rather clear that g XB is negative whereas g X ′ B is positive. Whereas our theoretical model meets the experimental data at the qualitative level, the energy ladder of the observed resonances requests further comments. Characteristic for the theoretical modelling is the fact that the energy distance, ∆(T ′ A , T A ), between the T ′ A and T A transitions is the same as the energy separation, ∆(T B , T ′ B ) between the T B and T ′ B transitions, both differences being determined by the amplitude ∆ c of the spin orbit splitting in the conduction band of the 2ML: This property is not seen in the experiment: we estimate that ∆(X ′ A , X A ) ≃ 70 meV, whereas ∆(X B , X ′ B ) ≃ 30 meV. The inconsistency between theory and experiment is likely due to different binding energies of excitons associated with different interband transitions: E(X) = E(T X ) − E b (X). In the first approximation we assume that binding en- ergies of our inter-and intralayer excitons are indeed different but that the excitons within each pair of indirect and direct resonances are the same: . Then, comparing the theoretical prediction with the experimental data one obtains that ]/2 ≃ 50 meV and that E b (intra) = E b (inter) + 20 meV. Larger binding energies of intralayer excitons with respect to those of interlayer excitons are logically expected since Coulomb attraction should be indeed stronger/weaker when the electron and hole charges are localized in the same or in the neighboring layers, as for the case of intra-or interlayer excitons, respectively. At first sight, the estimated amplitude of the spin orbit splitting in the conduction band of the MoS 2 bilayer appears to be surprisingly large. In the case of 1L MoS 2 , the commonly accepted results of the DFT calculations predict ∆ c in the range of few meV, though recent experimental works point out towards much higher values of about 15 meV [55]. According to our theoretical approach, the spin orbit interaction is sensitive to interlayer coupling and it is also affected by the interaction with the surrounding material (hBN). Thus the amplitudes of the ∆ c as well as ∆ v parameters are expected to evolve with the numbers of the stacked MoS 2 layers. Keeping our assumption about equal binding energies for X A and X B excitons we note that the expected energy separation between these excitonic resonance is given by |∆ c | + ∆ 2 v + 4t 2 (see Figs. 3 and 6), to be compared to the value of 170 meV estimated from the experiment. With the theoretically estimated value of t ≃ 40 meV [30] and the derived above |∆ c | ≃ 50 meV we find ∆ v ≃ 90 meV. Visibly, however, the "effective" spin orbit splitting in the valence band of the MoS 2 bilayer is larger, as given by ∆ 2 v + 4t 2 ≃ 120 meV (see Fig. 3).
Focusing now on the trilayer spectra (see bottom panel of Fig. 7) we observe a characteristic triplet of A-exciton resonances, which is in accordance with our theoretical expectations (see Fig. 6c). As expected, two strong X A and X o A resonances are characterised by negative g-factors whereas the g-factor of the interlayer X ′ A exciton is pretty much positive. The theoretically anticipated triplet structure of the ′′ X ′′ B -exciton is not resolved in the experiment but likely hidden within the observed broad spectrum of the B-exciton. Even with unresolved fine structure of the B-exciton for our trilayer MoS 2 , a rough estimation of ∆ c and ∆ v parameters can be done. Expecting that the energy separation between A , X A ) ≃ 70 meV from the experiment we conclude that |∆ c | ≃ 50 meV if the difference between binding energies of the intralayer X A and interlayer    binding energies of two X A and X o A intralayer excitons are the same we expect that they are separated in energy by ∆(X o A , X A ) = ( ∆ 2 v + 8t 2 − ∆ v )/2 and extract ∆ v ≃ 140 meV when reading ∆(X o A , X A ) ≃ 20 meV from the experiments and assuming again that t ≃ 40 meV.
With the above estimation, the band-edge structure of the 3L MoS 2 at the K ± points of the Brillouin zone is concluded to consist of two conduction band subbands split by 50 meV and 4 valence band subbands with outer two subbands split by 180 meV and the inner two sub-bands split by 140 meV.

IV. CONCLUSION
We have performed magnetooptical µ-reflectance measurements along with k · p theory based modelling on few-layer MoS 2 encapsulated in hBN structures, revealing the intralayer and interlayer nature of the newly discovered transitions. Such resonances form due to hybridisation of valence and conduction bands states when one adds new layers to the system.
The experiment i) manifests the new exciton resonances and their number in bi-and trilayers; ii) demonstrates that g-factors of intralayer and interlayer transitions have opposite signs; iii) the g-factor values of the latter ones are much larger (by magnitude) than of the former ones.
The symmetry based k · p description of the bi-and trilayers gives the natural explanation of the aforementioned experimental observations. Moreover, from the general symmetry point of view the theoretical model also predicts i) the renormalization of the band gap E g and spin splittings ∆ c , ∆ v ; ii) the deviation of the bi-and trilayer exciton g-factors from their monolayer analogs. These two phenomena appear as a result of the coupling between different layer of the system and effects of dielectric screening induced by hBN.
Finally, we point out some unique properties of the exciton states of the trilayer. First, we mention the existence of two groups of exciton resonances -"even" and "odd", which do not interact with each other because they belong to different irreps of the in-plane mirror symmetry of the system. However, electric field, applied perpendicularly to the crystal's plane, violates this symmetry and causes the controllable coupling of these states. Such feature can be used in future exciton based applications. We also note that the electrical positive and negative charges of interlayer excitons are separated between different layers. Such a property can be interesting in photovoltaic applications of S-TMD systems.

ACKNOWLEDGEMENTS
The work has been supported by the ATOMOPTO project (TEAM programme of the Foundation for Poli sh Science co-financed by the EU within the ERDFund), the EC Graphene Flagship project (no. 604391), the National Science Centre (grants no. DEC-2013/10/M/ST3/00791, UMO-2017/24/C/ST3/00119), the Nanofab facility of the Institut Néel, CNRS UGA and LNCMI-CNRS, a member of the European Magnetic Field Laboratory (EMFL).

Appendix A: Samples and experimental setups
In order to provide even and charged-defect-free support to van der Waals heterostructures studied in this work we assembled them on ultraflat unoxidized silicon (Si) substrates from Ted Pella Inc., whose surface roughness amounts to about one-tenth of that of standard Si/SiO 2 wafers. Each heterostructure consisted of three flakes obtained by means of mechanical exfoliation and successively stacked one on top of another using dry transfer techniques: a 60-70 nm thick bottom hBN flake, a central MoS 2 flake of 1ML, 2ML or 3ML thickness, and a 10-15 nm thick top hBN flake.
The optimum thicknesses of hBN flakes were estimated on the basis of transfer-matrix simulations of reflectance contrast spectra of complete heterostructures, aiming at maximizing the visibility of absorption resonances in MoS 2 flakes while keeping the thickness of bottom hBN flakes below 100 nm (a value up to which our exfoliation method yields a rich harvest of flat, defect-free flakes with large-area terraces of constant thickness).
Using high-purity ELP BT-150E-CM tape from Nitto Denko, bottom hBN flakes were exfoliated from highquality bulk crystals grown under high-pressure conditions [49] and then transferred in a non-deterministic way onto clean and freshly annealed at 200 • C Si substrates. The MoS 2 flakes were obtained by means of two-stage, tape-and polydimethylsiloxane (PDMS)-based exfoliation of a bulk crystal purchased from HQ Graphene and then stacked on the bottom hBN flakes using an all-dry deterministic stamping technique [48]. The same procedure was applied also to the top hBN flakes. While performing deterministic transfers, special attention was paid to mutual angular alignment of the flakes and to immediate 10-minute-long after-transfer annealing of the samples at 180 • C on a hot plate kept in clean ambient atmosphere.
Optical microscope images of investigated heterostructures which were fabricated in this way are shown in Fig. 10 (a)-(c). The spots of brown-to-blue colour visible in all images are up to 60-nm-high bubbles of air and/or water vapour (possibly with some amount of hydrocarbons) trapped between either the MoS 2 and the top hBN flake or the two hBN flakes. Their appearance directly results from the corrugation of thin, soft and flexible top hBN flakes supported by a PDMS stamp. Importantly, the areas between the bubbles, whose size reaches up to 5 by 10 micrometers, exhibit very flat and high-quality intefraces between the constitutent flakes, as revealed with optical measurements and atomic force microscopy (AFM) characterization.
An example of tapping-mode AFM imaging performed on finished heterostructures with the use of Digital Instruments Dimension 3100 microscope is shown in Fig. 10 (d)-(f). The false colour AFM images correspond to respective locations marked in panels (a)-(c) with orange circles. The grey arrows represent the line profiles drawn in Fig. 10 (g) and labelled with the same numbers as in images (d)-(f). As can be seen, the profiles unambiguously confirm the mono-, bi-, and trilayer thickness of MoS 2 flakes incorporated into the heterostructures shown in the upper row of Fig. 10.
Compared to based on neutron scattering measurements estimation equal to 0.615 nm [56], the value of 0.79 nm we got for the monolayer most probably indicates that the equilibrium distance between the hBN and MoS 2 layers differs from that between two neighboring MoS 2 layers in a bulk crystal. A smaller single-layer step hight obtained for the second and third layer in the biand trilayer MoS 2 flakes may on the other hand imply the existence of small uniaxial compressive strain along the c-axis in hBN/MoS 2 /hBN heterostructures.
Measurements at zero magnetic field were carried out with the aid of a continuous flow cryostat mounted on x − y motorized positioners. The samples were placed on a cold finger of the cryostat. The temperature of the samples was kept at T = 5 K. The excitation light was focused by means of a 50x long-working distance objective with a 0.5 numerical aperture producing a spot of about 1 µm. The signal was collected via the same microscope objective, sent through a 0.5 m monochromator, and then detected by a CCD camera.
Magneto-optical experiments were performed in the Faraday configuration using an optical-fiber-based insert placed in a resistive or a superconducting magnet producing magnetic fields up to 29 T or 14 T, respectively. The sample was mounted on top of an x − y − z piezo-stage kept in gaseous helium at T =4.2 K. The µ-RC experiments were performed with the use of 100 W tungsten halogen lamp. The excitation light was coupled to an optical fiber with a core of 50 µm diameter and focused on the sample by an aspheric lens (spot diameter around 10 µm). The signal was collected by the same lens, injected into a second optical fiber of the same diameter, and analyzed by a 0.5 m long monochromator equipped with a CCD camera. A combination of a quarter wave plate and a polarizer are used to analyse the circular polarization of signals.

Appendix B: Bilayer in magnetic field
The magnetic field correction to the valence band Hamiltonian at the K + point written in the basis {|Ψ Here µ B is the Bohr magneton, B is the magnetic field, applied perpendicularly to the bilayer plane, 1 is the unit matrix and The correction to the conduction band Hamiltonian written in the basis {|Ψ c ⊗ |s } has the same form δH The expressions for g v , g c , δg v , δg c are derived in [37]. The corrections to the Hamiltonians provide the energy shifts of excitons in magnetic field. The intralayer X A and X B excitons, constructed from the quasiparticles at the K + point of bilayer, are active in σ ± polarizations. The corresponding energy shifts are linear in magnetic field δE σ ± A,B = ±g We introduced g u = 2m 0 u 2 / 2 ∆ c parameter, which originates from k-dependent admixture of conduction bands. The mixing is negligibly small in the absence of magnetic field, but gives finite correction if B = 0. The interlayer excitons X ′ A and X ′ B at the bilayer K + point are also active in σ ± polarizations with the magnetic shifts A ′ ,B ′ µ B B respectively. The corresponding g-coefficients are g (2) Note that there are the following relations between gcoefficients of interlayer and intralayer exciton transitions One can mention that the corresponding transitions at the K − point in σ ± polarizations are characterized by the same values ±g (2) A,B and ±g (2) A ′ ,B ′ as in K + point. This is the consequences of time reversal and inversion symmetries of the crystal. As a result, we can restore the absorption spectra of the bilayer in σ − and σ + polarizations using our methodology of the experiment -by measuring the reflected light in fixed σ + polarization and changing the orientation of magnetic field from B to −B. The change of the direction of the magnetic field to −B in experiment mimics the transitions in σ − polarization of the reflected light.
We use the formula As a result the intralayer and interlayer A, B excitons have 2g (2) A,B and 2g (2) A ′ ,B ′ g-factors, respectively. The possible signs of intra-and interlayer exciton gfactors can be obtain from the expressions (B3), (B4), (B5) and (B6). Indeed, in the experiment we have 2g (2) A,B ≈ −4. It means that corrections g u , δg c and δg v are not significally large, and therefore very roughly 2g (2) A,B ≈ 2(g c − g v ) is nothing but the g-factors of monolayer A and B excitons, for which we know that g v > g c > 0. Substituting the positive g c and g v into Eqs. (B5) and (B6) one can see that 2g (2) A ′ ,B ′ > 0, which is confirmed from the experiment.
⊗ |s }, is δH Here n = v, c andḡ c andḡ v are the additional parameters which describe the magnetic dependent coupling between layers of the system [37]. In new basis, defined in the main part of the text, the full Hamiltonians H (C2) The 1 × 1 block corresponds to odd states |Υ (3) ns ⊗ |s with total energies with E v = 0 and E c = E g . The expressions coincide with the monolayer ones, but with the new g-coefficient g n − δg n +ḡ n . Moreover, the uk ± terms do not affect the odd states, and therefore such excitons have the same reduced masses as their monolayer analogs. Hence the corresponding trilayer exciton line has the same optoelectronic properties as it's monolayer analog. The odd A and B exciton transitions at the K + point are active only in σ + polarization. The corresponding energy shift in magnetic field is δE σ + Here we introduce the parameter From the experiment we know that g X o A ≈ −4.5, which is close to monolayer g XA = 2(g c − g v ) ≈ −4. Therefore, we can roughly estimate g ≈ 0.125.
The second block is 2 × 2 matrix, written in the basis of even states {|Υ (1) n ⊗ |s , |Υ (2) n ⊗ |s }. Note that the structure of this matrix does not coincide with the structure of G (2) n . Therefore the subsystem of even states of trilayer demonstrates another behavior in magnetic field than it's bilayer analog. The conduction band states in K + point are decoupled and have the energies The conduction band energies correspond to the states {|Υ The energy shifts for interlayer A ′ and B ′ excitons in both polarizations have the form δE σ ± A ′ ,B ′ = [±g In the absence of g the latter results coincide with the bilayer case up to redefinition of the parameters. The non-zero value of g shows remarkable difference between pure bilayer and effective bilayer cases. The corresponding energy shifts in σ ± polarizations at the K − point are δE σ ± A,B = [±g This non-equivalency makes the analysis of the gfactors of the system more complicated. Let us consider the results of the measurements presented on the Fig. 8 more carefully, focusing mainly on X A and X ′ A exciton transitions. Again, according to the methodology of our experiment we measure the g-factor using the formula g = [δE σ + (B) − δE σ + (−B)]/µ B B. Then the g-factors at the K + point of intralayer X A , X B and interlayer X ′ A , X ′ B excitons are g K + XA,B = 2g (3) A,B + 2g and g K + X A ′ ,B ′ = 2g (3) A ′ ,B ′ + 2g respectively. For K − point transitions we obtain g K − XA,B = 2g (3) A,B − 2g and g K − X A ′ ,B ′ = 2g (3) A ′ ,B ′ − 2g. Taking into account the relative smallness of g and absence of the results for magnetic field B larger than 14 T we suppose that the double g-factor structure of X A and X ′ A resonances is indistinguishable. Instead of this, probably, we observe only their average values g XA,B = 2g and g X A ′ ,B ′ = 2g (3) A ′ ,B ′ respectively. One can mention that these average g-factors surprisingly coincide with the ones we can get from the standard formula g = [E σ + (B)−E σ − (B)]/µ B B = [δE σ + (B)−δE σ − (B)]/µ B B. The analysis of the signs of intra-and interlayer exciton g-factors can be done in the same way as in bilayer case.