First-principles study of the magneto-Raman effect in van der Waals layered magnets

Introduction
Magnetism in two-dimensional (2D) van der Waals (vdW) layered materials has attracted much attention due to the fascinating physics and potential applications1,2,3,4. CrI3 is among the first 2D magnets to be discovered that exhibit spontaneous out-of-plane ferromagnetism in the monolayer limit1. A characteristic feature of 2D materials is that the intralayer bonding is strong while the interlayer coupling is weak and typically dominated by vdW interactions. Such a contrast gives rise to an interesting phenomenon that the intralayer magnetic order cannot be easily altered, but the spin orientations between the layers can be relatively easily manipulated between ferromagnetic (FM) and antiferromagnetic (AFM) orders by various methods such as interlayer stacking sliding/rotation5,6 and applying an electric/magnetic field7,8,9. The manipulation of interlayer magnetism in vdW layered magnets results in rich physical properties, including novel layered-dependent magneto-optical effects1,2, electric-field induced Kerr rotation7,8, and giant second-order harmonic generation10.
The long-range FM order can break the time-reversal symmetry, resulting in magneto-optical Kerr and Faraday effects. The Kerr effect was the one used for detecting the spontaneous magnetization of monolayer CrI3 and few-layer Cr2Ge2Te6 in experiments by magneto-optical Kerr effect microscopy1,2. The Kerr and Faraday rotation angles are closely related to the spin-dependent off-diagonal terms in the polarizability tensor (or dielectric tensor) relative to the spin-independent diagonal terms1,11,12. Such a spin-dependent polarizability tensor can be obtained via first-principles density functional theory (DFT) calculations by considering the spin-orbit coupling (SOC) effect in the magnetic materials as demonstrated in several prior theoretical works13,14,15. However, when the total magnetization of the system (e.g., bilayer CrI3) is zero due to the interlayer AFM order, Kerr and Faraday rotations vanish since the off-diagonal terms in the polarizability tensor become zero1. This calls for another technique to probe the AFM state more effectively.
The magneto-Raman effect has been recently demonstrated to provide such a technique for studying both the FM and AFM orders in 2D magnets16,17,18,19. Raman intensity of a phonon mode is determined by its Raman tensor, which corresponds to the derivative of the polarizability tensor with respect to the phonon mode’s vibration20,21,22. Since the polarizability tensor has spin-dependent anti-symmetric off-diagonal terms1,11,12, the Raman tensor has spin-dependent off-diagonal terms as well9,23,24,25,26. More importantly, a system only has one polarizability tensor, but can possess many Raman tensors depending on the number of phonon modes in the system. Therefore, unlike the polarizability tensor, the spin-dependent off-diagonal terms in Raman tensors do not simply vanish when the total magnetization is zero, as Raman tensors of some phonon modes retain the anti-symmetric off-diagonal terms in the AFM state. This allows magneto-Raman spectroscopy to probe both the FM and AFM orders in 2D materials, as recently demonstrated by several experimental works on CrI3 and VI39,23,24,25,26. Furthermore, the magneto-Raman effect naturally captures the intricate spin-phonon coupling27. Finally, compared to scanning magneto-optical Kerr effect microscopy and spin-polarized scanning tunneling microscopy, magneto-Raman spectroscopy is fast, nondestructive, and relatively cheap. Considering the significance of magneto-Raman scattering, it is imperative to develop first-principles simulation capability (i.e., digital twin of magneto-Raman spectroscopy) and a corresponding analytical model that can reveal the underlying mechanism. Although first-principles calculations of magneto-Raman scattering were attempted for monolayer and bulk CrI328,29, it has not yet been demonstrated for bilayer and few-layer CrI3 (the focus of most experimental works), let alone other 2D magnets. In addition, to the best of our knowledge, no generalized model derived from the theory of Raman scattering has been developed.
In this work, we have carried out systematic investigations of the magneto-Raman effect in CrI3 from single-layer (1L) to four-layer (4L) by performing first-principles simulations and developing a generalized polarizability model that can be broadly applied to 2D magnets of any thickness. The calculated polarized Raman spectra of 1L CrI3 with intralayer FM order and 2L-4L CrI3 with both interlayer FM and AFM orders agree well with existing experimental data9,23,24,25, capturing how Raman selection rules are changed due to the magnetic order, including Raman activation of odd-parity phonon modes. The model successfully reveals how the change of spin orientation in each layer is coupled to the layer’s vibration pattern to induce or eliminate the anti-symmetric off-diagonal terms in the Raman tensor and subsequently to alter the Raman intensity of a phonon mode. In addition, our simulations and modeling demonstrate that a similar magneto-Raman effect should be present in MnBi2Te4, a novel 2D topological magnetic material30,31,32, but it has not yet been experimentally observed. The quantitative analysis provided by DFT calculations sheds light on the underlying reason and reveals how to observe it experimentally. Furthermore, we show that the magneto-Raman effect exists in CrCl3 as well, but the change of the magnetic order from the out-of-plane to in-plane direction alters the location of the anti-symmetric off-diagonal terms in the Raman tensor and thus the Raman optical selection rules. Finally, to the best of our knowledge, our work uncovers the spin dependence of E modes for the first time in 2D magnets and predicts how they can be studied experimentally. We demonstrate that the correlation between phonon modes and magnetic orders is much more broadly present than the sole A2 mode of CrI3 studied in experiments.
Results
Theory of resonant and magneto-Raman scattering
Raman intensity of a phonon mode j can be expressed as (Ipropto frac{{n}_{j}+1}{{omega }_{j}}| {{bf{e}}}_{{rm{s}}}cdot {bf{R}}cdot {{bf{e}}}_{{rm{i}}}^{T}{| }^{2}), where es and ei are the electric polarization vectors of the scattered and incident light, respectively, and R is the Raman tensor associated with the phonon mode. ωj is the frequency of the phonon mode j, and ({n}_{j}={({e}^{hslash {omega }_{j}/{k}_{B}T}-1)}^{-1}) is the Boltzmann distribution function at the given temperature T. According to the Placzek approximation20,21,22,33,34, the matrix element of the 3 × 3 Raman tensor of the phonon mode j is given by
in which V is the volume of the primitive cell, χαβ is the electric polarizability tensor, rl(μ) is the position along the l direction of the μ-th atom (l = x, y, or z), (frac{partial {chi }_{alpha beta }}{partial {r}_{l}(mu )}) is the derivative of the polarizability tensor over the atomic displacement along the l direction, ({e}_{l}^{j}(mu )) is the atomic displacement along the l direction of the μ-th atom for the jth phonon mode (i.e., the eigenvector of the dynamic matrix), and Mμ is the atomic mass. Basically, one needs to calculate the derivatives of the polarizability tensor with respect to the atomic displacements for obtaining the Raman tensor. For non-resonant Raman scattering, we often ignore the dependence of the polarizability tensor on the incident photon energy and only calculate the static polarizability (or dielectric) constant at the zero photon energy20,21. However, to capture the photon energy dependence of Raman intensity in resonant Raman scattering, one needs to compute the frequency-dependent polarizability tensor at different atomic displacements and thus obtain the Raman tensor at the laser excitation energy EL using the finite difference method22,33,34. Here, we computed complex χαβ(EL) for each atomic displacement using the independent particle approximation, where the imaginary part is determined by a summation over all valence/conduction bands at every k-point in the Brillouin zone that satisfies the energy conservation, and the real part is obtained by the usual Kramers-Kronig transformation35. Finally, the derivatives of the polarizability tensor, combined with the phonon frequencies and phonon eigenvectors, yield the Raman tensor matrix and, subsequently, the Raman intensity for every phonon mode under a given laser polarization set-up. In the common experimental back-scattering geometry (light z in and z out), the polarization vectors of the scattered and incident light are in-plane (i.e., x−y). The Raman intensity can be simplified when the polarizations are along the crystalline axes (denoted as XX, XY, YY, YZ, etc.): I(XX) ∝ ∣R11∣2 when es = ei = (1, 0, 0); I(XY) ∝ ∣R12∣2 when es = (1, 0, 0), ei = (0, 1, 0); I(YY) ∝ ∣R22∣2 when es = ei = (0, 1, 0); I(YZ) ∝ ∣R23∣2 when es = (0, 1, 0), ei = (0, 0, 1).
The general form of the Raman tensor is given by
where the values of Raman tensor elements are complex and depend on the laser excitation energy36,37, which can be numerically computed by the first-principles DFT method outlined above. If the time-reversal symmetry is preserved in the system, the Raman tensor is symmetric, i.e., Rij = Rji. However, for electronic and vibrational Raman scattering of a magnetic material, the Raman tensor can be much more complicated since the magnetic order can break the time-reversal symmetry of the system9,23,24,25,26,28,29. It can be composed of a symmetric (S) spin-independent part (RS) where ({R}_{ij}^{{rm{S}}}={R}_{ji}^{{rm{S}}}) and an anti-symmetric (AS) spin-dependent part (RAS) where ({R}_{ij}^{{rm{AS}}}=-{R}_{ji}^{{rm{AS}}}). In other words, R = RS + RAS, where
With R directly calculated from DFT, the elements in RS and RAS can be subsequently obtained since ({R}_{ij}^{{rm{S}}}=({R}_{ij}+{R}_{ji})/2) and ({R}_{ij}^{{rm{AS}}}=({R}_{ij}-{R}_{ji})/2). It is important to note that the SOC effect with the noncollinear magnetic order has to be turned on during the calculations of the polarizability tensor corresponding to each atomic displacement, in order to fully capture the effects of magnetism on the polarizability tensor, the Raman tensor, and subsequently the Raman intensity28,29. In other words, the anti-symmetric spin-dependent terms in RAS require the inclusion of the SOC effect.
Magneto-Raman effect in CrI3
Monolayer CrI3 has a space group of P3 (No. 143), and its unit cell has eight atoms and, thus, 24 normal phonon modes at the Γ point. They are decomposed as ΓP3 = 8A + 16E, where the non-degenerate A symmetry modes and two-fold degenerate E symmetry modes are Raman active. As the A modes exhibit strong Raman signals in experiments9,23,24,25, we will focus on the analysis of these modes. In particular, the A mode around 127.4 cm−1 appears as the strongest Raman peak in experiments, which corresponds to the phonon mode at 124.7 cm−1 in our calculations, as shown in Fig. 1. If we do not consider the magnetic order, the point group of P3 is C3 and the Raman tensor of an A mode assumes the form on the left side of Eq. (3), where the off-diagonal elements are absent.
As mentioned above, the inclusion of the magnetic order (either collinear or noncollinear) but without considering the SOC effect does not yield the spin-dependent off-diagonal terms in the Raman tensor. Only after considering the SOC effect28,29, the Raman tensor develops into a spin-dependent anti-symmetric form, as shown on the right side of Eq. (3). The corresponding calculated values of the Raman tensors are listed in Supplementary Table 1. According to ref. 9, the spin-independent part RS is invariant under time-reversal symmetry and independent of the magnetic structure, while the spin-dependent part RAS is odd under time-reversal symmetry. This indicates that the reversal of the magnetization direction should change the sign of the anti-symmetric off-diagonal terms in the Raman tensor (i.e., c changed to − c while − c changed to c in Eq. (3)), which is confirmed by our calculations in Supplementary Table 1. In a first-order approximation, one can say that the anti-symmetric off-diagonal terms are proportional to the magnetization of the system or the applied magnetic field24,26. This can also explain why the sign is changed upon the reversal of the magnetization direction of the magnetic field24,26. The anti-symmetric Raman tensors in FM monolayer CrI3 are important for the analysis of Raman tensors and Raman intensities in the few-layer CrI3 below.

a Calculated Raman spectra of monolayer CrI3 under both parallel (XX) and cross (XY) polarizations. The vibration patterns of the two A modes are shown as insets. The red arrows indicate the out-of-plane FM order. b Calculated resonant Raman intensity versus the photon energy for A1 (green) and A2 (orange) modes. The intensity ratio between A2 and A1 is shown on the right (blue). Note that the photon energy has been scaled by a factor of 1.63.
Figure 1a shows the calculated Raman spectra of monolayer CrI3 in both parallel (XX) and cross (XY) polarization channels. According to prior experimental results and hybrid DFT (HSE06) calculations, monolayer CrI3’s band gap is around 1.52 eV, while bulk CrI3 has a band gap of about 1.24 eV38,39,40. Our conventional DFT calculations, with the SOC effect included, found the band gap as 0.79 eV for 1L and 0.71 eV for bulk CrI3. Since DFT underestimates the band gap, the theoretical laser excitation energy is lower than the experimental one, and thus a scaling factor is required to make calculated Raman results more comparable with experiments. By matching calculated Raman spectra of 1L CrI3 with experimental data produced by a laser line of 1.96 eV9, we determined the scaling factor as 1.63 and thus the DFT laser excitation energy as 1.20 eV. Although hybrid functional methods such as HSE06 yield better descriptions of the electronic properties of CrI3, it is computationally prohibitive to use HSE06 for magneto-Raman simulations with noncollinear magnetism and SOC considered. Conventional DFT calculations with a scaling factor provide a balanced approach for satisfactory simulations of magneto-Raman spectra, as will be discussed below. In Fig. 1a, besides the strong peak (A2) at 124.7 cm−1, there is another mode A1 at 73.9 cm−1 24 with much weaker intensity (see Supplementary Table 2 for the Raman tensor). Most experimental works focused on A2 mode due to its stronger signal, but we found out that their intensities can be comparable at certain laser energies. Figure 1b shows the Raman intensities of both modes as a function of the photon energy, called Raman excitation profiles (REPs). At the lower energy window, A2 shows a stronger Raman resonance effect than A1. The Raman intensity of A2 reaches a peak near 2 eV, which is an order of magnitude stronger than that of A1. With the energy further increasing, the intensity of A2 begins to decrease, but it is still favorable for experiments to observe (such as the commonly used 1.96 and 2.33 eV laser lines). Starting around 2.5 eV, however, the signal of A2 decreases, while the intensity of A1 begins to substantially increase and grows even to be larger than that of A2 starting from 2.7 eV. At excitation energy at 3 eV, for example, A1 shows stronger signals than A2. This is confirmed by the calculated Raman spectra shown in Supplementary Fig. 1. The calculated REPs in Fig. 1b reveal two key findings: first, A1 can be effectively observed by using a laser line with energy higher than 2.7 eV, but the two most commonly used laser energies are 1.96 and 2.33 eV that often miss this mode; second, the exciton-phonon couplings for A1 and A2 have both similarities and differences that arise from their vibration patterns as illustrated in the insets in Fig. 1a. A2 has atomic vibrations predominantly out-of-plane, while A1 has both out-of-plane and in-plane vibration components. When an exciton has the wavefunction pointed out-of-plane, it should couple more strongly with A2, such as at energies between 1.9 and 2.5 eV in Fig. 1b; when an exciton has the wavefunction spreading along both out-of-plane and in-plane directions, it could couple with both A2 and A1 so that they share a similar resonance window, e.g., between 3 and 3.5 eV. Such vibration-dependent exciton-phonon coupling behaviors have been previously reported for MoS241 and black phosphorus42,43.
Because both A1 and A2 have out-of-plane components in the vibration patterns, they exhibit strong spin-dependent Raman scattering as their out-of-plane vibration components are coupled with the out-of-plane magnetic order in CrI3. As shown in Eq. (3), above the Curie temperature where the magnetic order disappears, off-diagonal terms R12 = R21 = 0, and thus A symmetry Raman peaks disappear in the cross (XY) polarization since I(XY) ∝ ∣R12∣29. However, below the Curie temperature where the monolayer has an FM order, off-diagonal terms are no longer zero (R12 = −R21 = c), and hence A1 and A2 are still present in the XY polarization, although the intensity is weaker than that from the parallel (XX) polarization since ∣R12∣2 < ∣R11∣2 (Fig. 1a). The appearance of a A symmetry Raman peak in the cross-polarization is an important sign that the Raman mode is spin-dependent. Note that A1 mode is very weak in the XY polarization when the excitation energy is 1.96 or 2.33 eV, since it is already weak in the XX polarization. This explains why A1 in the XY polarization cannot be differentiated from the background noise and ignored for the magneto-Raman study9. However, as discussed above, with a higher excitation energy above 2.7 eV, A1 is significantly enhanced in both XX and XY polarizations and it shows strong spin dependence just like A2 (see Supplementary Fig. 1).
As for doubly degenerated E modes, the Raman intensity is basically the same for cross and parallel polarization setups (Fig. 1a), a common behavior for 2D materials in hexagonal lattices21. This may lead to a partially incorrect conclusion that E modes are not correlated with the magnetic order. In fact, we discover that their Raman tensors also have spin-dependent anti-symmetric off-diagonal terms, which correspond to R13 and R23 instead of R12. In particular, for the characteristic E modes at 105.2 cm−1, one of the E modes (E1) has large ({R}_{13}^{{rm{AS}}}) while the other (E2) has large ({R}_{23}^{{rm{AS}}}) (see Supplementary Tables 3 and 4, respectively). In the typical back-scattering geometry, R13 and R23 are not probed by Raman scattering, and thus, their spin-dependent behaviors are completely missed. To study the magneto-Raman scattering of E modes, therefore, it is important to change the travel direction of laser light to the x or y direction so the polarizations of light have z components for effectively capturing the spin-dependent R13 and R23 of E modes. To the best of our knowledge, our calculations uncover the spin dependence of E modes for the first time in CrI3 and predict how they can be studied experimentally. These results demonstrate that the correlation between phonon modes and magnetic orders is much more broadly present than the sole A2 mode studied in experiments.
Bilayer CrI3 has Davydov splitting of each phonon mode from the monolayer due to the vdW interactions between the layers. For instance, mode A in the monolayer is split into two modes: even-parity mode Ag and odd-parity mode Bu. A simple linear chain model can explain the Davydov splitting23. As shown in the insets in Fig. 2, Ag mode vibrates in phase between layers while Bu mode vibrates out of phase between layers. Without a magnetic order, the crystallographic space group of the bilayer belongs to C2/m (No. 12), and the point group is C2h(2/m). Ag is Raman active with a symmetric form of Raman tensor. However, with the inclusion of the magnetic order and SOC, it becomes much more complicated with anti-symmetric off-diagonal terms. In addition, the Raman inactive mode Bu also becomes active under the interlayer AFM order and shows up under the cross-polarization setup, as shown in Fig. 2. Note that the irreducible representations of phonon modes, such as Ag and Bu, originate from the crystallographic space group and strictly speaking should change with the magnetic order to reflect the change in the Raman tensors. However, for simplicity and consistency, we keep using the same symmetry notations for phonon modes regardless of the magnetism. Our first-principles simulated Raman spectra are well consistent with prior experimental results9,23,24,25. More importantly, our numerical simulations provide an atomic scale picture of the spin-phonon coupling and reveal how a magnetic order can introduce the anti-symmetric Raman tensor elements, which enables us to propose a unified model to explain the spin dependence of Raman modes in 2D magnets.

Calculated Raman spectra of bilayer CrI3 under both parallel (XX) and cross (XY) polarizations when the interlayer magnetic order is a AFM or b FM. The vibration patterns of Ag and its Davydov-split counterpart Bu are shown as insets. The red and green arrows indicate the out-of-plane magnetic order.
Since the interlayer coupling in CrI3 is weak, we can treat each layer approximately as a single object23,44 and rewrite the Raman tensor in Eq. (1) with each layer instead of each atom as a unit:
where L is the layer index, N is the total number of layers, ({P}_{alpha beta }^{L}) represents the contribution of layer L to the change of polarizability by a unit displacement of the layer, and ΔdL indicates the displacement of layer L in a phonon mode. ({P}_{alpha beta }^{L}) can also be divided into a symmetric part PS and an anti-symmetric part PAS since χαβ has both symmetric and anti-symmetric parts. Eq. (4) represents a generalized and simplified polarizability model for us to gain crucial insights into Raman scattering of 2D magnets, as will be discussed below.
For the bilayer, R ∝ P1Δd1 + P2Δd2. The symmetric parts of the polarizability tensor of layer 1 and layer 2 are equal regardless of the magnetic order between the layers, i.e., ({chi }_{{rm{S}}}^{1}={chi }_{{rm{S}}}^{2}={chi }_{{rm{S}}}); however, the anti-symmetric parts depend on the relative spin orientation between the layers. In the AFM state, two layers have the opposite spin alignments, and thus ({chi }_{{rm{AS}}}^{1}=-{chi }_{{rm{AS}}}^{2}={chi }_{{rm{AS}}}) (recall a similar behavior that the reversal of the spin direction in the monolayer changes the sign of the anti-symmetric terms in the Raman tensor); in the FM state, it is natural that ({chi }_{{rm{AS}}}^{1}={chi }_{{rm{AS}}}^{2}={chi }_{{rm{AS}}}). Therefore ({P}_{{rm{S}}}^{1}={P}_{{rm{S}}}^{2}={P}_{{rm{S}}}), ({P}_{{rm{AS}}}^{1}=-{P}_{{rm{AS}}}^{2}={P}_{{rm{AS}}}) for the AFM state and ({P}_{{rm{AS}}}^{1}={P}_{{rm{AS}}}^{2}={P}_{{rm{AS}}}) for the FM state. For the monolayer, R ∝ P1Δd1 = (PS + PAS)Δd1. Therefore, PS and PAS essentially correspond to the symmetric and anti-symmetric part of Raman tensor of A2 mode in monolayer CrI3, respectively, and they should assume the following simple forms according to Eq. (3):
For the AFM state, we have R ∝ PS(Δd1 + Δd2) + PAS(Δd1 − Δd2). For Ag mode (Fig. 2), Δd1 ≈ Δd2, and thus ({{bf{R}}}_{{A}_{g}}^{{rm{AFM}}}propto 2{P}_{{rm{S}}}Delta {d}_{1}) that only has the symmetric part; for Bu mode, Δd1 ≈ −Δd2, and thus ({{bf{R}}}_{{B}_{u}}^{{rm{AFM}}}propto 2{P}_{{rm{AS}}}Delta {d}_{1}) that only has the anti-symmetric part. This analysis is confirmed by the numerical Raman tensors from DFT calculations shown in Supplementary Table 5. Clearly, the Raman tensor of Ag (Bu) mode almost only has diagonal (off-diagonal) terms so that it almost only appears in the XX (XY) polarization set-up, as shown in Fig. 2a. For the FM state, on the other hand, we have R ∝ (PS + PAS)(Δd1 + Δd2). For Bu mode, Δd1 ≈ − Δd2, so its Raman tensor is basically zero, which is consistent with our DFT calculations shown in Supplementary Table 5. For Ag mode, Δd1 ≈ Δd2, ({{bf{R}}}_{{A}_{g}}^{{rm{FM}}}propto 2({P}_{{rm{S}}}+{P}_{{rm{AS}}})Delta {d}_{1}), so the Raman tensor contains both the symmetric PS and anti-symmetric PAS, indicating a much more complicated form similar to the one on the right side of Eq. (3). As shown in Fig. 2b, we can observe that Bu mode is inactive but Ag mode is active in both polarization setups. These trends are similar to the A2 mode of the monolayer shown in Fig. 1, since the bilayer here shares the same FM state with the monolayer. To summarize, our model successfully explains the experimental and DFT simulated magneto-Raman behaviors of Ag mode and the Davydov-split component Bu. It also reveals two key ingredients governing the abnormal behaviors of Ag and Bu modes under the AFM order: first, the opposite spin orientation between two layers gives rise to opposite anti-symmetric terms in the polarizability and Raman tensors between two layers; second, Ag and Bu modes have the opposite layer vibration patterns. Therefore, the spin-phonon coupling in CrI3 depends on both the magnetic order (reflected in PAS) and the vibration pattern (reflected in ΔdL).
Trilayer CrI3 has three Davydov splitting components for each phonon mode from the monolayer. For A2 mode in monolayer, it is split into three modes: two even-parity modes Ag and one odd-parity mode Bu, as shown in Fig. 3a. Their simulated Raman spectra exhibit strong dependence on the interlayer magnetic order as demonstrated in Fig. 3b, c, in agreement with the prior experimental work23. We continue to deploy our model for understanding both the experimental and our calculated magneto-Raman spectra. The Raman tensor of the trilayer can be written as R ∝ P1Δd1 + P2Δd2 + P3Δd3. In the AFM state, the polarizability tensor of each layer follows ({chi }_{{rm{S}}}^{1}={chi }_{{rm{S}}}^{2}={chi }_{{rm{S}}}^{3}={chi }_{{rm{S}}}), ({chi }_{{rm{AS}}}^{1}=-{chi }_{{rm{AS}}}^{2}={chi }_{{rm{AS}}}^{3}={chi }_{{rm{AS}}}). Therefore ({P}_{{rm{S}}}^{1}={P}_{{rm{S}}}^{2}={P}_{{rm{S}}}^{3}={P}_{{rm{S}}}), ({P}_{{rm{AS}}}^{1}=-{P}_{{rm{AS}}}^{2}={P}_{{rm{AS}}}^{3}={P}_{{rm{AS}}}), leading to R ∝ PS(Δd1 + Δd2 + Δd3) + PAS(Δd1 − Δd2 + Δd3). For Bu mode, Δd2 = 0, and Δd1 ≈ − Δd3, so R is essentially zero and Bu mode does not become active in the AFM order. However, for ({A}_{g}^{1}) mode, Δd1 ≈ Δd3 = Δd, while (Delta {d}_{2}=-Delta {d}^{{prime} }) which is in the opposite direction but larger than Δd. This yields ({{bf{R}}}_{{A}_{g}^{1}}^{{rm{AFM}}}propto {P}_{{rm{S}}}(2Delta d-Delta {d}^{{prime} })+{P}_{{rm{AS}}}(2Delta d+Delta {d}^{{prime} })), where (2Delta d-Delta {d}^{{prime} }) is smaller than (2Delta d+Delta {d}^{{prime} }). Recall Eq. (5), and we can conclude that the diagonal terms in the Raman tensor of ({A}_{g}^{1}) are smaller than the off-diagonal anti-symmetric terms. Therefore, ({A}_{g}^{1}) exhibits stronger signals in the XY polarization setup than that in the XX polarization setup as demonstrated in Fig. 3b. On the other hand, for ({A}_{g}^{2}) mode, (Delta{d_{1}} approxDelta{d_{3}}=Delta{d}^{primeprime}), while (Delta{d_{2}}=Delta{d}^{primeprimeprime}) that is in the same direction as Δd1 and Δd3, leading to ({{bf{R}}}_{{A}_{g}^{2}}^{{rm{AFM}}}propto {P}_{{rm{S}}}(2Delta {d}^{{primeprime} }+Delta {d}^{primeprimeprime})+{P}_{{rm{AS}}}(2Delta {d}^{{primeprime} }-Delta {d}^{primeprimeprime})). It is apparent that the diagonal terms in the Raman tensor of ({A}_{g}^{2}) should be larger, and hence it shows stronger signals in Fig. 3b in the parallel polarization setup than that in the cross-polarization setup, which is the opposite behavior of ({A}_{g}^{1}).

a The vibration patterns of three Davydov splitting components for the Raman mode around 125 cm−1 in trilayer CrI3. Calculated Raman spectra of trilayer CrI3 under both parallel (XX) and cross (XY) polarizations when the interlayer magnetic order is b AFM or c FM.
In the FM state, the Raman tensor is rewritten as R ∝ PS(Δd1 + Δd2 + Δd3) + PAS(Δd1 + Δd2 + Δd3). Bu mode still has its Raman tensor as zero and is Raman inactive. For ({A}_{g}^{1}) mode, Δd1 ≈ Δd3 = Δd, (Delta {d}_{2}=-Delta {d}^{{prime} }), and thus ({{bf{R}}}_{{A}_{g}^{1}}^{{rm{FM}}}propto ({P}_{{rm{S}}}+{P}_{{rm{AS}}})(2Delta d-Delta {d}^{{prime} })). As (2Delta d-Delta {d}^{{prime} }) is very small, the Raman intensity of ({A}_{g}^{1}) is weak in both polarization channels, as shown in Fig. 3c. For ({A}_{g}^{2}) mode, ({{bf{R}}}_{{A}_{g}^{2}}^{{rm{FM}}}propto ({P}_{{rm{S}}}+{P}_{{rm{AS}}})(2Delta {d}^{{primeprime} }+Delta {d}^{primeprimeprime})). This suggests that ({A}_{g}^{2}) should be stronger in both polarization channels compared to ({A}_{g}^{1}) since (2Delta {d}^{{primeprime} }+Delta {d}^{primeprimeprime}) is much larger than (2Delta d-Delta {d}^{{prime} }), as confirmed by calculated Raman spectra in Fig. 3c. Moreover, the Raman tensor of ({A}_{g}^{2}) assumes a similar form to that of A2 mode in monolayer CrI3, since the latter is R ∝ (PS + PAS)Δd1 for the monolayer. Therefore, similar to A2 mode shown in Fig. 1a, the Raman intensity of ({A}_{g}^{2}) in the XX channel is stronger than that in the XY channel. Raman intensities of the relevant phonon modes in 1L-3L CrI3 based on the generalized polarizability model are summarized below (more details in Supplementary Information):
which are consistent with ab initio calculated Raman spectra shown in Figs. 1, 2, 3, respectively. In short, the model can explain the magneto-Raman scattering data in 1L-3L CrI3, as well as the four-layer and bulk CrI3 systems (see more details in Supplementary Information and Supplementary Figs. 2, 3). The comparison between Raman spectra from first-principles calculations and the model is shown in Supplementary Fig. 4. Our calculated spin-dependent polarized Raman spectra of 1L-4L CrI3 are also consistent with existing experimental data (Supplementary Fig. 5). Finally, we note that the circularly polarized magneto-Raman spectra can also be predicted by our DFT simulations and explained by the model (more details in Supplementary Information).
Magneto-Raman effect in MnBi2Te4
MnBi2Te4 is another 2D vdW magnetic material that has attracted broad research interests due to its novel topological and magnetic properties30,31,32. The crystal structure of bulk MnBi2Te4 belongs to the space group (Roverline{3}m) (No. 166) with the point group D3d, and it has three doubly degenerate Eg Raman modes and three non-degenerate A1g Raman modes45,46. Similar to CrI3, the spin-dependent anti-symmetric off-diagonal matrix elements of MnBi2Te4 appear in R13 and R23 for Eg modes, while in R12 for A1g modes. Therefore, in the typical back-scattering geometry, only A1g modes can be probed for the spin-Raman correlation (R12 for Eg modes is spin-independent). Following the similar steps of bilayer CrI3 discussed above, we have R ∝ PS(Δd1 + Δd2) + PAS(Δd1 − Δd2) for the AFM state. For a A1g mode in Fig. 4a (in-phase layer vibrations), Δd1 ≈ Δd2, and thus ({{bf{R}}}_{{A}_{1g}}^{{rm{AFM}}}propto 2{P}_{{rm{S}}}Delta {d}_{1}); for a A2u mode (out-of-phase layer vibrations), Δd1 ≈ −Δd2, and thus ({{bf{R}}}_{{A}_{2u}}^{{rm{AFM}}}propto 2{P}_{{rm{AS}}}Delta {d}_{1}). Therefore, the Raman tensor of a A1g (A2u) mode only has diagonal (off-diagonal) terms in the AFM state so that it only appears in the XX (XY) polarization channel, which well explains the spin-dependent behaviors of three A1g and A2u modes in our calculated Raman spectra shown in Fig. 4b. For the FM state, we have R ∝ (PS + PAS)(Δd1 + Δd2). For a A2u mode, Δd1 ≈ − Δd2, so ({{bf{R}}}_{{A}_{2u}}^{{rm{FM}}}approx 0). For a A1g mode, Δd1 ≈ Δd2, and thus ({{bf{R}}}_{{A}_{1g}}^{{rm{FM}}}propto 2({P}_{{rm{S}}}+{P}_{{rm{AS}}})Delta {d}_{1}). These results corroborate our calculations in Fig. 4c that any A2u mode is inactive, but any A1g mode can show up in both linear polarization setups (see circularly polarized Raman spectra of MnBi2Te4 in Supplementary information). In short, MnBi2Te4 exhibits a similar magneto-Raman effect to CrI3. Our first-principles calculations and the generalized model demonstrate that the magneto-Raman scattering of phonon modes should broadly exist in 2D magnets. Finally, we note that compared to CrI3, the spin-dependent off-diagonal terms in MnBi2Te4 are obviously smaller (Supplementary Table 7) since the intensities of A2u modes in Fig. 4b and A1g modes in Fig. 4c in the XY channel are in general considerably weaker than those in the XX channel (even the intensities in XX are reduced by five times). Such a result predicted by the first-principles simulations may explain why the spin-dependent phenomenon of MnBi2Te4 has not yet been experimentally observed since the spin-related Raman signals could be too weak. Based on our calculations, exploring different laser excitation energies other than 1.96 and 2.33 eV, or applying a strong magnetic field should boost the signals for detection.

a Calculated vibration patterns of three A1g Raman modes and their Davydov-split counterparts in bulk MnBi2Te4. Calculated Raman spectra of bulk MnBi2Te4 under both parallel (XX) and cross (XY) polarizations when the interlayer magnetic order is b AFM or c FM. The Raman spectra in the XX channel are reduced by five times, so A2u modes in (b) and A1g modes in (c) from the XY channel are more visible.
We have also studied bilayer MnBi2Te4 (Fig. 5), where a notable difference from the bulk is the appearance of three additional A1g and A2u Davydov pairs besides the three bulk-like A1g and A2u pairs. This is due to the symmetry reduction in the bilayer as the exterior atoms in each layer are no longer equivalent to the interior atoms. These new Davydov pairs are denoted as ({A}_{1g}^{{prime} }) and ({A}_{2u}^{{prime} }), ({A}_{1g}^{{primeprime} }) and ({A}_{2u}^{{primeprime} }), and ({A}_{1g}^{primeprimeprime}) and ({A}_{2u}^{primeprimeprime}), respectively, to differentiate the bulk-like pairs (({A}_{1g}^{1}) and ({A}_{2u}^{1}), ({A}_{1g}^{2}) and ({A}_{2u}^{2}), and ({A}_{1g}^{3}) and ({A}_{2u}^{3})). All six Davydov pairs in 2L MnBi2Te4 exhibit similar magneto-Raman behaviors to bulk MnBi2Te4. Just like the bulk, the magneto-Raman effect of most pairs in 2L MnBi2Te4 is weaker than CrI3 since the spin-related Raman intensities in the XY channel are generally much weaker than intensities in the XX channel (even the latter are reduced by five times for most Raman peaks). However, it is interesting to note that for the ({A}_{1g}^{prime} ) and ({A}_{2u}^{{prime} }) pair that only shows up in 2L MnBi2Te4, the intensities in the cross channel are stronger than or comparable with those in the parallel channel, indicating a strong magneto-Raman effect. This is probably due to strong displacements of the Mn atoms and the neighboring Te atoms in the vibrations of ({A}_{1g}^{{prime} }) and ({A}_{2u}^{{prime} }) modes as these atoms host the majority of magnetic moments in the system. For the ({A}_{1g}^{primeprimeprime}) and ({A}_{2u}^{primeprimeprime}) pair involving strong displacements of the Mn atoms, the magneto-Raman effect is also relatively strong compared to many other modes. In short, 2L and few-layer MnBi2Te4 could present interesting magneto-Raman phenomena for future experimental exploration.

a Calculated vibration patterns of six A1g and A2u Davydov pairs in 2L MnBi2Te4. Calculated Raman spectra of 2L MnBi2Te4 under both parallel (XX) and cross (XY) polarizations when the interlayer magnetic order is b AFM or c FM. The intensities of four Raman peaks in the XX channel are reduced by five times so Raman peaks from the XY channel are more visible.
Magneto-Raman effect in CrCl3 with in-plane magnetization
CrCl3 is a 2D layered magnet that features an in-plane easy axis where magnetic orders are within the x−y plane47,48,49, compared to CrI3 and MnBi2Te4 that have an out-of-plane easy axis. We investigated bilayer CrCl3, where the spins are set along with the x axis without loss of generality. Similar magneto-Raman effect is found for CrCl3 (Supplementary Fig. 8); however, there are two notable differences between CrCl3 and CrI3. First, due to the in-plane magnetic orders, the major spin-dependent anti-symmetric terms in Raman tensors of Ag and Bu modes appear in R23 and R32 for CrCl3 (Supplementary Table 6) compared to R12 and R21 for CrI3 (Supplementary Table 5). Therefore, the magneto-Raman phenomena of CrCl3 are better observed when the travel direction of light is along the x axis (i.e., the laser geometry should be (bar{,text{X}}text{(YY)X},) or (bar{,text{X}}text{(YZ)X},) instead of the common back-scattering (bar{,text{Z}}text{(XX)Z},) or (bar{,text{Z}}text{(XY)Z},), see Supplementary Fig. 8). Second, because of its weaker SOC29,50, the spin-dependent off-diagonal terms in Raman tensors of CrCl3 are significantly smaller than the diagonal terms and hence the intensities in the cross channel are much weaker than those in the parallel channel (Supplementary Fig. 8), in contrast to CrI3 where the intensities in both the channels are comparable (Fig. 2). The comparative study between CrCl3 and CrI3 highlights that the spin orientation and the SOC strength are both crucial in governing magneto-Raman behaviors.
Finally, to demonstrate the robustness of our results, we tested different functionals (LDA and PBE) and vdW correction models (semi-empirical D3 corrections and non-local vdW-DF functional optB86b-vdW) on 2L CrI3 and 2L MnBi2Te4. These different calculation methods give rise to minor differences in the spin-dependent polarized Raman spectra while the magneto-Raman behaviors and the optical selection rules remain the same (Supplementary Figs. 9, 10), proving the validity of our magneto-Raman findings.
Discussion
In summary, we have carried out a comprehensive study of the magneto-Raman effect in 2D magnets by using first-principles DFT calculations and a simple yet generalized polarizability model derived from the theory of Raman scattering. For 1L CrI3, the out-of-plane FM magnetic order is found to introduce considerable anti-symmetric off-diagonal terms in the Raman tensor of the characteristic A symmetry mode near 125 cm−1 (denoted as A2 or Ag), which alters its optical selection rule so that it can be observed in the cross (XY) polarization channel. Such strong spin-phonon coupling also gives rise to a notable degree of circular polarization (i.e., Raman intensity is different in LL and RR circular polarization channels). Without the FM order, our simulations prove that such a magneto-Raman effect is completely absent. Our results are well consistent with prior experimental works, and more importantly, we discover that such spin-dependent off-diagonal terms are also present for another A symmetry mode near 74 cm−1 and several doubly degenerated E modes in 1L CrI3. DFT calculations provide quantitative evidence on why they are almost not studied and predict how they can be experimentally observed. This demonstrates that the correlation between phonon modes and magnetic orders is not limited to the sole A symmetry Raman mode studied in experiments so far, and it should be a universal phenomenon.
For 2L-4L CrI3, the calculated linearly and circularly polarized Raman spectra for both interlayer FM and AFM orders also agree well with the experimental data, highlighting the strong magnetism dependence of the Davydov-split components of Ag mode. The normally Raman inactive Davydov modes with Bu symmetry remain silent under the FM state, but can become Raman active under the AFM state as the anti-symmetric off-diagonal terms appear in the Raman tensor; the Davydov modes with Ag symmetry have anti-symmetric off-diagonal terms under the FM state, but can have these terms disappear under the AFM state. Such spin-dependent behaviors change the selection rules for both linear and circular polarizations and allow magneto-Raman spectroscopy to monitor the magnetic order. Our generalized polarizability model successfully explains the magneto-Raman spectra of 2L-4L CrI3 and can be applied to any thickness in principles. It reveals that the spin-phonon coupling in 2D magnets depends on both the spin orientation (reflected in PAS) and the vibration pattern (reflected in ΔdL) of each layer. Furthermore, we extend our simulations and modeling to MnBi2Te4 and reveal a similar magneto-Raman effect. The quantitative results by DFT calculations unveil why it has not yet been observed experimentally and present a path forward on how we can observe it in future experiments. Finally, we investigated 2L CrCl3 that has an in-plane easy axis, in contrast to CrI3 and MnBi2Te4 that have an out-of-plane easy axis. Although the magneto-Raman effect is still present, the in-plane magnetic order in CrCl3 changes the location of the major spin-dependent off-diagonal terms in the Raman tensor matrix, leading to different polarization selection rules. The weak SOC in CrCl3 renders its magneto-Raman effect much weaker than CrI3, and thus selecting different laser excitation energies and applying a strong magnetic field may be required for experimental detection.
Our first-principles simulations and analytical modeling are complementary to each other in a sense that the former provides quantitative Raman spectra guiding experiments to locate phonon modes showing magneto-Raman behaviors and the latter offers qualitative analysis facilitating the understanding of the phenomena. They could play an important role in accelerating the magneto-Raman research in the ever-growing family of 2D magnets. With the prediction power of first-principles simulations, in the future, we will generate a computational magneto-Raman library (i.e., digital twin of magneto-Raman spectroscopy) to facilitate rapid characterization of magnetic ordering in layered 2D magnets.
Methods
DFT parameters
First-principles spin-polarized DFT calculations are carried out using the Vienna ab initio simulation package (VASP) with the projected augmented wave (PAW) method51,52. Local density approximation (LDA) is adopted for the electron exchange-correlation functional for CrI3, since the phonon frequencies computed by LDA are well consistent with experimental data. The generalized gradient approximation (GGA) of Perdew, Burke, and Ernzerhof (PBE)53 is also used for the exchange-correlation functional, where both the DFT-D3 method54 (denoted as PBE + D3) and the optB86b-vdW functional55 (denoted as PBE + optB86b) are added to describe the van der Waals (vdW) interactions between the layers. The PBE + D3 and PBE + optB86b, however, notably underestimate the phonon frequencies of CrI3, and thus LDA is used for phonon and Raman intensity calculations of CrI3. On the other hand, PBE+D3 describes the lattice dynamics of MnBi2Te4 reasonably well46, and hence it is used for computing properties of MnBi2Te4. To account for the strong electron correlations effect from the localized d orbitals of Cr and Mn atoms, we employ Dudarev’s DFT + U approach56, where the effective U term is Ueff = 3.0 eV for CrI35 while 4.0 eV for MnBi2Te430. For both CrI3 and MnBi2Te4, each layer prefers the FM order with an out-of-plane easy axis1,30. The weak interlayer interactions, however, make it possible to switch interlayer magnetism between the FM and AFM order9,23,24,25,26. Therefore, we consider both interlayer AFM and FM magnetic orders in our calculations. All the structures are relaxed until the maximum force on each atom is less than 1 meV/Å. The cutoff energy of the plane waves is chosen as 350 eV. The Brillouin zone (BZ) is sampled using Γ-centered 12 × 12 × 1 for mono- and few-layer CrI3 and bilayer MnBi2Te4, while 12 × 12 × 3 uniform k-grid for bulk CrI3 and MnBi2Te4. For the 2D systems, a vacuum region of more than 20 Å in the direction is used to avoid spurious interactions with the neighboring cells. We also studied bilayer CrCl3 using the same setting as bilayer CrI3 with the PBE+D3 method; however, the magnetic orders (both FM and AFM) were set in-plane for CrCl3.
Phonon calculations
Based on the optimized unit cell, we perform phonon calculations using a finite difference scheme implemented in the Phonopy57. We use VASP to compute Hellmann-Feynman forces in a supercell. For mono- and few-layer CrI3 and bilayer CrCl3, the size of the supercell is 2 × 2 × 1 with the corresponding k-grid reduced to 6 × 6 × 1; for bulk CrI3, the size of the supercell is 2 × 2 × 1 with the corresponding k-grid reduced to 6 × 6 × 3; for bilayer MnBi2Te4, the supercell size is 3 × 3 × 1 with the corresponding k-grid reduced to 4 × 4 × 1; for bulk MnBi2Te4, the supercell size is 3 × 3 × 1 with the corresponding k-grid reduced to 4 × 4 × 3. The choice of the supercell size is to ensure at least 10 Å in every lattice direction to avoid interactions between the neighboring cells. We introduce both positive and negative atomic displacements (Δ = 0.03 Å) to the supercell for static calculations to obtain the forces, which are processed by Phonopy to construct the dynamic matrix whose diagonalization provides phonon frequencies and eigenvectors.
Raman intensity calculations
Finally, Raman intensities are computed within the Placzek approximation20,21,22,33,34, as discussed above in the section of “Theory of resonant and magneto-Raman scattering”. Basically, one needs to calculate the derivatives of the polarizability tensor with respect to the atomic displacements for obtaining the Raman tensor. For both positive and negative atomic displacements along each direction in the unit cell (σ = 0.03 Å), the frequency-dependent polarizability tensor is computed by VASP (LOPTICS =. TRUE.), and then the derivatives of the polarizability tensor are obtained via the finite difference approach. The SOC effect with the noncollinear magnetic order has to be turned on during the VASP calculations of the polarizability tensor corresponding to each atomic displacement, in order to fully capture the effects of magnetism on the polarizability tensor, the Raman tensor, and subsequently, the Raman intensity28,29. The derivatives of the polarizability tensor, combined with the phonon frequencies and phonon eigenvectors, yield the Raman tensor matrix and, subsequently, the Raman intensity for every phonon mode under a given laser polarization set-up. Since DFT tends to underestimate the band gap and excitation energy, the DFT excitation energy is chosen as 1.20 eV for 1L CrI3 and 2L-4L CrI3 under the interlayer AFM state unless mentioned otherwise. For 2L-4L CrI3 under the interlayer FM state, the electronic band structures are changed, as evidenced by the fact that the band gaps are different from the ones under the AFM state by about 0.1–0.14 eV, and hence the DFT excitation energy is correspondingly changed to 1.30 eV. The DFT excitation energy is set as 1.62 eV for MnBi2Te4.
Responses