Shubnikov-de Haas oscillations in optical conductivity of monolayer MoSe$_2$

We report polarization-resolved resonant reflection spectroscopy of a charge-tunable atomically-thin valley semiconductor hosting tightly bound excitons coupled to a dilute system of fully spin- and valley-polarized holes in the presence of a strong magnetic field. We find that exciton-hole interactions manifest themselves in hole-density dependent, Shubnikov-de Haas-like oscillations in the energy and line broadening of the excitonic resonances. These oscillations are evidenced to be precisely correlated with the occupation of Landau levels, thus demonstrating that strong interactions between the excitons and Landau-quantized itinerant carriers enable optical investigation of quantum-Hall physics in transition metal dichalcogenides.

Optical excitations of a semiconductor hosting a twodimensional electron (2DES) or hole (2DHS) system subjected to a strong magnetic field provide a direct tool for investigating quantum-Hall states arising from Landau-level (LL) quantization of the carrier orbital motion [1][2][3]. Over the last three decades, this approach has been widely applied to explore a plethora of fascinating many-body phenomena in III-V or II-VI quantum-well heterostructures, such as, modulation of the screening responsible for magneto-oscillations of the excitonic luminescence at integer [4,5] and fractional filling factors ν [6][7][8], or formation of collective spin excitations (Skyrmions) around ν = 1 [9,10].
Qualitatively new avenues for the magneto-optical studies of 2DES or 2DHS have emerged with the advent of atomically-thin transition metal dichalcogenides (TMD), which are direct band-gap semiconductors with two nonequivalent valleys of conduction and valence bands appearing at the K ± points of the hexagonal Brillouin zone [11][12][13][14]. Owing to heavy effective carrier masses and reduced screening, the cyclotron energy in a TMD monolayer is much smaller than the binding energies of the exciton (0.5 eV) or trion (30 meV) [15][16][17] even at high magnetic fields ∼10 T. In this regard, TMD monolayers remain in stark contrast to conventional semiconductors, such as GaAs. Moreover, spinvalley locking by the spin-orbit coupling [18] as well as the presence of the valley-contrasting π Berry phase [12] give rise to a unique ladder of spin-and valley-polarized LLs in TMD monolayers [19,20], the degeneracy of which is further lifted by the strong Zeeman effect [21,22]. Although the fingerprints of such LLs have been recently reported in several transport experiments [23][24][25][26][27][28], their optical signatures have been obtained only for the WSe 2 monolayer at high electron densities, where exciton binding is suppressed due to screening; in this limit, the optical excitation spectrum is similar to that of GaAs 2DES at moderate fields and shows band-to-band inter-LL transitions [29].
In this Letter, we report the optical response of a Landauquantized 2DHS in a MoSe 2 monolayer placed in a strong magnetic field of 8-16 T in the opposite limit of low densities p < ∼ 3 · 10 12 cm −2 , where interaction effects are manifest.
For such low hole densities (corresponding to a Wigner-Seitz radius of r s > ∼ 5), the excitons remain tightly bound and the itinerant holes are fully spin-and valley-polarized in the K + valley (for B > 0). As a consequence, the optical excitation spectrum is dominated by two cross-circularly-polarized resonances corresponding to the bare exciton in the K + valley and exciton-polaron in the K − valley [30][31][32][33]. We demonstrate that the transition energies and/or linewidths of both resonances exhibit Shubnikov-de Haas-like (SdH-like) oscillations with the hole density, which are correlated with LL filling, thus providing a first direct evidence for the influence of quantum-Hall states on the excitonic excitations in a TMD monolayer.
Our experiments have been carried out on three different van der Waals heterostructures, each consisting of a gate-controlled MoSe 2 monolayer encapsulated between two hexagonal boron-nitride (h-BN) layers. In the main text, we present the data acquired for one of these devices [ Fig. 1(a)], where density-dependent oscillations are manifest in the hole-doped regime; while all 3 devices show such oscillations in the presence of a 2DHS, the third device, described in the Supplemental Material (SM) [34], displays analogous oscillations for a 2DES. The valley-selective optical response of all samples has been analyzed by means of low-temperature (T ≈ 4 K), circular-polarization-resolved, white light reflection magneto-spectroscopy (for details, see SM [34]).
We first concentrate on the spectral reflectance contrast measured as a function of the gate voltage V g at B = 0, which is shown in Fig 1(b). For −7 V < ∼ V g < ∼ 4 V, when the sample is devoid of free carriers, the spectrum features only the bare exciton resonance. Once V g is increased (decreased) beyond this range, free electrons (holes) start to be injected to the monolayer. As demonstrated previously [30][31][32][33], the attractive interaction between these carriers and the excitons dress the latter into exciton-polarons, which in turn qualitatively alters the nature of the optical transitions. A prominent signature of this crossover is an emergence of a second, lower-energy resonance originating from an attractive exciton-polaron. Concurrently, the exciton resonance transforms into a repulsive exciton-polaron resonance, which ex- hibits a strong blue-shift and line broadening. Due to the transfer of oscillator strength to the attractive polaron, the repulsive polaron resonance becomes indiscernible for carrier densities > ∼ 1 · 10 12 cm −2 . Figs 2(a,b) present the gate-voltage-dependence of the reflectance contrast spectra detected in two circular polarizations upon application of the magnetic field of B = 16 T perpendicularly to the monolayer plane. Such a field lifts the degeneracy of the electronic states in K ± valleys due to the Zeeman splitting [see Fig. 2(c)], whose magnitude has been previously found to be significantly enhanced by the exchange interactions in the presence of free carriers [35]. Consequently, at appropriately low densities the electrons (holes) occupy only the states in K − (K + ) valley. Because the polaron dressing of the excitons arises predominantly due to the intervalley exciton-carrier interaction [31,33,35], at such densities the exciton-polaron resonances appear exclusively in σ + (σ − ) polarization, whereas the spectrum in the opposite polarization shows only the bare exciton resonance. As seen in Figs 2(a,b), this valley-polarized regime covers the whole experimentally accessible gate-voltage range −16 V < ∼ V g < ∼ −7 V on the hole-doped side, while for the electrons it holds for 5 V < ∼ V g < ∼ 10 V. Most importantly, there is a striking difference between these two cases: at electron doping the optical transitions evolve smoothly with V g , whereas at hole doping they exhibit a pronounced oscillatory behavior. As we show below, these oscillations are due to sequential filling of the hole LLs.
Although we have also observed similar oscillations on the electron-doped side for a better quality sample featuring narrower optical transitions (see SM [34]), they have been found to be much less pronounced than the oscillations at the hole doping for all studied devices. We speculate that this asymmetry may be, at least partially, a consequence of electronto-hole effective mass ratio that exceeds unity, in contrast to density functional theory (DFT) calculations predicting both masses to be similar [36][37][38]. This conjecture is supported by the fact that while the hole mass was shown by ARPES measurements [39][40][41][42] to remain in excellent agreement with DFT, recent transport measurements of MoS 2 and MoSe 2 monolayers indicated that the electron mass may be larger by a factor of ∼ 2 [26,28].
In the following we focus on the central finding of our work-the signatures of LL filling in the optical spectra of the exciton and attractive polaron transitions in the holedoped regime (the repulsive polaron transition is not examined, since it disappears at very low hole densities). As seen in Figs 3(a-c), both resonances exhibit aforementioned oscillatory behavior, which is particularly striking for the exciton, whose linewidth displays sharp, periodic minima [Fig. 3(a)] that appear to be correlated with cusp-like changes of the

Exciton
Attractive polaron 1660 1620 1660 Photon energy (meV) Photon energy (meV) Reflectance contrast Gate voltage (V)  slope of the energy increase with V g [Figs 3(b,c)]. Remarkably, this effect may be independently understood as SdH oscillations in MoSe 2 optical conductivity σ(E X ) at the excitonic energy [43], since Re[σ(E X )] is determined by the excitonic linewidth [44]. However, the origin of those oscillations is different to that of transport SdH oscillations in static conductivity. Specifically, we find that the oscillations in our system are due to the influence of the LL occupation on the strength of the interactions between the exciton and Landau-quantized 2DHS. This is particularly clear for the case of the exciton linewidth, the narrowing of which arises directly from the modulation of the efficiency of the intravalley exciton-hole coupling. This coupling becomes enhanced each time the lowest-energy LL is only partially filled (around half-integer ν), as in such a case the holes can be effectively scattered between the empty states belonging to this LL. Concurrently, when the Fermi level lies in the gap between the LLs (around integer ν), the possibility of such a hole scattering is substantially reduced due to lack of available final states. This leads to a suppression of the exciton-hole interaction, and thereby gives rise to a pronounced minima in the exciton linewidth.
The above reasoning is corroborated by the results of our theoretical simulations of the exciton absorption spectrum, which, for simplicity, were intended to reproduce the LLrelated oscillatory features on a qualitative level only. First, using a variational approach, we calculate the energy of the exciton embedded in the reservoir of valley-polarized holes, taking into account both the phase-space filling and hole-hole exchange interaction. Then, assuming contact-like repulsive intravalley exciton-hole interaction, we calculate (to second order in the interaction) the correlation energy and finite lifetime acquired by the exciton due to dressing with electronhole pairs from the Fermi sea (for details, see SM [34]). As shown in Figs 4(a-c), this simple model correctly captures periodic narrowing of the exciton transition around integer ν. Moreover, it also predicts the slope of exciton energy dependence to be altered for the same hole densities, thus explaining the presence of cusp-like features in Figs 3(b,c). We emphasize that although analogous oscillatory features appear in the model including only the phase-space filling and ignoring interactions [see the black curves in Figs 4(b,c)], in such a case the cusps have an opposite direction to that in the experiment: the slope of the exciton energy gets steeper around integer ν instead of becoming flatter. This finding unequivocally demonstrates that the correct description of the experimental data requires inclusion of the interactions, which in turn confirms their primary role in LL-filling-dependent modification of excitonic spectra.
Remarkably, the signatures of LL-filling emerge also for the K − exciton dressed into an attractive polaron, but are found to be much less pronounced. In fact, they are barely visible in the gate-voltage evolution of the polaron linewidth [ Fig. 3(a)] and appear only in the V g -dependence of the energy, taking form of a familiar cusp-like features around integer ν that, however, tend to vanish for larger hole densities [ Fig. 3(b,c)]. Although the reason behind this tendency remains not entirely clear, the presence of even faint polaron energy oscillations is by itself supportive of our explanation of the LL signatures to emerge in the optical spectra due to the interactions. This stems from the fact that, unlike the exciton, the polaron transition is not affected by the phase-space filling in the investigated valleypolarized regime [compare the schematics in Figs 3(e,f)].
As demonstrated above, the minima of the exciton linewidth occur at gate voltages corresponding to subsequent integer ν, which should imply these voltages V g (ν, B) to be equally spaced, as long as the hole density p remains proportional to V g . Although this prediction holds for ν > 3 [see Fig. 3(d)], at lower densities the inter-LL voltage gaps turn out to be significantly narrower. We attribute these deviations (that are found to be specific to the investigated device, see SM [34]) to the presence of a Schottky barrier at the contact to the MoSe 2 , which gives rise to initially non-linear increase of p with V g (for details, see SM [34]). This effect, however, becomes negligible at larger-density regime (i.e., for The data are presented only in a linear-response regime of the gate (i.e., for Vg < ∼ −8 V), in which the hole density depends linearly on Vg, and its value (marked on the right axis) can be obtained within the frame of the parallelplate capacitor model (for details, see SM [34]). Solid lines represent the fit to the experimental data with a set of linear dependencies Vg(ν, B) = V0 − δνB with two fitting parameters: V0 representing the common origin of these curves, and δ controlling their slopes.
needed to fill a LL in a unit magnetic field. On this basis, as well as based on the geometrical capacitance C geom of our device (obtained within a parallel-plate capacitor approximation, see SM [34]), we finally extract the number of states in each LL C geom δB/e that yields (1.0 ± 0.2)eB/h (where e is the electron charge, while h the Planck constant). This finding indicates complete lifting of the LL degeneracy, exactly as expected for LLs that are valley-and spin-polarized.
In conclusion, our experiments provide the first signatures of LL-quantization in the optical spectra of a TMD monolayer hosting a dilute system of spin-valley polarized holes. These signatures are evidenced to emerge due to the influence of LL occupation on the strength of interactions between itinerant holes and tightly bound excitons, which gives rise to prominent filling-factor-dependent SdH-like oscillations in the energy and linewidth of the excitonic transitions. The interaction-enabled optical access to the quantum-Hall physics demonstrated in our work constitutes the first step towards optical investigation of a rich field of strongly correlated phenomena at integer and fractional filling factors in atomically thin semiconductors.
The authors acknowledge discussions with Richard Schmidt. This work is supported by a European Re-

S1. SAMPLE AND EXPERIMENTAL SETUP
The device investigated in the main text consists of a MoSe 2 monolayer, which is contacted with a few-layer graphene (FLG) flake and fully encapsulated between two thin films of hexagonal boron-nitride (h-BN). These films are sandwiched between a pair of FLG flakes serving as top and back gates, out of which only the top gate was utilized to control the carrier density in the present study. All of the flakes were mechanically exfoliated from the bulk crystals (HQ Graphene MoSe 2 , NIMS h-BN, and natural graphite) using wafer dicing tape (Ultron) and then subsequently transferred onto Si substrates covered by a 285 nm thick SiO 2 layer. At this state the quality and thicknesses of the flakes were verified by optical contrast measurements [S1-S3] and/or atomic force microscopy (AFM). Moreover, the thicknesses of top and bottom h-BN layers, yielding respectively t t = (81 ± 5) nm and t b = (24 ± 3) nm, were deliberately chosen based on transfer-matrix simulations of the device reflectance spectrum [S4] in order to ensure that: (1) the excitonic resonance exhibits approximately Lorentzian lineshape, (2) the MoSe 2 monolayer is placed around a node of the optical field confined inside an effective cavity formed by the heterostucture, with the aim to prolong the radiative lifetime of the exciton [S5], and hence to reduce the linewidth of its optical transition.
The multilayer stack was assembled by means of a dry-transfer method [S6, S7] involving the use of a glass slide holding a hemispherical polydimethylsiloxane (PDMS) stamp [S8] covered with a thin polycarbonate (PC) layer to sequentially pick up the flakes. All stacking steps were performed inside the glove box in inert Ar atmosphere. The stacking was carried out in a high temperature of 120 • C, which, together with a slow stacking speed, allowed us to reduce the number of contamination pockets [S7]. In this regard, the utilized hemispherical shape of the PDMS stamp was a helpful aid, as it facilitated the control over the size of the contact front between the PC and the substrate. The finished stack was released onto the SiO 2 /Si substrate, and afterwards the residual PC film was removed from its surface by immersing in chloroform. No annealing was performed during the whole fabrication process. Finally, the FLG flakes were electrically contacted with metal electrodes prepared by standard electron-beam lithography and subsequent deposition of 105 nm thick gold layer on top of a 5 nm thick titanium sticking layer. The MoSe 2 monolayer was doped by applying a voltage V g between the FLG top gate and the FLG contact to the monolayer. The resulting carrier density was, however, found to be determined not only by the value of V g , as the device displayed a pronounced hysteresis when V g was swept in a loop [see Figs S1(a,b)], which presumably stems from a Schottky nature of the contacts. In order to avoid such a hysteretic behavior, in all of the experiments reported in the main text the gate voltage was always varied within a fixed range, at a constant rate, and in the same direction (i.e., from V g = 16 V to −16 V). In this way we were able to almost completely mitigate the impact of gate hysteresis on the measured quantities, as confirmed, e.g., by nearly perfect reproducibility of the V g values corresponding to the onset of filling the valence band with holes, which for all gate-voltage scans varied by less than 0.1 V. Our magneto-optical experiments were carried out in a high-resolution, confocal microscope setup schematically depicted in Fig. S1(c). The sample was mounted on x-y-z piezo-electric stages inside a stainless steel tube filled with 20 mbar helium exchange gas at T ≈ 4 K. The tube was immersed in a liquid helium bath cryostat equipped with a superconducting coil generating a magnetic field of up to 16 T in Faraday geometry. A free-space optical access to the sample was provided by a wedged window on top of the tube. The reflectance measurements were performed with the use of a single-mode-fiber-coupled broadband light emitting diode (LED) with a center wavelength of 760 nm and a 3 dB bandwidth of 20 nm. The excitation light, after exiting the fiber, was collimated before it entered the tube, where it was focused on the sample surface by a high numerical aperture aspheric lens (NA = 0.68). The excitation power was kept in the range of a few tens of nW. The light reflected off the sample was collected by the same lens, separated from incident light by a beam splitter, and finally coupled to a single-mode fiber that guided it to 0.5 m spectrometer equipped with a liquid nitrogen-cooled charge-coupled-device (CCD), which was utilized to analyze the reflectance spectra. A set of polarization optics, including linear polarizers and quarter wave-plates, was incorporated in both excitation and detection paths in order to linearly polarize the incident beam as well as to detect the reflected light in σ + or σ − circular polarizations.

S2. ANALYSIS OF THE REFLECTANCE SPECTRA
To obtain the reflectance contrast from the measured spectra of the LED light reflected off the sample it is necessary to determine the unperturbed reflection spectrum of the LED. In principle, such a reference spectrum may be obtained by moving the excitation spot off the MoSe 2 monolayer region on the sample surface. However, due to inevitable spatial inhomogeneity of the heterostructure, this procedure is often fraught with a systematic error. For this reason we follow a more reliable approach, in which the reference spectrum is constructed based on the reflectance spectra taken at a fixed spot, but for different gate voltages.  More specifically, from a given circular-polarization-resolved gate-voltage scan we select two spectra in such a way that each of them exhibits only one resonance corresponding either to the exciton or attractive polaron, as shown in Figs S2(a,b) for an example data set taken at B = 16 T. Owing to narrow linewidths and large energy splitting of these optical transitions, such a spectra feature partially overlapping, resonance-free spectral regions, which are combined together to obtain an unperturbed reference spectrum R 0 (E) [Fig. S2(c)]. On this basis we evaluate the reflectance contrast spectra R c (E) = [R(E) − R 0 (E)]/R 0 (E) from the reflectance spectra R(E) measured at each gate voltage (E represents the photon energy).
In order to determine the transition energy and the linewidth of the exciton and attractive polaron resonances in the holedoped regime we first recall that the reflectance contrast measured in our experiments is not only determined by the MoSe 2 monolayer optical susceptibility χ(E). Instead, it can be effectively described as Im[e iα(E) χ(E)], where α(E) stands for a wavelength-dependent phase-shift arising from the interferences of light reflected at different interfaces (e.g., h-BN/SiO 2 ) of the heterostructure [S9]. A more comprehensive analysis of the lineshape can be carried out using the transfer matrix method; however, since the absolute strength of the susceptibility is not of interest in this work, we carry out the outlined, simpler approach. Due to a suitable choice of the thicknesses of top and bottom h-BN films (as described in Sec. S1), the invoked phase-shift is essentially close to 180 • within the spectral region of interest in case of our device. In particular, α 180 • for the exciton resonance, as it features almost purely Lorentzian lineshape [see Figs S3(a-c)] that, except for forming a dip, exactly corresponds to the expected Lorentzian spectral profile of the imaginary part of the susceptibility. We remark here that, in general, the assumed Lorentzian lineshape is only valid in the absence of free carriers, whereas for non-zero hole densities the spectral profile of the susceptibility may be modified by the exciton-hole interactions [S4, S10]. However, we do not find any clear signatures of the exciton lineshape variation with the gate voltage. As seen in Figs S3(a-c), the attractive polaron resonance appearing in the opposite circular polarization, exhibits a slightly asymmetric lineshape. As revealed by our transfer-matrix simulations, this asymmetry originates from a combined effect of a reduced monolayer reflectivity, different wavelength and narrower linewidth of this resonance, and hence it can be accounted for by setting α 169 • for the attractive polaron instead of α 180 • as in the case of the exciton.
Based on the above considerations, we describe the reflectance contrast spectral profiles of both resonances with the following formula: where A, E 0 , and γ represent, respectively, the amplitude, energy, and linewidth of a given resonance, while the parameter C is introduced to capture any broad background contribution to the measured signal. Using this formula we are able to almost perfectly reproduce the experimental data for different carrier densities, as shown by the solid lines in Figs S3(a-c). From such fits we extract the energies and linewidths of the exciton/polaron resonances as a function of the gate voltage, which, in particular, are plotted in Figs 3(a,b) in the main text. To evaluate the derivatives of the energies E with respect to the voltage V g [presented in Fig. 3(c)], the extracted data are binned in 80 mV intervals, and |dE/dV g | is subsequently computed as an absolute value of the difference quotient between neighboring data points. The voltages V g (ν) corresponding to integer filling factors ν [shown in Figs 3(d) and 5] are in turn determined based on the oscillations of the exciton linewidth. To this end, the original data are first convolved with a Gaussian of standard deviation 30 mV in order to reduce the noise, and then the searched voltages V g (ν) are determined as those corresponding to the data points with locally minimal linewidth.

S3. CAPACITIVE MODEL FOR EVALUATION OF THE CARRIER DENSITY
During our experiments, the hole density p in the MoSe 2 monolayer was tuned by applying a gate voltage V g between the FLG top gate and the FLG contact to the monolayer. In order to establish a link between the voltage V g and the density p, we employ a simple capacitive model, in which the effective capacitance per unit area of the device is given by where C Q = e 2 D(E F ) denotes the quantum capacitance determined by the density D(E F ) of electronic states (DOS) in the MoSe 2 at the Fermi level E F , while C geom stands for a geometric contribution, which, in the parallel-plate approximation, can be described as Here t h−BN = (81 ± 5) nm represents the thickness of the top h-BN layer separating the top gate and the MoSe 2 monolayer, whereas ⊥ h−BN denotes the perpendicular, static dielectric constant of h-BN, which is estimated as (3.5 ± 0.5) based on the values provided in several previous reports [S11, S12].
In this model, the quantum capacitance plays an important role only if C Q C geom , that is when D(E F ) is negligibly small (e.g., in the band gap). In such a case the carrier density remains constant, while the change of the gate voltage results solely in the shift of the Fermi level by dE F /dV g = e. Concurrently, in the opposite regime of large D(E F ), which is realized when the Fermi level lies in the valence or conduction bands, the quantum capacitance C Q C geom can in turn be neglected, and the carrier density increases linearly with the gate voltage according to dp/dV g = −C geom /e. We note here that when the band states break up into a series of Landau Levels (LLs) upon application of an external magnetic field B, the condition C Q C geom does not necessarily have to be fulfilled for the Fermi level lying in the gap between the LLs. However, for a realistic device, owing to inevitable disorder-induced LL broadening, the DOS in the inter-LL gap remains presumably large enough for the invoked condition to hold. But even in an ideal case of delta-like energy spectrum of Landau-quantized DOS, the deviations from a linear dependence of the carrier density on the gate voltage still turn out to be negligible. More specifically, in such a scenario the invoked dependence would exhibit a step-like character, featuring a periodic linear-increase regions separated by plateaus, each extending over the gate voltage range ∆V inter = E c /e, in which the Fermi level moves across the inter-LL cyclotron energy gap E c =heB/m * , where m * is the carrier effective mass. Crucially, this voltage range remains orders of magnitude narrower than the range of ∆V LL = ∆p LL e/C geom , in which the Fermi level is pinned to a LL that is gradually filled with holes, where ∆p LL = f eB/h denotes the number of states in the LL. This can be readily seen by computing the ratio ∆V inter /∆V LL = 2πh 2 f C geom /e 2 m * , which for f = 1 and m * ≈ 0.5m 0 [S13, S14] yields 2·10 −3 , showing that nonlinearities in the dependence of p on V g related to DOS oscillations can indeed be safely neglected.
In light of the above theoretical analysis, the change of the gate voltage ∆V LL required to fill a LL with holes should be independent of the filling factor ν. However, according to our experimental results, this prediction remains valid only for appropriately large ν [e.g., ν > 3 for the 16 T data depicted in Fig. 3(d) in the main text], whereas at low hole densities (p < ∼ 1 · 10 12 cm −2 ) ∆V LL is clearly getting smaller with decreasing ν. We attribute these deviations to the presence of a considerable Schottky barrier at the interface of the MoSe 2 monolayer and the FLG contact, which is not taken into account in the above-described capacitive model, but turns out to play an important role for the main device. In particular, such a barrier prevents the holes from flowing into the monolayer, making it necessary to apply lower (i.e., more negative) gate voltages in order to reach the hole-doped regime. This, for example, entails larger than expected separation between the onsets of filling the valence and conduction bands, which yields approximately 10 V when the voltage is being swept towards negative values [see Fig. S1(a)]. Once the Schottky barrier is overcome (at V g ≈ −6.5 V), the holes are abruptly injected to the monolayer with the rate |dp/dV g | much larger than the one of C geom /e predicted by the capacitive model, resulting in a decreased value of ∆V LL as well as relatively rapid blue (red) shift of the exciton (attractive polaron) resonances in the reflectance spectra. As the gate voltage is further ramped down towards more negative values, the Schottky effects become less and less important, which in turn leads to a gradual reduction of the hole injection rate. Finally, at V g ≈ −8 V, we enter a linear-response regime, in which the injection rate saturates at a constant value, and further changes ∆p of the hole density can be described within the frame of the capacitive model as −C geom ∆V g /e. This is fully confirmed by a perfect agreement between the gate voltages V g (ν, B) corresponding to integer ν that were extracted from the data in the invoked range and a set of linear dependencies of the form V g (ν, B) = V 0 − δνB, as shown in Fig. 5 in the main text. Based on the value of the parameter δ = ∆V LL /B = f e 2 /C geom h controlling the slopes of these dependencies, we extract the LL degeneracy f = 1.0 ± 0.2, as stated in the main text. Moreover, by extrapolating these dependencies to zero magnetic field, we determine the gate voltage V 0 = −2.57 V that would correspond to the onset of filling the valence band with holes in the absence of the Schottky effects. This allows us to obtain the absolute value of the hole density for V g < ∼ −8 V [marked on the right vertical axis in Fig. 5 in the main text] using the following relation:

S4. REPRODUCIBILITY OF THE RESULTS ON A DIFFERENT DEVICE
As stated in the main text, our experiments have been repeated on a second van der Waals heterostructure yielding consistent results. This heterostructure was fabricated using very similar technique to that utilized for preparation of the device studied in the main text (see Sec. S1). It consisted of a charge-tunable MoSe 2 monolayer, which was electrically contacted with a FLG flake, encapsulated between two layers of h-BN, and finally sandwiched between another two FLG flakes serving as top and back gates, out of which only the back gate was utilized to tune the carrier density in the present experiments. Figs S4(a,b) show circular-polarization-resolved reflectance contrast spectra measured for this device as a function of the back gate voltage under a magnetic field of B = 14 T. The data are of clearly lower quality than those obtained for the main device. In particular, the exciton/polaron resonances are significantly broader and accompanied by several lower-intense side peaks, which most probably appear due to the contributions of distinct, simultaneously excited sub-micrometre sized regions lying within the excitation spot. Nonetheless, both the main and side exciton resonances demonstrate qualitatively similar behavior to that seen for the main sample, featuring a pronounced Shubnikov-de Haas-like oscillations on the hole-doped side (i.e., in σ + polarization for V g < ∼ −1.5 V) that are certainly much fainter at the electron doping (i.e., in σ − polarization for V g > ∼ 1.5 V). In order to quantitatively compare the behavior of the exciton/polaron resonances for the two devices, the spectral characteristics of the main (i.e., most intense) resonances extracted for the second sample on the hole-doped side are plotted in Figs S4(c-f) in the same format as in Figs 3(a-d) from the main text. As expected from the above analysis, the behavior is found to be indeed very similar, confirming the robustness of our conclusions drawn in the main text. In particular, the polaron linewidth remains almost constant with the gate voltage V g , whereas the exciton linewidth is clearly increasing for more negative V g . Additionally, it displays aforementioned, periodic oscillations related to sequential filling of the LLs, which are correlated with the oscillations of the exciton energy. However, in contrast to the main device, the signatures of such oscillations are not particularly evident in the gate-voltage dependence of the attractive polaron energy, which most probably stems from larger broadening of the resonances in the present device.
Most importantly, the gate voltages V g (ν) corresponding to integer filling factors ν (determined based on the positions of the local minima of the exciton linewidth) are almost equidistant, as revealed by a linear decrease of V g (ν) with ν [see S4(f)]. This finding shows that the hole density remains proportional to V g in the whole experimentally accessible voltage range on the holedoped side. This, in turn, indicates that the Schottky effects at the MoSe 2 /FLG interface play much less important role in case of the second sample, which is further corroborated by much lower hysteresis observed for this device when sweeping V g in a loop. By fitting the values V g (ν) with a linear curve, we obtain the voltage change ∆V LL = |V g (ν + 1) − V g (ν)| = f e 2 B/C geom h ≈ 1.00 V needed to fill a LL. On this basis, as well as based on the value of a geometric capacitance C geom evaluated using Eq. (S3) for t h−BN = (50 ± 5) nm revealed by the AFM measurements of the bottom h-BN flake, we finally extract the LL degeneracy for the second device f = 1.1 ± 0.3, which is found to be equal to 1 within the experimental uncertainty, similarly to the case of the main device.

S5. SIGNATURES OF LANDAU LEVELS FILLING IN THE OPTICAL SPECTRA AT ELECTRON DOPING
One of the most surprising results of our analysis of the magneto-optical spectra acquired for each of two investigated devices (either the main one or that studied in Sec. S4) is the fact that the signatures of LL filling were found to be prominent only on the hole-doped side, whereas at the electron doping they were much weaker, if present at all. As stated in the main text, we attribute this asymmetry to be at least partially due to a difference in effective masses of both carriers, with the electron mass being noticeably larger, as suggested by the results of recent transport experiments on Mo-based TMD monolayers [S15, S16]. This hypothesis would imply the electron LLs to feature lower cyclotron energy spacing, resulting in lower amplitude of the electron DOS modulation in the presence of inevitable disorder-induced LL broadening. Following this argument, one may anticipate analogous Shubnikov-de Haas-like oscillations to emerge also in the optical spectra on the electron-doped side for a betterquality device exhibiting lower LL broadening. We find that this is indeed the case, as revealed by the data obtained for the third device, whose optical micrograph is depicted in Fig. S5(a). This sample was fabricated in a similar way to the other two devices, following the procedure outlined in Sec. S1. It consisted of a charge-tunable MoSe 2 monolayer that was electrically contacted with two separate FLG flakes, encapsulated between two h-BN films, and finally sandwiched between top and bottom FLG gates, which were both utilized to control the carrier density in the monolayer. Due to very similar thicknesses t = (80 ± 10) nm of the two h-BN layers (extracted based on the optical contrast measurements), the geometrical capacitances C geom [given by Eq. (S3)] between each gate and the MoSe 2 were almost identical, implying the electron density to be determined by the sum V g = V top + V bottom of the top and back gate voltages. Fig. S5(b) shows linearly-polarized reflectance contrast spectra measured at B = 16 T on a particular spot on this device as a function of V g in the voltage range corresponding to either neutral or valley-polarized n-doped regime (in which the electrons occupy exclusively the states in K − valley). As seen, both the exciton and exciton-polaron resonances are remarkably sharp, featuring roughly two times narrower linewidth as compared to the corresponding resonances for the main device. It is most probably the reason for which the optical transitions observed for the present device display prominent Shubnikov-de Haas-like oscillations related to the filling of the electron LLs, which remained almost unresolved for previously investigated samples. Importantly, the character of these oscillations turns out to be almost identical to that established for the hole LLs in the main text, as proven by Figs S5(c-f) presenting the spectral characteristics of the exciton/polaron resonances for the present device in the same format as in Figs 3(a-d) for the case of the hole doping of the main device. First, the linewidth of the exciton transition undergoes familiar, periodic oscillations, whose sharp minima are perfectly correlated with cusp-like changes of the slopes of the exciton and attractive polaron energy dependencies on V g . Second, these cusp-like features become clearly weaker at larger electron densities for the polaron. Finally, the gate voltages V g (ν) associated with integer ν (i.e., with the minima of the exciton linewidth) are almost equidistant, with the voltage change needed to fill a LL ∆V LL = |V g (ν + 1) − V g (ν)| = f e 2 B/C geom h ≈ 1.85 V corresponding to a LL degeneracy f = 1.2 ± 0.3 that is equal to 1 within the experimental uncertainty, exactly as in the case of the hole LL. Altogether, the above-demonstrated similarity between the optical signatures of the electron and hole LLs clearly indicates that the interaction-related mechanisms, invoked in the main text to explain the hole-doped case, are also responsible for LL-induced oscillations at electron doping, which generalizes the conclusions drawn in our work to the case of both types of carriers.

S6. THEORETICAL MODEL
In this section we develop a theoretical model of the exciton in a TMD monolayer hosting a dilute hole system subjected to a strong magnetic field. Similarly to the experiment, we focus only on the spin-valley-polarized regime, in which the itinerant holes occupy only the states in the K + valley, and consider optical excitations in the same valley. In such a case, the exciton will experience two main effects due to the presence of the hole system: (1) its energy will be modified owing to phase-space filling; (2) it will acquire a finite correlation energy and lifetime owing to the interaction with the itinerant holes.

A. Single particle Hamiltonian
In the presence of a magnetic field, it is useful to use the basis of single particle eigenstates. In the symmetric gauge A = B(−y/2, x/2, 0), which is employed here, these states can be labelled with two quantum numbers |nl , where n ≥ 0 denotes the LL index, while l ≥ 0 is the canonical angular momentum of the state (with l − n corresponding to the actual angular momentum). The single particle states |nl can be generated using the ladder operators: starting from the ground state r|00 = 1 √ 2π e −zz * /4 , where r denotes the position of the particle. Notice that throughout this section we will use the complex representation of the position z = x + iy, and work in units of B = h/eB = 1 and h = 1, where B stands for the magnetic length, e is the electron charge, whileh the reduced Planck constant. In the above expression, a † = (Π x + iΠ y )/ √ 2 and b † = (Γ x − iΓ y )/ √ 2 , where Π = −i∇ r + eA(r) denotes the kinetic momentum, while Γ = −i∇ r − eA(r) denotes the magnetic momentum.
The single-particle kinetic energy Hamiltonian for a monolayer in magnetic field can be obtained using the Peierls substitution [S17, S18]. In case of the analyzed K + valley, this Hamiltonian can be described as: where c † cnm (c † vnm ) denotes the creation operator of an electron in the conduction (valence) band in the state |nm , ∆ g represents the semiconductor band gap, ω c = eB/m * stands for the cyclotron frequency (with m * representing the carrier effective mass, which assumed to be equal to 0.5m 0 for both valence and conductions bands), and N φ = A/(2π 2 B ) denotes the Landau level degeneracy with A representing the quantization area. Notice that due to the chirality of the K + valley, the Landau levels in this valley are shifted up in energy by ω c /2 (while those in K − valley are shifted by the same energy in the opposite direction). We remark that we neglect here a small corrections to the wavefunctions/energies arising due to the Dirac nature of the material. This assumption is equivalent to neglecting deviations from the parabolic dispersion of the electrons in the absence of the magnetic field [S17, S18].
It is useful to introduce the electron and hole operators, such that: which allow us to write down the electron-hole kinetic energy Hamiltonian:

B. Exciton energy in a truncated Hilbert space
In the absence of the electron, the holes will reside in a many-body ground state ψ h g . It is beyond the scope of this Letter to accurately calculate this state. Instead, in the following simplified analysis we will assume that this state is an eigenstate of the angular momentum operator, and that it does not mix different LLs, i.e., ψ h g h † n l h nl ψ h g = ν n δ n n δ l l , where ν n represents the filling of the nth hole LL. We further assume that ψ h g is a Gaussian state, which allows us to calculate any correlation functions using Wick's theorem. We expect this assumption to be justified due to the significant strength of Coulomb interactions in TMD monolayers.
We are now going to calculate the energy of the exciton in the presence of the hole system. To this end, apart from the kinetic energy Hamiltonian H 0 , we need also to take into account electron-hole Coulomb attraction as well as hole-hole Coulomb repulsion, which can be described using the following Hamiltonians: H hh = 1 2 where r e (r h ) denotes the electron (hole) position, while q represents scattered momentum (in units of 1/ B ), which will be further expressed in a complex representation q = q x + iq y . In the above formulas, V q denotes the Fourier transform of the Coulomb interaction, which in a two-dimensional system can be approximated as: where t = b ≈ 3.5 denotes the dielectric constant of the h-BN films on top and on the bottom of the TMD monolayer. The term of ρ 0 = 4πχ 2D /( t + b ) accounts for the stronger screening inside the monolayer, where χ 2D represents the 2D polarizability of the planar material, which for MoSe 2 yields χ 2D = 0.823 nm [S19]. In order to determine the exciton ground state, we attempt to diagonalize full Hamiltonian H = H 0 + H eh + H hh in a truncated Hilbert space of a single electron-hole excitation, which is spanned by the states of the form h † n l e † nl ψ h g . Owing to a significant size of this Hilbert space, a direct numerical calculation of the exciton binding energy is bound to fail. This can be readily seen by realizing that, due to strong Coulomb interactions, the exciton ground state is a superposition of many singleparticle electron-hole excitations between electron and hole LLs of indices ranging up to ∼ ( B /a * B ) 2 , which yields about 40 for B = 16 T and the exciton Bohr-radius of a * B ≈ 1 nm [S20]. Bearing in mind that we also need to include a similar number of angular momentum states in each LL, this finally leaves us with roughly 40 4 states that would be required to properly calculate the exciton binding.
Fortunately, the above-introduced problem can be greatly simplified by taking advantage of the fact that we can eliminate the angular degrees of freedom due to quenching of the kinetic energy. This can be most conveniently done by introducing the density-like operators: where G functions are defined as [S21]: where L b a represents the generalized Laguerre polynomials. These functions satisfy the following relations: ∞ l=0 G m l k1 G lm k2 = e −k * 1 k2/2 G m m k1+k2 , e −|k| 2 /2 G n n k * G l l k = n l | e −ikr |nl .