First Principles Prediction of Topological Phases in Thin Films of Pyrochlore Iridates

While the theoretical and experimental study of topological phases of matter has experienced rapid growth over the last few years, there remain a relatively small number of material classes that have been experimentally shown to host these phases. Most of these materials contain bismuth, and none so far are oxides. In this work we make materials-specific predictions for topological phases using density functional theory combined with Hartree-Fock theory that includes the full orbital structure of the relevant iridium d-orbitals and the strong but finite spin-orbit coupling strength. We find Y2Ir2O7 bilayer and trilayer films grown along the [111] direction can support topological metallic phases with a direct gap of up to 0.05 eV, which could potentially bring transition metal oxides to the fore as a new class of topological materials with potential applications in oxide electronics.


I. BULK DFT CALCULATION FOR Y2Ir2O7
The DFT calculation of bulk Y 2 Ir 2 O 7 was carried out with WIEN2k 1 and Quantum Espresso (QE) 2 . The results obtained from these two codes agree well. In QE, the pseudopotentials are generated by the included ATOMIC code. In the pseudopotential generation we used the fully relativistic PBESOL functional 3 . Because QE can only accept norm-conserving pseudopotentials for finite spin-orbit coupling (SOC) when Wannier projection is desired, we took all pseudopotentials to be norm-conserving. The valence shell of Y ([Kr]4d 1 5s 2 ) includes the 4s, 4p, 4d, 5s, 5p orbitals, and the valence shell of Ir ([Xe]4f 14 5d 7 6s 2 ) includes 5s, 5p, 5d, 6s, and 6p orbitals. Including the semi-core states, that is, 4s, 4p in Y, and 5s, 5p in Ir improves the transferability of the pseudopotentials. The valence states of O ([He]2s 2 2p 4 ) include the 2s and 2p states. The cutoff energy in QE calculation was selected to be 150 Rydberg (Ry) for plane waves, and 600 Ry for the charge densities. The structure information of bulk Y 2 Ir 2 O 7 is from Ref. [4]. The positions of the iridium ions in the bilayer and trilayer films grown along [111] are shown in Fig.S1. Because the tight-binding fitting method has been successfully applied to period four transition metal oxides, such as LaNiO 3 5,6 , we first tried to fit the DFT results for bulk Y 2 Ir 2 O 7 to a tight binding model that included both direct and indirect hopping between the local t 2g orbitals of the Ir 4+ ions. The hopping matrix elements were generated using Slater-Koster (SK) 7 integrals. As shown in Fig.S2 above, in order to obtain even a semi-quantitative fit for the j = 1/2 bands (upper 4) one must include at least out to third-neighbor hoppings. Even in this case, however, the lower j = 3/2 bands are fit rather poorly. The failure of the tight-binding fitting reveals the complexity of 5d orbits, and their much larger spatial extent compared to the 3d orbitals. In order to obtain an satisfactory model, one must include more SK hopping terms in the tight-binding fitting, which is not well motivated for a practical model Hamiltonian. Instead, we turn to the Wannier states, making use of the Wannier90 package 8 , which yields the fit shown in Fig.S3.

III. WANNIER PROJECTION OF THE BULK DFT CALCULATION FOR Y2Ir2O7
The spin-orbit coupling can be treated in two different ways based on Wannier projection. One route is to first perform a DFT calculation and Wannier projection without SOC, and then add an onsite SOC (∼ λ(−l) · s) term "by hand" whose strength is determined through least-square fitting of the Wannier+SOC bands to the DFT bands with the SOC turned on 8,9 . Alternatively, one can perform a DFT calculation with SOC turned on and then obtain a direct Wannier projection of the SOC coupled states 8 .
In those Wannier projections, the initial basis is selected as the local t 2g orbitals. Even though the trigonal distortion causes a mixing between those t 2g orbitals, the axis of the local t 2g orbitals on each of the 4 sites in the unit cell can be assigned as 8 (S1) The coordinates in Eqs.(S1) can be easily obtained from the rotational matrices R (0)−(3) in Ref. [10]. All the Wannier functions are centered on the iridium ions. The spin quantization axis for all 4 iridium ions in the unit cell is set to the the global z-axis. To maximumly preserve the bulk crystal symmetry, the Num Iter is set to 0 in the Wannierization process. However, even if a finite Num Iter is used, the Wannier functions change only slightly, so the trial wave functions are quite good. When the SOC is fit "by hand", we find a spin-orbit coupling strength of about 0.43 eV, corresponding roughly to the splitting of the j = 1/2 and j = 3/2 manifolds. As shown in Fig.S3, this leads to a reasonable fit to the fully relativistic (SOC on) DFT calculation. This approach should be able to capture the main character of the phases when interactions are further treated at the Hartree-Fock level. By comparison, a direct Wannier projection for the DFT calculation with SOC turned on captures the band feature more accurately (the blue dashed line in Fig.S3) and will be more accurate in determining the phase diagrams.
From the analysis of the Wannier projection ones sees appreciable hopping between iridium ions that are as far as 3 to 4 FCC basis vectors from each other. Our results agree with a previous study in Sr 2 IrO 4 , in which the hopping between Ir 4+ ions 1 nm (around 3 times the distance between our nearest neighbor hopping) away from each other still plays a role 11 . can also be used. The original FCC basis vectors for the bulk pyrochlore lattice are given by, a 1 = (0, a/2, a/2), a 2 = (a/2, 0, a/2), where a is the lattice constant.
To obtain the thin film structure along [111] direction, for the n = 2, m = 2 bilayer thin film, the basis vectors are selected as Here a and b are in plane, while c contains a component out of the plane spanned by a and b. The position of each ion can be expressed in terms of linear combinations of a, b, c. The size of the unit cell has been doubled compared to the bulk, with Ir 4+ ions in two adjacent atomic layers substituted with Hf 4+ ions to obtain the unrelaxed bilayer structure shown in Fig.S4. As an initial guess, the lattice constant a can be taken as the average of Y 2 Hf 2 O 7 and Y 2 Ir 2 O 7 . The unrelaxed trilayer film structure can be obtained in a similar way. Some package such as ASE Surface 13 can also help to achieve this.
The determination of the lattice structure of the sandwhich is carried out in a fully relaxed scheme with scalarrelativistic 14 pseudopotentials for both the bilayer superlattice (Y 2 Ir 2 O 7 ) 2 /(Y 2 Hf 2 O 7 ) 2 and the trilayer superlattice (Y 2 Ir 2 O 7 ) 3 /(Y 2 Hf 2 O 7 ) 3 . We have verified that increasing the number of layers of Y 2 Hf 2 O 7 does not have a significant effect on the final lattice structures and the resulting band structures. Once the lattice structure is determined, the band structures of the bilayer and trilayer Y 2 Ir 2 O 7 are calculated in a fully relativistic basis (i.e., with the spin-orbit coupling turned on). The pseudopotential of Hf (with atomic configuration [Xe]4f 14 5d 2 6s 2 ) is generated in the same way as the other elements, with the valence orbits as 5s, 5p, 5d, 6s, and 6p. The k-space is meshed by 7×7×1 with Monkhost-Pack Grid 15 . The Wannier fits shown in Fig. 2 of the main text are based on the Wannier fits obtained for the bulk system, which are still fairly accurate for states near the Fermi energy.

V. THIN FILM HARTREE-FOCK CALCULATION
The unrestricted Hartree-Fock calculation is carried out in the full t 2g subspace of the iridium orbitals with Wannier states obtained from the bulk band structure. The resulting fits are shown in Fig. 2 of the main text. While the quality of the fits are less good than they are for the bulk, Fig.S3, they are sufficient for producing reliable trends in terms of the phases expected as a function of the on-site interaction U . Our thin film calculations start with 21 randomly generated coupling constants m i αβ = ĉ i † αĉ i β per site, which are then iterated to convergence with a difference between m i αβ input and output of less than 10 −8 . The lattice size is selected by be 90 by 90 or 120 by 120. Around fifty sets of initial guesses are taken, and about two-thirds of them converge for most U values. The final convergent configuration with the minimum total energy is selected.
The magnetic moments are calculated in a way similar to the determination of the Landé-g factor. The g factor is obtained by where j = −l + s is the total angular momentum in local coordinates. Here −l is the effective orbital angular momentum in t 2g subspace obtained by P t2g LP t2g = −l l=1 , where P t2g is the projection operation into the t 2g subspace, and g l = 1, g s = 2. Then the total magnetic moment is determined by Most of the magnetic moment originates in the j = 1/2 subspace. Similar to the situation reported in bulk Y 2 Ir 2 O 7 16 , we find an energy difference (around 1-150 meV/site, depending on the distance from the magnetic phase transition) between states with different order. All of our time-reversal symmetry broken solutions have non-zero net magnetization.
To compute the Z 2 invariant in the non-magnetic insulators, we use the formulation of Ref. [17]. We compute the Chern number in the magnetic insulating phases in a similar way 18 .