Introduction

Long-range magnetic order is suppressed in two-dimensional (2D) isotropic systems at finite temperatures due to thermal fluctuations, according to the Mermin-Wagner theorem1. Therefore, recent discoveries of ferromagnetism in 2D materials, e.g., CrI32,3, Cr2Ge2Te64, Fe3GeTe25,6, and Fe4GeTe27 have attracted broad attention. The 2D ferromagnetic (2DFM) materials possess some fascinating properties in many aspects6,8,9, which also have great potential in device applications, such as 2D spintronics10,11,12,13,14.

However, only very few 2DFM materials have been experimentally synthesized so far, and the Curie temperatures (TCs) of these 2DFM materials are fairly low. The searching for 2DFM materials with high TC via first-principles calculations attracts more and more attention. Mounet et al. applied density functional theory (DFT) calculations to search for easily exfoliable magnetic compounds15, and they revealed a wealth of 2D magnetic systems including 37 ferromagnets. Zhu et al. found 15 2DFM materials16. Especially, Cr3Te4 was predicted to have an extremely high TC = 2057 K. Torelli et al. discovered 85 2DFM materials by high-throughput computation17. Kabiraj et al. found 26 2DFM materials with TC higher than 400 K from high-throughput scanning of 786 materials18.

In this work, we carry out a comprehensive high-throughput computational study to calculate the TCs of 2DFM materials. The magnetic exchange interactions are calculated via a first-principles linear response theory19,20,21,22. We simulate the magnetic phase diagrams and evaluate the TC by a replica-exchange Monte Carlo simulation23. We have identified 79 2DFM materials. We benchmark the calculated TCs with those of experimental synthesized materials, and also with the corresponding bulk materials. This suggests that the results are highly reliable. The compound with the highest TC we predicted is Co2F2, which has TC = 541 K. The 79 2DFM materials are classified into different structural prototypes according to their structural similarity. We perform sure independence screening and sparsifying operator (SISSO) analysis24,25 to explore the relation between TC and the material structures. The results suggest that the 2DFM materials with shorter distance between the magnetic atoms, larger local magnetic moments and more neighboring magnetic atoms are more likely to have higher TC.

Results

Two-dimensional magnetic atomic structural prototype

We select 198 2D magnetic materials which have high dynamical and thermodynamic stability in the Computational 2D Materials Database (C2DB)26,27,28. We also calculate 37 compounds from ref. 15, and 8 compounds from ref. 16. In addition to these materials, we also include several experimentally synthesized structures such as CrI32,3, Cr2Ge2Te64, Fe3GeTe25,6, and Fe4GeTe27. The workflow and calculation details are presented in the Methods Section. There are some 2D materials that show complicated magnetic ground states, which will be interesting for further studies. In this work, we focus on the FM materials, which have great potential in device applications.

We find 79 2D magnetic materials that have robust FM ground states. We analyze the structural characters of these 2DFM materials. For simplicity, we first consider only the magnetic atoms in the unit cell. The 79 2DFM structures can be divided into 11 categories according to the structural similarity of the magnetic atoms. The prototypical structures of the 11 categories are shown in Fig. 1. In Table 1, we list the chemical formula of materials as well as their space groups contained in each category. For convenience, we use “Chemical formula-Magnetic atom” to name the magnetic atomic structure prototypes. There are 6 categories, including Zr2I2S2-Zr, EuOI-Eu, V3N2O2-V, Cr3Te4-Cr, Fe3GeTe2-Fe, and Fe4GeTe2-Fe that each contains only one material. The remaining 73 materials are divided into five categories. Note that all the 2DFM materials studied in this work contain only one type of magnetic atom in each structure. The decoration of nonmagnetic atoms can change the structure dramatically. For example, if the nonmagnetic atoms are taken into account, the 29 original structures of the FeI2-Fe prototype can be further divided into 4 different structural types according to their structural similarity.

Fig. 1
figure 1

The structure prototypes of 2DFM materials with only magnetic atoms kept.

Full size image
Table 1 The structural prototypes of the 2DFM materials, and the chemical formula, space group of each compound.
Full size table

Among these 79 materials, some have multiple layers of magnetic atoms. Note that the layers we refer to here are those between the layers, the ions form strong chemical bonds, instead of weak van de Waals bonds. Among the 11 prototype categories, 4 have a single layer of magnetic atoms, 3 have double magnetic layers, 3 have triple magnetic layers, and 1 has quadruple layers. Obviously, the structures of different magnetic layers differ greatly.

The category of FeI2-Fe contains 29 2DFM materials, which is the most popular magnetic atomic structure prototype, and the magnetic atoms in this structure prototype reside on the same plane. The original structures have 3 different space groups P-6m2, P-3m1, and P3m1. The composition types of these structures are AB2 and ABC. The CrI3-Cr prototype with only one layer of magnetic atoms contains 14 2DFM materials which possess the P-62m and P-31m space groups, and adopt AB3 and ABC3 composition types. This category includes the famous CrI3 compound, which has been experimentally synthesized2, and investigated intensively.

For VO2-V prototype, the magnetic atoms also reside in one plane, and there are P4/nmm, P-4m2, and P4/mmm, three different space groups for the original structures. The composition types are AB and AB2. The Co2F2-Co prototype, which has double magnetic atom layers, has 6 2DFM materials and contains AB and ABC two different composition types. Their space groups are all P-3m1. There are 9 2DFM materials with ABC composition type and Pmmn space group belongs to Cr2I2S2-Cr prototype.

Curie temperatures

The calculated TCs for the 79 2DFM materials are also listed in Table 1. Since TCs of some of the 2DFM materials have been measured experimentally, we first compare the calculated TCs to the experimental values for these materials. The calculated TC of monolayer CrI3 is 26 K, which is slightly lower than 45 K of the experimental value2 and very close to the other theoretical calculation of 28 K17. The magnetic order at finite temperature can be stabilized when applying a small external magnetic field to Cr2Ge2Te64. The calculated TC of Cr2Ge2Te64 is 11 K with 0.25 meV/Cr magnetic anisotropic energy is very close to the estimated experimental value, and is slightly lower than the theoretical estimation of 20 K in ref. 4. It is worth mentioning that the calculated TC of Fe3GeTe2 is 158 K29, which is close to the experimental value of 130 K5. Our calculated TC of monolayer Fe4GeTe2 is about 207 K which is below 270 K of its bulk phase7. We also calculate the TC of monolayer Cr3Te4, and obtain TC = 133 K, which is also lower than the experimental TC = 316 K of the bulk material30, as expected. Note that in ref. 16, the TC of monolayer Cr3Te4 was severely overestimated to be 2057 K, much higher than the TC of the bulk Cr3Te4. All these results suggest that our calculated TCs of the 2DFM materials are highly reliable. The exchange interactions calculated by the linear response theory are somehow smaller than those calculated by the energy mapping method16,17,18, and therefore the calculated TCs are also lower.

The distributions of TCs of the first five categories of 2DFM materials are shown in Fig. 2. The average TC of the 73 materials is 71 K. The average TC for FeI2-Fe, CrI3-Cr, VO2-V, Co2F2-Co, and Cr2I2S2-Cr categories are 66 K, 38 K, 84 K, 132 K, and 88 K, respectively. The ratio of structures that have TC higher than the average value (71 K) for the above categories are 34%, 0%, 60%, 33%, and 78%, respectively. The last 6 prototypes in Table 1 have only one 2DFM material in each category. Among them, the TCs of V3N2O2, Cr3Te4, Fe3GeTe2, and Fe4GeTe2 are 71 K, 133 K, 158 K, and 207 K, respectively.

Fig. 2: The distribution of TC for 73 2DFM materials in different structure prototypes.
figure 2

The TCs are arranged in descending order. Blue dot line marks the average TC for the 73 structures. The percentages are the ratios of structures that have TC above the average Curie temperature for each category.

Full size image

Among the 79 2DFM materials, Co2F2 has the highest TC = 541 K. The TCs of the rest materials are all below room temperature, which is not very surprising. We find 3 materials whose TC are higher than 200 K, and 4 materials, whose TC are between 150 K and 200 K. There are 11 materials with TC between 100 K and 150 K. These results suggest that the TCs of 2DFM materials are relative low, compared to the 3D compounds, and finding room temperature 2DFM materials is quite difficult. The TCs of previous calculations16,18 are generally higher than the ones in this work. Especially, in ref. 16, heuristic factors 0.2–0.4 are multiplied to obtain reasonable TCs.

Discussion

Structure characteristics

It is important to understand the relation between the structure and TC. We have checked that the cation-anion-cation angles of all the 79 2DFM materials are close to 90∘, which obey the Goodenough-Kanamori rule31,32 for the FM exchange interactions. For example, the bond angles of Cr-O-Cr for CrO2 and Fe-Cl-Fe angle for FeCl2 are 97.53∘ and 85.49∘ respectively.

Naively, one could expect that having more neighboring magnetic atoms may lead to higher TC for a material. The coordination characterization function (CCF)33 defined in Eq. (6) of the Methods Section can reflect the coordination character of the structure, and the abscissa of the peak of CCF corresponds to the length of a atomic pair. The integral of CCF,

$${P}_{ij}(r)=\int\nolimits_{0}^{r}{{{{\rm{ccf}}}}}_{ij}(r^{\prime} )\,dr^{\prime} ,$$
(1)

characterizes the average number of (magnetic) atomic pairs within a certain range r, where \({{{{\rm{ccf}}}}}_{ij}(r^{\prime} )\) is the CCF for the i-th and j-th type of element of the structure.

The ccfij(r) (blue solid lines) and Pij(r) (red dashed lines) of four representative structures, CrI3, Co2F2, Fe3GeTe2, and Fe4GeTe2, from different structure prototypes are shown in Fig. 3. These structures have different numbers of magnetic atomic layers. CrI3 has quite low TC, which is about 26 K, whereas other components have much higher TC. Especially, the TC of Co2F2 is about 543 K, which is the highest among the 79 2DFM materials. Fe3GeTe2, and Fe4GeTe2 also have relative high TC. By comparing the ccfij(r) and Pij(r) of different structures, obviously CrI3 has much less neighboring magnetic atoms within the given range than other three components, due to two reasons. The first factor that affect ccfij(r) and Pij(r) is the length between the nearest magnetic atoms, \({d}_{\min }\). The smaller \({d}_{\min }\) will increase the neighboring magnetic atoms, and it may also lead to larger exchange interactions between magnetic atoms, which therefore leads to higher TC. The \({d}_{\min }\) of Co2F2 is very small (~2.48 Å), due to the small radii of F atoms. The \({d}_{\min }\)s of Fe3GeTe2 and Fe4GeTe2 are also quite small ~2.5 Å. In contrast, the \({d}_{\min }\) of CrI3 is about 3.97 Å, which is much larger than those of the other three compounds. The ccfij(r) and Pij(r) are also related to the number of magnetic atomic layers for the materials. For example, CrI3 has only a single layer of magnetic atoms, whereas Co2F2 has two magnetic atomic layers and Fe3GeTe2, Fe4GeTe2 have 3 and 4 magnetic atomic layers respectively, which therefore have more neighboring magnetic atoms. From Table 1, one may find that the 2DFM materials with multiple magnetic atomic layers are more likely to have higher TC. In fact, the TCs of the most compounds of a single magnetic atomic layer are relatively low, except CrO2 and FeCl2.

Fig. 3
figure 3

The ccfij(r) and their integral Pij(r) for a CrI3, b Co2F2, c Fe3GeTe2, and d Fe4GeTe2.

Full size image

The magnetization intensity is often used to evaluate bulk magnetic materials. For 2DFM materials, we define the areal magnetic moment density to measure their macroscopic magnetization,

$${{{{\bf{M}}}}}_{{{{\rm{2D}}}}}=\frac{{{{{\bf{m}}}}}_{{{{\rm{cell}}}}}}{{s}_{{{{\rm{cell}}}}}},$$
(2)

where mcell is the total magnetic moment of the unit cell, and scell is the area of the 2D unit cell. 2DFM materials with larger M2D have stronger magnetism and are more favorable for applications. Obviously, the 2DFM material with denser magnetic atomic mesh, multiple magnetic atomic layers, and larger local magnetic moments would have larger M2D. For example, the M2D for Co2F2, Fe4GeTe2, and Fe3GeTe2 are 0.74, 0.59, and 0.33 μB/Å2, respectively. They are much larger than that of CrI3 0.15 μB/Å2. This is because Co2F2 has two layers of magnetic atoms, which are densely packed. Fe4GeTe2 and Fe3GeTe2 have 4 and 3 layers of Fe atoms. In contrast, CrI3 has only one layer of magnetic atoms, with relatively larger Cr-Cr pair lengths. Interestingly, Co2F2 possesses both the highest TC and the highest M2D among the 79 2DFM materials, showing that it is a promising material.

SISSO analysis of Curie temperature

To dig out the factors that determine the TC of the 2DFM materials, we resort to the recently developed SISSO method which may capture the underlying mechanisms of materials’ properties using small data sets24,25. We use three parameters to construct the model, including the magnetic moment (M in μB per spin), the minimum magnetic atomic distance (\({d}_{\min }\), in Å), and the integral of CCF (P) as typical features to predict the TC of the above 2DFM materials. Among the three parameters, M is an important physical feature of the magnetic atoms, whereas \({d}_{\min }\) and P are used to characterize the structures of the magnetic atoms. The cutoff radius of Pij(r) is set to be 7.0 Å to calculate the P of the above 79 magnetic atomic structures.

Figure 4 illustrates the TC predicted by SISSO vs. the ab initio calculated TC. The one dimensional (1D) descriptor25 fitting gives,

$${T}_{{{{\rm{C}}}}}=27.88+110.03\times \frac{M\times P}{{d}_{\min }^{3}},$$
(3)

and 2D descriptor gives,

$${T}_{{{{\rm{C}}}}}=139.33+0.0078\times \frac{{M}^{3}\times {P}^{4}}{{d}_{\min }^{4}}-1.44\times \frac{{P}^{\frac{1}{4}}\times {d}_{\min }^{3}}{{M}^{\frac{1}{2}}}.$$
(4)

The root mean square error (RMSE) for 2D descriptor is 38 K, which is somehow smaller than that of the 1D descriptor of 46 K. According to the results of SISSO, it is clear that TC has a positive correlation with M and P, and negative correlation with \({d}_{\min }\), which is consistent with the analysis in the previous paragraphs.

Fig. 4: The TC predicted by SISSO vs. the TC calculated via the high-throughput first-principles calculations.
figure 4

The red circles are the results of 1D descriptor, whereas the green triangles are the results of 2D descriptor.

Full size image

The SISSO method correctly produces that Co2F2 has very high TC, because Co2F2 meets all the above beneficial conditions for high TC. The \({d}_{\min }\) for Co2F2 is only 2.48 Å, and the next nearest neighbor is as short as 2.73 Å. The exchange interaction between the nearest neighbor Co-Co pairs in Co2F2 is calculated to be −35.43 meV, and the exchange interaction between the next nearest neighbor Co-Co pair is −15.28 meV. Furthermore, there are about 20 Co-Co pairs with a mutual distance smaller than 7.0 Å for each Co atom, as shown in Fig. 3b. All these physical quantities suggest that it should have a high TC.

On the other hand, the TCs of the compounds in the CrI3-Cr structural prototype are quite low. For example, CrI3 has a low TC 26 K. It has only one magnetic atomic layer and for each Cr atom, there are only 5 Cr-Cr pairs within 7.0 Å. The \({d}_{\min }\) is as large as 3.97 Å. The situations for other structures of CrI3-Cr structural prototype are similar.

However, the magnetic interactions have very complicated mechanisms, and the TCs are determined by many factors. One would not expect that the simple SISSO models can predict very accurate TCs for the 2DFM materials. For example, while the SISSO method predicts the TCs of Fe3GeTe2 rather well, it somehow over estimates the TC of Fe4GeTe2. They also under estimate the TC of FeCl2 and CrO2. Nevertheless, it gives correct trends of TC respect to the structure of the 2DFM materials, which may provide useful guidance to further searching of high TC 2DFM materials.

Finally, we note that the PBE functional used in this work is suitable for large-scale, fast screening of the materials, but it is known to have some deficiencies. Once we find the promising candidate materials, we can use advanced techniques34,35,36,37 to obtain more accurate results for the electronic and magnetic structures. Furthermore, the 2D structures studied in this work are taken from existing databases. It is possible to investigate a wider range of 2D magnetic materials via crystal structure prediction methods38,39.

In this work, we carry out high-throughput first-principles computations to search for potential high TC 2DFM materials. We identify 79 2D materials that have robust FM ground state and calculate their TCs. Among the 79 2DFM materials, Co2F2 has the highest TC = 541K, well above room temperature. There are three 2D materials that also have relatively high TC, including CrO2 (226 K), FeCl2 (198 K), and Fe4GeTe2 (207 K). We perform SISSO method to analyze the relations between the TC and the structure of the 2DFM materials. The results suggest that the 2DFM materials with smaller distance between the magnetic atoms \({d}_{\min }\), larger local magnetic moments and more neighboring magnetic atoms are more likely to have higher TC. Our research is instructive to the discovery of new 2DFM materials with high TC.

Methods

Workflow

The work flow of our high-throughput computations is illustrated in Fig. 5. This high-throughput computing process is highly automated and almost no empirical parameters are introduced. Once the basic 2D magnetic crystals are collected, we run the Structure Prototype Analysis Package (SPAP) code33 to cluster structures according to structural similarity and label the structures with corresponding structure prototypes. We then perform self-consistent first-principle calculations, and some non-convergent structures will be removed at this step. The magnetic exchange interactions \({J}_{\alpha {{{\bf{R}}}},\beta {{{{\bf{R}}}}}^{\prime}}\) are calculated by a linear response method21,22. Based on the calculated magnetic exchange parameters, the magnetic ground state of the system is simulated by the Monte Carlo method at T ≈ 0 K. If the ground state is ferromagnetic, the TC will be calculated. We carry out further analysis based on these results.

Fig. 5
figure 5

The workflow of the high-throughput computation of the 2D magnetic materials.

Full size image

Structure prototype analysis

The properties of the materials are determined by their constituting elements and the residing structure. A group of materials having similar structures tend to possess similar properties. In order to find out the relationship between the TC (and other magnetic properties) of the 2DFM materials and their structures, we categorize the 2DFM materials according to their structural similarity. We use the SPAP code to classify these structures into corresponding structure prototypes, and details of the methods are given in ref. 33. Different from the similarity comparisons of bulk crystal structures, which needs to normalize two structures to the same volumetric atom number density, for 2D materials, we need to construct an areal atom number density normalization coefficient,

$${c}_{{{{\rm{2D}}}}}=\sqrt{\frac{{S}_{A}{N}_{B}}{{S}_{B}{N}_{A}}},$$
(5)

where SA, SB, NA, and NB denote the area of the 2D lattice of structure A and B, and the number of atoms in the 2D cell of structure A and B, respectively. We multiply the 3 lattice vectors of structure B by c2D. Note that the fractional coordinates of the atoms of the structure B should be kept unchanged. The manipulated structure B is adjusted to the same areal atom number density with that of structure A. The normalization helps us to identity similar structures with different scales. We use the CCF33 to measure the distance between normalized structures, which is defined as follows, for pairs of atomic types i and j,

$${{{{\rm{ccf}}}}}_{ij}(r)=\left(1-\frac{1}{2}{\delta }_{ij}\right)\frac{1}{N}\mathop{\sum}\limits_{{n}_{i}}\mathop{\sum}\limits_{{n}_{j}}f({r}_{{n}_{i}{n}_{j}})\sqrt{\frac{{a}_{{{{\rm{pw}}}}}}{\pi }}{e}^{-{a}_{{{{\rm{pw}}}}}{(r-{r}_{{n}_{i}{n}_{j}})}^{2}},$$
(6)

where N is the total number of atoms in the cell. ni runs over all atoms of the i-th type within the unit cell, whereas, nj runs over all atoms of the j-th type within the extended cell. \(f({r}_{{n}_{i}{n}_{j}})\) is the weighting function for different interatomic distances, whereas \({r}_{{n}_{i}{n}_{j}}\) is the interatomic distance within the cutoff radius. apw (usually 60.0 Å−2) is a parameter that controls the width of the normalized Gaussian smearing function.

We assign similar structures with the same prototype name with the help of three norms: (i) the same composition type and total number of atoms in the conventional cell; (ii) equivalent three-dimensional crystal symmetry; and (iii) the distance between two structures being below a threshold. In this work, we take 0.075 as the threshold. After the structural classification by SPAP, the structures are rechecked visually to guarantee that similar structures are labeled with the same prototype name.

Density functional calculation

The structural, electronic, and magnetic properties of the 2D magnetic materials are studied via density functional theory (DFT) by the WIEN2K code, which is based on a full-potential linearized augmented plane wave method (FP-LAPW)40. We adopt the Perdew-Burke-Ernzerhof (PBE) form generalized gradient approximation (GGA) of the exchange-correlation functional41. For all calculated materials, 1000 Brillouin k-points are used, and the number of k-points are automatically assigned according to the shape of the unit cell. The plane-wave cutoff RKmax = 7 is used.

Magnetic exchange interactions

In previous works16,17,18,42, the magnetic exchange interactions are usually calculated via an energy mapping method, i.e., the total energies of different spin configurations are calculated, and the exchange interactions are fitted to the energies of different spin configurations. There are several disadvantages of this method. First, this method is not user-friendly for the high-throughput calculations, since the spin configurations have to be chosen by hand for each material, with caution. Different choices of the spin configurations may lead to large discrepancies in the results. This is because, for some spin configurations, the electrons are forced to occupy higher energy orbitals, and in these cases, the total energies can not be fitted very well by the Heisenberg model, because an additional electron band (the so-called Stoner) energies have to be considered. The exchange interactions calculated by the energy mapping method may thus be overestimated and this might be one of the reasons that some of the previous works significantly overestimated the TC16. Furthermore, for some materials, the long-range exchange interactions are important for TC18. To take account of the long-range magnetic interactions, very large supercells are required to obtain the exchange interactions.

In this work, we calculate the magnetic interactions via a first-principle linear-response approach21,22, where the magnetic exchange interaction is calculated as:

$$\begin{array}{l}{J}_{i{{{\bf{R}}}},j{{{{\bf{R}}}}}^{\prime}}=\mathop{\sum}\limits_{{{{\bf{q}}}}}\mathop{\sum }\limits_{{{{\bf{k}}}}n{n}^{\prime}}\frac{{f}_{n{{{\bf{k}}}}}-{f}_{{n}^{\prime}{{{\bf{k}}}}+{{{\bf{q}}}}}}{{\epsilon }_{n{{{\bf{k}}}}}-{\epsilon }_{{n}^{\prime}{{{\bf{k}}}}+{{{\bf{q}}}}}}\left\langle {\psi }_{n{{{\bf{k}}}}}\left|\left[\sigma \times {{{{\bf{B}}}}}_{i}\right]\right|{\psi }_{{n}^{\prime}{{{\bf{k}}}}+{{{\bf{q}}}}}\right\rangle \\\qquad\quad \times\, \langle {\psi }_{{n}^{\prime}{{{\bf{k}}}}+{{{\bf{q}}}}}\left|\left[\sigma \times {{{{\bf{B}}}}}_{j}\right]\right|{\psi }_{n{{{\bf{k}}}}}\rangle {e}^{i{{{\bf{q}}}}\cdot \left({{{\bf{R}}}}-{{{{\bf{R}}}}}^{\prime}\right)},\end{array}$$
(7)

where, ϵnk is the one-electron band energy, σ are the Pauli matrices. ψnk is the Kohn-Sham wave function, Bi is the local magnetic field, R is the unit cell index, whereas i is the atom in the R-th unit cell. This technique has been successfully applied to a wide variety of complex magnetic materials22,29,43,44,45. Details about the method can be found in refs. 21,22. In this work, we calculate all the eligible exchange interactions \({J}_{i{{{\bf{R}}}},j{{{\bf{{R}}}^{\prime}}}}\), i.e., the exchange interactions between the i-th spin in the R-th unit cell and the j-th spin in the \({{{{\bf{R}}}}}^{\prime}\)-th unit cell, within the range of \(| {{{\bf{{R}}}_{\alpha }^{\prime}}}-{{{{\bf{R}}}}}_{\alpha }| \le 3\), where α = a, b.

Replica-exchange Monte Carlo simulation

The magnetic phase diagrams and the TCs are simulated via the following anisotropic Heisenberg-like Hamiltonian,

$$H=\mathop{\sum }\limits_{i{{{\bf{R}}}},j{{{{\bf{R}}}}}^{\prime}}{J}_{i{{{\bf{R}}}},j{{{{\bf{R}}}}}^{\prime}}{{{{\bf{S}}}}}_{i{{{\bf{R}}}}}\cdot {{{{\bf{S}}}}}_{j{{{{\bf{R}}}}}^{\prime}}+\mathop{\sum }\limits_{i}A{({{{{\bf{S}}}}}_{i}^{z})}^{2},$$
(8)

where SiR is the i-th spin in the R-th unit cell. The magnetic anisotropy energy A is important to stabilize the ferromagnetism in 2D. Because the magnetic anisotropy energy A is very sensitive to the calculation parameters, it is still quite difficult to obtain highly accurate and reliable values in the high-throughput DFT calculations. Fortunately, we find that TC is not very sensitive to the value of A, as demonstrated in our previous works29. For example, the TC of CrI3 with the magnetic anisotropic energy A = 0.25 meV for each Cr ion is 26 K, whereas the TC only increases to 30 K when A = 1.00 meV is used. Therefore, we set a reasonable value A = 0.25 meV per spin for all materials, which is enough for our purpose. We simulate the above spin Hamiltonian via the replica-exchange Monte Carlo method23,46,47, in which one simulates M replicas each at a different temperature T covering a range of interest, allowing configurational exchange between the replicas. The simulations are performed on a 36 × 36 lattice with periodic boundary conditions. We first perform the simulations at temperatures ranging from 5 to 640 K using 120 replicas. If no FM phase transition is found, we then perform simulation at temperatures ranging from 0.2 to 20 K. We discard the first 2 × 105 sweeps, when computing the equilibrium properties. Sample averages are accumulated over 2 × 105 sweeps.

SISSO analyses

SISSO is an approach for discovering descriptors for materials’ properties, which can tackle huge and possibly strongly correlated feature spaces, and converges to the optimal solution from a combination of features relevant to the materials’ target property. The outcome of SISSO is explicit, analytic functions consisted of basic physical quantities. A detailed description of the method can be found in refs. 24,25. In this work, the SISSO code from github.com/rouyang2017/SISSO is used to perform the analyses.