Introduction

The polaronic effect1, which describes the strong coupling between charge and lattice vibrations, has a key role in a broad class of novel quantum phenomena ranging from colossal magnetoresistance2 to anomalous photovoltaic effect3. In particular, the polaronic effect on excitons can profoundly modulate exciton dynamics upon photoexcitation and has been employed to describe intriguing optical and optoelectronic properties in materials such as hybrid organic-inorganic perovskite solar cells4,5,6. Compared with three-dimensional (3D) bulk systems, two-dimensional (2D) atomic crystals possess a couple of unique advantages in exploring the polaronic effect on exciton dynamics. First, the reduced dielectric screening in atomically thin samples enhances both the excitonic effect7 and the electron–phonon (e-ph) coupling8, which is expected to promote the polaronic effect of excitons. Second, unlike bulk materials in which the e-ph coupling is largely determined by intrinsic electronic and phonon band structures with limited tunability, 2D materials provide greater flexibility for engineering e-ph coupling through a number of approaches including carrier doping9,10 and interfacial coupling11,12,13,14,15, as well as dimensionality modulation8, and therefore hold high promise for the future development of optoelectronic devices.

The realization of a long-range magnetic order in 2D semiconducting CrI3 paves the way to engineer optical and optoelectronic properties of 2D magnetic semiconductors16,17,18,19,20,21,22,23. The large excitonic effect24 from the localized molecular orbitals, of neither Wannier-type in 2D TMDCs nor Frenkel-type in ionic crystals, can be considered as the microscopic origin of the giant magneto-optical Kerr effect25 and magnetic circular dichroism26 signals in 2D CrI3. Meanwhile, the strong e-ph coupling is suggested to cause the large Stokes shift, profound broadness, and skewed lineshape in the photoluminescence (PL) spectra of 2D CrI326. The coexistence of excitons and strong e-ph coupling in 2D CrI3 naturally leads to open experimental questions of whether polaronic character emerges in the exciton dynamics and whether they are affected by the long-range magnetic order.

One fingerprint for the polaronic effect is the development of phonon-dressed electronic bands that appear as satellite bands in proximity to the original undressed one. Such features manifest as multiple equally spaced replica bands in angle-resolved photoemission spectroscopy (ARPES)9,14,27,28,29,30,31,32 or as discrete absorption and emission lines in linear optical spectroscopy33. However, such signatures of the polaronic effect have not been revealed so far in 2D CrI3, as the sizable bandgap (~1.1 eV)26 and extreme surface sensitivity34 of CrI3 make ARPES measurements challenging, whereas the potential inhomogeneous broadening could largely smear out individual lines for phonon-dressed satellite bands in linear optical spectroscopy.

In this work, we exploit temperature and magnetic field-dependent resonant micro-Raman spectroscopy, to show the direct observation of the polaronic character of excitons in bilayer CrI3. The polaronic effect manifests in Raman spectra as a well-defined, periodic pattern of broad modes that is distinct from sharper phonon peaks. The profile of this periodic pattern and its temperature and magnetic field dependence further reveal essential information including the e-ph coupling strength and the tunability of polaronic effect by the magnetism in bilayer CrI3. We mainly focus on bilayer CrI3 because it features a single magnetic phase transition from the layered antiferromagnetic (AFM) to ferromagnetic (FM) order and briefly compare to the results on thicker CrI3 flakes afterwards.

Results

Excitonic transitions and strong electron–phonon coupling

We start by identifying excitonic transitions and e-ph coupling in bilayer CrI3 using temperature-dependent PL and linear absorption spectroscopy. Bilayer CrI3 was fully encapsulated between few-layer hexagonal BN (hBN) and placed on a sapphire substrate (for details, see “Methods”). Linear absorption spectroscopy measurements were then performed in a transmission geometry (see “Methods”). Figure 1 shows representative PL and absorbance spectra taken at 80 K, 40 K, and 10 K that correspond to well above, slightly below, and well below the magnetic critical temperature TC = 45 K, respectively25,26,34. A single PL mode at 1.11 eV and three prominent absorbance peaks at 1.51 eV, 1.96 eV, and 2.68 eV (denoted as A, B, and C, respectively) are observed across the entire temperature range. These three energies are in good agreement with the ligand-field electronic transitions assigned by differential reflectance measurements on monolayer CrI326 and bulk CrI335,36 and have been later revealed to be bright exciton states through sophisticated first principle GW and Bethe-Salpeter equation calculations24. The large Stokes shift (~400 meV) between the PL and A exciton absorption peak is consistent with previous report6 and indicates strong electron–phonon coupling in 2D CrI3. Although the absorbance spectra show little temperature dependence except for the appearance of a weak shoulder at 1.79 eV at 10 K (orange arrow), the PL spectra are clearly temperature dependent. In particular, the temperature dependence of the PL full width at half maximum, \({\Gamma}(T)\), is well fitted by the model functional form, \({\Gamma}\left( T \right) = {\Gamma}_0 + \frac{\gamma }{{\exp \left( {\frac{{\hbar \omega _{{\mathrm{LO}}}}}{{k_BT}}} \right) - 1}}\), with the first term for temperature-independent inhomogeneous broadening and the second term for homogeneous broadening from the exciton coupling with a longitudinal optical (LO) phonon at frequency \(\omega _{{\mathrm{LO}}}\). Taking \(\omega _{{\mathrm{LO}}} = 120.6\) cm−1 found later on in Fig. 2, we obtain \({\Gamma}_0 = 163.9 \pm 2.7\) meV and \(\gamma = 164.2 \pm 8.1\) meV, which suggests that the broadness of the exciton modes arises from both inhomogeneous broadening from disorders and homogeneous broadening from e-ph coupling. The large homogeneous broadening parameter (\(\gamma\)) indicates strong vibronic modes mixing in the PL spectra, which precludes the formation of well-resolved phonon sidebands6.

Fig. 1: Exciton transitions in bilayer CrI3.
figure 1

Photoluminescence (PL) (left, blue) and absorption (right, orange) spectra of a bilayer CrI3 encapsulated between two hBN flakes and placed on the sapphire substrate, at 10, 40, and 80 K. A, B, and C denote three main exciton transitions at 1.51, 1.96, and 2.68 eV, and the orange arrow marks a shoulder mode at 1.79 eV appearing only at low temperatures. Spectra at 10 K and 40 K are offset vertically for clarity. Inset shows the fitted full width at half maximum (FWHM) of PL spectra as a function of temperature, \({\Gamma}\left( T \right)\), (diamond symbols) and its fitting to the functional form \({\Gamma}\left( T \right) = {\Gamma}_0 + \frac{\gamma }{{\exp \left( {\frac{{\hbar \omega _{{\mathrm{LO}}}}}{{k_BT}}} \right) - 1}}\) with the first and second terms for impurity-related inhomogeneous broadening and e-ph coupling-induced homogeneous broadening, respectively. Error bars indicate one standard error in fitting the FWHM of PL spectra.

Full size image
Fig. 2: Polaronic character of exciton dynamics in bilayer CrI3.
figure 2

a Raman spectrum of bilayer CrI3 acquired in the linearly crossed polarization channel at 40 K using a 633 nm laser. Insets show 1st-order single-phonon modes (linear parallel/crossed, i.e., XX/XY, channel in blue/red), 2nd-order two-phonon modes, and 3rd order three-phonon modes. The red arrow indicates the Raman-inactive longitudinal optical (LO) phonon frequency of importance, \(\omega _{{\mathrm{LO}}}\). The spectral intensities in the frequency range above 118 cm−1 are scaled by a factor of 0.1. b Zoom-in of the orange shaded area in the spectrum in a. The solid orange line overlaid on the raw data is the fit to the periodic pattern in the Raman spectrum with a sum of eight Lorentzian profiles indexed from N = 1 to N = 8 and a constant background, i.e., \(\mathop {\sum }\nolimits_{N = 1}^8 \frac{{A_N\left( {\frac{{{\Gamma}_N}}{2}} \right)^2}}{{\left( {\omega - \omega _N} \right)^2 + \left( {\frac{{{\Gamma}_N}}{2}} \right)^2}} + C\). c Fit of the Raman spectrum over the 1st-order spectral range without (upper panel) and with (lower panel) the consideration of the N = 1 broad mode (shaded Lorentzian profile in the lower panel). d Plot of the fitted central frequency (\(\omega _N\)) of the Nth Lorentzian profile in data taken at 40 K (solid orange triangle) and 290 K (hollow blue triangle). Dashed line is a linear fit (\(\omega _N = N\omega _{{\mathrm{LO}}}\)) to the plot that yields a slot of \(\omega _{{\mathrm{LO}}} = 120.6 \pm 0.9\) cm−1. Inset shows the atomic displacement of the LO phonon mode.

Full size image

Polaronic character in Raman spectra

We next proceed to perform resonant micro-Raman spectroscopy measurements with an incident wavelength of 633 nm matching the energy of the B exciton on an encapsulated bilayer CrI3 flake placed on a SiO2/Si substrate (see “Methods”). Figure 2a displays a representative Raman spectrum acquired in the crossed linear polarization channel at 40 K (slightly below TC = 45 K). Note that this spectrum covers a much wider frequency range than earlier Raman studies on CrI334,37,38,39,40,41,42,43,44. The multiphonon scattering is visible up to the 3rd order, and their zoom-in Raman spectra are shown in the inset of Fig. 2a. The 1st-order single-phonon peaks appear in the relatively low frequency range of 50–150 cm−1, and are assigned to be of either Ag or Eg symmetries under the C3i point group (see Supplementary Note 1), which is consistent with earlier work34,37,38,39,40,41,42,43,44 and proves the high quality of our samples. The 2nd-order two-phonon and the 3rd-order three-phonon modes show up in slightly higher frequency ranges of 190–290 cm−1 and 310–410 cm−1, respectively, and show decreasing mode intensities at higher-order processes, same as typical multiphonon overtones under harmonic approximation45 or cascade model46. In addition to and distinct from these multiphonon features, we resolve a remarkable periodic modulation across a wide frequency range of 70–1100 cm−1 in the low intensity part of the Raman spectrum (highlighted by the orange shaded area in Fig. 2a). This low intensity periodic pattern consists of clean, individual Lorentzian profiles and survives up to the 8th order (Fig. 2b), well beyond the highest order (3rd order) of multiphonon overtones, and each order of it spans for ~50 cm−1 frequency range, much wider than the linewidth of any observed phonon modes (insets of Fig. 2a for phonons). Such a periodic pattern is also observed in the anti-Stoke’s side at higher temperatures in bilayer CrI3, for example, up to the 2nd order at 290 K (see Supplementary Note 2), which clearly supports its Raman origin instead of luminescence.

We fit this low intensity periodic pattern using a summation of Lorentzian profiles of the form \(\mathop {\sum }\nolimits_N \frac{{A_N\left( {\frac{{{\Gamma}_N}}{2}} \right)^2}}{{\left( {\omega - \omega _N} \right)^2 + \left( {\frac{{{\Gamma}_N}}{2}} \right)^2}} + C\) with central frequency \(\omega _N\), linewidth \({\Gamma}_N\), and peak intensity \(A_N\) of the Nth period and a constant background \(C\) (see fitting procedure in “Methods”). Among all eight orders (\(N = 1,\,2,\, \ldots ,8.\)) in Fig. 2b, the presence of the 1st-order broad mode is deliberately validated in Fig. 2c that fitting with this 1st-order broad mode (orange shaded broad peak in the bottom panel) is visibly better than without it (top panel). This improved fitting by involving the 1st-order broad mode is further rigorously confirmed by the bootstrap method47 (see Supplementary Note 3). Figure 2d shows a plot of the central frequency \({\upomega}_N\) as a function of the order N with data taken at 40 K (\(N = 1,\,2,\, \ldots ,\,8\)) and 290 K (\(N = - 2,\, - 1,\, \ldots ,\,3\)), from which a linear regression fit gives a periodicity of \(120.6 \pm 0.9\) cm−1 and an interception of 0 ± 0.2 cm−1. To the best of our knowledge, such a periodic pattern made of individual Lorentzian profiles previously has only been seen in multiphonon Raman spectra of Cd, Yb, and Eu monochalcogenides described by configuration-coordinate model48,49,50,51,52,53,54,55. However, the periodic pattern observed in bilayer CrI3 here differs from these monochalcogenide multiphonon modes, as the broad linewidth of 1st-order mode contradicts with the sharp 1st-order forbidden LO phonon in Cd and Yb monochalcogenides48,49,50,51 and the persistence (or even enhancement) of higher-order multiphonon below TC = 45 K is in stark contrast to the disappearance of paramagnetic spin disorder-induced multiphonon below magnetic phase transitions in Eu monochalcogenides52,53,54,55. Because no known multiphonon model can capture all characteristics of our observed periodic pattern as well as the broad linewidths of each mode, we are inspired to consider the electronic origin. Indeed, strikingly similar features have been seen in polaron systems through the energy dispersion curves (EDCs) of ARPES9,14,27,28,29,30,31,32 and linear absorption and PL spectroscopy5,33,56. In those cases, the periodic patterns in their energy spectra arise from the phonon-dressed electronic state replicas, or sometimes also referred as phonon-Floquet states57, and the periodicity is given by the frequency of the coupled phonon. Owing to the high resemblance between the lineshapes of our Raman spectrum and those polaron energy spectra5,14,27,28,29,30,31,32,56, we propose that this periodic pattern in Raman spectra of 2D CrI3 stems from inelastic light scattering between the phonon-dressed electronic states caused by the polaronic character of B excitons in 2D CrI3, whereby the B exciton at 1.96 eV, with the electron(hole) in the weakly dispersive conduction (highly dispersive valence) band of Cr 3d (I 5p) orbital character24,58, couples strongly to a phonon at 120.6 cm−159. It is worth noting that a recent theoretical work predicts magnetic polaronic states in 2D CrI3 because of charge-magnetism coupling60, whereas our work suggests polaronic exciton states due to charge-lattice coupling.

We then proceed to identify the source and the character of the phonon at 120.6 cm−1. We first rule out the possibility of this phonon arising from either the hBN encapsulation layers or the SiO2/Si substrate, as a similar periodic pattern in the Raman spectrum is also observed in bare bulk CrI3 crystals (see Supplementary Note 4). Compared with the calculated phonon band dispersion of monolayer CrI361, we then propose the LO phonon calculated to be at ~115 cm−1 as a promising candidate, whose slight energy difference from the experimental value of 120.6 cm−1 could result from the omission of e-ph coupling in calculations. This LO phonon mode belongs to the parity-odd Eu symmetry of the C3i point group, and its atomic displacement field transforms like an in-plane electronic field (Ex, Ey) (see inset of Fig. 2d)61. Its odd parity makes it Raman-inactive and absent in the 1st-order phonon spectra (Fig. 2a inset, top panel), whereas its polar displacement field allows for its strong coupling to electrons/holes and prompts the polaronic character of the charge-transfer B exciton (see Supplementary Note 5 for measurements with additional laser wavelengths). In addition, this LO phonon band is nearly dispersionless and has a large density of states, further increasing its potential for coupling with the B exciton in 2D CrI3.

Temperature dependence of the polaronic effect

Given the coexistence of a 2D long-range ferromagnetic order and polaronic effect of excitons below TC = 45 K in bilayer CrI3, it is natural to explore the interplay between the two. For this, we have performed careful temperature-dependent Raman spectroscopy measurements and fitted the periodic pattern in every spectrum with a sum of Lorentzian profiles. Figure 3a displays the periodic pattern in Raman spectra taken at 70 K and 10 K, well above and below TC, respectively. Comparing these spectra, not only do more high-order replica bands become visible at lower temperatures (i.e., from \(N = 6\) at 70 K to \(N = 8\) at 10 K), but also the spectral weight shifts toward the higher-order bands (i.e., from \(N = 1\) at 70 K for the strongest mode to between \(N =\) 3 and 4 at 10 K in Fig. 3b). The appearance of higher-order modes at lower temperatures possibly results from a combination of the narrow exciton linewidth (~50 cm−1) and the dispersionless nature of coupled LO phonon. More importantly, the spectral weight distribution (AN vs. N) quantifies the e-ph coupling strength, and its spectral shift across TC confirms the interplay between the polaronic effect and the magnetic order in bilayer CrI3. Theoretically, the polaron system consisting of dispersionless LO phonons and charges is one of the few exactly solvable models in many-body physics62, and the calculated polaron spectra can be well-described by a Poisson distribution function, \(A_N = A_0\frac{{e^{ - \alpha }\alpha ^N}}{{N!}}\)63,64, where \(A_0\) is the peak intensity of the original electronic band, \(A_N\) is the peak intensity for the Nth replica band with \(N\) phonon(s) dressed, and \(\alpha\) is a constant related to the e-ph coupling in 3D (i.e., \(\alpha _{3{\mathrm{D}}}\)) that can be scaled by a factor of \(3{\uppi}/4\) for 2D (i.e., \(\alpha _{2{\mathrm{D}}}\))65. By fitting the extracted Lorentzian peak intensity profile at every temperature to the Poisson distribution function (see fits of 10 K and 70 K data in Fig. 3b), we achieve a comparable fitting quality to that for ARPES EDCs in polaron systems9,31 at every temperature and eventually arrive at the temperature dependence of \(\alpha _{2{\mathrm{D}}}\), which remains nearly constant until the system is cooled to TC and then increases by almost 50% at the lowest available temperature 10 K of our setup (Fig. 3c). In addition to the anomalous enhancement of \(\alpha _{2{\mathrm{D}}}\) across TC, the value of \(\alpha _{2{\mathrm{D}}} = 1.5\) at 10 K is the highest among known 2D polaron systems including graphene/BN heterostructures (\(\alpha _{2{\mathrm{D}}} = 0.9\))14 and bare SrTiO3 surfaces (\(\alpha _{2{\mathrm{D}}} = 1.1\))29.

Fig. 3: Enhanced e-ph coupling across the magnetic onset TC in bilayer CrI3.
figure 3

a Raman spectra of bilayer CrI3 acquired at 10 K and 70 K, respectively. Solid orange lines are fits to the raw Raman spectra, using a sum of N Lorentzian profiles and a constant background, i.e., \(\mathop {\sum }\nolimits_N \frac{{A_N\left( {\frac{{{\Gamma}_N}}{2}} \right)^2}}{{\left( {\omega - \omega _N} \right)^2 + \left( {\frac{{{\Gamma}_N}}{2}} \right)^2}} + C\). b Histogram plot of the fitted Lorentzian mode intensity (\(A_N\)) as a function of order index (N) at 10 K and 70 K. Solid curves are fits of the peak intensity profiles to the Poisson distribution function, \(A_N = A_0\frac{{e^{ - \alpha }\alpha ^N}}{{N!}}\). Error bars represent one standard deviation that convolutes the standard errors from the fittings to the local individual and global multiple Lorentzian profiles. c Plot of two-dimensional electron–phonon coupling constant (\(\alpha _{2{\mathrm{D}}}\)) as a function of temperature. The dashed vertical line marks the magnetic onset TC = 45 K, and the solid line is the fit to the functional form \(\alpha _{2{\mathrm{D}}}\left( T \right) = \left\{ {\begin{array}{*{20}{c}} {A\sqrt {T_{\mathrm{C}} - T} + B;\,\, T {\,}< {\,}T_{\mathrm{C}}} \\ {B;{\hskip 60pt}T {\,}> {\,}T_{\mathrm{C}}} \end{array}} \right.\). Error bars in c represent one standard error in the Poisson fitting.

Full size image

Magnetic field dependence of the polaronic effect

It has been shown that bilayer CrI3 transitions from a layered AFM to FM with increasing out-of-plane magnetic field (\(B_ \bot\)) above the critical value BC of 0.7 T17,21,25,26. We then finally explore the evolution of the polaronic effect across this magnetic phase transition by performing magnetic field-dependent Raman spectroscopy measurements. Here, we choose circularly polarized light to perform magnetic field-dependent measurements in order to eliminate any Faraday effect from the optical components that are situated in close proximity to the strong magnetic field. Figure 4a shows Raman spectra taken at \(B_ \bot\) = 0 T and \(\pm{\!}\)1 T, below and above BC, respectively, in both RR and LL channels, where RR(LL) stands for the polarization channel selecting the right-handed (left-handed) circular polarization for both incident and scattered light (see Supplementary Note 6). At 0 T, the spectra are identical in the RR and LL channels, consistent with zero net magnetization in the layered AFM state for bilayer CrI3 at \(\left| {B_ \bot } \right| {\,}< {\,}B_{\mathrm{C}}\). At \(\pm{\!}\)1 T, the spectra in the RR and LL channels show opposite relative intensities under opposite magnetic field directions, owing to the fact that the net magnetization in the FM state for bilayer CrI3 at \(\left| {B_ \bot } \right| {\,}> {\,}B_{\mathrm{C}}\) breaks the equivalence between the RR and LL channels. To better quantify the magnetic field dependence of the spectra, we measured Raman spectra in the RR and LL channels at \(B_ \bot\) from \(-\)1.4 T to 1.4 T every 0.1 T. We fit the spectrum at every magnetic field to extract \(A_N\) first and then \(A_0\) and \(\alpha _{2{\mathrm{D}}}\). Figure 4b shows that \(A_0\) has abrupt changes at \(B_ \bot = \pm\)0.7 T in both RR and LL channels, consistent with the first order magnetic phase transition at BC. Furthermore, the magnetic field dependence of \(A_0\) shows an opposite trend in the RR channel from that in the LL channel, whereas the sum of \(A_0\) from both channels remain nearly constant to the varying magnetic field. This observation can be understood by that, under a time-reversal operation, the RR channel transforms into the LL channel and the direction of the net magnetization at \(\left| {B_ \bot } \right| {\,}> {\,}B_{\mathrm{C}}\) flips, resulting in that the Raman spectrum in the RR channel at \(B_ \bot {\,}> {\,}0.7\) T is equivalent to the spectrum in the LL channel at \(B_ \bot {\,}< -{\!}0.7\) T. Figure 4c shows that \(\alpha _{2{\mathrm{D}}}\) is magnetic field independent, suggesting that the interlayer magnetic order barely affects the e-ph coupling strength and that the in-plane long-range magnetic order is responsible for the strong enhancement of e-ph coupling at TC. This finding corroborates with the fact that the 120.6 cm−1 phonon has in-plane atomic displacement.

Fig. 4: Evolution of polaronic effect across the magnetic phase transition from the layered AFM to FM in bilayer CrI3.
figure 4

a Raman spectra of bilayer CrI3 acquired at 10 K in the circularly parallel polarization channels, RR and LL, where RR (LL) stands for that both incident and scattered polarizations are selected to be right-handed (left-handed) circularly polarized, with an applied out-of-plane magnetic field (\(B_ \bot\)) of 1 T (top), 0 T (middle), and −1 T (bottom), respectively. b, c Plots of the Poisson fit amplitude \(A_0\) (b) and electron–phonon coupling strength \(\alpha _{2{\mathrm{D}}}\) (c) as a function of the applied \(B_ \bot\) in RR (orange data points) and LL (blue) channels. Solid lines are step (orange and blue in b) and linear (gray in c) function fits to the magnetic field dependence of \(A_0\) and \(\alpha _{2{\mathrm{D}}}\), respectively. Error bars are defined as one standard error of the fitting parameters.

Full size image

Discussions

Our further Raman spectroscopy studies on tri-layer, four-layer, and five-layer CrI3 show qualitatively same findings as those in bilayer CrI3 (see Supplementary Note 7) and again echoes with the in-plane nature of the 120.6 cm−1 Eu phonon and the intralayer charge-transfer B exciton. Our data and analysis reveal the phonon-dressed electronic states and suggest the polaronic character of excitons in 2D CrI3, which arises from the strong coupling between the lattice and charge degrees of freedom and is dramatically modified by the spin degree of freedom of CrI3. The exceptionally high number of phonon-dressed electronic state replicas (up to \(N = 8\)) further suggests 2D CrI3 as an outstanding platform to explore nontrivial phases out of phonon-Floquet engineering, whereas the significant coupling to the spin degree of freedom adds an extra flavor whose impact on the phonon-Floquet states has not been studied. For example, one can imagine creating topological states through the band inversion between the phonon-dressed replicas of CrI3 and the electronic state of a material in close proximity.

Methods

Sample fabrication

CrI3 single crystals were grown by the chemical vapor transport method, as detailed in ref. 40. Bilayer CrI3 samples were exfoliated in a nitrogen-filled glove box. Using a polymer-stamping transfer technique inside the glove box, bilayer and few-layer CrI3 flakes were sandwiched between two few-layer hBN flakes and transferred onto SiO2/Si substrates and sapphire substrates for Raman spectroscopy and PL/linear absorption spectroscopy measurements, respectively.

Linear absorption spectroscopy

A bilayer CrI3 sample on a sapphire substrate was mounted in a closed-cycle cryostat for the temperature-dependent absorption spectroscopy measurements. A broadband tungsten lamp was focused onto the sample via a 50× long working distance objective. The transmitted light was collected by another objective and coupled to a spectrometer with a spectral resolution of 0.2 nm. The absorption spectra were determined by \(1 - \frac{{I_{{\mathrm{sample}}}(\lambda )}}{{I_{{\mathrm{substrate}}}(\lambda )}}\), where \(I_{{\mathrm{sample}}}(\lambda )\) and \(I_{{\mathrm{substrate}}}(\lambda )\) were the transmitted intensity through the combination of sample and substrate and through the bare substrate, respectively.

PL spectroscopy

PL spectra were acquired from the same bilayer CrI3 sample where we carried out linear absorption measurements. The sample was excited by a linearly polarized 633 nm laser focused to a ~2 μm spot. A power of 30 μW was used, which corresponds to a similar fluence reported in the literature26 (10 µW over a 1 μm-diameter spot). Transmitted right-handed circularly polarized PL signal was dispersed by a 600 grooves/mm, 750 nm blaze grating, and detected by an InGaAs camera.

Raman spectroscopy

Resonant micro-Raman spectroscopy measurements were carried out using a 633 nm excitation laser for the data in the main text and 473 nm, 532 nm, and 785 nm excitation lasers for data in Supplementary Note 5. The incident beam was focused by a 40× objective down to ~3 μm in diameter at the sample site, and the power was kept at ~ 80 μW. The scattered light was collected by the objective in a backscattering geometry, then dispersed by a Horiba LabRAM HR Evolution Raman spectrometer, and finally detected by a thermoelectric cooled CCD camera. A closed-cycle helium cryostat is interfaced with the micro-Raman system for the temperature-dependent measurements. All thermal cycles were performed at a base pressure lower than 7 × 10−7 mbar. In addition, a cryogen-free magnet is integrated with the low temperature cryostat for the magnetic field-dependent measurements. In this experiment, the magnetic field was applied along the out-of-plane direction and covered a range of \(- 1.4\) to \(+\)1.4 Tesla. In order to avoid the Faraday rotation of linearly polarized light as it transmits through the objective under the stray magnetic field, we used circularly polarized light to perform the magnetic field-dependent Raman measurements.

Fitting procedure

For every sample, we have taken temperature and magnetic field-dependent Raman spectra on the hBN/SiO2/Si substrate with the same experimental conditions as that on the CrI3 flakes. The Raman spectra from the substrate, an extremely gradual background with a Si phonon peak at ~525 cm−1, shows no dependence on temperature (over the range of 10–70 K) and magnetic field (0–2.2 T). To fit the periodic oscillations in Raman spectra of CrI3 flakes, we follow the procedure described below. (I) we fit the Si phonon peak at ~525 cm−1 in both spectra taken on the CrI3 thin flake and the bare substrate to extract the Si peak intensity, \(I_{{\mathrm{Si}}}^{{\mathrm{sample}}}\) and \(I_{{\mathrm{Si}}}^{{\mathrm{substrate}}}\). (II) we multiply the background spectrum by a factor of \(\frac{{I_{{\mathrm{Si}}}^{{\mathrm{sample}}}}}{{I_{{\mathrm{Si}}}^{{\mathrm{substrate}}}}}\), which is ~1, and then subtract off the factored background from the raw Raman spectrum of sample. This process leads to the pure Raman signal for CrI3 whose baselines are nearly identical over the temperature range of interest (10–70 K). (III) we fit the sharp CrI3 phonon peaks with Lorentzian functions and subtract their fitted functions from the background free Raman spectrum from step (II). This leads to a clean spectrum with only periodic broad modes for a global fitting. (IV) we fit the clean spectrum from step (III) with a sum of multiple Lorentzian functions, \(\mathop {\sum }\nolimits_N \frac{{A_N\left( {\frac{{{\Gamma}_N}}{2}} \right)^2}}{{\left( {\omega - \omega _N} \right)^2 + \left( {\frac{{{\Gamma}_N}}{2}} \right)^2}} + C\). For the neatness of the data presentation in Figs. 2 and 3, we only show the fitted line from step IV in the plots.