Distinctive magnetic properties of CrI₃ and CrBr₃ monolayers
caused by spin-orbit coupling

We start from the structural properties of monolayer CrX₃, where X stands for ligand I and Br atoms. CrX₃ crystallizes in the trigonal $P\overline{3}1m$ space group. The hexagonal planar lattice of Cr atoms is sandwiched between triangular planar lattices of ligand atoms as shown in Figs. 1 (a) and (d) for respective materials. One Cr atom bonds to 6 ligand atoms and each ligand atom bonds to 2 Cr atoms. The structures have 3-fold in-plane symmetry, where the triangular lattices of ligand atoms are 180^∘ rotated with respect to each other, which creates octahedral coordination around each Cr atom. As listed in Table 1, the Cr-I and Cr-Br bond lengths are found to be 2.78 Å and 2.55 Å, respectively, which is directly proportional with the radius of the bonding orbital of the ligand, I-5$p$ and Br-4$p$. Consequently, the lattice constants of monolayer CrI₃ and CrBr₃ are also different, 6.94 Å and 6.40 Å, respectively, consistent with the experimentally measured 6.95 Åli2020single and 6.50 Å.chen2019direct Note that both monolayers exhibit slightly larger lattice constant than their bulk or few-layer counterparts. The covalent bonding character consists of 1.10 $e^{-}$ and 1.33 $e^{-}$ donation of Cr and 0.37 $e^{-}$ and 0.44 $e^{-}$ gain of I and Br atoms, respectively. This bonding charge is visualized in Figs. 1 (b) and (e), as obtained by subtracting bare atom charge distributions from the charge distribution of the crystals. Charges obviously accumulate to the bonding sites and around Cr atoms. Observed charge accumulation around each ligand clearly shows three ’charge clouds’ that can be divided in two groups, representing two types of bonding between Cr and ligand atoms: the facial ($\sigma$) bonding, represented by two charge clouds facing the Cr atoms, and the lateral ($\pi$) bonding, represented by the third charge cloud, in direction orthogonal to the Cr-X-Cr-X plane.

To depict the energetic arrangements of the orbital states, we calculated the partial density of states for two considered materials in Figs. 1 (c) and (f). Both materials are semiconductors, in which the valence band maximum is dominated by $p$-orbital of the ligand. On the other hand, $d$-orbital of Cr resides much deeper in the valence band, with similar energetic delocalization for both CrI₃ and CrBr₃. $d_{x^{\prime}y^{\prime}}$, $d_{x^{\prime}z^{\prime}}$ and $d_{y^{\prime}z^{\prime}}$ are degenerate, and are plotted together. $d_{z^{\prime 2}}$ and $d_{x^{\prime 2}-y^{\prime 2}}$ orbitals are also degenerate and appear at the conduction band. $p_{x^{\prime}}$ and $p_{y^{\prime}}$ orbitals are degenerate as well, and dominate the valence band maximum. $p_{z^{\prime}}$ orbital mostly resides in the middle of the valence band. One should note that in order to obtain orbital states compatible with the octahedral coordination, the general coordinates are rotated as suggested by Rassekh et al.rassekh2020remarkably . Therefore we consider $xy$-plane and $z$-axis of the general Cartesian coordinates aligned with the local coordinates, where the Cr-X-Cr-X plane is considered as $x^{\prime}y^{\prime}$-plane, and $z^{\prime}$ axis is orthogonal to that plane. Orbital decomposition shows that both materials exhibit very similar orbital delocalizations. This is also visible in the charge density variations shown in Figs. 1 (b) and (e). Purple regions show the depletion of charges of isolated atoms through the charge density of the crystal. It is obvious that the orbitals laying along the Cr-X bonds exhibit most delocalization. We understand from the charge depletion and the orbital decomposition of DOS that $d_{x^{\prime}y^{\prime}}$, $d_{x^{\prime}z^{\prime}}$ and $d_{y^{\prime}z^{\prime}}$ orbitals are localized and do not show any variation from their single-atom form. $d_{z^{\prime 2}}$ and $d_{x^{\prime 2}-y^{\prime 2}}$ orbitals form a $spd^{2}$ hybridization with the $p_{x^{\prime}}$ and $p_{y^{\prime}}$ orbitals as expected for an octahedral coordination. Beside these similarities, the particular difference between CrI₃ and CrBr₃ is found between 4$p$ of Br and 5$p$ of I, leading to energetic differences of the $spd^{2}$ hybridization. Since 4$p$ orbital of Br is more confined as compared to 5$p$ of I, the states of CrBr₃ are shifted to higher energies. Therefore, the band gap values are different, found as 0.55 eV and 1.59 eV for CrI₃ and CrBr₃, respectively.

The responses of monolayer CrI₃ and CrBr₃ to biaxial strain are further revealing the influence of SOC on the magnetic exchange interaction. Principally, the structural changes upon straining are similar for CrI₃ and CrBr₃. Cr-Cr distance changes according to the degree of strain applied ($\pm 5\%$), however, the Cr-X bond length varies $\sim\pm 1\%$. Most of the changes therefore accumulate to the Cr-X-Cr (X = I, Br) angle, varying $\sim\pm 5\%$. Such changes in the angle cause distortion in the octahedral coordination, as illustrated in Fig. 4(a). Therefore, it is expected that strain induces the interaction of two $spd^{2}$ bonds of one ligand due to bong-angle modification and leads to modification of the SOC contribution in both Cr and the ligand.

We next calculated the total energy difference and the total SOC energy difference between out-of plane and in-plane spin alignments under biaxial strain. As shown in Figs. 4 (b) and (c), $E_{\text{SO}}$ of monolayer CrI₃ [orange curve in Fig. 4(b)] linearly increases under increasing tensile strain, which is similar to its behavior in CrBr₃ [green curve in Fig. 4(c)]. However, the behavior of the two materials under compressive strain is completely different. CrI₃ exhibits parabolic increase with increasing compressive strain, while CrBr₃ maintains linear behavior and decreases with increasing compressive strain! The atomic contributions reveal the source of these behaviors. As shown by magenta curve in Fig. 4(b), iodine contribution dominates the behavior for both compressive and tensile strain. For CrBr₃ in Fig. 4(c), in case of the compressive strain, the contributions from Cr and Br are almost equal but have opposite sign, meaning that Br atoms favor to have spin in out-of-plane direction while Cr atoms prefer in-plane spin direction. For tensile strain, Br contribution dominates and both Cr and Br prefer out-of-plane spin alignment. Since the structural changes of both monolayers under strain are similar, one expects similar variation of Cr energies, however, Cr of CrI₃ exhibits three times larger energy variation as compared to Cr of CrBr₃. This indicates that the anisotropy related with Cr also originates from the bonding electrons rather than Cr-only electrons. In Figs. 4 (d) and (e) we show the charge density variation for $5\%$ and $-5\%$ strain in both CrI₃ and CrBr₃. It is clearly seen that in case of compressive $-5\%$ strain the bonding charges of two different bonds of one ligand overlap. This reveals the origin of the interaction under compressive strain. For $5\%$ tensile strain, the bonding charges are clearly separated and exhibit no overlap.

Our results for behavior of MAE are generally in good agreement with results reported by Ref. PhysRevB.98.144411, . In case of CrBr₃ we also obtain that MAE is smallest in case of $-5\%$ strain and it is linearly growing, reaching the maximum for $+5\%$ strain. In case of CrI₃, however, our results agree with Ref. PhysRevB.98.144411, only for the compressive strain - with increased compressive strain, anisotropy is growing. On the other hand, with tensile strain, we report completely different behavior of MAE in CrI₃, compared to Ref. PhysRevB.98.144411, . There, anisotropy was decreasing with increased tensile strain, while in our study, the behavior is different - anisotropy grows with increased tensile strain.

The variations of exchange parameters and SIA under biaxial strain are also examined based on the considerations related to Table 2 and Fig. 3. The mean values $\langle J\rangle$ of the exchange matrices of three NNs are plotted as a function of strain in Fig. 5(a). The ferromagnetic exchange interaction increases with the tensile strain and decreases with compressive strain for both CrI₃ and CrBr₃, for all components. The variation on the curves of CrI₃ is much larger such that the compressive strain almost equalizes CrI₃ and CrBr₃ in terms of exchange interaction, and the difference between two monolayers increases with tensile strain. These results agree with Ref. PhysRevB.98.144411, only in case of compressive strain. Namely, both our study and mentioned previous work suggest that with compressive strain, two materials become less FM. However in case of tensile strain, our results suggest that materials are more FM than in pristine case, while in Ref. PhysRevB.98.144411, , the opposite was suggested. In Fig. 5(b), the anisotropy between in-plane and out-of-plane exchange parameters, $\Delta$, is plotted as a function of strain. $\Delta$ of CrI₃ increases for increasing either tensile or compressive strain. It is important to note that the behavior of $\Delta$ is consistent with the behavior of MAE in Fig. 4, since most of the contribution to MAE comes from the anisotropy of the exchange interaction of the first NN. For CrBr₃, the anisotropy slightly increases (decreases) under tensile (compressive) strain, which is also consistent with the behavior of MAE. In Fig. 5(c), SIA is presented as a function of strain. SIA in CrI₃ exhibits a gradually slowing decrease (in absolute value) when moving from compressive to tensile strain. In CrBr₃ on the other hand, SIA exhibits opposite behavior to that of CrI₃. It is significant that for compressive strain beyond $-3\%$ SIA parameter becomes positive, indicating preference for in-plane direction of magnetization. Notice that SIA is an order of magnitude smaller than the exchange parameter, therefore its contribution to overall magnetic properties is limited. However, MAE of CrBr₃ under high compressive strain is of the same order as SIA, indicating that the drop of MAE is due to the decrease of SIA parameter, while $\Delta$ stays almost constant.

The above-indicated differences in the strength of the magnetic exchange parameters, exchange anisotropy, and SIA parameter between two magnetic monolayers under investigation can be monitored via temperature-dependent magnetization and hysteretic behavior in applied magnetic field, which are both readily experimentally accessible. Having obtained all the parameters to construct the Heisenberg spin Hamiltonian in Eq. (1) for strained monolayers, we next calculate the temperature-dependent magnetization of CrI₃ and CrBr₃ and the corresponding $T_{\text{C}}$ where the magnetic phase transits from paramagnetic to ferromagnetic state. We used stochastic LLG simulation to obtain the temperature-dependent magnetization M_z/M_s($T$), where M_z is the out-of-plane magnetization and M_s is the saturation magnetization. As shown in Figs. 6 (a) and (b), the obtained $T_{\text{C}}$ of the unstrained monolayers of CrI₃ and CrBr₃ of 56 K and 38 K is reasonably close to the experimentally obtained values of 45 K and 21 K, respectively. These $T_{\text{C}}$ values are also consistent with those found in previous workstorelli2018calculating ; torelli2019high , obtained using renormalization spin-wave theory combined with classical MC calculations¹¹1Note that the feeding parameters in those works were obtained by DFT calculations with slightly different parametrization than used in the present work, such as $U=3.5$ eV instead of $U=4$ eV.. Further, as shown in Fig. 6(c), $T_{\text{C}}$ decreases (increases) under compressive (tensile) strain, mainly due to the previously described strong variation of the exchange parameters with strain. In case of CrBr₃, contrary to CrI₃, the influence of $\Delta$ and SIA is very small since the variation of those parameters under strain can be considered negligible. One should note that the behavior of $T_{\text{C}}$ as a function of strain presented here is completely different than that reported in Ref. PhysRevB.98.144411, where $T_{\text{C}}$ is calculated using mean-field theory. Such disagreement is hardly a surprise since in Ref. PhysRevB.98.144411, a single exchange parameter is used to determine $T_{\text{C}}$ in absence of information on anisotropy. Further we note that two monolayer materials have almost equal $T_{\text{C}}$ at compressive $-5\%$ strain. In the case of tensile strain on the other hand, the difference between $T_{\text{C}}$ of CrI₃ and CrBr₃ increases with strain. CrBr₃ reaches the saturation of $T_{\text{C}}\approx 45$ K at around 3% strain while CrI₃ exhibits $T_{\text{C}}\approx 74$ K at 5% strain and tends to further increased $T_{\text{C}}$ with further straining.

The behavior of M_z/M_s($T$) curves is also illustrative of the difference between the two materials. In general, all curves are well fitted by the functional behavior $(1-T/T_{\text{C}})^{\beta}$, where $\beta$ is the critical exponent. In Fig. 6(d), obtained exponents $\beta$ are plotted as a function of strain, together with $\beta$ values from different available models for comparison. In our results, the critical exponent of monolayer CrI₃ at all strains can be approximated by $\beta\approx 0.24$. This value suggests that strong out-of-plane anisotropy of CrI₃ separates its behavior from those expected by 3D models and sets it closer to the 2D limit. Our value of $\beta$ is also comparable with the value of $\sim 0.26$ obtained in Ref. PhysRevB.97.014420, for bulk CrI₃ since the intralayer interactions dominate the magnetic behavior even in bulk CrI₃.²²2It is interesting to mention that the critical exponent of the 2D XY model is estimated as $\beta\approx 0.23$bramwell1993magnetization from the experimental observations on several layered magnets such as BaNi₂(PO₄)₂, Rb₂CrCI₄ and K₂CuF₄, which is very close to $\beta$ of CrI₃. We believe this to be a coincidence, since 2D XY model is based on the easy-plane anisotropy contrary to the out-of-plane easy-axis of CrI₃. It is also known that 2D XY model does not exhibit standard spontaneous phase transition, which should deviate from the behavior we obtained in our simulations at around $T_{\text{C}}$. CrBr₃ with $\beta\approx 0.3$, on the other hand, is much closer to the 3D Ising value ($\beta=0.3226$) due to weak out-of-plane anisotropy. Our result is validated by (albeit at the lower margin of) $\beta=0.4\pm 0.1$ recently obtained from magnetization measurements by Kim et al.kim2019micromagnetometry .

Finally, motivated by its accessibility by advanced magnetometry,kim2019micromagnetometry we calculated the hysteretic behavior of the monolayers under magnetic field $B$ tilted by angle $\theta_{B}$ from the $z$-axis, for both CrI₃ and CrBr₃.

In the simulations, the spin system is initialized at high applied field, where all the spins are aligned to the magnetic field direction. The field is then looped in steps of $0.01$ T and $0.05$ T for CrBr₃ and CrI₃ respectively, where the magnetization is relaxed by $5\times 10^{4}$ time steps for each value of field. Dipole-dipole interactions have been included in the simulations. In Fig. 7(a,b), we show the obtained $\Delta\text{M}_{z}=|\text{M}_{z}(\rightarrow)-\text{M}_{z}(\leftarrow)|$ as a color-map plot where shades of red (CrI₃) and blue (CrBr₃) color show the deviation of $M_{z}(B)$ hysteresis loop from the rectangular shape. In Figs. 7 (c) and (d), we plot the corresponding hysteresis curves for selected tilt angles of the applied magnetic field. One clearly sees that hysteresis loops evolve from sharp rectangular to an oval form as the tilt angle is increased, as was recently shown experimentally by Kim et al. for the case of monolayer CrBr₃.

In absence of dipolar interactions, the behavior of magnetic spins in external field is captured by minimization of the energy $E_{TOT}(\theta)=-\rm{MAE}\cos(2\theta)-E_{B}\cos(\theta-\theta_{B})$, where $E_{B}$ is the energy associated with the magnetic field, and grows linearly with increasing $B$. The value of the critical, ’switching’ field for the given angle $\theta_{B}$, scaled to the switching field for $\theta_{B}=0^{\circ}$ [$\Gamma=B_{cr}(\theta_{B})/B_{cr}(0^{\circ})$] is then analytically obtained from the condition that the energy extrema coalesce, leading to equation

This functional dependence of critical field on the tilt angle $\theta_{B}$ perfectly reproduces the numerically calculated switching field, shown by dashed lines in Figs. 7(a,b). Notably, the switching field in absence of dipolar interactions ($\Gamma(\theta_{B})$) shows symmetric behavior with respect to $\theta_{B}=45^{\circ}$. However, with dipolar interactions included, that symmetry is broken in case of CrBr₃, and switching field for $\theta_{B}=0^{\circ}$ becomes significantly lower than for $\theta_{B}\rightarrow 90^{\circ}$. The latter was indeed validated experimentally, in Ref. kim2019micromagnetometry, . However, dipolar interactions cause no changes in the critical field of CrI₃ for any $\theta_{B}$, which is another important distinction between two materials that could be verified by Hall micromagnetometry.

One should however note that in our considerations we do not involve finite size effects nor demagnetization, or temperature fluctuations, likely playing an important role in experiment (next to the ever-present defects in the monolayers, that can facilitate magnetic reversal locally). Our simulations explore primarily the effect of the microscopic parameters on the apparent magnetic behavior, with a goal of capturing the intrinsic differences between CrI₃ and CrBr₃. As a consequence, the switching fields in our simulations are significantly larger than experimentally reported values of $\approx$ 0.15T huang2017layer for CrI₃ and $\approx$ 0.03Tkim2019micromagnetometry for CrBr₃, for $\theta_{B}=0$. Having said that, our simulations capture the large ratio between switching fields of the two monolayer materials (approximately 8 and 5 for simulation and experiment, respectively).

To feature yet another experimentally verifiable difference between monolayer CrI₃ and CrBr₃, we also calculated the hysteresis curves for the strained magnetic monolayers under out-of-plane magnetic field. We recall that MAE and SIA, plotted in Fig. 4(b,c) and Fig. 5(c), respectively, showed distinctively different behavior in two materials when under strain. As shown in Figs. 7 (e) and (f), the width of the hysteresis loop strongly changes with the strain. However, while the switching magnetic field of monolayer CrI₃ increases under both compressive and tensile strain, the switching field of monolayer CrBr₃ is reduced by compressive strain while the tensile strain increases the switching field stronger than was the case in CrI₃ (measured with respect to the unstrained case). This behavior can thus be mapped on the behavior of exchange anisotropy and SIA under strain, all rooted in the known difference in spin-orbit coupling between the two materials.