A complementary screening for quantum spin Hall insulators in 2D exfoliable materials

Quantum spin Hall insulators are a class of topological materials that has been extensively studied during the past decade. One of their distinctive features is the presence of a finite band gap in the bulk and gapless, topologically protected edge states that are spin-momentum locked. These materials are characterized by a $\mathbb{Z}_2$ topological order where, in the 2D case, a single topological invariant can be even or odd for a trivial or a topological material, respectively. Thanks to their interesting properties, such as the realization of dissipationless spin currents, spin pumping and spin filtering, they are of great interest in the field of electronics, spintronics and quantum computing. In this paper we perform an high-throughput screening of quantum spin Hall insulators starting from a set of 783 2D exfoliable materials, recently identified from a systematic screening of the ICSD, COD, and MPDS databases. We find a new $\mathbb{Z}_2$ topological insulator (Hg$_4$N$_4$S$_4$) as well as 3 already known ones and 7 direct gap metal that have the potential of becoming Quantum spin Hall insulators under a reasonable external perturbation.


I. INTRODUCTION
Topological states of matter 1 have been one of the main subject of studies for condensed matter physicists in the past decade. One of the driving aspects lies in the robustness of many properties of topological materials, whose existence can be directly related to a specific set of topological invariants [2][3][4] . In the case of insulators, the ground state of a topologically trivial material can be connected to the one of an atomic insulator or vacuum by a smooth adiabatic deformation of its Hamiltonian, while for a topologically non-trivial material, this cannot happen without closing a gap in its band structure. Among the possible classifications of topological insulators, quantum spin Hall insulators (QSHIs) represent a class of materials protected by time-reversal (TR) symmetry where the quantum spin Hall effect (QSHE) 5 can be realized. This is made possible by the presence of a finite bulk gap and spin-polarized gapless helical hedge states 6 that can sustain dissipationless spin currents 1 thanks to the absence of backscattering forbidden by TR. Other interesting properties arise from spin-momentum locking 7 , which allows for the possibility to realize spin filtering 8 , injection, detection 9 and pumping 3 , all of which make QSHIs of great interest for the realization of spintronic devices. It has also been shown that magnetic confinement can be used in order to realize qubits on the surface of a QSHI 10 , making these materials promising for applications in quantum computing.
QSHIs can be identified through a Z 2 invariant which can be either even (0) or odd (1), as first proposed by Kane and Mele, that showed that a Dirac cone, with a gap opened by spin-orbit coupling (SOC) can lead to a topologically non-trivial order 5,11 . This concept has later been extended as, in general, a QSHI can be ob-tained whenever a band inversion takes place, as in the paradigmatic cases of, e.g., HgTe quantum dots 12,13 and Bi 1−x Sb x 14 . A non-trivial Z 2 invariant can emerge naturally due to the material intrinsic band structure, or can be induced by a perturbation, such as the application of strain or electric fields or the interaction with a substrate [15][16][17] . Hence, a trivial insulator can be driven into a topological insulator and vice-versa, leaving open the possibility to realize a topological transistor 18,19 . While for two-dimensional (2D) materials a single invariant is sufficient in order to describe their topology with respect to the QSHE, it has been shown that for three-dimensional (3D) materials a set of 4 invariants is required, leading to the possibility of weak or strong Z 2 topological insulators 20 . It has also been shown that the topology of a system can often be related to its symmetries 21 . For example, in the case of QSHIs with inversion symmetry, the Z 2 topological invariant can be computed easily by calculating the eigenvalues of the parity operator at the time-reversal invariant momenta (TRIM) points of the Brillouin zone (BZ) 14 .
In this work we perform a complementary screening to that of Marrazzo et al. 22 , by using a set of newly discovered 2D materials from Ref 23, in which 1205 new monolayers, exfoliable from their layered (tipically vander-Waals bonded) parents, have been added to the initial 2D portfolio 24 used in Ref 22. We set out with this, the aim of identifying all the 2D Z 2 topological insulators present in this set, analyzing all the materials with less than 40 atoms per unit cell, with the exception of lanthanides due to their typically complex, correlated electronic structure 25 . We also identify direct gap metals (DGMs) which can serve as Z 2 candidates, by calculating the invariant for materials where the valence and conduction band manifolds are separated by a direct gap (DG), but where the indirect gap (IG) is negative.

II. METHODS
All the density-functional theory (DFT) calculations performed in this work use the Quantum ESPRESSO (QE) 26,27 distribution. We perform fully relativistic calculations with SOC, using the Perdew, Burke and Ernzerhof (PBE) pseudopotentials 28 from the ONCVPSP (v0.4) set 29 contained in the PseudoDojo library 30 . For every calculations we use a kinetic cutoff, high enough to ensure that the uncertainty on the energy eigenvalues is lower than 5 meV. In order to properly account for low-gap and metallic systems, a Fermi-Dirac smearing of 5 meV is used. For the initial self consistent calculations we use a k-point mesh of density 0.2Å −1 , which is refined to 0.1Å −1 for QSHI candidates. All the calculations have been performed assuming a non-magnetic ground state, magnetism has been later checked for the selected candidates following the same approach used in Ref. 24 .
The topological invariant are computed either from an analysis of the band structures at the TRIM points or using the Z2Pack package 31 , for which we set the maximum number of kpoints per line to 200, and the lowest acceptable distance between lines to 10 −7 . The threshold on the hybrid Wannier charge centers (HWCC) 32 position variation per step on a line (pos tol) is set to 0.01, the threshold on the relative movement of the HWCCs position between neighboring (move tol) to 0.3 and the threshold for the distance between the middle of the largest gap in a line and the position of the neighboring HWCC (gap tol) to 0.3.

A. Screening procedure
The original set of material of which we explore the properties is derived from the ICSD 33,34 , COD 35 and MPDS 36 databases, by identifying 2D exfoliable materials, starting from experimentally known bulk structures using vdW DF2-C09 DFT calculations 24 . First we exclude all materials containing lanthanides, due to the inaccuracy of DFT pertaining calculations with these atoms caused by the presence of strongly correlated electrons 25 where ν is the Z 2 topological invariant, ξ 2m is the eigenvalue of the parity operator calculated for a couple of TR-paired bands and Γ i are the TRIM points, of which there are 4 for a 2D material, defined as with G being any reciprocal lattice vector. For a 2D material, 4 TRIM points can be found at the coordinates (0, 0), (0, 0.5), (0.5, 0) and (0.5, 0.5) expressed in units of the reciprocal lattice vector. Given the formula in eq.(1), the most computationally efficient route for materials that have inversion symmetry is to calculate first the band structure at the TRIM points to determine the topological invariants. Materials with ν = 0 are discarded, while for those with ν = 1 we proceed with a band structure calculation along the high-symmetry lines. The result can either be a material with finite DG on the entire k-point path, or with a zero DG. The latter is discarded as a metal/semimetal, while the former can either be a QSHI if the IG, defined as the difference between the minimum of the conduction band and the maximum of the valence band, is positive, or a DGM if the IG is negative. We discard materials that are DGM with an IG lower than -0.15 eV, as attempting to drive such a material to a topological phase would be hardly feasible.
If inversion symmetry is not present, we follow a different route, by determining the DG and IG first through a band structure calculation and discarding the metals and semimetals without a finite DG. We then use the Z2Pack package, which tracks the positions of one-dimensional HWCCs across half of the BZ, in order to compute the topological invariant 37 . The package acts as an automation tool, which interfaces a DFT code (QE in our case) with Wannier90 38 in order to calculate the HWCC positions on a progressively denser grid of k-points. If one of the convergence criteria of Z2Pack, with the parameters mentioned before, is not satisfied, then the material is discarded. If the results of the calculation is a non trivial topological invariant, the material is then classified as a QSHI or DGM following the same criteria mentioned for the inversion-symmetric cases.
Finally, for the selected candidates we asses the dynamical stability by computing the full phonon dispersion by means of density-functional perturbation theory 39 and we check the non-magnetic ground state assumption comparing the energy of the non-magnetic ground state with different possible ferromagnetic, antiferromagnetic and ferrimagnetic configurations following the same procedure adopted in Ref. 24. A flowchart describing the entire high-throughput calculation is shown in Fig. 1. The entire process has been automated using AiiDA 40,41 , a workflow managing infrastructure, which enable us to keep track of the provenance at every computational step. The data produced in this work, including both the results and the entire provenance tree, has been made available on the Materials Cloud as Ref. 42.

III. RESULTS AND DISCUSSIONS
We start our screening from a set of 783 2D exfoliable structures taken from Ref 23, which are reduced to 603 after removing materials containing lanthanides. The set is further divided in two groups of 260/342 structures with/without inversion symmetry. For the inversion symmetric materials, we compute the eigenvalues of the parity operator ξ 2m , the calculation of which can sometimes fail due to limitations in the code to handle non-symmorphic space groups, or due to the accuracy threshold imposed for the calculation of the trace of the representation for each group of bands. 78 out of these 260 TRIM point calculation failed and were hence recalculated using Z2pack. Of the 182 successful calculations, 27 resulted in a non trivial topological invariant ν = 1. We then performed a band structure calculations along the high-symmetry path 43 , from which we excluded 22 metallic materials and identified 1 QSHI (Br 2 Hf 2 ) and 4 direct gap metals (DGMs). On the 343+78 remaining structures (the noninversion symmetric ones and those with failed parity calculations), the band structure calculations show that 207 materials have a DG greater than 0; for them we proceed with Z2Pack calculations that identify, 3 QSHIs (Zr 2 Te 10 , Hg 4 N 4 S 4 , Hf 2 Te 10 ) and 3 DGMs. The summary of the data for every material is reported in Table I, while a more detailed report is given in the supplementary materials, together with the crystal structure and band structure of the material identified and a plot of the HWCC for every Z2Pack calculation.
A summary for all the materials analyzed in the present work, together with those of the previous screening by Marrazzo et al. 22 (in absence of strain), is shown in Fig. 2, where the QSHIs are highlighted. In the figure, all materials with a negative IG have been put on the y = 0 axis in order to emphasize QSHIs and DGMs, especially with regard to the magnitude of the DG and IG, as these two parameters can be used to determine the range of temperatures at which a device could operate 44 .
The interest in DGMs with a non trivial ν is twofold. First it is possible that, given the limitations of DFT-PBE calculations, the band gap of the material is being underestimated and, using more appropriate methods (such as many-body perturbation theory in the GW approximation 45 ), a material that was computed to be a DGM with DFT-PBE could actually result to be a QSHI. It is worth highlighting that the opposite could also be possible, where a material that is estimated to be a QSHI within DFT-PBE could end up being a trivial insulator in a GW calculation (as discussed for the case of TiNI in Ref. 22 ). Second, a DGM could be driven into a Z 2 topological state through an external perturbation that could open a gap in the material without driving a band inversion. To this end, it is important to look at both the IG and DG of the material. The former is related to how strong a perturbation needs to be in order to drive the material in a semiconducting state, while the latter gives us an idea of how likely it is that a gap could be opened, without causing band inversion, or driving the material into a metallic state. Indeed, the higher the DG of a material, the more difficult it will be to close it. For this reason, DGMs with small negative IG and large DG are prime candidates for further studies; a visual summary of the direct and indirect gap of the DGMs is shown in Fig. 3. Among the various possible perturbations that could open the gap, the most interesting ones, from the standpoint of applications in a device, would be gating/electric fields 16 , the interaction with a substrate 46 and/or the application of strain 47 .
For all the materials found, we searched the current literature to determine which QSHIs were already known, and which one are novel to this study. We find that among the 4 QSHIs identified, 3 were already known, namely the two tellurides of hafnium/zirconium, and hafnium bromide. Br 2 Hf 2 48,49 was already known both for being an easily exfoliable material and for its topological properties in its 2D and 3D forms. This material has a DG of 233 meV and an IG of 48 meV making it a good candidate for experiments and possible device realization. It should also be noted that this material be- Hf 2 Te 10 (Zr 2 Te 10 ) 52 are two pentatellurides both of which have large DGs of 300(281) meV and IGs of 171(220) meV making them both very good candidates for QSHIs applications at room temperature. Interestingly, the inversion in these materials has been shown to be due to the non-symmorphic space group and not to SOC, the effect of which is only to open a gap in a band structure that would be metallic in a non "fullyrelativistic" calculation. These pentatellurides have also been proposed for the realization of the quantum anomalous Hall effect, when the transition metal Hf/Zr is substituted for a rare earth metal which could induce a magnetic ordering in the material, breaking time-reversal symmetry 52,53 .
Hg 4 N 4 S 4 is a new compound which has not been discussed in the literature for its topological properties. Its direct and indirect gaps are both of 7 meV. Such low gap, in between that of silicene (1.4 meV) and germanene (23.8 meV) precludes room temperature applications. Nevertheless, at the gap, situated along the Γ − Y direction, the material shows a linear band dispersion akin to that of germanene at K when SOC is taken into account. The similarity between the two systems, opens the possibility for a plethora of studies, such as the effect of an electric field on the gap and band inversion 16,17 or the possibilities for opto-electronic applications in THz devices 54  metal halaide and has been studied in Ref 48, in which the material is identified as a semi-metal for both PBE and PBE-SOC calculations, but as QSHI in a HSE06-SOC calculation. This reiterates that DGMs identified by DFT can or could actually be found to be QSHIs if studied with techniques, such as many body perturbation theory, that can correctly estimate the gap, or when exposed to a perturbation.

IV. CONCLUSIONS
We performed a systematic high-throughput study of 2D exfoliable materials in order identify novel QSHIs.
Starting from a set of 783 materials we found 4 QSHIs and 7 DGMs that could be driven into a QSHI state by material engineering. Of the 4 QSHIs identified, the 3 with the largest gaps had already been explored in the literature. In particular, the two pentatellurides Hf 2 Te 10 and Zr 2 Te 10 exhibit both large direct and indirect gaps which makes them promising candidates for room-temperature applications. The remaining material, namely Hg 4 N 4 S 4 , represents a newly identified QSHI with linear bands in proximity of its direct gap of 7 meV. The main feature of this material is the small gap and band linearity, akin to that of silicene or germanene. For this reason, Hg 4 N 4 S 4 could be a promising candidate for future electronic and photonic applications in the low energy regime.
All materials identified in this work are easily exfoliable, with a binding energy lower than 20 meV/Å 224 and for this reason they are prime candidates for possible experimental realization through exfoliation techniques starting from their experimentally known 3D bulk form.

V. ACKNOWLEDGMENTS
This work was supported by the National Centre for Computational Design and Discovery on Novel Materials (NCCR MARVEL) of the Swiss National Science Foundation. D.G. gratefully acknowledge support from the EU Centre of Excellence, MaX Materials design at the eXascale (Grant 824143). We acknowledge PRACE for awarding us access to Marconi at Cineca, Italy (project id. 2016163963).