Magnetism and Magneto-optical Effects in Bulk and Few-layer CrI$_3$: A Theoretical GGA + U Study

The latest discovery of ferromagnetism in atomically thin films of semiconductors Cr$_2$Ge$_2$Te$_6$ and CrI$_3$ has unleashed numerous opportunities for fundamental physics of magnetism in two-dimensional (2D) limit and also for technological applications based on 2D magnetic materials. In this paper, we present a comprehensive theoretical study of the magnetic, electronic, optical and magneto-optical(MO) properties of multilayers [monolayer(ML), bilayer and trilayer] and bulk CrI$_3$, based on the density functional theory with the generalized gradient approximation plus on-site Coulomb repulsion scheme. Interestingly, all the structures are found to be single-spin ferromagnetic(FM) semiconductors. They all have a large out-of-plane magnetic anisotropy energy(MAE) of $\sim$0.5 meV/Cr. These large MAEs suppress transverse spin fluctuations and thus stabilize long-range magnetic orders at finite temperatures down to the ML limit. They also exhibit strong MO effects with their Kerr and Faraday rotation angles being comparable to that of best-known bulk MO materials. The shape and position of the main features in the optical and MO spectra are found to be nearly thickness-independent although the magnitude of Kerr rotation angles increases monotonically with the film thickness. Magnetic transition temperatures estimated based on calculated exchange coupling parameters, calculated optical conductivity, MO Kerr and Faraday rotation angles agree quite well with available experimental data. The calculated MAE as well as optical and MO properties are analyzed in terms of the calculated orbital-decomposed densities of states, band state symmetries and dipole selection rules. Our findings of large out-of-plane MAEs and strong MO effects in these single-spin FM semiconducting CrI$_3$ ultrathin films suggest that they will find valuable applications in semiconductor MO and spintronic nanodevices.


Introduction
The recent discovery of intrinsic ferromagnetism in atomically thin films of semiconductors Cr 2 Ge 2 Te 6 [1] and CrI 3 [2] has opened numerous exciting opportunities in two-dimensional (2D) magnetism. Amalgamation of magnetism and 2D materials are highly desirable for at least two reasons, one for exploration of fundamental physics of magnetism in 2D limit and the other for fascinating technological applications ranging from magnetic memories, to sensing, to spintronics to magneto-optical devices based on 2D materials. Therefore, these atomically thin magnetic materials are currently subject to intensive investigations. In fact, bulk CrI 3 has already been well studied experimentally in the past [3][4][5] mainly because it is a layered ferromagnetic semiconductor with large Kerr and Faraday rotations which promise such device applications as optical insulators, magnetic sensors and rewritable optical memories [6]. Its 2D magnetic structures are presently attracting much renewed interest. For example, Seyler et al. [7] recently reported their observation of magneto-photoluminescence in monolayer and bilayer CrI 3 . Excitingly, manipulation of 2D magnetism (e.g., switching of magnetization direction and tuning antiferromagnetic to ferromagnetic transition) in BL CrI 3 by either applying a vertical electric field [8,9] or electrostatic doping [10] has also been recently demonstrated. Further, giant tunneling magnetoresistance has been observed in magnetic tunnel junctions made of atomically thin CrI 3 and othe van der Waals (vdW) materials. [11,12] Unprecedented control of spin and valley pseudospin in ultrathin CrI 3 and WSe 2 hetrostructures has also been demonstrated [13].
In order to exploit these emergent phenomena for various applications, the mechanisms that control the physical properties of these 2D materials should be thoroughly understood. As a step towards this goal, we present in this paper a comprehensive firstprinciples density functional study of the magnetic, electronic, optical and magneto-optical properties of multilayers [monolayer (ML), bilayer (BL) and trilayer (TL)] as well as bulk CrI 3 . Specifically, we focus the magnetic interactions, magnetic anisotropy energy and magneto-optical effects in these magnetic structures. Although magnetic exchange coupling and magnetic anisotropy in ML CrI 3 have been theoretically investigated by several groups based on density functional calculations [14][15][16][17], they have not been studied in BL, TL and bulk CrI 3 . Furthermore, the interlayer exchange coupling in bulk and TL CrI 3 has not been addressed although that in the BL is presently under intensive investigations [18][19][20][21].
Magnetic anisotropy energy (MAE) refers to the energy required to flip the magnetization from the easy to the hard axis, and is one of the principal specification parameters for a magnetic material. In fact, MAE is particularly important for 2D magnetic materials. According to the Mermin-Wagner theorem [22], a long-range magnetic order could not occur at any finite temperature in a 2D isotropic Heisenberg magnetic structure because of large thermal spin fluctuations in the 2D system. However, this Mermin-Wangner restriction [22] can be lifted if the 2D system has a significant outof-plane MAE which would suppress the thermal fluctuations and thus stabilize the long-range magnetic order at finite temperature even in the ML limit, as was observed recently in ML CrI 3 [2]. The MAE of a magnetic solid arises from two contributions, namely, the magnetocrystalline anisotropy energy (MCE) (∆E b ) due to the effect of electron relativistic spin-orbit coupling (SOC) on the band structure, and the magnetic dipolar anisotropy energy (MDE) (∆E d ) due to the magnetostatic dipole-dipole interaction in the magnetic solid. In a layered material, the MDE always prefers an in-plane magnetization while the MCE could favor either an in-plane or the out-ofplane magnetization depending on the band structure of the material. Although the MCE of ML CrI 3 has been calculated [14,17,23], the MDE was not considered in these previous reports. Although the MDE is negligibly small for the cubic and isotropic materials such as bcc Fe and fcc Ni, it could become significant in low-dimensional materials [24][25][26]. In fact, it was found recently that although the MCE in ML Cr 2 Ge 2 Te 6 favors the out-of-plane magnetization, the MDE is larger than the MCE such that ML Cr 2 Ge 2 Te 6 would have an in-plane magnetization [26], thereby explaining why the longe-range ferromagnetic order was not observed in the ML [1]. Furthermore, there are no reports on the MAE of bulk, BL and TL CrI 3 . Therefore, in this paper we present both calculated MCE and MDE for all the considered CrI 3 structures.
Magneto-optical (MO) effects are prominent manifestations of light-matter interactions in magnetic materials. When a linearly polarized light beam hits a magnetic material, the reflected and transmitted light beams would become elliptically polarized with the principal axis rotated with respect to the polarization direction of the incident light beam, as illustrated in Fig. 1(a). The former and latter are called MO Kerr (MOKE) and MO Faraday (MOFE) effects [27], respectively. MOKE has been widely used to probe the magnetic and electronic properties of solids, surface and thin films [27]. Indeed, long-range ferromagnetic orders in atomically thin films of Cr 2 Ge 2 Te 6 and CrI 3 were discovered by using the MOKE technique. [1,2] In contrast, MOFE has been less explored mainly because light can only transmit through very thin

films.
Magnetic materials with large Kerr and Faraday rotations are the promising candidates for valuable magneto-optical device applications [28,29], and thus they have continuously attracted attention in the past several decades. The latest discovery of 2D ferromagnetic semiconductors [1,2] provides especially exciting possibilities of scaling the optical and magneto-optical devices down to subnanometer scale. Therefore, here we carry out a systematic firstprinciples density functional study of the optical and MO properties of bulk and multilayer CrI 3 . Indeed, our study reveals that bulk and multilayer CrI 3 exhibit large MO effects in a wide optical frequency range with Kerr rotation angles being as large as ∼2.3 • and Faraday rotation angles being in the order of ∼200 • /µm. This suggests that 2D ferromagnetic CrI 3 semiconductor structures will provide an interesting material platform for further studies of novel magnetooptical phenomena and technological applications.

Computational methods
In this paper, we study the electronic, magnetic, and magneto-optical properties of bulk and multilayer CrI 3 structures. Bulk CrI 3 forms a layered structure with MLs separated by the van der Waals gap [ Fig.   1(c)]. Each CrI 3 ML consists of edge-sharing CrI 6 octahedra forming a planar network with Cr atoms in a honeycomb lattice [Figs 1(b)] and thus there are two formula units (f.u.) per lateral unit cell. These MLs are then stacked in an ABC sequence, thus resulting in a rhombohedral crystal with R3 symmetry and two f.u. per unit cell. This structure can also be regarded as an ABC-stacked hexagonal crystal with experimental lattice constants a = b = 6.867Å and c = 19.807Å [30]. We consider the experimental rhombohedral primitive cell in the bulk calculations. For the multilayer CrI 3 structures, the hexagonal unit cell [Figs 1(b)] with the experimental lattice constants is considered as the lateral unit cell. For BL and TL CrI 3 structures [see Fig. 1(c)], we consider the AB and ABC stackings, respectively. The slab-superlattice approach is used to construct the multilayer structures with the separations of neighboring slabs being at least 20Å.
First-principles density functional calculations are performed using the accurate projector augmented wave (PAW) method [31] as implemented in the Vienna ab-initio simulation package (VASP) [32,33]. The fully relativistic PAW potentials are adopted in order to include the SOC. The valence configurations of Cr and I atoms adopted in the calculations are 3d 5 4s 1 , 5s 2 5p 5 , respectively. A large plane-wave energy cutoff of 400 eV is used throughout the calculations. For the Brillouin zone integrations, k-point meshes of 16×16×16 and 20×20×1 are used for bulk and multilayers CrI 3 , respectively.
The generalized gradient approximation (GGA) to the exchange-correlation potential of Perdew-Burke-Ernzerhof parametrization [34] is used. However, an improved account of onsite electron correlation among 3d electrons is needed for a better description of the electronic and magnetic structures of 3d transition metal compounds (see, e.g., Ref. [35] and references therein). To better describe onsite Coulomb repulsion between Cr 3d electrons, we adopt the GGA+U scheme [36] in the present calculations. Feng et al. [26] recently carried out systematic GGA+U calculations for bulk and multilayer Cr 2 Ge 2 Te 6 with different effective Coulomb energy U values and concluded that the appropriate U value for Cr 3d electrons is 1.0 eV. Given the similarity between the Cr 2 Ge 2 Te 6 and CrI 3 , here we also use U = 1.0 eV. Nevertheless, we also perform a series of the GGA+U calculations with the U value ranging from 0 to 3 eV for all the considered CrI 3 systems and find that the calculated physical quantities do not depend significantly on the U value used [see Supplementary Note A, table S1 and Fig. S1 in the Supplementary Information (SI)]. This indicates that the results presented below would remain nearly the same even if a different U value were used.
We consider four intra-layer magnetic configurations, namely, one FM as well as three antiferromagnetic (AF) structures, labelled as AF-Néel, AF-zigzag, and AF-stripe in Fig. 2. For bulk, BL and TL CrI 3 , we also consider interlayer FM and AF magnetic configurations with intralayer FM CrI 3 MLs. To understand the magnetic interactions and also to estimate magnetic ordering temperature (T c ) for bulk and multilayers CrI 3 , we calculate the exchange coupling parameters by mapping the calculated total energies of these magnetic configurations (see Fig. 2) onto the classical Heisenberg Hamiltonian, where E 0 donates the nonmagnetic ground state energy; J ij the exchange coupling parameter between sites i and j;ê i , the unit vector representing the direction of the magnetic moment on site i. For one CrI 3 , we obtain a set of four linear equations of J 1 , J 2 and J 3 : , which can be solved to estimate J 1 , J 2 and J 3 . Similarly, one can estimate the effective interlayer coupling parameter J z based the calculated total energies per unit cell of two magnetic configurations by using expressions for the TL and J z = (E F iM − E F M )/8 for the bulk where E F iM denotes the total energy of the ferrimagnetic configuration with the Cr magnetic moments on one of the three layers in the hexagonal unit cell fliped.
We also calculate the MCE, MDE and hence MAE for all the considered structures in the FM state. We first perform two relativistic electronic structure calculations for the in-plane and out-ofplane magnetizations, and then obtain the MCE (∆E b ) as the total energy difference between the two calculations. Highly dense k-point meshes of 24×24×24 and 28×28×1 are used for bulk and multilayers CrI 3 , respectively. Further calculations by using different k-point meshes show that thus-obtained MCEs converge within 1 %. For a FM system, the MDE ∆E d in atomic Rydberg units is given by [24,25] where the speed of light c = 274.072 and the so-called magnetic dipolar Madelung constant where R are the lattice vectors, q are the atomic position vectors in the unit cell, and m q is the atomic magnetic moment (in units of µ B ) on site q. Using the calculated magnetic moments, the ∆E d is obtained as the difference in E d between the in-plane and out-ofplane magnetizations. In a 2D material, since all the R and q are in-plane, the second term in Eq. (2) would be zero for the out-of-plane magnetization, thereby resulting in the positive M qq ′ while for an in-plane magnetization the M qq ′ is negative. [26] Therefore, the MDE always prefers an in-plane magnetization in a 2D material. This is purely a geometric effect and consequently the MDE is also known as the magnetic shape anisotropy energy. In a FM solid with trigonal symmetry and magnetization rotational along z-axis, the optical conductivity tensor consists of three independent elements, namely, σ xx , σ zz and σ xy . Within the Kubo linear response theory [37,38], the adsorptive parts of these elements are given by where ω is the photon energy, and ǫ ki is the ith band energy at point k. Summations i and j are over the occupied and unoccupied bands, respectively. Dipole matrix elements p a ij = kj|p a |ki wherep a denotes Cartesian component a of the dipole operator, are obtained from the band structures within the PAW formalism [39], as implemented in the VASP package. The integration over the Brillouin zone is carried out by using the linear tetrahedron method (see Ref. [40] and references therein). The dispersive parts of the optical conductivity elements can be obtained using the Kramers-Kroing relations where P donates the principle value. The quasiparticle lifetime effects are accounted for by convoluting all the optical conductivity spectra with a Lorentzian of line width Γ. For layered van der Waals materials such as graphite, Γ is about 0.2 eV (see, e.g., Figs. 1(a) and 1(b) in Ref. [41]), which is thus used in this paper.
For a bulk magnetic material, the complex polar Kerr rotation angle is given by [42,43], However, for a magnetic thin film on a nonmagnetic substrate, the complex polar Kerr rotation angle is given by [44,45] where d stands for the thickness of the magnetic layer and ε s xx (σ s xx ) the diagonal part of the dielectric constant (optical conductivity) of the substrate. Experimentally, atomically thin CrI 3 films were prepared on SiO 2 /Si [2,7,13]. Thus, the dielectric constant (ε s xx = 2.25) of bulk SiO 2 is used here. Similarly, the complex Faraday rotation angle for a thin film can be written as [46] where n + and n − represent the refractive indices for left-and right-handed polarized lights, respectively, and are related to the optical conductivity (or dielectric function) via expressions n 2 Here the real parts of the optical conductivity σ ± can be written as , and thus σ xy would be nonzero only if σ + and σ − are different. In other words, magnetic circular dichroism is essential for the nonzero σ xy and hence the magneto-optical effects.  To find the ground state magnetic structure and also evaluate exchange coupling among the magnetic Cr atoms in all the considered CrI 3 structures, we first consider four possible intra-layer magnetic configurations, namely, the FM, Néel-type antiferromagnetic (AF-Néel), zigzag-antiferromagnetic (AF-zigzag) and stripe-like antiferromagnetic (AF-stripe) structures (see Fig. 2). The calculated total energies for these magnetic configurations indicate that for all the considered structures, the FM state is the most stable configuration. We then consider the interlayer ferromagnetic and antiferromagnetic [see Fig. 2(e)] structures for bulk, BL and TL CrI 3 , and again we find the interlayer FM state is more stable. Therefore, we list the calculated spin and orbital magnetic moments of bulk and multilayer CrI 3 in the FM state only in Table  1. The calculated Cr spin magnetic moment in all the CrI 3 structures is ∼3.2 µ B , being in good agreement with the experimental value of ∼3.0 µ B [30]. This is also consistent with three unpaired electrons in the Cr t 2g configuration in these structures. Further, the calculated Cr orbital magnetic moment for all the investigated structures is small (∼ 0.08 µ B ), thus lending support to the notion of the completely filled spin-up Cr t 2g subshell in these systems. Interestingly, there is also a small induced magnetic moment (-0.09 µ B ) on the I atom (Table 1), which is antipallel to that of the Cr magnetic moment, in all the systems. All these result in a total magnetic moment of 3.00 µ B /f.u. for all the structures. As mentioned in Sec. 2, using the total energies of the considered magnetic configurations, we estimate first three near-neighbor exchange coupling parameters (J 1 , J 2 , J 3 ) in each ML in all the considred systems and also the effective interlayer exchange coupling constant (J z ). The calculated exchange coupling constants are listed in Table 2. It is clear from Table  2 that the first near-neighbor exchange coupling (J 1 ) is ferromagnetic (i.e., positive) in all the considered CrI 3 structures. The second near-neighbor coupling (J 2 ) is also ferromagnetic (positive), although the J 2 value is about 3.5 times smaller than J 1 for all the structures. Interestingly, Table 2 shows that the magnitudes of both J 1 and J 2 increase monotonically with the number of layers, being consistent with the observed increase in magnetic transition temperature in these systemsi [2]. The third near-neighbor magnetic coupling (J 3 ) is also ferromagnetic (positive) in all the CrI 3 structures except ML CrI 3 where J 3 is negative (antiferromagnetic). Nonetheless, the calculated J 3 values are found to be 5 times smaller than J 1 and  only about half of J 2 for all the structures. For the CrI 3 ML, we find that all the calculated exchange coupling parameters are in good quantitive agreement with the recent GGA calculations [14][15][16], although the J 1 is slightly larger than that of the previous GGA+U calculation [17] with U = 2.7 eV and J = 0.7 eV. Note that only the first near-neighbor magnetic interaction was considered in Ref. [17] and also that the J values reported previously [14][15][16][17] should be multiplied by a factor of S 2 = 9/4 in order to compare with the present J values. Finally, we notice that the calculated interlayer exchange coupling constant J z is also ferromagnetic for bulk, BL and TL structures. This prediction agrees well with all previous experiments [2][3][4][5] on all the structures except BL CrI 3 which has been experimentally found to be interlayer antiferromagnetically coupled [2]. A microscopic explanation of this notable discrepancy between theory and experiment on BL CrI 3 is nontrivial. Indeed, very recently this problem has already been investigated theoretically by several groups [18][19][20][21]. These theoretical studies all showed that the interlayer magnetic coupling in BL CrI 3 is strongly stacking dependent and is also sensitive to the exchange-correlation functional used. In particular, the GGA + U calculation with U = 3 eV [18] indicates that in the AB ′ -stacked BL CrI 3 , the total energy of the AF configuration is slightly lower than that of the FM configuration (by merely ∼0.7 meV/unit cell). However, the AB-stacking order which is the global minimum and favors the FM state, still have a total energy significantly lower than the AB ′ -stacked structure (by ∼6 meV/unit cell). Consequently, it was proposed [18] that BL CrI 3 exfoliated at room-temperature is perhaps kinetically trapped in the metastable AB ′ -stacked structure, on cooling and thus the AF state was observed instead [2]. Furthermore, it was also found that even in the AB ′ -stacking order, the AF interlayer coupling becomes favorable only when U is larger than 2.5 eV, [19] and this result is confirmed by our own GGA + U calculations. Clearly, further measurements on the stacking order in BL CrI 3 are needed to resolve this interesting problem. To this end, we calculate the optical and magneto-optical properties of BL CrI 3 in both AB-and AB ′ -stacking orders. Unfortunately, we find that the magneto-optical effects of BL CrI 3 hardly depend on the stacking order, as presented in Figs. S4 and S5 in the supplementary information.
3.1.2. Magnetic anisotropy energy and ferromagnetic transition temperature: The calculated magnetic anisotropy energies (△E ma ) for bulk and multilayers CrI 3 are presented in Table 1. As mentioned before, the MAE consists of two competing contributions, namely, magnetocrystalline anisotropy energy (△E b ) due to the SOC effect on the band structure and magnetic dipolar (shape) anisotropy energy (△E d ). First of all, Table 1 shows that as could be expected from the 2D nature of their structures, the MDE per ML in all the considered systems is rather large and independent of the film thickness. Furthermore, the MDE is always negative and thus always prefers an in-plane magnetization. Remarkably, the MCE per ML is a few times larger than the MDE, and positive, meaning that it prefers the out-of-plane magnetization. The sum of these competing terms (△E b + △E d ) thus results in the large and positive MAE (Table 1), which is more or less independent of the film thickness. This is in strong contrast to the case of the Cr 2 Ge 2 Te 6 multilayers where the MAE is much smaller and depends significantly on the film thickness [26]. In particular, for ML Cr 2 Ge 2 Te 6 , the magnitude of MCE becomes smaller than that of MDE such that the MAE is negative, favoring an in-plane magnetization [26]. Importantly, since the significant out-of-plane MAE is needed to stabilize long-range magnetic orders in a 2D material by suppressing the thermal spin fluctuations, this explains why the FM state is observed in ML CrI 3 [2] but not in ML Cr 2 Ge 2 Te 6 [1]. These contrasting magnetic anisotropic properties between multilayers CrI 3 and Table 2. Intralayer first-, second-and third-neighbor exchange coupling parameters (J 1 , J 2 , J 3 ) as well as effective interlayer exchange coupling constant (Jz) in ML, BL, TL and bulk CrI 3 , derived from the calculated total energies for various magnetic configurations (see the text). Ferromagnetic transition temperatures estimated using the derived J values within the mean-field approximation (T m c ) and also within the spin wave theory (SWT) (T SW T c ) (see the text) are also included. For comparison, the available experimental Tc values are also listed in brackets. Cr 2 Ge 2 Te 6 thus lend support to the notion that the former is well described by the 2D Ising model while the latter is well described by the Heisenberg model. [47]. Note that the present MCE of ML CrI 3 is in good quantitative agreement with the theoretical value of 0.686 meV/f.u. reported recently in Ref. [14] in which the MDE was not calculated. Also, the calculated MAE of bulk CrI 3 agree reasonably well with the experimental value reported in Ref. [3] (see Table 1). Since the MAEs of the CrI 3 multilayers are large and comparable to heavy element magnetic alloys such as FePt and CoPt [48], which have the largest MAEs among magnetic transition metal alloys and compounds, these CrI 3 multilayers could have valuable applications in high density data storage devices. Now we consider the ferromagnetic transition temperature (T c ) of the CrI 3 systems, which is another important parameter for a magnetic material. A rough estimation of T c could be made within the meanfield approximation (MFA), i.e., k B T m c = 1 3 J 0 where J 0 = j J 0j and the present cases J 0 = 3J 1 + 6J 2 + 3J 3 + J z . [49] The estimated ferromagnetic transition temperatures (T m c ) are listed in Table 2. It is clear from Table 2 that the estimated T m c and measured T c agree in trend, i.e., the magnetic transition temperature increases monotonically from ML, to BL, to TL and to bulk CrI 3 . Nonetheless, the estimated T m c values are much higher. In particular, for ML CrI 3 , the estimated T m c of 109 K is more than two times larger than the experimental T c 45 K for ML CrI 3 , and in bulk CrI 3 , the T m c of 191 K is about three times larger than the measured T c 61 K. [2]. The reason for this is that the MFA neglects transverse spin fluctuations [50,51], which is especially strong in lowdimensional materials. Furthermore, the MFA violates Mermin-Wagner theorem [22] which says that longrange magnetic order at a finite temperature cannot exist in an isotropic 2D Heisenberg magnet. However, the out-of-plane MAE in 2D materials can establish long-range magnetic orders at finite temperatures by opening a gap in the spin wave spectrum which would suppress thermal spin fluctuations [50,51]. Therefore, a better estimation of the T c would be based on the spin wave theory (SWT) with the out-of-plane MAE taken into account. This leads to an approximate expression of T c as k B T SW T c = πJ 1 /2log(1 + 2πJ 1 /∆ 0 ) where ∆ 0 is the spin wave gap [17], which can be roughly set to ∆ 0 = ∆E ma . The T SW T c values evaluated using the calculated J 1 and ∆E ma are listed in Table 2. It can be seen that compared with T m c , the T SW T c values agree better with the corresponding experimental values although they are still too high ( Table 2).

Electronic structure
Better understanding of the electronic, magnetic and magneto-optical properties of the materials can be obtained from the detailed analysis of the electronic band structure. Therefore, fully relativistic electronic band structure of bulk and ML CrI 3 structures are displayed in Fig. 3, while that of the BL and TL are displayed in Fig. S2 in the supplementary information (SI). Notably, these figures show that all the multilayer structures are almost direct bandgap semiconductors with both valance band maximum (VBM) and conduction band minimum (CBM) located at Γ point. In contrast, the bulk structure is an indirect bandgap semiconductor with the VBM sitting at the Γ point and the CBM located at the T point. The obtained semiconducting band gaps for all the structures are listed in Table 1, which shows that the band gap increases slightly as bulk CrI 3 is thinned down to TL, to BL and finally to ML.
To analyze the spin characters of the bands at the VBM and CBM, we present the spin-polarized scalar-relativistic band structures of all the structures considered here in Fig. S3 in the SI. Figure S3 indicates that both the lower conduction bands and top valence bands of all the structures are of purely spin-up character. Interestingly, this shows that all the structures are single-spin semiconductors and thus would be promising candidates for efficient spintronic and spin photovoltaic devices (see Ref. [52] and references therein). We also calculate the total as well as site, orbital projected spin-up and spin-down densities of states (DOS) for all the CrI 3 structures. The calculated DOS spectra of bulk and monolayer CrI 3 are displayed in Figs. 4 and 5, respectively, as representatives. Figure 4 shows that the contributions to the valance and conduction bands in the energy ranges of -4.5 to -0.3 eV and 0.4 eV to 2.5 eV, respectively, come mostly from the Cr d orbitals with significant contributions from the I p orbitals as well. The contributions from both the Cr and I confirm the strong hybridization between Cr d and I p orbitals in CrI 3 structures. The spin-up valance band states are predominately arise from the Cr d orbitals, along with minor contribution with the I p orbitals at the valence band edge. Conversely, the DOS states in spindown valance band region are mostly contributed from I p orbitals. The states in the upper valance band In case of conduction bands ranging from 0.4 eV to 1.2 eV is mainly because of spin-up Cr d xz,yz orbitals. From these one could say that the band gap in the CrI 3 structures arise due to the crystal-field splitting of spin-up Cr d xy,x 2 −y 2 and d xz,yz bands. The spindown conduction bands appear only above 1.3 eV with pronounced peaks made up of mainly Cr d xy,x 2 −y 2 and d z 2 orbitals (Fig. 4). Furthermore, the calculated total as well as site, orbital projected spin-up and spin-down DOS spectra of ML CrI 3 (Fig. 5) are rather similar to that of bulk CrI 3 (Fig. 4). The main difference comes from the contributions of the I p orbital which is enhanced in the higher energy region in the ML.
As mentioned above, the insulating band gap in the CrI 3 structures arises due to the crystal-field splitting of Cr d xy,x 2 −y 2 and d xz,yz states. The symmetries of these states near the band gap edges are important to the magneto-crystalline anisotropy in the presence of SOC [24,53], according to the perturbation theory. In particular, the SOC ma- trix elements for the d orbital of d xz |H SO |d yz and d x 2 −y 2 |H SO |d xy would favor the perpendicular magnetic anisotropy, while other elements of d x 2 −y 2 |H SO |d yz , d xy |H SO |d yz and d z 2 |H SO |d yz would prefer an in-plane magnetic anisotropy [54]. The magnitude ratio of these d-orbital SOC matrix elements are d xz |H SO |d yz 2 : d x 2 −y 2 |H SO |d xy 2 : d x 2 −y 2 |H SO |d yz 2 : d xy |H SO |d yz 2 : d z 2 |H SO |d yz 2 = 1 : 4 : 1 : 1 : 3 [54]. Figures 4 and 5 show that in the upper valence and lower conduction band region, the d xz,yz and d x 2 −y 2 ,xy orbital-decomposed DOSs are large and this would result in large d xz |H SO |d yz and d x 2 −y 2 |H SO |d xy and hence large contributions to the perpendicular anisotropy. In contrast, the d z 2 DOS is almost zero in the lower conduction band region of 0.4∼1.5 eV in both bulk and ML CrI 3 , and thus this would give rise to the diminishing SOC matrix element of d z 2 |H SO |d yz 2 which favors an in-plane magnetization. All these together would lead to the magnetocrystalline anisotropy energes in the CrI 3 structures (see Table 1) which strongly favor the out-of-plane magnetization.

Optical conductivity
Now let us turn our attention to the calculated optical conductivity, which is the ingredient for evaluating the magneto-optical effects, as explained already in Sec. 2. The calculated optical conductivities of the investigated structures are displayed in Fig. 6 for bulk and ML CrI 3 as well as in Fig. S4 in the SI for BL and TL CrI 3 . Since the calculated optical conductivities for all the considered structures (see Figs. 6 and S4) are rather similar due to the weak van der Waals interlayer interaction, here we discuss only bulk and ML CrI 3 as examples. Notably, Fig. 6 shows that in both bulk and ML CrI 3 , the absorptive and dispersive parts of the diagonal optical conductivity elements for an in-plane electric polarization (E||ab) (σ xx ) differ significantly from that for the out-of-plane electric polarization (E||c) (σ zz ) above the absorption edge of about 1.5 eV. For example, the absorptive part σ 1 xx for E||ab of bulk CrI 3 is significantly larger than that (σ 1 zz ) for E||c in the energy range from 1.5 to 4.4 eV, while it becomes smaller than σ 1 zz for E||c above 4.4 eV. A similar profile is observed in ML CrI 3 in the energy region of 1.0 ∼ 4.3 eV, where higher value is found for E||ab compared to that for E||c, while above 4.3 eV σ 1 zz is higher than σ 1 xx . This manifests that these structures are highly optically anisotropic. The observed strong optical anisotropy is quite common in low dimensional materials and can be explained in terms of the orbital-decomposed DOS spectra reported in the previous subsection. In particular, Figure 5 The calculated real (σ 1 xy ) and imaginary (σ 2 xy ) parts of the off-diagonal optical conductivity element are displayed in Fig. 6(e) and Fig. 6(f) for bulk and ML CrI 3 , respectively. It is quite evident from these figures that the off-diagonal element spectra of bulk and ML CrI 3 are again rather similar, due to the weak van der Waals interlayer coupling in bulk CrI 3 . For example, the σ 1 xy shows a large positive peak at 4.2 eV for the bulk and at 4.4 eV for the ML. In addition, there is another prominant peak at ∼3.5 eV for bulk and ML structures. In the σ 2 xy spectra, there is a large positive peak in the vicinity of 4.8 eV for both structures. There are also negative peaks in both σ 1 xy and σ 2 xy spectra. For example, a negative peak occurs in the vicinity of 1.3 eV in the σ 1 xy spectrum and at 3.3 eV in the σ 2 xy spectrum for bulk and ML structures. The absorptive parts of the optical conductivity  Figure 6. are, respectively, the experimental differential reflection spectra (in abitrary units) for bulk and ML CrI 3 [7] which are proprotional to σ 1 xx . The red dashed curve in (a) is the σ 1 xx of bulk CrI 3 derived from the experimental refraction index and extinction coefficent. [2] The red diamond in (f) denotes the σ 2 xy of ML CrI 3 derived from the experimental photoluminescence [7] (see text). elements (σ 1 xx , σ 1 zz , σ 2 xy and σ 1 ± ) are directly related to the dipole allowed interband transitions, as Eqs. (3), (4) and (10) indicate. Therefore, we further analyze the origin of the main features in the σ 1 xx , σ 1 zz and σ 2 xy spectra by considering the symmetries of the involved band states and hence the dipole selection rules (see Supplementary Note B in the SI). Since the spectra of σ 1 xx , σ 1 zz and σ 2 xy for all the considered systems are rather similar, as mentioned above, here we perform the symmetry analysis for ML CrI 3 only (see the SI). The scalar-relativistic and relativistic band structures of ML CrI 3 with the symmetries of band states at Γ-point labelled, are displayed in Fig. S6 and Fig.  7, respectively. We then assign the main peaks of the σ 1 xx and σ 1 zz spectra [ Fig. 6(b)] to the interband transitions at the Γ-point using the dipole selection rules (see Table S3 in the SI), as presented in Fig.  S6. For example, we could relate the A 1 peak at ∼1.5 eV in σ 1 xx in Fig. 6(b) to the interband transition from the Γ + 3 state at -0.5 eV of the spin-down valence band to the conduction band state Γ − 3 at 1.2 eV. Of course, in addition to these dipole-allowed interband transitions at the Γ-point, there are also contributions from different interband transitions at other k-points. It should be noted that, without the SOC, the Γ + 3 and Γ − 3 band states are doubly degenerate (Fig. S6), and the absorption rates for left-and right-handed polarized lights are the same, thereby resulting in zero contribution to σ xy . In the presence of the SOC, the band states split (see Fig. 7) and this gives rise to optical magnetic circular dichroism. For example, we could relate the peak P 1 at 1.2 eV in the σ 2 xy of Fig.  6(f) to the interband transition from the Γ − 6 occupied states at -0.1 eV to the Γ + 4 unoccupied state at ∼1.1 eV (Fig. 7). For comparison with the available experimental data, we plot in Figs. 6 (a) and (b) the experimental differential reflection spectra (in abitrary units) of bulk and ML CrI 3 , respectively [7], which are proportional to σ 1 xx . Moreover, we also display the σ 1 xx of bulk CrI 3 derived from the experimental refraction index and extinction coefficent [2] in Fig. 6(a). It can be seen from Fig. 6(a) that twe peaks (A 2 and A 3 in the experimental σ 1 xx of bulk CrI 3 [2] agree well with the theoretical ones in both shape and position, although the magnitude of the experimental spectrum is slightly smaller. Moreover, three peaks (A 1 , A 2 and A 3 ) in the experimental differential reflection spectrum of bulk CrI 3 [7] also agree very well with the theoretical ones [see Fig. 6(a)]. This indicates that our GGA + U calculations could describe the optical properties of bulk CrI 3 quite well. Figure  6(b) shows that three peaks (A 1 , A 2 and A 3 ) in the experimental differential reflection spectrum of ML CrI 3 [7] are also quite well reproduced by our GGA + U calculations. In particular, the position and shape of peak A 1 in theoretical and experimental spectra are in good agreement, although theoretical peaks A 2 and A 3 are red-shifted by about 0.3 eV relative to the experimental ones. Assuming photoluminescence (PL) intensity (I) is proportional to photoabsorption (σ), PL circular polarization ρ = (I i.e., σ 2 xy = ρσ 1 xx . In Fig. 6(f), the red diamond denotes the σ 2 xy estimated this way based on the calculated σ 1 xx and experimental ρ for ML CrI 3 [7], which agree well with peak P 1 in the theoretical σ 2 xy spectrum. Table 1 indicates that the calculated band gap of bulk CrI 3 is about 50 % smaller than the experimental value [3]. As usual, we could attribute this significant discrepancy to the known underestimation of the band gaps by the GGA functional where many-body effects especially quasiparticle self-energy correction are neglected. Since Eqs. (3) and (4) indicate that correct band gaps would be important for obtaining accurate optical conductivity and hence magneto-optical effects, we further perform the band structure calculations using the hybrid Heyd-Scuseria-Ernzerhof (HSE) functional [55,56] which is known to produce improved band gaps for semiconductors (see Supplementary Note C in the SI for the details). The HSE band structures and band gaps are presented in Fig. S7 and Table S4 in the SI, respectively. Indeed, the HSE band gap of bulk CrI 3 agrees well with the experimental value [3]. Therefore, we also calculate the optical and magneto-optical properties of bulk and few-layer CrI 3 from the GGA + U band structures within the scissors correction (SC) scheme [57] using the differences between the HSE and GGA band gaps (Table S4). The optical conductivity spectra calculated with the SC are displayed in Fig. S8 in the SI. Surprisingly, the optical conductivity spectra of bulk and ML CrI 3 calculated with the SC are in worse agreement with the corresponding experimental values [2,7] [see Figs. S8(a), S8(b) and S8(j)]. This could be explained as follows. We notice that bulk CrI 3 has an indirect band gap (see Fig. 3). Consequently, the discrepancy in the band gap between theory and experiment in bulk CrI 3 could be attributed to the fact that the experimental band gap is the optical one due to direct interband transitions whose energy should be larger than the indirect band gap (see Fig. 7).

Magneto-optical Kerr and Faraday effects
Finally, let us move onto the MO Kerr and Faraday rotational angles calculated from the obtained optical conductivity spectra. To calculate the Kerr angles of the few-layer structures, we also need to know the dielectric constant of the substrate [see Eq. (8)]. Here we assume that the substrate is SiO 2 with dielectric constant ε s = 2.25 for all CrI 3 multilayers since SiO 2 was used in the MOKE experiments [2]. The calculated complex Kerr and Faraday rotation angles as a function of photon energy are presented, respectively, in Figs. 8 and 9 for all the investigated structures. Figures 8 and 9 show that the MO spectra for all the investigated structures except the bulk, are negligible below the band gap of ∼1.0 eV. Nonetheless, they all become pronounced above 1.0 eV and oscillate as the photon energy increases. We notice that the Kerr rotation (θ K ) and Kerr ellipticity (ε K ) of all the structures (Fig. 8) are, respectively, similar to that of the real part (σ 1 xy ) and imaginary part (σ 2 xy ) of the offdiagonal conductivity element (Fig. 6). This could be expected because the Kerr effect and the off-diagonal conductivity element are connected through Eqs. (7) and (8). In fact Eqs. 7 and 8 indicate that the complex Kerr rotation angle would be proportional to the offdiagonal optical conductivity of σ xy if the conductivity σ xx of the sample and substrate is roughly a constant of photon energy. For the CrI 3 multilayers, the dielectric constant of ε = 2.25 of the SiO 2 substrate is a constant and thus the Kerr rotation angle is proportional to the σ xy . For bulk CrI 3 , the Kerr rotation could become large if the σ xx , which is in the denominator of Eq. (7), becomes very small. This explains why the Kerr angles of the bulk are still pronounced below 0.5 eV [ Fig. 8(a)].
The overall behavior of the MO spectra for all the investigated structures (Figs. 8 and 9) is rather similar, which could be attributed to the weakness of the van der Waals interlayer interaction. Notably, the maximum Kerr rotation angles for all the systems are large. Specifically, the maximum θ K values for bulk, ML, BL and TL CrI 3 are, respectively, 1.5 • , 0.6 • , 0.8 • and 1.2 • at ∼1.3 eV. Kerr ellipticity shows a maximum value of 1.3 • at 1.6 eV for the bulk as well as 0.6 • , 0.9 • and 1.2 • at 3.3 eV for the ML, BL and TL, respectively. As for the calculated optical conductivities, we notice that the negative peaks also exist in both Kerr rotation and ellipticity spectra. The negative maximum of the Kerr rotation angle occurs at 4.3 eV for the bulk (-2.1 • ) and at 4. . Kerr rotation (θ K ) and ellipticity (ε K ) spectra of (a) bulk, (b) ML, (c) BL, and (d) TL CrI 3 . The filled red circle in each panel indicates the θ K value from the recent experiment [2] that the considered structures would have valuable applications in magneto-optical nanodevices.
To access the application potential of these CrI 3 materials for applications, let us now compare their MO properties with several well-known MO materials such as 3d transition metal alloys and compound semiconductors. Among the transition metal alloys, manganese-based pnictides are known to have remarkable MO properties. In particular, MnBi thin films were reported to have a large Kerr rotation angle of 2.3 • at 1.84 eV. [46,58] Pt alloys also possess large Kerr rotations in the range of 0.4-0.5 • in FePt and Co 2 Pt [43] and PtMnSb [59]. In these materials, the strong SOC of heavy Pt atoms was shown to play an important role [60]. For comparison, we note that 3d transition metal multilayers were found to have rather small Kerr rotation angles of 0.06 • at 2.2 eV for Fe/Cu multilayer, of 0.16 • at 2.5 eV for Fe/Au multilayer [44,61] and also of 0.06 • at 1.96 eV for Co/Au multilayer [62]. Among semiconductor MO materials, diluted magnetic semiconductors Ga 1−x Mn x As were reported to show significant Kerr rotation angles of ∼0.4 • near of 1.80 eV [63]. Y 3 Fe 5 O 12 also exhibits a significant Kerr rotation of 0.23 • at 2.95 eV [64]. In strong contrast, the measured Kerr rotation angles of bulk and atomically thin films of Cr 2 Ge 2 Te 6 are much smaller [1]. For example, the experimental Kerr angle at 0.8 eV is only ∼0.14 • for the bulk, and ranges from 0.0007 • in bilayer to 0.002 • in trilayer [1]. Overall, the Kerr rotation angles predicted for bulk and multilayer CrI 3 here are comparable to these important MO metals and semiconductors. Thus, because of their excellent MO properties, CrI 3 structures could have promising applications in, e.g., MO nanosensors and high density MO data-storage devices.
The calculated complex Faraday rotation angles for all the investigated structures are presented in Fig. 9. The main features of the Faraday rotation spectra of the considered structures resemble that of the Kerr rotation spectra (see Fig. 8). Notably, the large Faraday rotation maximum occurs at the photon energy of ∼2.3 eV for the bulk (36 • /µm), of ∼1.4 eV for the ML (50 • /µm), of 2.3 eV for the BL (75 • /µm) and TL (108 • /µm). Faraday ellipticity shows the large maximum near 3.3 eV for the bulk (48 • /µm), the ML (80 • /µm), the BL (130 • /µm) and the TL (170 • /µm). Also we notice that the large negative maximum exists in both Faraday rotation and ellipticity spectra. For example, the negative maximum of the Faraday rotation occurs at 4.3 eV for the bulk (60 • /µm) and the TL (-194 • /µm) as well as at 4.7 eV for the ML (-85 • /µm) and the BL (-155 • /µm). Faraday ellipticity shows the large negative maximum near 4.8 eV for all the structures with ε F = -47 • /µm, -85 • /µm, -148 • /µm and -155 • /µm for the bulk, ML, BL and TL, respectively. For comparison, we notice that the MnBi films are known to possess the largest Faraday rotation angles of ∼80 • /µm at 1.77 eV [46,58]. Ferromagnetic semiconductor Y 3 Fe 5 O 12 show only small Faraday rotation angles of ∼ 0.19 • /µm at 2.07 eV [65]. However, the replacement of Y with Bi considerably enhances the Faraday rotations to ∼35.0 • /µm at 2.76 eV [66], thus making Bi 3 Fe 5 O 12 an outstanding MO material. Clearly, the calculated Faraday angles of the CrI 3 multilayers are high compared to that of prominent bulk MO metals such as manganese pnictides.
The Kerr rotation angles θ K of bulk and atomically thin film CrI 3 have recently been measured at photon energy of 1.96 eV [2] and they are found to be large (see Fig. 8). If the calculated Kerr rotation spectra are red-shifted by merely 0.2 eV, the measured and calculated θ K values of bulk and ML CrI 3 would be in very good agreement, although the experimental θ K values of BL and TL CrI 3 are much larger than the corresponding theoretical θ K values (Fig. 8). Furthermore, the measured Faraday rotation angles θ F of bulk CrI 3 [4] are almost in perfect agreement with our theoretical prediction [see Fig. 9(a)]. We notice that in contrast, the photon energies of the experimental Faraday angles of bulk CrI 3 [4] differ significantly from that of the theoretical Faraday angles calculated with the scissors correction using the HSE band gap [see Fig. S10(a) in the SI].

Conclusion
Summarizing, we have systematically studied magnetism, electronic structure, optical and magnetooptical properties of bulk and atomically thin films of CrI 3 by performing extensive GGA + U calculations. Firstly, we find that both intralayer and interlayer magnetic couplings are strongly ferromagnetic in all the considered CrI 3 structures. Ferromagnetic ordering temperatures estimated based on calculated exchange coupling parameters agree quite well with the measured ones.
Secondly, we reveal that all the structures have a large out-of-plane MAE of ∼0.5 meV/Cr, in contrast to the significantly thickness-dependent MAE in atomically thin films of Cr 2 Ge 2 Te 6 . [26] These large MAEs suppress transverse spin fluctuations in atomically thin CrI 3 films and thus stabilize long-range magnetic orders at finite temperatures down to the ML limit.
Thirdly, all the structures also exhibit strong MO Kerr and Faraday effects. In particular, in these CrI 3 structures, MO Kerr rotations up to ∼2.6 • are found, being comparable to famous MO materials such as ferromagnetic semiconductor Y 3 Fe 5 O 12 and metal PbMnSb. Calculated MO Faraday rotation angles are also large, being up to ∼200 • /µm, and they are comparable to that of best known MO semiconductors such as Y 3 Fe 5 O 12 and Bi 3 Fe 5 O 12 . Calculated optical conductivity spectra, MO Kerr and Faraday rotation angles agree quite well with available experimental data. In BL CrI 3 , the interlayer magnetic coupling is found to be layer stacking dependent while the optical and magneto-optical properties are not. Therefore, further experimental determination of the structure of BL CrI 3 is needed to resolve why the antiferromagnetic state instead of the predicted ferromagnetic state, is observed. Fourthly, all the structures are found to be single-spin ferromagnetic semiconductors. Finally, calculated MAEs, optical and magneto-optical conductivity spectra of these structures are analyzed in terms of their underlying electronic band structures. Our findings of large outof-plane MAEs and strong MO effects in these singlespin ferromagnetic semiconducting CrI 3 ultrathin films suggest that they will find promising applications in high density semiconductor MO and low-power spintronic nanodevices as well as efficient spinphotovoltaics.