Full quantum search for high Tc two-dimensional van der Waals ferromagnetic semiconductors

Liang Liu ab, Zezhou Lin a, Jifan Hu c and Xi Zhang *a
aInstitute of Nanosurface Science and Engineering, Guangdong Provincial Key Laboratory of Micro/Nano Optomechatronics Engineering, Shenzhen University, Shenzhen 518060, China. E-mail: zh0005xi@szu.edu.cn
bKey Laboratory of Optoelectronic Devices and Systems of Ministry of Education and Guangdong Province, College of Optoelectronic Engineering, Shenzhen University, Shenzhen, 518060, China
cSchool of Physics, State key laboratory for crystal materials, Shandong University, Jinan 250100, China

Received 8th December 2020 , Accepted 5th March 2021

First published on 30th March 2021


Atomic thin two-dimensional (2D) ferromagnetic (FM) semiconductors with high Curie temperatures (Tc) are essential for future spintronic applications. However, reliable theoretical searching for 2D FM semiconductors is still hard due to the complexity of strong quantum fluctuations in 2D systems. We have proposed a full quantum search (FulQuanS) method to tackle the difficulty, and finally identified five 2D semiconductors of CrX3 (X = I, Br, Cl), CuCl3 and FeCl2 with FM order at finite temperature from the pool of 3721 potential 2D structures. Via the method of renormalized spin wave theory (SW) and quantum Monte Carlo simulations (QMC), we located the Tc for CrX3 (X = I, Br, Cl), CuCl3 and FeCl2 at 48 K, 31 K, 18 K, 74 K and 931 K respectively, which excellently agree with experiments for CrX3 and reveal the superior performances of the new predicted structures. Furthermore, our QMC results demonstrated that the systems with low-spin numbers and/or low anisotropies have much higher Tc than the estimations of classical models e.g., Monte Carlo simulations based on classical Heisenberg models. Our findings suggest excellent candidates for future room-temperature spintronics, and shed light on the quantum effects inherent in 2D magnetism.


2D van der Waals (vdW) layers are believed to be the cornerstones for future electronic devices, since the quasi-free-standing structures can protect them from the complexities of substrates and also make them easy to integrate.1–3 Recently, intrinsic magnetism was discovered in 2D vdW layered Cr2Ge2Te6 and CrI3 in which ferromagnetic orders survive even in a single quasi-free-standing layer, opening exciting prospects for 2D vdW materials.4,5 Particularly, the semiconducting nature of CrI3 helps it realize a series of fantastic functionalities e.g., electric field control of magnetism, magnetic topological phase, and valleytronics.6–10 However, except few vdW metallic layers such as Fe3GeTe2, CrTe2, MnSe2, VTe2 and VSe2, all experimentally observed 2D magnetic orders in semiconductors persist only at ultra-low temperature which is still far from practical applications.11–17 Hence, searching for 2D vdW magnetic semiconductors with high Tc is desperately demanded.

Recently, theoretical searching for 2D ferromagnetism has drawn intense attention, and significant progress has been made, providing essential references for the realization of superior 2D ferromagnets in laboratories.18–23 To date, researchers have proposed more than six hundred potential 2D ferromagnets via comparing the energies of ferromagnetic and some antiferromagnetic (AFM) configurations at the density functional theory (DFT) level.24,25 However, it is still highly difficult to provide accurate descriptions and predictions for 2D ferromagnets. Firstly, the complexities of AFM configurations e.g., C-type, G-type, A-type AFM and some long-range AFM configurations make it impossible to enumerate all possible AFM configurations and ensure the FM ground states within conventional methods. On the other hand, the Curie temperature which is one of the key properties in 2D magnets is very nontrivial to estimate via usual theoretical methods. In the literature, the locations of Tc were accomplished by the mean-field method, spin wave (SW) theory26 and Monte Carlo simulations (MC) of Ising model27 or classical Heisenberg model.19,28 The mean-field method is a linear mapping from exchange interactions J to Tc. It did not count the effects of magnetic anisotropy which is crucial in the 2D limits, failing to capture the phase transitions in 2D magnetism.19 So did the Ising model. SW and classical Heisenberg MC involve the considerations of magnetic anisotropies, but treat spins as classical sticks in which the approximation integrities were based on the large spin-number limit. To date, there is a surprising paucity of studies that seek to identify the effect of the quantum effect on Tc of 2D magnets which is supposed to be crucial in the case of a small spin-number such as in the famous CrX3 family. All these signify the imperatives of new and plausible methods that capture the ground state properties i.e. the exchange interactions and anisotropies of 2D FM systems, as well as going deep into the excited states near phase transitions with full considerations of quantum effects, predicting new 2D FM semiconductors with high efficiency and credibility.

In this work, we organized comprehensive investigations to find a series of 2D FM semiconductors with relatively high Tc. We proposed high throughput full quantum search (FulQuanS) to discover promising 2D ferromagnetic semiconductors and locate their Tc. The magnetic properties and effective parameters in the Heisenberg model were extracted from relativistic density functional theory (DFT) and we finally screened five candidates from 3712 structures in a computational 2D material database (C2DB), including the well-known CrX3 (X = Cl, Br, I), and the other two new structures of CuCl3 and FeCl2. We employed renormalized SW and a quantum Monte Carlo (QMC) method to locate the phase transitions. For CrX3 (X = Cl, Br, I), CuCl3 and FeCl2, Tc values predicted by FulQuanS were 48 K, 31 K, 18 K, 74 K and 931 K respectively, showing excellent agreement with experiments in the CrX3 family. More crucially, the results obtained using FulQuanS clarified the key roles of quantum effects in the phase transitions of 2D magnetism, including the frustrations in low spin systems with trigonal lattices, and the stabilizations of magnetic orders in systems with low anisotropies and low spin numbers. All these findings suggest promising candidates for future 2D magnetism study and spintronic application as well as clarifying the implications of quantum effects inherent in 2D magnetism.



Fig. 1(a) describes the main workflow of high-throughput FulQuanS, including three stages: (1) data mining and high-throughput calculations for searching the potential structures, (2) model mapping based on relativistic DFT calculations to obtain the anisotropic Heisenberg Hamiltonian, and (3) renormalized SW and QMC calculations to compute Tc and clarify the role of quantum effects inherent in 2D magnetic systems.
image file: d0nr08687h-f1.tif
Fig. 1 (a) Flowchart of FulQuanS utilized in this work for high throughput searching for 2D FM semiconductors and Tc estimations. Numbers in circles denote the quantity of input structures for each step while negative numbers in ellipse denote the quantity of ruled out structures. (b) Band gaps and energies above hull of every structure considered in this work. Grey, green, orange, and red dots represent the 772 stable 2D systems with band gap, 454 structures with open d/f shell, 128 structures with ferromagnetic ground state (verified in small supercell) and 8 ferromagnetic semiconductors with out-of-plane anisotropies.

Stage 1 is the data mining and high-throughput search for 2D FM semiconductors. In this stage we screened out all the potential high-Tc FM semiconductors by means of high throughput searching. We examined totally 3721 2D structures from C2DB,24 which is one of the largest 2D material databases. The considerations of stabilities (with formation energy above hull less than 0.2 eV) and band gaps (>0) screened 772 structures out of whole C2DB. In Fig. 1(b) we summarized their band gaps and formation energy above hull. Since all of them have relatively low formation energy, these structures thus are promising to be easily synthesized via bottom-up chemical methods. And the finite band gaps reveal the semiconducting nature, ensuring the integrities of utilizing Heisenberg models to capture the magnetic interactions in these structures.

Afterwards, we screened out structures with open d- or f-shells, consisting of 454 structures and denoted by orange dots in Fig. 1(b). The magnetism of s- and p-orbitals rarely exists in perfect crystals and nor does that of those with only closed shells, they are thus safely skipped. Then, the remaining structures were optimized and their ground electronic structures were computed with colinear-DFT. It was shown that only 128 of them converged to ferromagnetic ground states (orange dots in Fig. 1(b)).

Finally, we studied the magnetic anisotropies of the remaining structures based on noncollinear-DFT to screen out the FM structures with out-of-plane anisotropies, resulting in eight structures as shown by the red dots in Fig. 1(b). According to Mermin–Wagner theorem, breaking the continual symmetries via anisotropy is crucial for a 2D system to sustain magnetic order at finite temperature.29 Furthermore, most of the experimentally confirmed 2D magnetic semiconductors are with out-of-plane preferred magnetic anisotropy. Therefore, we only leave those with out-of-plane anisotropies for further considerations.

Stage 2 is the construction of the anisotropic Heisenberg model:

image file: d0nr08687h-t1.tif(1)
in which Szi is the z-component of the local spin operator acting on the ith atom. S±i = Sxi ± iSyi are the spin raising and lowering operators, respectively. Jij is the isotropic exchange interactions between ith and jth magnetic atoms. Bij is the exchange anisotropy and Az is the single-ion anisotropy. E0 is the energy unrelated to magnetism. The summation is over all inequal i, j pairs. All the magnetic parameters (Jij, Bij, Az) were deduced from the energy mapping method with a series of relativistic DFT calculations (see the ESI for more details).

Stage 3 is to solve the anisotropic Heisenberg model eqn (1) for the determinations of Tc based on Hatree–Fock renormalized spin wave theory (HF-SW) and path-integral QMC. In HF-SW, the low-energy behaviors of Heisenberg model were captured with boson (each spin wave is a boson) systems, via operator substitutions known as Dyson–Maleev transformation:30

image file: d0nr08687h-t2.tif(2)
where bi(bi) creates (annihilates) one spin wave boson on the ith magnetic ion, and ni = bibi counts the number of spin wave excitons on the ith ion. S is the maximally possible eigenvalue of Szi. Hence, the Heisenberg Hamiltonian in eqn (1) can be transformed into the Hamiltonian expression only consisting of bi and bi operators, which can be solved by the Hartree–Fock renormalization method (see the ESI for more details). The solutions of the Hamiltonian are the spin waves (eigenvector) and spin wave energy dispersions (eigenvalue). The ground state (E = 0, with ni = 0 for all sites) is exactly the ferromagnetic state of the Heisenberg model (with Sz = −S for all sites). The lowest SW eigenenergy is defined as spin-wave gap ΔSW of the HF-SW energy spectra. ΔSW qualifies the kinetic stability of ferromagnetism. ΔSW > 0 indicates that certain energy is required for the system to transmit into excited states (SW eigenstate), and the 2D ferromagnetism (ground state) is protected from small thermal fluctuations.

It may be noticed that the exact value of ΔSW depends on the occupations and the temperature once the interactions were taken into consideration. Fortunately, the property of ΔSW > 0 always remains in stable ferromagnetic systems until the phase transition (see the ESI for more discussions on this point). Therefore, utilizing the 0 K ΔSW as the criteria for screening is sufficient.

Besides, the averaged magnetic moment per site at finite temperature T is related to the quantity of spin wave excitons: M(T) = 2(Sn(T)), in which n(T) satisfies the standard Boson–Einstein distribution. TSWc is obtained at the temperature where magnetism vanishes, i.e. M(Tc) = 0.

The HF-SW method is analytical and costless, providing good descriptions for the large-scale magnetic configurations in the ground state and coarse description for the magnetic states near Tc. Meanwhile QMC is supposed to be numerically accurate near Tc but expensive. Therefore, FulQuanS was designed to integrate the advantages of both. We utilized HF-SW to have a coarse estimation for Tc, which is treated as the initial guess for QMC. Afterwards, we performed QMC31,32 on 16 × 16 and 32 × 32 supercells, using the scaling behavior of 4th order Binder cumulate B(T) = m22/m4 to extrapolate the numerically accurate Tc of infinitely sized systems.31,32 To investigate how the quantum effects impact phase transition and thermal dynamics of 2D magnetism, we also conducted the classical Monte Carlo (CMC)33,34 simulation in this stage for comparisons.

DFT settings

We performed the total-energy calculations via DFT with local density approximation (LDA) exchange–correlation functionals and projected augmented plane-wave (PAW) basis which is implemented in the Vienna Ab initio Simulation Package (VASP).35,36 The cutoff energy for the planewave basis is set to 500 eV. K-point meshes are sampled to make the actual spacing less than 0.02 Å−1 in each periodic direction. All the DFT calculations are based on optimized structures with residual Hellman–Feynman force <0.001 eV Å−1. We also employed the Hubbard U method37 (LDA + U) to describe the d-state interactions. As shown in Tables S1 and S2 in the ESI, up to the 0 eV ≤ Ueff ≤ 2 eV applied on d-states, the Tc order amongst all candidates did not change and all the conclusions are verified to be held well. To check the conclusions on levels beyond LDA functionals, we also verified the results based on the functionals of general gradient approximation (GGA) with a Perdew–Burke–Ernzerhof (PBE) form.36,38 As shown in the ESI, the distinctions between LDA and GGA are only numerical and all main conclusions were preserved. Also, the calculated Tc for experimentally known systems based on GGA and hybrid functionals39 seems to be more overestimated than LDA, we thus mainly present the LDA results. Finally, the dynamical and thermal stabilities of predicted candidates were also justified by the phonon vibrations with the method of linear response.40 And the formation energy diagrams were constructed with initial 2D structures extracted from C2DB (see the ESI for more details).

Results and discussion

High throughput screening

Fig. 2 shows the eight 2D FM semiconductor structures screened by stage 1, including three experimentally confirmed CrX3 (X = I, Br, Cl) and five new structures vis. CuCl3, FeCl2, TiX3 (X = Br, Cl) and Fe2Br2O2. Amongst them, CuCl3 has the same prototype structure as the CrX3 family, in which the magnetic cations comprise hexagonal magnetic crystals and every cation is surrounded by six anions which impose distorted octahedral crystal field on the magnetic ions. FeCl2 shares a similar structure with 2H-MoS2, and Fe2+ ions comprise triangular magnetic crystals and they are capsulated in triangular-prism crystal fields. The structures of the TiX3 family (X = Br, Cl) can be viewed as locally deformed CrX3 structures, and their magnetic ions also comprise hexagonal lattices but endure a triangular prism crystal field. Fe2O2Br2 has an orthorhombic lattice and Fe3+ ions endure distorted and unbalanced octahedral crystal fields of O2− and Br. Either an octahedral or triangular-prism local crystal field has the characteristics that the nearest magnetic cations are separated by the anions with nearly 90-degree angles which indicate that the ferromagnetic super-exchanges dominate the couplings between the nearest neighbored (1NN) magnetic ions. One can refer to the ESI for more detailed analyses of their electronic structures and magnetic structure.
image file: d0nr08687h-f2.tif
Fig. 2 Top and side views of eight 2D FM semiconductors with short-range ferromagnetism and out-of-plane anisotropies including three with experimental confirmations and five new predicted by DFT.

Anisotropic Heisenberg parameters obtained by DFT

In Table 1, we have listed the Heisenberg Hamiltonian parameters including net spin number S, exchange interaction J, exchange anisotropy B and single-ion anisotropy Az deduced from noncolinear DFT calculations for all eight structures. Our calculations show that the net spin number of Cr3+ is 1.5, indicating that the d-electrons in Cr3+ occupy a single spin-tunnel of t2g orbitals, in accordance with other studies.27,41 And the net spin of Cu3+ and Ti3+ is low with S = 1 and 0.5 respectively; such a low spin number may lead to significant quantum fluctuations which we will discuss later. Fe2+ ions in FeCl2 show a high spin state with S = 2. The Fe3+ ion in Fe2O2Br2 shows a low spin state with S = 1/2, which may be caused by the weak exchange splitting in Fe2O2Br2.
Table 1 Magnetic parameters for the candidates from noncolinear DFT calculations. S is the maximally possible net eigenvalues of Sz. J1J3 are the isotropic part of 1–3 NN exchange interactions. B are the anisotropic parts of exchange interactions. Az are the single-ion magnetic anisotropies. All energy units are meV
  S J1 (B1) J2 (B2) J3 (B3) A z
CrI3 1.5 −1.69 (−0.08) −0.51 (0.02) 0.46 (0.01) −0.272
CrBr3 1.5 −1.41 (−0.03) −0.34 (0.01) 0.44 (0.00) −0.08
CrCl3 1.5 −1.00 (−0.01) −0.18 (0.00) 0.29 (0.00) −0.031
CuCl3 1.0 −9.02 (−0.01) −1.84 (0.00) 3.20 (0.00) 0.01
TiBr3 0.5 −120.66 (−0.02) 39.43 (−0.03) −26.52 (0.00) 0.014
TiCl3 0.5 −114.60 (−0.02) 42.71 (−0.02) −29.96 (−0.01) −0.011
FeCl2 2.0 −15.40 (−0.02) −0.45 (−0.01) 1.52 (0.03) −0.049
Fe2O2Br2 0.5 −28.78 (−0.19) −8.88 (−0.63) −24.89 (0.36) −2.47

The calculated leading exchange couplings J1 in these nine structures are negative, indicating that the ferromagnetic order is favoured at least in the 1NN range, in agreement with the super-exchange mechanism discussed before. The 2NN exchange couplings J2 are much smaller than J1. And J2 values in all structures are also of the ferromagnetic type except TiCl3 and TiBr3. Most importantly, except FeCl2, exchange couplings of our candidates have non-negligible strengths up to the 3NN range, revealing the long-range nature of super-exchange. For instance, J3 values are shown to be 0.46 meV, 0.44 meV and 0.29 meV for CrI3, CrBr3 and CrCl3 respectively, nearly one third the amount of J1.

The strongest single-ion anisotropy was found in Fe2O2Br2 with Az = −2.47 meV, about nine times larger than that of CrI3. This stems from the symmetry-breaking in local octahedra crystal field caused by the large electrostatic differences between oxygen and bromine anions.

T c determined by SW and QMC

Fig. 3(a)–(h) show the spin-wave energy dispersion spectra for all eight structures. The spin-wave gap ΔSW in the HF-SW spectra was also marked. As we discussed before, ΔSW quantifies the kinetic stability of ferromagnetism. ΔSW > 0 indicates that 2D ferromagnetism (ground state) is stable and certain energy is required for the system to transmit into other magnetic states. Small negative eigenvalues are presented around the Γ point for TiBr3 and TiCl3, indicating that the FM configurations are not the most stable states for them and would spontaneously evolve to long-range AFM configurations at 0 K. Therefore, these two structures were ruled out, leaving the other six structures for further investigations.
image file: d0nr08687h-f3.tif
Fig. 3 (a)–(h) Spin-wave spectra of eight structures. The negative energies in the spectra of (g) TiBr3 and (h) TiCl3 are marked by arrows. (i) Illustration of the 2D projection of a particular spin wave eigenstate of TiBr3 with negative eigenvalues. Atoms at the left-bottom corner represent the lattice and the arrows represents the spin vectors of Cr3+ ions, and colours of arrows denote the phases of spins.

One of the particular HF-SW eigenvector of TiBr3 with negative eigenvalue and wave vector k = (2π/15a, 2π/15a) was illustrated in Fig. 3(i), in which part of the lattice (top view) was shown in the left-down corner. The arrows denote spins of Ti3+ ions and the colours represent the phases of spins. This kind of AFM configuration has energy lower than the FM states and may be missed by conventional considerations within small supercells due to its long-range nature (with a period of 15 times lattice constant).

Due to the out-of-plane magnetic anisotropies, finite energy gaps ΔSW at the Γ point are presented for the remaining six spectra. These energy gaps ΔSW are dependent on the temperature, and the evolutions of ΔSW are governed by the anisotropic properties including the exchange and single-ion parts. However, until the phase transition, all ΔSW values are preserved as positive to protect the magnetic orders (see the ESI for more details).

To describe the whole anisotropy quantitatively, including both effects from anisotropies and exchange interactions, we define a dimensionless quantity D = 103ΔSW/Em, where Em is the largest eigenvalue in SW spectra, reflecting the strengths of exchange interactions inherent in the system. In the CrX3 family, the largest anisotropy with D = 49.92 was found in CrI3 while the anisotropies in CrBr3 and CrCl3 are weaker with D = 22.66 and 12.22 respectively, in good agreement with experiments.5,13,14 The anisotropies of CuCl3, FeCl2 and Fe2O2Br2 were found to be 0.2387, 0.8565 and 24.29, respectively. Therefore, the 2D magnetism in Fe2O2Br2 has anisotropy comparative to CrBr3 while the other 2 systems show nearly isotropic behaviours.

Fig. 4 shows the net magnetic moment evolutions at finite temperatures for the six structures. To figure out the effects of quantum nature, both the results from classical Monte Carlo simulations and QMC were collected here. In the low temperature region, the quantum fluctuations have heavily suppressed magnetism in the low-spin system and even lead to the frustration in Fe2O2Br2. More importantly, quantum effects also play an essential role in high temperature regions and decide the actual positions of phase transitions.

image file: d0nr08687h-f4.tif
Fig. 4 Mean magnetic moment (per ion) vs. temperature, given by CMC (blue) and QMC (red) for each structure. Green stars correspond to experimentally observed Tc.

The exact Tc inherent in quantum models and classical models were obtained by scaling methods and marked in Fig. 4. Detailed scaling behaviours of these systems are shown in Fig. S1–S5 in the ESI. The experimental Tc of the CrX3 family (45 K, 28 K and 17 K for CrI3, CrBr3 and CrCl3 as reported in ref. 5, 13 and 14) were marked with green stars for comparison. The TQMCc predicted by QMC for CrI3, CrBr3 and CrCl3 are 48 K, 31 K and 18 K, respectively, in excellent agreement with experiments, demonstrating the high accuracy of FulQuamS. However TCMCc values predicted by classical Monte Carlo simulations for CrI3, CrBr3 and CrCl3 are 33 K, 20 K and 11 K, respectively, significantly below experimental observations, revealing the deficiencies of classical methods for the descriptions of 2D magnetism, also indicating the essential roles of quantum effects in the persisting of 2D magnetic orders near Tc.

Notably, the most robust ferromagnetism was found in FeCl2 with Tc = 931 K. Such high Tc is worthwhile for experimental verifications in future works, and the theoretical reasons for this result are two fold: (i) As shown in Table 1, 1NN ferromagnetic exchange couplings of FeCl2 are ∼60 meV (equals |J1 × S2|), almost 10 times larger than that in CrBr3 and (ii) the magnetic lattice of FeCl2 is trigonal, supporting six 1NN ferromagnetic couplings for each magnetic moment, which double the case of CrBr3. Therefore, it is straight to understand why FeCl2 has Tc of ∼30 times higher than CrBr3. Amongst the remaining candidates, Fe2O2Br2 was found to be frustrated in QMC due to the spin-1/2 nature and triangular lattice, and the Tc of new predicted structures CuCl3 turned out to be 74 K, well beyond the Tc in CrI3, worth future experimental investigations, too.

Role of quantum effects

To further quantify the quantum effects on Tc, we define the underestimation rate UR = (TQMCcTCMCc)/TQMCc. And URs for all six structures are shown in Fig. 5(a) with corresponding anisotropies D defined before. We found that URs are always positive which means that Monte Carlo simulations based on classical Heisenberg models always underestimate Tc. And it is apparent that the UR is inversely related to the out-of-plane anisotropies. In CrI3, which hosts the strongest anisotropies, the differences between classic and quantum are the smallest with UR = 0.31. However, in the case of CuCl3, which hosts nearly isotropic 2D ferromagnetism, the impact of quantum effects is huge with UR = 0.55.
image file: d0nr08687h-f5.tif
Fig. 5 (a) The underestimation rate of Tcvia classical models for the five candidates, relating to their anisotropies. Dashed lines are guide to eyes. (b) Temperature evolutions of spin structure factors for original CrCl3 (S = 3/2, blue), CrCl3 with reduced spin number (S = 1/2, red), CrCl3 with enlarged spin number (S = 9/2, yellow) and classical CrCl3 (S = , green). (c) Temperature evolutions of spin structure factors for original FeCl2 (S = 2, blue), FeCl2 with reduced spin number (S = 3/2, red), FeCl2 with enlarged spin number (S = 6, yellow), and classical FeCl2 (S = , green).

Another crucial factor related to quantum effects is the spin number. To figure out this, we performed QMC on a series of systems in which all the magnetic related energies except the spin numbers are the same. Fig. 5(b) shows the averaged spin structure factor (thermally averaged net spin divided by spin number S) of honeycomb ferromagnetic systems which share the same SW spectra with CrCl3 (their behaviours in CMC are the same, too) but different spin numbers of 1/2, 3/2 and 9/2, corresponding to the reduced spin, original and enlarged spin models. At temperatures much lower than Tc, the spin structure factor of the spin-1/2 system is slightly smaller than that of the spin-3/2 system due to the suppressions from quantum fluctuations embedded in spin-1/2. More importantly, systems with lower spin numbers always reveal higher Tc. As a result, at high temperatures, the spin structure factors of spin-1/2, 3/2, 9/2 and infinite (corresponding to classical case) systems are in descending order.

In triangular ferromagnetic systems which share the same SW spectra with FeCl2, we found similar behaviours, as shown in Fig. 5(c). In the low temperature region, the averaged spin structure factors are positively related to the spin numbers due to the strong quantum fluctuations. On the other hand, at high temperatures where phase transition occurs, the spin structure factors become negatively related to the temperature. Hence, the quantum effects from a low spin number also play important roles in the phase transition and always lead to the enhancements of Tc in trigonal systems, too.


In summary, we proposed a high throughput process FulQuanS to screen 2D ferromagnetic semiconductors from a large database with high reliabilities and predict Tc with high accuracies. In total, we have found five candidates including CrX3 (X = I, Br, Cl), CuCl3 and FeCl2 which are supposed to host stable ferromagnetism at finite temperatures, and their Tc values are located at 48 K, 31 K, 18 K, 74 K and 931 K respectively, with full considerations of quantum effects. The results for the CrX3 family agree with experiments excellently, indicating the superior accuracies of FulQuanS. And the two new predicted candidates all have Tc beyond CrI3, especially that the Tc in FeCl2 is far beyond room-temperature, revealing promising applications in future spintronics. Furthermore, we found that the quantum effects play essential roles in 2D magnetic systems for not only the low temperature region but also the high temperature region near phase transition, and classical models either falsely predict the magnetic orders or underestimate Tc for systems with low spin numbers and/or low anisotropies. All these provide us an excellent process for future studies on 2D magnetism, and enriched our understanding in the intrinsic quantum nature of 2D systems.

Author contributions

Liang Liu developed high throughput codes and analysed the data. Zezhou Lin checked and verified all results and constructed diagrams. Jifan Hu analysed the data. Xi Zhang proposed the theme, analysed the results, and supervised the whole work. All authors contributed to the manuscript writing.

Conflicts of interest

There are no conflicts to declare.


The authors thank Dr Xiaolin Li for the help in editing figures. This work was supported by the National Natural Science Foundation (No. 51605306) of China, and Shenzhen Overseas High-level Talents Innovation and Entrepreneurship Plan (No. KQJSCX20180328094853770).

Notes and references

  1. C. Gong and X. Zhang, Science, 2019, 363, eaav4450 CrossRef CAS PubMed.
  2. K. S. Burch, D. Mandrus and J.-G. Park, Nature, 2018, 563, 47–52 CrossRef CAS PubMed.
  3. X. Li, B. Dong, X. Sun, H. Wang, T. Yang, G. Yu and Z. V. Han, J. Semicond., 2019, 40, 081508 CrossRef CAS.
  4. 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 CrossRef CAS PubMed.
  5. 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 CrossRef CAS PubMed.
  6. S. Jiang, J. Shan and K. F. Mak, Nat. Mater., 2018, 17, 406–410 CrossRef CAS PubMed.
  7. N. Samarth, Nature, 2017, 546, 216 CrossRef CAS PubMed.
  8. Nat. Nanotechnol., 2018, 13, 269 Search PubMed.
  9. D. Ghazaryan, M. T. Greenaway, Z. Wang, V. H. Guarochicomoreira, I. J. Veramarun, J. Yin, Y. Liao, S. V. Morozov, O. Kristanovski and A. I. Lichtenstein, arXiv: Mesoscale and Nanoscale Physics, 2018, 1, 344–349.
  10. O. Johansen, V. Risinggard, A. Sudbo, J. Linder and A. Brataas, Phys. Rev. Lett., 2019, 122, 217203 CrossRef CAS PubMed.
  11. M. Bonilla, S. Kolekar, Y. Ma, H. C. Diaz, V. Kalappattil, R. Das, T. Eggers, H. R. Gutierrez, M. H. Phan and M. Batzill, Nat. Nanotechnol., 2018, 13, 289–293 CrossRef CAS PubMed.
  12. Y. Deng, Y. Yu, Y. Song, J. Zhang, N. Z. Wang, Z. Sun, Y. Yi, Y. Z. Wu, S. Wu, J. Zhu, J. Wang, X. H. Chen and Y. Zhang, Nature, 2018, 563, 94–99 CrossRef CAS PubMed.
  13. X. Cai, T. Song, N. P. Wilson, G. Clark, M. He, X. Zhang, T. Taniguchi, K. Watanabe, W. Yao, D. Xiao, M. A. McGuire, D. H. Cobden and X. Xu, Nano Lett., 2019, 19, 3993–3998 CrossRef CAS PubMed.
  14. Z. Zhang, J. Shang, C. Jiang, A. Rasmita, W. Gao and T. Yu, Nano Lett., 2019, 19, 3138–3142 CrossRef CAS PubMed.
  15. A. O. Fumega, M. Gobbi, P. Dreher, W. Wan, C. González-Orellana, M. Peña-Díaz, C. Rogero, J. Herrero-Martín, P. Gargiani, M. Ilyn, M. M. Ugeda, V. Pardo and S. Blanco-Canosa, J. Phys. Chem. C, 2019, 123, 27802–27810 CrossRef CAS.
  16. A. O. Fumega, J. Phillips and V. Pardo, J. Phys. Chem. C, 2020, 124, 21047–21053 CrossRef.
  17. J. Li, B. Zhao, P. Chen, R. Wu, B. Li, Q. Xia, G. Guo, J. Luo, K. Zang and Z. Zhang, Adv. Mater., 2018, 30, 1801043 CrossRef PubMed.
  18. A. Kabiraj, M. Kumar and S. Mahapatra, npj Comput. Mater., 2020, 6, 1–9 CrossRef.
  19. D. Torelli, K. S. Thygesen and T. Olsen, 2D Mater., 2019, 6, 045018 CrossRef CAS.
  20. C. Huang, J. Feng, F. Wu, D. Ahmed, B. Huang, H. Xiang, K. Deng and E. Kan, J. Am. Chem. Soc., 2018, 140, 11519–11525 CrossRef CAS PubMed.
  21. S. Chen, F. Wu, Q. Li, H. Sun, J. Ding, C. Huang and E. Kan, Nanoscale, 2020, 12, 15670–15676 RSC.
  22. S. Zheng, C. Huang, T. Yu, M. Xu, S. Zhang, H. Xu, Y. Liu, E. Kan, Y. Wang and G. Yang, J. Phys. Chem. Lett., 2019, 10, 2733–2738 CrossRef CAS PubMed.
  23. C. Huang, J. Feng, J. Zhou, H. Xiang, K. Deng and E. Kan, J. Am. Chem. Soc., 2019, 141, 12413–12418 CrossRef CAS PubMed.
  24. S. Haastrup, M. Strange, M. Pandey, T. Deilmann, P. S. Schmidt, N. F. Hinsche, M. N. Gjerding, D. Torelli, P. M. Larsen, A. C. Riis-Jensen, J. Gath, K. W. Jacobsen, J. J. Mortensen, T. Olsen and K. S. Thygesen, 2D Mater., 2018, 5, 042002 CrossRef CAS.
  25. N. Mounet, M. Gibertini, P. Schwaller, D. Campi, A. Merkys, A. Marrazzo, T. Sohier, I. E. Castelli, A. Cepellotti, G. Pizzi and N. Marzari, Nat. Nanotechnol., 2018, 13, 246–252 CrossRef CAS PubMed.
  26. D. Torelli and T. Olsen, 2D Mater., 2018, 6, 015028 CrossRef.
  27. J. Liu, Q. Sun, Y. Kawazoe and P. Jena, Phys. Chem. Chem. Phys., 2016, 18, 8777–8784 RSC.
  28. S. Lu, Q. Zhou, Y. Guo, Y. Zhang, Y. Wu and J. Wang, Adv. Mater., 2020, 32, e2002658 CrossRef PubMed.
  29. N. D. Mermin and H. Wagner, Phys. Rev. Lett., 1966, 17, 1133 CrossRef CAS.
  30. F. J. Dyson, Phys. Rev., 1956, 102, 1217 CrossRef.
  31. B. Bauer, L. D. Carr, H. G. Evertz, A. Feiguin, J. Freire, S. Fuchs, L. Gamper, J. Gukelberger, E. Gull, S. Guertler, A. Hehn, R. Igarashi, S. V. Isakov, D. Koop, P. N. Ma, P. Mates, H. Matsuo, O. Parcollet, G. Pawłowski, J. D. Picon, L. Pollet, E. Santos, V. W. Scarola, U. Schollwöck, C. Silva, B. Surer, S. Todo, S. Trebst, M. Troyer, M. L. Wall, P. Werner and S. Wessel, J. Stat. Mech.: Theory Exp., 2011, 2011, P05001 Search PubMed.
  32. A. F. Albuquerque, F. Alet, P. Corboz, P. Dayal, A. Feiguin, S. Fuchs, L. Gamper, E. Gull, S. Gürtler, A. Honecker, R. Igarashi, M. Körner, A. Kozhevnikov, A. Läuchli, S. R. Manmana, M. Matsumoto, I. P. McCulloch, F. Michel, R. M. Noack, G. Pawłowski, L. Pollet, T. Pruschke, U. Schollwöck, S. Todo, S. Trebst, M. Troyer, P. Werner and S. Wessel, J. Magn. Magn. Mater., 2007, 310, 1187–1193 CrossRef CAS.
  33. L. Liu, S. Chen, Z. Lin and X. Zhang, J. Phys. Chem. Lett., 2020, 11, 7893–7900 CrossRef CAS PubMed.
  34. L. Liu, X. Ren, J. Xie, B. Cheng, W. Liu, T. An, H. Qin and J. Hu, Appl. Surf. Sci., 2019, 480, 300–307 CrossRef CAS.
  35. G. Kresse and J. Furthmuller, Comput. Mater. Sci., 1996, 6, 15–50 CrossRef CAS.
  36. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868 CrossRef CAS PubMed.
  37. V. I. Anisimov, J. Zaanen and O. K. Andersen, Phys. Rev. B: Condens. Matter Mater. Phys., 1991, 44, 943–954 CrossRef CAS PubMed.
  38. J. P. Perdew and Y. Wang, Phys. Rev. B: Condens. Matter Mater. Phys., 1992, 45, 13244 CrossRef PubMed.
  39. J. Heyd, G. E. Scuseria and M. Ernzerhof, J. Chem. Phys., 2003, 118, 8207 CrossRef CAS.
  40. X. Wu, D. Vanderbilt and D. R. Hamann, Phys. Rev. B: Condens. Matter Mater. Phys., 2005, 72, 035105 CrossRef.
  41. J. L. Lado and J. Fernández-Rossier, 2D Mater., 2017, 4, 035002 CrossRef.


Electronic supplementary information (ESI) available: The detailed energy mapping method to extract Heisenberg parameters from DFT calculations, the spin wave theory with interactions solved by Hartree–Fock approximations, the electronic and magnetic structures of five candidates and the impacts of Hubbard U on the d-shell, and verifications on the GGA-level were presented. See DOI: 10.1039/d0nr08687h

This journal is © The Royal Society of Chemistry 2021