Engineering Kitaev exchange in stacked iridate layers: impact of inter-layer species on in-plane magnetism

Electrostatic effects involving the inter-layer species are important: the largest Kitaev interactions comewith a more isotropic distribution of inter-layer cations around a given ligand.


Introduction
The prospect of realizing spin-liquid (SL) ground states in layered honeycomb materials with strong spin-orbit interactions 1,2 has triggered intense research activity in relation to these lattice systems. Quantum SLs are of particular interest in connection with properties such as protection of quantum information and the emergence of Majorana fermions. On a honeycomb lattice (Fig. 1), the essential ingredient for the formation of a quantum SL state is the so-called Kitaev coupling (K) between nearest-neighbor (NN) magnetic sites, a bonddependent Ising-like exchange 1,2 that must be large enough as compared to the more conventional NN Heisenberg J. It reaches quite robust values for d 5 electron congurations in iridium honeycomb oxides such as Na 2 IrO 3 (ref. 3 and 4) but also in the ruthenium halide RuCl 3 . 5,6 In the latter, a SL phase is realized by applying an external magnetic eld. 7,8 One peculiar prediction on the computational side is an enhancement of the Kitaev coupling K at large Ir-O-Ir bond angles. 9 The Ir-O-Ir bond angles are 90 for cubic edge-sharing octahedra, but in most honeycomb compounds they reach larger values due to trigonal compression of the oxygen cages (see Fig. 2). The largest Ir-O-Ir bond angles so far have been reported for H 3 LiIr 2 O 6 , nearly 100 . 10 Interestingly, Kitagawa et al. inferred a SL ground state for this material. 11 We examined in this context the Kitaev interactions of H 3 LiIr 2 O 6 but for ideal stacking of the honeycomb layers found rather modest K values as compared to, e.g., Na 2 IrO 3 (ref. 3 and 4) and earlier predictions for 100 Ir-O-Ir angles. 9 In an attempt to reconcile these apparently contradicting sets of computational results for large Ir-O-Ir bond angles, we addressed in detail the effect of having a single adjacent H site for each O ion and only 'vertical' O-H-O paths for the simplest stacking pattern. 10 We establish that the axial potential created through this kind of O-H coordination polarizes and bends towards the O-H link the O 2p orbitals bridging between Ir t 2g components orthogonal to the Ir 2 O 2 plaquette (see Fig. 3). Such polarization effects are absent in Li 2 IrO 3 and Na 2 IrO 3 , for coordination with several inter-layer cations of the O sites, but for ideal stacking in H 3 LiIr 2 O 6 (ref. 10) they disrupt Ir-O-Ir electron hopping and consequently reduce the Kitaev exchange.
Numerical tests in which the two H ions next to an Ir 2 O 2 plaquette are simply removed yield an impressively large ferromagnetic (FM) |K| value of 40 meV. Given the experimental indications for a SL ground state in H 3 LiIr 2 O 6 , 11 these computational ndings provide additional support for the existence of stacking faults 10 and H-ion disorder in this system, since larger K's should in principle make the quantum SL more likely. Moreover, our results provide valuable guidelines for the rational design of Kitaev quantum magnets, indicating that electrostatic effects involving inter-layer species are very important and that the largest K values come with a more isotropic distribution of inter-layer cations around a given ligand.

Quantum chemistry exchange couplings
Each IrO 6 octahedron shares edges with three NN IrO 6 octahedra such that the Ir sites frame a two-dimensional (2D) honeycomb lattice in H 3 LiIr 2 O 6 . Similar to its parent compound a-Li 2 IrO 3 , a Li ion is present at the center of each Ir 6 hexagon; an essential structural difference is the replacement of Li species between adjacent honeycomb layers by H ions. The octahedral ligand coordination of the Ir sites gives rise to a large gap between the t 2g and e g 5d levels. The leading ground-state conguration is therefore 2 T 2g (t 5 2g ), with all ve valence electrons in the t 2g orbitals. This corresponds to an L ¼ 1 orbital angular momentum 12 and via strong spin-orbit coupling yields a magnetically active j eff ¼ 1/2 ground-state doublet and lowerlying, occupied j eff ¼ 3/2 states. 2,12 The remarkable feature displayed by NN t 5 2g ions with strong spin-orbit interactions and edge-sharing connectivity of the encircling ligand cages is a large Ising-like interaction KS g iS g j between spin components perpendicular to a given M 2 L 2 plaquette (L stands for a ligand, see Fig. 2 and 3 for more details). On a honeycomb network of ML 6 octahedra this Ising-like coupling is bond dependent, i.e., the index g(g˛{x,y,z}) is different for each of the three M-M links emerging out of a given metal site M. Following Kitaev's conceptualization and initial analysis, 1 strong interactions of this type were suggested to be realized in Ir oxide compounds such as Na 2 IrO 3 and Li 2 IrO 3 . 2 To derive the strength of K in the related material H 3 LiIr 2 O 6 , and also of additional symmetric anisotropic terms and of the isotropic Heisenberg component, we rely on ab initio many-body computational schemes from quantum chemistry. The NN superexchange is analyzed on clusters of two edge-sharing IrO 6 octahedra, embedded in an effective eld that models the remaining part of the crystalline lattice.
Pairs of adjacent IrO 6 octahedra of two slightly different types were reported on the basis of X-ray diffraction data, 10 with Ir-O-Ir bond angles of either 99.0 or 99.8 degrees. Since the difference between these two values is rather small, we consider in our calculations a slightly idealized crystal structure with 'averaged' Ir-O-Ir bond angles of 99.4 and Ir-Ir bond lengths of 3.08Å. For two NN octahedra, the C 2h point-group symmetry allows two extra, symmetric off-diagonal exchange terms in addition to the isotropic Heisenberg and anisotropic Kitaev components. 4 The effective spin Hamiltonian for a pair of pseudospins at NN Ir sites i and j can then be written as where a, b˛{x, y, z}. We use in the following a local Kitaev reference frame 2,4 in which the z axis is perpendicular to the Ir 2 O 2 plaquette (i.e., g ¼ z); given the C 2h symmetry, G zx ¼ ÀG yz  in this setting. The NN magnetic couplings discussed in the following are derived by mapping 5,13 the ab intio quantum chemistry data onto such an effective spin Hamiltonian. Multicongurational wavefunctions were obtained to this end by complete-active-space self-consistent-eld (CASSCF) calculations, 14 using an active space consisting of six t 2g orbitals at the two NN Ir sites (see the Methods section for additional details). Post-CASSCF, we also performed multireference conguration-interaction (MRCI) computations 14 accounting for single and double excitations out of the Ir 5d t 2g and bridging-ligand 2p orbitals. The reference CASSCF wavefunctions were variationally optimized for the lowest nine singlets and nine triplets, which subsequently entered the spinorbit treatment to yield 36 spin-orbit states. The lowest four of these states dene the actual magnetic problem of two interacting pseudospin-1/2 sites, as pointed out by Jackeli and Khaliullin, 2 and were mapped onto the effective spin Hamiltonian (1), in order to derive the NN effective exchange couplings. The other 32 levels lie at signicantly higher energy, S0.5 eV, and are associated with j z 3/2 to j z 1/2 excitations, 2,12 as shown by joint resonant inelastic X-ray scattering (RIXS) and quantum chemistry investigations. 15,16 This large gap between the two sets of spin-orbit states, Ir j z 1/2 and j z 3/2, ensures that decoupling these two sectors in the mapping procedure is a safe approximation. 2 Such a strategy is also widely applied in the context of f electron superexchange. 17 All calculations were performed using the MOLPRO quantum chemistry package. 18 Using the same methodology, MRCI exchange couplings in good agreement with experimental data were reported earlier for square-lattice iridates, 13,19 the pyrochlore Sm 2 Ir 2 O 7 , 20 and the perovskite CaIrO 3 . 21,22

Results and discussion
The M-L-M angle is one of the key factors in tuning the magnitude of the Kitaev and Heisenberg components: as pointed out in ref. 9 and 23, larger angles lead to larger anisotropic interactions in honeycomb t 5 2g oxides. Since in comparison to the related iridates Na 2 IrO 3 and a-Li 2 IrO 3 , the Ir-O-Ir angles are on the larger side in H 3 LiIr 2 O 6 (99-100 ), 10,11 one would expect Kitaev couplings of larger magnitude in this system. For the other iridates, FM values in the range of 15-20 meV were found by MRCI; 4,9 however, for angles close to 100 and ideal stacking of successive honeycomb layers, we here compute a FM K value of only 10.9 meV (see Table 1). This suggests some subtle differences between interactions in H 3 LiIr 2 O 6 and in, e.g., Na 2 IrO 3 , the identication of which constitutes a main purpose of this paper. The actual splittings among the lowest four 'magnetic' levels along with the other effective coupling constants are also listed in Table 1, as obtained by both CASSCF and MRCI calculations.
We note that results similar to those provided in Table 1 are obtained by spin-orbit MRCI 24 when considering the presence of two, structurally different Ir-Ir links (referred to as B1 and B2) of the type proposed in ref. 10 and 11, with somewhat different bond angles and bond lengths between the B1 and B2 blocks of NN IrO 6 octahedra. Estimates for J, K and the off-diagonal couplings were also derived in very recent studies on the basis of density-functional computations. 25,26 K, for example, becomes signicantly stronger in the latter investigations, by a few meV in ref. 25 and by a factor of nearly 2 in ref. 26.

Interlayer electrostatics, impact on Kitaev exchange
From a structural point of view, two groups of compounds can be identied within the family of 5d 5 honeycomb iridates: Na 2 IrO 3 and a-Li 2 IrO 3 , where each cation in-between the honeycomb-like sheets has six adjacent oxygen sites, and We nd that removal of two H's next to the bridging ligands results in a nearly four-fold increase of the Kitaev exchange between in-plane NNS ¼ 1/2 sites: from z À11 meV (see Table  1), it now reaches À40 meV in the spin-orbit MRCI calculations. The proximity of the positive H ions is associated with two different effects on the in-plane spin-spin interactions: (i) the 'bare' effect of the H-ion Coulomb potential on on-site orbital energies and intersite hopping matrix elements and (ii) O 2p orbital polarization effects that can additionally affect the Ir-O-Ir orbital overlaps and therefore, once again, the intersite hoppings. To determine which is the dominant mechanism, we performed extra computations at the CASSCF level. NN CASSCF magnetic couplings obtained for a cluster where each of the H ions next to a bridging ligand is removed and the associated ionic charge is redistributed within the embedding are listed on the rst line in Table 2. This K value, À27.4 meV, corresponds to the CASSCF states used as a reference in the conguration interaction calculations leading to the MRCI result K ¼ À40 meV. In a second step, we modeled those two H ions as simple point charges (PCs) but did not allow relaxation of the cluster Table 1 Splittings among the lowest four spin-orbit states, mapped onto the eigenstates of the effective model defined by (1), and the resulting effective exchange couplings for two edge-sharing IrO 6 octahedra (all values in meV); a slightly idealized crystal structure with averaged bond lengths and bond angles was used (see the text) orbitals. In other words, multicongurational computations were performed without orbital reoptimization, which are also referred to as frozen-orbital, CASFO calculations. The exchange interactions are somewhat suppressed due to the presence of the nearby positive charge (see the second line in Table 2). However, allowing the orbitals to fully relax, i.e., to react to the axial potential generated by the adjacent unit PCs, results in a much more drastic reduction of the NN magnetic couplings (see the third row in Table 2). This step by step analysis makes it clear that orbital polarization in response to the electrostatic potential induced by the inter-layer H cations is the primary cause of the lower NN interaction constants listed in Table 1.

Magnetic splittings CASSCF + SOC MRCI + SOC
A large amount of stacking faults was evidenced in H 3 LiIr 2 O 6 , 10 most probably related to the rare situation in which H is bridging two adjacent O sheets. Having the hydroxyl bond in mind, it has been pointed out that an alternative way of writing the chemical formula of this compound is LiIr 2 O 3 (OH) 3 . 10 An idealized picture arising from this formula is then that of alternating [LiIr 2 O 6 ] 3À and [LiIr 2 (OH) 6 ] 3+ honeycomb-like layers (or slabs), the latter with all bridging O's replaced by hydroxyl groups, as in the related material Li 2 Pt(OH) 6 . 28 The weak bonding between layers and the inherent stacking disorder is even better highlighted in such a representation, the frail hydrogen bonds O-H/O being more apparent. In this context, 10,28 our results strongly suggest the existence of both 'ideally stacked' 11 (i.e., weak, see Table 1) and 'fault-present' 10 (i.e., strong, see Table 2) exchange couplings in this system, which then makes the modelling of the extended magnetic lattice more complicated.
While the role of inter-layer ionic species is analyzed here for H 3 LiIr 2 O 6 , ongoing work 29 yields similar results for the copper iridate Cu 2 IrO 3 , 27 displaying similar O-M 0 -O interlayer contacts. These data nicely complement earlier ndings concerning the sensitivity of various effective magnetic couplings to the position and charge of secondary/tertiary cations in oxide compounds. 13,30 Quantities addressed in the earlier investigations were on-site parameters such as g factors 13 and zero-eld splittings. 30 Here it is explicitly shown that also the intersite magnetic couplings can be adjusted by using electrostatic effects involving ionic species beyond the crystalline region (i.e., 'beyond' the bridging ligands) that is commonly assumed to be of relevance.

Angle dependence, the Kitaev limit
Looking at tuning structural parameters having to do with the honeycomb-like [LiIr 2 O 6 ] 3À slab, we further determined the evolution of both K and J with gradually modifying the Ir-O-Ir bond angles. We focused on atomic congurations where the two H cations next to the bridging ligands are removed since such arrangements are found to yield very large K's and are additionally likely to occur in the actual material. To maintain overall neutrality, the formal ionic charge associated with the two H's was again redistributed within the embedding; the Ir-Ir distance was xed at 3.08Å while a variable Ir-O-Ir bond angle was achieved through gradual trigonal distortion of the IrO 6 octahedra (see Fig. 2). The resulting K's and J's are shown in Fig. 4. The remarkable feature is that for bond angles close to 98 , J / 0. That is, a purely anisotropic effective magnetic model can be realized according to our calculations for z98 , with FM Kitaev coupling constants as large as 31 meV. This critical angle dening on the computational side the Kitaev 'limit' is actually not so far from the bond angles reported for H 3 LiIr 2 O 6 (ref. 11)-according to the quantum chemistry results, only small structural modications would be required for reaching the regime of vanishing J. Interestingly, it turns out that not only the M-L-M angle constitutes here a tuning knob but also the interatomic distances; especially in the vicinity of the critical value q c , J can be reduced towards zero by also varying the bond lengths, via tensile strain for instance. 31 Exchange between honeycomb planes A recent theoretical model 32 attempting to explain the experimentally observed magnetic properties of H 3 LiIr 2 O 6 (ref. 11) assumes interlayer isotropic couplings as large as 10 meV. To verify this assumption we carried out an additional set of calculations, on a cluster having as the active region two IrO 6 octahedra that belong to adjacent honeycomb layers and are connected through double O-H-O paths as proposed in ref. 11 (see Fig. 2). A rather similar kind of connectivity is in fact also    34,35 and it becomes apparent that this research area holds much potential for engineering magnetic couplings in Kitaev-Heisenberg systems.

Methods
The magnetic exchange couplings between NN Ir sites were derived by calculations on embedded clusters having two edgesharing octahedra (Ir 2 O 10 units) as the central region. To properly describe multiorbital physics within the Ir t 2g sector, we rely on a CASSCF scheme 14 as the starting point. In this frame, the most obvious choice for the active multicongurational space is that based on having six t 2g orbitals (three t 2g orbitals per Ir site) and ten electrons (two holes in the Ir t 2g channel). Since all possible electron congurations are here accounted for, we naturally describe in this way superexchange processes involving virtual excited states of t 4 2g -t 6 2g type. Additional intersite excitations, both M-M (t 2g -e g ) and L-M, enter our correlation treatment in the subsequent MRCI calculations. 36,37 Spin-orbit couplings were computed according to the methodology described in ref. 38. To model the nite charge distribution in the immediate neighborhood, the adjacent four octahedra were also explicitly included in the calculations. Energy-consistent relativistic pseudopotentials along with quadruple-zeta basis functions 39 were used for the Ir ions of the central unit. All-electron basis sets of quintuple-zeta quality 40 were employed for the bridging O ligands while all-electron basis sets of triple-zeta quality 40 were used for the other O anions within the two-octahedron central region. Ir 4+ sites belonging to octahedra adjacent to the reference unit were described as closed-shell Pt 4+ t 6 2g species, using relativistic pseudopotentials and valence triple-zeta basis functions. 39 Ligands of these adjacent octahedra that are not shared with the central reference unit were modeled with minimal all-electron atomic-natural-orbital basis sets. 41 For the Li NNs we employed total-ion effective potentials with a single s valence basis function 42 while two s and one p valence basis functions were used for the H NNs. 43 To extract the effective magnetic couplings, the lowest four spin-orbit states associated with the t 5 2g -t 5 2g manifold in the quantum chemical treatment were mapped onto the eigenstates of the effective spin Hamiltonian (1). To evaluate all symmetryallowed coupling constants in (1), we additionally considered the Zeeman coupling termĤ For the latter, all required matrix elements are available in the MOLPRO output data, i.e., the expectation values of the orbital angular momentum (L i , L j ) and spin (S i , S j ) operators (see also ref. 5 and 13 for additional details, in particular, the matrix elements listed in Tables 5 and 6 in ref. 5). The one to one correspondence between Hamiltonian matrix elements obtained at the ab initio and effective-model levels allows us to evaluate all coupling constants involved in (1). 5,13,20 Additional quantum chemistry computations were performed to investigate the strength of magnetic exchange between Ir sites belonging to adjacent honeycomb layers, connected via double O-H-O pathways. For this set of calculations, an Ir 2 O 12 H 2 unit was used as the central region. The six NN IrO 6 octahedra directly coordinating the two reference octahedra (three NN octahedra around each 'central' octahedron) plus eight H NNs lying close to the central fragment were also included in the calculations, for better representation of the immediate neighborhood. The two Ir ions in the central unit were represented by energy-consistent relativistic pseudopotentials along with triple-zeta basis functions. 39 All-electron basis sets of quintuple-zeta quality 40 were employed for the four O ligands on the O-H-O contacts while all-electron basis sets of triple-zeta quality 40 were used for the remaining O ions within the central region; triple-zeta basis sets 40 were applied for the two bridging H's. Ir 4+ sites belonging to octahedra adjacent to the central region were described as closed-shell Pt 4+ t 6 2g species, using relativistic pseudopotentials and valence double-zeta basis functions. 39 Ligands of these adjacent octahedra that are not shared with the central reference unit were modeled with minimal all-electron atomic-natural-orbital basis sets. 41 For the Li NNs we employed total-ion effective potentials with a single s valence basis function 42 while double-zeta basis functions were utilized for the remaining H ions. 40 PC embeddings were used in all calculations, as in earlier quantum chemistry investigations of honeycomb iridates. 4,9,15

Conflicts of interest
There are no conicts to declare.