A combined first principles study of the structural, magnetic, and phonon properties of monolayer CrI₃

2D materials represent an exciting new frontier for materials research.¹ Due to their reduced dimensionality, these materials tend to exhibit stronger and longer-range electron correlation than their 3D counterparts that give rise to exotic new physics and phase behavior, including Moiré patterns,^2,3 two-dimensional superconductivity,⁴ and exotic spin and charge density waves.^5–7 The properties of 2D materials can also be tuned just by stacking,^8–10 crinkling,¹¹ straining,^12–14 and twisting^15–18 them. This versatility facilitates the layer-by-layer construction of designer materials with unique properties through the careful selection and ordering of their constituent layers.¹⁹ 

An exciting recent development in this regard is the discovery of new magnetic 2D materials.²⁰ While the Mermin–Wagner theorem^21,22 prohibits finite-temperature magnetism for the isotropic Heisenberg model in 2D, magnetic 2D materials such as CrI₃,²³ Fe₃GeTe₂,²⁴ and VSe₂²⁵ have large anisotropies arising from strong spin–orbit couplings (SOCs) or structural anisotropies that break continuous symmetries and allow for their finite-temperature long-range magnetic order. 2D magnetic materials are of great interest as they may accelerate the development of next-generation spintronic data transmission²⁶ and storage²⁷ technologies. Furthermore, if integrated into 2D heterostructures, magnetic 2D materials may additionally be able to induce magnetism in nearby nonmagnetic layers through proximity effects, enabling the realization of magnetic graphene.²⁸ Few-layer magnetic materials that have particularly strong magneto-optical responses and/or electron–phonon coupling, e.g., MPSe₃/TMD heterointerfaces,²⁹ with M = Mn, Fe, WSe₂/CrI₃ heterostructures,³⁰ or bilayer (BL) CrI₃,³¹ are also potentially promising for applications that bridge spintronics, photonics, and phononics.

However, in order to realize the potentials of 2D materials, their structural, electronic, magnetic, and optical properties have to be well characterized and understood. Experimentally, this is very challenging for mono- and few-layer systems. Monolayers (MLs) are often grown on substrates that can distort their intrinsic geometries.^32–34 Due to their inherently small thickness, monolayers are also not immediately accessible to neutron and other scattering experiments, which typically require a critical thickness to yield meaningful scattering patterns.³⁵ Given these experimental challenges, it is imperative to predictively and accurately model these materials’ atomic geometries, spin densities, and magnetic moments using first principles methods in order to elucidate the origins of their magnetic anisotropy and also to inform models upon which more macroscopic predictions can be built.^36,37 Recent advances in optical response-based methods, such as fluorescence and magneto-optical Kerr effect (MOKE) microscopy, can resolve magnetic moments down to sub-micrometer scale,³⁷ and single-spin microscopy can improve this resolution down to nanometer scale,³⁸ complementing such first principles-based predictions.

To date, the overwhelming majority of first principles modeling of these materials has relied on Density Functional Theory (DFT).³⁹ Even though DFT has led to transformative changes in our understanding of materials and chemistry and typically yields good results for ground state properties in 3D, especially lattice constants, it can yield widely varying results for the properties of 2D materials, such as lattice parameters and bandgaps, depending upon the functional employed. Recent studies of 2D materials have shown that DFT consistently overestimates interlayer binding energies between 2D monolayers⁴⁰ and predicts lattice parameters that routinely differ from experimental ones by several percent.^41,42 In contrast, Diffusion Monte Carlo (DMC), a real-space, many-body quantum Monte Carlo (QMC) method, routinely predicts lattice constants to within 1% and bandgaps to within 10% of their experimental values for the same materials.^41–43 It is worth emphasizing that DMC obtains such high-accuracy results for all properties within a single, consistent framework with few approximations; moreover, the approximations can be systematically improved upon because of the variational nature of DMC. In contrast, obtaining accurate DFT results of different properties often necessitates using a different functional for each and every property, which substantially reduces the predictive power of such calculations.

In this work, we use a combination of DFT and DMC⁴⁴ methods to construct a full atomistic picture of the physics and properties of the CrI₃ monolayer (ML). ML CrI₃ is a strongly correlated, magnetically tunable Mott insulator with a hexagonal lattice structure in which each Cr atom is coordinated by six I atoms to form a distorted octahedron (see Fig. 1).⁴⁵ Atomic force microscopy and scanning tunneling microscopy point to a lattice constant of roughly 7 ų² and a monolayer thickness of 0.7 nm.^23,32 What makes ML CrI₃ particularly intriguing as a material is its magnetism. In its monolayer form, CrI₃ has been demonstrated by MOKE microscopy²³ and single-spin nitrogen vacancy microscopy³⁸ experiments to be a ferromagnet. The strong magneto-crystalline anisotropy with an out-of-plane easy axis makes the Ising model a natural starting point for modeling the magnetic properties of CrI₃ as it allows for finite-temperature long-range order in 2D.²¹ However, several DFT studies have suggested that the material may be better described by the Heisenberg and Kitaev models.^46,47 These studies concluded that the ligand superexchange within slightly distorted Cr–I crystal environments⁴⁶ and spin–orbit coupling⁴⁷ ultimately stabilize ML CrI₃’s ferromagnetism. However, one DFT study of bilayer (BL) CrI₃ found non-negligible magnetic moments present on the iodine atoms,⁴⁸ suggesting that a more comprehensive magnetic model may be needed to quantify all of the relevant exchange interactions. In addition to its magnetism, CrI₃ exhibits strong electron–phonon coupling and magnetoelasticity, as evidenced by its alternating ferromagnetic (FM) odd-layered structures/antiferromagnetic (AFM) even-layered structures⁹ and magnetic ordering-dependent Raman spectra.⁴⁹ Although computational studies of ML CrI₃’s properties have previously been performed using a variety of DFT functionals, including the generalized gradient approximation (GGA),^45,50,51 hybrid functionals such as Heyd–Scuseria–Ernzerhof (HSE06),⁵¹ and the Local Density Approximation (LDA) with added spin–orbit coupling (SOC) corrections (LDA_SOC),⁵² they have yielded conflicting information regarding magnetic moments, lattice constants, and exchange parameters because of the strong functional dependence of these properties (see Table SI of the supplementary material). These discrepancies underscore the need to employ many-body approaches beyond DFT with fewer or no uncontrolled approximations in order to faithfully resolve ML CrI₃’s remaining controversies and uncertainties.

We utilize a combination of DFT with an on-site Hubbard-U correction (DFT + U) as well as many-body DMC simulations to resolve the structural, magnetic, and phonon properties of monolayer CrI₃. Our results show that previous DFT predictions of ML CrI₃’s magnetic spin moments⁴⁸ are correct on average across a unit cell but miss critical local spatial fluctuations and anisotropies in the spin density. Our DMC calculations predict a magnetic moment of 3.62 μ_B per chromium atom and −0.145 μ_B per iodine atom, which is considerably larger than past estimates in the literature. We furthermore exploit a recently developed surrogate Hessian DMC line search technique⁴² to determine CrI₃’s monolayer geometry with DMC accuracy, yielding high-accuracy bond lengths and bond angles that resolve previous structural ambiguities. Finally, our calculations also reveal that ML CrI₃ possesses a substantial spin–phonon coupling, 3.32 cm⁻¹, in line with couplings recently observed in other magnetic 2D materials. These findings indicate a far more anisotropic spin density involving more substantial ligand magnetism than previously thought.

In Sec. II, we describe the DFT and quantum Monte Carlo methods used, and in Sec. III, we present results on the monolayer’s structure, spin density, magnetic moments, and phonon spectra. Finally, Sec. IV contains a summary and conclusions.

The ground state properties of the CrI₃ monolayer, including its lattice constant and magnetic moments, were modeled using DFT + U, Variational Monte Carlo (VMC), and DMC. DFT simulations were first performed to obtain reference structural data, guide our surrogate Hessian line search, and generate DFT + U wave functions that were subsequently used as inputs into progressively more accurate VMC and DMC simulations. VMC and DMC calculations were then performed to obtain CrI₃’s geometry, spin density, and magnetic moments. In addition, DFT + U simulations with a self-consistent Hubbard-U correction⁵³ were employed to investigate the effects of long-range magnetic ordering on lattice phonons where the range of U values considered was guided by our DMC findings.

DFT calculations were performed using the Vienna Ab Initio Simulation Package (VASP).^55–57 The ion–electron interaction was described with the projector augmented wave (PAW) method.^58,59 A cutoff energy of 500 eV was used for the plane-wave basis set. The Brillouin zone was sampled with an 8 × 8 × 1 Monkhorst–Pack k-point mesh for the unit cell of CrI₃, which includes two Cr and six I atoms. A wide range of exchange correlation (XC) functionals were used to investigate the electronic and magnetic properties of the CrI₃ monolayer. The Local Density Approximation (LDA), Perdew–Burke–Ernzerhof (PBE), PBE + U, and PBEsol⁶⁰ functionals were specifically employed in this study. All calculations were spin-polarized. In its ferromagnetic (FM) configuration, all of CrI₃’s magnetic moments were initialized in the same out-of-plane direction, while in its antiferromagnetic (AFM) configuration, the two Cr atoms per unit cell were set to have antiparallel spins. A large vacuum of more than 25 Å along the z-direction was employed to avoid artificial interactions between images. The energies were converged with a 1 × 10⁻⁸ eV tolerance. In addition to screening various XC functionals, a self-consistent Hubbard U (U_LR) was determined from first principles by using the linear response approach proposed by Cococcioni and De Gironcoli,⁵³ in which U was determined by the difference between the screened and the bare second derivative of the energy with respect to localized state occupations. The linear-response calculation was also performed using Quantum ESPRESSO, using the same high-quality (but hard) pseudopotentials as were used to generate trial wave functions for the Quantum Monte Carlo simulations (discussed below). We obtained a linear-response U (U_LR) value of ∼3.3 eV for Cr with the LDA functional.

VMC and DMC calculations were undertaken using QMCPACK,⁴⁴ a process that was facilitated by a Nexus workflow.⁶¹ All QMC calculations were performed using a norm-conserving scalar-relativistic Opium-generated Cr pseudopotential and a norm-conserving non-relativistic iodine pseudopotential that was developed by Burkatzki, Fillipi and Dolg.⁶² The comparison of DMC energies obtained using nonrelativistic and relativistic iodine pseudopotentials yielded an energy difference of 0.009 meV/f.u., which was deemed negligible for this study. Both pseudopotentials were of the Troullier–Martins flavor. Trial wave functions used as starting points for the VMC calculations were produced using LDA + U with U = 2 eV within Quantum ESPRESSO. In the DFT calculations, a plane-wave energy cutoff of 300 Ry was necessary to converge the total energy due to the high-quality (but hard) pseudopotentials. The trial Jastrow factor consisted of inhomogeneous one- and homogeneous two-body terms, each represented by sums of 1D B-spline-based correlation functions in electron–ion or electron–electron pair distances. The trial Jastrow factor was optimized using the linear method. The VMC energy and variance were converged at a B-spline meshfactor of 1.2 using the fixed optimized Jastrow factor.

To determine a time step that is a reasonable compromise between speed and accuracy, as well as the U in LDA + U trial wave functions that minimizes fixed-node errors, preliminary calculations were performed on a 2 × 2 × 1 supercell. These calculations demonstrated that a time step of 0.01 Ha⁻¹ is capable of converging the DMC energies to within 0.01 Ha per formula unit (see Fig. 2). Moreover, as depicted in Fig. 3, the lowest DMC energies were achieved using a LDA + U trial wave function with U = 2 eV. As shown in Fig. SI, U = 2 eV also minimized the DMC energy when a smaller time step of 0.001 Ha⁻¹ was employed, demonstrating that this U is robust to time step errors when sufficiently small time steps are used. A time step of 0.01 Ha⁻¹ and LDA + U trial wave functions with U = 2 eV were thus used in all of our subsequent DMC production runs. Note that the obtained optimal value of U from DMC is close to that of the U_LR ∼ 3.3 eV value obtained from linear-response calculations. Furthermore, both values agree with statistically indistinguishable total energies from a recent DMC study of bulk CrI₃ (see Fig. 1 in Ref. 63). All of our calculations also utilized a 3 × 3 × 1 QMC twist grid and a meshfactor of 1.2.

VMC calculations employing twist-averaged boundary conditions were performed for increasingly larger twist grids and were shown to converge with a 3 × 3 × 1 twist grid; this size twist grid was thus used for twist-averaging of all subsequent VMC and DMC calculations to minimize one-body finite size effects.⁶⁴ We further employed extrapolation to the thermodynamic limit of Model Periodic Coulomb-corrected and uncorrected DMC energies for increasingly larger supercells with 1 × 2 × 1 (16 atoms), 2 × 2 × 1 (32 atoms), 2 × 3 × 1 (48 atoms), and 3 × 3 × 1 (72 atoms) tilings to minimize the two-body finite-size errors (Fig. 4).^44,65 All DMC runs used the T-moves scheme for pseudopotential evaluation to minimize localization errors.^66,67

To predict the site-averaged atomic magnetic moment per chromium, M_Cr, and iodine, M_I, the spin densities, ρ_s(r), obtained from both our DFT and DMC simulations were integrated up to a cutoff radius r_cut, defined as the zero-recrossing radius of the sign of the spin density (see Fig. SVI). In particular, we sum over the spherically interpolated spin densities within a 24-atom unit cell to predict the magnetic moment per atom (M_A),

where the r_i denote the distances from the center of the atom to the given point on the grid. To facilitate a direct comparison between the DFT and DMC spin densities, the spin densities obtained from DFT simulations were interpolated onto the dimensions of the DMC spin density grid, before being mapped onto a spherical grid.

In order to predict highly accurate lattice parameters for ML CrI₃, we used a surrogate Hessian-accelerated DMC energy-based structural optimization method which has recently demonstrated robust performance when applied to two-dimensional GeSe.⁴² This method leverages the DFT energy Hessian, which is substantially cheaper to compute than the DMC energy Hessian, to direct where to compute DMC energies based on optimal statistical sampling to ultimately determine the lowest-energy material geometry. This line search method resolves structural parameters to an accuracy higher than that obtainable via DFT energy gradient-based structural optimization techniques; this resolution is particularly important for materials which exhibit sensitive coupling between electronic, magnetic, and structural degrees of freedom, such as 2D CrI₃.

Our line search begins by first constructing an approximate energy Hessian (H_p), or a force constant matrix, using a cheaper theory, which in our case is DFT. This step is meant to (1) address the expense that would result from exploring an arbitrarily complex PES from scratch using only QMC energies and (2) minimize the noise that would result from numerous stochastic energy evaluations by allowing for a few optimally placed QMC energy calculations. We expanded the DFT PES to the second order in the Wyckoff parameter space, p, as in Ref. 68,

and diagonalized the Hessian to obtain search directions conjugate to the PES isosurfaces,

where the columns of U^T form an optimal basis of parameter directions for the line search.

With these search directions, a set of structures containing an equilibrium structure and a total of 18 “strained” structures are generated, which comprise the structure population for the first line search iteration. More specifically, three increasingly “positively” (tensile) strained structures and three increasingly “negatively” (compressive) strained structures are generated along each search direction by systematically varying the corresponding parameters of the equilibrium structure. Next, DMC energies are obtained for the entire structure population. The shape of the resultant DMC PES is then used to locate a candidate minimum energy structure. The next iteration of the line search is performed identically to the first but starting from this new candidate minimum structure. This process is repeated until the statistical noise on the parameter error bars is within the desired limit.⁶⁸ Like all conjugate direction methods, a property of the method is the potential to converge after only a single iteration within a suitably quadratic region of the PES and using ideal search directions. In this work, converged DMC-based structural parameters were obtained after only three iterations and were averaged over the last two iterations to account for remnant statistical fluctuations in the lattice constant.

Phonon calculations meant to examine monolayer CrI₃’s spin–phonon coupling were performed using the frozen phonon method as implemented in the PHONOPY code⁶⁹ based upon DFT calculations performed in VASP.^57,59 We explored how CrI₃’s phonons change with different magnetic orderings, including nonmagnetic, ferromagnetic, and antiferromagnetic orderings, using both the LDA and LDA + U functionals. Our choice of U = 3 eV is further motivated below. 2 × 2 × 1 supercell structures with displacements were created from a unit cell fully preserving the material’s crystal symmetry. To properly analyze how the phonons vary in different magnetic phases, the supercell studied must be commensurate with the magnetic ordering. The 2 × 2 × 1 supercell selected was large enough to be commensurate with the FM and in-plane AFM ordering considered in this paper. Force constants were calculated using the structure files from the computed forces on the atoms. A part of the dynamical matrix was built from the force constants. Phonon frequencies and eigenvectors were calculated from the dynamical matrices with the specified q-points. Monolayer CrI₃ possesses D_3d point group symmetry, and hence, the phonon modes at the Γ point can be decomposed as G_D3d = 2A_1g + 2A_2g + 4E_g + 2A_1u + 2A_2u + 4E_u. Excluding the three acoustic modes (the doubly-degenerate E_u and A_1u modes) and noticing that each of the E_g and E_u modes is doubly degenerate, there are 21 modes in total.

As a starting point for our subsequent DMC calculations, we began by modeling the monolayer’s properties using standard density functionals, including the LDA, the GGA implementation of Perdew–Burke–Ernzenhof (PBE),⁷² PBEsol,⁷³ and a series of PBE + U functionals. Previous work has shown that the material properties can depend strongly on the functional employed, and thus, we initially sought a functional that could simultaneously predict the lattice constant, magnetic moment, and electronic bandgap of the monolayer. Since relatively few ab initio and experimental studies have been able to ascertain the properties of the monolayer, we fixed the internal geometry of the monolayer to the previously determined bulk values. As illustrated in Fig. 5, no single functional came close to simultaneously reproducing all three bulk properties: PBE + U = 1 eV most closely approximated the bulk lattice constant, PBE + U = 6 eV most closely approximated the bulk magnetic moment, and PBEsol most closely approximated the bulk bandgap. Even more glaringly, all three properties continuously increased with the U employed in our PBE + U calculations, strongly suggesting that the U’s employed could be continuously tuned to reproduce any individual quantity desired. This is illustrative of the sensitivity of CrI₃ properties to the variability in DFT functionals, which motivates the use of DMC throughout the rest of this work.

To select a DFT functional that could serve as a meaningful reference, we thus turned to DMC simulations. Based on Fig. 3, LDA + U wave functions with U’s of 2–3 eV minimize the DMC energy, so in all subsequent DFT analyses, we have employed LDA + U wave functions with U = 2–3 eV. This range of U values is consistent with our linear-response estimate of U_LR = 3.3 eV and also that employed in a wide range of materials and a recent DMC study of bulk CrI₃.⁶³ 

With a U = 2 eV value, we used LDA + U to relax the CrI₃ structure to determine its equilibrium geometry. As shown in Tables I and SI, LDA + U yields a lattice constant of 6.695 Å and an axial iodine–chromium–iodine bond angle (θ_I−Cr–_I) of 175.72°. The predicted lattice constant is within the range of lattice constants previously obtained using other density functionals, which range from 6.686 Å from previous DFT + U calculations⁴⁶ to 6.978 Å from GGA + U calculations⁷¹ (a variation of roughly 5%). In contrast, recent STM experiments point to an experimental monolayer lattice parameter of 6.84 Å, which is greater than our best LDA + U estimate (a point to which we will return below).

Given this variability in DFT-derived lattice constants, we employed DMC guided by a surrogate Hessian method to determine monolayer CrI₃’s equilibrium geometry with DMC-level accuracy. For monolayer CrI₃, the surrogate Hessian method yielded converged structural parameters to within the desired accuracy of 0.5% of the lattice parameter within three iterations.

As evidenced by the tightly bunched contour lines in Fig. 6, the potential energy surface around the DMC structural minimum is steep and highly harmonic. In terms of the Wyckoff parameters x and a (note that a ≡ a₀), the DMC PES is more shallow along the direction in which x and a simultaneously increase and steeper in the direction in which x increases but a decreases. In contrast, the PES is steeper along the direction in which z and a simultaneously increase and more shallow along the direction in which z decreases but a increases. Finally, although the relationship between z and x is slightly more anharmonic than the previous two relationships, the PES roughly tends to be more shallow as both z and x increase, but steeper in the direction in which z decreases, while x increases.

As bond angles and lengths consist of contributions from both the x and z Wyckoff parameters, a physical understanding of the DMC PES is enhanced by considering Fig. 7 of the main text and Figs. SII and SIII of the supplementary material. In these figures, slices of the DMC PES near the true minima are taken so as to illustrate (i) how the CrI₃ energy changes with the lattice constant, a₀; Cr–I bond length, d_Cr−I; and CrI₃ bond angles, θ₁ and θ₂ (see Fig. 1 for a visualization of how these quantities are defined) and (ii) that these structural quantities differ significantly between LDA + U and DMC. These figures reveal that the energy changes most rapidly as the lattice constant is varied. The energy minima become increasingly more shallow along the d_Cr−I, θ₁, and θ₂ directions, signifying that it is more energetically costly to isotropically strain the lattice than to more locally vary bond lengths and angles. The relatively low barriers to changing the bond distances and angles are what, moreover, make it challenging to accurately resolve CrI₃’s structure, a feature that CrI₃ holds in common with many 2D materials, as we have illustrated using DMC in our previous works.⁴² 

Ultimately, our line search converges upon an LS-DMC lattice parameter of a₀ = 6.87 Å, as presented in Fig. 7 and Table I. This lies within the range of previous DFT-derived estimates and is 2% larger than the DFT + U lattice constant discussed earlier. Nevertheless, what adds particular confidence to these results is that this independently derived lattice parameter is just 0.4% off from that recently obtained for a CrI₃ monolayer grown using molecular beam epitaxy and analyzed using scanning tunneling microscopy (see Li et al. in Table I).³² The fact that our surrogate Hessian approach was able to so closely reproduce an experimental value underscores both DMC and the surrogate Hessian approach’s accuracy for CrI₃. In addition to obtaining the monolayer lattice parameter, these calculations also yielded estimates of the chromium–iodine bond distance (d_Cr−I = 2.723 Å), as also presented in Fig. 7, and the monolayer bond angles (θ₁ = 90.4° and θ₂ = 175.4°), as tabulated in Table SI and Figs. SII and SIII of the supplementary material.

Although bond angles and distances have not been measured experimentally in the monolayer and are sparsely reported in DFT studies, our line search-predicted values all fall within 0.5% of experimental values for bulk CrI₃. The other bond angles and distances obtained using different DFT functionals tabulated in Table SI of the supplementary material differ from the experimental bulk values by up to 3%, further highlighting the robustness of our surrogate Hessian line search method for predicting structural parameters.

Using our DMC-optimized CrI₃ structure, we then proceeded to compute the magnetic moments on the Cr and I atoms in the hopes of fully and accurately resolving the magnetic structure of the monolayer for the first time. By integrating over the DMC spin density, as described in Sec. II C, we obtained a magnetic moment of 3.62 μ_B on each Cr and a moment of −0.15 μ_B on each I. Interestingly, this indicates that each Cr possesses a magnetic moment significantly larger than that of the +3 charge that one would naively assume based upon chromium’s typical oxidation state. Many previous DFT studies obtained an average magnetic moment of 3 μ_B across the entire CrI₃ unit cell, leading many researchers to conclude that all of the unit cell’s magnetism could be attributed to that from the Cr site.

While the sum of our magnetic moments for the Cr and three I atoms that constitute the unit cell also totals to 3 μ_B, we instead find that the iodine atoms additionally carry a significant magnetic moment, suggesting that CrI₃’s magnetism is more complex—and nonlocal—than previously assumed. As a check on these moments, we also employed LDA + U and DMC to compute the moments of a monolayer cleaved from the experimental bulk structure.⁴⁵ The moments obtained using this structure were different by up to 0.06 μ_B, supporting the assertion that DMC is capturing correlations missed by DFT (Fig. 8). Furthermore, the DMC moments obtained using this structure were within 0.005 μ_B of those from our line search-optimized structure (see Table I), a statistically insignificant discrepancy that strongly corroborates our findings.

In contrast, integration over the spin densities of the LDA + U-optimized structure yielded magnetic moments of 3.497 μ_B on each Cr and −0.099 μ_B on each I, meaning that the correlation accounted for by DMC calculations leads to an increased spin anisotropy in the unit cell. That correlation leads to an increase in spin anisotropy is further substantiated by the fact that we see an increase in the magnitudes of the moments in the DMC-optimized structure to 3.655 μ_B on each Cr and −0.153 μ_B on each I after applying VMC to our DFT trial wave function, followed by a further tuning to their final values upon applying DMC to our VMC trial wave function (see Fig. SIV of the supplementary material).

Finally, we used LDA + U = 3 eV in VASP to calculate atomic magnetic moments with and without spin–orbit coupling (SOC) to gauge the effect of SOC on the atomic magnetic moments. When SOC was turned off, we obtained moments of 2.848 μ_B and −0.073 μ_B on the Cr and I, respectively. When the SOC was turned on, the magnetic moment on the Cr grew to 2.856 μ_B, while the magnetic moment on the I remained −0.073 μ_B. This suggests that leaving SOC out of our DMC calculations should not drastically affect our magnetic moment results.

The discrepancies in the magnetic moments that we observe with varying monolayer geometries point to appreciable spin–phonon and spin–lattice couplings. Spin–phonon and spin–lattice couplings are associated with changes in the magnetic exchange interaction with changes in ionic motion or strain, respectively, and together can give rise to a large magneto-elastic effect. A material with a large spin–phonon coupling can be leveraged for magneto-caloric device applications.

Up to the lowest order, the shift in the phonon frequency due to coupling to the lattice can be attributed to ν = ν₀ + λ⟨S_i · S_j⟩, where ν₀ is the frequency in the paramagnetic case and λ is the spin–phonon coupling constant with i and j as site indices. For CrI₃, if we assume a nominal +3 oxidation state for Cr, the expectation value of this spin–spin correlation function equates to $∼9/4$⁠. Hence, for a given phonon frequency shift, we can estimate λ. The spin–spin correlation function is a constant for all modes in a given material. Hence, we simply report the frequency shifts in much of our discussion below. Nevertheless, we would like to note that for a given frequency shift, the effects of spin-anisotropy due to strong correlation, as revealed by our DMC calculations, should significantly increase the spin–spin correlation function, thereby effectively reducing the spin–phonon coupling.

Experimentally, the largest spin–phonon coupling is observed in a 5d perovskite oxide, corresponding to a frequency shift of 40 cm⁻¹.⁷⁴ In comparison, a frequency shift of ∼2.7 cm⁻¹ is observed for the E_g mode of the chromium-based 2D material Cr₂Ge₂Te₆ (CGT) around its T_c. This corresponds to a λ of ∼1.2 cm⁻¹. Spin–phonon couplings in the range of 0.27–0.4 cm⁻¹ have subsequently been observed in other transition-metal trihalides, again arising from the E_g modes.⁷⁵ 

In our phonon calculations, we initially fixed CrI₃’s lattice parameter to its bulk value. Considering AFM ordering as one specific realization of the paramagnetic phase, we find that the largest change in phonon frequency between the AFM and FM phases is for the E_g symmetric mode and is ∼4 cm⁻¹ (see Table II, modes 8 and 9). This is consistent with CrI₃ having a sizable spin–phonon coupling at finite temperatures, comparable to that of CGT and other 2D layered magnetic materials.

We next performed full geometry relaxations (using LDA + U = 3 eV in VASP) for the different magnetic configurations. The DMC calculations depicted in Fig. 3 helped us to narrow the range of U’s to employ in our phonon calculations down to U = 2–3. We ultimately employed U = 3 to obtain the data in Table II. Phonon frequencies obtained using U = 2 differed by at most 2.5 cm⁻¹ from those reported in the table, supporting this decision (see Table SII). We find the optimal lattice parameters for the FM, AFM, and NM phases of the CrI₃ monolayer to be 6.677, 6.656, and 6.617 Å, respectively (see Table II). A large, ∼1%, difference in lattice-parameters due to magnetic ordering is indicative of a non-trivial spin–lattice coupling and is consistent with unstable in-plane acoustic modes when ML-CrI₃ is forced to be in a non-magnetic state (see Fig. S5). Unstable acoustic modes have been observed in certain Heusler compounds.⁷⁶ A strong spin–lattice coupling, particularly involving soft in-plane acoustic modes, should lead to structural transitions across the Curie temperature. Indeed, bulk CrI₃ has been experimentally shown to transition from a rhombohedral to a monoclinic structure across the paramagnetic transition,⁴⁵ consistent with these findings.

We also find that the phonon frequencies shift to larger values across all of the different magnetic orderings when the lattice parameter is allowed to relax, as shown in Table II. Specifically, the phonon frequency change between the FM and AFM phases increases to ∼10 cm⁻¹ for the E_g mode (modes 8 and 9). This particular mode resembles shearing of the iodine planes and is the one with the largest spin–phonon coupling in CGT as well as in other Cr trihalides. Assuming a spin-eigenvalue S of 3.47ℏ/2 based upon our DMC calculations of the magnetic moment, this corresponds to a λ of ∼3.32 cm⁻¹. In comparison, the theoretically predicted spin–phonon coupling for CGT using a rigorous perturbation theory approach was 3.19 cm⁻¹.⁷⁷ 

These results demonstrate that even in the monolayer limit, CrI₃ possesses strong spin–phonon and spin–lattice couplings. While recent studies have demonstrated a magneto-optic effect in few-layer CrI₃, the large coupling of the magnetism to the Raman active phonon frequencies and lattice-strain that we observe in this study suggest that one could potentially use magnetic fields to modulate inelastically scattered light. Indeed, a very recent experiment demonstrates that an out-of-plane magnetic field change from −2.5 to 2.5 T leads to a rotation in the plane of the polarization of inelastically scattered light from −20° to +60°.⁷⁸ 

In summary, we have used a potent combination of first principles DFT and DMC calculations to produce some of the most accurate estimates of the electronic, magnetic, and structural properties of monolayer CrI₃ to date. Using a surrogate Hessian line search optimization technique combined with DMC, we were able to independently resolve the lattice and other structural parameters of CrI₃ to within a fraction of a percent of recently published STM measurements, an accomplishment given the up to 10% variability in previous DFT-derived estimates of the lattice parameter, depending upon the functional employed. Based upon the DMC-quality structure we obtained, we were then able to acquire a high-resolution monolayer spin density that showed each Cr atom to possess a magnetic moment of 3.62 μ_B and each I atom to have a moment of −0.145 μ_B, substantially larger moments than previously reported values that suggest that CrI₃ manifests substantial ligand magnetism. Given the Cr–I–Cr angle of 90°, this would indicate a superexchange stabilization of the ferromagnetic ordering. In conjunction with the expected large spin–orbit coupling on the I atom,⁴⁶ this could also give rise to a large ligand superexchange-dominated magnetic-anisotropy,^79,80 explaining recent observations of spin-waves and magnons in thin-film CrI₃.⁸¹ We, moreover, demonstrated that CrI₃ possesses remarkably strong spin–phonon coupling, with a predicted value as large as 3.32 cm⁻¹. This work, thus, demonstrates CrI₃’s promise for magnon-based spintronic applications,⁸² potentially with optical controls,⁷⁸ while also demonstrating the capability of DMC to accurately model the structural and magnetic properties of 2D materials.

See the supplementary material for additional figures and tables describing the time step extrapolations performed, lattice parameters, and magnetic moment calculations.

The authors thank Paul Kent, Kemp Plumb, Anand Bhattacharya, Ho Nyung Lee, Fernando Reboredo, Nikhil Sivadas, and the QMCPACK team for thoughtful conversations. The work by G.H., J.T., R.N., J.K., M.C.B., O.H., P.G., and B.R.; modeling and analysis efforts by D.S.; and the scientific applications of QMCPACK were supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division, as part of the Computational Materials Sciences Program and the Center for Predictive Simulation of Functional Materials. The preparation of this manuscript by D.S. was supported by the NASA Rhode Island Space Grant Consortium. This research was conducted using computational resources and services at the Center for Computation and Visualization, Brown University, and the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231.

The authors have no conflicts to disclose.

D.S. and G.H. made equal contributions to the paper.

The data that support the findings of this study are openly available in Blaiszik https://doi.org/10.1007/s11837-016-2001-3, Ref. 83, Blaiszik https://doi.org/10.1557/mrc.2019.118, Ref. 84 and the supplementary materials.