The vicinity of hyper-honeycomb β-Li2IrO3 to a three-dimensional Kitaev spin liquid state

Due to the combination of a substantial spin-orbit coupling and correlation effects, iridium oxides hold a prominent place in the search for novel quantum states of matter, including, e.g., Kitaev spin liquids and topological Weyl states. We establish the promise of the very recently synthesized hyper-honeycomb iridate β-Li2IrO3 in this regard. A detailed theoretical analysis reveals the presence of large ferromagnetic first-neighbor Kitaev interactions, while a second-neighbor antiferromagnetic Heisenberg exchange drives the ground state from ferro to zigzag order via a three-dimensional Kitaev spin liquid and an incommensurate phase. Experiment puts the system in the latter regime but the Kitaev spin liquid is very close and reachable by a slight modification of the ratio between the second- and first-neighbor couplings, for instance via strain.

In magnetism, frustration refers to the existence of competing exchange interactions that cannot be simultaneously satisfied. Such effects can spawn new states of matter with quite exotic physical properties. Most famous in this regard are the different kinds of quantum spin liquids (QSL's) that emerge from frustrated spin couplings 1 . In these collective states of matter quantum fluctuations are so strong that they disorder the spins even at the lowest temperatures. The types of QSL states that then emerge range from chiral ones 2,3 to Z 2 topological spin liquids [4][5][6] carrying fractionalized excitations. Both experimentally and theoretically such QSL's have been observed and intensely studied in two-dimensional (2D) systems [1][2][3][4][5][6][7][8][9] . How this situation carries over to three spatial dimensions (3D), in which tendencies towards formation of long-range ordered magnetic states are in principle stronger and the disordering effect of quantum fluctuations therefore less potent, is largely unexplored. This is not only due to the limitations of theoretical and numerical approaches in 3D but also to the the sparsity of relevant candidate materials 10 . Very recently the latter however fundamentally changed through the synthesis of insulating Li 2 IrO 3 polymorphs [11][12][13] in which the magnetic moments of Ir 4+ ions form 3D honeycomb structures with threefold coordination. Here we concentrate on the β-Li 2 IrO 3 polymorph, which forms a so-called hyper-honeycomb lattice, see Fig. 1. Such a lattice might in principle support a 3D Kitaev spin liquid [14][15][16][17] , a direct counterpart of its lower-dimensional, 2D equivalent [18][19][20] .
The 2D Kitaev-Heisenberg model on the honeycomb lattice is characterised by the presence of large uniaxial symmetric magnetic couplings that cyclically permute on the bonds of a given hexagonal ring [18][19][20] . A QSL phase is present in this model if the ratio between the Kitaev interaction K and Heisenberg coupling J is larger than 8 20 . Quasi-2D honeycomb compounds initially put forward for the experimental realization of the Kitaev-Heisenberg Hamiltonian are 5d 5 and 4d 5 j ≈ 1/2 systems 19,21 such as Na 2 IrO 3 , α-Li 2 IrO 3 and Li 2 RhO 3 . Subsequent measurements evidenced, however, either antiferromagnetically ordered [22][23][24][25] or spin-glass 26 ground states in these materials.
The three factors that complicate a straightforward materialisation of the Kitaev QSL ground state in the quasi-2D honeycomb compounds are the presence of (i) appreciable additional exchange anisotropies [27][28][29] , (ii) two crystallographically inequivalent Ir-Ir bonds and (iii) longer-range magnetic interactions between second-and third-neighbor iridium moments 22,23,27,30,31 . These additional interactions push quasi-2D Na 2 IrO 3 and α-Li 2 IrO 3 towards the formation of long-range antiferromagnetic (AF) order at temperatures below 15 K. Also the 3D honeycomb system β-Li 2 IrO 3 orders magnetically: at 38 K the spins form an incommensurate (IC) Scientific RepoRts | 6:29585 | DOI: 10.1038/srep29585 ordering pattern 32 with strong ferromagnetic (FM) correlations 11 . Apparently additional interactions beyond only the nearest-neighbor (NN) Kitaev and Heisenberg ones are relevant also in the 3D system. This leaves two main challenges: first, one would like to precisely quantify the different magnetic exchange interactions between the Ir moments and second, one should like to determine how far away the magnetic ground state is from a Kitaev-type 3D QSL. Here we meet these challenges through a combination of ab initio quantum chemistry calculations by which we determine the NN magnetic couplings in β-Li 2 IrO 3 and exact diagonalization (ED) of the resulting effective spin Hamiltonian, on large clusters, to determine how far β-Li 2 IrO 3 is situated from the QSL ground state in the magnetic phase diagram.
The ab initio results show that the NN exchange in β-Li 2 IrO 3 is mostly FM, with relatively weak FM Heisenberg couplings of a few meV, large FM Kitaev interactions in the range of 10-15 meV, and additional anisotropies not included in the plain Kitaev-Heisenberg model. The sign and magnitude of second-neighbor Heisenberg couplings we determine from fits of the ED calculations to the experimental magnetization data. This second-neighbor effective coupling comes out as J 2 ≈ 0.2-0.3 meV and is thus small and AF. Remarkably, this AF J 2 stabilizes an IC magnetic structure that puts the system to be only a jot apart from the transition to a QSL ground state. Our findings provide strong theoretical motivation for further investigations on the material preparation side. The Kitaev QSL phase might be achieved by for instance epitaxial strain and relaxation in β-Li 2 IrO 3 thin films, slightly modifying the J 2 /K ratio.

Results
Quantum Chemistry Calculations. Quantum chemistry calculations were first performed for the on-site d-d excitations, on embedded clusters consisting of one central octahedron and the three adjacent octahedra (for technical details, see Supplementary Information (SI) and ref. 33). Reference complete-active-space (CAS) multiconfigurational wave functions 34 were in this case generated with an active orbital space defined by the five 5d functions at the central Ir site. While all possible occupations are allowed within the set of Ir 5d orbitals, double occupancy is imposed in the CAS calculations on the O 2p levels and other lower-energy orbitals. The self-consistent optimization was here carried out for an average of four states, i.e., T g  Table 1.
Due to slight distortion of the O cage 11 and possibly anisotropic fields associated with the extended surroundings, the degeneracy of the Ir t 2g levels is lifted. Without SOC, the Ir t g 2 5 states are spread over an energy window of   Table 1. Ir 4+ 5d 5 multiplet structure in β-Li 2 IrO 3 , all numbers in eV. Due to the noncubic environment, the T 2g /T 1g (and spin-orbit coupled j = 3/2) states are split appart. We still use however notations corresponding to O h symmetry. Only the lowest and highest Kramers doublets are shown for each set of higher-lying spin-orbit states.
Scientific RepoRts | 6:29585 | DOI: 10.1038/srep29585 ≈ 0.1 eV (see Table 1). Similar results were earlier reported for the quasi-2D honeycomb iridates 33 . The low-symmetry fields additionaly remove the degeneracy of the j = 3/2 spin-orbit quartet. With orbitals optimized for an average of 5d 5 states, i.e., T g  Table 1). If the reference active space in the prior CAS self-consistent-field (CASSCF) calculation 34 is restricted to only three (t 2g ) orbitals and five electrons, the relative energies of the j ≈ 3/2 components in the subsequent MRCI + SOC treatment are somewhat lower, 0.69 and 0.73 eV. The Ir t 2g to e g transitions require excitation energies of at least 3 eV according to the MRCI data in Table 1, similar to values computed for α-Li 2 IrO 3 33 .
While the quantum chemistry results for the on-site excitations in β-Li 2 IrO 3 resemble very much the data for the quasi-2D honeycomb iridates, the computed intersite effective interactions show significant differences. The latter were estimated by MRCI + SOC calculations for embedded fragments having two edge-sharing IrO 6 octahedra in the active region. As detailed in earlier work 27,35,36 , the ab initio quantum chemistry data for the lowest four spin-orbit states describing the magnetic spectrum of two NN octahedra is mapped in our scheme onto an effective spin Hamiltonian including both isotropic Heisenberg exchange and symmetric anisotropies. Yet the spin-orbit calculations, CASSCF or MRCI, incorporate all nine triplet and nine singlet states that arise from the two-Ir-site t g MRCI + SOC results for the NN effective couplings are listed in Table 2. The two, structurally different sets of Ir-Ir links are labeled B1 and B2, see Fig. 1. For each of those, the O ions are distributed around the Ir sites such that the Ir-O-Ir bond angles deviate significantly from 90°. While the B1 links display effective D 2 point-group symmetry (the effective symmetry of a block of two NN octahedra is dictated not only by the precise arrangement of the O ions coordinating the two magnetically active Ir sites but also by the symmetry of the extended surroundings), the B2 bonds possess C i symmetry, slightly away from C 2h due to small differences between the Ir-O bond lengths on the Ir 2 O 2 plaquette of two Ir ions and two bridging ligands (2.025 vs 2.023 Å 11 ). The absence of an inversion center allows a nonzero antisymmetric exchange on the B1 links. However, our analysis shows this antisymmetric Dzyaloshinskii-Moriya coupling is the smallest effective parameter in the problem -two orders of magnitude smaller than the dominant NN interactions, i.e., the Kitaev exchange. On this basis and further symmetry considerations (see the discussion in refs 27, 35-37), the effective spin Hamiltonian for the B1 links is assumed D 2h -like and in the local Kitaev reference frame (with the z axis perpendicular to the Ir 2 O 2 plaquette and x, y within the plane of the plaquette 19,27 ) it reads We find for the B2 links that slight distortions lowering the bond symmetry from C 2h to C i have minor effects on the computed wave functions and the quantum chemistry data can be safely mapped onto a C 2h model. For C 2h symmetry, the elements of the symmetric anisotropic tensor are such that Γ zx = − Γ yz .
The wave functions for the low-lying four states in the two-Ir-site problem can be conveniently expressed in terms of 1/2 pseudospins as in Table 2. In D 2 symmetry (B1 links) these pseudospin wave functions, singlet Φ S and triplet Φ 1 , Φ 2 , Φ 3 , transform according to the A u , B 2 , B 1 and A u irreducible representations, respectively. For (nearly) C 2h symmetry (B2 links), Φ S , Φ 1 , Φ 2 and Φ 3 transform according to A g , B u , B u and A u , respectively. The amount of Φ S -Φ 3 (B1) and Φ 1 -Φ 2 (B2) mixing (see Table 2) is determined by analysis of the "full" spin-orbit wave functions obtained in the quantum chemistry calculations.
As seen in Table 2, for each set of Ir-Ir links in β-Li 2 IrO 3 , B1 and B2, both J and K are FM. In contrast, J is AF for all pairs of Ir NN's in honeycomb Na 2 IrO 3 27 and features different signs for the two types of Ir-Ir links in α-Li 2 IrO 3 35 . The Kitaev exchange, on the other hand, is found to be large and FM in all 213 compounds, see Table 2 and refs 27 and 35. In addition to the Kitaev coupling, sizable off-diagonal symmetric anisotropic interactions are predicted. In β-Li 2 IrO 3 , these are FM for the B1 bonds and show up with both + and − signs for the B2 links (the sign of these terms is with respect to the local Kitaev reference frame), see Table 2.
Magnetic Phase Diagram. Having established the nature and the magnitude of the NN effective spin couplings, we now turn to the magnetic phase diagram of β-Li 2 IrO 3 . In addition to the NN MRCI + SOC data of Table 2, we have to take into account explicitly the second-neighbor Heisenberg interactions. Due to the 3D nature of the iridium lattice, with alternate rotation of two adjacent B2 bonds around the B1 link with which both share an Ir ion, one can safely assume that the third-neighbor exchange is vanishingly small. Results of ED calculations for an extended (pseudo)spin Hamiltonian including the MRCI NN interactions and a variable second-neighbor Heisenberg coupling parameter J 2 are shown in Fig. 2. Different types of clusters were considered, with either 16, 20 or 24 Ir sites. The 24-site cluster used in ED calculations with periodic boundary conditions is displayed in Fig. 2(a) while the structure of the smaller clusters is detailed in SI.
In order to investigate the magnetic properties of β-Li 2 IrO 3 , we calculated the static spin-structure factor ij i j i j along two paths denoted as θ (bc-diagonal) and φ (ab-diagonal) in Fig. 2(a), where the distance between neighboring B1 bonds is taken as 1. The results for several J 2 values with the 24-site cluster are plotted in Fig. 2(b). The propagation vector for each path θ φ q q ( / ) m m , determined as the wave number q providing a maximum of S(q), is plotted in Fig. 2(c). For J 2 = 0 the ground state is characterized by long-range FM order, i.e., , consistent with a previous classical Monte Carlo study 38,39 . Given the strong FM character of the NN exchange, ground states different from FM order are only obtained for finite AF J 2 . With increasing strength of the AF J 2 , q θ develops finite values starting at J 2 = J 2,c1 and reaches π at J 2 = J 2,c2 whereas q φ is finite but small in the range J 2,c1 < J 2 < J 2,c2 and zero otherwise. This evidences two magnetic phase transitions, from FM to IC order and further to a commensurate ground state. The latter commensurate structure corresponds to zigzag AF order, a schematic picture of which is shown in Fig. 2(d). The ED results for the four different types of periodic clusters are here in good overall agreement, as shown in Fig. 2(e). Some differences arise only with respect to the precise position of the critical points.
An intriguing feature is the appearance of a SL state in between the FM and IC phases. Since the total spin 2 S/N falls off rapidly and continuously near J 2 = J 2,c1 [see Fig. 2(c)], the FM ground state is expected to change into SL before reaching the IC regime. It can be confirmed by a structureless static spin-structure factor, like nearly flat q-dependence of S(q) at J 2 = 0.65 in Fig. 2(b). In Fig. 2(e) we also provide the critical values marking the transition between the FM and SL states. This was estimated as the point where any of the 〈 ⋅ 〉   S S i j expectation values turn negative, which implies a collapse of long-range FM order. Importantly, we find that the SL phase shows up in each of the four different types of periodic clusters. A more detailed analysis of the spin-spin correlations is provided in SI.

Discussion
Typically, a commensurate-to-IC transition critical point tends to be overestimated by using periodicity. For estimating more precisely the critical J 2 values we therefore additionally studied clusters with open boundary conditions along the c direction. Also, for a direct comparison between our ED results and the experimentally observed magnetic structure, we introduce an additional path δ (ac-diagonal), sketched in Fig. 3(a). The size of the cluster along a and b has insignificant effect on the computed critical J 2 values because φ q m (ab-diagonal) is either zero, around the critical points (periodic 24-site cluster), or very small, in the IC phase (periodic 16-and 20-site clusters), as seen in Fig. 2(c).
The value of the propagation vector along the δ-path δ q ( ) m is shown in Fig. 3(b) as function of J 2 for various cluster "lengths" in the c direction. The inset displays a finite-size scaling analysis for the critical values. In the infinite-length limit, we find J 2,c1 = 0.02 and J 2,c2 = 1.43 meV. The corresponding phase diagram is provided in Fig. 3(c). Similar critical points, i.e., J 2,c1 = 0.02 and J 2,c2 = 1.48 meV, are obtained for θ q m (see SI). As shown in Fig. 3(b), the dropdown of 2 S/N near J 2 = J 2,c1 is more clearly seen than in the case of periodic clusters because the formation of IC order is not hindereded for open clusters. Defining the FM-SL J 2,c1 critical value as the point where 〈 ⋅ 〉   S S i j turns negative for any (i, j) pair, the SL phase in the vicinity of J 2 ≈ J 2,c1 = 0.02 meV would have a width of about 0.01 J 2 . In other words, a very tiny FM J 2 coupling may drive the system from FM order to a SL state. With further increasing J 2 , the system goes through an IC phase to AF zigzag order at J 2 = 1.43 meV.
To finally determine the value of J 2 in β-Li 2 IrO 3 , we fitted the magnetization curve obtained by ED calculations at T = 0 K [see Fig. 3(d)] to the experimental data at T = 5 K 11 . Such an exercise yields J 2 = 0.2-0.3 meV, i.e., J 2 ≈ 0.1 J 2,c2 , so that the system is relatively far from the instability to zigzag order but very close to the transition to the SL ground state. Since with increasing J 2 the propagation vector δ q m of the IC phase increases smoothly from that of the SL = δ q ( 0) m to that of the zigzag state π = δ q ( ) m , long-wavelength IC order with a small propagation vector is expected for β-Li 2 IrO 3 . By performing a finite-size scaling analysis of δ q m at J 2 = J 2,c1 (N) + 0.28 meV, we obtain π = . ± .  sonably well our theoretical estimate. The stabilization of an IC state by J 2 couplings has been previously discussed for 1D zigzag chains like the path we label here as δ in ref. 40.
The value extracted for J 2 from our fit of the magnetization data is thus within our theoretical framework fully consistent with the experimentally observed IC magnetic order in β-Li 2 IrO 3 . Nevertheless we find that the system is remarkably close to a three-dimensional spin-liquid ground state, which can be reached by a minute change of ~0.25 meV, an energy scale that corresponds to about 3 K, in the second-neighbour exchange parameter J 2 . Changes of this order of magnitude can easily be induced by pressure or strain.