Magnetism and magneto-optical effects in bulk and few-layer CrI3: a theoretical GGA + U study

Vijay Kumar Gudelli; Guang-Yu Guo

doi:10.1088/1367-2630/ab1ae9

1. Introduction

The recent discovery of intrinsic ferromagnetism in atomically thin films of semiconductors Cr₂Ge₂Te₆ [1] and CrI₃ [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 and to magneto-optical (MO) devices based on 2D materials. Therefore, these atomically thin magnetic materials are currently subject to intensive investigations. In fact, bulk CrI₃ has already been well studied experimentally in the past [3–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₃. Excitingly, manipulation of 2D magnetism (e.g. switching of magnetization direction and tuning antiferromagnetic to ferromagnetic transition) in BL CrI₃ 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₃ and other van der Waals (vdW) materials [11, 12]. Unprecedented control of spin and valley pseudospin in ultrathin CrI₃ and WSe₂ 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 first-principles density functional study of the magnetic, electronic, optical and MO properties of multilayers (monolayer (ML), bilayer (BL) and trilayer (TL)) as well as bulk CrI₃. Specifically, we focus the magnetic interactions, magnetic anisotropy energy (MAE) and MO effects in these magnetic structures. Although magnetic exchange coupling and magnetic anisotropy in ML CrI₃ have been theoretically investigated by several groups based on density functional calculations [14–17], they have not been studied in BL, TL and bulk CrI₃. Furthermore, the interlayer exchange coupling in bulk and TL CrI₃ has not been addressed although that in the BL is presently under intensive investigations [18–21].

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 out-of-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₃ [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-of-plane magnetization depending on the band structure of the material. Although the MCE of ML CrI₃ 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–26]. In fact, it was found recently that although the MCE in ML Cr₂Ge₂Te₆ favors the out-of-plane magnetization, the MDE is larger than the MCE such that ML Cr₂Ge₂Te₆ 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₃. Therefore, in this paper we present both calculated MCE and MDE for all the considered CrI₃ structures.

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 figure 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₂Ge₂Te₆ and CrI₃ 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 MO 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 MO devices down to subnanometer scale. Therefore, here we carry out a systematic first-principles density functional study of the optical and MO properties of bulk and multilayer CrI₃. Indeed, our study reveals that bulk and multilayer CrI₃ 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₃ semiconductor structures will provide an interesting material platform for further studies of novel MO phenomena and technological applications.

**Figure 1.** (a) Schematic illustration of magneto-optical Kerr and Faraday effects in CrI₃. (b) ML CrI₃ top-view and (c) Crystalline structure of TL CrI₃ along with the primitive bulk structure (solid black lines) and the blue dashed rectangle indicates ML and BL CrI₃ (side-view).
Download figure:
Standard image High-resolution image

2. Computational methods

In this paper, we study the electronic, magnetic, and MO properties of bulk and multilayer CrI₃ structures. Bulk CrI₃ forms a layered structure with MLs separated by the van der Waals gap (figure 1(c)). Each CrI₃ ML consists of edge-sharing CrI₆ octahedra forming a planar network with Cr atoms in a honeycomb lattice (figure 1(b)) and thus there are two formula units (fu) per lateral unit cell. These MLs are then stacked in an ABC sequence, thus resulting in a rhombohedral crystal with $R\bar{3}$ symmetry and two fu 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₃ structures, the hexagonal unit cell (figure 1(b)) with the experimental lattice constants is considered as the lateral unit cell. For BL and TL CrI₃ structures (see figure 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⁵ 4s¹, 5s² 5p⁵, 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₃, 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. [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₂Ge₂Te₆ 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₂Ge₂Te₆ and CrI₃, 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₃ systems and find that the calculated physical quantities do not depend significantly on the U value used (see supplementary note A, table S1 and figure S1 in the supplementary information (SI) available online at stacks.iop.org/NJP/21/053012/mmedia). 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, labeled as AF-Néel, AF-zigzag, and AF-stripe in figure 2. For bulk, BL and TL CrI₃, we also consider interlayer FM and AF magnetic configurations with intralayer FM CrI₃ MLs. To understand the magnetic interactions and also to estimate magnetic ordering temperature (T_c) for bulk and multilayers CrI₃, we calculate the exchange coupling parameters by mapping the calculated total energies of these magnetic configurations (see figure 2) onto the classical Heisenberg Hamiltonian, ${\text{}}E={E}_{0}-{\sum }_{i,j}{J}_{{ij}}{\hat{{\bf{e}}}}_{i}\ {\boldsymbol{\cdot }}\ {\hat{{\bf{e}}}}_{j},$ where E₀ donates the nonmagnetic ground state energy; J_ij the exchange coupling parameter between sites i and j; ${\hat{{\bf{e}}}}_{i}$ , the unit vector representing the direction of the magnetic moment on site i. For one CrI₃, we obtain a set of four linear equations of J₁, J₂ and J₃: ${E}_{{\rm{FM}}}={E}_{0}-3{J}_{1}-6{J}_{2}-3{J}_{3}$ , E_AF-Néel = ${E}_{0}+3{J}_{1}-6{J}_{2}+3{J}_{3}$ , ${E}_{{\rm{AF}} \mbox{-} {\rm{zigzag}}}={E}_{0}-{J}_{1}+2{J}_{2}+3{J}_{3}$ and ${E}_{{\rm{AF}} \mbox{-} {\rm{stripe}}}={E}_{0}+{J}_{1}+2{J}_{2}-3{J}_{3}$ , which can be solved to estimate J₁, J₂ and J₃. 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 ${J}_{z}=({E}_{{\rm{AF}}}-{E}_{{\rm{FM}}})/4$ for the BL, ${J}_{z}=({E}_{{\rm{AF}}}-{E}_{{\rm{FM}}})/8$ for the TL and ${J}_{z}=({E}_{{\rm{FiM}}}-{E}_{{\rm{FM}}})/8$ for the bulk where E_FiM denotes the total energy of the ferrimagnetic configuration with the Cr magnetic moments on one of the three layers in the hexagonal unit cell flipped.

**Figure 2.** Four considered intralayer magnetic configurations: (a) FM, (b) AF-Néel, (c) AF-zigzag, (d) AF-stripe types. The magnetic Cr atoms are shown with red (blue) arrows indicating up (down) spins. Intralayer exchange coupling parameters J₁, J₂, and J₃ are indicated by green arrows.
Download figure:
Standard image High-resolution image

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-of-plane 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₃, 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]

$\begin{eqnarray}&&{E}_{d}=\displaystyle \sum _{{qq}^{\prime} }\displaystyle \frac{2{m}_{q}{m}_{q^{\prime} }}{{c}^{2}}{M}_{{qq}^{\prime} },\end{eqnarray} \tag{ 1 }$

where the speed of light c = 274.072 and the so-called magnetic dipolar Madelung constant

$\begin{eqnarray}&&{M}_{{qq}^{\prime} }=\displaystyle \sum _{{\bf{R}}}^{^{\prime} }\displaystyle \frac{1}{| {\bf{R}}+{\bf{q}}+{\bf{q}}^{\prime} {| }^{3}}\left\{1-3\displaystyle \frac{{\left[({\bf{R}}+{\bf{q}}+{\bf{q}}^{\prime} ).{\hat{m}}_{q}\right]}^{2}}{| {\bf{R}}+{\bf{q}}+{\bf{q}}^{\prime} {| }^{2}}\right\},\end{eqnarray} \tag{ 2 }$

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 ${\mu }_{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-of-plane magnetizations. In a 2D material, since all the R and q are in-plane, the second term in equation (2) would be zero for the out-of-plane magnetization, thereby resulting in the positive ${M}_{{qq}^{\prime} }$ while for an in-plane magnetization the ${M}_{{qq}^{\prime} }$ 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

$\begin{eqnarray}&&{\sigma }_{{aa}}^{1}(\omega )=\displaystyle \frac{\pi {e}^{2}}{{\hslash }\omega {m}^{2}}\displaystyle \sum _{i,j}{\int }_{{\rm{BZ}}}\displaystyle \frac{{\rm{d}}{\bf{k}}}{{\left(2\pi \right)}^{3}}| {p}_{{ij}}^{a}{| }^{2}\delta ({\epsilon }_{{\bf{k}}j}-{\epsilon }_{{\bf{k}}i}-{\hslash }\omega ),\end{eqnarray} \tag{ 3 }$

$\begin{eqnarray}&&{\sigma }_{{xy}}^{2}(\omega )=\displaystyle \frac{\pi {e}^{2}}{{\hslash }\omega {m}^{2}}\displaystyle \sum _{i,j}{\int }_{{\rm{BZ}}}\displaystyle \frac{{\rm{d}}{\bf{k}}}{{\left(2\pi \right)}^{3}}\mathrm{Im}[{p}_{{ij}}^{x}{p}_{{ji}}^{y}]\delta ({\epsilon }_{{\bf{k}}j}-{\epsilon }_{{\bf{k}}i}-{\hslash }\omega ),\end{eqnarray} \tag{ 4 }$

where ℏω is the photon energy, and ${\epsilon }_{{\bf{k}}i}$ is the ith band energy at point k. Summations i and j are over the occupied and unoccupied bands, respectively. Dipole matrix elements ${p}_{{ij}}^{a}=\langle {\bf{k}}{\text{}}j| {\hat{p}}_{a}| {\bf{k}}i\rangle$ where ${\hat{p}}_{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 [40] and references therein). The dispersive parts of the optical conductivity elements can be obtained using the Kramers–Kroing relations

$\begin{eqnarray}&&{\sigma }_{{aa}}^{2}(\omega )=-\displaystyle \frac{2\omega }{\pi }P{\int }_{0}^{\infty }\displaystyle \frac{{\sigma }_{{aa}}^{1}(\omega ^{\prime} )}{{\omega }^{{\prime} 2}-{\omega }^{2}}{\rm{d}}{\omega }^{{\prime} },\end{eqnarray} \tag{ 5 }$

$\begin{eqnarray}&&{\sigma }_{{xy}}^{1}(\omega )=\displaystyle \frac{2}{\pi }P{\int }_{0}^{\infty }\displaystyle \frac{{\omega }^{{\prime} }{\sigma }_{{xy}}^{2}(\omega ^{\prime} )}{{\omega }^{{\prime} 2}-{\omega }^{2}}{\rm{d}}{\omega }^{{\prime} }\end{eqnarray} \tag{ 6 }$

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.1 eV (see, e.g. figures 1(a) and (b) in [41]), which is thus used in this paper.

For a bulk magnetic material, the complex polar Kerr rotation angle is given by [42, 43]

$\begin{eqnarray}&&{\theta }_{{\rm{K}}}+{\rm{i}}{\epsilon }_{{\rm{K}}}=\displaystyle \frac{-{\sigma }_{{xy}}}{{\sigma }_{{xx}}\sqrt{1+{\rm{i}}(4\pi /\omega ){\sigma }_{{xx}}}}.\end{eqnarray} \tag{ 7 }$

However, for a magnetic thin film on a nonmagnetic substrate, the complex polar Kerr rotation angle is given by [44, 45]

$\begin{eqnarray}&&{\theta }_{{\rm{K}}}+{\rm{i}}{\epsilon }_{{\rm{K}}}={\rm{i}}\displaystyle \frac{2\omega d}{c}\displaystyle \frac{{\sigma }_{{xy}}}{{\sigma }_{{xx}}^{s}}=\displaystyle \frac{8\pi d}{c}\displaystyle \frac{{\sigma }_{{xy}}}{(1-{\varepsilon }_{{xx}}^{s})},\end{eqnarray} \tag{ 8 }$

where d stands for the thickness of the magnetic layer and ${\varepsilon }_{{xx}}^{s}({\sigma }_{{xx}}^{s})$ the diagonal part of the dielectric constant (optical conductivity) of the substrate. Experimentally, atomically thin CrI₃ films were prepared on SiO₂/Si [2, 7, 13]. Thus, the dielectric constant ( ${\varepsilon }_{{xx}}^{s}=2.25$ ) of bulk SiO₂ is used here.

Similarly, the complex Faraday rotation angle for a thin film can be written as [46]

$\begin{eqnarray}&&{\theta }_{{\rm{F}}}+{\rm{i}}{\epsilon }_{{\rm{F}}}=\displaystyle \frac{\omega d}{2c}({n}_{+}-{n}_{-}),\end{eqnarray} \tag{ 9 }$

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}_{\pm }^{2}={\varepsilon }_{\pm }=1+\tfrac{4\pi {\rm{i}}}{\omega }{\sigma }_{\pm }=1+\tfrac{4\pi {\rm{i}}}{\omega }({\sigma }_{{xx}}\pm {\rm{i}}{\sigma }_{{xy}})$ . Here the real parts of the optical conductivity σ_±can be written as

$\begin{eqnarray}&&{\sigma }_{\pm }^{1}(\omega )=\displaystyle \frac{\pi {e}^{2}}{{\hslash }\omega {m}^{2}}\displaystyle \sum _{{i},j}{\int }_{{\rm{BZ}}}\displaystyle \frac{{\rm{d}}{\bf{k}}}{{\left(2\pi \right)}^{3}}| {{\rm{\Pi }}}_{{ij}}^{\pm }{| }^{2}\delta ({\epsilon }_{{\bf{k}}j}-{\epsilon }_{{\bf{k}}{i}}-{\hslash }\omega ),\end{eqnarray} \tag{ 10 }$

where ${{\rm{\Pi }}}_{{ij}}^{\pm }=\langle {\bf{k}}{\text{}}j| \tfrac{1}{\sqrt{2}}({\hat{p}}_{x}\pm {\text{}}{\rm{i}}{\hat{p}}_{y})| {\bf{k}}i\rangle$ . Clearly, ${\sigma }_{{xy}}=\tfrac{1}{2{\rm{i}}}({\sigma }_{+}-{\sigma }_{-})$ , 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 MO effects.

3. Results and discussion

3.1. Magnetic properties

3.1.1. Magnetic structure and exchange coupling

To find the ground state magnetic structure and also evaluate exchange coupling among the magnetic Cr atoms in all the considered CrI₃ 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 figure 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 figure 2(e)) structures for bulk, BL and TL CrI₃, 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₃ in the FM state only in table 1. The calculated Cr spin magnetic moment in all the CrI₃ 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/fu for all the structures.

Table 1. Total spin magnetic moment (m_s^t), atomic (averaged) spin (m_s^Cr, m_s^I) and orbital (m_o^Cr, m_o^I) magnetic moments as well as band gap (E_g), magnetocrystalline anisotropy energy (ΔE_b), dipolar anisotropy energy (ΔE_d) and (total) magnetic anisotropy energy (ΔE_ma = ΔE_b + ΔE_d) of bulk and few-layer CrI₃ in the ferromagnetic ground state (magnetization being perpendicular to the layers) calculated using the GGA + U scheme with the spin–orbit coupling included. Also listed are the experimental magnetic anisotropy energy ( ${\rm{\Delta }}{E}_{{\rm{ma}}}^{{\rm{\exp }}}$ ) and band gap (E_g^exp) for comparison.

	m_s^t	m_s^Cr (m_o^Cr)	m_s^I (m_o^I)	ΔE_b (ΔE_d)	ΔE_ma ( ${\rm{\Delta }}{E}_{{\rm{ma}}}^{{\rm{\exp }}}$ )	E_g (E_g^exp)
	(μ_B/fu)	(μ_B/atom)	(μ_B/atom)	(meV/fu)	(meV/fu)	(eV)
ML	6.00	3.20 (0.083)	−0.092 (−0.014)	0.678 (−0.152)	0.526	0.81
BL	6.00	3.21 (0.083)	−0.094 (−0.014)	0.620 (−0.152)	0.468	0.71
TL	6.00	3.21 (0.083)	−0.093 (−0.014)	0.586 (−0.152)	0.434	0.68
Bulk	6.00	3.21 (0.082)	−0.094 (−0.014)	0.545 (−0.064)	0.481 (0.26^a)	0.62 (1.2^a)

Note.

^aReference[3].

As mentioned in section 2, using the total energies of the considered magnetic configurations, we estimate first three near-neighbor exchange coupling parameters (J₁, J₂, J₃) in each ML in all the considered 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₁) is ferromagnetic (i.e. positive) in all the considered CrI₃ structures. The second near-neighbor coupling (J₂) is also ferromagnetic (positive), although the J₂ value is about 3.5 times smaller than J₁ for all the structures. Interestingly, table 2 shows that the magnitudes of both J₁ and J₂ increase monotonically with the number of layers, being consistent with the observed increase in magnetic transition temperature in these systems [2]. The third near-neighbor magnetic coupling (J₃) is also ferromagnetic (positive) in all the CrI₃ structures except ML CrI₃ where J₃ is negative (antiferromagnetic). Nonetheless, the calculated J₃ values are found to be 5 times smaller than J₁ and only about half of J₂ for all the structures. For the CrI₃ ML, we find that all the calculated exchange coupling parameters are in good quantitative agreement with the recent GGA calculations [14–16], although the J₁ 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 [17] and also that the J values reported previously [14–17] should be multiplied by a factor of S² = 9/4 in order to compare with the present J values.

Table 2. Intralayer first-, second- and third-neighbor exchange coupling parameters (J₁, J₂, J₃) as well as effective interlayer exchange coupling constant (J_z) in ML, BL, TL and bulk CrI₃, 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}_{c}^{m}$ ) and also within the spin wave theory (SWT) ( ${T}_{c}^{{\rm{SWT}}}$ ) (see the text) are also included. For comparison, the available experimental T_c values are also listed in brackets.

	J₁	J₂	J₃	J_z	${T}_{c}^{m}$	${T}_{c}^{{\rm{SWT}}}$ ( ${T}_{c}^{{\rm{\exp }}}$ )
	(meV)	(meV)	(meV)	(meV)	(K)	(K)
ML	6.91	1.48	−0.46	—	109	66 (45^a)
BL	7.59	2.30	1.02	2.97	165	69
TL	8.05	2.45	0.91	2.85	172	71
Bulk	8.08	2.89	1.62	2.95	191	73 (61^b)

Notes.

^aReference [2]. ^bReference [14].

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–5] on all the structures except BL CrI₃ 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₃ is nontrivial. Indeed, very recently this problem has already been investigated theoretically by several groups [18–21]. These theoretical studies all showed that the interlayer magnetic coupling in BL CrI₃ 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₃, 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₃ 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₃ are needed to resolve this interesting problem. To this end, we calculate the optical and MO properties of BL CrI₃ in both AB- and AB'-stacking orders. Nonetheless, we find that the MO effects of BL CrI₃ hardly depend on the stacking order, as presented in figures S4 and S5 in the supplementary information.

3.1.2. MAE and ferromagnetic transition temperature:

The calculated magnetic anisotropy energies ( $\bigtriangleup {E}_{{\rm{ma}}}$ ) for bulk and multilayers CrI₃ are presented in table 1. As mentioned before, the MAE consists of two competing contributions, namely, MCE ( $\bigtriangleup {E}_{b}$ ) due to the SOC effect on the band structure and magnetic dipolar (shape) anisotropy energy ( $\bigtriangleup {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 ( $\bigtriangleup {E}_{b}+\bigtriangleup {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₂Ge₂Te₆ multilayers where the MAE is much smaller and depends significantly on the film thickness [26]. In particular, for ML Cr₂Ge₂Te₆, 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₃ [2] but not in ML Cr₂Ge₂Te₆ [1]. These contrasting magnetic anisotropic properties between multilayers CrI₃ and Cr₂Ge₂Te₆ 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₃ is in good quantitative agreement with the theoretical value of 0.686 meV/fu reported recently in [14] in which the MDE was not calculated. Also, the calculated MAE of bulk CrI₃ agree reasonably well with the experimental value reported in [3] (see table 1). Interestingly, the MAEs of the CrI₃ 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. This appears to be somewhat surprising since the 3d transition metals such as Cr usually have small SOC. The large MAEs could be attributed to two factors, namely, the 2D and quasi-2D structures of these materials and also the large SOC of the heavy ligand ion (I) [49]. As a result, these CrI₃ multilayers could have valuable applications in high density data storage devices.

Now we consider the ferromagnetic transition temperature (T_c) of the CrI₃ systems, which is another important parameter for a magnetic material. A rough estimation of T_c could be made within the mean-field approximation (MFA), i.e. ${k}_{B}{T}_{c}^{m}=\tfrac{1}{3}{J}_{0}$ where ${J}_{0}={\sum }_{j}{J}_{0j}$ and the present cases ${J}_{0}=3{J}_{1}+6{J}_{2}+3{J}_{3}+{J}_{z}$ [50]. The estimated ferromagnetic transition temperatures ( ${T}_{c}^{m}$ ) are listed in table 2. It is clear from table 2 that the estimated ${T}_{c}^{m}$ 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₃. Nonetheless, the estimated ${T}_{c}^{m}$ values are much higher. In particular, for ML CrI₃, the estimated ${T}_{c}^{m}$ of 109 K is more than two times larger than the experimental T_c of 45 K for ML CrI₃, and in bulk CrI₃, the ${T}_{c}^{m}$ of 191 K is about three times larger than the measured T_c of 61 K [2]. The reason for this is that the MFA neglects transverse spin fluctuations [51, 52], which is especially strong in low-dimensional materials. Furthermore, the MFA violates Mermin–Wagner theorem [22] which says that long-range 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 [51, 52]. 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}_{c}^{{\rm{SWT}}}=\pi {J}_{1}/2{\rm{log}}(1+2\pi {J}_{1}/{{\rm{\Delta }}}_{0})$ where Δ₀ is the spin wave gap [17], which can be roughly set to Δ₀ = ΔE_ma. The ${T}_{c}^{{\rm{SWT}}}$ values evaluated using the calculated J₁ and ΔE_ma are listed in table 2. It can be seen that compared with T_c^m, the ${T}_{c}^{{\rm{SWT}}}$ values agree better with the corresponding experimental values although they are still too high (table 2).

3.2. Electronic structure

Better understanding of the electronic, magnetic and MO 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₃ structures are displayed in figure 3, while that of the BL and TL are displayed in figure S2 in the supplementary information (SI). Notably, these figures show that all the multilayer structures are almost direct band-gap 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₃ is thinned down to TL, to BL and finally to ML.

**Figure 3.** Relativistic band structures of (a) bulk and (b) ML CrI₃, respectively. Horizontal dashed lines denote the top of valance band.
Download figure:
Standard image High-resolution image

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 figure 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 [53] 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₃ structures. The calculated DOS spectra of bulk and monolayer CrI₃ are displayed in figures 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₃ 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 spectrum in spin-down valance band region are mostly contributed from I p orbitals. The states in the upper valance band region of −4.5 eV to −0.3 eV are made up of spin-up Cr ${d}_{{xy},{x}^{2}-{y}^{2}}$ orbitals, with a broad peak at -1.3 eV. The conduction bands in the range from 0.4 eV to 1.2 eV are mainly of spin-up Cr d_xz,yz orbitals. From these one could say that the band gap in the CrI₃ structures arise due to the crystal-field splitting of spin-up Cr ${d}_{{xy},{x}^{2}-{y}^{2}}$ and d_xz,yz bands. The spin-down 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² orbitals (figure 4). Furthermore, the calculated total as well as site, orbital projected spin-up and spin-down DOS spectra of ML CrI₃ (figure 5) are rather similar to that of bulk CrI₃ (figure 4).The main difference comes from the contributions of the I p orbital which is enhanced in the higher energy region in the ML.

**Figure 4.** Site-, orbital-, and spin-projected scalar-relativistic densities of states (DOS) of bulk CrI₃.
Download figure:
Standard image High-resolution image

**Figure 5.** Site-, orbital-, and spin-projected scalar-relativistic densities of states (DOS) of ML CrI₃.
Download figure:
Standard image High-resolution image

As mentioned above, the insulating band gap in the CrI₃ 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, 54], according to the perturbation theory. In particular, the SOC matrix elements for the d orbital of $\langle {d}_{{xz}}| {H}_{{\rm{SO}}}| {d}_{{yz}}\rangle$ and $\langle {d}_{{x}^{2}-{y}^{2}}| {H}_{{\rm{SO}}}| {d}_{{xy}}\rangle$ would favor the perpendicular magnetic anisotropy, while other elements of $\langle {d}_{{x}^{2}-{y}^{2}}| {H}_{{\rm{SO}}}| {d}_{{yz}}\rangle$ , $\langle {d}_{{xy}}| {H}_{{\rm{SO}}}| {d}_{{yz}}\rangle$ and $\langle {d}_{{z}^{2}}| {H}_{{\rm{SO}}}| {d}_{{yz}}\rangle$ would prefer an in-plane magnetic anisotropy [55]. The magnitude ratio of these d-orbital SOC matrix elements are $\langle {d}_{{xz}}| {H}_{{\rm{SO}}}| {d}_{{yz}}{\rangle }^{2}$ : $\langle {d}_{{x}^{2}-{y}^{2}}| {H}_{{\rm{SO}}}| {d}_{{xy}}{\rangle }^{2}$ : $\langle {d}_{{x}^{2}-{y}^{2}}| {H}_{{\rm{SO}}}| {d}_{{yz}}{\rangle }^{2}$ : $\langle {d}_{{xy}}| {H}_{{\rm{SO}}}| {d}_{{yz}}{\rangle }^{2}$ : $\langle {d}_{{z}^{2}}| {H}_{{\rm{SO}}}| {d}_{{yz}}{\rangle }^{2}$ =1:4:1:1:3 [55]. 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 $\langle {d}_{{xz}}| {H}_{{\rm{SO}}}| {d}_{{yz}}\rangle$ and $\langle {d}_{{x}^{2}-{y}^{2}}| {H}_{{\rm{SO}}}| {d}_{{xy}}\rangle$ and hence large contributions to the perpendicular anisotropy. In contrast, the d_z² DOS is almost zero in the lower conduction band region of 0.4–1.5eV in both bulk and ML CrI₃, and thus this would give rise to the diminishing SOC matrix element of $\langle {d}_{{z}^{2}}| {H}_{{\rm{SO}}}| {d}_{{yz}}{\rangle }^{2}$ which favors an in-plane magnetization. All these together would lead to the magneto-crystalline anisotropy energies in the CrI₃ structures (see table 1) which strongly favor the out-of-plane magnetization.

3.3. Optical conductivity

Now let us turn our attention to the calculated optical conductivity, which is the ingredient for evaluating the MO effects, as explained already in section 2. The calculated optical conductivities of the investigated structures are displayed in figure 6 for bulk and ML CrI₃ as well as in figure S4 in the SI for BL and TL CrI₃. Since the calculated optical conductivities for all the considered structures (see figures 6 and S4) are rather similar due to the weak van der Waals interlayer interaction, here we discuss only bulk and ML CrI₃ as examples. Notably, figure 6 shows that in both bulk and ML CrI₃, 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 σ_xx¹ for $E| | {ab}$ of bulk CrI₃ is significantly larger than that ( ${\sigma }_{{zz}}^{1}$ ) for $E| | c$ in the energy range from 1.5 to 4.4 eV, while it becomes smaller than ${\sigma }_{{zz}}^{1}$ for $E| | c$ above 4.4 eV. A similar profile is observed in ML CrI₃ 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 ${\sigma }_{{zz}}^{1}$ is higher than ${\sigma }_{{xx}}^{1}$ . 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(b) shows that the upper valance bands in the energy range of −4.5 to −0.3 eV are dominated by Cr ${d}_{{xy},{x}^{2}-{y}^{2}}$ orbitals. Given that ${d}_{{xy},{x}^{2}-{y}^{2}}$ (d_z²) states can be excited by only ${E}\perp c$ (E $\parallel c$ ) polarized light while Cr d_xz,yz orbitals can be excited by light of both polarizations, consequently ${\sigma }_{{xx}}^{1}$ would be greater than ${\sigma }_{{zz}}^{1}$ in the low photon energy region up to 4.3 eV (see figures 6(a) and (b)).

**Figure 6.** (a) and (b) Real, (c) and (d) imaginary diagonal components and (e) , (f) off-diagonal components of the optical conductivity tensor of bulk (left panels) and ML (right panels) CrI₃. The pink dashed–dotted lines in (a) and (b) are, respectively, the experimental differential reflection spectra (in arbitrary units) for bulk and ML CrI₃ [7] which are proportional to ${\sigma }_{{xx}}^{1}$ . The red dashed curve in (a) is the ${\sigma }_{{xx}}^{1}$ of bulk CrI₃ derived from the experimental refraction index and extinction coefficient [2]. The red diamond in (f) denotes the ${\sigma }_{{xy}}^{2}$ of ML CrI₃ derived from the experimental photoluminescence [7] (see text).
Download figure:
Standard image High-resolution image

The calculated real ( ${\sigma }_{{xy}}^{1}$ ) and imaginary ( ${\sigma }_{{xy}}^{2}$ ) parts of the off-diagonal optical conductivity element are displayed in figures 6(e) and (f) for bulk and ML CrI₃, respectively. It is quite evident from these figures that the off-diagonal element spectra of bulk and ML CrI₃ are again rather similar, due to the weak van der Waals interlayer coupling in bulk CrI₃. For example, the ${\sigma }_{{xy}}^{1}$ 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 prominent peak at ∼3.5 eV for bulk and ML structures. In the ${\sigma }_{{xy}}^{2}$ spectra, there is a large positive peak in the vicinity of 4.8 eV for both structures. There are also negative peaks in both ${\sigma }_{{xy}}^{1}$ and ${\sigma }_{{xy}}^{2}$ spectra. For example, a negative peak occurs in the vicinity of 1.3 eV in the ${\sigma }_{{xy}}^{1}$ spectrum and at 3.3 eV in the ${\sigma }_{{xy}}^{2}$ spectrum for bulk and ML structures.

The absorptive parts of the optical conductivity elements ( ${\sigma }_{{xx}}^{1}$ , ${\sigma }_{{zz}}^{1}$ , ${\sigma }_{{xy}}^{2}$ and ${\sigma }_{\pm }^{1}$ ) are directly related to the dipole allowed interband transitions, as equations (3), (4) and (10) indicate. Therefore, we further analyze the origin of the main features in the ${\sigma }_{{xx}}^{1}$ , ${\sigma }_{{zz}}^{1}$ and ${\sigma }_{{xy}}^{2}$ 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 ${\sigma }_{{xx}}^{1}$ , ${\sigma }_{{zz}}^{1}$ and ${\sigma }_{{xy}}^{2}$ for all the considered systems are rather similar, as mentioned above, here we perform the symmetry analysis for ML CrI₃ only (see the SI). The scalar-relativistic and relativistic band structures of ML CrI₃ with the symmetries of band states at Γ-point labeled, are displayed in figure S6 and figure 7, respectively. We then assign the main peaks of the ${\sigma }_{{xx}}^{1}$ and ${\sigma }_{{zz}}^{1}$ spectra (figure 6(b)) to the interband transitions at the Γ-point using the dipole selection rules (see table S3 in the SI), as presented in figure S6. For example, we could relate the A₁ peak at ∼1.5 eV in ${\sigma }_{{xx}}^{1}$ in figure 6(b) to the interband transition from the ${{\rm{\Gamma }}}_{3}^{+}$ state at −0.5 eV of the spin-down valence band to the conduction band state ${{\rm{\Gamma }}}_{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 ${{\rm{\Gamma }}}_{3}^{+}$ and ${{\rm{\Gamma }}}_{3}^{-}$ band states are doubly degenerate (figure S6), and the absorption rates for left- and right-handed polarized lights are the same, thereby resulting in zero contribution to ${\sigma }_{{xy}}$ . In the presence of the SOC, the band states split (see figure 7) and this gives rise to optical magnetic circular dichroism. For example, we could relate the peak P₁ at 1.2 eV in the ${\sigma }_{{xy}}^{2}$ of figure 6(f) to the interband transition from the ${{\rm{\Gamma }}}_{6}^{-}$ occupied states at −0.1 eV to the ${{\rm{\Gamma }}}_{4}^{+}$ unoccupied state at ∼1.1 eV (figure 7).

**Figure 7.** Relativistic band structures of ML CrI₃. The horizontal dashed lines denote the top of valance band. The states at the Γ-point in the band structure are labeled according to the irreducible representations of the C_3i double point group. The principal interband transitions and the corresponding peaks in the σ_xy in figure 6 (f) are indicated by red and green arrows.
Download figure:
Standard image High-resolution image

For comparison with the available experimental data, we plot in figures 6 (a) and (b) the experimental differential reflection spectra (in arbitrary units) of bulk and ML CrI₃, respectively [7], which are proportional to ${\sigma }_{{xx}}^{1}$ . Moreover, we also display the ${\sigma }_{{xx}}^{1}$ of bulk CrI₃ derived from the experimental refraction index and extinction coefficient [2] in figure 6(a). It can be seen from figure 6(a) that twe peaks (A₂ and A₃) in the experimental ${\sigma }_{{xx}}^{1}$ of bulk CrI₃ [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₁, A₂ and A₃) in the experimental differential reflection spectrum of bulk CrI₃ [7] also agree very well with the theoretical ones (see figure 6(a)). This indicates that our GGA + U calculations could describe the optical properties of bulk CrI₃ quite well. Figure 6(b) shows that three peaks (A₁, A₂ and A₃) in the experimental differential reflection spectrum of ML CrI₃ [7] are also quite well reproduced by our GGA + U calculations. In particular, the position and shape of peak A₁ in theoretical and experimental spectra are in good agreement, although theoretical peaks A₂ and A₃ 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 $\rho =({I}_{+}-{I}_{-})/({I}_{+}+{I}_{-})\approx ({\sigma }_{+}^{1}-{\sigma }_{-}^{1})$ / $({\sigma }_{+}^{1}+{\sigma }_{-}^{1})$ = ${\sigma }_{{xy}}^{2}/{\sigma }_{{xx}}^{1}$ , i.e. ${\sigma }_{{xy}}^{2}=\rho {\sigma }_{{xx}}^{1}$ . In figure 6(f), the red diamond denotes the ${\sigma }_{{xy}}^{2}$ estimated this way based on the calculated ${\sigma }_{{xx}}^{1}$ and experimental ρ for ML CrI₃ [7], which agree well with peak P₁ in the theoretical ${\sigma }_{{xy}}^{2}$ spectrum.

Table 1 indicates that the calculated band gap of bulk CrI₃ 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 equations (3) and (4)indicate that correct band gaps would be important for obtaining accurate optical conductivity and hence MO effects, we further perform the band structure calculations using the hybrid Heyd–Scuseria–Ernzerhof (HSE) functional [56, 57] 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 figure S7 and table S4 in the SI, respectively. Indeed, the HSE band gap of bulk CrI₃ agrees well with the experimental value [3]. Therefore, we also calculate the optical and MO properties of bulk and few-layer CrI₃ from the GGA + U band structures within the scissors correction (SC) scheme [58] using the differences between the HSE and GGA band gaps (table S4). The optical conductivity spectra calculated with the SC are displayed in figure S8 in the SI. Surprisingly, the optical conductivity spectra of bulk and ML CrI₃ calculated with the SC are in worse agreement with the corresponding experimental values [2, 7] (see figures S8(a), S8(b) and S8(j)). This could be explained as follows. We notice that bulk CrI₃ has an indirect band gap (see figure 3). Consequently, the discrepancy in the band gap between theory and experiment in bulk CrI₃ 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 figure 7).

3.4. MO 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 equation (8)). Here we assume that the substrate is SiO₂ with dielectric constant ε_s = 2.25 for all CrI₃ multilayers since SiO₂ 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 figures 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 (figure 8) are, respectively, similar to that of the real part ( ${\sigma }_{{xy}}^{1}$ ) and imaginary part ( ${\sigma }_{{xy}}^{2}$ ) of the off-diagonal conductivity element (figure 6). This could be expected because the Kerr effect and the off-diagonal conductivity element are connected through equations (7) and (8) In fact equations (7) and (8) indicate that the complex Kerr rotation angle would be proportional to the off-diagonal optical conductivity of σ_xy if the conductivity σ_xx of the sample and substrate is roughly a constant of photon energy. For the CrI₃ multilayers, the dielectric constant of ε = 2.25 of the SiO₂ substrate is a constant and thus the Kerr rotation angle is proportional to the σ_xy. For bulk CrI₃, the Kerr rotation could become large if theσ_xx, which is in the denominator of equation (7), becomes very small. This explains why the Kerr angles of the bulk are still pronounced below 0.5 eV (figure 8(a)).

**Figure 8.** Kerr rotation (θ_K) and ellipticity (ε_K) spectra of (a) bulk, (b) ML, (c) BL, and (d) TL CrI₃. The filled red circle in each panel indicates the θ_K value from the recent experiment [2].
Download figure:
Standard image High-resolution image

**Figure 9.** Faraday rotation (θ_F) amd ellipticity (ε_F) spectra of (a) bulk, (b) ML, (c) BL, and (d) TL CrI₃. The filled red circles in (a) denote the measured θ_F values of bulk CrI₃ [4].
Download figure:
Standard image High-resolution image

The overall behavior of the MO spectra for all the investigated structures (figures 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₃ 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.2 eV for the TL (−1.3°) as well as at 3.3 eV for the ML (−0.5°) and 4.2 eV for the BL (−1.0°). Similarly, Kerr ellipticity shows the negative maximum near 4.8 eV for the bulk (−1.9°) as well as at 4.9 eV for the ML (−0.6°) and at 4.8 eV for BL (−1.2°) and TL (−1.3°). These large Kerr angles and ellipticities indicate that the considered structures would have valuable applications in MO nanodevices.

To access the application potential of these CrI₃ 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, 59]. Pt alloys also possess large Kerr rotations in the range of 0.4°–0.5° in FePt and Co₂Pt [43] and PtMnSb [60]. In these materials, the strong SOC of heavy Pt atoms was shown to play an important role [61]. 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, 62] and also of 0.06° at 1.96 eV for Co/Au multilayer [63]. Among semiconductor MO materials, diluted magnetic semiconductors Ga_1−xMn_xAs were reported to show significant Kerr rotation angles of ∼0.4° near 1.80 eV [64]. Y₃Fe₅O₁₂ also exhibits a significant Kerr rotation of 0.23° at 2.95 eV [65]. In strong contrast, the measured Kerr rotation angles of bulk and atomically thin films of Cr₂Ge₂Te₆ 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₃ here are comparable to these important MO metals and semiconductors. Thus, because of their excellent MO properties, CrI₃ 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 figure 9. The main features of the Faraday rotation spectra of the considered structures resemble that of the Kerr rotation spectra (see figure 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, 59]. Ferromagnetic semiconductor Y₃Fe₅O₁₂ show only small Faraday rotation angles of ∼0.19°/μm at 2.07 eV [66]. However, the replacement of Y with Bi considerably enhances the Faraday rotations to ∼35.0°/μm at 2.76 eV [67], thus making Bi₃Fe₅O₁₂ an outstanding MO material. Clearly, the calculated Faraday angles of the CrI₃ 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₃ have recently been measured at photon energy of 1.96 eV [2] and they are found to be large (see figure 8). If the calculated Kerr rotation spectra are red-shifted by merely 0.2eV, the measured and calculated θ_K values of bulk and ML CrI₃ would be in very good agreement, although the experimental θ_K values of BL and TL CrI₃ are much larger than the corresponding theoretical θ_K values (figure 8). Furthermore, the measured Faraday rotation angles θ_F of bulk CrI₃ [4] are almost in perfect agreement with our theoretical prediction (see figure 9(a)). We notice that in contrast, the photon energies of the experimental Faraday angles of bulk CrI₃ [4] differ significantly from that of the theoretical Faraday angles calculated with the SC using the HSE band gap (see figure S10(a) in the SI).

4. Conclusion

Summarizing, we have systematically studied magnetism, electronic structure, optical and MO properties of bulk and atomically thin films of CrI₃ by performing extensive GGA + U calculations. Firstly, we find that both intralayer and interlayer magnetic couplings are strongly ferromagnetic in all the considered CrI₃ 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₂Ge₂Te₆ [26]. These large MAEs suppress transverse spin fluctuations in atomically thin CrI₃ 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₃ structures, MO Kerr rotations up to ∼2.6° are found, being comparable to famous MO materials such as ferromagnetic semiconductor Y₃Fe₅O₁₂ 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₃Fe₅O₁₂ and Bi₃Fe₅O₁₂. Calculated optical conductivity spectra, MO Kerr and Faraday rotation angles agree quite well with available experimental data. In BL CrI₃, the interlayer magnetic coupling is found to be layer stacking dependent while the optical and MO properties are not. Therefore, further experimental determination of the structure of BL CrI₃ 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 MO conductivity spectra of these structures are analyzed in terms of their underlying electronic band structures. Our findings of large out-of-plane MAEs and strong MO effects in these single-spin ferromagnetic semiconducting CrI₃ ultrathin films suggest that they will find promising applications in high density semiconductor MO and low-power spintronic nanodevices as well as efficient spin-photovoltaics.

Acknowledgments

VKG and GYG acknowledges the support by the Ministry of Science and Technology and the National Center for Theoretical Sciences, Taiwan. GYG also thanks the Academia Sinica Thematic Research Program (AS-TP-106-M07) and the Kenda Foundation of Taiwan for the supports.

Magnetism and magneto-optical effects in bulk and few-layer CrI₃: a theoretical GGA + U study

Article metrics

Author e-mails

Author affiliations

ORCID iDs

Dates

Abstract

1. Introduction

2. Computational methods