Discovering new two-dimensional topological insulators from computational screening

TL;DR

This study computationally screened 3331 2D materials to identify new topological insulators, revealing 46 quantum spin Hall, 7 quantum anomalous Hall, and 9 crystalline topological insulators, with detailed analysis on their stability and electronic properties.

Contribution

Introduced a comprehensive b initio scheme to identify topological phases in 2D materials from a large database, discovering several novel topological insulators.

Findings

Identified 46 quantum spin Hall insulators, 7 quantum anomalous Hall insulators, and 9 crystalline topological insulators.

Found that topological indices are highly sensitive to exchange-correlation functional approximations.

Calculated a 0.65 eV gap for PdSe, larger than known 2D topological insulators.

Abstract

We have performed a computational screening of topological two-dimensional (2D) materials from the Computational 2D Materials Database (C2DB) employing density functional theory. A full \textit{ab initio} scheme for calculating hybrid Wannier functions directly from the Kohn-Sham orbitals has been implemented and the method was used to extract $\mathbb{Z}_2$ indices, Chern numbers and Mirror Chern numbers of 3331 2D systems including both experimentally known and hypothetical 2D materials. We have found a total of 46 quantum spin Hall insulators, 7 quantum anomalous Hall insulators and 9 crystalline topological insulators that are all predicted to be dynamically stable. Roughly one third of these were known prior to the screening. The most interesting of the novel topological insulators are investigated in more detail. We show that the calculated topological indices of the quantum anomalous Hall insulators are highly sensitive to the approximation used for the exchange-correlation functional and reliable predictions of the topological properties of these materials thus require methods beyond density functional theory. We also performed $GW$ calculations, which yield a gap of 0.65 eV for the quantum spin Hall insulator PdSe$_2$ in the MoS$_2$ crystal structure. This is significantly higher than any known 2D topological insulator and three times larger than the Kohn-Sham gap.

Figures

Tables

Table 1. Table 1: Overview of known topological materials found by computational screening. The topology is specified by either the Chern number C 𝐶 C , the ℤ 2 subscript ℤ 2 \mathbb{Z}_{2} index ν 𝜈 \nu or the mirror Chern number C M subscript 𝐶 𝑀 C_{M} . We also state the calculated Kohn-Sham gap (KS gap).

Material	Prototype	Topology	KS gap [meV]
C2Kane2005a	C2	$\nu =1$, $C_M=1$	0.3
Si2Ezawa2015	C2	$\nu =1$	1.6
Ge2Ezawa2015	C2	$\nu =1$	25
Sn2Ezawa2015	C2	$\nu =1$	65
SnFXu2013	CH	$\nu =1$	316
HgSeLi2015	GeSe	$\nu =1$	90
HgTeLi2015	GeSe	$\nu =1$	156
MoS2Qian2014	WTe2	$\nu =1$	51
MoSe2Qian2014	WTe2	$\nu =1$	41
WSe2Qian2014	WTe2	$\nu =1$	32
OsCl3Sheng2017	BiI3	$C=1$	64
GeTe2Zhang2018a	MnS2	$\nu =1$	32
SnSLiu2015c	PbS	$C_M=2$	67
SnSeLiu2015c	PbS	$C_M=2$	84
SnTeLiu2015c	PbS	$C_M=2$	28
PbSLiu2015c	PbS	$C_M=2$	422
PbSeLiu2015c	PbS	$C_M=2$	478
PbTeLiu2015c	PbS	$C_M=2$	271

Table 2. Table 2: Overview of new QSHIs without mirror symmetry. All the materials have a ℤ 2 subscript ℤ 2 \mathbb{Z}_{2} index of ν = 1 𝜈 1 \nu=1 . We also state the calculated Kohn-Sham gap (KS gap), the heat of formation (HOF) and the energy above the convex hull (EACH). The dynamically stable materials are shown in bold face.

Material	Prototype	Topology	KS gap [meV]	HOF [eV]	EACH [eV]
AuCl	FeSe	$\nu =1$	20	0.10	0.30
CrAsBi	BiTeI	$\nu =1$	35	0.31	0.46
IrSe	GaSe	$\nu =1$	134	0.18	0.45
TiTe	GaSe	$\nu =1$	109	-0.08	0.63
ZrTe	GaSe	$\nu =1$	207	-0.28	0.65
AuI3	BiI3	$\nu =1$	109	0.10	0.10
TiIN	FeOCl	$\nu =1$	62	-1.18	-0.26
TlClSe	FeOCl	$\nu =1$	27	-0.36	0.32
TiS	CuI	$\nu =1$	54	-1.13	0.31
TiCl	CuI	$\nu =1$	13	-0.64	0.45
ZrS	CuI	$\nu =1$	132	-1.16	0.26
ZrSe	CuI	$\nu =1$	20	-0.91	0.25
ZrCl	CuI	$\nu =1$	37	-0.59	0.73
ZrBr	CuI	$\nu =1$	45	-0.34	0.68
SbCl	CuI	$\nu =1$	434	-0.46	0.13
SbBr	CuI	$\nu =1$	442	-0.31	0.10
SbI	CuI	$\nu =1$	584	-0.12	0.13
HfS	CuI	$\nu =1$	158	-0.89	0.43
HfSe	CuI	$\nu =1$	42	-0.64	0.54
ReS	CuI	$\nu =1$	309	0.10	0.54
HgCl	CuI	$\nu =1$	129	-0.37	0.23
HgBr	CuI	$\nu =1$	188	-0.25	0.22
PbF	CuI	$\nu =1$	116	-1.43	0.45

Table 3. Table 3: Overview of new crystalline topological insulators protected by mirror symmetry. The topology is specified the mirror Chern number C M subscript 𝐶 𝑀 C_{M} . All the materials are invariant under time-reversal symmetry and the associated ℤ 2 subscript ℤ 2 \mathbb{Z}_{2} index thus becomes ν = 1 𝜈 1 \nu=1 if C M subscript 𝐶 𝑀 C_{M} is odd and ν = 0 𝜈 0 \nu=0 if C M subscript 𝐶 𝑀 C_{M} is even. We also state the calculated Kohn-Sham gap (KS gap), the heat of formation (HOF) and the energy above the convex hull (EACH). The dynamically stable materials are shown in bold face.

Material	Prototype	Topology	KS gap [meV]	HOF [eV]	EACH [eV]
HgO	BN	$\nu =1$, $C_M=1$	302	-0.15	0.21
PdSe2	MoS2	$\nu =1,C_M=1$	229	-0.02	0.27
AuTe	GaS	$\nu =1,C_M=1$	37	0.03	0.11
WO	GaS	$\nu =1,C_M=1$	53	-1.01	0.49
RhO	GaS	$\nu =1,C_M=3$	67	-0.38	0.25
IrO	GaS	$\nu =1,C_M=3$	122	-0.16	0.51
ReI3	AgBr3	$\nu =1,C_M=1$	141	0.37	0.43
ReCl3	TiCl33	$\nu =1,C_M=1$	220	-0.18	0.52
WI3	TiCl33	$\nu =1,C_M=1$	222	0.19	0.19
GeTe2	MnS2	$\nu =1,C_M=1$	32	0.27	0.33
RuTe2	MnS2	$\nu =1,C_M=1$	157	0.34	0.68
OsS2	MnS2	$\nu =0,C_M=2$	111	0.25	0.63
OsSe2	MnS2	$\nu =0,C_M=2$	144	0.51	0.68
OsTe2	MnS2	$\nu =0,C_M=2$	117	0.70	0.77

Table 4. Table 4: Overview of new QAHIs found by computational screening. The topology is specified by either the Chern number C 𝐶 C . We also state the calculated Kohn-Sham gap (KS gap), the heat of formation (HOF) and the energy above the convex hull (EACH). The dynamically stable materials are shown in bold face.

Material	Prototype	Topology	KS gap [meV]	HOF [eV]	EACH [eV]
PdF2	GeS2	$C=-2$	27	-0.83	0.38
FeCl3	BiI3	$C=1$	13	-0.66	-0.08
FeBr3	BiI3	$C=1$	42	-0.38	-0.04
PdI3	BiI3	$C=-1$	36	-0.13	0.06
CoBr3	AgBr3	$C=-2$	27	-0.21	-0.16
MoS2	MnS2	$C=2$	133	0.10	1.05

Equations

$$\mathbf{k}=k_1{\mathbf{b}}_1+k_2{\mathbf{b}}_2,$$

$$|W_{nj{k}_2}⟩={\int }_0^{1}d{k}_1{e}^{-i\mathbf{k}\cdot (\hat{\mathbf{r}}+{\hat{\mathbf{R}}}_j)}|{\tilde{u}}_{n\mathbf{k}}⟩.$$

$$x_{1,n}(k_2)$$

$${\gamma }_{1,n}(k_2)=i{\int }_0^{1}d{k}_1⟨{\tilde{u}}_{n\mathbf{k}}|{\partial }_{k_1}{\tilde{u}}_{n\mathbf{k}}⟩.$$

$$S_n^{(z)}(k_2)\equiv ⟨W_{n0k_2}|{\hat{S}}^{(z)}|W_{n0k_2}⟩,$$

$${\sigma }_{xy}=-\frac{e^2}{h}{\int }_{BZ}\frac{d^{2}k}{2\pi }\Omega (\mathbf{k}),$$

$$\Omega (\mathbf{k})=i\sum _{ijn}{f}_{n\mathbf{k}}{\epsilon }_{ij}{\partial }_{k_i}⟨u_{n\mathbf{k}}|{\partial }_{k_j}{u}_{n\mathbf{k}}⟩.$$

$${\sigma }_{xy}=-C\frac{e^2}{h},C\in \mathbb{Z}.$$

$$M_{mn}(\mathbf{k},\mathbf{k}+\delta \mathbf{k})$$

$$=⟨{\psi }_m(\mathbf{k})|e^{-i\delta \mathbf{k}\cdot \hat{\mathbf{r}}}|{\psi }_n(\mathbf{k}+\delta \mathbf{k})⟩,$$

$$\displaystyle|\psi_{n}\rangle=|\tilde{\psi}_{n}\rangle+\sum_{ai}\langle\tilde{p}_{i}^{a}|\tilde{\psi}_{n}\rangle\Big{[}|\phi^{a}_{i}\rangle-|\tilde{\phi}^{a}_{i}\rangle\Big{]},$$

$$M_{mn}(\mathbf{k},\mathbf{k}+\delta \mathbf{k})$$

$$+\sum _{aij}{e}^{-i\delta \mathbf{k}\cdot {\mathbf{r}}_a}⟨{\tilde{\psi }}_m(\mathbf{k})|{\tilde{p}}_i^a⟩$$

$$\displaystyle\times\Big{[}\langle\phi_{i}^{a}|\phi_{j}^{a}\rangle-\langle\tilde{\phi}_{i}^{a}|\tilde{\phi}_{j}^{a}\rangle\Big{]}\langle\tilde{p}_{j}^{a}|\tilde{\psi}_{n}(\mathbf{k}+\delta\mathbf{k})\rangle,$$

$$H=-\frac{1}{2}J\sum _{⟨ij⟩}{\mathbf{S}}_i\cdot {\mathbf{S}}_i-A\sum _i(S_j^z{)}^2-\frac{1}{2}B\sum _{⟨ij⟩}{S}_i^z{S}_j^z,$$

Full text

I Introduction

The concept of topological band theory was initially developed in order to explain the quantum Hall effect, which was observed experimentally in 1980.Klitzing The measurements were soon interpreted as a topological effect arising from the phases of Bloch states winding around the boundary of the magnetic Brillouin zoneThouless1982 and is thus closely related to the $k$-space Berry phase.Berry1984a In 1988, Haldane proposed a model system that exhibited the quantum Hall effect without an external magnetic field, but with intrinsically broken time-reversal symmetry.Haldane1988 Such materials are referred to as quantum anomalous Hall insulators (QAHI) and the first experimental demonstration of the effect was reported in 2013 - 25 years after it was proposed. In the mean time, Kane and Mele had showed that any time-reversal invariant 2D insulator can be characterized by a $\mathbb{Z}_{2}$ topological index $\nu$. Time-reversal invariant materials with non-trivial topology ($\nu=1$) are known as quantum spin Hall insulators (QSHI) and the effect was observed immediately after a theoretical prediction of the effect in HgTe quantum wells.Mosfets2013 ; Konig2007 Subsequently, the concepts have been generalized to bulk 3D systemsQi2008 ; Teo2008 and several well-known materials have been shown to comprise examples of topological insulators.Hasan2010 Most notably, Sb2Te3, Bi2Se3 and Bi2Te3,Hsieh2009 ; Zhang2009 but also several 2D materials have been shown to exhibit a non-trivial band topology. In fact, graphene comprised the first theoretical prediction for a QSHI and while it is still believed that graphene has a non-trivial band topology it is practically impossible to verify experimentally due to the small band gap. However, several other 2D materials have been shown to comprise examples of QSHIs. For example, the graphene-like materials silicene,Ezawa2012 ; Ezawa2013 germanene,Zhang2016a ; Amlaki2016 and staneneTang2014 are all predicted to be QSHIsMolle2017 and so are several of the transition metal dichalcogenides in the 1T’ phase.Qian2014 ; Olsen2016a

While the quantum spin Hall effect has been observed in a wide range of both 2D and 3D materials, the quantum anomalous Hall effect has proven more elusive and has so far only been observed in a few $\mathbb{Z}_{2}$ topological insulators, where time-reversal symmetry is broken by introducing magnetic impurities.Chang2013a ; Mogi2017a There has been a few proposals for pristine 2D materials that are predicted to be intrinsic QAHIs by first principles calculations,Tse2011 ; Garrity2014 ; Sheng2017 ; Wu2017 ; Chen2017a but the effect has not yet been confirmed experimentally for any of the materials and the topological properties seem to be somewhat sensitive to the details of the calculations. Moreover, first principles calculations typically only pertain to the case of zero Kelvin, but for 2D materials magnetic order is highly fragile to finite temperature effects and can only be stabilized in the presence of magnetic anisotropy.Mermin ; Torelli2018 Realistic theoretical predictions of novel 2D QAHIs thus have to take into account that the magnetic order must persist at experimentally relevant temperatures - preferably room temperature.

Since the discovery of topological classifications of solids, the field has witnessed a tremendous development of the theoretical concepts, which have been extended to include topological semi-metals,Armitage2018 topological crystalline insulatorsFu2011 ; Hsieh2012 and higher order topological effects.Benalcazar2017a ; Benalcazar2017b ; Schindler2018a ; Schindler2018 In addition to QAHIs and QSHIs we will focus on a particular class of topological insulators in the present work, namely, the topological crystalline insulators where the topology is protected by mirror symmetry in a plane parallel to the 2D material.Liu2015c ; Niu2015 It is particularly easy to understand the topology in this case since the topological index is simply defined as the difference between the quantum anomalous Hall conductance of the two eigenspaces of the Mirror operator. Moreover, the fact that the mirror plane coincides with the plane of the materials implies that any edge that conserves the mirror symmetry will host gapless states that are protected from back-scattering by the topology.

Although topological properties are in principle derivable from the ground state of bulk materials, the observable consequences are very limited in bulk. It is, for example, not possible to distinguish a non-trivial topological state from a trivial one by looking at the band structure alone. However, any interface between materials belonging to different topological classes is guaranteed to host gapless states that are localized at the boundary. Since the boundary states are protected by topology, they are extremely stable and has been proposed as candidates for novel dissipationless electronics circuits. Moreover, since the spin of boundary states are locked to the direction of propagation in QSHIs, such materials are promising for spintronics applications.Pesin2012

A crucial ingredient for operating topological insulators at room temperature is a sizable band gap. Typically, the band gap in topological insulators is determined by the strength of spin-orbit coupling and common values of the gap are on the order of 0.1-0.2 eV. With such small band gaps it becomes hard to maintain full control of the gapless boundary states and the applicability of the materials becomes questionable. It would thus be highly desirable to find new topological insulators with large band gaps. Large-scale screening studies based on first principles computations have previously been performed for 3D materialsYang2012b ; Zhang2018b and very recently for 2D materials based on experimentally known van der Waals bonded crystals.Mounet2018 The latter study predicted a novel 2D material in the Jacutingaite prototype to be a QSHI with a band gap of 0.5 eV.Marrazzo2017a Remarkably, the PBE gap of the material was only 0.15 eV, but a quasi-particle gap 0.5 eV was obtained with the $G_{0}W_{0}$ approximation. It is thus far from obvious that the simple Kohn-Sham gap provides a good estimate of the quasi-particle gap in 2D topological insulators.

Finally, we emphasize that it has not yet been possible to demonstrate the quantum anomalous Hall effect in a pristine 2D material. This is perhaps not so surprising since a magnetically ordered ground state is a minimal requirement for the effect and until very recently magnetic order had not been observed in 2D.Huang2017a The discovery of the first pristine 2D QAHI, which exhibits magnetic order at reasonable temperatures thus comprises a tremendous challenge and it is highly likely that theoretical predictions may aid this quest by significantly decreasing the number of relevant materials to investigate.

In the present work we have screened more than 3000 hypothetical 2D materials using first principles simulations and identified novel as well as well-known 2D topological insulators. We have focused on three topological classes: 1) QAHIs, which requires a magnetically ordered ground state and are characterized by the Chern number $C$ that may take any integer value. 2) QSHIs, which require time-reversal symmetry and are classified by the binary $\mathbb{Z}_{2}$ index. 3) Mirror crystalline topological insulators, where the topology is protected by mirror symmetry and the ground state is classified according to the mirror Chern number $C_{M}$, which may take any integer value. In order to extract the topological properties we have implemented a full calculation of $k$-space Berry phases that allow us to extract the topological indices in a semi-automated way and does not depend on a mapping to tight binding models through Wannier functions.Gresch2017

The paper is organized as follows. In Sec. II , we describe the Berry-phase implementation and exemplify how the topological indices are are extracted for the three cases described above. In Sec. III, we present the computational details and provide a comprehensive list of all the 2D materials with non-trivial topology that have been found in the screening. We then analyze the topological properties of a few representative materials in more details and investigate the effect of the approximations for exchange-correlation energy and $G_{0}W_{0}$ calculations. In Sec. IV we provide a discussion of the results.

II Berry phases and Hybrid Wannier functions

In this section we will briefly introduce the notion of parallel transport and show how it can be applied to obtain the Berry phase matrix of a closed path in $k$-space. We will closely follow the discussions in Refs. Marzari1997, , Taherinejad2014, , and Gresch2017,

We consider a minimal unit cell in reciprocal space spanned by the reciprocal lattice vectors $\mathbf{b}_{1}$ and $\mathbf{b}_{2}$. A generic point in the reciprocal unit cell can then be written as

[TABLE]

where $0\leq k_{i}<1$ are the reciprocal fractional coordinates. We wish to calculate the Berry phase obtained by transporting the occupied Bloch states along $k_{1}$ through the reciprocal unit cell at a fixed value of $k_{2}$ (See Fig. 1). In a numerical treatment of the Bloch Hamiltonian $H(\mathbf{k})$, one obtains a set of occupied eigenstates $|u_{n\mathbf{k}}\rangle$ at different $\mathbf{k}$-points and each set of eigenstates come with an arbitrary set of phases. In fact, any unitary rotation in the space of occupied states leaves the ground state energy invariant and in order to evaluate the phases picked up along a closed path in $k$-space one needs to construct eigenstates with phases that are smooth along the path. This can be accomplished by the so-called parallel transport gauge. For a single occupied state we fix the phase along the path by requiring that $\langle u_{\mathbf{k}}|\partial_{k_{1}}u_{\mathbf{k}}\rangle=0$, which enforces that the change in the state along the path is orthogonal to the state itself. When the Bloch states are calculated at a string of $N$ $\mathbf{k}$-points with fixed $k_{2}$ and $k_{1}=0,1/N,2/N,\ldots$, the parallel transport condition can be implemented by requiring that $\langle u_{\mathbf{k}}|u_{\mathbf{k}+\mathbf{b}_{1}/N}\rangle$ is real at any point along the path. Thus for an initial state at $\mathbf{k}$ on the path we can fix all phases on the path sequentially by the parallel transport condition, which will result in a smooth phase along the path. The phase at $k_{1}=1$ can then be obtained by imposing the periodic gauge such that $|u_{\mathbf{k}+\mathbf{b}_{1}}\rangle=e^{-i\mathbf{\hat{r}}\cdot\mathbf{b}_{1}}|u_{\mathbf{k}}\rangle$ and the Berry phase is obtained as the phase difference between $k_{1}=0$ and $k_{1}=1$.

In the case of multiple occupied bands, the condition is generalized by requiring that the matrix $M_{mn\mathbf{k}}=\langle u_{m\mathbf{k}}|u_{n\mathbf{k}+\mathbf{b}_{1}/N}\rangle$ is Hermitian. This uniquely fixes the unitary rotation among the occupied states at all points along the path in terms of an initial set of occupied states at particular point $\mathbf{k}$ on the path. This is due to the fact that a single value decomposition yields $M=V\Sigma W^{\dagger}=V\Sigma V^{\dagger}VW^{\dagger}$, where $V$ and $W$ are unitary matrices and $\Sigma$ is diagonal. We can thus take $|u_{n\mathbf{k}+\mathbf{b}_{1}/N}\rangle\rightarrow WV^{\dagger}|u_{n\mathbf{k}+\mathbf{b}_{1}/N}\rangle$, which renders $M$ Hermitian and completely fixes the gauge at $\mathbf{k}+\mathbf{b}_{1}/N$. This procedure is continued along the path until the states at $\mathbf{k}+\mathbf{b}_{1}$ are obtained using the periodic gauge where $|u_{n\mathbf{k}+\mathbf{b}_{1}}\rangle=e^{-i\mathbf{\hat{r}}\cdot\mathbf{b}_{1}}|u_{n\mathbf{k}}\rangle$. Finally, the eigenvalues of the unitary matrix relating the states at $\mathbf{k}$ and $\mathbf{k}+\mathbf{b}_{1}$ are the Berry phases acquired by the individual bands.

The method can also be used to obtain the individual states $|\tilde{u}_{n\mathbf{k}}\rangle$ that are parallel transported without mixing and thus acquire the distinct eigenvalues of the Berry phase matrix.Marzari1997 Since these are smooth, one may construct hybrid Wannier functions (HWFs) localized along the direction parallel to $\mathbf{a}_{1}$ as

[TABLE]

Writing $\mathbf{r}=x_{1}\mathbf{a}_{1}+x_{2}\mathbf{a}_{2}$ with $\mathbf{a}_{i}\cdot\mathbf{b}_{j}=2\pi\delta_{ij}$, one can show thatMarzari2012

[TABLE]

with

[TABLE]

Except for a factor of $2\pi$, the Berry phases obtained from the parallel transport gauge are thus the charge centers of the HWFs. This construction allows one to calculate various properties associated with the individual phases. For example the spin expectation value

[TABLE]

which will serve as a useful tool for analyzing the topological properties of 2D materials below.

The Berry phases will be smooth functions of $k_{2}$, which implies that one can track the evolution of the phases while $k_{2}$ is cycled through the reciprocal space unit cell. The dispersion of the spectrum of Berry phases gives rise to topological classifications as explained below. We will provide examples of ab initio calculations of the Berry phase spectrum for the three topological classes considered in the present work, but postpone a compilation of the computational details until Sec. III.

II.1 Quantum Anomalous Hall Insulators

The Hall conductance $\sigma_{xy}$ relates an electric current in the $x$-direction to a uniform field in the $y$-direction by $J_{x}=\sigma_{xy}E_{y}$. A finite Hall conductance requires broken time-reversal symmetry and it follows from the Kubo formula that it can be written in terms of the $k$-space Berry curvature $\Omega_{z}(\mathbf{k})$ asNagaosa2010

[TABLE]

with

[TABLE]

Here $i,j$ runs over $x$ and $y$, $\varepsilon$ is the two-dimensional Levi-Civita symbol, and $f_{n\mathbf{k}}$ are occupation factors.

For insulators, the integral can be shown to yield an integer known as the Chern number $C$ and the Hall conductance becomes

[TABLE]

In 2D metals, a gap can be opened due to the Landau levels emerging when an external magnetic field is introduced and the value of the Chern number can be controlled by the magnitude of magnetic field. That is the quantum Hall effect. Moreover, as shown by HaldaneHaldane1988 materials with spontaneously broken time-reversal symmetry can exhibit intrinsic quantum Hall effect without an external magnetic field. Such materials are known as quantum anomalous Hall insulators (QAHI) (or Chern insulators) and have a finite Chern number and a non-trivial band topology.

We now briefly discuss how the $k$-space Berry phase calculations outlined above can be related to the Chern number in QAHIs. A constant electric field in the $y$-direction of magnitude $E_{0}$ can be included in the Bloch Hamiltonian by the substitution $k_{y}\rightarrow k_{y}-eE_{0}t/\hbar$. Clearly the physical properties of the Bloch Hamiltonian are restored after a period $T=2\pi\hbar/eE_{0}$. If $E_{0}$ is sufficiently small, the system will evolve adiabatically and the charge transported in the $x$-direction in the time interval is $Q=\int_{0}^{T}J_{x}dt=T\sigma_{xy}E_{0}=h\sigma_{xy}/e$. Expressing the Hall conductance as $\sigma_{xy}=-Ce^{2}/h$ we see that the transferred charge is $Q=-eC$. One can thus obtain the Chern number as the number of Wannier charge centers that are transported by a unit cell in the $x$-direction, while $k_{y}$ is cycled through the reciprocal space unit cell in the backward direction. This argument comprises a modified versionQi2008 of the explanation originally provided by Laughlin to account for the quantum Hall effect.Laughlin1981 This also explains the appearance of gapless edge states in QAHIs. For a bulk system an adiabatic cycling of $k_{2}$ through a reciprocal unit vector will return the system to itself by transferring $C$ charges by one lattice vector. However, in the presence of an edge charges will pile up at the edge and the argument breaks down unless there is a different mechanism that can remove the charges from the edge. One may thus conclude that any edge has to host chiral gapless states that connect the valence band with the conduction band, such that $C$ units of charge are transferred out of the valence bands at the edge while $k_{2}$ is cycled by a unit.

In Fig. 2 we show the Berry phases of FeBr3 in the BiI3 crystal structure. Due to Eq. (3), the vertical axis in the plot can be regarded as the unit cell in the $x$-direction and it is clear that at any horizontal line there will be a total of one state crossing in the downward direction. We can thus calculate the Chern number as the total number of chiral crossings - that is, the number of crossing points with negative slope minus the number of crossing points with positive slope at any horizontal line and we conclude that $C=1$. We have confirmed this result by a direct integration of the Berry curvature $\Omega_{z}$ over the Brillouin zone.

II.2 Quantum Spin Hall Insulators

Materials with time-reversal symmetry must have a vanishing Hall conductance and thus cannot have a non-zero Chern number. From the Berry phase perspective this is due to the fact that any phase at $\mathbf{k}$ will be accompanied by an equal phase at $-\mathbf{k}$, which excludes the possibility of having a finite number of chiral crossings at a horizontal line. Nevertheless, as shown by Kane and MeleKane2005 ; Kane2005a all 2D time-reversal invariant insulators belong to one of two topological classes and can thus be characterized by a $\mathbb{Z}_{2}$ topological index. The simplest way to understand the $\mathbb{Z}_{2}$ index is by considering a system where $S_{z}$ is a good quantum number. Then one can calculate the Chern numbers of the two spin states and obtain $C_{\uparrow}=-C_{\downarrow}$, since $C=C_{\uparrow}+C_{\downarrow}=0$. Due to time-reversal symmetry all Berry phases have Kramers degenerate partners at time-reversal invariant points such as $\mathrm{M}$ and $\Gamma$. If one introduces spin-mixing (but time-reversal conserving) perturbations in the Hamiltonian, the two Chern numbers are no longer well-defined but two distinct ways of connecting phases between time-reversal invariant points remains, which gives rise to a $\mathbb{Z}_{2}$ topological classification. In particular if $C_{\uparrow}$ is even in the case where $S_{z}$ is a good quantum number, there will be an even number of Kramers pairs at $\Gamma$ and $M$ that may to hybridize and open a gap in the Berry phase spectrum, once spin-mixing terms are introduced. In contrast, if $C_{\uparrow}$ is odd, there will be an odd number of Kramers pairs at $\mathrm{M}$ and $\Gamma$ and the Berry phase spectrum must remain gapless when spin-mixing terms are included. In general the $\mathbb{Z}_{2}$ index distinguishes whether there is an odd or even number of Berry phases crossing any horizontal line in half the Brillouin zone.

In Fig. 3 we show the band structure and Berry phase spectrum for stanene (Sn), which comprises an example of a quantum spin Hall insulator. Due to inversion symmetry in stanene, all bands are doubly degenerate and cannot be colored according to spin. However, the degeneracy of the Berry phases are split due to spin-orbit coupling, which will allow the spectrum to exhibit a single crossing in half the Brillouin zone. Qualitatively, the spin-down electrons are transported upwards, while the spin-up electrons are transported downwards indicating that one would be able to assign spin Chern numbers of $C_{\uparrow}=-C_{\downarrow}=1$ if the Hamiltonian of the system could be continuously connected to a Hamiltonian that commutes with $S_{z}$ without closing the gap.

II.3 Mirror Crystalline Topological Insulators

As shown by Fu,Fu2011 crystal symmetries alone may give rise to a topologically non-trivial band structure. However, the consequences (gapless boundary states) are only observable at edges or surfaces that conserve the crystal symmetry. In this respect, mirror symmetry comprises a particularly simple type of crystal symmetry that gives rise to an integer topological classification. Whenever a material has mirror symmetry all the occupied states may be labeled according to their mirror eigenvalues $\pm i$ and one can define Chern numbers $C_{\pm}$ for these subsets of bands. The total Chern number is then $C=C_{+}+C_{-}$, but we might as well consider the $\mathbb{Z}\times\mathbb{Z}$ classification based on $C_{+}$ and $C_{-}$. In the case where the total Chern number vanishes the topology can be specified by the mirror Chern number $C_{M}=(C_{+}-C_{-})/2$, which is readily verified to be an integer.

For 2D materials, mirror symmetry in the plane of the material plays a special role, since any clean edge of the material will conserve the symmetry. A 2D material with non-vanishing mirror Chern number is thus guaranteed to host gapless edge states and gapless boundary states if the material is interfaced with another material that exhibits mirror symmetry in the plane. Strictly speaking, the interface or edge itself may break the symmetry by reconstruction of adsorption of atoms or molecules and as such the gapless edge states are more fragile than in the case of QSHI where the topology is defined by time-reversal symmetry. If $C=0$ an edge will host $C_{M}$ edge states of positive chirality and $C_{M}$ edge states of negative chirality. In Fig. 4 we show the Berry phase spectrum of SnTe, which is an example of a mirror crystalline topological insulator with $C_{M}=2$.Liu2015c We note that spinorial part of the mirror operator in the $z$-plane is represented as $i\sigma_{z}$ and the mirror eigenvalues thus closely follow the spin of the hybrid Wannier functions. This is evident in Fig. 4, where the spin up states are largely transported downward ($C_{+}=2$) and the spin down states ($C_{-}=-2$) are largely transported upward. This is of course only a qualitative argument, since the mirror operator also affects the orbital part, but the spin structure of the hybrid Wannier functions does allow one to obtain an intuitive picture of the (mirror-resolved) charge transport.

It is interesting to note that in the presence of time-reversal symmetry, the $\mathbb{Z}_{2}$ index can be obtained as $C_{M}$ mod 2. This is due to the fact that time-reversal symmetry enforces mirror symmetry on the Berry phase spectrum around the vertical axis at the $\Gamma$-point. In particular, we can represent the time-reversal symmetry operator as $T=\sigma_{y}K$, where $K$ denotes complex conjugation and the mirror symmetry operator by $M_{z}=i\sigma_{z}$. Due to the fact that $\sigma_{y}$ and $\sigma_{z}$ anti-commute, it follows that if $|u_{+}\rangle$ is an eigenstate of $M_{z}$ with eigenvalue $+i$ then $T|u_{+}\rangle$ will be an eigenstate of $M_{z}$ with eigenvalue $-i$. The Berry phase spectrum of the negative eigenvalue sector can thus be obtained from that of the positive mirror eigenvalues by reflection through the vertical line at $\Gamma$. This results in an odd number of crossings in half the Brillouin zone if $C_{M}$ is odd and an even number of crossings if $C_{M}$ is even. Materials with odd mirror Chern number and time-reversal symmetry thus have a dual topological character such that the edge states will be remain protected if either mirror symmetry or time-reversal symmetry is broken. This can be regarded as a 2D analogue of the dual topological character in Bi2Se3.Rauch2014 ; Olsen2016c In that case the strong topological index can be calculated from the Chern-Simons axion coupling $\theta$, which is restricted to the values of [math] and $\pi$ in the presence of either mirror symmetry or time-reversal symmetry.

III Results

III.1 Computational details

All the calculations presented in the this work were based on the Computational 2D Materials Database (C2DB), which currently contains DFT calculations for 3331 2D materials.Haastrup2018 The calculations were performed with the electronic structure software package GPAW,Enkovaara2010a which is based on the projector augmented wave methodblochl combined with the Atomic Simulation Environment (ASE).Larsen2017 All materials have been fully relaxed with the PBE functionalpbe and treated according to a strict work flow and a wide range of properties are calculated for materials that are both dynamically and thermodynamically stable. We refer to Ref. Haastrup2018, for details on the calculations. The database can be browsed online at https://c2db.fysik.dtu.dk or the full database can be downloaded from the repository.

As a first screening for 2D materials with non-trivial topology, we have sorted the C2DB for insulators with direct PBE band gaps below 0.7 eV. This primary criterion is based on the fact that the topological gaps are driven by spin-orbit coupling and we do not expect spin-orbit coupling to open band gaps by more than 0.7 eV. For all these materials, we carried out the parallel transport described in Sec. II and identified topological insulators as the materials with a gapless Berry phase spectrum. The procedure is highly convenient for automated screening, because no reference to the type of topological insulator is needed for the calculations and one can simply sort out whether a given materials is a QAHI, QSHI or a crystalline topological insulator (protected by mirror symmetry) afterwards.

The only non-standard ingredient in the procedure is the calculations of the matrix elements

[TABLE]

which is needed for the parallel transport algorithm. Within the PAW formalism the all-electron wavefunctions are written as

[TABLE]

where $|\tilde{\psi}_{n}\rangle$ are soft pseudo-wavefunctions, $|\phi_{i}^{a}\rangle$ are atomic orbital of atom $a$ and $|\tilde{p}_{i}^{a}\rangle$ are projector functions satisfying $\langle\tilde{p}_{i}^{a}|\tilde{\phi}_{j}^{a}\rangle=\delta_{ij}$. Denoting the position of atom $a$ by $\mathbf{r}_{a}$ the matrix elements can thus be written as

[TABLE]

where we assumed that $\langle\phi_{i}^{a}|e^{-i\delta\mathbf{k}\cdot\mathbf{\hat{r}}}|\phi_{j}^{a}\rangle=e^{-i\delta\mathbf{k}\cdot\mathbf{r}_{a}}\langle\phi_{i}^{a}|\phi_{j}^{a}\rangle$, since the partial waves $|\phi_{i}^{a}\rangle$ are localized at the atom $a$ and $e^{-i\delta\mathbf{k}\cdot\mathbf{r}}$ is a slowly varying function when $\delta\mathbf{k}$ is small. All quantities entering in Eq. (III.1) are calculated during standard DFT calculations with GPAW and are thus readily available for Berry phase calculations.

III.2 Overview of topological insulators in C2DB

In Tabs. 1 and 2 we provide an overview of all the topological insulators found in the screening. We emphasize again that all these were simply identified by looking for materials with a gapless Berry phase spectrum. The topological indices $\nu$, $C$, and $C_{M}$, relevant for QSHIs, QAHIs and MCTIs respectively, were then assigned to the different materials afterwards.

Tab. I contains known topological 2D materials. We find most of the materials that have previously been predicted to exhibit a non-trivial band topology. For example, grapheneKane2005a and its derivatives silicene, germanene, and staneneEzawa2015 as well as the transition metal dichalcogenides MoS2, MoSe2, and WSe2 in the 1T’ phase (WTe2 crystal structure).Qian2014 However, some well-known 2D topological insulators are missing from this table. For example, WTe2 and WS2 in the 1T’ phase. These materials are present in C2DB, but are semi-metals in the PBE approximationQian2014 and are therefore excluded from the present compilation. In addition CoBr2 in the CdI2 prototype has previously been predicted to be a QAHI based on PBEsol,Chen2017a which we have confirmed but the material is metallic within PBE and is therefore not included here. On the other hand, we find both HgSe and HgTe to be QSHIs although these materials have previously been reported to be trivial insulators based on calculations with a modified Becke-Johnson LDA functional.Li2015

Tabs. II-IV contains all the topological insulators that to our knowledge have not been reported prior to this work. The dynamically stable compounds (shown in bold face) are likely to be the only ones that are experimentally relevant, but we include all materials that we have found for completeness.

We find 27 new $\mathbb{Z}_{2}$ topological insulators that are dynamically stable (Tabs. II and III). We start by noting that the list is dominated by a 15 materials in the CuI crystal structure, with the antimony halides exhibiting band gaps exceeding 0.4 eV. Among the remaining materials, PdSe2 in the MoS2 crystal structure is the one with the largest band gap of 0.23 eV. However, it is situated 0.25 eV/atom above the convex hull with two other materials of the same stoichiometry being more stable. Nevertheless, the difference in stability is similar to the difference between MoS2 in the 2H phase and the 1T’ phase, which are both accessible by modern synthesis techniques. To test the reliability of the topological gap obtained with PBE, we have calculated LDA and $GW$ band structures fully including spin-orbit coupling in the self-energy. The result is shown in Fig. 5. The LDA gap is 0.225 eV, which is very similar to the PBE value. However, the $GW$ gap is 0.65 eV, which nearly comprises a threefold increase of the Kohn-Sham gap. In fact, to our knowledge this is one of the largest gaps reported for a two-dimensional topological insulator. A similar dramatic increase was recently reported for a Jacutingaite crystal structure where PBE yielded a gap of 0.15 eV and GW predicted a gap of 0.5 eV.Marrazzo2017a Such a dramatic increase in band gap has not been reported for three-dimensional topological insulators and could be related to the reduced screening in 2D. In principle, the predicted topology of the materials may depend on the approximation used to obtain the eigenstates. For example, in Ref. Vidal2011, it was that shown that DFT can lead to a False-positive conclusion for the non-trivial topology of certain 3D materials and GW calculatians may reverse the band inversion leading to the predicted non-trivial topology in DFT. However, we will not perform full GW calculations (with spin-orbit coupling) for all the materials in the present work, but note that false-positives could be a caveat for the present method. We discuss this issue further below in the context of QAHIs.

In addition to the 6 well-known mirror crystalline topological insulators displayed in Tab. 1, we also find 3 osmium dichalcogenides in the MnS2 crystal structure with a mirror Chern number of $C_{M}=2$, which are displayed in Tab. III. Only one of these - OsO2 - is stable though. As noted previously, any $\mathbb{Z}_{2}$ odd topological insulator with mirror symmetry must have an odd Chern number. In those cases the $\mathbb{Z}_{2}$ index does not exhaust the topological properties and we may distinguish the topological classes corresponding to $C_{M}=1$ and $C_{M}=3$ for example as in the case of RhO and AuTe in the GaS crystal structure. Again the physical consequences only emerge when considering an edge where the difference in Chern numbers would yield the number of protected gapless edge states. In general, one would not expect an interface between two $\mathbb{Z}_{2}$ topological insulators ($\nu=1$) to exhibit gapless interface modes. However, any interface between RhO and AuTe any interface that conserves the mirror symmetry would host four topologically protected gapless edge modes in - two for each mirror sector.

We would also like to emphasize that the screening has resulted in 6 new candidates for quantum anomalous Hall insulators, which are displayed in Tab. IV. This is of particular interest, since an experimental demonstration of the quantum anomalous Hall effect in a pristine 2D material is still lacking. Specifically, the materials FeCl3 and FeBr3 are highly stable and situated less than 0.1 eV/atom above the convex hull. However, the band structure of FeBr3 shown in Fig. 2 exhibits rather flat bands and indicates that the electrons in these materials are strongly correlated. It is thus likely that the PBE band gap provides a poor estimate of the fundamental gap of the material and even the topological properties could be wrong if PBE does not describe the band inversion correctly. In order to test the reliability of PBE we have tested the topological properties with various semi-local functionals and with PBE+U. Using LDA, RPBE, and PBEsol yield a trivial topology whereas PBE and revPBE predicts a QAHI with $C=1$. Using PBE+U we obtain a QAHI for values of U below 0.18 eV and a trivial insulator for U$>0.18$ eV. In all cases the geometry was optimized with the given functional. The topological properties of FeBr3 are thus highly sensitive to the method used and with the methodology applied here it is not possible to determine whether or not the material is a QAHI.

III.3 Magnetic anisotropy

The quantum anomalous Hall effect is driven by spin-orbit mediated band inversion in magnetic materials. Moreover, as a consequence of the Mermin-Wagner theorem magnetic order cannot exist without magnetic anisotropy in 2D materials and spin-orbit effects thus have another crucial role to play for these materials. In fact, the Curie temperature in a 2D ferromagnet is strongly dependent on the anisotropy as well as the exchange coupling constants and the first example of 2D ferromagnetic order was only observed very recently in the trivial insulator CrI3. Some of the present authors have shown that the Curie temperatures in 2D can be obtained from Monte Carlo simulationsTorelli2018 based on the classical Heisenberg model

[TABLE]

where $\langle ij\rangle$ denotes sum over nearest neighbors. The exchange, and anisotropy constants $J$, $A$, and $B$ can be obtained from first principles calculations combined with an energy mapping scheme that assumes magnetic moments to be located on transition metal atoms.Olsen2017 ; Torelli2018 Since the magnetic structure is typically isotropic in the plane of the material a minimal requirement for magnetic order at finite temperatures is an out-of-plane easy axis. For the magnetic materials in Tab. IV we find $T_{c}=2K$ and $T_{c}=274K$ respectively for FeBr3, CoBr3. The remaining magnetic materials are predicted to lack magnetic order at any finite temperature due to an in-plane easy axis or, more precisely, a negative spin-wave gap.Torelli2018 From these calculations CoBr6 appears highly promising. However, additional calculations shows that this material is more stable in the BiI3 crystal structure (similar to CrI3), which is non-magnetic and a trivial insulator. We note that the materials with an in-plane easy axis that are excluded here could in general give rise to a finite critical temperature if there is additional small in-plane anisotropy such that the rotational symmetry is explicitly broken. It is, however, rather difficult to predict the critical temperature in those cases and we expect critical temperatures to be low due to the approximate in-plane anisotropy. One exception to this may be provided by OsCl3Sheng2017 , which has a rather large magnetic anisotropy of $\sim 40$ meV per Os atom.

It is interesting to note that small gap magnetic insulators with large spin-orbit coupling may undergo topological phase transitions upon a rotation of the spin structure. As an example we take OsO2 in the CdI2 crystal structure, which has an in-plane easy axis and is not found in Tab. IV because it is a trivial insulator. However, if the magnetic moments are rotated out of plane to align with the axis perpendicular to the plane of the material it becomes a QAHI with a Chern number of $C=-2$. This is only possible if the band gap closes at some point when the magnetic moments are rotated from the in-plane to the out-of-plane configuration. In Fig. 6 we show the band gap as a function of the polar angle $\theta$ that the magnetic moments form with the $z$-axis. The band gap is seen to close in the vicinity of $\theta=\pi/4$ where the system undergoes a topological phase transition from $C=0$ to $C=2$. Such a rotation can be accomplished by applying a magnetic field and comprises a mechanism under which the gapless edge states can switched on or off by external means.

IV Discussion

We have implemented and performed an automated search for topologically non-trivial 2D materials in the Computational 2D Database. The method is based on a direct evaluation of Berry phases from the Kohn-Sham states and circumvents the common mapping to tight binding models via Wannier functionsGresch2017 that can sometimes make automation cumbersome.

In addition to several well-known topological insulators we have found 45 new materials of which, 18 are predicted to be stable. Of particular interest are the 6 magnetic QAHI. The experimental demonstration of the quantum anomalous Hall effect in a pristine 2D material would constitute a major breakthrough in the field of topological materials science. However, even if any of these materials could be synthesized the experimental verification of the effect will be highly tricky, since all of the experimentally relevant QAHIs have band gaps below 0.1 eV. Moreover, we have shown that the theoretical prediction of the topological properties is non-trivial due to strong correlation and the predictions based on PBE calculations may not be reliable as we exemplified in the case of FeBr3. Nevertheless, the PBE predictions for QAHIs presented here provide indications that some of these materials could be highly interesting to put under experimental scrutiny. Most notable CoBr3 in the AgBr3 crystal structure, which we predict to have a Curie temperature of 274 K, but which is also predicted to be less stable compared to the similar BiI3 crystal structure.

The 10 stable QSHIs found have PBE band gaps between 0.05 and 0.23 eV. We have only performed full spinorial GW calculations for PdSe2 where we found a threefold increase of the band gap. It would be highly interesting to perform carefully converged spinorial GW calculations for all the topological insulators presented in Tabs. I-IV, but this is beyond the scope of the present work. However, based on the cases of PdSe2 and jacutegaiteMarrazzo2017a we believe it is likely that several of the materials could have significantly larger gabs than predicted by PBE. None of the materials in Tabs. II-IV have been synthesized yet. But considering the rapid pace at which experimental techniques are currently evolving, we expect an experimental realization of one or several of the predicted topological insulators should be within reach in near future.

V Acknowledgments

This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 773122, LIMA). T. Olsen and D. Torelli was supported by the Independent Research Fund Denmark, grant number 6108-00464A. T. Okugawa acknowledges support from the Scandinavia-Japan Sasakawa Foundation. T. Deilmann acknowledges financial support from the Villum foundation.