Minimal models of $\alpha$-Li$_2$IrO$_3$: On the range of the interactions, ground state properties, and magnetization processes

TL;DR

This study investigates minimal theoretical models for the magnetic properties of $ ext{Li}_2 ext{IrO}_3$, identifying two key models that explain experimental observations and predicting their magnetization behaviors.

Contribution

The paper narrows down the range of minimal models for $ ext{Li}_2 ext{IrO}_3$ using Monte Carlo simulations, highlighting the importance of further neighbor interactions and anisotropies.

Findings

Two minimal models reproduce experimental results.

Strong bond anisotropies stabilize observed spiral states.

Predicted magnetization behaviors for future experiments.

Abstract

In recent years, a lot of effort has been devoted to the quest for experimental realizations of Kitaev interactions in spin systems. Recently, many materials have been synthesized which seem to realize extended Kitaev models, where Kitaev interactions are supplemented by Heisenberg and other bond dependent terms. The crystal and electronic structure of these materials renders the determination of a minimal model a non trivial pursuit. In this work, we will concentrate on one of these particular materials, $\alpha$-Li$_2$IrO$_3$, for which various minimal models have been proposed. Employing large scale Monte Carlo simulations we show how the number of prospective models can be reduced. We study in detail six models with different range of the interactions, and show how only two of those reproduce the most recent experimental results for this material. We obtain two possible minimal models, one of them with nearest neighbour interactions, while the other includes interactions up to third neighbours. Furthermore, we show that strong bond anisotropies and further neighbour interactions are crucial to stabilize the tilted counterotating spirals found in $\alpha$-Li$_2$IrO$_3$. We further clarify the picture, and distinguish these two models by studying the magnetization processes. We predict the magnetization behaviour of these models, and propose future experimental directions.

Tables

Table 1. Table 1: Values that reproduce the experimental results with the corresponding bond anisotropies for the anisotropic third neighbour model. All values are given in terms of | K 1 | subscript 𝐾 1 |K_{1}| .

	${\delta }_W$	Bond average W	Anisotropic interactions
$J_1$	0.14	0.2	$J_1^{XY}=0.06J_1^Z=0.34$
$K_1$	0.34	-1	$K_1^{XY}=-1.34K_1^Z=-0.66$
${\Gamma }_1$	0.195	1	${\Gamma }_1^{XY}=0.805{\Gamma }_1^Z=1.195$
$K_2$	0	-0.275	$K_2^{XY}=-0.275K_2^Z=-0.275$
${\Gamma }_2$	-0.06	0.275	${\Gamma }_2^{XY}=0.335{\Gamma }_2^Z=0.215$
$J_3$	0	0.3	$J_3^{XY}=0.3J_3^Z=0.3$

Table 2. Table 2: Values of the anisotropic interactions as obtained in Ref.  Winter  et al. ( 2016 ) . All interactions are given in meV. For the study of the anisotropic model we have selected those interactions which are greater than one (bold).

Bond	$J_n$	$K_n$	${\Gamma }_n$
$X_1$, $Y_1$	$-1.0$	$-13.0$	$+6.6$
$Z_1$	$-4.6$	$-4.2$	$+11.6$
$X_2$, $Y_2$	$+0.9$	$-2.9$	$+3.0$
$Z_2$	$-0.9$	$+0.1$	$+1.5$
$X_3$, $Y_3$	$+4.7$	$-0.2$	$0$
$Z_3$	$+4.4$	$+0.4$	$-0.1$

Table 3. Table 3: Values of the anisotropy constants for the anisotropic third neighbour model, extracted from Ref.  Winter  et al. ( 2016 ) following the process described in the text. All values are given in terms of | K 1 | subscript 𝐾 1 |K_{1}| .

	$\delta$
$J_1$	0.14
$K_1$	0.34
${\Gamma }_1$	0.195
$K_2$	0
${\Gamma }_2$	-0.06
$J_3$	0

Equations

$$\mathcal{H}=J\sum _{ij}{\mathbf{S}}_i{\mathbf{S}}_j+\sum _{ij\in \gamma -\mathrm{bonds}}\left(K{S}_i^{\gamma }{S}_j^{\gamma }+\Gamma (S_i^{\alpha }{S}_j^{\beta }+S_i^{\beta }{S}_j^{\alpha })\right),$$

$$\mathcal{H}(J_n,K_n,I_c,I_d)$$

$$\displaystyle+J_{n}\sum_{<ij>_{n}}\mathbf{S}_{i}\cdot\mathbf{S}_{j}+K_{n}\sum_{<ij>_{n}}\sum_{\gamma}S^{\gamma}_{i}S^{\gamma}_{j}\big{]}$$

$$\mathcal{H}(J_{1,3},K_{1,2},{\Gamma }_{1,2})=J_1\sum _{<ij{>}_1}{\mathbf{S}}_i\cdot {\mathbf{S}}_j+K_1\sum _{<ij{>}_1}{S}_i^{\gamma }{S}_j^{\gamma }$$

$$+{\Gamma }_1\sum _{<ij{>}_1}(S_i^{\alpha }{S}_j^{\beta }+S_i^{\beta }{S}_j^{\alpha })+K_2\sum _{<ij{>}_2}{S}_i^{\gamma }{S}_j^{\gamma }$$

$$+{\Gamma }_2\sum _{<ij{>}_2}(S_i^{\alpha }{S}_j^{\beta }+S_i^{\beta }{S}_j^{\alpha })+J_3\sum _{<ij{>}_3}{\mathbf{S}}_i\cdot {\mathbf{S}}_j,$$

$$J_1^{XY}=J_1-\delta J_1^Z=J_1+\delta$$

$$K_1^{XY}=K_1-\delta K_1^Z=K_1+\delta$$

$${\Gamma }_1^{XY}={\Gamma }_1-\delta {\Gamma }_1^Z={\Gamma }_1+\delta$$

$${\Gamma }_2^{XY}={\Gamma }_1-\delta {\Gamma }_1^Z={\Gamma }_1+\delta ,$$

$$C_{ij}=⟨{\mathbf{S}}_i\cdot {\mathbf{S}}_j⟩-⟨{\mathbf{S}}_i⟩⟨{\mathbf{S}}_j⟩,$$

$$|{\mathbf{OP}}_{FM}|=\frac{1}{N}|\sum _i⟨{\mathbf{S}}_i⟩|,$$

$$O{P}_{M\cdot H}=\frac{1}{N}(O{P}_{FMx}{H}_x+O{P}_{FMy}{H}_y+O{P}_{FMz}{H}_z),$$

$$\chi =\frac{d{\mathbf{FM}}_{OP}}{dH}=\frac{1}{N}\frac{1}{T}\sum _i\left(⟨{\mathbf{S}}_i^2⟩-⟨{\mathbf{S}}_i{⟩}^2\right)$$

$${\alpha }_{XY}={\alpha }_m-\delta {\alpha }_Z={\alpha }_m+\delta$$

$$I_c\sum _{<ij>}{S}_i^{r_{ij}}{S}_j^{r_{ij}}$$

$$=I_c(S_i^x{S}_j^x+S_i^y{S}_j^y)+I_c{\Gamma }^{xy},$$

$${\Theta }_{xx}={\Theta }_{yy}=\frac{S(S+1)}{3K_b}(-3J-I_c-I_d-K),$$

$${\Theta }_{zz}=\frac{S(S+1)}{3K_b}(-3J-2I_d-K).$$

$${\Theta }_{xx}={\Theta }_{yy}=\frac{S(S+1)}{3K_b}(-3J_1^{XY}-K_1^{XY}-2K_2-3J_3),$$

$${\Theta }_{zz}=\frac{S(S+1)}{3K_b}(-3J_1^Z-K_1^Z-2K_2-3J_3),$$

$$\mathcal{H}=J\sum _{<ij>}{\mathbf{S}}_i\cdot {\mathbf{S}}_j+K\sum _{<ij>}\sum _{\gamma }{S}_i^{\gamma }{S}_j^{\gamma },$$

$$\mathcal{H}(J_1,K_1,I_c\ne 0,I_d=0)$$

$$+I_c\sum _{<ij{>}_{zz}}{S}_i^{r_{ij}}{S}_j^{r_{ij}},$$

$$\mathcal{H}(J_{1,2},$$

$$+K_1\sum _{<ij>}\sum _{\gamma }{S}_i^{\gamma }{S}_j^{\gamma }+J_2\sum _{<<ij>>}{\mathbf{S}}_i\cdot {\mathbf{S}}_j$$

$$+K_2\sum _{<<ij>>}\sum _{\gamma }{S}_i^{\gamma }{S}_j^{\gamma }$$

$$\mathcal{H}=$$

$$+K_2\sum _{<<ij>>}\sum _{\gamma }{S}_i^{\gamma }{S}_j^{\gamma }+I_c\sum _{<ij>}{S}_i^{r_{ij}}{S}_j^{r_{ij}}$$

Full text

I Introduction

The theoretical study of Kitaev materials is a flourishing part of the frustrated magnetism field Jackeli and Khaliullin (2009); Chaloupka et al. (2010); Song et al. (2016); Lee et al. (2016); Winter et al. (2017). Since Kitaev proposed his exactly solvable honeycomb model Kitaev (2006), which belongs to a wider range of quantum compass Hamiltonians Nussinov and van den Brink (2015), a great number of resources have been dedicated at studying this system. Kitaev’s proposal kick-started a quest for experimental realizations of Kitaev spin liquids, and since then many materials have been synthesized that have been proposed to realize Kitaev interactions ($\mathrm{Na}_{2}\mathrm{IrO}_{3}$Singh and Gegenwart (2010), $\mathrm{Li}_{2}\mathrm{IrO}_{3}$Singh et al. (2012); Takayama et al. (2015); Modic et al. (2014), $\mathrm{Ru}\mathrm{Cl}_{3}$Plumb et al. (2014), $\mathrm{Ba}_{3}\mathrm{IrTi}_{2}\mathrm{O}_{9}$Dey et al. (2012), among others).

Until now, though, these materials not only present Kitaev interactions but are also accompanied by Heisenberg and other bond dependent exchangesJackeli and Khaliullin (2009); Chaloupka et al. (2010); Rau et al. (2014); Winter et al. (2017), the so called extended Kitaev models. These perturbing interactions have a strong effect on the ground state of the these materials. Furthermore, even though Kitaev interactions give rise to spin liquid states, the perturbations are so strong that the synthesized materials exhibit classical magnetic order.

While the quest for materials which only realize Kitaev interactions is still ongoingNasu et al. (2016); Kitagawa et al. (2018), the theoretical understanding of extended Kitaev models remains challenging. The complexity of their crystal structure and orbital hybridization Winter et al. (2017) has provided significant challenges to their modelling. So far, we know that the minimal models for these materials can be built as modifications of the same basic HamiltonianJackeli and Khaliullin (2009); Chaloupka et al. (2010); Rau et al. (2014); Winter et al. (2017), but where different variations of the model seem to be realized in different materials Singh and Gegenwart (2010); Singh et al. (2012); Takayama et al. (2015); Modic et al. (2014); Plumb et al. (2014); Winter et al. (2017), as for example in $\mathrm{Na}_{2}\mathrm{IrO}_{3}$Chaloupka et al. (2010); Singh et al. (2012); Choi et al. (2012); Yamaji et al. (2014); Hwan Chun et al. (2015); Winter et al. (2016) and $\alpha-\mathrm{Li}_{2}\mathrm{IrO}_{3}$ Singh et al. (2012); Biffin et al. (2014); Reuther et al. (2014); Kimchi et al. (2015); Williams et al. (2016); Winter et al. (2017) . In particular, different materials present different bond anisotropies on the interactions and further neighbour exchanges. Since different variations of the same model can reproduce experimental results, the question that arises is related to how can one distinguish between these models further, and which experiments can be performed to asses their validity. In this work, we will attempt to answer this question in the context of $\alpha$-$\mathrm{Li}_{2}\mathrm{IrO}_{3}$.

$\alpha-\mathrm{Li}_{2}\mathrm{IrO}_{3}$ (isostructural to $\mathrm{Na}_{2}\mathrm{IrO}_{3}$) possesses a layered crystal structure where the Ir4+ ions, surrounded by an octahedral cage of oxygens, form a honeycomb lattice Freund et al. (2016). This lattice resides in the $[111]$ Cartesian plane. Both compounds, $\mathrm{Na}_{2}\mathrm{IrO}_{3}$ and $\alpha-\mathrm{Li}_{2}\mathrm{IrO}_{3}$, present long range order, showing anomalies in the specific heat and in the magnetic susceptibility at a critical temperature $T_{c}\sim 15K$, with a Curie-Weiss temperature of $\Theta=-125(6)K$ and $\Theta=-33(3)K$ respectively Singh et al. (2012). $\mathrm{Na}_{2}\mathrm{IrO}_{3}$ presents a zig-zag ordered state and it would be expected that, given the similarities in the thermodynamic anomalies, $\alpha-\mathrm{Li}_{2}\mathrm{IrO}_{3}$ would also present a zig-zag order. While the thermodynamics show similarities in the behavior of the Na and Li compound, recent studies performed on single crystals and powder samples of $\alpha-\mathrm{Li}_{2}\mathrm{IrO}_{3}$ have shown that the magnetic order present in this material is not of zig-zag type, but of spiral nature. Magnetic resonant X-ray diffraction (MRXD) together with magnetic powder neutron diffraction determined a magnetic structure composed of counterrotating incommensurate coplanar spin spirals Williams et al. (2016), with a propagation wavevector $\mathbf{q}=(0.32(1),0,0)$. At the same time, the authors of Ref. Williams et al. (2016) were able to determine that the plane of rotation of the spirals is uniform between the different sublattices and tilted with respect to the lattice plane by $80^{o}$. In their studies they see that the MRXD, at a temperature of $5K$, presents satellite peaks at positions $\mathbf{\tau}\pm\mathbf{q}$, where $\mathbf{\tau}$ are the positions of allowed structural reflections $\mathbf{\tau}=(h,k,l)$, with $h+k=$even.

Similarities have been found between the ground state of $\alpha-\mathrm{Li}_{2}\mathrm{IrO}_{3}$ and the two structural polytypes, $\beta-\mathrm{Li}_{2}\mathrm{IrO}_{3}$ and $\gamma-\mathrm{Li}_{2}\mathrm{IrO}_{3}$, which correspond to hyper-honeycomb and stripy-honeycomb magnetic lattices respectively. All three of these polytypes are members of the “harmonic honeycomb” structural series Modic et al. (2014), and the similarities in their magnetic ordering have lead to proposals of universality between the members of the family of harmonic honeycomb Iridates Kimchi et al. (2015).

It has been shown previously that the magnetic ordering of the $\beta$ and $\gamma$ structures are well described by a dominant ferromagnetic Kitaev interaction, combined with other smaller exchange terms of the form of Heisenberg and bond dependent terms, Biffin et al. (2014); Lee and Kim (2015); Kimchi et al. (2015); Lee et al. (2016).

Given the similarities between the different polytypes, it is expected that the model which correctly describes the magnetic structure of $\alpha$–$\mathrm{Li}_{2}\mathrm{IrO}_{3}$ will be composed of dominant Kitaev interactions, supplemented by weaker Heisenberg and bond dependent terms. While in principle this proposal is in agreement with what is known for the $\beta$ and $\gamma$ polytypes, the particularities of the minimal model have not yet been defined. Many models have been proposed which reproduce some characteristics of the $\alpha$ polytype, but which have been tested only on toy models, or by a Luttinger-Tisza (LT) approximation, which are not ensured to succeed at the detection of incommensurate states.

Regarding $\alpha$–$\mathrm{Li}_{2}\mathrm{IrO}_{3}$, the range of its interactions is still an open question. There have been proposals where this material is modeled as a nearest neighbour Singh et al. (2012); Kimchi et al. (2015); Williams et al. (2016), second neighbour Reuther et al. (2014), or third neighbour Winter et al. (2017) extended Kitaev model, and in this work we perform a comparative analysis between all those models, via Monte Carlo simulations. We aim at exploring the possible ground states of the different models via a method which is not biased towards any type of magnetic order, and which can deal with the incommensurate nature of the expected ground state appropriately. We find that among the six proposed models up to date, only two, a nearest and a third neighbour model, can reproduce the experimental features correctly. The nearest neighbour model we study is a variation of the model proposed by I. Kimchi et. al. Kimchi et al. (2015); Williams et al. (2016). In the model of Ref.Kimchi et al. (2015), bond dependent interaction are allowed by the crystallographic symmetry, and was studied by the authors in a toy model consisting of 1D zig-zag chains. We will extend this model to include further bond dependent interactions, and we will study it employing large scale Monte Carlo simulations. On the other hand, we also find that the third neighbour model proposed by Winter et. al. via DFT calculationsWinter et al. (2017) can also accurately reproduce the experimental results, but provided that bond anisotropies are included.

The discrepancy between the range of the interactions of these two models raises the question of how we can further distinguish them to ascertain which one corresponds to $\alpha$–$\mathrm{Li}_{2}\mathrm{IrO}_{3}$. We propose to answer this question studying the magnetization behaviour of both models. Our aim is that the differences in the magnetization processes with different applied field directions can provide an efficient experimental route to probe which of the proposed models is more feasible. In the presence of a magnetic field, we find that both models present different magnetization processes depending on the direction in which the field is applied. This indicates an experimentally realizable way of determining if one of the proposed models is correct, by studying the low temperature magnetization of single crystals of $\alpha$–$\mathrm{Li}_{2}\mathrm{IrO}_{3}$.

The paper is structured as follows: in section II we show the fundamental characteristics of the studied models, and briefly describe the numerical method employed. In section III we show the ground state properties for those models which reproduce the experimental results, and in section IV we study in detail the magnetization processes for those same models. Finally in section V we discuss our results in the context of the available experimental evidence, and propose further theoretical and experimental studies which can help determine which minimal model corresponds to $\alpha$–$\mathrm{Li}_{2}\mathrm{IrO}_{3}$.

II Models and method

II.1 Models

We classify the different models according to the range of their interactions, into nearest, second, and third neighbour models. While here we mention the common characteristics of all of them, in the main text we will only show the results for those models which reproduce the experimental features of $\alpha-\mathrm{Li}_{2}\mathrm{IrO}_{3}$. In Appendices B, C, D, and E we show the ground state properties for the rest of the models.

The models treated here are variations of the extended Kitaev model

[TABLE]

where $\gamma=\{x,y,z\}$, and where $\alpha$ and $\beta$ indicate the two spin components perpendicular to $\gamma$. $J$ represents a Heisenberg coupling, while the bond dependent terms contain Kitaev interactions coupled by $K$. The $\Gamma$ exchange couples two orthogonal spin components, $\alpha$ and $\beta$, along the bond with Kitaev interactions in the $\gamma$ spin component. A diagram of the lattice and the Kitaev interactions is shown in Fig. 1.

The nearest ($\mathcal{H}(J_{1},K_{1},I_{c},I_{d})$) and second neighbor ($\mathcal{H}(J_{1,2},K_{1,2},I_{c},I_{d})$) models we study are an anisotropic version of Eq.1, based on the model proposed in Refs.  Kimchi et al. (2015); Williams et al. (2016):

[TABLE]

where $<ij>_{n}$ denote a sum over $n$-th neighbours, and $\mathbf{S}_{i}$ = $(S^{x}_{i},S^{y}_{i},S^{z}_{i})$ denotes the classical spin acting on site $i$. The terms containing the couplings $I_{c}$ and $I_{d}$ are Ising terms that couple the spins components parallel to the bond orientation, i.e $\mathbf{S}^{r_{ij}}=\mathbf{S}\cdot\mathbf{\hat{r}}_{ij}$, where $\mathbf{\hat{r}}_{ij}$ is the unit vector connecting the spins at sites $i$ and $j$. Please note that the Kitaev model has a particular symmetry in the bond isotropic case, where a $60^{o}$ rotation in real and spin space leave the ground state invariant, and this symmetry is preserved in the model of Eq. 2 when $I_{c}=I_{d}$. In the real material, the octahedral cage enclosing the Ir4+ atoms is not perfect, presenting deformations. These deformations induce a bond anisotropy on the interactions, where the couplings of the Ising terms is not the same on the $zz$-, $xx$- and $yy$-bonds. For this, we choose $I_{c}$ to be active only on the zz-bonds, while $I_{d}$ acts on the rest of them (zig-zag bonds).

We studied two nearest neighbour models, the model given by $\mathcal{H}(J_{1},K_{1},I_{c}\neq 0,I_{d}=0)$ contains dominant ferromagnetic Kitaev interactions $K_{1}<0$ as well as small antiferromagnetic Heisenberg terms $J_{1}>0$, the Ising term is ferromagnetic as well, $I_{c}<0$. The second nearest neighbour model also includes terms with $I_{d}<0$, $\mathcal{H}(J_{1},K_{1},I_{c}\neq 0,I_{d}\neq 0)$. For the second neighbour models we studied an isotropic model, $\mathcal{H}(J_{1,2},K_{1,2},I_{c}=0,I_{d}=0)$, with nearest and second neighbour Heisenberg interactions, $J_{1}>0$ and $J_{2}<0$ respectively, as well as nearest and second neighbour Kitaev exchanges, $K_{1}<0$ and $K_{2}>0$. Finally we study an anisotropic second neighbour model, $\mathcal{H}(J_{1,2},K_{1,2},I_{c}\neq 0,I_{d}=0)$, including terms with $I_{c}<0$.

For the third neighbor case, $\mathcal{H}(J_{1,3},K_{1,2},\Gamma_{1,2})$, we implemented the model proposed by Winter et.al. based on DFT calculations Winter et al. (2016). In their work, they propose a model with Heisenberg, Kitaev, and other bond dependent interactions, some of them ranging up to third neighbours. The effective Hamiltonian obtained results in

[TABLE]

where $X_{n}$, with $X=J$, $K$, or $\Gamma$, represent the exchange coupling for an interaction between $n$th-neighbours, $\sum_{<ij>_{n}}$ represents a sum over $n$th-neighbours, and $\{\alpha,\beta,\gamma\}=\{x,y,z\}$ indicate the spin component. A particularity of this model is that the Kitaev and $\Gamma$ exchanges are equal in magnitude but opposite in sign, $|K_{n}|=-|\Gamma_{n}|$.

Ref. Winter et al. (2016) reports strong bond anisotropies in the model, an observation which coincides with previous proposals for this material Kimchi et al. (2015). To gain an understanding of such complex Hamiltonian we first studied the bond isotropic case $\mathcal{H}_{I}(J_{1,3},K_{1,2},\Gamma_{1,2})$, to afterwards introduce anisotropies, giving rise to the anisotropic third neighbour model $\mathcal{H}_{A}(J_{1,3},K_{1,2},\Gamma_{1,2})$. The following bond anisotropies based on Ref. Winter et al. (2016) are introduced,

[TABLE]

where the quantities $J_{1}$, $K_{1}$, $\Gamma_{1}$, and $\Gamma_{2}$ are understood as bond averaged interactions, and the superscript indicate on which bond these interactions act.

In Table 1 we show the bond average value of the different couplings together with the anisotropic component $\delta$ that reproduce the experimental results.

II.2 Method

The numerical solution consist on Monte Carlo simulations for classical Heisenberg spins implementing a Metropolis-Hastings algorithm. Even though we want to study ground state properties of the bulk, some of these ground states are incommensurate phases which, to the effect of the algorithm, means that we implement free edge boundary conditions (FEBs), where the spins on the edges are exposed to a putative vacuum. The presence of FEBs carries some added effects, and as such the system will exhibit edge modes that are, in principle, not relevant to the study of the bulk physics. We show in Appendix A that the edge modes do not affect the bulk properties of the data, provided the lattice sizes are sufficiently big. We achieve this via a benchmark of our code implementing FEBs against the results of Ref. Price and Perkins (2013) for the Kitaev-Heisenberg model ($J_{1}\neq 0$, $K_{1}\neq 0$) with periodic boundary conditions. To minimize the finite size effects and further enhance equilibration we implement an iterative minimization and parallel tempering algorithms respectively, with system sizes ranging from 2400 up to 5400 sites.

To identify the different states in the phase diagram we will rely on the study of the Fourier transform of the spin-spin correlation function

[TABLE]

where the average is taken over different Monte Carlo sweeps. For the ground state of $\alpha$–$\mathrm{Li}_{2}\mathrm{IrO}_{3}$ we expect to find peaks at positions indicated by the red stars in Fig. 2 in accordance with the results of Ref.Williams et al. (2016). We can determine the propagation wavevector of the magnetic spirals from the maxima of the Fourier transform of Eq. 5, as they will be located at the points $\tau\pm\mathbf{q}$, with $\mathbf{q}=(\pm q,0)$ the propagation wavevector. At the same time, and since Monte Carlo produces the spin pattern of the state, we confirm our results by calculating the relative angle between nearest neighbour spins in one spiral.

To study the magnetization processes we will be interested in the ferromagnetic order parameter,

[TABLE]

where the average is taken over different Monte Carlo sweeps, and the magnetization in the direction of an externally applied field,

[TABLE]

where N is the number of sites, $\mathbf{S}_{i}$ represent the spin at site $i$, and $OP_{FMi}$ is the $i$-th component of the ferromagnetic order parameter. To probe the existence of a ferromagnetic order with a particular polarization the ferromagnetic order parameter has to be projected onto this direction. As we will include magnetic fields, we choose to study the projection of $\mathbf{OP}_{FM}$ in the direction of the applied field. To estimate critical fields we will also employ the associated response function of the ferromagnetic order parameter, the magnetic susceptibility

[TABLE]

where T is the temperature of the simulation.

When everything is taken into consideration, our simulations consist of system sizes ranging from 2400 to 5400 Heisenberg spins, with a temperature consistently set at $T=0.001$ in units of the dominant coupling of the model. We employ $2$x$10^{6}$ Monte Carlo sweeps, of which $10^{5}$ are used as equilibration steps, and the rest are employed to calculate observables. Furthermore, we average over 10 independent runs per data point. We interlace these sweeps with parallel tempering swaps, for 32 replicas, every 100 sweeps. Once the low temperature state was equilibrated we employ an iterative minimization algorithm with a threshold $10^{-20}$, to minimize the energy further.

III Ground state properties

III.1 Nearest neighbour model, $\mathcal{H}(J_{1},K_{1},I_{c}\neq 0,I_{d}\neq 0)$

In the following, we will show the results for the ground state properties of the model described by $\mathcal{H}(J_{1},K_{1},I_{c}\neq 0,I_{d}\neq 0)$ (Eq. 2). The limit $I_{d}=0$ corresponds to the $\mathcal{H}(J_{1},K_{1},I_{c}\neq 0,I_{d}=0)$ model (Eq. 2), which we show in Appendix B.

Phase diagram

We study the model described in Eq. 2 employing as a starting point the suggested values of the exchange couplings given in Ref. Williams et al. (2016). We will set fixed values for the Heisenberg and Kitaev interactions $J/|K|=0.2$, $K/|K|=-1$. To analyse the effect of bond dependent interactions we allow the couplings $I_{c}/|K|$ and $I_{d}/|K|$ to move in the range $\{-1,...0\}$. Studying both the real space configuration of the spins, as well as the correlation function in reciprocal space we are able to map the phase diagram for finite $I_{c}$ and $I_{d}$. We show this phase diagram in Fig. 3(right), where we observe a rich behaviour, with different commensurate-incommensurate transitions. The green squares corresponds to the regions of the phase diagram where an incommensurate spin spiral state is found.

The regime $I_{c}=0$ and for moderate values of $I_{d}$ ($I_{d}>-0.25$) we observe a degenerate stripy phase. This state arise from the presence of magnetic domains exhibiting stripy phases polarized in the X (st-X phase) and Y (st-Y phase) directions. We denote this state by st-XY (orange triangles). For smaller $I_{d}$ we obtain a distorted $120^{o}$ order (red diamonds). This latter state was previously studied via a soft spin approximationWilliams et al. (2016), and identified as an incommensurate spin spiral propagating in the vertical direction of Fig. 1. The Fourier transform of the correlation function for this phase presents maxima close to the first Brillouin zone, Fig. 3(left), while the spin pattern in real space shows that the state is that of a distorted $120^{o}$. For small enough $I_{d}$, the distortion in the $120^{o}$ is more pronounced, and we expect that in the limit $I_{d}\rightarrow-\infty$ we recover the perfect $120^{o}$ order.

The degeneracy shown in the st-XY phase is expected since the Heisenberg-Kitaev model stabilizes a triple degenerated stripy phase Price and Perkins (2013). The inclusion of a small bond dependent interaction strengthening one particular bond, breaks the degeneracy of the stripy phase. When both $I_{c}$ and $I_{d}$ are non-zero, we see a clear separation of phases through the line $I_{c}=I_{d}$. When $|I_{c}|>|I_{d}|$ a counterrotating spiral state dominates the phase diagram. This phase reproduces the experimental results for $\alpha-\mathrm{Li}_{2}\mathrm{IrO}_{3}$ and will be studied in detail in the next section.

Finally, we note that, at high enough values of $|I_{c}|$ and $|I_{d}|$ the model adopts a ferromagnetic order with polarization in the lattice plane. This order lives on both sides of the $I_{c}=I_{d}$ line, presenting a ferromagnetic order with a net magnetization in the direction of the $zz$-bonds when $|I_{c}|>|I_{d}|$, and perpendicular to it in the $|I_{c}|<|I_{d}|$ case. The line in the phase diagram separating both polarizations is special given that the low energy fluctuations in this regime are of a different nature. In this ferromagnetic state, and over the line $I_{c}=I_{d}$, the low energy fluctuations are vortex-like, and appear in pairs of vortex-antivortex fluctuations. An example of a spin pattern presenting this behaviour is shown, for a calculation over 2400 sites, in Fig. 4.

These vortex fluctuations where studied in systems ranging from 24 to 5400 lattice sites, and they consistently appear in all the studied system sizes over this particular line in the phase diagram, which rules out this behaviour as a finite size effect. A detailed study of these vortices is beyond the scope of this paper, and is left for future study.

Spiral properties

A big part of the phase diagram on Fig. 3 is dominated by an incommensurate phase (green squares). This state represents an incommensurate counterrotating spiral which propagates in the horizontal direction according to Fig. 1 (the direction perpendicular to the zz-bonds). In the regime $I_{d}=0$ the wavevector varies between $0.5$ and $0.4$ (in units of $2\pi$) and the plane of rotation is tilted with respect to the lattice plane by $54^{o}$, i.e, the rotation plane is oriented parallel to the XY-Cartesian plane.

In the regime where both $I_{c}\neq 0$ and $I_{d}\neq 0$ the commensurate phases survive for values of $I_{d}$ down to $-0.4$ and for values of $I_{c}$ such that $|I_{c}|>|I_{d}|-0.4$ for $|I_{d}|<0.2$, and $|I_{c}|>|I_{d}|$ for $0.2<|I_{d}|<0.4$. In this regime, some properties of the spiral phase found for the $I_{d}=0$ case are modified. As $I_{c}$ and $I_{d}$ are varied, the wavevector varies between $0.5$ and [math] (wavevector [math] correspond to the onset of ferromagnetic order). The rotation plane’s tilt now also varies along the phase diagram, between $\sim 50^{o}$ (consistent with the spin spiral known to appear at $I_{d}=0$) and $\sim 90^{o}$. We show in Fig. 5 a real space pattern of the spin spiral, at a point in the phase diagram which reproduces the experimental results.

In Fig.5(right)we show two heat maps, one for the wavevector and another for the tilt angle of the rotation plane, respectively. We can see the variation of these quantities in the phase diagram. In the heat map corresponding to the wavevector (Fig. 5(top right)) we indicate with a black line the zone of wavevector $q=0.32$ (in units of $2\pi$) and superimpose this line over the tilt angle heat map (Fig. 5(bottom right)). Please note that this mark is a guide to the eye, it does not arise from a fit to the data.

Further confirming the existence of counterrotating spin spirals reproducing the experimental results, we show the Fourier transform of the correlation function in Fig. 6. This figure was calculated for the parameters $I_{c}=-0.5$ and $I_{d}=-0.35$, which reproduce the wavevector and rotation plane tilt found in experiments. The correlation function presents maxima at the expected positions. They appear as satellites of the $\Gamma$ point, with positions $\tau+\mathbf{q}$, where $\tau$ indicates the location of the $\Gamma$ points and $\mathbf{q}=(\pm 0.32,0)$ in units of $2\pi$, which coincide with the experimental results mentioned in the introduction.

III.2 Anisotropic third neighbour model, $\mathcal{H}_{A}(J_{1,3},K_{1,2},\Gamma_{1,2})$

The third neighbour model, $\mathcal{H}(J_{1,3},K_{1,2},\Gamma_{1,2})$, was proposed by Winter et. al. in Ref. Winter et al. (2016) where it was determined (via a combination of DFT and exact diagonalization on a small cluster of hexagons) that this model presents large anisotropies in various parameters. The isotropic model, $\mathcal{H}_{I}(J_{1,3},K_{1,2},\Gamma_{1,2})$, presents incommensurate spin spirals in the ground state (results for the isotropic model are shown in Appendix. C). However, since the model is bond isotropic, the spin spirals can propagate in three symmetry allowed directions. The powder MRXD measurements performed on $\alpha$-Li2IrO3 indicate that the material does not seem to show any type of degeneracy of the ground state. While this degeneracy could be broken by order by disorder effects, the fact that the model proposed by Winter et. al. Winter et al. (2016) exhibits strong bond anisotropies seem to indicate that these anisotropies need to be included in the model to break the degeneracy. In the following we will treat the model shown in Eq. 3 including these anisotropies and study what their effect is on the phase diagram and the spiral properties of the system.

The values of the exchange parameters, as obtained in Ref.Winter et al. (2016), are shown in table 2. Studying this table it becomes clear that an analysis considering that the third neighbour model is bond isotropic is an excessive simplification. We will introduce the anisotropies in the model as shown in Eq. 3 in the following way: for a given coupling $\alpha$ we have a bond anisotropy which differentiates the X and Y bond ($\alpha_{XY}$) from the Z bond ($\alpha_{Z}$). We will define $\alpha_{m}$ as the bond average of the exchange coupling, $\alpha_{XY}+\alpha_{z}/2$, and $\delta$ as the anisotropy constant (with an appropriate sign) such that (for a direct comparison please look at Eqs. 4)

[TABLE]

This way, calculating the bond average from Table 2 we can determine what the anisotropy for each exchange is. Please note that after this process is performed all couplings and anisotropy constants are expressed in terms of $|K^{XY}_{1}|=13$. The values obtained for the anisotropy constants are given in Table. 3.

When we study the bond isotropic model with the exchange constants obtained from Table. 2 we obtain a zig-zag order, which agrees with the results from Ref. Winter et al. (2016). In the bond anisotropic case, where the interaction couplings reduce to those shown in Table. 2, we again obtain a zig-zag state 111Please note that in the phase diagrams shown in this paper no zig-zag order is shown, as we have concentrated in mapping a part of the diagram which exhibits incommensurate spirals. If we were to map the full phase diagram we would see (and we have confirmed this via numerical simulations) that, as in Ref. Winter et al. (2016) for $J_{3}\gtrsim 0.4$ a zig-zag phase is present.. By changing the bond averages but maintaining the anisotropic parameters $\delta$ constant, we can map a phase diagram including anisotropies. Changing the bond averages is not a radical idea, since the values obtained by Winter et. al. Winter et al. (2016) have been obtained via exact diagonalization on small clusters. This, combined with the uncertainty in the crystal structure which has been resolved until this point indicates that, while the nature of the interactions might not change, their coupling strength could.

Phase diagram

The phase diagram of the anisotropic third neighbour model is shown in Fig. 7 top. We map the phase diagram for the following coupling strengths: a dominant nearest neighbour Kitaev coupling $K_{1}=-1$ supported by (here and in the following, all exchange couplings are given in units of $|K_{1}|$) $K_{2}=-0.275$, $\Gamma_{1}=1$, $\Gamma_{2}=0.275$, $J_{1}\in(-0.1,...,0.4)$, and $J_{3}\in(0,...,0.4)$. The phase diagram presents two dominant phases, a ferromagnetic state (black dots) and an incommensurate state (green squares). At the bottom right corner a small stripy phase is observed (blue squares). A comparison with the phase diagram for the isotropic model (Fig. 25) indicates that both the ferromagnetic and the spin spiral states are displacing the rest of the phases.

The ferromagnetic phase present in this model exhibits an in plane net magnetization, in the direction parallel to the $zz$-bonds. Domain walls separate two domains exhibiting the two possible orientations of the polarization, where these domain walls are realized by spins aligned antiferromagnetically. As $J_{1}$ and $J_{3}$ are increased the domains multiply, until the system enters an incommensurate phase.

This commensurate-incommensurate transition can be seen as the generation of magnetic domains overpowering this ferromagnetic state. In Fig. 8 we show a sketch of the projected spin pattern on the [111] plane (the lattice plane) for a cut through a fixed value of $J_{3}$, where the system transitions from a ferromagnetic to a spiral state. To simplify the argument we have assumed for this discussion that the angles between the spins in the resulting spirals are $45^{o}$, and there is no tilt angle in the rotation plane 222The reason for choosing this particular example is that of convenience, given that at this angle and plane tilt, the spin spiral can be seen as alternating ferromagnetic/antiferromagnetic domains of size two. For another wavevector there would also be ferromagnetic and antiferromagnetic domains but the size of these regions would not be the same.. For this particular example, at small values of $J_{1}$ the state presents no domain walls (Fig. 8 (a)). As $J_{1}$ is increased domain walls start to span the length of the system, separating big domains of ferromagnetic order (Fig. 8 (b)) with opposing polarization vectors. For even bigger $J_{1}$ the system now contains ferromagnetic and antiferromagnetic domains spanning two sites in the vertical direction, each. This is seen in Fig. 8 (c), which corresponds to spin spirals of a wavevector such that the angle between spins is $45^{o}$.

Spiral properties

As can be seen from Fig. 7, a big part of the phase diagram is dominated by an incommensurate state. In Fig. 9 we show the Fourier transform of the spin pattern obtained from this model (the real space spin pattern is shown in Fig. 5(left)). From this we observe that the maxima are located inside the Brillouin zone as satellites of the $\Gamma$ point, and secondary maxima appear as satellite points of the $K$ points. From the Fourier transform we can further extract a wavevector $(q,0)$, $q=0.32$ in units of $2\pi$, which coincides with the analysis performed on the real space spin pattern which we show below. Overall, this is consistent with the MRXD results.

The incommensurate state are spin spirals of the same nature as those found for the nearest neighbour model studied in the previous section. In essence, the state is such that planar spirals propagate in the direction perpendicular to the $zz$-bonds. Furthermore, the spirals counterotate, with those formed by spins on sublattice $a$ and $c$ rotating with opposite chirality to those formed by the spins of sublattices $b$ and $d$. A real space pattern of the spirals can be observed in Fig. 5(left), where we show a spiral which coincides with the experimental results for $\alpha$-Li2IrO3, exhibiting a wavevector $(q,0)$ where $q=0.32$ in units of $2\pi$, and a tilt of the rotation plane of $\sim 80^{0}$.

Comparing this spiral with that found for the bond isotropic model, it is not surprising to notice that the effect of the anisotropies was that of destroying the degeneracy of the spiral phase. Recall that in the bond isotropic model our spirals were planar spirals, but also degenerate, this degeneracy arising from the fact that the Hamiltonian retains the discrete Kitaev symmetry. Thus the spirals were free to propagate in three possible directions.

Throughout the phase diagram, we notice that the wavevector and tilt angle change. In Fig. 10 we show a heat map of the variation of the wavevector (Fig. 10(top)) and of the tilt angle (Fig. 10(bottom)) as the exchanges $J_{1}$ and $J_{3}$ are modified.

IV Magnetization processes

Up to this point, we have performed an analysis of the possible minimal models for $\alpha$-Li2IrO3. All these models where proposed either on the basis of symmetry allowed interactions Singh et al. (2012); Rau et al. (2016, 2014); Khaliullin (2005); Jackeli and Khaliullin (2009), or on DFT studies Winter et al. (2016, 2017). Comparing the results for the models shown previously, and those shown in the appendices, we have determined that only two models reproduce the experimental features of the material. While all models present incommensurate spin spiral phases, only the nearest neighbour ($\mathcal{H}(J_{1},K_{1},I_{c}\neq 0,I_{d}\neq 0)$) and anisotropic third neighbour models reproduce not only the counterotating nature and wavevector of the spirals, but also the tilt of the rotation plane.

Now we want to put forward a prediction regarding the magnetization behaviour of these models in a way which can be verifiable experimentally. The fact that the models which reproduce the experimental results have different range of interactions mean that more studies are needed to reduce the number of possible models further. For this we have chosen to study the magnetization processes of the different models by applying an external magnetic field in different directions. Given the different interactions and bond anisotropy of both models, the behaviour will be different for both models, depending on the direction of the applied field. As we will see, the magnetization processes of the nearest neighbour and anisotropic third neighbour models are radically different, which indicates an experimentally feasible way of probing whether one of this models is in fact a good representation of $\alpha$-Li2IrO3.

IV.1 Nearest neighbour model, $\mathcal{H}(J_{1},K_{1},I_{c}\neq 0,I_{d}\neq 0)$

We will study the magnetization processes of the model given by Eq. 2, for magnetic fields applied in three different global directions: $H\parallel[111]$, $H\parallel[\bar{1}10]$, and $H\parallel[11\bar{2}]$. Of these three directions, $[111]$ corresponds to the direction perpendicular to the lattice plane, while $[\bar{1}10]$ ($[11\bar{2}]$) is the direction parallel (perpendicular) to the $zz$- bonds (see Fig. 1).

We will concentrate on the point of the phase diagram that reproduces the experimental results: $I_{c}=-0.5$ and $I_{d}=-0.35$ (in units of the Kitaev coupling $|K|$). Fig. 11 shows the magnetization curves for the different directions of the field. While the magnetization process is consistently the same for all field directions, ferromagnetic domains are created, which increase in size until they dominate the system, we observe different critical fields for the different field directions. For field directions in the lattice plane ($H\parallel[\bar{1}10]$, $H\parallel[11\bar{2}]$) the critical field is much lower than for the out of plane direction ($H\parallel[111]$). This indicates the existence of an easy plane anisotropy.

The origin of this easy plane can be seen from the position of the experimentally relevant point in the phase diagram, as this point is close to the phase boundary between the spiral and the ferromagnetic phase. We have seen before that the ferromagnetic phase is such that the polarization vector is in the lattice plane, in a direction that depends on the values of $I_{c}$ and $I_{d}$. Applying a magnetic field in the direction of one of the bonds then is equivalent to reinforcing the ferromagnetic Ising terms that act on those bonds. Since these bonds reinforce the tendency to order in plane, then the stabilization of ferromagnetic order in the lattice plane is enhanced when an in plane field is applied, and the system transitions towards and in plane ferromagnet easily. Furthermore, since at these values of $I_{c}$ and $I_{d}$ the ferromagnetic phase has a polarization vector parallel to the zz-bonds (the $[\bar{1}10]$ direction), then the magnetic field in this direction possesses the smallest critical field.

We will also compare the magnetization in the direction of the field ($OP_{M\cdot H}=\boldmath{M}.\boldmath{H}$) with the ferromagnetic order parameter ($OP_{FM}$) and analyse the susceptibility, $\chi$, of $OP_{FM}$, to determine the critical field at which ferromagnetic order is realized. The ferromagnetic order parameter will probe the system for the presence of a ferromagnetic state, while $OP_{M\cdot H}$ will probe whether the polarization vector is in the direction of the field or not.

For both in plane direction, $[11\bar{2}]$ and $[\bar{1}10]$ (Figs. 12(left) and 12(center)), the behaviour in $OP_{FM}$ and $OP_{M\cdot H}$ is the same. The magnetization monotonically increases until saturation is reached at rather small critical fields (in comparison with the relevant couplings of the model). The susceptibility curves show a peak that indicates a phase transition towards the fully polarized state.

In the case of the field in the $[111]$ direction (Fig 12(right)) the behaviour of $OP_{M\cdot H}$ and $OP_{FM}$ is rather different. Both curves increase monotonically up to $H=0.5$ (unless stated otherwise, the magnetic fields are given in units of $|K|$), at this point they separate into two different behaviours. While $\boldmath{M}.\boldmath{H}$ keeps slowly growing until saturation at approximately $H=0.8$, the ferromagnetic order parameter suddenly reaches saturation at $H=0.6$. This indicates a stable intermediate state, a ferromagnetically ordered phase, in which the ferromagnetic state is realized in a direction not parallel to the applied field. The resulting ferromagnet is a state where the spins are out of plane. The polarization vector of this state can be partitioned into two components. A component parallel to the external field, and one perpendicular to it. In the intermediate state, the perpendicular component is non-zero and aligned along the zz-bonds.

We determine the critical fields by analysing the susceptibility response. We obtain a critical field (in units of $|K|$) $H_{crit}=0.05$ for a field direction $[\bar{1}10]$, $H_{crit}=0.1$ for the magnetic field in the $[11\bar{2}]$ direction, and $H_{crit}=0.8$ for the field in the $[111]$ direction.

IV.2 Anisotropic third neighbour model, $\mathcal{H}_{A}(J_{1,3},K_{1,2},\Gamma_{1,2})$

In this section we will study the magnetization processes for the bond anisotropic version of the model shown in Eq. 3. We will employ the same values of the exchange couplings as those in table 1 and the same field directions as for the nearest neighbour model, $H\parallel[111]$, $H\parallel[\bar{1}10]$, and $H\parallel[11\bar{2}]$. As we shall see, the magnetization processes change drastically from those in the nearest neighbour model, given the strong off-diagonal and further neighbour interactions present in the third neighbour model. In Fig. 13 we show the magnetization in the direction of the field for the three field directions studied. We observe that the polarized state is reached at $H\sim 0.25$ (in units of $|K_{1}|$) in the $[\bar{1}10]$ direction, while for the magnetic field in the $[11\bar{2}]$ direction the critical field is $H\sim 1.5$ and $H\sim 4$ for the field in the $[111]$ direction. This behaviour points towards the existence of an easy axis anisotropy in the $[\bar{1}10]$ direction (please remember that the $[\bar{1}10]$ direction is the direction parallel to the $zz$-bonds). Note that we have dedicated a big part of the numerical effort to fields $H\lesssim 1$, which correspond to experimentally realizable fields.

We show in Fig. 14(top) that the magnetization for the field in the direction $H\parallel[11\bar{2}]$ grows monotonically with a constant slope, reaching saturation at $H\sim 1.5$. At $H=0$ the state is a planar counterrotating spiral as shown in the previous section, but for $H>0$ a continuous transition between a counterotating planar spiral and a counterotating conical spiral is realized. As the name indicates, a conic spiral is a helimagnetic state in which the spiral does not rotate in a plane but in a cone around a certain common direction. Since the magnetic field is applied in the direction parallel to the propagation direction, at moderate fields the spins cant in that direction, which transform the plane of rotation into a cone. We show in Fig. 14(bottom) a scheme of this process. We indicate the direction of propagation (which coincides with the direction of the applied field, $[11\bar{2}]$) as well as the direction perpendicular to it ($[\bar{1}10]$). At low fields the canting in the spins induces the transition from a rotation plane to a rotation cone. As the field increases the cone gets narrower, until at high fields the spins point in the direction of the field, thus reaching saturation.

For the other in plane direction, $H\parallel[\bar{110}]$ (the direction parallel to the $zz$-bonds), the behaviour at small fields is different. We show in Fig. 15(top) the magnetization curve for an applied magnetic field in the $[\bar{1}10]$ direction. We observe that the magnetization monotonically increases with an increasing slope at small fields up to $H\sim 0.3$, at which point it suddenly reaches saturation. This transition towards saturation follows the same mechanism as that developed in section  III.2 for the commensurate-incommensurate transition.

The spiral phase in the $H=0$ case can be seen, projected over the $[111]$ plane, as that of alternating ferromagnetic and antiferromagnetic domains (as in Fig. 7), where the ferromagnetic domains are oriented in the direction parallel to the $zz$-bonds (the $[\bar{1}10]$ direction). If we introduce a small magnetic field in this direction, the wavevector of the magnetic spiral decreases, which is evidenced in the projection on the $[111]$ plane as the increase of the ferromagnetic domains’ size, until they overpower the system at a critical value of the field. In Fig. 15(bottom) we show the change in the spiral wavevector as a function of field. It can be seen that at small fields, $H<0.15$ the change in the wavevector is not pronounced, ranging from $q=0.31$ to $q=0.25$ (in units of $2\pi$). On the other hand, for $H>0.15$, the wavevector decreases rapidly, reaching $q=0$ at $H\sim 0.3$. In this case, a vanishing wavevector indicates the onset of a ferromagnetic state, which is confirmed by the real space pattern and the correlation functions obtained from the simulations. It is also at $H\sim 0.3$ that we observe the maxima in the susceptibility of the ferromagnetic parameter (inset of Fig. 15(top)) indicating that the critical field for this field direction is $H_{\text{crit}}\sim 0.3$.

Since the model contains strong in plane interactions, saturation is only achieved at strong magnetic fields for out of plane directions. In the case of a magnetic field in the $[111]$ direction (direction perpendicular to the lattice plane) the saturation is reached at $H\sim 4$ as can be observed in Fig. 16(left). In this case, as well as in the nearest neighbour model, both order parameters saturate at different field intensities. Both the ferromagnetic order parameter and the magnetization in the direction of the field have to saturate at high enough fields, but while the ferromagnetic order parameter saturates at $H\sim 4$, the magnetization in the direction of the field does not completely saturate up to the biggest calculated fields.

A magnetic field in the $[111]$ direction tilts the spins in this direction, producing a net magnetization, up to the biggest calculated field, which is deviated from the magnetic field direction. The magnetization process can be understood recalling Fig. 7 bottom. Here we mentioned that the spin spiral state can be considered as arrangements of ferromagnetic/antiferromagnetic domains of spins in the $[111]$ plane. When a magnetic field is applied perpendicular to the direction of the ferromagnetic domain’s polarization (as is the case here), the spins of the ferromagnetic domains will be tilted in the direction of the field. The same will happen to the spins in the antiferromagnetic domains which are not already pointing parallel to the field, but a different tilt angle. In Fig. 16(right) we show the spin pattern at $H=0.35$, where the tilt in the direction of the field can be observed.

V Discussion

We have studied a variety of extended Kitaev Hamiltonians as minimal models of $\alpha$-Li2IrO3. The experimental studies of this material indicate a strong bond anisotropy, which led us to study models which contain strong anisotropic bond dependent interactions beyond Kitaev interactions. Furthermore we have considered interactions ranging from nearest to third neighbours. While all models share similarities with each other, we have shown that the range of the interactions can produce radical differences in the behaviour of the system.

Our nearest neighbour models (those given by $\mathcal{H}(J_{1},K_{1},I_{c}\neq 0,I_{d}=0)$ and $\mathcal{H}(J_{1},K_{1},I_{c}\neq 0,I_{d}\neq 0)$ Eq. 2) exhibit a phase diagram where a big part of it presents an incommensurate counterrotating coplanar spiral state. For the model with non zero $I_{c}$ and $I_{d}$ we encounter a spiral state whose wavevector coincides with that of $\alpha$-Li2IrO3 on a section of the phase diagram. However, while in the Iridate material the plane of rotation is tilted away from the lattice plane by $\sim 80^{o}$, in this model the plane of rotation is parallel to the $XY$ Cartesian plane ($54^{o}$ tilt with respect to the lattice plane). This can be understood realizing that a Hamiltonian only containing the $I_{c}$ term is equivalent to an XY model, where in this case the model posses a U(1) symmetry around the Cartesian $z$-axis. The same is true for the $I_{d}$ terms when only one of them is present. When the Kitaev and Heisenberg terms are taken into account together with the $I_{c}$ term, the U(1) symmetry is broken, but there is a residue of this symmetry, which is exhibited in the strong anisotropy of the state, where the rotation plane of the spin spiral state is perpendicular to the Cartesian $z$-axis. The introduction of a small term of the same nature as the $I_{c}$ but over the $xx$- and $yy$-bonds (that is, the $I_{d}$ term in the nearest neighbour model) will induce a tilt of the rotation plane, induced by the symmetries of these terms. The tilt angle will thus be a function of the ratio $I_{c}/I_{d}$, allowing us to reproduce not only the wavevector but also the tilt of the rotation plane characteristic of $\alpha$-Li2IrO3.

While the introduction of second neighbor interactions do not reproduce the experimental results, as shown in Appendices D and E, recent DFT calculations Winter et al. (2016) have put forward a third neighbour model which indicate that further neighbour interactions are needed. We tested this model, in both its bond isotropic and anisotropic forms. Given that the experiments performed on $\alpha$-Li2IrO3 indicate strong bond anisotropy, is no surprise that the model which reproduced the experimental features of the material is the bond anisotropic third neighbour model, $\mathcal{H}_{I}(J_{1,3},K_{1,2},\Gamma_{1,2})$ (Eq. 3).

In the third neighbour models, the interactions are slightly different than in the nearest neighbour models previously studied. We modelled the bond dependent interactions in the nearest and second neighbour models as $I_{c/d}\sum_{<ij>}S^{r_{ij}}_{i}S^{r_{ij}}_{j}$, with $S^{r_{ij}}_{i}=\mathbf{S}_{i}\cdot\mathbf{r}_{ij}$, where $r_{ij}$ is the bond between sites $i$ and $j$. This in turn takes the form, for the $zz$-bonds (the expression for the remaining bonds is analogous)

[TABLE]

which has the form of the $\Gamma$ terms included in the third neighbour model, supplemented by two Heisenberg terms. This expression indeed shows a clear connection between the nearest neighbour model and the nearest neighbour terms of the third neighbour model, where those terms can be deformed into one another by finely tuning the bond anisotropy parameters. While the nearest neighbour terms are connected on both models, the presence of further neighbour exchanges stabilize the properties of the magnetic order. For the nearest neighbour model, the parameters which reproduce the experimental signatures of $\alpha$-Li2IrO3 are constrained to a small part of the phase diagram, close to a phase boundary, which makes the behaviour of the spin spirals highly susceptible to small parameter variation and magnetic fields. On the other hand, the third neighbour model reproduces these signatures over a broader region, making the spiral state robust to parameter variations and applied magnetic fields. Turning to the parameters studied here, it is worth pointing out that we have only mapped the phase diagram as a function of a two parameters per model, the $I_{c}$ and $I_{d}$ terms in the nearest neighbour model, and $J_{1}$ and $J_{3}$ in the third neighbour model, since we have tried to remain as faithful as possible to the exchange couplings proposed in Refs. Winter et al. (2016); Williams et al. (2016) However, while it is expected that the dominant exchange interaction in these materials is of Kitaev type and ferromagnetic, it cannot be discarded that different combination of parameters (in particular regarding the Heisenberg exchange in the nearest neighbour model, or the second neighbour Kitaev and $\Gamma$ exchanges in the third neighbour model) could also reproduce the experimental results.

Since both nearest neighbour and anisotropic third neighbour models reproduce the experimental features of the specific heat, we proceed to calculate the Curie-Weiss temperature within a mean field approach. For the nearest neighbour model we find an anisotropic susceptibility $\mathbf{\chi}=(\chi_{xx},\chi_{yy},\chi_{zz})$ arising from the $I_{c}$ and $I_{d}$ terms, where the associated Curie-Weiss temperatures are given by

[TABLE]

[TABLE]

At the point of interest ($I_{c}=-0.5$ and $I_{d}=-0.35$), where experimental results are reproduced, we obtain positive Curie-Weiss temperatures: $\Theta_{xx}=\Theta_{yy}=16.35K$ and $\Theta_{zz}=14.38K$ in units of the Kitaev coupling, eploying the values of the couplings in meV from Ref. Williams et al. (2016).

In the case of the third neighbour model we also find an anisotropic susceptibility, in this case arising from the anisotropy of the nearest neighbour Heisenberg and Kitaev interactions

[TABLE]

[TABLE]

Evaluated at the experimentally relevant point shown in Table. 1 we obtain $\Theta_{xx}=\Theta_{yy}=20.15K$ and $\Theta_{zz}=-17.17K$.

The experimentally measured Curie-Weiss temperature Singh et al. (2012) was obtained from a fit to the high temperature magnetic susceptibility for polycrystalline samples, $\Theta_{\text{exp}}=-33(3)K$. Since the experiments have been performed on a polycrystalline sample, the Curie-Weiss temperature obtained is the average of the anisotropic Curie-Weiss temperatures. In the case of the nearest neighbour and anisotropic third neighbour models we obtain an averaged Curie-Weiss temperature $\Theta^{I_{c}I_{d}}_{av}=15.7K$ and $\Theta^{W}_{av}=7.58K$. As previously shown, a considerable part of the phase diagram for the anisotropic third neighbour model reproduces the experimental results. This would indicate that a different selection of coupling strengths could modify the averaged Curie-Weiss temperature of the model as to obtain the experimental value, while maintaining the overall behaviour of the system to be the same as in the neutron diffraction experiments. While we could do a similar analysis for the nearest neighbour model, the region of the phase diagram that reproduces the results is much more reduced, and as such the variation of the exchange parameters is not enough to change the sign of the Curie-Weiss temperatures. Extra thermodynamic studies on single crystals would prove useful at determining the anisotropic susceptibilities and the respective signs of the associated Curie-Weiss temperatures. Given that the experiments where performed at high temperature, the comparison between our ground state calculations and the experimental results are not sufficient to determine if the anisotropic third neighbour model is the minimal model of $\alpha$-Li2IrO3.

We studied the magnetization processes for both nearest neighbour and third neighbour models. In both cases we find a strong anisotropy in the magnetic response, with different behaviors for different directions of the applied field. For the nearest neighbour model we find a strong easy plane anisotropy, arising from the off diagonal terms of the Hamiltonian. Since these terms favor in plane orderings, the tendency to order ferromagnetically when an in-plane magnetic field is applied is strong, which presents in our study as a lower critical field for in-plane rather than out of plane directions of the external field. The tendency to in-plane ordering is also observed in the magnetization behavior when a field in the $[111]$ plane is applied, as there exist an intermediate state where the magnetization has a net component in the direction of the field, but also a component in the lattice plane. At intermediate fields, the three field directions studied show similar magnetic behavior, which resembles that of a ferromagnetic/paramagnetic transition, where ferromagnetic domains with polarization parallel to the direction of the external field are induced, and which grow in size as the field is increased. We mention that, since the magnetic fields are expressed in terms of $K$, we can extract the value in Tesla since $H=1$ in units of $|K|$ is equivalent to $g\mu_{B}H/|K|=1$. The gyromagnetic factor is believed to be anisotropic, given the deformations of the octahedral cage of Oxygens, but the exact value is not known. Recent calculations Winter et al. (2017) show that the g factor can drastically change depending on the octahedral deformation, with values ranging from $\sim 1.5$ to $\sim 4$ for the parallel component of the gyromagnetic factor (the component in the lattice plane), and from $\sim 2.5$ to $\sim 0.5$ for the out of plane component. Given the strong variation of this parameter and the lack of experimental data, we have employed an isotropic factor $g=2$. Finally, if we assume the interaction strength of the Kitaev interaction to be $-4.5$meV (as proposed in Ref. Williams et al. (2016)), we obtain the critical fields of $\sim 31$T for the $[111]$ field direction, $\sim 6$T for $[11\bar{2}]$, and $\sim 2.7$T for $[\bar{1}10]$.

For the case of the anisotropic third neighbour model the behaviours is quite different. While the system also exhibits an anisotropic magnetic response, in this case it possesses an easy axis anisotropy. Additionally, assuming that the overall energy scale of the couplings determined by DFT is correct, we find that the critical field for an external field in the $[111]$ direction results $\sim 200$T, while it is $\sim 111$T for the $[11\bar{2}]$ direction, and $\sim 22$T for $[\bar{1}10]$. While such high fields are hard to achieve in experiments, the low field regime already exhibits a highly anisotropic behaviour which can be probed experimentally. For fields in the $[\bar{1}10]$ direction, the magnetization seems to suddenly increase at $|H|\sim 0.25$ (in units of the averaged nearest neighbour Kitaev coupling), exhibiting a change in the spiral wavevector. For fields in the $[11\bar{2}]$ and $[111]$ directions the transition is continuous, presenting a continuous transformation of the spin pattern (in the case of the $[11\bar{2}]$ field direction, from a planar spiral towards a conical spiral).

VI Conclusions

The similarities in the ground states found for both systems indicate that one of them could prove to be the minimal model of $\alpha$-Li2IrO3. The differences in specific heat as well as different magnetization behaviors indicate that further experimental studies are needed to decide which one, if any, is the model corresponding to this material. While the Curie-Weiss temperature for the nearest neighbour model does not reproduce the sign found from the calorimetric experiments, the nature of the experiments (the powder average and the high temperature measures) could obscure further details which could clarify the discrepancy. Studies on single crystals would prove valuable as in these cases the different crystallographic directions could be probed to assert the existence of anisotropic susceptibilities. While Ref. Winter et al. (2016) proposes a set of values for all couplings in the third neighbour model (which were chosen as the starting point for our study), in both isotropic and anisotropic cases these values do not correspond to spiral phases. By changing the couplings of the nearest and third neighbour Heisenberg exchanges we find a range of values where the experimental results are obtained. This indicates that further studies regarding the quantum chemistry as well as the crystal structure of the material are necessary to further determine the strength of the couplings.

To further distinguish the models, magnetization measurements can be fruitful. While the critical fields for the third neighbour model are far beyond current capabilities, interesting studies can be performed at low field. In particular, for the nearest neighbour model the calculated critical fields are within the possibility of experimental realization, while for the third neighbour model, the anisotropy present at low fields can be measured experimentally (we remind the reader that the critical field measured for a field in the $[\bar{1}10]$ direction in third neighbour model is $22$T). In particular, the low field behaviour, which is different for the different field directions, could be observed in single crystals. In this context it is worth pointing out that even if the critical fields obtained experimentally were to results smaller than the ones reported in this work, that would not indicate a failure of the proposed models. It has been shown previously that quantum fluctuations can further reduce the critical fields Winter et al. (2018) for $\alpha$-RuCl3 with respect to those obtained from classical solutions, but that the overall magnetization behaviour is not necessarily modified by these fluctuations.

We conclude then, pointing out that questions remain open regarding the minimal model of $\alpha$-Li2IrO3. We have reduced the number of possible models and explored the different magnetization behavior of those models which reproduce the experimental signatures of the material. We expect that magnetization measurements can point in the direction of one these models being correct in the low temperature limit. However, we recommend further studies: We expect that electronic structure calculations could clarify the current situation in which the obtained exchange couplings do not lead to the experimentally measured spin pattern. On the other hand, and with respect to the third neighbour model, it was pointed out in Ref. Winter et al. (2016) that the crystal structure of the material is not well understood. In their work, the authors of Ref. Winter et al. (2016) also studied a relaxed crystal structure, which coincided with the couplings later obtained in Ref. Freund et al. (2016). We do expect that a study of the kind performed in this work employing the relaxed crystal structure will give a better agreement with experiments for this particular model. However, the importance of further neighbour interactions and anisotropies to the stabilization of the spin spirals will not be modified by this study.

The recent growth of single crystals could help refine the crystal structure. Furthermore, our study of the third neighbour model suggest that modest long-range interactions can stabilize the counterotating spirals found in $\alpha$-Li2IrO3, but that anisotropy is crucial to reproduce the experimental results. Since materials realizing Kitaev interactions show strong bond anisotropy, it cannot be discarded that perhaps a different combination of interactions with a different combination of anisotropies could also reproduce the experimental results and be relevant in real materials. We further mention that also more Monte Carlo studies can be beneficial. By studying the finite temperature magnetization curves, we might encounter lower critical magnetizations and interesting intermediate states which could be easier to realize experimentally (as was shown to be the case for $\beta$-Li2IrO3 Ruiz et al. (2017)). Furthermore, a study of the full phase diagram in the presence of different external magnetic fields could help model future materials which could be realized. Finally, we mentioned that both models could also be distinguished by their behaviour in the presence of dilution. The study performed in Ref. Manni et al. (2014) probe the thermodynamic response of the material when the magnetic ions have been removed, as well as the variation of the spin glass critical temperature as a function of dilution. They conclude that their results indicate that the further neighbour interactions are needed in the minimal model of $\alpha$-Li2IrO3. A detailed study of the models proposed here in the presence of dilution could clarify the picture regarding this material, and it will be left for future studies.

VII Acknowledgements

We would like to thank Roser Valenti, Simon Trebst, Johannes Reuther, and Marek Gluza for fruitful discussions and help with this manuscript. M. L. B is supported by the Freie Universität Berlin within the Excellence Initiative of the German Research Foundation .

Appendix A Benchmark: Heisenberg-Kitaev model

We will benchmark our code by studying the Heisenberg-Kitaev model on the honeycomb lattice. The Hamiltonian for this models is given by,

[TABLE]

where $\gamma=\{x,y,z\}$ and $\mathbf{S}_{i}$ are classical spins on sites $i$ of the honeycomb lattice. This model contains exactly solvable points and has been studied in detail by Price et. al  Price and Perkins (2013) via Monte Carlo simulations employing periodic boundary conditions (PBCs). We reproduce some of the results in Ref. Price and Perkins (2013), and show how (in the case of FEBs), at big enough system sizes we recover the bulk behaviour expected from a simulation employing PBCs.

We perform Monte Carlo simulations on the Hamiltonian shown in Eq. 15, for sizes ranging from 24 to 5400 sites. The temperature of the simulations is consistently chosen as $T=0.001$ in units of $|K|$. We employ $2\times 10^{5}$ Monte Carlo sweeps from which $1\times 10^{5}$ are used as equilibration steps. We concentrate here on the behavior of a simple Metropolis-Hastings algorithm, without recurring to parallel tempering or iterative minimization schemes, to test the effectiveness of the Monte Carlo code. For the results in the main text we employ further modifications of the algorithm to reduce the effects of domain walls and of rough energy landscapes as mentioned in section II.2.

We have mapped the phase diagram for the Heisenberg-Kitaev model. We parametrized the exchange couplings as $K=2\alpha$ and $J=1-\alpha$ and calculated the ground state for values between $\alpha=0$ and $\alpha=1$. We show the phase diagram in Fig. 17.

The phase diagram exhibits two phases, at $\alpha<1/3$ we obtain a ferromagnetic state, while for $\alpha>1/3$ the state is a triple degenerated stripy phase, due to the symmetries of the Kitaev interaction, allowing for the three possible stripy phases, st-X, st-Y, and st-Z.

A.1 Space pattern and correlation functions

For this study we will select the lattice plane as the Cartesian XY-plane, and observe the real space spin pattern on the YZ plane where both Néel and stripy phases are easy to distinguish. We will choose to show phases where the spins are oriented maximally in the $z$-direction to ease the comparison, but states where spins are aligned in other directions are also possible, and have been also obtained within our approach.

We study the real space configurations for different values of $\alpha$ to compare with the low temperature results of  Price and Perkins (2013). We show the results for $\alpha=0$, $\alpha=0.5$, and $\alpha=0.75$ obtained from a calculation employing 216 sites.

For the case $\alpha=0$ we recover the antiferromagnetic Heisenberg model. In Fig 18 we show one of the simulations performed, in which the system exhibits a Neel order, where the spins are ordered along the Cartesian $z$-direction.

For the stripy phases, we show the spin pattern for $\alpha=0.75$ (Fig. 19). In this case we observe that in the center of the system, the spins point in two directions forming stripes that span the system. In this case we separate the system in four sublattices, two with spins pointing in the $+z$ directions (red and pink spins), and two in the $-z$ (blue and green). The finite size effects are noticeable here, where the spins deviate from the $\pm z$ orientation, and where this deviation is more pronounced at the boundaries of the system. Even for this small system size, where the finite size effects are considerable, the simulation already reproduces the known results for the Heisenberg-Kitaev model for these parameters.

For the stripy phase, an important point needs to be mentioned. The point $\alpha=0.5$ is special, since at this point the obtained stripy phase becomes and exact ground state of the system, which can be seen as a ferromagnetic state in a rotated basis. It can be proven that in this case the system exhibits a SU(2) symmetry Price and Perkins (2013). In Fig.20 we show a spin pattern for this case. The ground state realizes a stripy phase, but in this particular case the spins are not arranged along only one direction, due to the emergent symmetry of the state. While in the $\alpha=0.5$ case some disorder in the state can be seen, these finite size effects are not as pronounced as in the rest of the stripy state.

These findings coincide with what was reported in Ref. Price and Perkins (2013) for the full phase diagram. We have purposefully showed results for a relatively small system size (216 sites), to exemplify the effect of FEBs.

We proceed to show in Fig. 21 the Fourier transform of the correlation function for both the antiferromagnetic and stripy phases, obtained from simulations on 2400 sites.

Fig. 21(left) shows the Fourier transform for the antiferromagnetic state. We observe the distinct features of the Néel phase, presenting maxima at the corners of the extended Brillouin zone. On the other hand, the stripy phase (Fig. 21, right) presents two maxima as well as secondary maxima on the sides of the extended Brillouin zone. This maxima indicate that, in this particular simulation, the state is dominated by a st-Y phase. The secondary maxima on the sides of the extended Brillouin are effects arising from the FEBs, where small regions of the system realize the other two degenerate stripy phase.

A.2 FEBs effect on big system sizes

Even though the effect of FEBs is still noticeable in system sizes as big as 2400 sites, the behavior of the system is not radically affected by them. While the presence of FEBs has the tendency to generate domain walls when the low temperature state might be degenerate, the simulations consistently converge to states that correspond the expected behavior of the system.

In the following, we show the results of two identical simulations, run with different seeds, for the case $\alpha=0.75$ and 5400 sites. In this case, both simulations converge to the stripy phase on the $z$-direction. Fig. 22 shows the stripy phase for a simulation in which no domain walls have been generated. In this case, away from the bulk the spins deviate from their $\pm z$ orientation, but it is only close to the edges of the system where the effects of FEBs destroy the stripy order. Furthermore, these deviations do not take an arbitrary form. When the spins deviate from the $\pm z$ direction, they do it by inducing a non-zero $x$- and/or $y$-spin components arranged according to the st-X and st-Y phases.

On the other hand, Fig. 23 shows a spin pattern where a domain wall has been generated in the middle of the system. In this case, both sides of the domain wall present a stripy order, but the orientation of the sublattices is flipped. In this case, the domain wall also generate non-zero x- y-spin components which are arranged accordingly to the st-X and st-Y phases.

While in Fig 23 the effects of FEBs seem to be more noticeable than in Fig.22, the correlations for both states are the same, i.e, the Fourier transform of the correlation function presents the same characteristics as Fig. 21.

In this work we implement further algorithms to minimize the effect of FEBs, but while their effect can be minimized, they cannot be eliminated entirely. With this in mind, and the comparison of our results to the work in Ref. Price and Perkins (2013), we are ensured that while FEBs will induce domains (when the ground state is degenerated) and edge defects in the systems, these effects will not modify the nature of the studied system nor their correlations. If the simulations is able to converge to a low temperature state (please note that here we are not considering other effects such as local minima or critical slowing down) then we are confident that this state will be a low temperature state of the studied Hamiltonian.

Appendix B $\mathcal{H}(J_{1},K_{1},I_{c}\neq 0,I_{d}=0)$ model

Kimchi et. al. proposed Kimchi et al. (2015) a nearest neighbour Hamiltonian that is capable of capturing the counterotating features of the ground state properties common to these three compounds. The Hamiltonian takes the form of Eq. 2 with $I_{d}=0$, explicitly

[TABLE]

where $<ij>_{\gamma}$ indicates the bonds in which the interactions acts. This Hamiltonian contains terms that are symmetry allowed by the microscopic structure of the materials, and in their paperKimchi et al. (2015) Kimchi and collaborators studied it via a reduction to a 1D chain model and a subsequent solution employing an LT approximation. The 1D toy model relies on the reduction of the Hamiltonian in Eq. 16 to that of decoupled zig-zag chains by taking the $I_{c}$ term to zero and introducing a second neighbour Heisenberg interactions. This Hamiltonian is solved by proposing an ansatz for the spin components and then the interchain coupling arising from the $I_{c}$ terms is slowly introduced. They subsequently solve the perturbed Hamiltonian by an LT approximation. They find that even when the second neighbour interactions is zero and the $I_{c}$ term is completely introduced, the spin spirals survive in the phase diagram. We study the model in Eq. 16 exactly, via Monte Carlo simulations to test the accuracy of our methodology at detecting incommensurate states.

Phase diagram

Employing large scale Monte Carlo simulations as detailed in section II ,we are able to map the phase diagram of this model ($\mathcal{H}(J_{1},K_{1},I_{c}\neq 0,I_{d}=0)$) in the phase space of the $J$ and $I_{c}$ couplings, where the Kitaev exchange is consistently set to $-1$. The obtained phase diagram can be seen in Fig. 24. This diagram shows evidence of the strong Kitaev interactions present, displaying a broad region where the stripy phase lives. At the boundary $I_{c}=0$ ($J\neq 0$) we recover the Kitaev-Heisenberg model, which presents a stripy phase for these values of $J$ and $K$. Furthermore, and as expected, this state is affected by a spin locking effect. As soon as we set $I_{c}\neq 0$ the degeneracy of the stripy phase is broken, and a stripy in the $z$ direction is chosen (we will name this phase st-Z), which we show in the phase diagram by blue squares.

On the other limit, when $J=0$ and $I_{c}\neq 0$ the model in Eq. 16 reduces to two ferromagnetic couplings, which produces a large ferromagnetic phase which survives up to finite $J$ (black dots in Fig. 24). The case $I_{c}=0$ and $J=0$ is special, as this point reduces to the Kitaev model, which is a macroscopically degenerate state without LRO. For finite $J$ and $I_{c}$ we find an incommensurate spiral phase (green squares in Fig. 24). It is worth noting that Kimchi et.al report the existence of a small area where the st-X and st-Y phases should be present, and they locate it at the intersection of the spiral and st-Z phase.

The incommensurate order present in the phase diagram (green squares in Fig. 24) exhibits signatures of incommensurate counterotating spirals. This spiral phase propagates in the horizontal direction according to Fig. 1 (the direction perpendicular to the zz-bonds), with a wavevector that varies from $0.5$ to $0.30$ in units of $2\pi$. The phase of rotation tilted with respect to the lattice plane remains constant throughout the phase diagram at $\theta=54^{o}$, i.e, the rotation plane is oriented parallel to the XY-Cartesian plane. This phase mostly reproduces the experimental structure determined from neutron diffraction except for the tilt of the rotation plane (we find a tilt of $54^{o}$ and the expected value for the tilt angle is $80^{o}$).

Appendix C Isotropic third neighbour model, $\mathcal{H}_{I}(J_{1,3},K_{1,2},\Gamma_{1,2})$

Phase diagram

We will start by mapping part of the phase diagram of the isotropic case for the model shown in Eq. 3. We will implement a dominant nearest neighbour Kitaev coupling $K_{1}=-1$ supported by (here and in the following, all exchange couplings are given in units of $|K_{1}|$) $K_{2}=-0.275$, $\Gamma_{1}=1$, $\Gamma_{2}=0.275$, $J_{1}\in(-0.1,...,0.4)$, and $J_{3}\in(0,...,0.4)$.

The phase diagram is shown in Fig. 25. In this model, the boundary between the different phases was not clearly distinguished from the Monte Carlo simulations, and future work will be needed to obtain reliable bounds. However, as we are interested in the nature of the phases arising in these models, we leave the study of the precise location of the boundaries for future work. In the phase diagram shown in Fig. 25 the phase boundaries will be located over the white spaces.

We observe three commensurate phases, and an incommensurate one. In the limit $J_{3}=0$ a transition from a ferromagnetic to a stripy phase at $J_{1}\sim 0.15$ is realized. At a critical value of the third neighbour coupling, $J_{3}\sim 0.5$, a commensurate-incommensurate transition takes place and an intermediate incommensurate phase emerges. At a critical value of $J_{1}$ (which depends on the value of $J_{3}\gtrsim 0.05$) the system enters into a spiral phase (green squares in Fig. 25). If $J_{1}$ is further increased the system can enter a stripy phase ($J_{3}\lesssim 0.1$), remains in the spiral phase, or enter an antiferromagnetic state ($J_{3}\gtrsim 0.25$). Both antiferromagnetic and stripy phases have spins oriented perpendicular to the lattice plane, which allows us to classify the stripy phase as st-Z order. Both phases are two fold degenerated, and this arises in the spin pattern as domains separating these degenerate states. We have confirmed the existence of these phases via a study of the Fourier transform of the correlation functions obtained from Monte Carlo, as well as with a LT minimization. These observations confirm the results of Ref. Winter et al. (2016).

Spiral properties

The incommensurate order present in this model exhibits the signatures of a co-planar counterrotating spin spiral. The fundamental difference between these spirals and the ones appearing in $\alpha$-Li2IrO3 is that in this case, the spirals are degenerate. Since the model is bond isotropic, the spin spirals can propagate in three symmetry allowed directions, i.e. the propagation direction can be perpendicular to the $x$-, $y$- or $z$ bonds. In our simulations this implies that our state can contain domain walls separating domains where the spin spirals propagate in different directions. This can be seen in Fig. 26 the degeneracy is evident in the Fourier transform of the correlation function. Here we see that indeed, we obtain maxima inside the first Brillouin zone, and that these are six fold degenerated. The study of the wavevector and nature of the spirals confirm that these spirals are of the same type as those obtained in the nearest neighbour model shown in the main text, with the caveat that the state in these models is degenerated.

Appendix D $\mathcal{H}(J_{1,2},K_{1,2},I_{c}=0,I_{d}=0)$ model

We proceed now to study the first of the second neighbour models, which consists of first and second Heisenberg and Kitaev interactions. We show in Fig. 1 the exchange couplings in the honeycomb lattice for this model. The Hamiltonian corresponding to this model is given by Eq. 17.

[TABLE]

This model was previously studied in the quantum limit, and in the context of the $\mathrm{Li}_{2}\mathrm{IrO}_{3}$ family, via the pseudofermionic functional renormalization group method (PFFRG) Reuther et al. (2014). In their paper, Reuther et. al parametrize the different couplings via two angles, $P_{1}$ and $P_{2}$, as $J_{1}=\mathrm{cos}(\pi P_{1}/2)$, $K_{1}=-\mathrm{sin}(\pi P_{1}/2)$, $J_{2}=-g\mathrm{cos}(\pi P_{2}/2)$, $K_{2}=g\mathrm{sin}(\pi P_{2}/2)$, and map the phase diagram of 17 for $P_{1}\in(0,1)$ and $P_{2}\in(0,1)$.

They find two incommensurate spiral phases for $P_{2}\gtrsim 0.5$, the spirals SP1 and SP2, while for values below $0.5$ they find a ferro and antiferromagnet. Studying the maximum of the susceptibility they find that the state SP1 corresponds to maxima outside of the first Brillouin zone, while the SP2 has maxima inside the first zone. We will employ Monte Carlo simulations to study the classical equivalent of this model, and we will restrict ourselves to the value $g=0.8$, as according to the evidence in Ref.Reuther et al. (2014) the phase diagram does not change drastically for different values of $g$.

Phase diagram

The phase diagram for the this model, as obtained from our Monte Carlo simulations, is shown in Fig. 27. Please note that since the aim of this study is to identify the nature of the spiral phases, we have not gone to great lengths mapping the boundary in between the phases in this model. In the phase diagram shown in Fig. 27 the phase boundaries will be located somewhere over the white spaces.

The phase diagram presents two commensurate and two incommensurate phases. Pink dots represent antiferromagnetic order, black ones indicate the onset of ferromagnetic order. We indicate by red diamonds the incommensurate spiral phase SP1, and by green squares the phase SP2. We will study these phases in detail in the next section.

The ferro and antiferromagnetic phases exhibit clear features in the correlation function which confirms their nature. However, when we study the real space pattern of the spins we find that these phases come together with domain walls as well as with vortex-like defects.

Spiral properties

The two incommensurate spin spiral phases we find correspond to the red diamonds in Fig. 27 (SP1), and to the green squares (SP2). Upon examination of the correlation function (Fig. 28) we see that the SP1 phase present maxima outside of the first Brillouin zone, while the maxima corresponding to the SP2 phase are contained within the first Brillouin zone.

The Fourier transform of the correlation function for the SP2 phase presents maxima as satellites around the $\Gamma$ point, but no secondary maxima around the $K$ points are observed, unlike the expected signal for $\alpha$-Li2IrO3. Furthermore a degeneracy is present, in which the maxima appear in the three symmetry related positions, which would indicate that the spin spirals propagate in the three directions allowed by the Kitaev symmetry. On the other hand, the spiral phase SP1 presents only maxima as satellite peaks around the $K$ points.

Upon a closer inspection of the real space spin pattern we see that for both the SP1 and SP2 phases, the actual spiral present in this model does not coincide with that of $\alpha-\mathrm{Li}_{2}\mathrm{IrO}_{3}$. In the experimental case one observes co-planar spirals, where the plane of rotation can be defined by $(\mathbf{S}_{1}\times\mathbf{S}_{2})\cdot\mathbf{S}_{3}$, where $S_{x}$ (with $x=\{1,2,3\}$) are consecutive spins on a single spiral. However, this model exhibits a non-coplanar helimagnet, where a plane of rotation cannot be defined.

While both spirals are non-coplanar, we can further distinguish them by studying their nearest neighbor correlations. Calculating the correlation function up to nearest neighbours, and Fourier transforming it, we obtain Fig. 30. From here we see that the main distinguishing feature between both non-coplanar spirals resides in the nature of the nearest neighbor correlations. While the maxima are broadened, we see that for phase SP1 we find maxima in the corners of the extended Brillouin, which indicates antiferromagnetic nearest neighbor correlations. On the other hand, SP2 has a maximum in the $\Gamma$ point, which coincides with ferromagnetic correlations.

Appendix E $\mathcal{H}(J_{1,2},K_{1,2},I_{c}\neq 0,I_{d}=0)$ model

Since $\alpha-\mathrm{Li}_{2}\mathrm{IrO}_{3}$ presents bond dependent interactions, a question that arises after studying the $2N$-model is what the effect of other symmetry allowed interactions are in this particular model. It is not absurd to think that an $I_{c}$-like term could break the degeneracy and perhaps induce a coplanar spiral. For this purpose we will now study a modification of the $2N$ model, where we introduce an $I_{c}$ term of the same form as used in the nearest neighbour models. The Hamiltonian results.

[TABLE]

As before we parametrize the different couplings via two angles, $P_{1}$ and $P_{2}$, as $J_{1}=\mathrm{cos}(\pi P_{1}/2)$, $K_{1}=-\mathrm{sin}(\pi P_{1}/2)$, $J_{2}=-g\mathrm{cos}(\pi P_{2}/2)$, $K_{2}=g\mathrm{sin}(\pi P_{2}/2)$, and set $g=0.8$. The coupling $I_{c}$ will take values in the range $\{-1,-0.9,...,0\}$. We run Monte Carlo simulations for the pairs $\{P_{1}=0.5,P_{2}=0.5\}$, $\{P_{1}=0.5,P_{2}=1\}$, $\{P_{1}=1,P_{2}=0.5\}$, $\{P_{1}=1,P_{2}=1\}$, and $\{P_{1}=0.9,P_{2}=0.9\}$ as a way to probe the effect of the $I_{c}$ term at different points in the phase diagram.

$p_{1}=0.5$, $p_{2}=0.5$

For the limit $I_{c}=0$ an antiferromagnetic state is realized, as shown in Fig. 27. This state is degenerate, with spins able to point in different symmetry allowed directions. For $I_{c}>-0.5$ the state remains an antiferromagnet, but now the degeneracy disappears and the ground state selects the state with spins pointing in plane and perpendicular to the $zz$-bonds, as shown in Fig. 31. For $I_{c}<-0.5$ the state changes slightly, maintaining its antiferromagnetic nature, but with the staggered magnetization in a direction perpendicular to the lattice plane.

$p_{1}=1$, $p_{2}=0.5$

For these set of parameters, a ferromagnetic state is realized, where vortex-like defects appear. When the value of $I_{c}$ is non zero, the vortex defects disappear, and we find a ferromagnetic state with a net magnetization in the direction of the zz-bonds. This configuration is reached for the smallest values of $I_{c}$ studied and remains unchanged through the whole range. Here, two possible orientations of the net magnetization are possible, and they appear through the simulation separated by extended domain walls spanning through the system.

$p_{1}=0.5$, $p_{2}=1$

The case for $I_{c}=0$ realized a helimagnetic state for these parameter values.The inclusion of a finite $I_{c}>-0.2$ breaks the degeneracy of the spin spiral states, and spirals only propagate in the direction perpendicular to the $zz$-bonds. The nature of the spin spiral also changes with respect to the case $I_{c}=0$, the spirals are still non-coplanar, but the nearest neighbor correlations are not purely antiferromagnetic, as the $I_{c}$ term introduces a ferromagnetic binding in the $zz$- bonds. For $I_{c}=-0.3$ a stripy phase dominates the phase diagram, with domain walls separating different orientations of the spins. In this stripy phase, the spins are in plane, aligned in the direction of the zz-bonds. We show in Fig. 32 the resulting stripy phase for the case $I_{c}=1$.

$p_{1}=1$, $p_{2}=1$, and $p_{1}=0.9$, $p_{2}=0.9$

For the limit $I_{c}=0$ we observed a spin spiral state which we analysed in the previous section. When $I_{c}>0$ the non-coplanar nature of the spin spiral disappears, the anisotropy along the $zz$-bonds introduced by the ferromagnetic $I_{c}$ term enforces a ferromagnetic in plane alignment of the spins in this bond, which transforms the spin spiral into an antiferromagnetic state. As $I_{c}$ decreases the state remains unchanged.

Since one would expect that the introduction of an $I_{d}$ term would modify the plane of rotation, but not the wavevector, as we have shown in the nearest neighbor model that only an $I_{c}$ term is enough to generate incommensurate counterotating spirals, a further modification of the model shown in this appendix was not attempted. We can conclude then, that the second neighbour models are not minimal models for $\alpha$-Li2IrO3.