Physical origin of giant excitonic and magneto-optical responses in two-dimensional ferromagnetic insulators

The recent discovery of magnetism in atomically thin layers of van der Waals crystals has created great opportunities for exploring light–matter interactions and magneto-optical phenomena in the two-dimensional limit. Optical and magneto-optical experiments have provided insights into these topics, revealing strong magnetic circular dichroism and giant Kerr signals in atomically thin ferromagnetic insulators. However, the nature of the giant magneto-optical responses and their microscopic mechanism remain unclear. Here, by performing first-principles GW and Bethe-Salpeter equation calculations, we show that excitonic effects dominate the optical and magneto-optical responses in the prototypical two-dimensional ferromagnetic insulator, CrI3. We simulate the Kerr and Faraday effects in realistic experimental setups, and based on which we predict the sensitive frequency- and substrate-dependence of magneto-optical responses. These findings provide physical understanding of the phenomena as well as potential design principles for engineering magneto-optical and optoelectronic devices using two-dimensional magnets.

(solid red curve), Cr minor-spin 3d orbitals (dashed blue curve), I major-spin 5p orbitals (green dotted curve) and I minor-spin 5p orbitals (black solid curve). The energy of the VBM is set to zero.

C. Numerical evaluation of exciton radius
We use the exciton wave functions shown in Fig. 2 and Fig. 3 in the main text to evaluate the arithmetic mean radius ⟨| |⟩ and the root mean square radius 4⟨ 6

D. Spin-decomposed exciton probability amplitudes
Since the two-component spinor wave functions are used in our calculations, it is possible for us to choose certain spinor components of wave functions and calculate the spin configuration of constituent carriers for each exciton state. To be specific, an exciton wave function in the spinor formalism can be written as, The fractions for each spin configuration are listed in Supplementary

E. Computational details for Si and SiO2
The GW (at G0W0 level) and GW-BSE calculations of bulk Si are performed with the BerkeleyGW code 6 . The experimental lattice constant of a = 5.43 Å at 300 K is adopted in the calculations. The resulting frequency-dependent dielectric function of Si is shown in Supplementary Gaussian broadening is employed.
The reflective index for bulk Si is calculated as, where the static value Re[ ( → 0)] = 11.1 is in excellent agreement with experimental value of 11.68 or 11.4 7,8 .
For bulk SiO2, because the band gap (8.9 eV) is much larger than the energy range of interest in this problem (< 3.5 eV), we will only consider a static reflective index (with experimental value) with a very small imaginary part (~0.01i) to account for other dissipation channels 9 ,

F. Normal Modes
In the P-MOKE configuration with at least C3 rotational symmetry along the magnetization direction ( = l^), the dielectric tensor takes the following form, The Fresnel equation is given by, where n is the complex refractory index, = .
After solving the Fresnel equations with the dielectric function in Supplementary Eq. 5 with ∥ l^, we get the normal modes as the y and z circularly polarized plane waves, with distinct refractive indices, where the +(−) in ± denotes the circularly polarized light with the complex electric field amplitude along the direction of the spherical basis:

G. Kerr signals and Faraday signals
Following the convention in previous works 10 , we get the ratio of complex amplitude of ± circularly polarized reflection light through the ratio of the corresponding complex reflectivity ̃( ±) as, Suppose that an incident linearly polarized light is along the x-axis, the relative complex amplitudes of the reflection light • ‚ (at fixed and t) along x-and y-axis are then given by, which defines an ellipse oriented slightly away from the x-axis as shown in Supplementary  We can calculate the Kerr angle † and Kerr ellipticity † as 11 , and

H. Multi-interface P-MOKE Setup
Here we model the multi-interface P-MOKE setup in a systematic way 12,13 . The goal is to calculate the complex reflection coefficients for ± circularly polarized light, ̃( ±) = • ‚ (±) / • • (±) at the interface. The resulting complex reflection coefficients (at the topmost interface) for the three-interface model (shown in Fig. 4a in the main text) are given by, where the one-interface complex reflection coefficient between the i-th layer and the j-th layer is defined , and the light path within the i-th layer is defined as The twointerface model can be achieved by taking •ž ¢ → 0. The one-interface model can be achieved by taking •ž Ÿ → 0 and •ž ¢ → 0.

I. Effects of substrate thickness on MO signals
Here we investigate the influence of substrate thickness on MO signals. As described in the main text, we construct a three-interface P-MOKE setup with the order of vacuum-CrI3-SiO2-Si. We vary the thickness of the SiO2 layer ( (¨© ¢ ) and assume a semi-infinitely thick Si layer. We find that the interference of light reflected on each interface will be very sensitive to (

J. In-plane ferromagnetic monolayer CrI3
Our LSDA+U and G0W0 calculations have reproduced the reported magneto band structure effect where relevant bands are degenerate at the Γ point for the case of an in-plane ferromagnetic monolayer CrI3 14 .
We have obtained an indirect bandgap of 2.64 eV at the G0W0 level in the in-plane polarized structure, which is only slightly smaller than the direct bandgap at the Γ point (2.69 eV). The rotated magnetization has a strong impact on the polar MO signals, which can be understood from the symmetry. The dielectric function tensor is an axial tensor, which is in analogue to a dyad of two vectors. The broken C3 rotational

L. Crystal structure and structure relaxation
We use the experimental structure for both ferromagnetic bulk and monolayer CrI3. Bulk CrI3 belongs to the space group R 3 ® (148) 15 . Both bulk and monolayer CrI3 belong to the point group S6. In fact, we have checked the validity of the employed pseudopotentials with first-principles structure relaxation. During the relaxation, we used a kinetic energy cutoff of 120 Ry. The van der Waals interaction is included in two ways: the rVV10 nonlocal density functional 16 and the semiempirical Grimme's DFT-D3 method 17 , both of which have been implemented in the Quantum ESPRESSO package 18 . The structures have been fully relaxed until the force on each atom is less than 0.02 eV/Å. The spin-orbit coupling effects are not considered in the relaxation. We find that there is little deviation between the relaxed bulk structure and the experimental bulk structure. In addition, the lattice constants and internal coordinates barely change from bulk to monolayer structures, indicating much stronger intralayer bondings than interlayer bondings.
For these reasons, we use the experimental structure for both bulk and monolayer calculations. The detailed structure parameters are listed in Supplementary Table 3.

M. Dielectric function for quasi-2D materials
As an extensive physical quantity, the dielectric function is ill-defined for quasi-2D materials. In this work, we rescale the calculated dielectric function in a slab model by the thickness of a monolayer CrI3 (d=cbulk/3=6.6 Å), and where

N. Effect of a hexagonal boron nitride (hBN) substrate on exciton energies
To study the effects of an insulating substrate, we perform GW-BSE calculations of a ferromagnetic monolayer CrI3 on top of a monolayer hBN with an interlayer distance of 3.42 Å. The interlayer distance was determined with the van der Waals interaction through the DFT-D3 method 17 . It is expected that the substrate-induced exciton excitation energy redshift, albeit small, will quickly saturate with increasing thickness and additional layers of substrate will introduce negligible deviations [19][20][21] To simplify our following discussions, we now assume J = 0. The energy correction from +U can be estimated as, where ¿À is the occupation number for the diagonalized and spin-polarized local basis | ⟩, and the Bloch states | ¶ ⟩ are assumed to be fixed before and after +U. To get ¿À and | ⟩, we first use atomic wave functions (five spin-up d orbitals and five spin-down d orbitals for each Cr atom) as projectors to build the local orbital-resolved spin density ÅÅ AE À by summing over all the contributions from occupied Bloch states, with the quantization axis chosen along the z-axis. Here , ′ = 1, … ,5 labels the d orbitals and =↑, ↓ labels the spin polarization. We then diagonalize this spin density and get the eigenvalues ¿ À and eigenvectors | ⟩ as a linear combination of the five d orbitals with spin polarization . In this way, the eigenvalue correction should incorporate an extra factor to quantify how much the projection of the Bloch state is onto the localized basis. This is in contrast to the atomic limit, Δ • yº = ( ¡ 6 − • ), which indicates that the occupied and unoccupied states at the LSDA level are further split by U. But there is strong p-d hybridization in the CrI3 systems, especially for eg states which form bonds between Cr d orbitals and I p orbitals. In this way, the major-spin eg states become delocalized, which means the overlap between the first set of conduction Bloch states and the localized atomic orbital | ⟩ used in LSDA+U is small. To be explicit, Δ ¶ yº should be smaller for the eg states than for the t2g states, which has also been verified in our calculations. Moreover, the strong p-d hybridization makes the local occupation ¿À deviate from the atomic limit: that is, 3 major-spin d orbitals are completely occupied, while all the remaining d orbitals are empty. Here we used U = 2.0 eV in a case study. As listed in Supplementary   Table 4, the spin density has a more uniform distribution than that in the atomic limit. Also note that the local spin density is completely determined from the occupied Bloch states, while the overlap between unoccupied Bloch states and occupied diagonal local basis ( ÅÀ > 0.5) is not necessarily zero. Indeed, our numerical results agree with our intuition that the eg states will be more delocalized and have smaller overlap with localized diagonal orbitals. Also, the empty major-spin eg states have small but noticeable overlap between the occupied major-spin diagonal orbitals, giving rise to the negative shift of CBM energy with increasing U, as confirmed from both our Janak analysis and DFT calculations (last two columns in Supplementary Table 4). Moreover, the p-dominant VBM state (v1) at the Γ point has little overlap with those d-like localized diagonal orbitals, leading to negligible energy shift upon +U. For these reasons, the bandgap slightly decreases with increasing U.