Opto-valleytronic imaging of atomically thin semiconductors

TL;DR

This paper introduces two-dimensional polarimetry for direct imaging of valley pseudospin in monolayer transition metal dichalcogenides, enabling detailed analysis of valley polarization, coherence, and disorder in atomically thin semiconductors.

Contribution

It presents a novel polarimetric imaging technique for visualizing valley pseudospin in monolayer TMDs, with quantitative analysis and applications to device inspection.

Findings

Valley polarization and coherence are sensitive to crystalline disorder.

Site-dependent thermal and non-thermal regimes of valley excitons are identified.

Wide-field polarimetry enables rapid inspection of opto-valleytronic devices.

Abstract

Transition metal dichalcogenide semiconductors represent elementary components of layered heterostructures for emergent technologies beyond conventional opto-electronics. In their monolayer form they host electrons with quantized circular motion and associated valley polarization and valley coherence as key elements of opto-valleytronic functionality. Here, we introduce two-dimensional polarimetry as means of direct imaging of the valley pseudospin degree of freedom in monolayer transition metal dichalcogenides. Using MoS$_2$ as a representative material with valley-selective optical transitions, we establish quantitative image analysis for polarimetric maps of extended crystals, and identify valley polarization and valley coherence as sensitive probes of crystalline disorder. Moreover, we find site-dependent thermal and non-thermal regimes of valley-polarized excitons in perpendicular magnetic fields. Finally, we demonstrate the potential of wide-field polarimetry for rapid inspection of opto-valleytronic devices based on atomically thin semiconductors and heterostructures.

Equations

$$\begin{array}[]{llllrcl}\textrm{ground state}&\lvert 0\rangle&:&&\dot{\rho_{0}}&=&\rho_{1}\gamma_{0}+\rho_{2}\gamma_{0}-\rho_{0}\gamma_{03}\\ \textrm{$K$ state}&\lvert 1\rangle&:&&\dot{\rho_{1}}&=&\rho_{3}\gamma_{31}+\rho_{2}\gamma_{21}-\rho_{1}\gamma_{12}-\rho_{1}\gamma_{0}\\ \textrm{$K^{\prime}$ state}&\lvert 2\rangle&:&&\dot{\rho_{2}}&=&\rho_{3}\gamma_{32}-\rho_{2}\gamma_{21}+\rho_{1}\gamma_{12}-\rho_{2}\gamma_{0}\\ \textrm{excited state}&\lvert 3\rangle&:&&\dot{\rho_{3}}&=&\rho_{0}\gamma_{03}-\rho_{3}\gamma_{31}-\rho_{3}\gamma_{32},\end{array}$$

$$P_c=\kappa \frac{{\rho }_1-{\rho }_2}{{\rho }_1+{\rho }_2},$$

$$P_c=\kappa \left(\frac{{\gamma }_0}{{\gamma }_0+{\gamma }_{12}+{\gamma }_{21}}\cdot \frac{{\gamma }_{31}-{\gamma }_{32}}{{\gamma }_{31}+{\gamma }_{32}}+\frac{{\gamma }_{21}-{\gamma }_{12}}{{\gamma }_0+{\gamma }_{12}+{\gamma }_{21}}\right).$$

$$P_c=\kappa \left(\frac{{\gamma }_0}{{\gamma }_0+{\gamma }_{12}+{\gamma }_{21}}\cdot \frac{1-b_{\kappa }}{1+b_{\kappa }}+\frac{{\gamma }_{21}-{\gamma }_{12}}{{\gamma }_0+{\gamma }_{12}+{\gamma }_{21}}\right),$$

$$b_{\kappa }={\left(\frac{{\gamma }_{32}}{{\gamma }_{31}}\right)}^{\kappa }.$$

$$P_c^0=\frac{1}{1+2{\gamma }_l/{\gamma }_0}=\frac{1}{1+2r_0},$$

$$P_c^0=\frac{1}{1+{\tau }_0/{\tau }_l},$$

$$P_c^0=P_i^0\cdot \frac{1}{1+2r_0}=\frac{1-b_{\kappa }}{1+b_{\kappa }}\cdot \frac{1}{1+2r_0},$$

$$b_{\kappa }(B)=b_0+\kappa \delta B,$$

$$\begin{array}[]{rcl}\gamma_{12}&=&r_{0}\gamma_{0}\exp(-\beta)\\ \gamma_{21}&=&r_{0}\gamma_{0}\exp(+\beta),\end{array}$$

$$P_c=\frac{1}{1+2r_0\cosh (\beta )}\cdot \frac{1-b_{\kappa }}{1+b_{\kappa }}+\kappa \cdot \frac{2r_0\sinh (\beta )}{1+2r_0\cosh (\beta )}.$$

Full text

Opto-valleytronic imaging of atomically thin semiconductors

Andre Neumann1, Jessica Lindlau1, Léo Colombier1, Manuel Nutz1,

Sina Najmaei3, Jun Lou3, Aditya D. Mohite2, Hisato Yamaguchi2, and Alexander Högele1

1Fakultät für Physik, Munich Quantum Center, and Center for NanoScience (CeNS),

Ludwig-Maximilians-Universität München, D-80539 München, Germany

2MPA-11 Materials Synthesis and Integrated Devices, Materials Physics and Applications Division,

Los Alamos National Laboratory (LANL), Los Alamos, NM 87545, U.S.A.

3Department of Materials Science and NanoEngineering, Rice University, Houston, TX 77005, U.S.A.

Transition metal dichalcogenide semiconductors represent elementary components of layered heterostructures for emergent technologies beyond conventional opto-electronics. In their monolayer form they host electrons with quantized circular motion and associated valley polarization and valley coherence as key elements of opto-valleytronic functionality. Here, we introduce two-dimensional polarimetry as means of direct imaging of the valley pseudospin degree of freedom in monolayer transition metal dichalcogenides. Using MoS2 as a representative material with valley-selective optical transitions, we establish quantitative image analysis for polarimetric maps of extended crystals, and identify valley polarization and valley coherence as sensitive probes of crystalline disorder. Moreover, we find site-dependent thermal and non-thermal regimes of valley-polarized excitons in perpendicular magnetic fields. Finally, we demonstrate the potential of wide-field polarimetry for rapid inspection of opto-valleytronic devices based on atomically thin semiconductors and heterostructures.

Published in Nature Nanotechnology 12, 329–334 (2017).

The version of record is available online at https://doi.org/10.1038/nnano.2016.282.

The valley quantum degree of freedom manifests itself as the quantized angular motion of crystal electrons near the band edges. In analogy to spin it represents a resource for quantum information in conventional semiconductors such as aluminum arsenide 1 and silicon 2, or in atomically thin materials including graphene 3 and TMDs 4, 5. Direct band-gap monolayer TMDs 6, 7 are particularly viable for practical realizations of valleytronic concepts as they enable initialization 4, manipulation 8 and detection 9, 10, 11 of the valley pseudospin by all-optical means. Despite successful realization of first opto-valleytronic devices 12, 13, controversy prevails with respect to possible intrinsic and extrinsic origins of significant variations in the degrees of valley polarization 9, 10, 11, 14, 15 and valley coherence 16, 17 observed for different material representatives of monolayer TMD semiconductors.

The valley pseudospin in TMDs is most conveniently accessed with photoluminescence (PL) polarimetry of band-edge excitons 5. The degrees of circular and linear PL polarization, $P_{c}$ and $P_{l}$, defined as the ratio of emission intensities $P=(I_{co}-I_{cr})/(I_{co}+I_{cr})$ detected in co-polarized ($I_{co}$) and cross-polarized ($I_{cr}$) configurations with a circularly ($\sigma$) or linearly ($\pi$) polarized excitation laser, are direct measures of valley polarization 9, 10, 11 and valley coherence 16, 17. Most values reported for circular and linear PL polarizations in monolayer TMDs are well below unity, and they vary significantly with material quality or the underlying substrate 9, 10, 11, 14, 15, 18, 19, 20. A detailed understanding of the variations in $P_{c}$ and $P_{l}$ has been elusive to date and partly attributed to different ratios of the exciton and valley dynamics in different samples. Analogous to optical spin orientation 21, ideal initialization of the valley polarization yields in a steady-state PL measurement $P_{c}=1/(1+\tau_{0}/\tau_{l})$ close to unity if the exciton lifetime $\tau_{0}$ is short as compared to the longitudinal valley decay time $\tau_{l}$. The same argument holds for $P_{l}$ as a measure of valley coherence with the transverse valley decay time $\tau_{t}$. For TMD monolayers with long-lived excitons, one therefore expects a sizable reduction in the degrees of circular and linear polarizations due to exchange-mediated valley decay and dephasing 22, 23, 24 active during the exciton lifetime.

In the following we demonstrate that the notion of the degrees of circular and linear polarization as being determined simply by the exciton and valley lifetimes is of limited validity, and therefore previous quantitative conclusions drawn from this simplistic picture should be critically revised. However, the model is helpful for a qualitative interpretation of varying degrees of valley polarization and coherence across single TMD monolayers. A more quantitative analysis can be provided by taking into account the optical valley initialization processes. This insight is based on our experiments on extended MoS2 flakes grown by chemical vapor deposition (CVD) and transferred onto SiO2/Si substrates as a representative TMD material system (see the Supplementary Information for sample details). The valley pseudospin physics in single- and poly-crystalline MoS2 monolayers were addressed with steady-state two-dimensional circular and linear PL polarimetry.

Photoluminescence spectroscopy of extended monolayers

We begin our studies by characterizing individual MoS2 crystals with cryogenic PL. The experimental setup for confocal PL spectroscopy and raster-scan imaging is shown schematically in Fig. 1a. The sample was cooled to $3.1$ K in a closed-cycle cryostat and positioned within the confocal excitation and detection spots of a low-temperature apochromatic objective. The polarization of the excitation and detection pathways was set independently for circular or linear polarimetry by according combinations of linear polarizers, half- and quarter-waveplates. A superconducting solenoid was used to apply magnetic fields up to $9$ T perpendicular to the sample.

A representative cryogenic spectrum of monolayer MoS2 recorded with a non-resonant excitation laser at $532$ nm is shown in Fig. 1b. It features $A$ and $B$ exciton PL around $1.9$ eV and $2.0$ eV characteristic of monolayer MoS2 on SiO2 6, 7. Additionally, we observed a red-shifted PL of localized excitons ($L$) which exhibit saturation as a function of the excitation laser power 25 in contrast to $A$ excitons with linear response (Fig. 1c). The maps in Fig. 1d - g, constructed by color-coding normalized raster-scan PL within exciton-specific bands indicated in Fig. 1b, show the $A$ and $L$ exciton PL profiles for two representative MoS2 flakes of our studies. The triangular shape (Fig. 1d and e) is typical for CVD-grown single-crystal flakes, while the star-shaped geometry of the poly-crystalline flake (Fig. 1f and g) reflects a cluster of single-crystal domains separated by grain boundaries 26, 27.

The PL intensity profiles show spatial variations across the flakes due to the presence of spatial inhomogeneities in the crystal quality. The edges and the center of the single-crystal triangle exhibit a more intense $A$ exciton PL than the rest of the flake (Fig. 1d), while the $L$ exciton emission is most intense in a ’puddle’ at the right edge of the triangle (Fig. 1e). Based on previous studies 26, 27, Raman spectroscopy (see the Supplementary Information) and optical inspection of other flakes on our sample, we identify the central region as a bilayer triangle formed on top of the monolayer single-crystal MoS2. The puddle is likely a collection of point defect contaminants responsible for exciton localization 25. The poly-crystalline star shows homogeneous domains of both $A$ and $L$ emission separated by lines of enhanced intensity (Fig. 1f and g) at the grain boundaries 26, 27 which also seem to favor exciton localization. At low excitation powers, emission hot-spots of cryogenic quantum dots 28, 29, 30, 31 dominate the PL intensity profiles of $L$ excitons (see the Supplementary Information).

Raster-scan opto-valleytronic imaging

The polarization-resolving feature of our setup allows us to identify characteristic signatures of crystal defects with circular and linear polarimetry. The two-dimensional maps of the degrees of circular and linear PL polarizations in Fig. 2a - d were recorded for the flakes of Fig. 1 with a laser at $637$ nm in resonance with the blue shoulder of the $A$ exciton PL which typically results in high near-resonant polarization values 10, 16, 32, 15. In the absence of an external magnetic field, time-reversal symmetry implies identical degrees of $K$ and $K^{\prime}$ valley polarization and thus equal $P_{c}$ recorded with $\sigma^{+}$ and $\sigma^{-}$ polarimetry 4. Equivalently, the valley coherence is independent of the choice of the linear basis 16 (see the Supplementary Information for the equivalence of complementary bases at zero magnetic field). First, we focus our analysis on the fundamental $A$ exciton to discuss the polarization in the $K$ valley, and the valley coherence in terms of $K$-$K^{\prime}$ superpositions.

The feature of highest contrast in the $P_{c}$ maps of the single- and poly-crystalline flakes is located at the triangle center (Fig. 2a) and coincides with the bilayer region discussed earlier. Optical selection rules are different for TMD bilayers 33 and provide direct means to single out bilayer regions in polarimetric imaging. Away from bilayer regions and edges, the degree of $P_{c}$ in Fig. 2a and b is remarkably high and homogeneous. Most surprisingly, the grain boundaries of the poly-crystalline star are almost invisible in the map of Fig. 2b. This observation of homogeneous valley polarization is in stark contrast to the sizable variations in the PL intensity profiles (Fig. 1d and f). It is also strongly contrasted by the pronounced inhomogeneities of the respective $P_{l}$ profiles in Fig. 2c and d, where both the puddle and the grain boundaries appear most prominently as sites of low valley coherence. To discuss these differences qualitatively, we recall the simple model where valley depolarization and dephasing are responsible for the reduction of $P_{c}$ and $P_{l}$ during the exciton lifetime. High $P_{c}$ values in the maps of Fig. 2a and b imply that the degree of optical valley pseudospin initialization is high in our experiments – a necessary condition for the optical generation of valley coherence. With the exception of the bilayer region, the homogenous $P_{c}$ profiles also suggest that the longitudinal valley relaxation is slow on the timescale of the exciton lifetime even in the presence of disorder. On the same timescale, however, rapid transverse valley relaxation results in locally reduced $P_{l}$ values as in the defective regions of the flakes in Fig. 2c and d.

In order to visualize the regions of rapid valley dephasing signified by large differences in $P_{c}$ and $P_{l}$, we present differential polarization profiles for both representative flakes. Fig. 2e and f were computed with a scaling factor $\eta$ of $0.65$ and $0.91$, respectively, to enhance the map contrast and thus the visibility of monolayer defects by highlighting bilayers in green and sites of valley decoherence (crystal edges, grain boundaries and the puddle of surface contaminants) in orange. Orange regions of disorder are characterized by rapid dephasing of optically generated quantum coherent superpositions of $K$ and $K^{\prime}$ excitons without signatures of low valley polarization. The white areas of the differential profiles in Fig. 2e and f identify the least defective monolayer crystals most favorable for opto-valleytronic applications.

Opto-valleytronic imaging in magnetic field

Our analysis is complemented by circular polarimetric profiling in non-zero magnetic fields. A positive magnetic field oriented orthogonal to the TMD monolayer plane breaks the time-reversal symmetry and lifts the valley degeneracy by decreasing (increasing) the exciton energy in the $K$ ($K^{\prime}$) valley 34, 35, 36, 37, 38, 39. Moreover, it changes the degree of circular polarization 34, 35, 36, 37, 38, which we monitored for our flakes with magneto-optical polarimetry. The polarimetric maps recorded in a magnetic field of $+9$ T within the $A$ and $L$ exciton bands are shown in Fig. 3a - d for the representative single-crystal triangle (see the Supplementary Information for complementary measurements in a magnetic field of $-9$ T). The main signature of the valley Zeeman effect in opto-valleytronic imaging is the field-induced difference in $P_{c}$ measured with $\sigma^{+}$ and $\sigma^{-}$ excitations (compare Fig. 3a and c with b and d). In a positive field of $9$ T the $P_{c}$ within both $A$ and $L$ spectral bands exhibits a decrease (increase) under $\sigma^{+}$ ($\sigma^{-}$) excitation in regular regions of the triangular flake.

To analyze the field-induced changes in the valley polarization of $A$ excitons we removed the contribution of the trion PL from the $A$ exciton band by spectral deconvolution (see the Supplementary Information). The resulting evolution of $P_{c}$ with magnetic field on a regular site of the single-crystal triangle is presented in the upper and central panels of Fig. 3e for the deconvolved $A$ exciton and for the $L$ exciton band, respectively. The data in the upper panel of Fig. 3e exhibit an evolution similar to $A$ excitons in WSe2 monolayer 36 with an $X$-shaped pattern that can be much less pronounced on different samples and under less resonant excitation 38. This $X$-shape is not a generic feature of TMD magnetoluminescence as it is distinct from the $P_{c}$ evolutions of MoSe2 neutral excitons in magnetic field37, 38. It rather reflects a counter-intuitive steady-state distribution of excitons among the two Zeeman-split $K$ and $K^{\prime}$ valleys 36: at positive magnetic fields the population in the energetically higher valley Zeeman branch $K^{\prime}$ becomes progressively protected from the relaxation into the lower branch $K$ with increasing magnetic field. This ’hot’ valley exciton population on a regular position is contrasted by the ’thermal’ distribution of $L$ excitons in the puddle (Fig. 3e, lower panel) with a sign reversal of the $X$-pattern as one would expect for a population redistribution that favors the exciton state of lowest energy.

The scheme of non-thermal, optically induced valley population imbalance is the opto-valleytronic counterpart of optical spin orientation in conventional semiconductors 21. For the fundamental $A$ exciton, rate-equation analysis 36 suggests that out-of-equilibrium $K$ and $K^{\prime}$ valley populations result from both a finite branching of the optically excited valley polarization into valley flipping and conserving relaxation channels and a slow valley depolarization in the exciton ground state. In fact, valley depolarization times much longer than the exciton decay time underpin the non-thermal valley population regime 36. Our analysis that accounts for inter-valley thermalization in addition to finite branching of the photoexcited population yields best fits to the data (solid lines in Fig. 3e) with $\tau_{l}/\tau_{0}\simeq 17$ for $A$ and $50$ for $L$ and $L_{P}$ excitons (see the Supplementary Information for model details and other relevant fit parameters). In contrast, the model with ideal initial polarization of the luminescent exciton states predicts $\tau_{l}/\tau_{0}\simeq 9$ for $A$ and $0.1$ for $L$ excitons, respectively. The discrepancies stem from the ignorance of the idealized model to the details of the valley pseudospin initialization in $K$ and $K^{\prime}$ exciton ground states. Actually it holds only if the valley conserving relaxation significantly outcompetes the valley flipping relaxation during the formation process of the $K$ and $K^{\prime}$ excitons out of the selectively excited valley as in the case of $A$ excitons in our experiment with near-resonant excitation. This does not apply to localized states, where significant population branching occurs upon relaxation, and different branching scenarios determine the non-thermal and thermal populations of $L$ and $L_{P}$ excitons: for the former, alike for $A$ excitons, the valley polarization is increasingly accumulated in the photoexcited upper valley due to a reduction of branching with magnetic field, while for the the latter defect-assisted relaxation likely renders the branching field-independent.

The regime of $\tau_{l}/\tau_{0}\gg 1$ identified by our analysis for both $A$ and $L$ excitons in monolayer MoS2 yields upper bounds on the longitudinal valley depolarization times of $\tau_{l}\simeq 80$ ps for $A$ and $230$ ps for $L$ and $L_{P}$ excitons if we take the PL decay time of $\tau_{0}=4.5$ ps for both the $A$ excitons in MoS2 15 and the dominant decay timescale of localized excitons in our sample (see the Supplementary Information for time-resolved PL data). In absolute terms, the longitudinal depolarization time of $A$ excitons might be shorter if scaled to a more rapid exciton decay potentially present in our CVD-grown MoS2 flakes but inaccessible in our experiments due to limited temporal resolution. In relative terms, however, our finding of $\tau_{l}\gg\tau_{0}$ is in accord with previous results on WSe2 36. At the same time it is contrasted by the theoretical estimate $\tau_{l}\simeq\tau_{0}\simeq 1$ ps for the valley depolarization dynamics in the presence of long-range exchange 23, 24 which in turn is an integral part of the finite-branching model 36. This caveat is qualitatively resolved by the notion of exchange-mediated valley depolarization timescales beyond $1$ ps for luminescent excitons with small center-of-mass momentum away from the light-cone edges 23, 24. For quantitative consistency, however, efforts in theory and experiment on the details of exciton valley dynamics in the presence of valley conserving and flipping relaxation channels and external magnetic fields are required beyond the scope of this work.

Wide-field opto-valleytronic imaging

We return to polarimetric mapping to demonstrate that the main valleytronic signatures of layered TMDs discussed above can be obtained in direct two-dimensional imaging. To this end we defocussed the excitation laser to illuminate a spot of $\sim 100~{}\mu$m diameter and replaced the single-mode fiber in the detection path with an imaging lens and a CCD. The resulting optical system images a sample area of $\sim(0.1\times 0.1)~{}\mu$m2 onto a single pixel of a standard room-temperature CCD array. A tunable band-pass filter placed before the CCD imaging lens was used to select exciton-specific bands, and polarimetric imaging was performed either in the circular or linear basis.

Polarimetric images of the two representatives flakes recorded in cross-linear configuration (orthogonal excitation and detection polarizations) within the $A$ exciton band are shown in Fig. 4a and b. In the given configuration, the high intensity features correspond to bilayer and dephasing regions of the flakes. Apart from a scaling factor, the images are equivalent to gray-scale versions of the false-color maps in Fig. 2e and f. They demonstrate qualitatively that disorder in monolayer TMDs can be visualized directly with a rather simple and efficient imaging technique. A quantitative analysis of the spatial distribution of $P_{l}$ can be carried out to yield maps analogous to Fig. 2c and d (see the Supplementary Information).

It is also straight forward to apply the technique in the circular basis. The images of Fig. 4c and d show the PL intensity of $A$ excitons in $\sigma^{-}$ co-circular and cross-circular polarimetry for a magnetic field of $+9$ T. The normalized difference of the images yields the $P_{c}$ profile equivalent to Fig. 3b (see the Supplementary Information). However, the non-thermalized $A$ exciton valley population of the upper Zeeman branch $K^{\prime}$ addressed by $\sigma^{-}$ excitation can be readily deduced from Fig. 4c and d; for a fully thermalized valley population the flake in Fig. 4c would be less intense than in Fig. 4d. The analogous set of measurements for the $L$ exciton band in $\sigma^{+}$ co- and cross-polarized configurations (Fig. 4e and f) illustrates the negative degree of circular polarization: the cross-polarized image of the flake is brighter (except for the puddle and a few hot-spots of point-like quantum dot emission) in accord with the raster-scan polarization image of Fig. 3c.

Conclusions

Our study identifies both raster-scan and wide-field polarimetric imaging as viable tools to explore the valley pseudospin physics in layered TMD semiconductors. Experiments that do not require the full spectral information but are meaningful within a limited PL bandwidth will benefit from a decrease of the measurement time required to achieve the same signal-to-noise performance in wide-field and hyperspectral raster-scan polarimetry. For $(20\times 20)~{}\mu$m2 sample areas like in Fig. 4 we estimate a speed-up by at least three orders of magnitude for wide-field imaging as compared to raster-scanning (see the Supplementary Information) which enables rapid large-scale monitoring of the sample (see the Supplementary Video). Clearly, our polarimetric studies are not limited to MoS2 crystals but establish novel analytic means for the entire class of layered TMD semiconductors and heterostructures, and can be extended to spin-polarization imaging of conventional semiconductor quantum wells 21. It will perform most efficiently when applied on TMD materials with spectrally separated exciton and trion emission bands, a condition currently accessible with MoS2 only by chemical treatment 40. By virtue of simplicity and efficiency our technique will facilitate the analysis of TMD materials in terms of crystalline and environmental disorder and, combined with a microscopic model for the optical valley initialization, it will enable quantitative access to site-dependent valley pseudospin dynamics. Moreover, our findings encourage the experimental exploration of topological exciton-polaritons by identifying sufficiently large regions of low valley depolarization and decoherence in monolayer TMDs for the realization of periodic exciton potential arrays that are key to the implementation of topolaritonic devices 41.

1 Sample fabrication and basic characterization

The chemical vapor deposition (CVD) of monolayer MoS2 crystals was adopted from Najmaei and coworkers 1. Briefly, SiO2/Si substrates with MoO3 seeding particles were placed in a quartz furnace in the presence of sulfur powder. The furnace was heated to $\sim 900$ ∘C for $15$ min with a flow of inert gases (N2 and Ar) under atmospheric pressure. After cooling down to room temperature, synthesized MoS2 crystals were transferred onto p-doped SiO2/Si substrates with a conventional polymethyl methacrylate (PMMA) transfer method 2. The PMMA was removed by rinsing the resulting sample with MoS2 on SiO2/Si in an acetone bath for three cycles of $15$ min. Optical and atomic force microscopy (AFM) images of typical MoS2 monolayer flakes after transfer onto SiO2/Si are shown in Supplementary Fig. 1.

2 Experimental setup and settings

Cryogenic confocal spectroscopy, raster-scan and wide-field imaging were performed in an ultra-low vibration closed-cycle cryostat (attocube systems, attoDRY1000) with a base temperature of $3.1$ K and a superconducting magnet with fields of up to $\pm 9$ T. The sample was positioned with nanopositioners (attocube systems, ANP101 series) into the focal plane of a low-temperature apochromat with a numerical aperture of $0.65$ (attocube systems, LT-APO/VIS/0.65) and confocal excitation and detection spots with $0.7~{}\mu$m full-width at half-maximum (FWHM) diameters. Linear polarizers (Thorlabs, LPVIS and LPVISB), half- and quarter-waveplates (B. Halle, RAC 3 series) mounted on piezo-rotators (attocube systems, ANR240) were used to control the photon polarization in the excitation and detection pathways. PL spectroscopy was performed with continuous wave (cw) excitation lasers at $532$ nm (CNI, MLL-III-532-50-1) or $637$ nm (New Focus, Velocity TLB-6704), and a standard spectrometer (PI, Acton SP-2558) with a liquid nitrogen cooled silicon CCD (PI, Spec-10:100BR/LN, gain setting of $4$ e-/count). The spectral resolution of the system was $\sim 0.35$ meV for confocal PL and $\sim 0.6$ cm*-1* for Raman measurements. A streak camera (Hamamatsu, C5680-24S) and a femtosecond (fs) optical parametric oscillator (OPO) operated at $630$ nm or $637$ nm (Coherent, Mira-OPO) were used for time-resolved PL studies. The fiber-coupled system exhibited dispersion-limited temporal resolution down to $\sim 10$ ps depending on the spectral bandwidth and the length of the single-mode fiber. For wide-field imaging, the OPO was tuned to $637$ nm and defocussed to illuminate a sample area of $\sim 4.4\cdot 10^{-3}$ mm2 at an average excitation power density of $\sim 200$ W/cm2. PL images with $50\times$ magnification were acquired in the spectral bands of tunable band-pass filters (Semrock, VersaChrome) with a silicon CCD (Point Grey, GS3-U3-14S5M-C, pixel size of $\left(6.45\times 6.45\right)~{}\mu$m2 and saturation capacity of $17$ ke-). All wide-field images were recorded with an integration time of $500$ s ($300$ s in Supplementary Fig. 10b) and a CCD gain setting of [math] dB (except for $5$ dB in Fig. 4a of the main text and $3$ dB in Supplementary Fig. 10a - d). The images digitized to a $16$ bit format with $14$ bit conversion were acquired with linear gamma encoding (gamma encoding of $0.7$ in Fig. 4 of the main text). The intensities were multiplied with scaling factors ranging from $1.0$ to $1.7$ to span the full dynamic range of the grayscale (identical factors were used for co- and cross-polarized images of the same exciton bands).

3 Spectroscopy and polarimetry experiments

3.1 Basic spectral characteristics

The PL of monolayer MoS2 crystals on SiO2/Si was excited either non-resonantly at $532$ nm ($2.33$ eV) or at $637$ nm ($1.95$ eV) in resonance with the blue shoulder of the $A$ exciton. A non-resonant spectrum recorded on a regular position of the triangular flake is shown in Supplementary Fig. 2a (same as in Fig. 1b of the main text). The non-resonant PL exhibits contributions of both neutral and charged excitons, $A$ and $A^{-}$ due to unintentional doping of MoS2 crystals on p-doped SiO2/Si substrates 3. The asymmetric lineshape of the total PL (spectrum in Supplementary Fig. 2a) can be decomposed into two Lorentzians with equal FWHM linewidths of $\sim 62$ meV separated by the trion binding energy of $\sim 30$ meV 4. The near-resonant spectrum recorded with a long-pass filter at $652$ nm ($1.90$ eV) in Supplementary Fig. 2b (gray and blue traces show for the same flake position non-resonant and near-resonant spectra, respectively) exhibited additional sharp features identified as Raman scattered photons by the substrate (Si) and the MoS2 monolayer, with assignments given in Supplementary Fig. 2c adopted from Ref. 5. The Raman spectra were useful to identify bilayer regions (as in the center of the single-crystal triangle) where they exhibited a characteristic splitting of the $A_{1g}(\Gamma)$ mode shown in Supplementary Fig. 2d under resonant excitation with the fundamental exciton 6. To avoid contamination of the opto-valleytronic properties of $A$ excitons by Raman photons, we selected a spectral band of $3$ meV width centered at $1.879$ eV near the PL maximum of the $A$ exciton emission and away from Raman resonances (the red spectral interval in Supplementary Fig. 2a - c is the same as in Fig. 1b of the main text ). Moreover, for the quantitative analysis of the degree of circular polarization of $A$ excitons presented in Fig. 3e of the main text, the contribution of $A^{-}$ was removed by fitting the total PL with two Lorentzians and subtracting the contribution of the trion PL from the spectral band chosen for the opto-valleytronic analysis as described above.

3.2 Time-resolved photoluminescence measurements

The PL decay dynamics of the monolayer MoS2 triangle from the main text are presented in Supplementary Fig. 3. The $A$ exciton PL, measured in a spectral band of $35$ meV, exhibited radiative decay with a time constant below the resolution limit of $\sim 14$ ps given by the dispersion-limited instrument response function (Supplementary Fig. 3a). Localized excitons, selected by a bandpass filter of $50$ meV width, exhibited two decay timescales (Supplementary Fig. 3b): the dominant fast component was resolution-limited (with partial contribution from $A$ excitons), and the slow component exhibited monoexponential decay with a lifetime of $174\pm 12$ ps in agreement with previous results 7. The slow decay channel accounted for $8$% of the total PL intensity on a regular position ($L$) and for $12$% in the puddle ($L_{P}$) of the triangular crystal.

3.3 Spectral characteristics of polarization-resolved photoluminescence

Supplementary Fig. 4a and c show co-polarized ($I_{co}$) and cross-polarized ($I_{cr}$) PL spectra recorded with an excitation laser at $637$ nm ($1.95$ eV) in circular and linear bases, respectively. The spectral characteristics of $P_{c}$ and $P_{l}$, respectively shown in Supplementary Fig. 4b and d, were calculated as the normalized differences between co- and cross-polarized PL intensities according to $P=(I_{co}-I_{cr})/(I_{co}+I_{cr})$. Data in Supplementary Fig. 4e and f confirm that the degree of linear polarization is independent of the choice of the linear basis (the linear PL polarization is parallel to the axis of the excitation laser set along $\pi^{H}$ and $\pi^{D}$ in Supplementary Fig. 4e and f, respectively). The red bar in Supplementary Fig. 4a - d was used for confocal opto-valleytronic imaging of $A$ excitons (same band as in Supplementary Fig. 2a - c and Fig. 1b of the main text).

3.4 Circular and linear raster-scan polarimetry

To construct polarimetric maps we raster-scanned the sample with respect to fixed confocal excitation and detection spots, and performed spectral acquisition of co- and cross-polarized PL at each raster pixel with an excitation power density of $\sim 2.5\cdot 10^{4}$ W/cm2. An averaged background spectrum was subtracted from all image pixels, and co- and cross-polarized PL intensities of $A$ and $L$ excitons were integrated within the spectral bands shown in Fig. 1b of the main text. Pixels with a standard deviation of $P$ above 0.05 (stemming from vanishingly small PL intensities away from the flake) were set to zero. In Supplementary Fig. 5a - d we reproduce the maps of Fig. 2a - d of the main text to indicate specific positions (numbered from $1$ to $6$) for a quantitative comparison of site-to-site variations of $P_{c}$ (blue bars) and $P_{l}$ (purple bars) summarized in Supplementary Fig. 5e for the of $A$ exciton band. Similar polarimetric measurements were performed on other MoS2 flakes. Supplementary Fig. 6a, c and Supplementary Fig. 6b, d show profiles of $P_{c}$ and $P_{l}$ for a poly-crystalline and a single-crystalline MoS2 monolayer, respectively. The bar chart in Supplementary Fig. 6e summarizes the statistics of $P_{c}$ and $P_{l}$ acquired on five different monolayer flakes.

3.5 Equivalence of polarization bases at zero magnetic field

In the absence of an external magnetic field, time-reversal symmetry implies the equivalence of the degrees of circular polarization measured in $\sigma^{+}$ and $\sigma^{-}$ configurations. The set of data in Supplementary Fig. 7a - c demonstrates that the degree of the circular polarization is independent of the choice of the circular basis: $P_{c}$ is identical (within the precision of our measurement) for $\sigma^{+}$ and $\sigma^{-}$ polarization bases (compare Supplementary Fig. 7a and b, respectively). This fact is also reflected by the vanishing difference profile $\Delta P_{c}=P_{c}(\sigma^{+})-P_{c}(\sigma^{-})$ shown in Supplementary Fig. 7c. The analogous set of data is respectively shown for the degrees of linear polarization recorded under horizontal ($\pi^{H}$) and diagonal ($\pi^{D}$) linearly polarized excitations in Supplementary Fig. 7d and e; their vanishing difference $\Delta P_{l}=P_{l}(\pi^{H})-P_{l}(\pi^{D})$ is plotted in Supplementary Fig. 7f. The data demonstrate the independence of $P_{l}$ of the choice of the linear basis.

3.6 Magnetic field dependence of the degree of circular polarization

A finite external magnetic field applied in Faraday geometry (perpendicular to the TMD monolayer surface) lifts the valley degeneracy and gives rise to a Zeeman shift of opposite sign for excitons in the $K$ and $K^{\prime}$ valleys. A positive magnetic field decreases (increases) the energy of the $K$ ($K^{\prime}$) valley exciton. Consequently, an applied magnetic field changes the degrees of circular polarization of the valley exciton emission. The changes are discussed in the main text for spectrally selected $A$ and $L$ excitons under $\sigma^{+}$ and $\sigma^{-}$ excitation at $+9$ T (the corresponding data of Fig. 3 are reproduced in Supplementary Fig. 8a, b, d, and e). In Supplementary Fig. 8c and f we plot for both types of excitons the change in the degree of circular polarization, quantified as $\Delta P_{c}=P_{c}({\sigma^{+}})-P_{c}({\sigma^{-}})$. For $A$ excitons at $B=+9$ T, $\Delta P_{c}$ is negative (magenta-colored) throughout the MoS2 flake (Supplementary Fig. 8c). This corresponds to a decrease (increase) of the $\sigma^{+}$ ($\sigma^{-}$) degree of circular polarization for the lower (upper) $K$ ($K^{\prime}$) Zeeman branch of $A$ excitons. The $L$ excitons show the same trend away from the puddle and reversed features in the puddle (cyan-colored region in Supplementary Fig. 8f) in accord with thermal population distribution discussed in the main text. The right panel of Supplementary Fig. 8 shows the same set of data but for $B=-9$ T. The sign reversal of the magnetic field results in interchanged roles of $\sigma^{+}$ and $\sigma^{-}$ polarizations (Supplementary Fig. 8g, h, j, and k) and thus in a sign reversal of $\Delta P_{c}$ (Supplementary Fig. 8i and l).

3.7 Wide-field imaging and polarimetry

As described previously the confocal setup was modified for wide-field imaging by defocusing the excitation laser to illuminate a sample area of $\sim 4.4\cdot 10^{-3}$ mm2 and by replacing the single-mode fiber in the detection path with an imaging lens and a silicon CCD. A tunable band-pass filter (with $15$ nm FWHM bandwidth) was inserted before the imaging lens to select exciton-specific bands. The corresponding setup schematics are shown in Supplementary Fig. 9a. A fs-OPO operated at $637$ nm was used to excite the PL at an average power density of $\sim 200$ W/cm2. This average power density was a factor of $\sim 150$ and $\sim 500$ lower than in confocal measurements with the excitation laser at $637$ nm (cw power density of $\sim 2.5\cdot 10^{4}$ W/cm2) and $532$ nm (cw power density of $\sim 9.4\cdot 10^{4}$ W/cm2), respectively. In the following, we discuss the main implications of the reduced power density for the observations in wide-field imaging and polarimetry.

To compare the confocal and wide-field imaging modes we recorded additional PL spectra and intensity maps in the confocal setup configuration at reduced power densities. In Supplementary Fig. 9b - g we present data recorded in the regime of low power excitation for the single-crystal MoS2 triangle discussed in the main text. A confocal PL spectrum recorded at an excitation power density comparable to that of the wide-field imaging mode is shown in Supplementary Fig. 9b. Colored bands indicate the spectral intervals that were used for the evaluation of confocal PL intensity maps in Supplementary Fig. 9f and g for $A$ and $L$ excitons, respectively. The bands were selected to match the setting of the band-pass filter in wide-field imaging.

At low excitation powers the PL intensity of $L$ excitons is not saturated (Supplementary Fig. 9c). This results in an increased intensity ratio of $L$ to $A$ exciton PL in the confocal spectrum of Supplementary Fig. 9b as compared to the data in Fig. 1b of the main text. Moreover, since the saturation responses of $L$ excitons away from the puddle and in the puddle are different (the respective data are denoted as $L$ and $L_{P}$ in Supplementary Fig. 9c), one expects a crossover in the relative PL intensities of defect-bound excitons in the puddle and on a regular position of the flake. This effect is observed in PL imaging at low excitation power densities: both wide-field (Supplementary Fig. 9e) and confocal (Supplementary Fig. 9g) images of the monolayer triangle show a more intense $L$ exciton PL away from the puddle, whereas at high excitation power densities the $L$ exciton PL in the puddle is more intense (Fig. 1e of the main text). The more pronounced appearance of hot-spots in the PL intensity profiles of $A$ excitons under low-power illumination (Supplementary Fig. 9d and f) as compared to the PL intensity map under high-power illumination (Fig. 1d of the main text) is also attributed to saturation effects.

As in confocal polarimetry, polarization-resolved wide-field imaging can be used to construct polarimetric profiles of valley polarization and valley coherence. To demonstrate the quantitative character of the technique, we show in Supplementary Fig. 10a - f the $P_{l}$ and $P_{c}$ profiles of the single-crystal MoS2 triangle of the main text obtained with wide-field polarimetry. Note that the opto-valleytronic profiles are in excellent quantitative agreement with confocal data shown in Fig. 2 and Fig. 3 of the main text if one takes into account the reduction of $P_{l}$ and $P_{c}$ due to wider spectral bands of the band-pass filter used in wide-field polarimetry.

The advantage of the wide-field imaging technique is that it reduces the integration time for opto-valleytronic profiling whenever full spectral information at each pixel of the map is not required. Ultimately, for an image size of $\sim\left(20\times 20\right)~{}\mu$m2 with a pixel size of $\sim\left(0.4\times 0.4\right)~{}\mu$m2 an acquisition acceleration of at least $\sim 2.5\cdot 10^{3}$ can be obtained with polarimetric profiling in wide-field imaging as compared to confocal raster-scanning with the same spatial resolution, spectral width, noise level, and transmission characteristics. For this estimation we neglected the confocal scanning time and assumed equal excitation power densities of the imaging techniques as well as shot-noise limited signals. We note that for our confocal raster-scan polarimetry the scanning time was a factor of two longer than the integration time, and the transmission properties were a factor of ten lower than for wide-field imaging.

4 Theoretical modeling

In the following we model the degree of circular polarization of excitons in monolayer MoS2 in the presence of an external magnetic field applied in Faraday geometry (the magnetic field is oriented perpendicular to the TMD monolayer surface). To this end we consider a four level system with the crystal ground state $\lvert 0\rangle$, the $K$ exciton state $\lvert 1\rangle$, the $K^{\prime}$ exciton state $\lvert 2\rangle$, and an excited state $\lvert 3\rangle$ in one of the two valleys. All levels are denoted in Supplementary Fig. 11, where without loss of generality the excited state $\lvert 3\rangle$ is in the $K$ valley. In the following, exciton states in $K$ ($K^{\prime}$) that couple to $\sigma^{+}$ ($\sigma^{-}$) polarized optical transitions are tagged with the valley index $\kappa=+1$ ($\kappa=-1$).

For a $\sigma^{+}$ ($\sigma^{-}$) polarized optical transition that is assumed to create population exclusively in the excited $K$ ($K^{\prime}$) state, the temporal evolution of the populations in each state is given by the following set of rate equations:

[TABLE]

where $\rho_{i}$ is the population of the $i$th state, $\dot{\rho_{i}}$ is the temporal derivative of $\rho_{i}$, and $\gamma_{ij}$ is the transition rate from state $i$ to $j$. The population dynamics are governed by the following rates: $\gamma_{03}$ denotes the absorption rate from the crystal ground state into an excited state in one of the two valleys, $\gamma_{31}$ and $\gamma_{32}$ are the valley conserving and valley flipping relaxation rates (or vice versa), $\gamma_{12}$ and $\gamma_{21}$ are the valley flipping rates among the $K$ and $K^{\prime}$ states, and $\gamma_{0}$ is the $K$ and $K^{\prime}$ exciton decay rate (including both radiative and non-radiative channels). The total population of the system is normalized to $\rho_{0}+\rho_{1}+\rho_{2}+\rho_{3}=1$, and steady-state solutions are obtained for $\dot{\rho_{i}}=0$.

The degree of circular polarization $P_{c}=(I_{co}-I_{cr})/(I_{co}+I_{cr})$ is obtained from steady-state populations $\rho_{1}$ and $\rho_{2}$ of the $K$ and $K^{\prime}$ valleys:

[TABLE]

with the explicit expression equivalent to the one derived in the Supplementary Information of Ref. 8:

[TABLE]

We rewrite this expression as:

[TABLE]

with a branching parameter $b_{\kappa}$ given by the ratio of the valley flipping to the valley conserving relaxation rates:

[TABLE]

The definition of the branching parameter implies $b_{\kappa}=\gamma_{32}/\gamma_{31}$ $\left(b_{\kappa}=\gamma_{31}/\gamma_{32}\right)$ for $\kappa=+1$ ($\kappa=-1$).

First, we examine the the degree of circular polarization given by Eqs. 3 and 4 at zero magnetic field. In the presence of time-reversal symmetry, the valley flipping processes are symmetric and thus $\gamma_{12}=\gamma_{21}$. Moreover, in the limit of ideal initial polarization, i. e. for $b_{\kappa}=0$, Eqs. 3 and 4 simplify to:

[TABLE]

where $r_{0}=\gamma_{12}/\gamma_{0}=\gamma_{21}/\gamma_{0}$ is the ratio of the zero-field longitudinal valley depolarization rate $\gamma_{12}=\gamma_{21}=\gamma_{l}$ to the exciton decay rate $\gamma_{0}$. This expression is equivalent to the one derived in Ref. 9, and it corresponds to the steady-state degree of circular polarization in optical spin orientation with ideal initial polarization 10:

[TABLE]

if we identify $\tau_{0}=1/\gamma_{0}$ as the exciton lifetime and $\tau_{l}=1/(2\gamma_{l})$ as the valley depolarization time. This limit where $r_{0}$ and thus $\tau_{0}/\tau_{l}$ are fixed by the degree of circular polarization $P^{0}_{c}$ was used in the main text to discuss the data at zero magnetic field.

If the initial polarization of the fundamental exciton populations in $K$ and $K^{\prime}$ valleys is non-ideal because of valley flipping events upon relaxation from the optically excited state, we have analogous to imperfect spin orientation 10:

[TABLE]

with the initial polarization $P_{i}^{0}=(1-b_{\kappa})/(1+b_{\kappa})$ given by the yield of optical valley polarization of $K$ and $K^{\prime}$ populations at zero magnetic field in the presence of finite branching $b_{\kappa}>0$.

In a finite magnetic field $B$ applied perpendicular to the sample, the valley Zeeman splitting $\Delta_{\mathrm{v}}$ 11, 12, 8, 13, 14, 15 introduces an imbalance between $\gamma_{12}$ and $\gamma_{21}$. In the presence of non-zero branching, we return to Eqs. 3 and 4 to recapitulate the findings of Ref. 8 with respect to the magnetic field evolution of $P_{c}$. First we note qualitatively that due to the valley Zeeman splitting, the imbalance of $\gamma_{12}$ and $\gamma_{21}$ increases (decreases) the population transfer rate from the upper (lower) to the lower (upper) Zeeman branch. Thus, provided sufficiently rapid valley depolarization on the timescale of the exciton lifetime, we expect the population to relax into the energetically lower valley. Consequently, the $P_{c}$ of the upper (lower) Zeeman valley should decrease (increase) with magnetic field. The opposite trend is observed in our study of $A$ and $L$ excitons in monolayer MoS2 away from defects in agreement with the $A$ exciton response in monolayer WSe2 reported in Ref. 8. This counter-intuitive population distribution was attributed to (i) slow inter-valley scattering and (ii) polarization protecting branching ratio $b_{\kappa}$ that decreases (increases) for the upper (lower) valley with magnetic field 8. Since no explicit expression for the functional form of the evolution of $b_{\kappa}$ with magnetic field was given in Ref. 8, and to account for different materials and experimental conditions, we approximate $b_{\kappa}$ with a linear function in $B$:

[TABLE]

where $b_{0}$ is the zero-field branching ratio and $\delta$ is a proportionality factor. It is worth noting that Eq. 9 should be interpreted as a low-field approximation since in it would yield unphysical negative values for $b_{\kappa}(B)$ in the limit of sufficiently high fields. Moreover, a difference in $\delta$ for the $K$ and $K^{\prime}$ valleys could account for broken time-reversal symmetry in absence of a magnetic field as pointed out in Ref. 8.

To proceed with the analysis of our data using the model of Ref. 8 where $\gamma_{l}\ll\gamma_{0}$ (or correspondingly $r_{0}\simeq 0$) was assumed, we fit the $P_{c}$ evolutions of the $A$ and $L$ excitons on a regular position of the triangular flake with $r_{0}$, $b_{0}$ and $\delta$ as fitting parameters. The results of best fits to the data, obtained from least $\chi^{2}_{\mathrm{red}}$ deviation (defined by the unweighted sum of squared deviations between fits and data and divided by the number of degrees of freedom), are presented in Supplementary Fig. 12.

The non-thermal population distribution of both $A$ and $L$ excitons is reflected by the counter-intuitive $X$-pattern of higher $P_{c}$ values for the Zeeman branches with higher energy. It is consistent with slow longitudinal valley depolarization on the timescale of the exciton lifetime, and best fits actually yield $r_{0}=0$ for both $A$ and $L$ excitons. The zero-field branching ratio for the selectively initialized population is small for $A$ excitons ($b_{0}=0.07$) and large for $L$ excitons ($b_{0}=0.90$) as expected from the less resonant excitation of the latter. In this framework, the evolution of the branching ratio with magnetic field critically determines the respective $P_{c}$ pattern: both $A$ and $L$ excitons exhibit an increasing protection of the optical valley polarization in the upper Zeeman valley as the branching ratio between the valley flipping and conserving relaxation rates decreases ($\delta>0$ for both $A$ and $L$) due to exchange-modified exciton dispersions 8.

The model of Ref. 8 within the framework of Eqs. 4 and 9 requires vanishingly small values of $r_{0}$ without quantifying its smallness. To obtain an estimate for the ratio of the longitudinal valley depolarization time to the exciton lifetime, we extend the model of Ref. 8 by assigning thermal imbalances to the longitudinal inter-valley flipping rates $\gamma_{12}$ and $\gamma_{21}$ explicitly using the Boltzmann factor $\exp(\beta)$:

[TABLE]

where, at a temperature $T$, the exponent $\beta=(\Delta_{\mathrm{v}}/2)/(k_{\mathrm{B}}T)$ is determined by the Boltzmann constant $k_{\mathrm{B}}$ and the valley Zeeman splitting $\Delta_{\mathrm{v}}=g\mu_{\mathrm{B}}B$ linear in the magnetic field $B$ (with $g$ being the effective exciton $g$-factor, and $\mu_{\mathrm{B}}$ the Bohr magneton) 11, 12, 8, 13, 14, 15. The zero-field ratio $r_{0}$, defined by the Eq. 6 above, is likely to exhibit a magnetic field dependence but short of knowledge of its functional form is assumed here as constant. The exponential imbalance of $\gamma_{12}$ and $\gamma_{21}$ is responsible for an effective unidirectional thermalization of population from the upper to the lower Zeeman branch during the exciton lifetime, and modifies Eq. 4 to:

[TABLE]

The degree of circular polarization is now sensitive to the thermal imbalance of the inter-valley scattering rates through the Boltzmann factor, and it accounts for the non-zero branching in the polarization initialization with its characteristic magnetic field dependence via $b_{\kappa}(B)$ given by Eq. 9. The second term drives the thermal population distribution among the $K$ and $K^{\prime}$ valleys, while the first term counteracts this thermalization via the functional dependence of the branching ratio $b_{\kappa}$ on the magnetic field in favor of a ’hot’ population distribution.

With this expression for the field-evolution of $P_{c}$, we performed a simultaneous fit of our set of data for $A$ and $L$ excitons away from the puddle and $L_{P}$ excitons. To minimize the number of free fitting parameters, we assumed no field-dependence of $b_{\kappa}$ for the case of $L_{P}$ excitons that exhibit a thermal population distribution, and set both the $g$-factors and the ratios $r_{0}$ equal for both $L$ and $L_{P}$. Moreover, we fixed the $A$ exciton $g$-factor to $4.0$ 14, 15, and assumed a global temperature $T$ for the entire set of data. With these parameters fixed, best simultaneous fit was obtained with $b_{0}=0.03$, $0.90$, and $0.82$ for $A$, $L$, and $L_{P}$; $r_{0}=0.03$ for $A$ and $0.01$ for $L$ and $L_{P}$; $\delta=0.012$ T*-1* for $A$ and $0.027$ T*-1* for $L$; $g=4.4$ for $L$ and $L_{P}$; $T=9$ K. The corresponding evolutions of $P_{c}$ for $A$, $L$ and $L_{P}$ excitons are shown in Supplementary Fig. 13a, b and c, respectively. We note that the best fit yielded an effective temperature of $T=9$ K consistent with local laser heating, and $g=4.4$ for $L$ and $L_{P}$ excitons in agreement with elevated values reported for excitons localized in quantum dots 16, 17, 18, 19.

The main effect of the model of Eq. 11 as compared to the model of Eq. 4 is the increased field dependence of the $A$ exciton branching ($\delta$ increased by a factor of two from $0.006$ T*-1* to $0.012$ T*-1*), whereas the values of the zero-field branching ratio of $b_{0}$ for the valley flipping to the valley conserving rate are comparable for the two models. Most importantly, the model yields an estimate for $r_{0}$ for both $A$ and localized excitons: bet fit values of $r_{0}=0.03$ for $A$, and $0.01$ for $L$ and $L_{P}$ correspond to $\tau_{l}/\tau_{0}\simeq 17$ and $50$, respectively. Taking the PL decay time of $4.5$ ps for the $A$ excitons in MoS2 7 that we do not resolve in our time-correlated PL experiments, we thus obtain an estimate for the longitudinal valley depolarization time of $\tau_{l}\simeq 80$ ps for $A$ excitons. With the same scaling we estimate $\tau_{l}\simeq 230$ ps for localized excitons that exhibit PL decay dynamics dominated by the fast component below the resolution limit of our fiber-based setup with dispersion (as shown in Supplementary Fig. 3b, only $\sim 10\%$ of the total PL intensity of both $L$ and $L_{P}$ excitons contribute to the slow decay component). The main difference for $L$ and $L_{p}$ populations stems from different branching during relaxation: while the upper valley polarization of $L$ excitons is weakly yet increasingly protected with magnetic field (with a slope of $0.027$ T*-1*), the branching ratio of $L_{P}$ excitons is independent of magnetic field. Given similar PL decay dynamics in the $L$ and $L_{P}$ bands, and similar branching ratios at zero magnetic field, we conclude that hot and thermal valley populations of localized excitons arise from field-dependent and field-independent branching ratios, respectively.