Strongly frustrated triangular spin lattice emerging from triplet dimer formation in honeycomb Li2IrO3

Iridium oxides with a honeycomb lattice have been identified as platforms for the much anticipated Kitaev topological spin liquid: the spin-orbit entangled states of Ir4+ in principle generate precisely the required type of anisotropic exchange. However, other magnetic couplings can drive the system away from the spin-liquid phase. With this in mind, here we disentangle the different magnetic interactions in Li2IrO3, a honeycomb iridate with two crystallographically inequivalent sets of adjacent Ir sites. Our ab initio many-body calculations show that, while both Heisenberg and Kitaev nearest-neighbour couplings are present, on one set of Ir–Ir bonds the former dominates, resulting in the formation of spin-triplet dimers. The triplet dimers frame a strongly frustrated triangular lattice and by exact cluster diagonalization we show that they remain protected in a wide region of the phase diagram.

Supplementary Tables: Supplementary Table 1: Exchange coupling parameters as obtained at the MRCI level of theory. Units of meV are used.   In the ab initio quantum chemistry calculations we used energy-consistent relativistic pseudopotentials with quadruple-zeta basis sets for the valence shells of the two reference Ir ions [1], all-electron quintuple-zeta basis sets for the bridging ligands [2] and triple-zeta basis functions for the other O's of the two reference octahedra [2]. We also applied polarization functions for the central Ir ions and the two bridging ligands, two Ir f and four O d functions [1,2]. For the NN Ir 4+ and Li + ions we employed total-ion effective potentials supplemented with [2s2p2d] and [1s] valence basis functions, respectively (Stoll 2013) [3]. Identical basis sets were used for Ir and O in the calculations for the idealized structures of Na213 (Fig. 2 in main text), with the same type of Na pseudopotentials and basis functions as for the Li ions in Li213.
The effective magnetic Hamiltonian for two adjacent Ir sites is most conveniently written in the local reference frame {X b ,Y b ,Z b } with the X b axis along the Ir-Ir bond and Z b perpendicular to the Ir 2 O 2 plaquette (see Methods in the main text). For C 2h point-group symmetry, it reads For simplicity, we omitted the bond index b. A straightforward diagonalization of H ij yields the following eigenvalues and eigenfunctions: Here Φ S is the total spin singlet and Φ 1−3 are combinations of the three triplet components. The latter are degenerate in the plain Heisenberg model. The above diagonalization procedure is equivalent with a rotation of the coordinate system {X,Y,Z} around X by an angle α to a new frame {X ,Y ,Z } in which the symmetric anisotropic exchange matrix is diagonal [15]. {X ,Y ,Z } are also referred to as principal axes and the angle α is given by In C 2h symmetry, the Φ S , Φ 1 , Φ 2 and Φ 3 spin-orbit wave functions transform according to the A g , B u , B u and A u irreducible representations, respectively. Since states Φ 1 and Φ 2 belong to the same irreducible representation B u , they are in general admixed, i.e., in the reference frame {X,Y,Z} the corresponding eigenfunctions should be written as The mixing parameter ζ = sin α is given by Since the quantum chemistry calculations were actually performed in C 1 symmetry, to determine the nature of each of the lowest four spin-orbit states we explicitly computed the dipole and quadrupole transition matrix elements within that manifold. Standard selection rules and the nonzero dipole and quadrupole matrix elements in the quantum chemistry outputs then clearly indicate which state is which. We also carried out the transformation of the spin-orbit wave functions from the usual {L 1 ,M L1 ,L 2 ,M L2 ,S,M S } basis in standard quantum chemistry programs to the {S 1 ,S 2 ,M S1 ,M S2 } basis. This allows the study of Φ 1 -Φ 2 mixing due to the offdiagonal Γ yz and Γ zx couplings.

Supplementary Note 2: Lattice spin model, conventions and ab initio effective couplings
For each type of Ir-Ir link, we used two different local axes frames, namely Let us choose our global frame to be {x, y, z} = {x 1 , y 1 , z 1 }. The relation between the local bond axes {X b , Y b , Z b } and the global coordinate system is the following (see Supplementary Figure 1): We note that the system is invariant under two-fold rotations around the vector X 1 . This is why we must choose X 3 to be the rotated version of X 2 , i.e., C 2 · X 2 = X 3 . If we adopt the opposite direction for X 3 then we must change the sign of the coupling C.
Let us now write down the NN terms of the Hamiltonian for the three types of bonds b = 1, 2, 3. In the reference frame where . For simplicity, we replace J b=1 by J, J b=2,3 by J and similarly for the remaining parameters. Supplementary Equations (7) and (6) then yield . With the conventions introduced in Supplementary Figure 1, the spin Hamiltonian now reads: where R = na + mb labels the unit cells and j = 1−2 labels the sublattice index.
The numerical values of the above coupling parameters are provided in Supplementary Tables 1 and 2.
Supplementary Note 3: Effective spin-1 description in momentum space As mentioned in the main text, the effective spin-1 model takes the following form in momentum space: where T R = 1 √ N k e −ik·R T k , N is the number of B1 bonds on the lattice, the matrix Λ is given by The classical phase diagram of H eff can be investigated by a numerical minimization of the eigenvalues of Λ(k) over the entire Brillouin zone (BZ) of the triangular lattice (shown in Supplementary Figure 3), which is essentially the Luttinger-Tisza method [16] for a Hamiltonian with anisotropic terms in spin space. In contrast to isotropic Hamiltonians, here not all directions in spin space are equivalent. This is an important aspect reflected in the form of the eigenvectors of Λ.
When the minimum sits at a commensurate wavevector Q -as in the ferromagnetic (FM), diagonal zigzag and stripy phases -the ground state is given by the corresponding commensurate modulation (see section 1 below) and the actual directions in spin space are set by the eigenvector corresponding to the minimum eigenvalue. If Q is an incommensurate point and the minimum eigenvalue at Q is at least two-fold degenerate then the ground state is given by a spiral configuration in the plane formed by the corresponding eigenvectors. The situation is more complex when Q is an incommensurate point and the minimum eigenvalue is non-degenerate, as in the regions ICx and ICy. In this case we cannot write down a spiral solution that satisfies the spin length constraint |T R | = 1 for all R and the actual nature of the ground state goes beyond the standard helical configuration (see also main text).
The classical phase diagram for the case of C = D = C = D = 0 is provided in Supplementary Figure 4(a), while the one for finite C, D, C and D is shown in Figure 3 of the main text. In the same figure we show the quantum phase diagram of the original spin-1/2 model, as obtained from exact diagonalization (ED) calculations on the 24-site cluster. As for the case of finite C, D, C and D , we do not find any traces of the incommensurate phase ICy in the ED calculations. This is either a finite-size effect or due to the fact that the classical ICy region is very narrow, meaning that the competing nearby FM and zigzag phases are very close in energy and are more likely to be favored by quantum mechanical fluctuations. The surprising result is that the quantum ICx region is extremely narrow for small J 2 and J 3 , drastically different from what we found for finite C, D, C and D . Since the classical minimization does not show such a drastic dependence of the ICx region on C, D, C and D , we can attribute this to finite-size effects in the ED calculations.
In the following two subsections we shall focus on the zigzag phase to demonstrate what happens for the simplest, commensurate cases.
Zigzag phase -limit of C = C = D = D = 0 In this limit, Λ(k) becomes diagonal: Now, the minimum of f (k) is at the two M points M 1 = (π, π 2 is minimum at both M 1 and M 3 , with the eigenvalue λ z = λ 0 − |K| 2 . Since |K|/2 > |K |/4, the lowest eigenvalue is λ z . This means that spins will prefer to point along the z axis and not within the xy-plane. Still, there are two M points available and we therefore have two possible solutions: T(na + mb) = e iM1·r e z = (−1) n+m e z , T(na + mb) = e iM3·r e z = (−1) n e z , where R = na + mb and n, m are integers.

Remark 1:
The direction in spin space basically boils down to the competition between the two anisotropic Kitaev terms |K|/2 and |K |/4. The difference in the two prefactors by a factor of 2 goes back to the value of ξ = 1/2 in the equivalent operator for the onsite quadrupolar terms. This difference is a quantum mechanical effect since, as we show below, in the classical treatment of the original spin-1/2 model, the direction in spin space is actually fixed by the competition between |K|/4 and |K |/4.

Remark 2:
If the couplings K and K were such that |K|/2 < |K |/4, λ x,y would be the lowest eigenvalues. Then the two solutions from the two M points would give T(na + mb) = e iM1·r e x = (−1) n+m e x , T(na + mb) = e iM1·r e y = (−1) n e y .
In this case we can in fact combine the two solutions to get a one-parameter family of degenerate states that contains coplanar as well as collinear states. The quantum order-by-disorder effect would most likely select the collinear states α = 0 or π, as they do in the original spin-1/2 model, see below.
Zigzag phase -finite C, D, C and D Diagonalizing the matrix Λ(k) we find that the minimum eigenvalue occurs again at the points M 1 and M 3 . The difference with the previous case is that now the principal eigenvectors differ from the (1, 0, 0), (0, 0, 1) or (0, 0, 1) directions. For example, for J 2 = J 3 = 3 meV, we obtain the following eigensolutions: This means that the spatial arrangement of the spins is that of the diagonal stripy phase but the overall direction in spin space is now along v 1 or v 3 , which are tilted away from the z axis: Supplementary Note 4: Back to the original spin-1/2 model: Classical limit of large −J > 0 and J 3 > 0 Here we discuss a semiclassical description of the diagonal stripy phase in the original spin-1/2 model. We begin by taking the limit J = K = K = J 2 = D = D = C = C = 0. The coupling J forces all pairs of spins on B1 bonds to be aligned, see Supplementary Figure 2(b). Together with the large antiferromagnetic (AF) coupling J 3 , these interactions give rise to two decoupled sublattices, denoted in Supplementary Figure 2(b) by blue arrows (±S b ) and red arrows (±S r ). Neither J nor J 3 can fix the relative orientation φ between blue and red spins. So, let us see what the remaining couplings do at the classical level. For simplicity we replace S b → b and S r → r.
• The contributions from the J couplings cancel out.
• The contributions from the J 2 couplings between blue and red cancel out as well.
• The contributions from the J 2 couplings between blue with blue and red with red spins give an overall constant.
• The energy contribution from the K and K couplings (where N stands for the number of B1 bonds) is • The energy contribution from D and D is • The energy contribution from C and C is Altogether, in spherical coordinates : The highly frustrated point C = D = C = D = 0: In this limit we have Since the energy does not depend on the angle difference φ b − φ r , this limit features a ground state manifold with an accidental one-parameter family of degenerate states. If we require that the solutions have the property θ b = θ r (confirmed by numerical minimization), then The extremum value of the second term occurs for | cos(φ b + φ r )| = 1. For cos(φ b + φ r ) = 1, we obtain and the minimum is at θ b = θ r = 0 with E 1 = −|K|. On the other hand, for cos(φ b + φ r ) = −1, we find and, since |K| < |K |, the minimum occurs for θ b = θ r = π/2 with E 2 = −|K |. Comparing E 1 and E 2 we find that the minimum energy is E 2 , i.e., when the spins lie within the xy-plane.
Quantum order-by-disorder -The above degeneracy of the classical ground state manifold is lifted by quantum fluctuations at the harmonic spin-wave level. This can be shown by performing a linear spin-wave expansion around the classical ground state manifold, which is parametrized by the angle α between the two sublattices. We follow the standard procedure and introduce two bosons a R,j=1,2 per unit cell of the honeycomb lattice, followed by a Holstein-Primakoff expansion of the spin operators. The final form of the bosonic Hamiltonian reads H S(S + 1) E clas S 2 + S 2 k f 1r (k)(a k,1 a −k,1 + a k,2 a −k,2 ) + f 2 (k)a k,1 a −k,2 + f 2 (−k)a −k,1 a k,2 + g 2 (k)a k,1 a + k,2 where , g 1 (k) = g 1 (k) + B clas and g 2 (k) = J + J (sin 2 α 2 e ik·c + cos 2 α 2 e −ik·b ) + K 2 − K 2 (cos 2 α 2 e ik·c + sin 2 α 2 e −ik·b ) .
We may now group all terms in matrix form as with , A + k = a + k,1 , a + k,2 , a −k,1 , a −k,2 .
We next define the commutator matrix and search for a transformation A k = S k ·Ã k that conserves the bosonic commutation relations (i.e.,g = g) and at the same time diagonalizes the Hamiltonian. Specifically, where Ω M k is diagonal. Ω M k is found using the eigenvalue equation One can show [17] that if M k is semidefinite positive then the eigenvectors of gM k come in pairs of opposite eigenvalues and lead to : The ground state of this Hamiltonian is the vacuum ofã k,1 ,ã k,2 , and so the quadratic energy correction is δE = S 2 k ω 1 (k) + ω 2 (k) .
Integrating over the BZ of the lattice, we find the quadratic energy correction as a function of the angle α between the two sublattices. Results are shown in Supplementary Figure 5 (energy per site) for the case of J 2 = J 3 = 3 meV and demonstrate that spin waves favor the collinear confingurations α = 0 or π.