Single-layer CrI₂ as a magnetic semiconductor: a detailed first-principles study

Abstract

The study of magnetic monolayer (ML) structures has garnered increasing attention due to their potential for unveiling diverse physical phenomena in two-dimensional (2D) systems, as well as their promising applications in spintronics, optoelectronics, and magnonics. In this work, we present a comprehensive density functional theory investigation of the physical properties of monolayer chromium diiodide (CrI₂). As a starting point, we revisit the bulk CrI₂ system and show that its structural and vibrational properties, including Raman and infrared (IR) spectra, are highly sensitive to the treatment of van der Waals interactions, as modeled by different dispersion correction schemes. We then examine in detail the structural, electronic, magnetic, vibrational, and thermodynamic properties of the CrI₂ ML. Our results indicate that the antiferromagnetic (AFM) configuration corresponds to the ground state, with a magnetic moment of 3.8μ_B per Cr atom. The ML exhibits semiconducting behavior with an indirect band gap of 0.64 eV, and effective masses of −1.25m_e and 0.24m_e for the valence band maximum and the conduction band minimum, respectively. Additionally, we evaluate the magnetocrystalline anisotropy energy and identify key energy extrema in the AFM ML. A complete characterization of the vibrational modes at the Γ point is also provided, detailing their Raman and IR activity and atomic displacement patterns across different frequency ranges. The interplay of low-dispersive valence bands, semiconducting behavior, and robust antiferromagnetism suggests that CrI₂ MLs are promising candidates for next-generation spintronic and optoelectronic applications, offering both fundamental insights and technological potential.

---

1 Introduction

The study of low-dimensional materials has revealed a wide variety of exotic phenomena and promising applications across spintronics, electronics, and energy storage. Within this context, prior research has explored how impurities affect electronic properties and bonding configurations in graphene, highlighting the sensitivity of two-dimensional (2D) systems to atomic-scale modifications.¹ Similar insights have emerged from the study of silicon sulfide (SiS) and carbon monosulfide (CS) chains and monolayers (ML), where electronic bandgaps and topological features can be tuned through structural control and mechanical deformation.^2,3 Other geometries as SiS nanotubes⁴ and TiO nanocages close to the bulk chemical stoichiometry⁵ further expand this family of materials, exhibiting novel structural and electronic behaviors that enrich the landscape of functional materials at the nanoscale. From an application-oriented perspective, low-dimensional systems have also been explored for hydrogen storage. For example, borophene decorated with lithium atoms was found to be a promising candidate when van der Waals (vdW) corrections were properly considered in DFT calculations.⁶

Moreover, due to the exotic phenomena they exhibit, such as strong electron correlations and unusual phase behaviors, 2D systems have emerged as an exciting and active area of research compared to three-dimensional (3D) materials.^7–9 In fact, many 2D samples exhibit a wide variety of novel physical properties, including Moiré patterns,^10,11 superconductivity,¹² exotic spin and charge density waves,^13,14 magnetic anisotropy and spin polarization,¹⁵ as well as unique electrical behaviors¹⁶ and unexpected magnetic properties.¹⁷ These attributes have sparked considerable interest due to their potential applications in data storage, quantum computing devices, optoelectronics, and the development of ultra-thin spintronic devices.^18–20

Following the discovery of graphene, 2D materials exhibiting magnetic order have become the focus of intense investigation.^21–23 In particular, phenomena such as polarized valley states, topological quantum states, and topological superconductivity have been reported to arise from magnetic ordering in these systems.^24–26 Also, in the search for intrinsic magnetic semiconductors, transition-metal-based 2D oxides have drawn considerable attention due to role of structural distortions in tailoring their magnetic properties.²⁷

In contrast to the significant progress made in the study of 2D ferromagnets, research on low-dimensional antiferromagnetic (AFM) materials has been largely confined to exploring magnetic couplings between ferromagnetic (FM) vdW MLs.^21,23,28 In these cases, precise control over the number of MLs is crucial, as an even number of layers results in a global FM behavior, whereas an odd number yields a global AFM state. Moreover, since lattice vibrations (phonons) play a fundamental role in electrical, magnetic, thermal, and optical phenomena, 2D systems offer unique opportunities to tune the phonon energy spectrum and the associated material properties for a range of specific applications.

Among the growing array of available 2D systems, chromium diiodide (CrI₂) AFM MLs have recently begun to attract significant attention.²⁹ CrI₂ MLs can be synthesized either by direct growth at high temperatures^30–32 or through the decomposition of CrI₃,^30,33,34 emerging as a promising material for investigating strongly correlated magnetic systems. Very recently, experimental realization of CrI₂ in the single-layer limit has been reported, revealing phenomena such as degeneracy-lifted Hubbard bands, the characterization of its Mott insulating behavior via spectral weight transfer upon charge doping, and orbital-dependent Hubbard band splitting associated with a single charge defect.³⁰ These findings underscore the critical role of accurate theoretical modeling in predicting realistic functionalities.

Despite the burgeoning interest in the phenomenology of CrI₂ MLs, systematic studies that comprehensively explore their physical properties remain scarce. In particular, a recent contribution by Bo et al.²³ reported a comparative analysis of the magnetic ground states in a variety of transition metal dihalide MLs, including CrI₂, based on DFT + U and Monte Carlo simulations. While their study provides valuable insight into broader material trends, the scope, methodology, and focus differ substantially from our work. For instance, Bo et al. adopted a uniform Hubbard U = 2 eV across all systems, whereas we show that U = 0 yields structural parameters in closer agreement with experimental values for CrI₂, thus supporting a different methodological approach. Furthermore, they reported an FM ground state for the ML, while our total energy calculations identify the AFM configuration as more stable. Their assessment of magnetic anisotropy is limited to two directions, whereas we perform a broader directional analysis and reveal that the easy axis lies outside those high-symmetry directions. Additionally, Bo et al. do not explore vibrational or thermodynamic properties, which are central elements in our investigation. Therefore, while both works include CrI₂ as a system of interest, our study presents a more detailed and system-specific analysis, addressing structural, electronic, magnetic, vibrational, and thermal behavior in a unified framework.

Accordingly, in this study, we present a detailed first-principles investigation of the physical behavior of the CrI₂ system, with a special focus on the ML.

The paper is organized as follows: In sec. 2, we describe the computational procedure. In sec. 3.1, we validate our computational parameters by presenting results for the CrI₂ bulk phase and comparing them with available references, focusing on the monoclinic, orthorhombic, and hexagonal phases. Our analysis aims to identify the lowest energy states, which depend on the treatment of vdW interactions; hence, particular attention is given to the influence of different vdW approaches. In sec. 3.2, we present the structural, electronic, vibrational, and magnetic properties of the CrI₂ ML, including an analysis of its magnetic anisotropy energy and its vibrational and thermal properties. Finally, sec. 4 presents some final remarks regarding practical challenges and prospective solutions and sec. 5 summarizes our findings.

2 Systems, methods and computational details

Both the bulk CrI₂ and its ML form are considered in this study. A comparison was made among the bulk monoclinic (M-phase), orthorhombic (O-phase), and hexagonal (H-phase) structures. Fig. 1 displays schematic representations of the crystal structures for the M-, O-, and H-phases in bulk, along with top views for structural characterization.

Fig. 1 Unit cell representation of the fully relaxed monoclinic (a), orthorhombic (b), and hexagonal (c) bulk CrI₂ systems. The blue spheres represent Cr atoms, while the purple balls correspond to I atoms. The unit cells are delineated by black lines. Top views of the monoclinic (d), orthorhombic (e), and hexagonal (f) structures are also shown to highlight the different Cr–I distances.

Calculations were carried out within the framework of Density Functional Theory (DFT)³⁵ using a plane-wave basis set, as implemented in the Vienna Ab initio Simulation Package (VASP).^36–39 The Projector-Augmented Wave (PAW) method^40,41 was employed with a plane-wave energy cutoff of 500 eV. Calculations were performed by increasing the energy cutoff to 600 eV. The total energy difference between E_cutoff = 600 and E_cutoff = 500 (the latter of which was used) was 0.096 meV per atom. To treat the exchange–correlation energy and facilitate comparisons, two different formulations of the Generalized Gradient Approximation (GGA)⁴² were used in the bulk case: the Perdew–Burke–Ernzerhof (PBE)⁴³ and PBEsol⁴⁴ functionals. Also, aiming to the same purpose in the ML case, calculations were performed by using the PBE and HSE06 hybrid functionals.⁴⁵

In the bulk case, van der Waals interactions were accounted for by incorporating various dispersion correction schemes available in VASP;^46,47 namely: Nonlocal vdW-DF functional (vdW-DF),^48,49 Tkatchenko–Scheffler method (TS),⁵⁰ Tkatchenko–Scheffler method with iterative Hirshfeld partitioning (TS-Hi),^50–52 Many-body dispersion energy method (MBDE),⁵³ Density-dependent energy correction method (DDEC),^54,55 and Grimme's method with Becke–Johnson damping function (GBJDF).⁵⁶ a detailed evaluation of these vdW methods, in terms of both energy and structural parameters, is provided in the following section.

The Monkhorst–Pack scheme was employed for integrations in the Brillouin zone (BZ)⁵⁷ with meshes of 7 × 13 × 7 (for the M-phase), 13 × 7 × 4 (for the O-phase) and 13 × 13 × 5 (for the H-phase). A 7 × 13 × 1 mesh was used for the ML case. These meshes correspond to sets of 175, 56, 63 and 28 special k-points in the irreducible BZ, respectively. Denser meshes were used for the electronic structure calculations.

The structures were fully relaxed until the total energy change was below 10⁻⁹ eV and the forces on each atom were less than 0.001 eV Å⁻¹ in the bulk and 0.0002 eV Å⁻¹ in the ML case, along each Cartesian direction. Such high convergence in the forces is essential for calculating the dynamical matrix via the direct force constant approach,^58,59 which enabled the identification of the irreducible representations and the character of the phonon modes at the Γ point. For the ML CrI₂ calculations, periodic images were separated by 20 Å of vacuum to eliminate spurious interactions.

The PHONOPY⁶⁰ code was employed to assign the Raman modes at the Γ point, while the VASP_RAMAN⁶¹ code was used to compute the Raman intensities, which are related to the derivative of the dielectric susceptibility with respect to the normal mode coordinate. For the CrI₂ ML, phonon dispersions were calculated using the finite displacement method on 2 × 2 × 1 supercells.

The anisotropy constants were determined from non-collinear calculations of the system's energy. In these calculations, the Cr spins were varied over the polar angle θ in the [0,π] range and the azimuthal angle ϕ in the [0,2π] range with 10-degree increments. An expansion of the energy function E(θ,ϕ) in spherical harmonics was then used to extract the anisotropy constants.

All DFT calculations were performed with tight convergence criteria to ensure the accuracy and reliability of the computed structural, electronic, magnetic, and vibrational properties.

3 Results and discussion

3.1 Bulk structure

Fig. 1 depicts the unit cell representations for the fully relaxed bulk structures considering their M-, O-, and H-phases. These structures consist of stacked atomic monolayers weakly bonded by van der Waals (vdW) interactions. The monolayers exhibit a similar triple-layer stacking of I–Cr–I atoms, where the Cr atoms form a hexagonal lattice.

As a primary result, the antiferromagnetic (AFM) state is found to be the lowest-energy magnetic configuration for the M- and O-phases, which is in agreement with previous reports.⁶² In contrast, the H-phase exhibits a ferromagnetic (FM) ground state, with energy above that of the M- and O-phases.

Notable effects of the vdW scheme on the structural parameters of these phases are observed. Table 1 summarizes the bulk properties of the AFM states for the M- and O-phases, as well as those of the FM H-phase, computed within the PBE approximation for different vdW schemes.

Table 1 Structural parameters of the bulk AFM monoclinic (M-), AFM orthorhombic (O-), and FM hexagonal (H-) phases computed for each vdW scheme. ΔE represents the energy difference relative to the lowest-energy phase in each scheme. a, b, and c denote the lattice parameters, β is the angle between the b and c cell vectors, and γ is the angle between a and c. Cr1–I and Cr2–I refer to the distinct Cr–I bond lengths. The experimental values are listed in the bottom section. The lowest-energy phase for each vdW scheme is marked with an asterisk (*), with its energy set to 0.0 meV per atom

	ΔE (meV per atom)	a (Å)	b (Å)	c (Å)	β (°)	γ (°)	Cr1–I (Å)	Cr2–I (Å)
vdW-DF
M-phase	1.05	8.09	4.04	7.91	117.9		2.81	3.51
O-phase (*)	0.00	4.03	8.11	13.92	90.0		2.81	3.52
H-phase	2146	4.24	4.24	13.82		120.0	2.95	2.96
TS
M-phase (*)	0.00	7.48	3.95	7.56	115.1		2.74	3.24
O-phase	0.01	3.95	7.49	13.67	90.0		2.74	3.24
H-phase	14.10	4.09	4.09	13.47		120.0	2.85	2.86
TS-Hi
M-phase (*)	0.00	7.54	3.92	7.74	115.4		2.73	3.29
O-phase (*)	0.00	3.92	7.56	13.96	90.0		2.73	3.29
H-phase	19.10	4.06	4.06	13.80		120.0	2.85	2.85
MBDE
M-phase	0.68	7.31	3.90	7.44	114.0		2.73	3.16
O-phase (*)	0.00	3.89	7.34	13.57	90.0		2.73	3.17
H-phase	6.88	4.01	4.01	13.34		120.0	2.83	2.84
DDEC
M-phase	0.79	7.44	3.91	7.49	114.5		2.73	3.20
O-phase (*)	0.00	3.91	7.46	13.61	90.0		2.73	3.20
H-phase	4.29	4.04	4.04	13.38		120.0	2.83	2.84
GBJDF
M-phase	0.21	7.26	3.90	7.28	114.0		2.73	3.13
O-phase (*)	0.00	3.89	7.30	13.26	90.0		2.72	3.13
H-phase	3.82	4.01	4.01	13.03		120.0	2.82	2.83
Experiment^63
M-phase		7.54	3.93	7.50	115.5		2.74	3.24
O-phase		3.92	7.56	13.55	90.0		2.72, 2.75	3.22, 3.24

---

A general good agreement is observed between the obtained results and experimental lattice parameters, with the TS scheme providing the most accurate values. Furthermore, the lowest-energy phase varies depending on the vdW scheme: either the M- or O-phase can be the ground state, and in the case of the TS-Hi scheme, they are nearly degenerate. This strong dependence on the vdW scheme highlights the importance of careful selection when analyzing the quantitative results and conclusions regarding these systems.

Our results show that including different dispersion correction schemes produces variations in the c parameter of up to ∼13%, which are expected since corrections for van der Waals interactions primarily impact the out-of-plane direction (c-axis), which is typically associated with inter-laminar stacking. Interestingly, non-negligible variations of up to 2.5% have been found in the b parameter. Although b corresponds to a periodic internal direction of the ML, van der Waals interactions can affect directions within the plane due to at least three interrelated factors: (i) global geometric coupling: when different dispersion corrections are applied, total relaxation of the structure occurs, affecting not just the c direction, but also the internal atomic positions to reach a new energy minimum. This can induce small modifications in the periodic in-plane directions, even if they are not directly related to van der Waals interactions. (ii) Angular distortion: the significant change observed in β (from 113.7° to 117.9°) suggests generalized distortion of the cell. This angular change alters the actual projections of the cell vectors onto the basal plane, also affecting the in-plane dimension. (iii) Changes in Cr–I distances: these distances also adjust slightly for different functionals, suggesting a redistribution of the bonding environment (slight charge redistribution). Although subtle, these changes affect the balance of internal forces and, consequently, the entire geometry of the cell. Together, these effects account for the sensitivity of the b parameter to different dispersion corrections.

From the other side, a recent study focused on this system⁶² and reported an energy difference of 0.13 meV per atom between the O- and M-phases when using the Tkatchenko–Scheffler method. The authors attributed the experimental difficulty in distinguishing these phases to this small energy difference. In our higher-accuracy study, an even smaller energy difference is found for the same method, with similarly small energy differences across other vdW functionals.

To further explore the conclusions of ref. 62, we computed the Raman and infrared spectra of bulk CrI₂ using different vdW schemes, as shown in Fig. 2. Some differences between the O- and M-phases are noted. In particular, the Raman spectrum of the O-phase exhibits three peaks below 100 cm⁻¹, whereas only a single peak is observed for the M-phase in the same region. Small shifts in peak positions occur depending on the vdW scheme. Additionally, a new peak appears at approximately 225 cm⁻¹ for the O-phase with the DDEC scheme and for the M-phase with the TS scheme. Conversely, the infrared spectra are largely similar, remaining within experimental uncertainty, except for the DDEC scheme, which exhibits a peak at around 290 cm⁻¹ for the O-phase and 225 cm⁻¹ for the M-phase.

Fig. 2 Raman spectra for (a) the O- and (b) the M-phases with no vdW interaction (green dashed line) and for selected vdW schemes (blue, orange, and red signals). Infrared spectra for (c) the O- and (d) the M-phases obtained using the same selected vdW schemes.

These results suggest that, despite typical experimental uncertainties of approximately 50 cm⁻¹, Raman spectroscopy could serve as a means of distinguishing between the M- and O-phases in bulk CrI₂ systems, provided the experimental spectrum matches one of the computed ones and the best vdW scheme is identified. However, determining the lowest-energy state via this approach remains challenging unless the experimental uncertainty is sufficiently low to resolve the small peak shifts. The infrared spectrum could also aid in phase identification if the experimental data align best with the computed spectra obtained using the DDEC scheme.

3.2 CrI₂ monolayer

In this section, we focus on the primary objective of this study: to explore the properties of the CrI₂ ML. To this end, we extracted a single layer from each of the M- and O-bulk phases. As an initial observation, the atomic structures in the ML limit corresponding to phases M- and O- are structurally equivalent. Indeed, they consist of alternating I–Cr–I layers with identical Cr1–I and Cr2–I distances. This is evident from an analysis of the radial and angular distribution functions, which are identical for both MLs, indicating equivalent local atomic environments. This fact is also evident in the geometry; we observe that it is possible to map the atomic positions of one ML onto the other by means of a rigid translation within the basal plane. These findings suggest that, in the bulk, the difference between the M- and O-phases lies mainly in how these equivalent MLs are stacked to form the 3D structure (different spatial symmetries and inter-laminar orderings). At the monolayer limit, however, both structures are equivalent. Consequently, from this point onward, we will refer to the ML without distinguishing between the M- and O-phases.

Moreover, the vdW interaction does not significantly affect the ML properties, since the bonds exhibit a mixed ionic-covalent character (as discussed below); hence, the ML results are presented without any vdW corrections.

Fig. 3 shows a schematic representation of the fully relaxed unit cell of the ML system, which comprises 4 I atoms and 2 Cr atoms. In this structure, the Cr layer is sandwiched between two I layers, with lattice parameters a = 7.81 Å and b = 3.94 Å, and a Cr–I distance of 2.74 Å. Furthermore, as shown in Fig. 3(a), an isosceles triangle is formed by the three closest iodine atoms in the top layer, with distances d₁ = 3.93 Å and d₂ = d₃ = 4.33 Å. These values are in good agreement with the experimental measurements reported by P. Li et al.⁶⁴ (d₁ = 4.0 Å and d₂ = d₃ = 4.2 Å), lending confidence to the DFT methodology used to compute the ML properties.

Fig. 3 Unit cell representation for the fully relaxed single-layer system. The blue spheres represent Cr atoms, while the purple balls represent I atoms. (a), (b), and (c) show different views of the structure. The green lines in (a) indicate the various I–I distances obtained in the relaxed ML, with d₂ = d₃ = 4.33 Å.

The structural stability of the ML was confirmed by analyzing the frequencies of its vibrational modes to ascertain whether the relaxed structure resides at a minimum or at a saddle point on the potential energy surface. Fig. 4 presents the vibrational mode frequencies; all 18 normal modes (corresponding to the system's degrees of freedom) exhibit positive (or zero) frequencies, confirming the dynamic stability of the ML.

Fig. 4 Vibrational modes for the AFM semiconducting single-layer CrI₂. All 18 normal modes are greater than or equal to zero, indicating structural stability.

Bader charge analysis provides insights into the charge distribution and magnetic moments in the CrI₂ ML, revealing the key aspects of its electronic and magnetic properties. The computed charges, valence values, and magnetic moments for each atom are summarized in Table 2. A net charge transfer from the Cr atoms to the I atoms is observed. This suggests that the Cr–I bond has a partially ionic nature, which is consistent with iodine's higher electronegativity. Specifically, the Cr atoms have a Bader charge of approximately 5.1e, indicating a loss of about 0.9e from their nominal six valence electrons (3d⁵4s¹). The I atoms, on the other hand, have a charge of about 7.5e, meaning that they have gained about 0.5e with respect to their seven valence electrons (5s²5p⁵).

Table 2 Bader analysis for the AFM single-layer system. The table shows the charge, valence, and magnetic moment for each atom

Atom	Charge (e)	Valence (e)	Magnetization (μB)
Cr1	5.1	0.9	−3.8 ± 0.1
Cr2	5.1	0.9	3.8 ± 0.1
I1	7.5	−0.5	0.1 ± 0.1
I2	7.5	−0.5	0.1 ± 0.1
I3	7.5	−0.5	−0.1 ± 0.1
I4	7.5	−0.5	−0.1 ± 0.1

---

This charge distribution suggests that the Cr–I bond has a mixed ionic–covalent character. The ionic component is evident in the charge transfer from Cr to I, while the covalent component arises from the incomplete electron sharing, as indicated by the partial charge on the I atoms. This interpretation is consistent with the charge density visualization in Fig. 5, which shows significant charge accumulation around the Cr and I atoms, and some charge distribution present along the Cr–I bonds.

Fig. 5 (a) 3D charge density for the single-layer CrI₂ within the unit cell. (b) 2D projection of the charge density. (c) Schematic representation of the magnetic moment directions in the AFM configuration, indicating an out-of-plane easy axis for the magnetic anisotropy energy of the CrI₂ single layer (green arrows). The light-blue plane represents the plane formed by the Cr atoms and their I neighbors, with the easy axis normal to this plane.

Regarding the magnetic structure, the minimum-energy configuration corresponds to the (stripy) AFM state, which is 6.12 meV lower in energy than the FM configuration. The nature of the AFM ordering was determined by evaluating different magnetic configurations, namely Néel, Stripy-y, Zigzag-y and Zigzag-x (equivalent to Stripy-x). The results indicate that the Stripy-y configuration is the most stable, corresponding to the state with the lowest total energy. The energy differences with respect to this configuration are 13 meV for Néel, 52.5 meV for Zigzag-y, and 33.2 meV for Zigzag-x, confirming the energy preference for Stripy AFM ordering along the y direction. The related energy difference (6.12 meV) relative to the FM state corresponds to a thermal energy of approximately 71 K. This suggests that at elevated temperatures, magnetic transitions or phase coexistence could take place. This scenario is consistent with the estimated Néel temperature (about 40 K) obtained from our calculated exchange constants, reinforcing the possibility that thermal fluctuations may stabilize competing magnetic configurations near or above this temperature. Furthermore, external factors such as mechanical strain, applied magnetic fields, or electron doping could alter the magnetic exchange couplings (J) between Cr atoms, potentially favoring configurations other than the predicted AFM ground state. It is also worth noting that the existence of non-collinear magnetic solutions is conceptually possible due to the presence of complex magnetic interactions, such as isotropic exchange, anisotropies, and, potentially, the Dzyaloshinskii–Moriya interaction, which can favor non-collinear configurations in systems with suitable symmetry and magnetic frustration. However, such calculations are significantly more computationally demanding and complex. Therefore, at this stage of the study, we focused on collinear magnetic solutions as a reasonable and efficient approximation of the system. The exploration of non-collinear states remains an open avenue for future work.

Furthermore, the computed atomic magnetic moments indicate that the Cr atoms dominate the magnetism in the ML, each exhibiting a magnetic moment of 3.8(±0.1)μ_B, whereas the I atoms exhibit much smaller moments (less than 0.1μ_B). This pronounced difference in magnetic moments suggests a relative strong localization of magnetism around the Cr atoms. Although the iodine atoms do not significantly contribute to the overall magnetism, they play a crucial role as ligands by surrounding the Cr magnetic centers and contribute to the stabilization of the crystal structure by mediating magnetic interactions among Cr atoms via superexchange mechanisms.

The density of states (DOS) and band structure of the CrI₂ ML provide crucial insights into its electronic and magnetic properties. Fig. 6(a) shows the total DOS along with the contributions from the s, p, and d orbitals, exhibiting a symmetric distribution of up and down spins confirming the AFM ordering in the system. A band gap of 0.64 eV is obtained, indicating the semiconducting nature of the ML. Regarding orbital contributions, the s states dominate at energies below −12 eV, while the p orbitals are significant between about −5 and −1.7 eV, and the valence and conduction bands are predominantly governed by the Cr d orbitals (as observed from the PDOS of the Cr atoms represented by the yellow and violet regions in Fig. 6(a)). Consequently, these d orbitals, which underpin the material's magnetic and electronic properties, are likely to play a key role in its transport behavior.

Fig. 6 (a) Total density of states (DOS) for the CrI₂ monolayer (light green), together with the partial density of states (PDOS) of the d orbitals of the two Cr atoms (orange and violet). The inset shows the partial DOS for the s (green), p (blue), and d (red) orbitals on a reduced scale. (b) Band structure for the CrI₂ monolayer, highlighting the valence band maximum (VBM) at the Γ point, the conduction band minimum (CBM) at the midpoint of the Γ–S pathway, and the Fermi level (red dotted line).

These results confirm that the CrI₂ ML exhibits AFM semiconductor behavior, characterized by antiparallel spin polarization in the two Cr atoms and a finite band gap in the majority spin channel of each site. This classification as an AFM semiconductor is particularly interesting for spintronic applications since these materials can provide magnetic order without stray fields and greater thermal stability, which makes them suitable for high-density, low-power logic and memory devices.

Fig. 6(b) presents the band structure, which corroborates the semiconducting character of the CrI₂ ML by exhibiting an indirect band gap. The VBM is located at the Γ point, while the CBM is found along the Γ–S pathway. The effective masses determined for the VBM and CBM are −1.25m_e and 0.24m_e, respectively.

To evaluate the influence of the exchange–correlation functional on the electronic properties of the CrI₂ ML, we performed calculations using the HSE06 hybrid functional on a geometry optimized with the PBE functional. Although a full relaxation with the HSE06 hybrid functional was not performed due to its high computational cost, the remaining residual forces (below 0.10 eV Å⁻¹) and internal pressures (in the sub-kilobar range) ensure reliable analysis.

The band gap obtained with the HSE06 hybrid functional is significantly larger (2.4 eV) than the PBE value (0.64 eV) value, reflecting the known gap underestimation by semi-local functionals. The PDOS around the Fermi level is qualitatively similar for both functionals (see the SI), suggesting minimal changes in the effective mass. Furthermore, the local magnetic moment of Cr remains at ∼3.8μ_B with the HSE06 hybrid functional, and the AFM order remains robust. Furthermore, we also considered the DFT + U approach during our initial methodological assessment. However, as briefly mentioned in the introduction, we found that for CrI₂, a value of U = 0 eV provided structural parameters in better agreement with available experimental data. Therefore, given the computational expense of the HSE06 hybrid functional, particularly for dense k-point sampling or full relaxations, the PBE functional was selected as the most appropriate choice for the consistent structural and electronic comparison presented in this work. Thus, the remainder of the study was conducted with the PBE functional, which reliably captures the semiconducting and AFM character of the system, as mentioned above.

The presence of a low-dispersive valence band, coupled with pronounced peaks in the DOS (indicative of localized states), along with the semiconducting and AFM properties, make the CrI₂ ML a promising candidate for applications in both spintronics and optoelectronics.

3.3 Magnetic anisotropy energy of the ML

A detailed analysis of the magnetic anisotropy reveals complex features that are fundamental to understanding the magnetic properties. The Magnetocrystalline Anisotropy Energy (MAE), obtained by fitting spherical harmonics (with l_max = 2) to the DFT results, can be expressed as:

E_MAE(θ, ϕ) = K₀ + K₁ cos² θ + K₂ sin(2θ) cos ϕ + K₃ sin² θ cos(2ϕ), (1)

where the coefficients are K₀ = 1.266 meV, K₁ = −0.082 meV, K₂ = 0.695 meV, and K₃ = −0.048 meV. This equation describes how the anisotropy energy varies with the magnetic orientation angles θ (polar) and ϕ (azimuthal). Here, the polar angle θ was measured with respect to an axis perpendicular to the monolayer, while the azimuthal angle ϕ was measured within the plane defined by the Cr atoms relative to the direction of the lattice parameter a.

The easy magnetization axes were identified in the directions θ = 59.45°, ϕ = π and θ = 120.56°, ϕ = 0, which are perpendicular to the plane formed by a Cr atom and its nearest I neighbors, as illustrated in Fig. 5(c). The out-of-plane orientation is crucial for maintaining stable magnetic order in two-dimensional materials.

The 3D and 2D representations of the MAE, as shown in Fig. 7, reveal pronounced minima corresponding to the easy magnetization axes, along with complex energetic variations that indicate the presence of multiple metastable states. This feature could imply rich magnetic domain dynamics and non-trivial magnetic states. Table 3 provides additional details on the critical points (CPs), including easy-axis minima, saddle points, and local maxima.

Fig. 7 (a) 3D energy anisotropy landscape and (b) its 2D projection for the CrI₂ monolayer.

Table 3 Critical points of the CrI₂ anisotropy energy. The first column indicates the type of critical point, the second column shows its location (in degrees), and the third column gives the energy in meV

Type	CP (°)	Energy (meV)
Saddle point	(0, 0)	1.184
Saddle point	(30.57, 0)	1.595
Local maximum	(90, −90)	2.114
Local minimum	(120.56, 0)	0.007

---

Identification of these critical points suggests the existence of multiple pathways for magnetic reorientation, which has important implications. Indeed, flexibility in manipulating the magnetic direction allows for more precise control, enabling the selection of reorientation paths that may be faster, more stable, or more energy efficient. Furthermore, the complex anisotropic behavior in the monolayer, evidenced by higher-order terms and angular dependence in eqn (1), can give rise to non-trivial magnetic states that can manifest as non-collinear spin configurations, complex magnetic textures, or multiple metastable states. This could be promising for developing multivalent magnetic memory technologies or for designing devices with reconfigurable magnetic behavior. These findings open the door for future studies to explore such aspects from a dynamic perspective, including the possibility of non-collinear configurations, complex magnetic textures, or even topological effects under specific conditions.

3.4 Vibrational and thermodynamic properties

According to the harmonic approximation, a system is dynamically stable only if all its phonon frequencies are positive. Thus, it is essential to determine the phonon dispersion relationships across the entire Brillouin zone and compute the corresponding density of states to confirm the dynamic stability of the structure. Initially, the phonon frequencies of the monolayer were characterized using group theory. Table 4 summarizes the key results, providing information on how the phonon modes, representing the atomic vibrations, transform under the symmetry operations of the crystal's point group at the Γ point. The results indicate that the structure belongs to the space group Pm and the point group m.⁶⁵ In addition, two types of irreducible representations, A′ and A′′, are identified; these differ according to how the phonon modes transform under the monolayer's symmetry operations, thereby helping to distinguish the Raman-active modes from the infrared (IR) modes.

Table 4 Phonon frequencies at the Γ point (in THz). “R” and “IR” indicate the Raman-active and infrared-active modes, respectively. Raman modes are labeled as A′, while A′′ represents the IR modes. The mode labels correspond to those in Fig. 4

Space group: Pm
Point group: m
Mode	Frequency	Irrep
[001]	0.499	A′ (R)
[002]	0.929	A′ (R)
[003]	0.986	A′′ (IR)
[004]	1.248	A′ (R)
[005]	1.381	A′ (R)
[006]	2.159	A′ (R)
[007]	2.202	A′ (R)
[008]	2.672	A′′ (IR)
[009]	2.691	A′′ (IR)
[010]	3.432	A′ (R)
[011]	3.514	A′ (R)
[012]	6.144	A′′ (IR)
[013]	6.161	A′′ (IR)
[014]	6.736	A′ (R)
[015]	6.877	A′ (R)

---

To identify the Raman-active modes, one must select the modes that exhibit a nonzero character under symmetry operations excluding the inversion center (i.e., those with parity different from −1).⁶⁶ At the Γ point, these symmetry operations are defined by the rotation matrices and the corresponding transformation matrix. Modes satisfying this condition are labeled with the irreducible representation A′. Conversely, the infrared-active modes are those that exhibit a nonzero character under symmetry operations corresponding to the inversion center (parity equal to −1), and these are labeled as A′′. Table 4 lists each monolayer mode with its respective frequency and irreducible representation.

Once the phonon frequencies are characterized, the vibrational patterns for each mode can be determined. As shown in Fig. 8, the Raman-active modes (Irrep A′) exhibit vibrations both within the plane (x–y) and out-of-the-plane, indicating that these modes can induce changes in the system's dipole moment.^67,68

Fig. 8 Raman-active vibrational modes. The blue spheres represent Cr atoms, while the purple spheres represent I atoms. The red arrows indicate the vibration direction. The labels correspond to the modes listed in Table 4. Modes [016], [017], and [018] have zero frequency, as shown in Fig. 4.

In contrast, the IR-active modes are predominantly confined to vibrations within the (x–y) plane. As illustrated in Fig. 9, during these vibrational modes, the system's dipole moment remains unchanged.

Fig. 9 Infrared-active vibrational modes. The blue spheres represent Cr atoms, while the purple spheres represent I atoms. The red arrows indicate the direction of vibration. The labels correspond to the modes listed in Table 4.

Fig. 10 displays the phonon dispersion relationship across the Brillouin zone, where three acoustic branches are identified: two transverse acoustic (TA) branches and one longitudinal acoustic (LA) branch. These acoustic branches appear at very low frequencies on the order of 0.0196 cm⁻¹. A gap in the phonon spectrum is observed between 115.31 and 202.14 cm⁻¹, which can be attributed to the distinct contributions and transitions of the vibrational modes associated with each atom.⁶⁹ Furthermore, at low frequencies (w < 115.31 cm⁻¹), the vibrational contributions from the I atoms dominate, whereas at higher frequencies (w > 202.14 cm⁻¹), the Cr atoms contribute more significantly to the normal modes.

Fig. 10 Phonon dispersion relationship. The green triangles and blue diamonds at the Γ point correspond to the frequencies of the Raman and IR modes, respectively. The acoustic branches—LA (longitudinal) and TA (transverse)—are explicitly indicated.

Fig. 11 presents the total and partial density of states (PDOS) for the phonon frequencies. The PDOS reveals that the I atoms contribute predominantly to the low-frequency modes (w < 3.457 THz), while the Cr atoms dominate the high-frequency region, consistent with previous observations. Moreover, the absence of negative frequencies confirms the dynamic stability of the structure.

Fig. 11 Partial density of states for phonon modes: Cr (green solid line) and I (red dashed line with yellow shading). The inset shows the total phonon DOS.

Thermodynamic properties were evaluated by calculating the specific heat (c_V) and entropy (S) as functions of temperature, as shown in Fig. 12. These calculations highlight the phonon contribution to these properties. It is evident that above 300 K, c_V saturates at approximately 150 J K⁻¹ mol⁻¹, indicating that all vibrational modes are fully excited at high temperatures, while phonon contributions are found to be negligible as the temperature decreases. At around 291 K, the c_V begins to decrease slightly; this temperature corresponds to a frequency of 202.14 cm⁻¹, where the phonon density of states is low (see Fig. 11), resulting in a reduced phonon contribution. Conversely, at temperatures below 166 K (corresponding to frequencies below 115.08 cm⁻¹), the c_V decreases more substantially, possibly due to the higher density of low-frequency phonon states.

Fig. 12 Thermodynamic properties of the CrI₂ monolayer: entropy S (red continuous line), specific heat c_V (dot-dashed blue line), harmonic phonon energy (HPE; dashed green line), and Helmholtz free energy (HFE; yellow circles).

The entropy exhibits a more pronounced decrease at low temperatures. This behavior is expected because entropy, which reflects the number of accessible excitation modes, diminishes as the system cools and fewer vibrational states are thermally populated, indicating reduced disorder and fewer available configurations.

4 Final insights

Although this study does not address experimental aspects directly, the technological potential of the CrI₂ ML merits some final remarks regarding practical challenges and prospective solutions.

First, synthesizing a stable, free-standing CrI₂ ML faces several challenges. These include the high chemical reactivity of iodine, susceptibility to environmental degradation, and the need for precise control over structural quality to ensure reproducibility of the predicted properties. Such challenges are common in 2D halide-based materials and may require encapsulation strategies or inert environments during fabrication and characterization. Regarding the AFM ordering, it is well known that substrate interactions can significantly alter the magnetic behavior of 2D materials. Substrates may induce mechanical strain, charge transfer, or electronic hybridization, which could suppress or destabilize AFM ordering, particularly in the case of metallic or reactive substrates. Therefore, careful substrate selection is critical for preserving the intrinsic properties observed in the free-standing case. Ideal substrates include hexagonal boron nitride (h-BN), SiO₂ with suitable surface passivation, or other inert van der Waals materials, which can minimize electronic interference while providing mechanical support. Although our calculations for the free CrI₂ ML provide a fundamental basis, selecting the appropriate substrate is essential for stabilizing and preserving magnetic behavior in experimental practice.

Second, understanding the stability of the AFM ground state of the CrI₂ ML in the face of external perturbations, such as mechanical strain, electric fields, and interaction with substrates, is crucial for grasping and applying the material in practice. The application of strain is known to modify magnetic exchange constants and anisotropy, which can stabilize or alter AFM ordering depending on the nature and magnitude of the applied deformation. Similarly, external electric fields can influence electronic distribution and magnetic anisotropy. However, their effects are typically more subtle in systems with strong intrinsic magnetic coupling, such as CrI₂. Although the present study analyzed the system under ideal, disturbance-free conditions to establish a reliable baseline, evaluating the impact of these variables is an important challenge and the subject of future work. Such studies would improve our understanding of the robustness of the AFM state and optimize the conditions for device integration.

Finally, the ability of the CrI₂ ML to support magnon transport and applications in magnonic devices is of great interest for developing spintronic technologies based on 2D materials. Preliminary calculations on magnonic excitation suggest the existence of well-defined spin-wave excitation with compatible energies and dispersion for stable magnon transport at room temperature. These results suggest that the CrI₂ ML could be a promising material for magnonic applications. Additionally, the observed AFM order and anisotropic properties favor the coherence and manipulation of magnons, which are essential for designing efficient magnonic devices. These results motivated our group to conduct a new study focused on characterizing the magnonic properties, dynamics, and response to external perturbations. The goal is to evaluate the feasibility of CrI₂ as a platform for magnon-based devices.

5 Conclusions

In summary, we conducted high-precision DFT calculations to investigate the structural, electronic, magnetic, vibrational, and thermodynamic properties of CrI₂. For the bulk, we found that its physical properties and vibrational spectra strongly depend on the treatment of van der Waals interactions. Our results highlight the potential of Raman spectroscopy to detect subtle structural variations, especially when compared to spectra computed using density-dependent dispersion corrections.^54,55

In the monolayer, the antiferromagnetic ground state features a magnetic moment of 3.8μ_B per Cr atom. The system is a semiconductor with an indirect band gap of ∼0.64 eV and effective masses of −1.25m_e and 0.24m_e for the valence and conduction band edges, respectively. The MAE landscape reveals multiple critical points, suggesting nontrivial magnetic behavior. Vibrational mode analysis identifies Raman- and IR-active branches, while the phonon dispersion confirms dynamic stability.

Overall, the interplay of antiferromagnetism, semiconducting behavior, and rich vibrational features positions the CrI₂ monolayer as a promising material for next-generation spintronic and optoelectronic applications.

Author contributions

F. A. Núñez-Murillo: data curation, formal analysis, investigation, methodology, validation, visualization, and writing – review & editing. L. A. López: investigation, methodology, validation, visualization, and writing – review & editing. D. Mejía-Burgos: investigation, methodology, validation, visualization, and writing – review & editing. J. Mazo-Zuluaga: formal analysis, investigation, methodology, validation, visualization, writing – original draft, and writing – review & editing. J. Mejía-López: conceptualization, data curation, formal analysis, investigation, methodology, project administration, software, resources, supervision, validation, visualization, writing – original draft, and writing – review & editing.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). The Supplementary Information file includes (i) CONTCAR files of the studied system, (ii) details on the van der Waals functional used, (iii) Radial and Angular Distribution functions for the O- and M- phases of the ML, (iv) details on the energy convergence, (v) comparisons between PBE and HSE06 functionals, (vi) comparisons with other 2D magnetic semiconductors, (vii) effect if spin-orbit coupling, (viii) potential strategies to improve CrI2 ML optoelectronic performance. See DOI: https://doi.org/10.1039/d5nr01077b.

Acknowledgements

The authors gratefully acknowledge support from ANID/CONICYT FONDECYT grant no. 1210193 and the Millenium Institute on Green Ammonia as Energy Vector MIGA-ANID/Millennium Science Initiative Program/ICN2021_023 (Chile). J. M.-Z. acknowledges the “Año Sabático”, CODI 2023-63352, “Estrategia de Sostenibilidad GICM 2023-ES84230019” and “Dedicación Exclusiva 2024-2025” grants at the Universidad de Antioquia (Colombia). The unit cell renders featured in the graphical abstract and figures were created with VESTA (K. Momma and F. Izumi “VESTA 3 for three-dimensional visualization of crystal, volumetric and morphology data, J. Appl. Crystallogr. 2011, 44, 1272-1276)”.

References

1. T. Alonso-Lanza, A. Ayuela and F. Aguilera-Granja, Phys. Chem. Chem. Phys., 2016, 18, 21913–21920.
2. T. Alonso-Lanza, A. Ayuela and F. Aguilera-Granja, Phys. Rev. B, 2016, 94, 245441.
3. T. Alonso-Lanza, A. Ayuela and F. Aguilera-Granja, Phys. Rev. Mater., 2017, 1, 024001.
4. T. Alonso-Lanza, A. Ayuela and F. Aguilera-Granja, Nanomaterials, 2023, 13, 3033.
5. F. Aguilera-Granja, J. Montejano-Carrizales, E. Vogel and J. Escrig, J. Phys. Chem. Solids, 2022, 164, 110646.
6. A. Lebon, J. A. Reyes-Retana and A. D. Hernández-Nieves, Int. J. Hydrogen Energy, 2019, 44, 1021–1033.
7. K. S. Novoselov, A. Mishchenko, A. Carvalho and A. H. C. Neto, Science, 2016, 353, 1234–1240.
8. G. Wang, A. Chernikov, M. M. Glazov, T. F. Heinz, X. Marie, T. Amand and B. Urbaszek, Rev. Mod. Phys., 2018, 90, 021001.
9. D. Staros, G. Hu, J. Tiihonen, R. Nanguneri, J. Krogel, M. C. Bennett, O. Heinonen, P. Ganesh and B. Rubenstein, J. Chem. Phys., 2022, 156, 014707.
10. K. Tang and W. Qi, Adv. Funct. Mater., 2020, 30, 2002672.
11. F. He, Y. Zhou, Z. Ye, S. H. Cho, J. Jeong, X. Meng and Y. Wang, ACS Nano, 2021, 15, 5944–5958.
12. D. Qiu, C. Gong, S. S. Wang, M. Zhang, C. Yang, X. Wang and J. Xiong, Adv. Mater., 2021, 33, 2006124.
13. X. Chen, J. L. Schmehr, Z. Islam, Z. Porter, E. Zoghlin, K. Finkelstein, J. P. Ruff and S. D. Wilson, Nat. Commun., 2018, 9, 103.
14. H. Huang, S. J. Lee, Y. Ikeda, T. Taniguchi, M. Takahama, C. C. Kao, M. Fujita and J. S. Lee, Phys. Rev. Lett., 2021, 126, 167001.
15. R. H. A. del Toro, T. Alonso-Lanza and A. Ayuela, Nanoscale Adv., 2023, 5, 2813–2819.
16. S. Manzeli, D. Ovchinnikov, D. Pasquier, O. V. Yazyev and A. Kis, Nat. Rev. Mater., 2017, 2, 1–15.
17. L. Zhang, J. Zhou, H. Li, L. Shen and Y. P. Feng, Appl. Phys. Rev., 2021, 8, 021308.
18. Y. Liu, S. Zhang, J. He, Z. M. Wang and Z. Liu, Nano-Micro Lett., 2019, 11, 13.
19. X. Zhou, X. Hu, J. Yu, S. Liu, Z. Shu, Q. Zhang, H. Li, Y. Ma, H. Xu and T. Zhai, Adv. Funct. Mater., 2018, 28, 1706587.
20. S. M. Hus and A.-P. Li, Prog. Surf. Sci., 2017, 92, 176–201.
21. B. Huang, G. Clark, E. Navarro-Moratalla, D. R. Klein, R. Cheng, K. L. Seyler, D. Zhong, E. Schmidgall, M. A. McGuire, D. H. Cobden, W. Yao, D. Xiao, P. Jarillo-Herrero and X. Xu, Nature, 2017, 546, 270–273.
22. C. Gong, L. Li, Z. Li, H. Ji, A. Stern, Y. Xia, T. Cao, W. Bao, C. Wang, Y. Wang, Z. Q. Qiu, R. J. Cava, S. G. Louie, J. Xia and X. Zhang, Nature, 2017, 546, 265–269.
23. X. Bo, L. Fu, X. Wan, S. Li and Y. Pu, Phys. Rev. B, 2024, 109, 014405.
24. H. Zhang, W. Yang, Y. Ning and X. Xu, Phys. Rev. B, 2020, 101, 205404.
25. Y. Hou, J. Kim and R. Wu, Sci. Adv., 2019, 5, eaaw1874.
26. S. Kezilebieke, M. N. Huda, V. Vaňo, M. Aapro, S. C. Ganguli, O. J. Silveira, S. Głodzik, A. S. Foster, T. Ojanen and P. Liljeroth, Nature, 2020, 588, 424–428.
27. A. Ayuela, J. Phys. Chem. C, 2020, 124, 2634–2643.
28. S. J. Gong, C. Gong, Y. Y. Sun, W. Y. Tong, C. G. Duan, J. H. Chu and X. Zhang, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, 8511–8516.
29. J. Zhang, J. Yang, L. Lin and J. Zhu, J. Semicond., 2020, 41, 122502.
30. L. Peng, J. Zhao, M. Cai, G. Y. Hua, Z. Y. Liu, H. N. Xia, Y. Yuan, W. H. Zhang, G. Xu, L. X. Zhao, Z. W. Zhu, T. Xiang and Y. S. Fu, Phys. Rev. Res., 2020, 2, 023264.
31. P. Li, C. Wang, J. Zhang, S. Chen, D. Guo, W. Ji and D. Zhong, Sci. Bull., 2020, 65, 1064–1071.
32. X. Cai, Z. Xu, S. H. Ji, N. Li and X. Chen, Chin. Phys. B, 2021, 30, 028102.
33. H. Liu, J. T. Sun, M. Liu and S. Meng, J. Phys. Chem. Lett., 2018, 9, 6709–6715.
34. Y. Zhao, H. Liu, J. Gao and J. Zhao, Phys. Chem. Chem. Phys., 2021, 23, 25291–25297.
35. R. O. Jones, Rev. Mod. Phys., 2015, 87, 897–923.
36. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558–561.
37. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 49, 14251–14269.
38. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–150.
39. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
40. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
41. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758.
42. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
43. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
44. J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou and K. Burke, Phys. Rev. Lett., 2008, 100, 136406.
45. J. Heyd, G. E. Scuseria and M. Ernzerhof, J. Chem. Phys., 2003, 118, 8207–8215.
46. J. Klimeš, D. R. Bowler and A. Michaelides, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 83, 195131.
47. M. Dion, H. Rydberg, E. Schröder, D. C. Langreth and B. I. Lundqvist, Phys. Rev. Lett., 2004, 92, 246401.
48. K. Lee, E. D. Murray, L. Kong, B. I. Lundqvist and D. C. Langreth, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 081101.
49. G. Roman-Perez and J. M. Soler, Phys. Rev. Lett., 2009, 103, 096102.
50. A. Tkatchenko and M. Scheffler, Phys. Rev. Lett., 2009, 102, 073005.
51. T. Bučko, S. Lebegue, J. Hafner and J. G. Angyan, J. Chem. Theory Comput., 2013, 9, 4293–4299.
52. T. Bučko, S. Lebegue, J. G. Angyan and J. Hafner, J. Chem. Phys., 2014, 141, 034114.
53. A. Tkatchenko, R. A. DiStasio, R. Car and M. Scheffler, Phys. Rev. Lett., 2012, 108, 236402.
54. S. N. Steinmann and C. Corminboeuf, J. Chem. Theory Comput., 2011, 7, 3567–3577.
55. S. N. Steinmann and C. Corminboeuf, J. Chem. Phys., 2011, 134, 044117.
56. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456–1465.
57. H. J. Monkhorst and J. D. Pack, Phys. Rev. B, 1976, 13, 5188–5192.
58. K. Parlinski, Computer Code PHONON. See: https://wolf.ifj.edu.pl/phonon/, 2008.
59. S. López-Moreno, A. H. Romero, J. Mejía-López and A. Muñoz, Phys. Chem. Chem. Phys., 2016, 18, 33250–33263.
60. A. Togo and I. Tanaka, Scr. Mater., 2015, 108, 1–5.
61. A. Fonari and S. Stauffer, vasp_raman.py, 2013, https://github.com/raman-sc/VASP/.
62. S. Zhang, F. Tang, X. Song and X. Zhang, Phys. Rev. B, 2022, 105, 104105.
63. F. Besrest and S. Jaulmes, Acta Crystallogr., Sect. B, 1973, 29, 1560–1563.
64. P. Li, N. Liu, J. Zhang, S. Chen, X. Zhou, D. Guo, C. Wang, W. Ji and D. Zhong, ACS Appl. Mater. Interfaces, 2023, 15, 19574–19581.
65. F. A. Cotton, Chemical applications of group theory, John Wiley & Sons, 1991.
66. M. S. Dresselhaus, G. Dresselhaus and A. Jorio, Group theory: Application to the physics of condensed matter, Springer Science & Business Media, 2007.
67. R. Merlin, A. Pinczuk and W. Weber, Raman scattering in materials science, Springer, 2000, pp. 1–29.
68. Y. Peter and M. Cardona, Fundamentals of semiconductors: Physics and materials properties, Springer Science & Business Media, 2010.
69. J.-W. Jiang, B.-S. Wang and H. S. Park, Nanoscale, 2018, 10, 13913–13923.

---