Computational Search for Magnetic and Non-magnetic 2D Topological Materials using Unified Spin-orbit Spillage Screening

Two-dimensional topological materials (TMs) have a variety of properties that make them attractive for applications including spintronics and quantum computation. However, there are only a few such experimentally known materials. To help discover new 2D TMs, we develop a unified and computationally inexpensive approach to identify magnetic and non-magnetic 2D TMs, including gapped and semi-metallic topological classifications, in a high-throughput way using density functional theory-based spin-orbit spillage, Wannier-interpolation, and related techniques. We first compute the spin-orbit spillage for the ~1000 2D materials in the JARVIS-DFT dataset (https://www.ctcms.nist.gov/~knc6/JVASP.html ), resulting in 122 materials with high-spillage values. Then, we use Wannier-interpolation to carry-out Z2, Chern-number, anomalous Hall conductivity, Curie temperature, and edge state calculations to further support the predictions. We identify various topologically non-trivial classes such as quantum spin-hall insulators (QSHI), quantum anomalous-hall insulators (QAHI), and semimetals. For a few predicted materials, we run G0W0+SOC and DFT+U calculations. We find that as we introduce many-body effects, only a few materials retain non-trivial band-topology, suggesting the importance of high-level DFT methods in predicting 2D topological materials. However, as an initial step, the automated spillage screening and Wannier-approach provide useful predictions for finding new topological materials and to narrow down candidates for experimental synthesis and characterization.


Introduction
In recent years, there has been a huge upsurge in topological materials research, following the predictions and observations of Dirac, Weyl, and Majorana fermions in condensed matter systems 1,2 . Several classes of topological materials have been proposed for applications in errorreduced quantum computing [3][4][5][6] , or as high mobility and dissipationless conductors. While there have been several recent detailed screening efforts for 3D non-magnetic topological materials 7-17 , such systematic searches for 2D materials are still developing [18][19][20] , especially for magnetic systems.
Nevertheless, 2D materials could be more important than 3D ones because of their unique potential as miniaturized devices and their tunability via layering or functionalization 21 .
2D topological insulating phases can be classified into two primary types: quantum spin Hall insulators (QSHI 22 ), which have time-reversal symmetry (TRS), and quantum anomalous Hall insulators (QAHI 23,24 ), which lack TRS. QSHIs, characterized by a Z2 invariant, have an insulating bulk and Dirac cone edge features. Examples include graphene, silicene 25 , germanene 26 , stanine 27 , and 1T' metal dichalcogenides 28,29 . Quantum anomalous Hall insulators (QAHI), characterized by a Z invariant known as the Chern number, are magnetic materials with a bulk gap and quantized conducting edge channels, even in the absence of an external magnetic field. Thus far, QAHI-like behavior has been observed experimentally in very few systems: Cr,V doped (Bi,Sb )2Te3 23 and MnBi2Te4 [30][31][32][33] , under highly controlled conditions. While there have been many theoretical works predicting QAHIs in 2D materials, surfaces, or interfaces, a systematic investigation of monolayer 2D materials is lacking. Some previously explored examples include transition metal halides such as CoBr2, FeCl3, NiRuCl6, V2O3, FeCl3, RuCl3 34-39 . Compared to gapped topological systems, semimetals in two-dimensions with spin-orbit coupling are relatively under-explored. For example, graphene is a well-known example of a twodimensional semimetal; however, the Dirac point in graphene is not robust to the addition of a finite amount of spin-orbit coupling 22 . This splitting of symmetry-protected band crossings by spin-orbit coupling occurs for most crossings protected by two-dimensional point group symmetries. Also, unlike in three dimensions, where Weyl points can occur at generic points in the Brillion zone, Weyl points will not occur generically in two dimensions without the tuning of some external parameter 40 . However, it is possible for non-symmorphic symmetries to protect Dirac or Weyl points in two-dimensions 41 , and we find several examples of this in our work. While 2D insulting topological materials are relatively easy to classify using topological indices such as Z2 and the Chern number, finding topological semimetals is non-trivial. Hence, having a universal strategy for screening both topological insulators and semimetals is highly desirable.
In our previous work 8 , we discovered several classes of non-magnetic topological 3D materials using spin-orbit spillage 42 to perform the initial screening step. The spillage technique is based on comparing density functional theory (DFT) wave-functions obtained with-and without spin-orbit coupling. Materials with high spillage value (discussed later) are considered non-trivial. In this paper, we extend this approach to screen the JARVIS-DFT 2D database to search for topological insulators and semimetals, with and without TRS. The JARVIS-DFT database is a part of Materials Genome Initiative (MGI) at National Institute of Standards and Technology (NIST) and contains about 40000 3D and 1000 2D materials with their DFT-computed structural, energetics 43  For metals, we search for band crossings near the Fermi level. We also predict surface (2D) and edge (1D) band-structures, Curie temperatures, and anomalous Hall conductivity. For a subset of predicted materials, we run G0W0+SOC and DFT+U calculations. We find that as we introduce improved treatments of many-body effects, only a few materials retain a non-trivial band-topology, suggesting the importance of higher-level electronic structure methods in predicting 2D topological materials. However, as an initial step, the automated spillage screening and Wannier analysis provide useful predictions of potential topological materials, narrowing down candidate materials for experimental synthesis and characterization.

Results and discussion
As mentioned above, the spin-orbit spillage criterion is applicable for both magnetic/non-magnetic and metallic/insulating classes of materials. A flow-chart describing our computational search using the spillage as well as the traditional Wannier-based method is shown in Fig. 1. Starting from 963 2D materials in the JARVIS-DFT database, we first screen for materials with OptB88vdW bandgaps < 1.5 eV, because SOC-induced band-inversion is limited to the magnitude of SOC. This leads to 506 materials. Then we carry out spin-polarized calculations with and without SOC, and compute spillage using Eq.1. Setting a threshold of 0.5, we find 122 materials with high-spillage. Note that for magnetic materials, we also screen in terms of the magnetic moment, selecting only cases with magnetic moment > 0.5 B. As a computational note, for the spin-polarized calculations, we set the initial magnetic moment to a high initial value of 6 B per atom, to search for high spin configurations. After this initial screening, we find 19 magneticinsulating, 40 magnetic-metallic, 10 non-magnetic insulating, and 53 non-magnetic metallic materials. In Table. 1, we present selected examples of each class of topological materials that we consider. The full list is given in the supplementary information. Clearly, metallic topological candidates outnumber insulators. Note that previous 2D topological material searches were mostly limited to insulators, but our approach can be extended to semimetals as well. Some of the predicted materials have already been experimentally synthesized, including Bi2TeI 50 , RuCl3 51 , FeTe 52 etc., but the experimental confirmation of their topological properties is still ongoing.
In the remaining part of the discussion, we analyze the overall distribution of topological materials, and then we focus on individual topological classes, exploring a few examples. Information for other similar materials is provided in the JAVRIS-DFT database. Finally, we discuss a few cases of DFT+U and G0W0+SOC calculations as a way to understand the limits of semi-local DFT calculations as screening tools for topological properties.

b) magnetic moment distribution for high-spillage materials, c) bandgap-distribution, d) chemical prototype distribution. Zero bandgap materials should be topological-semimetals., e) space-group distribution for high spillage materials, f) spillage of bulk (3D) vs monolayers (2D) for nonmagnetic systems.
In Figure 2 we show the spillage-based distribution for two-dimensional materials. We find that 122 materials have high spillage (with a threshold 0.5). The spillage criterion does, therefore, eliminate many materials, as shown in Fig. 2a. The 122 selected materials include both magnetic and non-magnetic materials, as well as both metals and insulators, and contains examples with a variety of chemical prototypes and crystal structures. The magnetic moments of topological materials range from 0 to 6 B, as shown in Fig. 2b. The most common topological 2D material types are AB2, ABC, AB3, and AB. Note that in experimental synthesis it might be easier to synthesize simple chemical prototypes such as AB or 1:1 structures. Hence, having a variety of chemical prototypes provides a vast amount of opportunity for synthesis. We find that most of the high-spillage materials belong to highly symmetric space groups, as shown in Fig. 2e. This is consistent with our previous 3D topological materials search. A comparison of 3D and 2D spillage in 2f indicates that the 2D-monolayer spillage values are generally lower than their 3D-bulk counterparts, but the trend is rather weak. This may be explained in part by quantum confinement in 2D-monolayers tending to increase band gaps.  Fig. 3. We observe high spillage peaks for all three materials (Fig. 3a,d,g), which occur near the Z and Gamma points. It is difficult to observe bandinversion in Fig. c because all the states near Ef are derived from Bi orbitals, so the orbitalprojected bandstructure doesn't show any obvious band-inversion. However, the spillage method suggests such inversion as in Fig. 3e. In all of these cases, we find that the bulk (Fig. 2b,e,h) is insulating while the edge states (Fig. 2c,f, (Fig 4 c, d), which carries the quantized AHC (Fig 4 d, h).

Fig. 5 Examples of non-magnetic (a,b,c,d) and magnetic (e,f,g,h) topological semimetals with spillage (a,c,e,g) and surface bandstructure (b,d,f,g) plots. a-b) Ti2TeP, c-d) AuTe2 e-f) TiCl3 g-h) MnSe
Next, we identify the non-magnetic topological semimetals. This is more challenging than analyzing insulators because there are no specific Z2/Chern-like indices for such materials. Rather

Magnetic semimetals
Finally, we look at ferromagnetic phases with band crossings, which are all Weyl crossings. We Fermi-level as shown in Fig. 5g and h. Note that the bandgaps and magnetic moments of the above systems are heavily dependent on the calculation method. This limitation of the methodology is discussed later.

Magnetic ordering
In addition to the electronic structure, we report the magnetic ordering of the structures and Curie temperature obtained using the method described in Ref. 54  but with spins in-plane, but for the remaining materials, which we predict to be FM at finite temperature, we present their Curie temperatures and Chern numbers in Table 1. We find that Tc of FeTe is unusually high because of the high J parameter resulting from the AFM configuration.
The FePS3 and FeTe have Tc higher than CrI3, which suggests they can be used for above liquid nitrogen temperature Tc topological applications. On the other hand, the Tc of ZrFeCl6 is very low, which is consistent with the large separation between Fe atoms in that crystal structure. As a next step, it would be interesting to study the layer dependence of Tc, which is beyond the scope of present work.

Limitations of GGA-based Calculations
Next, we investigate the reliability of our calculations for a few systems using GGA+U 55,56 and G0W0+SOC 57-59 methods. We note that fully ab initio beyond-DFT methods like G0W0 are very time-consuming, and hence unfeasible for a high-throughput search, while GGA+U is computationally fast, but has an adjustable parameter that limits predictive power. In Fig. 6a we show the GGA+U dependence of bandgaps and Chern numbers of ZrFeCl6. We find that as we increase the U parameter, ZrFeCl6 has a non-zero Chern number only for U < 0.2. This indicates that Chern behavior can be quite sensitive to U. Similar behavior was also observed for FeX3 (X=Cl,Br,I) by Li 35 , where the critical value of U varies from 0.43 to 0.80 eV.In Fig. 6b we show how G0W0+SOC can change the band gaps for VAgP(Se3)2, an example QAHI. In this case, G0W0+SOC increases the bandgap (0.018 eV vs 0.103 eV), but the band-shapes remain very similar, and the topology is non-trivial even in the G0W0+SOC calculation. However, for the other cases, we found that the band gap increase was substantial (Table. 2), and caused the bands to become un-inverted, resulting in trivial materials. We note that G0W0+SOC can also predict incorrect band-structures, especially for correlated transition metal compounds 60 . Here we are using only single-shot GW (G0W0), and more accurate results could be obtained by using fully self-consistent (sc)-GW 63 , which is not carried out here due to the very high computational cost.
Also, note that the bandgaps, magnetic moments can also depend on the selection of pseudopotentials. Fully accurate ab initio calculations of 2D magnetic materials remain very challenging, even for a single material, with techniques like dynamical mean field theory (DMFT) and quantum Monte-Carlo (QMC) as possible approaches for future work.

Conclusions
We have presented a comprehensive search of 2D topological materials, including both magnetic and non-magnetic materials, and considering insulating and metallic phases. Using the JARVIS-DFT 2D material dataset, we first identify materials with high spin-orbit spillage, resulting in 122 materials. Then we use Wannier-interpolation to carry-out Z2, Chern-number, anomalous Hall conductivity, Curie temperature, and surface and edge state calculations to identify topological insulators and semimetals. For a subset of predicted QAHI materials, we run G0W0+SOC and GGA+U calculations. We find that as we introduce many-body effects, only a few materials retain a non-trivial band-topology, suggesting that higher-level electronic structure methods will be necessary to fully analyze the topological properties of 2D materials. However, we believe that as an initial step, the automated spillage screening and Wannier-approach provides useful predictions for finding new topological materials.
Please note commercial software is identified to specify procedures. Such identification does not imply recommendation by National Institute of Standards and Technology (NIST). We use the OptB88vdW functional 63 , which gives accurate lattice parameters for both vdW and non-vdW (3D-bulk) solids 43 . We optimize the crystal-structures of the bulk and monolayer phases using VASP with OptB88vdW. The initial screening step for <1.5 eV bandgap materials is done with OptB88vdW bandgaps from the JARVIS-DFT database. Because SOC is not currently implemented for OptB88vdW in VASP, we carry out spin-polarized PBE and spin-orbit PBE calculations in order to calculate the spillage for each material. Such an approach has been validated by Refs. 8,64 . The crystal structure was optimized until the forces on the ions were less than 0.01 eV/Å and energy less than 10 -6 eV. We also use G0W0+SOC and DFT+U on selected structures. We use Wannier90 65 and Wannier-tools 66 to perform the Wannier-based evaluation of topological invariants.
First, out of 963 2D monolayer materials, we identify materials with bandgap < 1.5 eV, and heavy elements (atomic weight ≥ 65). We calculate the exfoliation energy of a 2D material as the difference in energy per atom for bulk and monolayer counterparts 43 . Then we use the spillage technique to quickly narrow down the number of possible materials. As introduced in Ref. 42 , we calculate the spin-orbit spillage, ( ), given by the following equation: where, is the projector onto the occupied wavefunctions without SOC, and ̃ is the same projector with SOC for band n and k-point k. We use a k-dependent occupancy ( ) of the non-spin-orbit calculation so that we can treat metals, which have varying number of occupied electrons at each k-point 8 . Here, 'Tr' denotes trace over the occupied bands. We can write the spillage equivalently as: where ( ) = 〈 |̃〉 is the overlap between occupied Bloch functions with and without SOC at the same wave vector k. If the SOC does not change the character of the occupied wavefunctions, the spillage will be near zero, while band inversion will result in a large spillage.
In the case of insulating topological materials driven by spin-orbit based band inversion, the spillage will be at least 1.0 for Chern insulators or 2.0 for TRS-invariant topological insulators at the k-point(s) where band inversion occurs 42 . In other words, the spillage can be viewed as the number of band-inverted electrons at each k-point. For topological metals and semimetals, or cases where the inclusion of SOC opens a gap, the spillage is not required to be an integer, but we find empirically that a high spillage value is an indication of a change in the band structure due to SOC that can indicate topological behavior. We choose a threshold value of = 0.5 at any k-point for our screening procedure, based on comparison with known topological semimetals 8 . Our screening method can also detect topological materials with small SOC, like the small band gap that is opened at the Dirac point in graphene 22 (JVASP-667). However, the screening method is not expected to work for semimetals with topological features that are not caused or modified by spinorbit coupling.
After spillage calculations, we run Wannier based Chern and Z2-index calculations for these materials.
The Chern number, C is calculated over the Brillouin zone, BZ, as: Here, is the Berry curvature, unk being the periodic part of the Bloch wave in the nth band, = ћꞷ , vx and vy are velocity operators. The Berry curvature as a function of k is given by: Then, the intrinsic anomalous Hall conductivity (AHC) is given by: In addition to searching for gapped phases, we also search for Dirac and Weyl semimetals by numerically searching for band crossings between the highest occupied and lowest unoccupied band, using the algorithm from WannierTools 66 . This search for crossings can be performed efficiently because it takes advantage of Wannier-based band interpolation. In an ideal case, the band crossings will be the only points at the Fermi level; however, in most cases, we find additional trivial metallic states at the Fermi level. The surface spectrum was calculated by using the Wannier functions and the iterative Green's function method 67,68 .
For magnetic systems, we primarily screen the ferromagnetic phase, with spins oriented out of the plane, which we expect is possible to achieve experimentally in most cases, using an external field if necessary. We initialize the magnetic moment with 6 for spin-polarized calculations. For a subset of interesting compounds, we perform a set of three additional calculations: ferromagnetic with spins in-plane, and antiferromagnetic with spins in-and out-of-plane. With these energies, we can fit a minimal magnetic model and calculate an estimated Curie temperature for ferromagnetic materials with out-of-plane spins, using the method of Ref. 54 , which takes into account exchange constants and magnetic anisotropy. Magnetic anisotropy that favors out-ofplane spin is crucial to enable long-range magnetic ordering in two-dimensions following the Mermin-Wagner theorem 69  To better account for correlation effects and self-interaction error in describing transition metal elements, we apply the DFT+U method to a subset of materials 55,56 . For a few systems, we also carry out a systematic U-scan by varying the U parameter from 0 to 3 eV and monitor changes in the band structure 70 . We also perform G0W0+SOC 57-59 calculations with an ENCUTGW parameter (energy cutoff for response function) of 333.3 eV for a few materials to analyze the impact of many-body effects.