High-throughput computational screening for two-dimensional magnetic materials based on experimental databases of three-dimensional compounds

We perform a computational screening for two-dimensional (2D) magnetic materials based on experimental bulk compounds present in the Inorganic Crystal Structure Database and Crystallography Open Database. A recently proposed geometric descriptor is used to extract materials that are exfoliable into 2D derivatives and we find 85 ferromagnetic and 61 antiferromagnetic materials for which we obtain magnetic exchange and anisotropy parameters using density functional theory. For the easy-axis ferromagnetic insulators we calculate the Curie temperature based on a fit to classical Monte Carlo simulations of anisotropic Heisenberg models. We find good agreement with the experimentally reported Curie temperatures of known 2D ferromagnets and identify 10 potentially exfoliable 2D ferromagnets that have not been reported previously. In addition, we find 18 easy-axis antiferromagnetic insulators with several compounds exhibiting very strong exchange coupling and magnetic anisotropy.

We perform a computational screening for two-dimensional magnetic materials based on experimental bulk compounds present in the Inorganic Crystal Structure Database and Crystallography Open Database.A recently proposed geometric descriptor is used to extract materials that are exfoliable into two-dimensional derivatives and we find 85 ferromagnetic and 61 anti-ferromagnetic materials for which we obtain magnetic exchange and anisotropy parameters using density functional theory.For the easy-axis ferromagnetic insulators we calculate the Curie temperature based on classical Monte Carlo simulations of anisotropic Heisenberg models.We find good agreement with the experimentally reported Curie temperatures of known 2D ferromagnets and identify 10 potentially exfoliable two-dimensional ferromagnets that have not been reported previously.In addition, we find 18 easy-axis anti-ferromagnetic insulators with several compounds exhibiting very strong exchange coupling and magnetic anisotropy.

I. INTRODUCTION
The discovery of two-dimensional (2D) ferromagnetism in 2017 1,2 has initiated a vast interest in the field of the field.The origin of magnetic order in 2D is fundamentally different from the spontaneously broken continuous symmetry that is responsible for magnetism in three-dimensional materials.In particular, the Mermin-Wagner theorem states that a continuous symmetry cannot be broken at finite temperatures in 2D and magnetic anisotropy therefore becomes a crucial ingredient for magnetic order in 2D.The first report on 2D ferromagnetism involved a monolayer of CrI 3 , 1 which has a strong easy-axis orthogonal to the plane and has a Curie temperature of 45 K.In addition, few-layer structures of CrGeTe 3 was reported to exhibit ferromagnetic order down to the bilayer limit. 2 However, for the case of a monolayer of CrGeTe 3 magnetic order is lost due to the presence of an easy-plane, which comprises a continuous symmetry that cannot be broken spontaneously.Since then several materials have joined the family of 2D magnets.Most notably, CrBr 3 , 3 which have properties very similar to CrI 3 but with lower Curie temperatures of 34 K due to smaller magnetic anisotropy, Fe 3 GeTe 2 , which is metallic and has a Curie temperature of 130 K 4 , FePS 3 5 which is anti-ferromagnetic with an ordering temperature of 118 K, and VSe 2 where some evidence has been provided for ferromagnetic order at room temperature 6 although the presence of magnetism is being debated 7 .10][11][12][13][14][15][16] Although the handful of known magnetic 2D materials have been shown to exhibit a wide variety of interesting physics, there is a dire need for discovering new materials with better stability at ambient conditions and higher critical temperatures for magnetic order.Such conditions are not only crucial for technological applications of 2D magnets, but could also serve as a boost for the experimental progress.In addition, the theoretical efforts in the field are largely limited by the few materials that are available for comparison between measurements and calculations.An important step towards discovery of novel 2D materials were taken by Mounet et al. 17 where Density Functional Theory (DFT) was applied to search for potentially exfoliable 2D materials in the Inorganic Crystal Structure Database (ICSD) and the Crystallography Open Database (COD).More than 1000 potential 2D materials were identified and 56 of these were predicted to have a magnetically ordered ground state.Another approach towards 2D materials discovery were based on the Computational 2D Materials Database (C2DB), [18][19][20] which comprises more than 3700 2D materials that have been computationally scrutinized based on lattice decoration of existing prototypes of 2D materials.The C2DB presently contains 152 ferromagnets and 50 anti-ferromagnets that are predicted to be stable by DFT.In addition to these high throughput screening studies there are several reports on particular 2D materials that are predicted to exhibit magnetic order in the ground state by DFT, [21][22][23][24][25][26] as well as a compilation of known van der Waals bonded magnetic materials that might serve as a good starting point for discovering new 2D magnets. 27ue to the Mermin-Wagner theorem a magnetically ordered ground state does not necessarily imply magnetic order at finite temperatures and the 2D magnets discovered by high throughput screening studies mentioned above may not represent materials with observable magnetic properties.In three-dimensional bulk compounds the critical temperature for magnetic order is set by the magnetic exchange coupling between magnetic moments in the compound and a rough estimate of critical temperatures can be obtained from mean field theory. 28In 2D materials, however, this is no longer true since magnetic order cannot exist with magnetic anisotropy and mean field theory is always bound to fail.The critical temperature thus has to be evaluated from either classical Monte Carlo simulations or renormalized spin-wave theory of an anisotropic Heisenberg model derived from first principles 2,[29][30][31] The former approach neglects quantum effects whereas the latter approximates correlation effects at the mean field level.Monte Carlo simulations are not well suited to high-throughput studies, but it has recently been shown that such calculations can be fitted to an analytical expression that is easily evaluated for a given material once the exchange and anisotropy parameters have been computed. 30,32This approach has been applied to the C2DB resulting in the discovery of 11 new 2D ferromagnetic insulators that are predicted to be stable. 33In addition 26 (unstable) ferromagnetic materials with Curie temperatures esceeding 400 K have been identified from the C2DB. 34However, it is far from obvious that any of these materials can be synthesised in the lab even if DFT predicts them to be stable since they are not derived from experimentally known van der Waals bonded bulk compounds.
In the present work we have performed a full computational screening for magnetic 2D materials based on experimentally known van der Waals bonded materials present in the ICSD and COD.In contrast to previous high throughput screening of these databases we evaluate exchange and magnetic anisotropy constants for all materials with a magnetic ground state and use these to predict the Curie temperature from an expression fitted to Monte Carlo simulation of the anisotropic Heisenberg model.

II. METHODOLOGY
The first step in the computational screening is to identify potentially exfoliable 2D structures from the bulk materials present in ICSD and COD.In Ref. 17 this was accomplished by identifying layered chemically bonded sub-units followed by a calculation of the exfoliation energy from van der Waals corrected DFT.Here we will instead use a recently proposed purely geometrical method that quantifies the amount of zero-dimensional (0D), one-dimensional (1D), two-dimensional (2D) and three-dimensional (3D) components present in a given material. 35The method thus assigns a 0D, 1D, 2D, and 3D score to all materials and thus quantifies the 0D, 1D, 2D, and 3D character.The scores are defined such that they sum to unity and taking the 2D score > 0.5 thus provides a conservative measure of a material being (mostly) composed of 2D components that are likely to be exfoliable.
The magnetic properties of possible candidate 2D materials are then investigated using first principles Heisenberg models derived from DFT. 2,[29][30][31]36 In particular, if a 2D candidate material has a magnetic ground state we model the magnetic properties by the Hamiltonian where J is the nearest neighbor exchange coupling, λ is the nearest neighbor anisotropic exchange coupling, A is the single-ion anisotropy, and ij denotes sum over nearest neighbors.J may be positive(negative) signifying a ferromagnetic(anti-ferromagnetic) ground state and we have assumed that the z-direction is orthogonal to the atomic plane and that there is in-plane magnetic isotropy.This model obviously does not exhaust the possible magnetic interactions in a material, 37 but has previously been shown to provide good estimates of the Curie temperature of CrI 3 29,30 and provides a good starting point for computational screening studies.
The thermal properties can then be investigated from either renormalized spin-wave calculations [28][29][30]38,39 or classical Monte Carlo simulations, 30,40 based on the model (1) Due to the Mermin-Wagner theorem the magnetic anisotropy constants are crucial for having magnetic order at finite temperatures and for ferromagnetic compounds the amount of anisotropy can be quantified by the spin-wave gap where S is the maximum eigenvalue of S z i and N nn is the number of nearest neighbors.This expression was calculated by assuming out-of-plane magnetic order and in the present context a negative spin-wave gap signals that the ground state favors in-plane alignment of spins in the model (1) and implies that the assumption leading to Eq. (2) breaks down.Nevertheless, the sign of the spinwave gap comprises an efficient descriptor for the presence of magnetic order at finite temperatures in 2D, since a positive value is equivalent to having a fully broken rotational symmetry in spin-space.
For bipartite lattices with anti-ferromagnetic ordering (J < 0) the spinwave analysis based on Eq. (1) (with out-of-plane easy-axis) yields a spinwave gap of It is straightforward to show that ∆ AFM is real and positive if (2S − 1)A > N nn Sλ, real and negative if (2S − 1)A < N nn S(2J + λ) and imaginary otherwise.The latter case corresponds to favouring of in-plane antiferromagnetic order and negative real values correspond to favouring of ferromagnetic order (this may happen if λ is a large positive number even if J < 0).∆ AFM thus only represents the physical spinwave gap in the case where it is positive and real.However, in the case of an imaginary spinwave gap the norm of the gap may be used to quantify the strength of confinement to the plane.In the case of non-bipartite lattices we use the expression (3) as an approximate measure of the anisotropy.More details on this can be found in Sec.IV A.
In Ref. 30 it was shown that the critical temperature for ferromagnetic order (J > 0) can be accurately obtained by classical Monte Carlo simulations of the model (1) and for S > 1/2 the result can be fitted to the function where and γ = 0.033.T Ising C is the critical temperature of the corresponding Ising model (in units of JS 2 /k B ).The expression ( 5) is readily evaluated for any 2D material with a ferromagnetic ground state once the Heisenberg parameters J, λ and A have been determined.This can be accomplished with four DFT calculations of ferromagnetic and anti-ferromagnetic spin configurations including spin-orbit coupling.Specifically, for S = 1/2 the exchange and anisotropy constants are determined by 33,41 where ∆E FM(AFM) = E FM(AFM) − E ⊥ FM(AFM) are the energy differences between in-plane and out-of-plane magnetization for ferromagnetic(anti-ferromagnetic) spin configurations and N FM(AFM) is the number of nearest neighbors with aligned(anti-aligned) spins in the antiferromagnetic configuration.For bipartite magnetic lattices (square and honeycomb) N FM = 0.However, several of the candidate magnetic materials found below contain a triangular lattice of transition metal atoms and in that case there is no natural anti-ferromagnetic collinear structure to compare with and we have chosen to extract the Heisenberg parameters using a striped anti-ferromagnetic configurations with N FM = 2 and N AFM = 4. Finally the factor of (1 + β/2S) in the denominator of Eq. ( 9) accounts for quantum corrections to anti-ferromagnetic states of the Heisenberg model where β is given by 0.202 and 0.158 for N AFM = 3 (honeycomb lattice) and N AFM = 4 (square and triangular lattices) respectively. 41For S = 1/2 we take A = 0 and λ = ∆E FM /N S 2 for J > 0 and λ = −∆E AFM /(N AFM − N FM )S 2 for J < 0.More details on the energy mapping analysis is provided in appendix A. All DFT calculations were performed with the electronic structure package GPAW 42,43 including non-selfconsistent spinorbit coupling 44 and the Perdew-Burke-Ernzerhof 45 (PBE) functional.

A. Computational screening of COD and ICSD
The ICSD and COD databases combined count more than 500.000materials, but removing corrupted or incomplete entries and duplicates, reduces the number to 167767 bulk materials. 35Of these, a subset of 4264 are predicted to have a 2D score higher than 0.5 and these materials are the starting point of the present study.We then perform a computational exfoliation by isolating the 2D component and performing a full relaxation of the resulting 2D material with DFT.We restrict ourselves to materials that have a 2D component with less than five different elements and less than a total of 20 atoms in the minimal unit cell.This reduces the number of candidate 2D materials to 651 compounds.We find 85 materials with a ferromagnetic ground state and 61 materials with an anti-ferromagnetic ground state.A schematic illustration of the workflow is shown in Fig. 1.
For all of the magnetic materials we calculate the exchange coupling J and the spinwave gap ∆ according to the energy mapping approach. 31,33,46The results are shown in Fig. 2 and all the materials along with the calculated parameters can be found in Tabs.V-VIII.The spinwave gap is on the order of 0-4 meV for all materials.The exchange couplings fall in the range of 0-10 meV for the insulators but can acquire somewhat larger values for the metals.However, the energy mapping analysis is somewhat ill-defined for metals, since the electronic structure may change significantly when comparing energy differences between ferromagnetic and anti-ferromagnetic configurations.In particular, for insulators the value of S is a well-defined integer that can be extracted from the ferromagnetic ground state without spin-orbit coupling.But for metals it is not clear what value to use in the model (1).In addition the Heisenberg model itself is likely to be unsuitable for a description of the magnetic properties of metals and we restrict ourselves to insulators in the following and then subsequently comment on promising metallic compounds.

B. Insulating 2D ferromagnets
In Tab.I we display the calculated exchange coupling constants and spinwave gaps for ferromagnetic insulators with ∆ > 0. Assuming in-plane magnetic isotropy these are the only insulators that will exhibit magnetic order at finite temperatures.For the compounds with S = 1/2 we calculate the Curie temperatures according to Eq. (5).
It is reassuring that the well-known Ising type 2D ferromagnets CrBr 3 3 and CrI 3 1 are reproduced by the screening.In addition, CrClO, CrCl 3 , MnO 2 , CoCl 2 , and NiI 2 have previously been predicted to be ferromagnetic 2D insulators by DFT. 18,33,46Multi-layered CrSiTe 3 has been reported to exhibit a large magnetic anisotropy in the direction perpendicular to the layers and a ferromagnetic phase transition has been observed at 33 K. 47 In addition, strained CrSiTe 3 has very recently been predicted to comprise an ideal candidate for a 2D Kitaev  I. List of 2D ferromagnetic insulators (J > 0) with out-of-plane easy axis (∆ > 0).The Curie temperature for materials with S = 1/2 was calculated from Eq. ( 5).
spin spin-liquid. 48e also find 10 novel 2D ferromagnetic insulators - have not been studied prior to the present work.Of particular interest is the compound CoCa 2 O 3 , which is predicted to be ferromagnetic up to 57 K.However, it exhibits a rather small band gap of 40 meV, which may imply that the electronic structure could be sensitive to the choice of exchange-correlation functional.Such ambiguities have indeed been reported for FeCl 3 and FeBr 3 , which are both predicted to be small-gap quantum anomalous Hall insulators by PBE, but trivial insulators by PBE+U as well as other GGA functionals. 19he largest exchange coupling constant in Tab.I of 11 meV is found for MnNa 2 P 2 F 3 O 7 , which appears highly promising.However, we do not have a reliable estimate for the critical temperature due to large in-plane anisotropy (only two nearest neighbors per Mn atom), which renders the inclusion of second nearest neighbors crucial.A faithful estimation of the critical temperature would thus require a full Monte Carlo simulation of an extended Heisenberg model including in-plane anisotropy and exchange couplings.This is, however, beyond the scope of the present screening study.
The materials NiRe 2 O 8 and CoCl 2 O 8 are interesting variants of the common CdI 2 prototype (for example NiI 2 ) where the halide atom is replaced by units of ReO 4 and ClO 4 respectively.For 2D materials discovery based on computational lattice decoration such compounds opens the possibility of a wide range of new materials, since the number of possible ligands in the CdI 2 prototype is dramatically increased.
We also wish to mention the compound CuC 6 H 4 N 6 O 2 , which is an example of a 2D metal-organic framework (MOF).It is composed of a rectangular lattice of Cu atoms connected by pyrazine (C 4 H 4 N 2 ) and C 2 N 4 O 2 units.Such 2D MOFs have recently attracted an increasing amount of attention and it has been shown that the quasi-2D MOF CrCl 2 (pyrazine) 2 exhibits ferrimagnetic order below 55 K. 49 Due to the spin-1/2 nature of the magnetic lattice we cannot obtain a reliable estimate of the critical temperature of this material.Moreover, the material have large in-plane anisotropy and the second nearest neighbors must play a crucial role since the nearest neighbor approximation gives rise to chains that cannot order themselves at finite temperatures.Nevertheless the sizable value of the intrachain exchange coupling (3.04 meV ) could imply a critical temperature comparable to that of CrI 3 .
It should be stressed that the results of a screening study like the present one should be taken as a preliminary prediction.The first principles description of magnetic insulators is challenging for DFT since many of these exhibit strong correlation of the Mott-Hubbard type and the calculated Heisenberg parameters may be rather sensitive to the choice of functional. 31,33A detailed study of the functional dependence or inclusion of Hubbard corrections is required in order to support the theoretical prediction of these 2D materials being ferromagnetic.

C. Itinerant 2D ferromagnets
For metallic materials the prediction of thermodynamical properties is more challenging since it is not obvious that the Heisenberg Hamiltonian (1) comprises a good starting point for the analysis.Nevertheless, the exchange coupling J and spin-wave gap ∆ still provides a rough measure of the magnetic interactions and magnetic anisotropy respectively.Alternatively, one could specify the energy difference per magnetic atom in ferromagnetic and anti-ferromagnetic configurations as well as the energy cost of rotating the magnetic moments from the out-of-plane direction to the atomic plane.However, for the sake of comparison we have chosen to report the values of J and ∆ resulting from the energy mapping analysis although it comprises a rather naive approach for metals.The value of S is obtained by rounding off the total magnetic moment per atom to nearest half integer and we then evaluate the critical temperature from Eq. ( 5), which is the prediction obtained by assuming a Heisenberg model description using the calculated parameters.The results are shown in Tab.II, but it should be kept in mind that the exchange coupling constants and predicted critical temperatures in this case only provides a qualitative measure of the magnetic interactions.
Again, we rediscover a few materials (FeTe and VBrO) that were previously predicted to be ferromagnetic from  II.List of 2D itinerant ferromagnets (J > 0 and EGap = 0) with out-of plane easy axis (∆ > 0).The Curie temperature for materials with S = 1/2 was calculated from Eq. ( 5).
computational screening of the C2DB.FeClO has recently been exfoliated to bilayer nanoflakes and were shown to retain the anti-ferromagnetic ordering known from the bulk material. 50The discrepancy with our prediction of ferromagnetic order could either be due to an inaccurate description by PBE or due to the fact that the true anti-ferromagnetic structure of bulk FeClO is strongly non-collinear, 51 which is not taken into account in the present simplistic calculations.We find a few materials with two nearest neighbors, implying a strongly anisotropic in-plane magnetic lattice.For example, VFC 4 O 4 (H 2 O) 2 is a MOF with hydrated alternating linear chains of V and F atoms interconnected by cyclobutanetetrone (C 4 O 4 ) units.The intra-chain exchange coupling is significant (22.3 meV), but a reliable estimate of the critical temperature requires inclusion of the inter-chain exchange, which is not addressed in the present study.We also find a few materials with 9 nearest neighbors, which originates from a strongly buckled lattice of magnetic atoms and the analysis based on nearest neighbor interactions is expected to be insufficient in this case as well.We observe that several materials have predicted exchange couplings on the order of 10-50 meV, which far exceeds the values found for the insulators.But it should be emphasized that the comparison is not necessarily fair since the electronic structure of the antiferromagnetic state may be significantly different compared to the ferromagnetic state.Such differences will lead to large predictions for J that do not originate from magnetic interactions.Nevertheless, Tab.II provides a promising starting point for the discovery of new 2D itinerant ferromagnets, but there is a dire need for a better theoretical framework that can quantitatively deal with the thermodynamical properties of itinerant magnetism in 2D.
We finally note that certain known itinerant 2D ferromagnets (VSe 2 6 and CrGeTe 3 2 ) are not present in Tabs.I and II due to in-plane magnetization, which results in a negative spinwave gap in the present study.For the case of CrGeTe 3 this is in accordance with the experimentally observed loss of magnetism in the monolayer limit whereas for VSe 2 the origin of magnetic order is still unresolved. 7In addition, we do not find the itinerant 2D ferromagnet Fe 3 GeTe 2 , 4 which cannot be found in a bulk parent form in either the COD or ICSD.

D. Insulating 2D anti-ferromagnets
In the case of anti-ferromagnetic insulators we do not have a quantitative estimate of the Néel temperature given the nearest neighbor exchange coupling and spinwave gap.However, it is clear that an easy-axis (positive spinwave gap) is required to escape the Mermin-Wagner theorem for materials with isotropic in-plane magnetic lattices.Moreover, although the formula for the critical temperature Eq. ( 5) was fitted to Monte Carlo simulations we expect that a rather similar expression must be valid for the Néel temperature of anti-ferromagnets.This is partly based on the fact that mean field theory yields similar critical temperatures for ferromagnetic and antiferromagnetic interactions in the nearest neighbor model and we thus use the expression (5) as a very rough estimate of the critical temperatures for the anti-ferromagnet candidates found in the present work.In Tab.III we thus display a list of the anti-ferromagnetic insulators with positive spin-wave gap.In addition to the exchange coupling and spin-wave gap we also report the critical temperatures calculated from Eq. ( 5).
The most conspicuous result is the exchange coupling of VPS 3 , which exceeds 0.1 eV.However, while the use of the energy mapping analysis seems to be justified by the gapped anti-ferromagnetic ground state, the ferromagnetic configuration entering the analysis is metallic and may thus imply that the energy difference is not solely due to magnetic interactions.Nevertheless, the local magnetic moments in the ferromagnetic and antiferromagnet states are almost identical, which indicates that the large energy difference between the ferromagnetic and anti-ferromagnetic states originates in magnetic interactions.
We also observe that the V and Mn halides are predicted to be anti-ferromagnetic insulators with large exchange coupling constants.However, these compounds exhibits the CdI 2 prototype where the magnetic atoms form a triangular lattice.In the present study we have only considered collinear spin configurations, but the true ground state of a triangular lattice with antiferromagnetic nearest neighbor exchange has to exhibit a frustrated non-collinear spin structure. 52Second-nearest neighbors may complicate this picture and the true ground state of these materials could have a complicated structure.Moreover, it has previously been shown that the Mn halides are predicted to be ferromagnetic with the PBE+U functional, which underlines the importance of further investigating the predictions of the present work with respect to exchange-correlation functional, second nearest neighbor interactions etc.
In analogy with the ferromagnetic insulators NiRe ) 2 ) units.Again, the material exhibits strong nearest neighbor interactions (across oxalate units), but the second nearest interactions (mediated by ethylenediamine units) will play a crucial role in determining the critical temperature, which is predicted to vanish in the present study, which is only based on nearest neighbor interactions.
Finally, we remark that MnBi 2 Te 4 in 3D bulk form has recently attracted significant attention as it has been demonstrated to comprise the first example of a magnetic Z 2 topological insulator. 53,54The bulk material is comprised of ferromagnetic layers with anti-ferromagnetic interlayer coupling.In contrast we predict that the individual layers exhibit anti-ferromagnetic order.Like the case of the Mn halides the sign of the exchange coupling constant changes upon inclusion of Hubbard corrections to the DFT description.We have tested that PBE+U calculations yields ferromagnetic ordering for U > 2.0 eV.In addition, we do not find the Ising anti-ferromagnet FePS 3 , 5 , since PBE without Hubbard corrections predicts this material to be non-magnetic.This could imply that PBE+U is likely to be a more accurate framework for the present type of calculations, but we leave it to future work to unravel the sensitivity to the choice of xc-functional used for the DFT calculations.

E. Itinerant 2D anti-ferromagnets
For completeness we also display all the predicted antiferromagnetic metals in Tab.(IV).For S = 1/2, we have provided rough estimates of the critical temperatures based on Eq. ( 5), but in this case it should be regarded as a simple descriptor combining the effect of exchange and anisotropy rather than an actual prediction for the critical temperature.Neither the energy mapping analysis or the Heisenberg model is expected to comprise good approximations for these materials.However, DFT  IV.List of 2D itinerant anti-ferromagnets (J < 0) with out-of-plane easy axis (∆ > 0).The Curie temperature for materials with S = 1/2 was calculated from Eq. ( 5).
(with the PBE functional) certainly predicts that these materials exhibit ferromagnetic order at some finite temperature and Tab (IV) may provide a good starting point for further investigation or prediction of itinerant antiferromagnetism in 2D.

IV. DISCUSSION
We have performed a computational screening for 2D magnetic materials based on 3D bulk materials present in the ICSD and COD.We find a total of 85 ferromagnetic and 61 anti-ferromagnetic materials, which are listed in Tabs.V-VIII.The strength of magnetic interactions in the materials have been quantified by the nearest neighbor exchange coupling constants and the magnetic anisotropy has been quantified by the spinwave gap derived from the anisotropic Heisenberg model (1).Due to the Mermin-Wagner theorem only materials exhibiting an easy-axis (positive spinwave gap) will give rise to magnetic order at finite temperatures and these materials have been presented in Tabs.I-IV.For these we have also estimated the critical temperature for magnetic order from an expression that were fitted to classical Monte Carlo simulations of the anisotropic Heisenberg model.
The insulating materials are expected to be well described by the Heisenberg model and for S = 1/2 we have evaluated the critical temperatures from an analytical expression fitted to classical Monte Carlo simulations.However, for simplicity this expression was based on a Heisenberg model with in-plane isotropy and nearest neighbor interactions only.This may introduce errors in the prediction of critical temperatures, but for any given material the approach is easily generalized to include other interactions and in-plane anisotropy, which will yield more accurate predictions for critical temperatures.
A more crucial challenge is related to the determination of Heisenberg parameters from DFT.We have already seen that PBE+U can modify the predictions significantly 33 and even change the sign of the exchange coupling.Is is, however, not obvious that PBE+U will always provide a more accurate prediction compared to PBE (or other exchange-correlation functional for that matter) and benchmarking of such calculations is currently limited by the scarceness of experimental observations.
For anti-ferromagnetic insulators, we expect that classical Monte Carlo simulations combined with the energy mapping analysis will provide an accurate framework for predicting critical temperatures.In the present work we have simply used the expression (5) as a crude descriptor and leave the Monte Carlo simulations for antiferromagnets to future work.In general, the phase diagrams for anti-ferromagnets will be more complicated compared to ferromagnets 52 and there may be several critical temperatures associated with transitions between different magnetic phases.
The case of itinerant magnets are far more tricky to handle by first principles methods.It is not expected that the applied energy mapping analysis comprises a good approximation for metallic materials and it is not even clear if the Heisenberg description and associated Monte Carlo simulations is the proper framework for such systems.A much better approach would be to use Greens function methods 55,56 or frozen magnon calculations to access J(q) ∼ i J 0i e iq•R01 directly from which the magnon dispersion can be evaluated directly.It may then be possible to estimate critical temperatures based on renormalized spinwave theory 28 or spin fluctuation theory. 57espite the inaccuracies in the predicted critical temperatures of the present work, all of the 146 reported magnetic materials constitute interesting candidates for further scrutiny of 2D magnetism.All materials are likely to be exfoliable from bulk structures and contains mag-netic correlation in some form.Even the materials with an isotropic magnetic easy-plane that cannot host strict long-range order according to the Mermin-Wagner theorem, may be good candidates for studying KosterLitz-Thouless physics 58 Moreover, such materials exhibit algebraic decay of correlations below the Kosterlitz-Thouless transition, which may give rise to finite magnetization for macroscopic flakes. 31,59PENDIX A. Energy mapping analysis Here we provide the details of Eqs. ( 7)-( 9) used to extract the Heisenberg parameters from first principles.The energy mapping analysis is based on ferromagnetic and anti-ferromagnetic configurations.We only consider nearest neighbor interactions and in the number of nearest neighbors in the ferromagnetic configurations is denoted by N .Only bipartite lattices allow for antiferromagnetic configurations where all magnetic atoms have anti-parallel spin alignments with all nearest neighbors.For non-bipartite lattices we thus consider frustrated configurations where each atom has N FM nearest neighbors with parallel spin alignment and N AFM nearest neighbors with anti-parallel spin alignment.Assuming a classical Heisenberg description represented by the model ( 1), the ferromagnetic (FM) and antiferromagnetic (AFM) DFT energies per magnetic atom with in-plane ( ) and perpendicular spin configurations are written as where E 0 represents a reference energy that is independent of the magnetic configuration.The Heisenberg parameters can then be calculated as where ∆E FM(AFM) = E FM(AFM) − E ⊥ FM(AFM) are the energy differences between in-plane and out-of-plane magnetization for ferromagnetic(anti-ferromagnetic) spin configurations.
However, we wish to base the energy mapping on the quantum mechanical Heisenberg model, which is less trivial.If we start with the anisotropic Heisenberg model where spin-orbit coupling is neglected the ferromagnetic configuration with energy E FM corresponds to an eigenstate with energy −J/2N AFM per magnetic atom, which is the same as the classical Heisenberg model.However the anti-ferromagnetic configuration does not correspond to a simple eigenstate of the Heisenberg model.In particular, for bipartite lattices the Neel state where all sites host spin that are eigenstates of S z is not the eigenstate of lowest(highest) energy of the Heisenberg Hamiltonian model with J < 0(J > 0).Rather the classical energy corresponds to the expectation value of the Heisenberg Hamiltonian with respect to this state whereas the true ground state has lower(higher energy) leading to an overestimation of J if the energy mapping is based on the classical Heisenberg model.We have recently shown how to include quantum corrections to J for bipartite lattices using a correlated state, which has an energy in close proximity to the true anti-ferromagnetic ground state. 41We note that the magnetic moments obtained with DFT support the fact that the DFT energy of the anti-ferromagnetic configuration represents a proper eigenstate of the Heisenberg model rather than the classical state.The result is the factor of (1 + β/2S) in equation ( 9).
Including spin-orbit coupling and magnetic anisotropy in the energy mapping complicates the picture since only one of the states E FM , E ⊥ FM represents an eigenstate of the anisotropic Heisenberg model.On the DFT side this is reflected by the fact that only one of these configurations would be obtainable as a self-consistent solution and we have to calculate these energies by including spin-orbit coupling non-self-consistently.We thus retain the classical expression for the anisotropy constants, but retain the quantum correction for the exchange constants.Is is, however, clear that the single-ion anisotropy term becomes a constant for any system with S = 1/2.Since A does not have any physical significance it cannot influence the values of E FM(AFM) and E ⊥ FM(AFM) and we take A = 0 and λ = ∆E FM /N S 2 for J > 0 and λ = −∆E AFM /(N AFM − N FM )S 2 for J < 0. In principle, the two choices for λ should be equivalent and we have tested that they yield nearly the same value for a few spin-1/2 insulators.But in order to obtain full consistency with the spinwave gap we use different expressions depending on the sign of J.In addition for S = 1/2 the classical analysis leads to an inconsistency since the spinwave gap (2) is not guaranteed to yield the same sign as −∆E FM .This can be fixed by taking 2S → (2S − 1)S in Eq. ( 14), which leads to Eq. (7).Finally, the antiferromagnetic spinwave gap Eq. ( 3) was derived for bipartite lattices and it is not possible to derive a gap for non-bipartite lattices in a collinear spin configuration, since such a state will not represent the ground state leading to an instability in the gap.However, we will apply the expression naively to non-bipartite lattices as well but taking N nn → N AFM − N FM to ensure that the sign of the gap corresponds to the sign of −∆E AFM .

B. List of predicted magnetic materials
In Tab.V we list all the predicted ferromagnetic insulators containing two elements and tn Tab.VI we list the ferromagnetic materials containing three, four or five elements.For all materials we provide the COD/ICSD identifier for the bulk parent compound from which the 2D material was derived.We also stats the spin S, the number of nearest neighbors N nn , the exchange coupling J, the spinwave gap ∆, and Kohn-Sham band gap E Gap .For materials with S = 1/2 and N nn = 2 we have calculated an estimated critical temperature from Eq. ( 5).In Tab.VII we show all the anti-ferromagnetic compounds found in the computational screening.In addition, we found 11 materials for which we were not able to evaluate exchange coupling constants.This was either due to problems converging the anti-ferromagnetic spin configuration (converged to ferromagnetic state), more than two magnetic atoms in the unit cell, or that the two magnetic atoms in the unit cell form a vertical dimer.All of the materials are, however, predicted to be magnetic and could comprise interesting magnetic 2D materials that are exfoliable from 3D parent compounds.

FIG. 1 .
FIG. 1. Schematic workflow of the computational discovery of 2D magnets performed in the present work.

FIG. 2 .
FIG.2.Exchange coupling J and spinwave gap ∆ calculated for the magnetic 2D materials obtained from computational screening of ICSD and COD.
2 O 8 and CoCl 2 O 8 the anti-ferromagnetic insulator CoRe 2 O 8 comprises a variant of the CdI 2 prototype (represented by the V and Mn halides in Tab.III) where the halide atom has been replaced by ReO 4 .NiC 2 O 4 C 2 H 8 N 2 , constitutes an anti-ferromagnetic example of a MOF with a rectangular lattice of Ni atoms connected by a network of oxalate (C 2 O 4 ) and ethylenediamine (C 2 H 4 (NH 2

TABLE V .
List of 2D materials with a ferromagnetic ground state (within the PBE approximation) containing two elements.ID denotes the unique ICSD/COD identifier (materials from ICSD have ID < 10 6 ) for the bulk parent material and J is the nearest neighbor exchange interaction obtained from the energy mapping.EGap denotes the electronic (Kohn-Sham) band gap.∆ is the spin wave gap obtained from the anisotropy constants and positive values indicate an out-of-plane easy axis.

TABLE VIII .
List of 2D ferromagnetic compounds, which did not allow for a simple estimation of a nearest neighbor exchange coupling constant.