Multiple spin-orbit excitons in α-RuCl3 from bulk to atomically thin layers

The van der Waals Kitaev magnet α-RuCl3 has recently garnered considerable attention due to its possible realization of topological spin liquids. Combining Raman spectroscopy with numerical calculations, we report here the thickness dependence of electronic structure and ensuing low-energy excitations for exfoliated α-RuCl3. We observe two pronounced peaks at A1 = 249 meV and A2 = 454 meV, which are assigned to single and double spin-orbit (SO) excitons, respectively. Our numerical calculations support this interpretation by reproducing their spectral energy and shape with the electronic parameters: SO coupling λ = 140 meV, Hund’s coupling JH = 350 meV, and on-site Coulomb interaction U = 2.35 eV. The multiple SO excitons persist down to a single layer, whereas their peaks shift slightly to lower energy. For frequencies below 350 cm−1, both a magnetic continuum and phonons show noticeable thickness dependence. These results demonstrate that a SO entangled jeff = 1/2 picture remains valid in a monolayer limit despite the presence of lattice distortions.


INTRODUCTION
Two-dimensional van der Waals (2D vdW) materials have offered an excellent platform to explore electronic, optical, and magnetic phenomena emergent in a dimensional crossover from bulk to the 2D limit as well as to enable interface engineering by forming heterostructures 1,2 . When exfoliated down to the monolayer limit, a particular class of 2D materials shows prominent layerdependent physical properties, ranging from an-indirect-todirect bandgap transition, charge density wave, superconductivity, and a paramagnetic-to-ferromagnetic transition [3][4][5][6] . On the other hand, another class of 2D vdW materials retains their bulk properties in the ultrathin limit. The 2D Kitaev candidate α-RuCl 3a graphene-like quantum magnet is a valuable addition to the family of 2D vdW magnetic materials, in which the edge-sharing RuCl 6 octahedra constitute RuCl 3 layers connected by weak interlayer vdW interactions [7][8][9][10] . However, little is known about the electronic stability, the nature of exchange interactions, and the magnetic transition temperature in 2D atomic layers.
The Kitaev honeycomb lattice is an analytically solvable spin model known to harbor quantum spin liquids and Majorana quasiparticles as a consequence of spin fractionalization 11 . The application of a magnetic field creates a topological mass term, thereby generating a topological spin liquid with non-Abelian anyonic excitations. Remarkably, non-Abelian anyons are key elements of the implementation of a topological quantum computer. This exciting prospect has been vigorously tested in the spin-orbit (SO) coupled Mott insulator α-RuCl 3 at both zero and finite magnetic fields.
The bulk α-RuCl 3 has a magnetically ordered ground state that adjoins a proximate Kitaev spin-liquid phase. Although the zigzag magnetic order occurs at T N ≈ 7 K, there are thermodynamic and spectroscopic signatures of emergent Majorana fermions at elevated temperatures and finite magnetic fields [12][13][14][15][16] . On applying an in-plane external field, the zigzag order is suppressed at B C ≈ 7 T [17][18][19] . Before entering a partial spin-polarized state, which is topologically trivial, the narrow intermediate-field phase appears at B = 7-10 T, in which the half-integer quantized thermal Hall conductance is observed 20 . This is taken as evidence for a fieldinduced topological spin liquid with chiral Majorana edge current. At the respective field range, quasiparticle excitations are governed by the coexistence of Majorana and magnon bound states 21,22 . Despite extensive works, however, it is still controversial whether the bulk α-RuCl 3 hosts the Majorana quasiparticles predicted in the pure Kitaev model system. In the pursuit of the pure Kitaev spin liquid, it has been proposed that interfacial effects in α-RuCl 3 /graphene render the original Mott insulator α-RuCl 3 electron-doped and that an interface strain enhances Kitaev interactions [23][24][25][26] . To verify the tunability of Kitaev magnetism in α-RuCl 3 /graphene or related heterostructures, electronic and magnetic properties should be first investigated in a single layer limit.
Previous studies on atomically thin or nanoflake α-RuCl 3 , mainly focusing on low-energy magnetic and phonon excitations, showed that fractionalized magnetic excitations remain robust in the 2D limit [27][28][29][30][31][32] . Besides, there exists a tantalizing indication of enhanced spin-phonon coupling and weak symmetry reduction, compared with the bulk sample. Overall, reduced dimensionality seems to exert a minor impact on the low-energy magnetism, while preserving a proximate Kitaev phase. However, no systematic investigations have been made for electronic structures and excitations 16,[33][34][35] . Given that a SO-entangled state in α-RuCl 3 is delicately conditioned by the interplay of spin-orbit coupling (SOC), Coulomb interaction, and trigonal crystal field, a layer-dependent electronic structure of atomically thin α-RuCl 3 should be elucidated for better understanding the vdW Kitaev magnet in its 2D limit.
In the present work, we combine Raman spectroscopy with exact diagonalization calculations to identify layer-dependent lowlying excitations in the mechanically exfoliated α-RuCl 3 . We observe single and double SO excitations at A1 = 249 and 1 Department of Physics, Chung-Ang University, Seoul, Republic of Korea. 2   = 454 meV, respectively, which persist down to a single layer while undergoing a slight softening. Our theoretical calculations reproduce the gross feature with the SOC λ = 140 meV, the Hund's coupling J H = 350 meV, and the on-site Coulomb interaction U = 2.35 eV. Our results attest that lattice distortions modify the electronic structure and magnetic interactions to some degree, yet a SO-entangled picture holds in a monolayer limit.

RESULTS AND DISCUSSION
Thickness dependence of the room-temperature Raman spectra We first examine the thickness dependence of the low-frequency Raman spectra measured at T = 300 K in the circular RL polarization, allowing for probing the E g symmetry channel. In an ideal Kitaev honeycomb model, a magnetic Raman signal is exclusively allowed in this channel. In addition, the RL configuration minimizes a contribution from stray light scattering. The thickness of the exfoliated α-RuCl 3 flakes, ranging from 0.8 nm (=1 L) to 50 nm, was determined by atomic force microscope (AFM) images and their height profiles (see Fig. 1a, b and Supplementary Fig. 3). As exhibited in Fig. 1c, we are able to resolve four E g phonon modes, which are consistent with factor group analysis for the D 3d space group. The asterisk at 220 cm −1 appears owing to symmetry reduction induced by layer stacking. Our spectra agree nicely with the previously reported data with respect to the spectral shape, energy, and symmetry 13,[28][29][30][31][32][33]36,37 . We note that all the presented Raman spectra are properly corrected using the previously reported optical constant 13,37 .
As the sample is exfoliated down to 0.8 nm (=monolayer), the Raman spectra show neither a drastic change nor apparent anomalies except for the enhancement of quasielastic scattering and the spectral broadening of phonon peaks. This indicates that the samples are hardly degraded. Before proceeding, we admit that prior Raman studies have addressed a thickness dependence of phonons and magnetic excitations at low energies [28][29][30][31][32][33] . Nonetheless, probing low-and high-energy excitations at the same spot of the sample surface is essential to obtain the consistent thickness dependence of magnetic and electronic excitations. For this reason, we repeat the analysis of the lowenergy spectra.
As shown in Fig. 1d, the low-energy spectrum is decomposed to a sum of two Gaussian functions and two Fano profiles Here, the asymmetry parameter q measures the strength of coupling between a phonon and a continuum and ε is the reduced energy. This fitting procedure nicely reproduces the data. The quasielastic scattering centered at zero energy, described by a Gaussian-like spectral function, stems * 10 20  c Thickness dependence of the Raman spectra measured at T = 300 K in RL polarization. Four major phonons are assigned to E g symmetry. The asterisk denotes a symmetry-forbidden mode. The spectra are vertically shifted for clarity. d Representative fit of the low-energy Raman spectrum to a sum of two Gaussian and two Fano resonances. The green shading denotes a background magnetic continuum, which is thermally damped. e Thickness dependence of the integrated intensities of the E g (2) Fano phonon (cyan balls), the quasielastic response (open triangles), and the thermally damped magnetic continuum (green squares). f, g Thickness dependence of the Fano asymmetry and linewidth for the E g (2) Fano mode. The solid colored lines represent guides to the eye.
from diffusive spin fluctuations 38 . The broad background (green shading) mainly arises from the continuum of thermally damped fractional and paramagnon excitations, as alluded by the E g (1) and E g (2) Fano resonances. As the sample thickness is decreased through t = 3 nm, the intensity of the magnetic continuum and quasielastic scattering is gradually decreased and then shows an upturn (see Fig. 1e). Overall, this trend stands out in the quasielastic scattering. In contrast, the Fano E g (2) mode shows a persistent reduction in its intensity, while the E g (4) phonon hardly varies with thickness (not shown here). Based on the observed anomalies around t=3 nm, we infer a detectable change of magnetic behaviors in its 2D limit. As seen in Fig. 1f, g, both the asymmetry parameter (1/|q|) and linewidth of the Fano resonance increase steeply on approaching a monolayer limit. This is opposite to the continuing suppression of the Fano E g (2) mode (see Fig. 1e). This thickness dependence points to the presence of an extra channel for phonon decay, which may be associated with stronger coupling to the magnetic continuum. In addition, as α-RuCl 3 is dominated by several anisotropic exchange interactions, it is far from clear to what extent reduced dimensionality affects Kitaev magnetism in atomically thin α-RuCl 3 .
Temperature dependence of the Raman spectra of bilayer α-RuCl 3 Next, we turn to a temperature-dependent Raman study of bilayer α-RuCl 3 . The data were recorded over a wide temperature range between 4 and 300 K on the slow warming process. As shown in Fig. 2a, the magnetic continuum evolves systematically into a quasielastic response with increasing temperature. Unlike the A 1g mode 32 , involving mainly the out-of-plane displacements of the Ru and Cl atoms, the E g phonons show no apparent anomalies through the monoclinic-to-rhombohedral transition around T = 150 K. In Fig. 2b, we compare the temperature dependence of the Fano asymmetry between the bulk and bilayer samples. For both samples, 1/|q| shows a common trend. To be more specific, on cooling down to T* = 130 K, 1/|q| is temperature-independent and then increases quasilinearly. In the Kitaev related system, the increasing 1/|q | is taken as a Raman spectroscopic signature of spin fractionalization, reflecting the enhanced coupling between the phonon and the fractionalized continuum 37,38 . In this vein, the onset temperature T*, giving an energy scale of the Kitaev exchange interaction J K , is associated with a crossover from a simple to a Kitaev paramagnetic state. The Fermi statistics of fractionalized excitations gives another spectroscopic diagnostic for the spin fractionalization. In the intermediate-energy window ω ≈ (0.5-1.25) J K (Kitaev exchange interaction J K = 7-15 meV) a magnetic Raman scattering process is dictated by the creation or annihilation of pairs of Majorana fermions [36][37][38][39] . As such, the integrated Raman intensity I mid (ω) obeys the asymptotic Both inelastic neutron and Raman scattering data disclose that anharmonic magnon excitations or non-Kitaev interactions prevail for energies below 6 meV (≈48 cm −1 ) 15,37 . In consideration of this, we choose the middle energy of ω = 8.4-16.1 meV (=70-130 cm −1 ) to single out fermionic excitations. Even in the selected energy range, a bosonic contribution (e.g., multimagnon scatterings) is not completely removed. Thus, the bosonic terms should also be taken into account to give a complete description of I mid (T). This phenomenology is applied to our bulk and bilayer samples. As plotted in Fig. 2c, I mid (T) of the bulk and bilayer samples provides evidence for fermionic excitations, namely, a slight upturn below T = 50 K. Compared with the bulk, the low-T upturn feature becomes less pronounced for the bilayer sample. This seems to indicate a small suppression of the fermionic excitations in the atomic limit. Considering no essential difference of the T-dependence of 1/|q| between the bulk and the bilayer α-RuCl 3 , however, we conclude that both the atomically thin and bulk samples possess a similar degree of the Kitaev magnetism.

Thickness dependence of single and double SO excitons
Having established the thickness dependence of the low-energy lattice and magnetic excitations, we turn to excitations at highenergy transfers measured up to 4500 cm −1 (= 558 meV), which are the main focus of the present study. Exhibited in Fig. 3a are the Raman spectra with decreasing thickness from bulk to 0.8 nm (= 1 L). At room temperature, we observe two prominent peaks A1 and A2, which are located at ω = 249 meV and 454 meV, respectively, in the bulk limit. The peak energies are largely consistent with previous studies: A1 = 231 ± 3 meV, A2 = 524 ± 10 meV, and A3 = 745 ± 10 meV extracted by resonant inelastic x-ray scattering (RIXS) 34 , A1 ≈ 300 meV, A2 ≈ 530 meV, and A3 ≈ 750 meV by optical spectroscopy, and A1 = 248 meV and A2 = 450 meV by Raman spectroscopy at T = 5 K 35 (see Table 1 for comparison). Several comments on the spectral features are in order.
First, the A3 excitation is not detectable by Raman spectroscopy (Fig. 3b for the Raman spectrum taken up to 1 eV). In the case that the A3 excitation is of the triple SO exciton, it is natural that the third-order scattering has a negligible scattering intensity. Second, we observe two sharp peaks at 142 and 151 meV at low temperatures, which correspond to the A 0 (≈ 145 meV) peak (see Supplementary Fig. 2). The origin of the A 0 peak was previously interpreted in terms of a SO exciton, which corresponds to chargeneutral and spin-dressed excitations between the SOC-split levels j eff = 1/2 and j eff = 3/2 33 . Our data indicate that the 142 and 151 meV peaks are the overtones of the 37 meV E g (4) mode. This four-phonon scattering may be facilitated by a sideband of electronic excitations (SO exciton) together with the structural phase transition 38 . Third, the Raman and RIXS data show the nearly identical position of the A1 peak. However, the A1 peak energy extracted using optical spectroscopy is substantially larger than that by the Raman and RIXS data. This difference is owing to the distinct selection rule between Raman and optical spectroscopy. To be specific, the optical spectroscopy detects the SO exciton (not carrying a dipole moment) through a phononassisted process (for example, 37 meV phonon mode). Consequently, the infrared peak is blue-shifted by the phonon energy relative to the Raman peak. Fourth, the peak A2 detected by Raman spectroscopy is red-shifted by~80 meV in comparison to the RIXS and optical data.
Previously, the A1-A3 peaks were erroneously assigned to transitions between SOC-split e g states 12,33 . According to quantum chemistry calculations, the dd excitations from j eff = 3/2 to j eff = 1/ 2 are expected at 195 meV and 234 meV with a trigonal splitting 40 . In this vein, the A1 peak is attributed to a SO exciton. This kind of excitonic quasiparticle has been observed in a SO assisted Mott insulator Sr 2 IrO 4 41 . In contrast, as compared in Table 1, two scenarios have been proposed for the origin of the A2 and A3 peaks. Warzanowski et al. 35 assigned these peaks as double and triple SO excitons. As the RIXS and optical spectroscopy show the identical energies of the A2 and A3 peaks, multiple excitations should be directly infrared-active. However, the A2 and A3 optical peaks obey distinct temperature dependence, suggesting that the A2 peak is a phonon-assisted transition. This is at odds with no energy shift of the A2 peak between the RIXS and optical spectroscopy data. Lebert et al. 34 tentatively interpreted these peaks as charge-transfer type excitations from Cl 3p to Ru 4d. Fig. 3 Multiple spin-orbit excitons and their thickness dependence. a Thickness dependence of spin-orbit exciton excitations measured at T = 300 K in RL polarization. The spectra are vertically shifted for clarity. The asterisks are single-, two, and three-phonon excitations arising from the Si substrate. b Room-temperature Raman spectrum of bulk α-RuCl 3 measured up to Raman shift of 1 eV. c Representative fit of the spin-orbit excitons to four Gaussian profiles A1, A'1, A2, and A'2. False-color map of the Raman intensity in the frequency-thickness plane. d Frequency of the peaks A1 and A2 as a function of thickness. e Thickness dependence of the peak intensities normalized to the value of t = 50 nm as well as of the intensity ratio of the peak A1 to A2. The dashed lines are guides to the eye.
As such, identifying the origin of these peaks is crucial in unraveling the intricate electronic structures in α-RuCl 3 .
With this in mind, we analyze the thickness dependence of the Raman-active A1 and A2 peaks. Both peaks are modeled using Gaussian line shapes, as marked by color shadings in Fig. 3c. The best fit reveals the additional A'1 and A'2 peaks that may be associated with either the splitting of j eff = 3/2 levels owing to trigonal distortion or the intrinsic spectral features of SO excitons, see Fig. 4a. The energy difference between the A1 and A1' peaks amounts to 34 meV, which may correspond to the trigonal splitting energy. This value falls in the 10-40 meV range reported in the prior studies 34,40,42 . Nonetheless, this should not be taken literally because the characteristic two-peak structure is not apparently resolved, hindering their unambiguous modeling.
Fitting results are plotted in Fig. 3d, e. The peak positions of A1 and A2 undergo a quasilinear red-shift with reducing thickness down to t=7 nm and then show a steep drop below 3 nm. As seen from Fig. 3c, the SO excitons undergo a qualitative change in their scattering intensity. Our numerical calculations enable us to identify the A1 and A2 peaks to single and double SO excitons, respectively (as discussed later). As the A1 peak energy is related to the SOC constant λ ≈ (2/3)•A1, the SOC constant is evaluated to λ = 166 meV in bulk. The softening rate of the A2 peak (dω A2 /dt = 0.395 meV nm −1 ) is slightly faster than twice that of the A1 peak (dω A1 /dt = 0.184 meV nm −1 ). Besides, the energy of the double SO exciton is smaller than twice that of the single exciton, namely, 2 • A1-A2 ≈ 44 meV. The energy renormalization of the A2 peak suggests a strong effect of many-body interactions. Here we note that the Raman A2 peak is lower in energy than the RIXS A2 peak. This may be because the double SO exciton evolves the entire Brillioun zone in the Raman scattering process, or excitonic excitations are mediated through different intermediate states between Raman and RIXS. Nonetheless, we cannot exclude the possibility that the RIXS A2 peak is of distinct origin. As clearly seen from Fig. 3c, e, with reducing thickness, the A2 peak is much rapidly reduced in intensity than the A1 peak. In contrast to the monotonic decrease of the A2 peak, the A1 peak becomes stronger below t = 7 nm. We stress that a similar trend is observed for the thickness dependence of the low-energy magnetic and lattice excitations, see Fig. 1e. Despite the enhanced softening and varying intensity of the SO excitons, we can see that the SO excitons still survive.
Taken together, we conclude that vdW interlayer interactions modify the electronic structure possibly owing to lattice distortions in a few nm thickness range, yet are not strong enough to demolish a Mott insulating j eff = 1/2 state even in the monolayer.
Theoretically calculated Raman spectra For a microscopic understanding of the Raman spectrum in the high-energy region, we performed numerical calculations based on the Hubbard model of t 2g orbitals. Figure 4a presents the theoretical Raman spectra in RL polarization calculated with the physical parameters listed in Table 2. The theoretical A1 and A2 peaks appear at 244 meV and 456 meV, respectively. Our numerical calculations of resonant Raman scattering reproduce the experimental spectrum of the bulk sample in its energy. However, they fail to describe the relative intensity between the A1 and A2 peaks. In our model, e g orbital contributions are omitted, which can play a role in the intermediate states in the resonant Raman process with laser photon energy (~2.33 eV). This limitation may be responsible for the wrong intensity ratio.
To identify the origin of the A1 and A2 peaks, we calculate the distributions of single and double SO-exciton states from the ground state (see Methods for details). As plotted in Fig. 4b, the single and double SO excitons are distributed with the main peaks at 244 meV (almost the same as the A1 peak) and 485 meV (29 meV higher than the A2 peak), respectively. The resemblance Table 1. Low-energy electronic excitations probed by different spectroscopies and their assignments. between the distributions of the single and double SO excitons and the Raman spectrum corroborates an excitonic origin of both the A1 and A2 peaks. However, a close look unveils that the Raman peaks differ from the SO-exciton peaks. This discrepancy gives a rationale for why the A2 peak is lower than twice the A1 peak. The creation of a double SO exciton leads to a renormalization of its peak energy owing to many-body effects. In addition, the theoretical Raman spectrum shows a shoulder feature at 176 meV. This energy is somewhat larger than the A 0 peak energy of 145 meV observed in the low-T Raman spectra 33 . Further, this peak is reminiscent of the RIXS peak at 410 meV in Na 2 IrO 3 43 . Its origin was interpreted as the exciton-like peak stemming from an unusual coupling between the SO exciton and electron-hole excitations 43 . In α-RuCl 3, the electronic bandgap of 1.1 eV is well separated from the SO excitons, and thus this process is improbable. Rather, the theoretical 176 meV peak may be related to the experimental A'1 shoulder feature. Noteworthy is that the theoretical A1 peak of 210 meV is evaluated from 3λ /2 with the bare SOC λ ¼ 140 meV. The calculated A1 peak is somewhat smaller than the experimental value of 249 meV. Many-body effects of electronic correlations and kinetics certainly enhance the bare SOC splitting in the lattice.
We next discuss the thickness dependence of the single and double SO excitons (see Fig. 3). To quantify the red-shift trend toward the 2D limit, we calculated the Raman A1 and A2 peaks as functions of the trigonal distortion, the Coulomb repulsion, the bare SOC strength, and the scaled hopping integrals. Here, the parameter Δ t represents the trigonal distortion, which is defined as the energy splitting of t 2g orbitals Δ t ¼ E a1g À E e π g . The calculation results are summarized in Fig. 4c, d as well as in supplementary Figs. 6 and 7. We find a monotonic decrease of the A1 and A2 peaks as the trigonal distortion Δ t , the bare SOC, or the hopping integrals decrease. In addition, we can identify the discontinuous change at the specific values as marked with gray dotted lines in Fig. 4c, d. Our six-site cluster calculations show that a magnetic phase transition occurs from a disordered to a magnetically ordered phase when the trigonal distortion and the hopping integrals become larger than the critical values corresponding to the gray dotted lines. As the six-site is insufficient to describe a thermodynamic transition, this feature is ascribed to finite-size effects.
When approaching the 2D limit, local lattice distortions are apt to take place, which can alter an electronic structure and the magnitude and nature of exchange interactions. In this light, we scrutinize the Δ t dependence of the A1 and A2 peaks exhibited in Fig. 4c. With decreasing Δ t j j, the energies of the A1 and A2 peaks decrease quasilinearly and the A2 peak becomes narrower regardless of the sign of Δ t . We stress that the reduced trigonal distortion captures the softening of the single and double SO excitons with decreasing thickness (see Fig. 3). The trigonal compression or elongation leads to the variation of the Ru-Cl-Ru bond angle, which in turn determines the strength of exchange  Table 2. The light green line refers to the measured Raman spectra of the bulk sample. b Distributions of single and double spin-orbit (SO) exciton states. Arrows in a, b indicate the theoretical positions of A1 and A2 peaks. c, d Evolution of Raman intensities and the A1 and A2 peaks as functions of trigonal distortion parameter Δ t (¼ E a1g À E e π g ) and scaled hopping integral strengths. The parameter set is detailed in Table 2. The dotted lines in c, d are twice the A1 peak. Horizontal gray dotted lines in c, d indicate the specific values of the trigonal distortion and scaled hopping integral strengths, which show a discontinuous variation of the peak positions. interactions. Quantum chemistry calculations show that the Kitaev interaction becomes a maximum at~94°, being close to the experimental value of 94.09 •40,44 . In this nearly optimal bonding geometry, the reduced trigonal distortion generically increases the strength of non-Kitaev terms, while suppressing the Kitaev interaction 45 . Furthermore, recent exact diagonalization calculations on the K-Γ model show an enhanced low-energy Raman response when the ratio of Γ/K increases 46 . This theoretical result is consistent with the observed increase of the low-energy magnetic excitations (see Fig. 1c). However, we cannot totally exclude the possibility that a concerted interplay of the structural distortions, bare SOC, hopping integrals, or Coulomb repulsion comes into play. The thickness dependence of the electronic bandgap may shed further light on understanding the variation of local electron distribution in atomically thin α-RuCl 3 .
A combined study of exact diagonalization calculations and Raman scattering measurements allows determining the parameters of the local electronic structure: the SOC constant λ = 140 meV, the effective Hubbard parameter U = 2.35 eV, the Hund's coupling J H = 350 meV, and the four hopping integrals for the bulk sample (see Supplementary Discussion). While being exfoliated to the monolayer α-RuCl 3 , the magnetic and SO excitonic excitations and the Fano asymmetry show a noticeable thickness dependence. The comparison between experimental data and calculations suggests that the red-shift of the SO excitons is associated with the reduction of a trigonal crystal field. This in turn enhances the low-energy magnetic excitation and the Fano asymmetry. However, the fermionic contribution to the Raman intensity shows little change (see Fig. 2c). A small increase of the non-Kitaev terms against the Kitaev interactions does not alter the Fermi statistics substantially 37 .
At last, we will discuss the implication of our results. α-RuCl 3 is a vdW Kitaev magnet known thus far in the family of 2D magnetic crystals. This layered material realizes a proximate Kitaev spin liquid, possessing exciting application perspectives in topological quantum computation. According to prior experimental and theoretical works, there is a further need for suppressing non-Kitaev interactions to stabilize the Majorana fermions and Z 2 fluxes. The exfoliated α-RuCl 3 on SiO2/Si shows an opposite trend. In this situation, the interfacial engineering of a few-layer α-RuCl 3 with graphene attracts our attention as the rule of thumb in controlling Kitaev magnetism. According to theoretical results 25,26 , a lattice strain of the α-RuCl 3 /graphene heterostructure places the atomically thin α-RuCl 3 in the Kitaev spin-liquid phase. At the same time, interface charge transfer may alter the trigonal distortion, hopping integrals, or Coulomb repulsion as well as induce a band renormalization, which in turn modifies the exchange parameters.
Given these multifaceted interfacial effects, experimental techniques that can characterize the electronic structure and magnetic excitations in atomically thin samples should be established to better control the physical properties of heterostructures. As demonstrated above, the multiple SO excitons in conjunction with the low-frequency magnetic continuum can be used as a road map for gauging Kitaev interaction, crystal structure, SOC, and Mottness. By taking full advantage of these spectroscopic indicators, we will be able to weigh the heterostructure constituents rationally for the material-level realization of vdW Kitaev magnets.
In summary, combining Raman spectroscopy with exact diagonalization calculations, we observe a reduced dimensionality effect of α-RuCl 3 on its magnetic, lattice, and electronic properties in the few nm thickness range. Further, we identify single and double SO excitons, allowing for measuring Kitaev magnetism as well as a local electronic structure. This lays a firm foundation for controlling a Kitaev spin liquid in a few-layer α-RuCl 3 through strain, current, and heterostructure.

Crystal growth and characterization
Single crystals of α-RuCl 3 were synthesized by a vacuum sublimation method. A commercial compound of RuCl 3 (Alfa Aesar) was thoroughly ground, and dried in a quartz tube under vacuum to fully dehydrate. The evacuated quartz ampoule was sealed and placed in a temperature gradient furnace. A powder of RuCl 3 was heated at 1080°C for 24 h and then slowly cooled down to 600°C at a rate of −2°C per hour after dwelling for 5 hours. We obtained millimeter-sized crystals with a shiny black surface. Their structural, thermodynamic, and spectroscopicnetic characterizations were extensively characterized in our previous works 15,17,38 . Our single crystals showed a magnetic ordering at T~6.2 K, as shown in Supplementary Information.

Raman scattering measurements
The α-RuCl 3 samples of various thicknesses were mechanically exfoliated from its bulk crystal using Scotch transparent tape (3 M) and then were transferred onto Si/SiO 2 substrates. The thickness of the exfoliated nanoflakes, determined by AFM, ranged from 0.8 nm to 50 nm. To minimize exposure to air, the exfoliated samples were mounted on a liquid-He-cooled continuous cryostat right after the exfoliation and then were cooled down to T ¼ 4 K.
For Raman scattering experiments, we used a single-grating spectrometer (Princeton Instruments, SP-2500i) with a focal length of 50 cm equipped with a liquid-nitrogen cooled CCD detector (Princeton Instruments, Spec-10). The Raman spectra were recorded in exact backscattering geometry by using a 532 nm laser as an excitation source. The laser power was kept to P = 300 µW with a spot size of~1 µm. Local heating of the sample did not exceed 3 K. We calibrated the Raman spectra using a standard Neon calibration lamb (Newport 6032). All of the experiments were carried out in a circularly polarized configuration to minimize the effect of low-frequency stray light. Besides, the penetration depth effects of Raman spectra were corrected using the optical constants.

Numerical calculations
We employed the Hubbard model of t 2g orbitals with periodic six-site clusters incorporating the SOC, Kanomori-type Coulomb interactions, and hopping integrals between nearest-neighboring Ru ions (see the detail in ref. 47 ). The physical parameters are listed in Table 2. The value of J H is adapted from the previous literature 48 . U is determined for the main optical peak position (U À 3J H ) to be consistent with the experimental value of 1.2 eV 12,47 . Hopping integrals are characterized by four parameters (t 1 ; t 2 ; t 3 ; t 4 ) (see the details in Supplementary Information and ref. 48 ).
To keep the C 3 rotation, they are determined by averaging three nearestneighboring results of the monoclinic α-RuCl 3 as in ref. 45 . We note that t 2 is chosen with a smaller value than −120.1 meV to fit the Raman spectra well (see Supplementary Information). To calculate the theoretical resonant Raman spectra, we exploited the Kramers-Heisenberg-Dirac equation of Stokes scattering 49 as following where Ψ g and E g are the ground state and its energy, Ψ f j i and E f are the excited state and its energy, ω in and δ in are the energy and broadening of an incident laser, j is the current operator, and ε and ε 0 are the polarization vectors of the incident and scattered photons. We set ω in ¼ 2:33 eV and δ in ¼ 0:5 eV.
To calculate the distribution of sinlge and double SO excitions, we consider sixfold t 2g orbitals that are split into twofold j eff ¼ 1 2 and fourfold j eff ¼ 3 2 orbitals. Let a y iσ and b y im be the creation operators of σð¼ ± 1 2 Þ state among j eff ¼ 1 2 orbitals and mð¼ ± 1 2 ; ± 3 2 Þ state among j eff ¼ 3 2 orbitals at the i site, respectively. When single SO exciton is created at the i site, unoccupied j eff ¼ 1 2 hole and occupied j eff ¼ 3 2 electron exchange each other. The operator b im jF i ihF i ja y iσ creates the single SO exciton at the i site. Here, F i j i is the state, in which six t 2g orbitals are fully occupied at the i site. Thus, the single SO-exciton state at the i site Ψ i mσ can be defined as jΨ i mσ i ¼ b im jF i ihF i ja y iσ jΨ g , where Ψ g is the ground state. Its distribution is calculated with the following relation where E g is the ground energy and δ is the broadening parameter (¼ 200  meV). Similarly, double SO-exciton state at the i and j sites jΨ ij mi σi mj σj i can be defined as jΨ ij mi σi mj σj i ¼ b imi jF i ihF i ja y iσi b jmj jF j ihF j ja y jσj jΨ g i. Its distribution is also calculated with the following relation