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

Magnetism in two-dimensional (2D) van der Waals (vdW) layered materials has attracted much attention due to the fascinating physics and potential applications^1,2,3,4. CrI₃ is among the first 2D magnets to be discovered that exhibit spontaneous out-of-plane ferromagnetism in the monolayer limit¹. 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/rotation^5,6 and applying an electric/magnetic field^7,8,9. The manipulation of interlayer magnetism in vdW layered magnets results in rich physical properties, including novel layered-dependent magneto-optical effects^1,2, electric-field induced Kerr rotation^7,8, and giant second-order harmonic generation¹⁰.

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 CrI₃ and few-layer Cr₂Ge₂Te₆ in experiments by magneto-optical Kerr effect microscopy^1,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 terms^1,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 works^13,14,15. However, when the total magnetization of the system (e.g., bilayer CrI₃) is zero due to the interlayer AFM order, Kerr and Faraday rotations vanish since the off-diagonal terms in the polarizability tensor become zero¹. 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 magnets^16,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 vibration^20,21,22. Since the polarizability tensor has spin-dependent anti-symmetric off-diagonal terms^1,11,12, the Raman tensor has spin-dependent off-diagonal terms as well^{9,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 CrI₃ and VI₃^{9,23,24,25,26}. Furthermore, the magneto-Raman effect naturally captures the intricate spin-phonon coupling²⁷. 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 CrI₃^28,29, it has not yet been demonstrated for bilayer and few-layer CrI₃ (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 CrI₃ 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 CrI₃ with intralayer FM order and 2L-4L CrI₃ with both interlayer FM and AFM orders agree well with existing experimental data^9,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 MnBi₂Te₄, a novel 2D topological magnetic material^30,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 CrCl₃ 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 A² mode of CrI₃ studied in experiments.

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 e_s and e_i 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 approximation^{20,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, r_l(μ) 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 energy^20,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 E_L using the finite difference method^22,33,34. Here, we computed complex χ_αβ(E_L) 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 transformation³⁵. 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) ∝ ∣R₁₁∣² when e_s = e_i = (1, 0, 0); I(XY) ∝ ∣R₁₂∣² when e_s = (1, 0, 0), e_i = (0, 1, 0); I(YY) ∝ ∣R₂₂∣² when e_s = e_i = (0, 1, 0); I(YZ) ∝ ∣R₂₃∣² when e_s = (0, 1, 0), e_i = (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 energy^36,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., R_ij = R_ji. 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 system^{9,23,24,25,26,28,29}. It can be composed of a symmetric (S) spin-independent part (R^S) where ({R}_{ij}^{{rm{S}}}={R}_{ji}^{{rm{S}}}) and an anti-symmetric (AS) spin-dependent part (R^AS) where ({R}_{ij}^{{rm{AS}}}=-{R}_{ji}^{{rm{AS}}}). In other words, R = R^S + R^AS, where

With R directly calculated from DFT, the elements in R^S and R^AS 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 intensity^28,29. In other words, the anti-symmetric spin-dependent terms in R^AS require the inclusion of the SOC effect.

Magneto-Raman effect in CrI₃

Monolayer CrI₃ 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 experiments^9,23,24,25, we will focus on the analysis of these modes. In particular, the A mode around 127.4 cm⁻¹ appears as the strongest Raman peak in experiments, which corresponds to the phonon mode at 124.7 cm⁻¹ in our calculations, as shown in Fig. 1. If we do not consider the magnetic order, the point group of P3 is C₃ 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 effect^28,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. ⁹, the spin-independent part R_S is invariant under time-reversal symmetry and independent of the magnetic structure, while the spin-dependent part R_AS 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 field^24,26. This can also explain why the sign is changed upon the reversal of the magnetization direction of the magnetic field^24,26. The anti-symmetric Raman tensors in FM monolayer CrI₃ are important for the analysis of Raman tensors and Raman intensities in the few-layer CrI₃ below.

a Calculated Raman spectra of monolayer CrI₃ 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 A¹ (green) and A² (orange) modes. The intensity ratio between A² and A¹ is shown on the right (blue). Note that the photon energy has been scaled by a factor of 1.63.

Full size image

Figure 1a shows the calculated Raman spectra of monolayer CrI₃ in both parallel (XX) and cross (XY) polarization channels. According to prior experimental results and hybrid DFT (HSE06) calculations, monolayer CrI₃’s band gap is around 1.52 eV, while bulk CrI₃ has a band gap of about 1.24 eV^38,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 CrI₃. 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 CrI₃ with experimental data produced by a laser line of 1.96 eV⁹, 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 CrI₃, 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 (A²) at 124.7 cm⁻¹, there is another mode A¹ at 73.9 cm^−1 24 with much weaker intensity (see Supplementary Table 2 for the Raman tensor). Most experimental works focused on A² 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, A² shows a stronger Raman resonance effect than A¹. The Raman intensity of A² reaches a peak near 2 eV, which is an order of magnitude stronger than that of A¹. With the energy further increasing, the intensity of A² 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 A² decreases, while the intensity of A¹ begins to substantially increase and grows even to be larger than that of A² starting from 2.7 eV. At excitation energy at 3 eV, for example, A¹ shows stronger signals than A². This is confirmed by the calculated Raman spectra shown in Supplementary Fig. 1. The calculated REPs in Fig. 1b reveal two key findings: first, A¹ 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 A¹ and A² have both similarities and differences that arise from their vibration patterns as illustrated in the insets in Fig. 1a. A² has atomic vibrations predominantly out-of-plane, while A¹ 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 A², 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 A² and A¹ 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 MoS₂⁴¹ and black phosphorus^42,43.

Because both A¹ and A² 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 CrI₃. As shown in Eq. (3), above the Curie temperature where the magnetic order disappears, off-diagonal terms R₁₂ = R₂₁ = 0, and thus A symmetry Raman peaks disappear in the cross (XY) polarization since I(XY) ∝ ∣R₁₂∣²⁹. However, below the Curie temperature where the monolayer has an FM order, off-diagonal terms are no longer zero (R₁₂ = −R₂₁ = c), and hence A¹ and A² are still present in the XY polarization, although the intensity is weaker than that from the parallel (XX) polarization since ∣R₁₂∣² < ∣R₁₁∣² (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 A¹ 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 A¹ in the XY polarization cannot be differentiated from the background noise and ignored for the magneto-Raman study⁹. However, as discussed above, with a higher excitation energy above 2.7 eV, A¹ is significantly enhanced in both XX and XY polarizations and it shows strong spin dependence just like A² (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 lattices²¹. 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 R₁₃ and R₂₃ instead of R₁₂. In particular, for the characteristic E modes at 105.2 cm⁻¹, one of the E modes (E¹) has large ({R}_{13}^{{rm{AS}}}) while the other (E²) has large ({R}_{23}^{{rm{AS}}}) (see Supplementary Tables 3 and 4, respectively). In the typical back-scattering geometry, R₁₃ and R₂₃ 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 R₁₃ and R₂₃ of E modes. To the best of our knowledge, our calculations uncover the spin dependence of E modes for the first time in CrI₃ 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 A² mode studied in experiments.

Bilayer CrI₃ 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 A_g and odd-parity mode B_u. A simple linear chain model can explain the Davydov splitting²³. As shown in the insets in Fig. 2, A_g mode vibrates in phase between layers while B_u 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 C_2h(2/m). A_g 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 B_u 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 A_g and B_u, 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 results^9,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 CrI₃ under both parallel (XX) and cross (XY) polarizations when the interlayer magnetic order is a AFM or b FM. The vibration patterns of A_g and its Davydov-split counterpart B_u are shown as insets. The red and green arrows indicate the out-of-plane magnetic order.

Full size image

Since the interlayer coupling in CrI₃ is weak, we can treat each layer approximately as a single object^23,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 Δd_L indicates the displacement of layer L in a phonon mode. ({P}_{alpha beta }^{L}) can also be divided into a symmetric part P_S and an anti-symmetric part P_AS 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 ∝ P¹Δd₁ + P²Δd₂. 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 ∝ P¹Δd₁ = (P_S + P_AS)Δd₁. Therefore, P_S and P_AS essentially correspond to the symmetric and anti-symmetric part of Raman tensor of A² mode in monolayer CrI₃, respectively, and they should assume the following simple forms according to Eq. (3):

For the AFM state, we have R ∝ P_S(Δd₁ + Δd₂) + P_AS(Δd₁ − Δd₂). For A_g mode (Fig. 2), Δd₁ ≈ Δd₂, and thus ({{bf{R}}}_{{A}_{g}}^{{rm{AFM}}}propto 2{P}_{{rm{S}}}Delta {d}_{1}) that only has the symmetric part; for B_u mode, Δd₁ ≈ −Δd₂, 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 A_g (B_u) 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 ∝ (P_S + P_AS)(Δd₁ + Δd₂). For B_u mode, Δd₁ ≈ − Δd₂, so its Raman tensor is basically zero, which is consistent with our DFT calculations shown in Supplementary Table 5. For A_g mode, Δd₁ ≈ Δd₂, ({{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 P_S and anti-symmetric P_AS, 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 B_u mode is inactive but A_g mode is active in both polarization setups. These trends are similar to the A² 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 A_g mode and the Davydov-split component B_u. It also reveals two key ingredients governing the abnormal behaviors of A_g and B_u 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, A_g and B_u modes have the opposite layer vibration patterns. Therefore, the spin-phonon coupling in CrI₃ depends on both the magnetic order (reflected in P_AS) and the vibration pattern (reflected in Δd_L).

Trilayer CrI₃ has three Davydov splitting components for each phonon mode from the monolayer. For A² mode in monolayer, it is split into three modes: two even-parity modes A_g and one odd-parity mode B_u, 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 work²³. 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 ∝ P¹Δd₁ + P²Δd₂ + P³Δd₃. 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 ∝ P_S(Δd₁ + Δd₂ + Δd₃) + P_AS(Δd₁ − Δd₂ + Δd₃). For B_u mode, Δd₂ = 0, and Δd₁ ≈ − Δd₃, so R is essentially zero and B_u mode does not become active in the AFM order. However, for ({A}_{g}^{1}) mode, Δd₁ ≈ Δd₃ = Δ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 Δd₁ and Δd₃, 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⁻¹ in trilayer CrI₃. Calculated Raman spectra of trilayer CrI₃ under both parallel (XX) and cross (XY) polarizations when the interlayer magnetic order is b AFM or c FM.

Full size image

In the FM state, the Raman tensor is rewritten as R ∝ P_S(Δd₁ + Δd₂ + Δd₃) + P_AS(Δd₁ + Δd₂ + Δd₃). B_u mode still has its Raman tensor as zero and is Raman inactive. For ({A}_{g}^{1}) mode, Δd₁ ≈ Δd₃ = Δ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 A² mode in monolayer CrI₃, since the latter is R ∝ (P_S + P_AS)Δd₁ for the monolayer. Therefore, similar to A² 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 CrI₃ 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 CrI₃, as well as the four-layer and bulk CrI₃ 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 CrI₃ 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 MnBi₂Te₄

MnBi₂Te₄ is another 2D vdW magnetic material that has attracted broad research interests due to its novel topological and magnetic properties^30,31,32. The crystal structure of bulk MnBi₂Te₄ belongs to the space group (Roverline{3}m) (No. 166) with the point group D_3d, and it has three doubly degenerate E_g Raman modes and three non-degenerate A_1g Raman modes^45,46. Similar to CrI₃, the spin-dependent anti-symmetric off-diagonal matrix elements of MnBi₂Te₄ appear in R₁₃ and R₂₃ for E_g modes, while in R₁₂ for A_1g modes. Therefore, in the typical back-scattering geometry, only A_1g modes can be probed for the spin-Raman correlation (R₁₂ for E_g modes is spin-independent). Following the similar steps of bilayer CrI₃ discussed above, we have R ∝ P_S(Δd₁ + Δd₂) + P_AS(Δd₁ − Δd₂) for the AFM state. For a A_1g mode in Fig. 4a (in-phase layer vibrations), Δd₁ ≈ Δd₂, and thus ({{bf{R}}}_{{A}_{1g}}^{{rm{AFM}}}propto 2{P}_{{rm{S}}}Delta {d}_{1}); for a A_2u mode (out-of-phase layer vibrations), Δd₁ ≈ −Δd₂, and thus ({{bf{R}}}_{{A}_{2u}}^{{rm{AFM}}}propto 2{P}_{{rm{AS}}}Delta {d}_{1}). Therefore, the Raman tensor of a A_1g (A_2u) 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 A_1g and A_2u modes in our calculated Raman spectra shown in Fig. 4b. For the FM state, we have R ∝ (P_S + P_AS)(Δd₁ + Δd₂). For a A_2u mode, Δd₁ ≈ − Δd₂, so ({{bf{R}}}_{{A}_{2u}}^{{rm{FM}}}approx 0). For a A_1g mode, Δd₁ ≈ Δd₂, 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 A_2u mode is inactive, but any A_1g mode can show up in both linear polarization setups (see circularly polarized Raman spectra of MnBi₂Te₄ in Supplementary information). In short, MnBi₂Te₄ exhibits a similar magneto-Raman effect to CrI₃. 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 CrI₃, the spin-dependent off-diagonal terms in MnBi₂Te₄ are obviously smaller (Supplementary Table 7) since the intensities of A_2u modes in Fig. 4b and A_1g 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 MnBi₂Te₄ 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 A_1g Raman modes and their Davydov-split counterparts in bulk MnBi₂Te₄. Calculated Raman spectra of bulk MnBi₂Te₄ 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 A_2u modes in (b) and A_1g modes in (c) from the XY channel are more visible.

Full size image

We have also studied bilayer MnBi₂Te₄ (Fig. 5), where a notable difference from the bulk is the appearance of three additional A_1g and A_2u Davydov pairs besides the three bulk-like A_1g and A_2u 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 MnBi₂Te₄ exhibit similar magneto-Raman behaviors to bulk MnBi₂Te₄. Just like the bulk, the magneto-Raman effect of most pairs in 2L MnBi₂Te₄ is weaker than CrI₃ 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 MnBi₂Te₄, 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 MnBi₂Te₄ could present interesting magneto-Raman phenomena for future experimental exploration.

a Calculated vibration patterns of six A_1g and A_2u Davydov pairs in 2L MnBi₂Te₄. Calculated Raman spectra of 2L MnBi₂Te₄ 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.

Full size image

Magneto-Raman effect in CrCl₃ with in-plane magnetization

CrCl₃ is a 2D layered magnet that features an in-plane easy axis where magnetic orders are within the x−y plane^47,48,49, compared to CrI₃ and MnBi₂Te₄ that have an out-of-plane easy axis. We investigated bilayer CrCl₃, where the spins are set along with the x axis without loss of generality. Similar magneto-Raman effect is found for CrCl₃ (Supplementary Fig. 8); however, there are two notable differences between CrCl₃ and CrI₃. First, due to the in-plane magnetic orders, the major spin-dependent anti-symmetric terms in Raman tensors of A_g and B_u modes appear in R₂₃ and R₃₂ for CrCl₃ (Supplementary Table 6) compared to R₁₂ and R₂₁ for CrI₃ (Supplementary Table 5). Therefore, the magneto-Raman phenomena of CrCl₃ 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 SOC^29,50, the spin-dependent off-diagonal terms in Raman tensors of CrCl₃ 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 CrI₃ where the intensities in both the channels are comparable (Fig. 2). The comparative study between CrCl₃ and CrI₃ 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 CrI₃ and 2L MnBi₂Te₄. 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.

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 CrI₃, 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⁻¹ (denoted as A² or A_g), 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⁻¹ and several doubly degenerated E modes in 1L CrI₃. 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 CrI₃, 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 A_g mode. The normally Raman inactive Davydov modes with B_u 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 A_g 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 CrI₃ 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 P_AS) and the vibration pattern (reflected in Δd_L) of each layer. Furthermore, we extend our simulations and modeling to MnBi₂Te₄ 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 CrCl₃ that has an in-plane easy axis, in contrast to CrI₃ and MnBi₂Te₄ that have an out-of-plane easy axis. Although the magneto-Raman effect is still present, the in-plane magnetic order in CrCl₃ 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 CrCl₃ renders its magneto-Raman effect much weaker than CrI₃, 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.

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) method^51,52. Local density approximation (LDA) is adopted for the electron exchange-correlation functional for CrI₃, since the phonon frequencies computed by LDA are well consistent with experimental data. The generalized gradient approximation (GGA) of Perdew, Burke, and Ernzerhof (PBE)⁵³ is also used for the exchange-correlation functional, where both the DFT-D3 method⁵⁴ (denoted as PBE + D3) and the optB86b-vdW functional⁵⁵ (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 CrI₃, and thus LDA is used for phonon and Raman intensity calculations of CrI₃. On the other hand, PBE+D3 describes the lattice dynamics of MnBi₂Te₄ reasonably well⁴⁶, and hence it is used for computing properties of MnBi₂Te₄. To account for the strong electron correlations effect from the localized d orbitals of Cr and Mn atoms, we employ Dudarev’s DFT + U approach⁵⁶, where the effective U term is U_eff = 3.0 eV for CrI₃⁵ while 4.0 eV for MnBi₂Te₄³⁰. For both CrI₃ and MnBi₂Te₄, each layer prefers the FM order with an out-of-plane easy axis^1,30. The weak interlayer interactions, however, make it possible to switch interlayer magnetism between the FM and AFM order^{9,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 CrI₃ and bilayer MnBi₂Te₄, while 12 × 12 × 3 uniform k-grid for bulk CrI₃ and MnBi₂Te₄. 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 CrCl₃ using the same setting as bilayer CrI₃ with the PBE+D3 method; however, the magnetic orders (both FM and AFM) were set in-plane for CrCl₃.

Phonon calculations

Based on the optimized unit cell, we perform phonon calculations using a finite difference scheme implemented in the Phonopy⁵⁷. We use VASP to compute Hellmann-Feynman forces in a supercell. For mono- and few-layer CrI₃ and bilayer CrCl₃, the size of the supercell is 2 × 2 × 1 with the corresponding k-grid reduced to 6 × 6 × 1; for bulk CrI₃, the size of the supercell is 2 × 2 × 1 with the corresponding k-grid reduced to 6 × 6 × 3; for bilayer MnBi₂Te₄, the supercell size is 3 × 3 × 1 with the corresponding k-grid reduced to 4 × 4 × 1; for bulk MnBi₂Te₄, 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 approximation^{20,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 intensity^28,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 CrI₃ and 2L-4L CrI₃ under the interlayer AFM state unless mentioned otherwise. For 2L-4L CrI₃ 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 MnBi₂Te₄.