Prospect of quantum anomalous Hall and quantum spin Hall effect in doped kagome lattice Mott insulators

Electronic states with non-trivial topology host a number of novel phenomena with potential for revolutionizing information technology. The quantum anomalous Hall effect provides spin-polarized dissipation-free transport of electrons, while the quantum spin Hall effect in combination with superconductivity has been proposed as the basis for realizing decoherence-free quantum computing. We introduce a new strategy for realizing these effects, namely by hole and electron doping kagome lattice Mott insulators through, for instance, chemical substitution. As an example, we apply this new approach to the natural mineral herbertsmithite. We prove the feasibility of the proposed modifications by performing ab-initio density functional theory calculations and demonstrate the occurrence of the predicted effects using realistic models. Our results herald a new family of quantum anomalous Hall and quantum spin Hall insulators at affordable energy/temperature scales based on kagome lattices of transition metal ions.


I. FORMATION AND DOPING ENERGIES
Total energies were evaluated using ab-initio density functional theory (DFT) calculations within an allelectron full-potential local orbital (FPLO) [1] basis. For the exchange-correlation functional we employ the generalized gradient approximation (GGA) [2]. Total energies were extracted from calculations converged using 8 3 kpoint grids.
Experimental and hypothetical crystal structures were fully relaxed using the projector augmented wave (PAW) method [3] implemented in GPAW [4] with a plane-wave cutoff of 1000 eV and the GGA exchange-correlation functional [2]. We optimized the stoichiometric structures of herbertsmithite-based compounds using 6 3 kpoints (4 3 k-points for non-stoichiometric structures) until forces were below 10 meV/Å.
Hypothetical crystal structures were prepared starting from the experimental crystal structure of herbertsmithite [5], substituting zinc (Zn) atoms between the copper kagome layers by monovalent A=Li, Na (holedoping) and trivalent Al, Ga, In, Sc, Y (electron-doping). We refer to these compounds as A-herbertsmithite, ACu 3 (OH) 6 Cl 2 (case 1) and estimate their stability by comparison to clinoatacamite [Cu 2 (OH) 3 Cl]. We directly compare energies of A-substituted materials, where the dopant sits in the interlayer site versus structures where it occupied a kagome site (case 2, see Fig. 3 in main paper). To analyze the phase stability, we also evaluate total energies of 3 In order to determine whether the presence of excess copper in the synthesis process makes the formation of clinoatacamite favorable compared to ACu 3 (OH) 6 Cl 2 , we compare the total energies of Cu in a copper crystal plus the energy of ACu 3 (OH) 6 Cl 2 to the total energy of A in a crystal of substance A plus the total energy of clinoatacamite. Proper normalization of total energies to formula units is taken into account. The calculated energy difference determines whether the formation of metallic patches of substance A is energetically favourable compared to forming A-substituted herbertsmithite. The formation energy for this case is defined as where E A refers to the pure crystalline metal A, E AH to A-herbertsmithite and E C to clinoatacamite. As all total energies are negative, a negative formation energy signals stability of the target compound ACu 3 (OH) 6 Cl 2 . As a consistency check, we also prepared a structure with copper on the interlayer site so that the chemical formula is identical to clinoatacamite. For this structure, we find a positive formation energy. Therefore, hypothetical Cuherbertsmithite is unstable with respect to the formation of clinoatacamite, as expected.
For the formation energy of impurity or vacancy structures, we use where E AHM is the total energy of the compound Formation energies for modifications of herbertsmithite are shown in Table I. As all calculated energy differences are negative, formation of A-herbertsmithite is always favorable compared to formation of clinoatacamite. All proposed modifications are robust against formation of vacancies or copper impurities on the (interlayer) A-site.
The doping energies shown in Fig. 3 of the main paper are computed from energy differences between structures with dopant atoms occupying the interlayer sites and structures with dopant atoms on a kagome site (case 2). The values for the doping energies are given in Table I, together with ionic radii taken from Ref. [6].
Having established the stability of stoichiometric compounds, we investigate the solid solution of Zn and Ga TABLE I. Calculated formation and doping energies for substituted herbertsmithite. Energies are given in eV (electron volts) per formula unit of A-substituted herbertsmithite ACu3(OH)6Cl2 (A = Li, Na, Zn, Cd, Sc, Y, Al, Ga, In). Negative values signal stability of the stoichiometric compound against the formation of clinoatacamite (case 1), doping into the kagome layer (case 2, see Fig. 3 in main paper), vacancies on the interlayer site (case 3) and impurities on the interlayer site (case 4). Ionic radii in coordination number 6 are taken from Ref. [6]. In this configuration Cu 2+ has an ionic radius of 72 pm.  dopants. To accomodate non-integer dopant concentrations, we use 2 × 1 × 1 and 3 × 1 × 1 super cells. The formation energy of the target compound having total energy E T is defined as where E ZnH refers to herbertsmithite and E GaH to Gallium-substituted herbertsmithite. Note that as a reference energy here we use (1 − x) · E ZnH + x · E GaH , which is the energy that one obtains if the solid solution were to dissociate spontaneously into herbertsmithite and Gasubstituted herbertsmithite. Therefore, formation energies are negative if the solid solution is stable against phase separation. Our results show that the solid solution Ga x Zn 1−x Cu 3 (OH) 6 Cl 2 should exist in a broad range of doping ratios (Fig. 1).

II. DETAILS ON PREPARATION OF HYPOTHETICAL CRYSTAL STRUCTURES
For all cases discussed in the previous section, crystal structures had to be prepared from a stoichiometric structure of herbertsmithite [5]. We started from the rhombohedral unit cell of space group R3m, containing one formula unit. After reducing symmetry to space group P 1, herbertsmithite allows for only one possibility to arrange dopant atoms or defects for all compositions investigated here, if the smallest possible supercell is considered.  Fully relativistic bandstructure of Gaherbertsmithite. Plus and minus signs denote the parities of the three bands closest to the Fermi level at eight inversion-symmetric points (Γ, 3 × F , 3 × L, Z) in the Brillouin zone. In Ga-herbertsmithite the bandstructure is more three-dimensional than in Li-herbertsmithite, which leads to a substantial displacement of the Dirac point from the high-symmetry point K. In turn the apparent gap at K is exaggerated.

III. BANDSTRUCTURE OF HERBERTSMITHITE IN FERROMAGNETIC SPIN-CONFIGURATION
Herbertsmithite is an antiferromagnetic Mott insulator with possible quantum spin liquid ground state. Although density functional theory cannot describe such a paramagnetic ground state with large fluctuating moments, a DFT+U calculation for the most simple ordered magnetic configuration (ferromagnetic) already reproduces the insulating ground state of herbertsmithite with a large band gap of about 2 eV (see Fig. 3).

IV. BANDSTRUCTURE AND TOPOLOGICAL PROPERTIES OF ELECTRON-DOPED HERBERTSMITHITE
We also analysed the bandstructure and topological properties of electron-doped herbertsmithite. As an example, in Fig. 4 we show the fully relativistic non-spin-polarized electronic bandstructure of Gaherbertsmithite, where the parities of the three (dominantly Cu d x 2 −y 2 ) bands closest to the Fermi level are indicated. Ga-herbertsmithite is more three-dimensional than Li-herbertsmithite. Therefore, the Dirac point is displace from K, which exaggerates the gap at K. At the Dirac point the gap magnitude is the same as in Liherbertsmithite.
The paritiy eigenvalues of the bands closest to the Fermi level are the same as in all other compounds derived from herbertsmithite that we predict to be stable. That means, topological numbers of the bands below the Dirac point are ν 0 ; (ν 1 , ν 2 , ν 3 ) = 0; (111).

V. FORMALISM FOR COMPUTING SURFACE STATES
We use the Green's function method for the semiinfinite system [7][8][9] to calculate the states on the (001) surface of interlayer substituted herbertsmithite. In this scheme, the crystal must be partitioned into so-called principal layers parallel to the surface of interest. The Hamiltonian within the principal layer is denoted as H 0 , while the coupling between neighboring principal layers is denoted as C. Note that only hoppings with components in direction either negative or positive perpendicular to the surface must be included into C. The corresponding couplings in the reverse direction are taken into account in the formalism by the adjoint matrix C † .
The surface Green's function of a system consisting of N principal layers can be written as The initial condition is Indices i, j denote combined site, orbital and spin variables. The Green's function of the dual surface can be calculated by iterating As we use Green's functions on the real-frequency axis, an artificial imaginary part must be added to ensure numerical convergence. To this end, we transform ω → ω + iη with η = 10 −5 eV in the equations above. Although the method in principle allows for the treatment of semi-infinite systems, the numerical iteration of Eq. (4) or (6) has to be stopped with a finite number of layers. We use N = 10 5 .

VI. (001) SURFACE OF HERBERTSMITHITE
In the main paper we presented topological invariants ν 0 ; (ν 1 , ν 2 , ν 3 ) = 0; (111) for the herbertsmithite system. According to Ref. [10], a material with those topological invariants displays conducting states on a (001) surface. The orientation of this lattice plane refers to the unit cell in rhombohedral setting of space group R3m. The (001) surface of Li-herbertsmithite is shown in Fig. 5. Note how the Cu ions at the surface are arranged in chains, with all kagome layers terminated equivalently.
The unit cell of (substituted) herbertsmithite contains in total three Cu ions. These three ions comprise the minimal model of the herbertsmithite material. In the (001) surface geometry we use, two of the ions are at the surface and one is below the surface. The minimal principal layer therefore comprises Cu ions at different heights within the unit cell. This complicates the partitioning of hopping elements into matrices H 0 and C needed for the Green's function method. The situation is depicted in Fig. 6.
As mentioned in the main paper, and obvious from Fig. 6, the semi-infinite system can be built up either in positive or negative direction perpendicular to the surface. In the kagome lattice this results in an edge termination with either chains or triangles. The case of chains is presented in the main paper. For the dual surface (triangles termination) one obtains a spectral function with equivalent features (see Fig. 7) from Eq. (6).
Both in the main paper and in the supplement we show the spectral function A(k, ω) only in positive k-direction within the k x -k y -plane. As can bee seen from Fig. 5, the surface unit cell of (001) herbertsmithite contains two chains of Cu ions. The surface is mirror symmetric with respect to the k x and the k y -direction. Therefore, the symmetry equivalent edge states of the second kagome layer appear in negative k-direction.

VII. DETERMINATION OF HEISENBERG HAMILTONIAN PARAMETERS
Based on the theoretically predicted structures of Liherbertsmithite (see Table II, Fig. 8a), we use density functional theory calculations to determine the most important couplings of the Heisenberg Hamiltonian. We use the all electron full potential local orbital (FPLO) code [1] with a generalized gradient approximation [2] exchange and correlation functional and correct for the strong correlations on the Cu 2+ 3d orbitals using GGA+U [11].
The results for LiCu 3 (OH) 6 Cl 2 are given in Table III    allows for 171 out of 4096 unique spin configurations, and the data in Table III are based on fits to total energies for about 30 of these configurations. The exchange couplings are given assuming S = 1/2. However, due to the changed filling of the Cu bands, the moments in the calculation are reduced from 1 µ B to, on average, about 0.67 µ B . The relevant exchange paths are visualized in Fig. 8b. Based on previous work [12], we consider the GGA+U functional with U = 6 eV, J H = 1 eV most relevant for Cu in the square planar environment as realized in herbertsmithite type crystals.

VIII. ESTIMATE FOR THE CURIE TEMPERATURE
We estimate the Curie temperature T C from a simple Weiss mean-field formula [13]. For 3d transition metals, in this approximation the Curie temperature is given by Eq. (7), where z i is the coordination number and J i are the exchange couplings in Kelvin (U = 6 eV), taken from Table III.
We obtain T C = 1160 K, where we took into account i = {1, 3, 5} with z i = {4, 4, 6}. Instead of the spin-1/2 model used here, it is also conceivable to parametrize the Heisenberg model with DFT effective moments for S, which would result in a higher value for T C .