Magnetic Field Induced Competing Phases in Spin-Orbital Entangled Kitaev Magnets

TL;DR

This paper investigates the complex magnetic phases in Kitaev magnets under magnetic fields, revealing large unit cell orders and quantum fluctuations that influence thermal Hall effects, advancing understanding of spin liquids.

Contribution

It provides a semiclassical analysis of large system sizes, uncovering competing magnetic orders missed by previous quantum studies, and links these to thermal Hall conductivity.

Findings

Large unit cell magnetic orders at intermediate fields

Strong quantum fluctuations in these orders

Potential for spin liquid states and thermal Hall effects

Abstract

There has been a great interest in magnetic field induced quantum spin liquids in Kitaev magnets after the discovery of neutron scattering continuum and half quantized thermal Hall conductivity in the material $\alpha$-RuCl$_3$. In this work, we provide a semiclassical analysis of the relevant theoretical models on large system sizes, and compare the results to previous studies on quantum models with small system sizes. We find a series of competing magnetic orders with fairly large unit cells at intermediate magnetic fields, which are most likely missed by previous approaches. We show that quantum fluctuations are typically strong in these large unit cell orders, while their magnetic excitations may resemble a scattering continuum and give rise to a large thermal Hall conductivity. Our work provides an important basis for a thorough investigation of emergent spin liquids and competing phases in Kitaev magnets.

Equations

$$H=\sum _{\lambda =x,y,z}\sum _{⟨ij⟩\in \lambda }[K{S}_i^{\lambda }{S}_j^{\lambda }+\Gamma (S_i^{\mu }{S}_j^{\nu }+S_i^{\nu }{S}_j^{\mu })$$

$$+{\Gamma }^{\prime }(S_i^{\mu }{S}_j^{\lambda }+S_i^{\lambda }{S}_j^{\mu }+S_i^{\nu }{S}_j^{\lambda }+S_i^{\lambda }{S}_j^{\nu })]-\mathbf{h}\cdot \sum _i{\mathbf{S}}_i,$$

$$\frac{⟨S_{\mathrm{LSWT}}⟩}{S}=1-\frac{{\sum }_i⟨b_i^{†}{b}_i⟩/N_{\mathrm{site}}}{S},$$

$${\kappa }_{xy}=-\frac{k_{\mathrm{B}}^{2}T}{\hslash V}\sum _{\mathbf{k}}\sum _{n=1}^{\mathcal{N}}\left\{c_2\left[f_{\mathrm{BE}}({\epsilon }_{\mathbf{k},n})\right]-\frac{{\pi }^2}{3}\right\}{\Omega }_{\mathbf{k},n}.$$

$$H=N{\mathbf{S}}^{\mathrm{T}}\left(H_K+H_{\Gamma }\right)\mathbf{S}-2N\mathbf{h}\cdot \mathbf{S},$$

$$H_K=\left(\begin{array}{ccc}K& 0& 0\\ 0& K& 0\\ 0& 0& K\end{array}\right),H_{\Gamma }=\left(\begin{array}{ccc}0& \Gamma & \Gamma \\ \Gamma & 0& \Gamma \\ \Gamma & \Gamma & 0\end{array}\right),$$

$$\begin{gathered} \tilde{\mathbf{S}}=\left(\begin{array}{c}{S}^a\\ S^b\\ S^c\end{array}\right) \\ =\left(\begin{array}{ccc}\frac{1}{\sqrt{6}}& \frac{1}{\sqrt{6}}& -\sqrt{\frac{2}{3}}\\ -\frac{1}{\sqrt{2}}& \frac{1}{\sqrt{2}}& 0\\ \frac{1}{\sqrt{3}}& \frac{1}{\sqrt{3}}& \frac{1}{\sqrt{3}}\end{array}\right)\left(\begin{array}{c}{S}^x\\ S^y\\ S^z\end{array}\right) \\ =R\mathbf{S}. \end{gathered}$$

$${\tilde{H}}_{\Gamma }=R{H}_{\Gamma }{R}^{\mathrm{T}}=\left(\begin{array}{ccc}-\Gamma & 0& 0\\ 0& -\Gamma & 0\\ 0& 0& 2\Gamma \end{array}\right).$$

$${\mathbf{B}}_i=-\sum _j{H}_{ij}{\mathbf{S}}_j+\mathbf{h},$$

$$R_i=\left(\begin{array}{ccc}\cos {\theta }_i\cos {\varphi }_i& -\sin {\varphi }_i& \sin {\theta }_i\cos {\varphi }_i\\ \cos {\theta }_i\sin {\varphi }_i& \cos {\varphi }_i& \sin {\theta }_i\sin {\varphi }_i\\ -\sin {\theta }_i& 0& \cos {\theta }_i\end{array}\right),$$

$${\tilde{S}}_i^z=S-b_i^{†}{b}_i=S-n_i,$$

$${\tilde{S}}_i^x=\frac{\sqrt{2S-n_i}{b}_i+b_i^{†}\sqrt{2S-n_i}}{2}\approx \sqrt{\frac{S}{2}}(b_i+b_i^{†}),$$

$${\tilde{S}}_i^y=\frac{\sqrt{2S-n_i}{b}_i-b_i^{†}\sqrt{2S-n_i}}{2i}\approx -i\sqrt{\frac{S}{2}}(b_i-b_i^{†}),$$

$$H=\sum _{ij}{\mathbf{S}}_i^{\mathrm{T}}{H}_{ij}{\mathbf{S}}_j-{\mathbf{h}}^{\mathrm{T}}\sum _i{\mathbf{S}}_i=\sum _{ij}{\tilde{\mathbf{S}}}_i^{\mathrm{T}}{\tilde{H}}_{ij}{\tilde{\mathbf{S}}}_j-\sum _i{\tilde{\mathbf{h}}}_i^{\mathrm{T}}{\tilde{\mathbf{S}}}_i,$$

$$b_{\mathbf{k},s}=\frac{1}{\sqrt{N}}\sum _i{b}_{i,s}{e}^{\mathbf{k}\cdot {\mathbf{R}}_i},$$

$$H=\sum _{\mathbf{k}}{\Psi }_{\mathbf{k}}^{†}{\mathrm{D}}_{\mathbf{k}}{\Psi }_{\mathbf{k}},$$

$${\mathrm{D}}_{\mathbf{k}}=\left(\begin{array}{cc}{\mathrm{A}}_{\mathbf{k}}& {\mathrm{B}}_{\mathbf{k}}\\ {\mathrm{B}}_{-\mathbf{k}}^{\ast }& {\mathrm{A}}_{-\mathbf{k}}^{\mathrm{T}}\end{array}\right),$$

$${\mathrm{T}}_{\mathbf{k}}^{†}{\mathrm{D}}_{\mathbf{k}}{\mathrm{T}}_{\mathbf{k}}={\mathcal{E}}_{\mathbf{k}},{\mathrm{T}}_{\mathbf{k}}{\sigma }^3{\mathrm{T}}_{\mathbf{k}}^{†}={\sigma }^3,$$

$$\frac{1}{N_{\mathrm{site}}}\sum _{is}⟨b_{is}^{†}{b}_{is}⟩=\frac{1}{N\mathcal{N}}\sum _{\mathbf{k}}\sum _{m,n=1}^{\mathcal{N}}\{{\mathrm{T}}_{\mathbf{k}}^{\ast }(m,n){\mathrm{T}}_{\mathbf{k}}(m,n)f_{\mathrm{BE}}({\epsilon }_{\mathbf{k},n})$$

$$+{\mathrm{T}}_{\mathbf{k}}^{\ast }(m,n+\mathcal{N}){\mathrm{T}}_{\mathbf{k}}(m,n+\mathcal{N})\left[1+f_{\mathrm{BE}}({\epsilon }_{-\mathbf{k},n})\right]\},$$

$$f_{\mathrm{BE}}(\epsilon )=\frac{1}{e^{\epsilon /T}-1}.$$

$$c_2(x)$$

$$=(1+x){\left(\ln \frac{1+x}{x}\right)}^2-(\ln x{)}^2-2{\mathrm{Li}}_2(-x),$$

$${\Omega }_{\mathbf{k},n}=i{ϵ}_{\mu \nu }{\left[{\sigma }^3\frac{\partial {\mathrm{T}}_{\mathbf{k}}^{†}}{\partial k_{\mu }}{\sigma }^3\frac{\partial {\mathrm{T}}_{\mathbf{k}}}{\partial k_{\nu }}\right]}_{nn},$$

$$C_n=\frac{1}{2\pi }\sum _{\mathbf{k}}\frac{(2\pi {)}^2}{A}{\Omega }_{\mathbf{k},n},$$

$$S_{\mathbf{k}}=\frac{1}{N_{\mathrm{site}}}\sum _{ij}{\mathbf{S}}_i\cdot {\mathbf{S}}_j{e}^{i\mathbf{k}\cdot ({\mathbf{r}}_i-{\mathbf{r}}_j)},$$

$${\mathbf{S}}_i\stackrel{C_3}{⟶}{C}_3^{-1}{\mathbf{S}}_{C_3(i)}.$$

Full text

Discovery of quantum spin liquidsZhou et al. (2017); Takagi et al. with emergent quasiparticles has been an important subject in modern condensed matter physics. This serves as an ultimate test of our understanding of highly quantum entangled phases in interacting electron systems. Recent research has invested tremendous effort on a number of materials with strong spin-orbit couplingWitczak-Krempa et al. (2014); Rau et al. (2016), which leads to intriguing bond dependent exchange interactions between spin-orbital entangled pseudospin-$1/2$ moments. These studies are largely motivated by the exact solution of the Kitaev honeycomb modelKitaev (2006). The Kitaev interaction is naturally present in the systems with $4d/5d$ transition metal elementsJackeli and Khaliullin (2009), such as honeycomb/hyperhoneycomb iridatesChaloupka et al. (2010); Katukuri et al. (2014); Takayama et al. (2015) and $\alpha$-RuCl3Plumb et al. (2014). However, other exchange interactions are present tooRau et al. (2014), which often lead to magnetically ordered ground states instead of the desired quantum spin liquidLiu et al. (2011); Ye et al. (2012); Sears et al. (2015); Johnson et al. (2015). Hence much effort has been spent to suppress the magnetic orders and gain access to the possible spin liquid phases.

Over the past few years, great experimental progress has been achieved in $\alpha$-RuCl3. At zero magnetic field, this material orders magnetically in the zigzag (ZZ) orderSears et al. (2015); Johnson et al. (2015). Upon the application of an external field, neutron scattering experimentsWinter et al. (2018); Banerjee et al. (2018); Balz et al. find an intermediate window of fields before the system enters the polarized state, where sharp magnon modes are absent, but a scattering continuum appears instead. Under a $[111]$ field (perpendicular to the honeycomb plane), the measured thermal Hall conductivity above the ordering temperature $T_{\mathrm{N}}\approx 7\,\mathrm{K}$ follows the predicted trend of itinerant Majorana fermions in the pure Kitaev modelKasahara et al. (2018a). When the field is tilted away from the $[111]$ direction by $45$°  and $60$°, half quantized thermal Hall conductivity is observedKasahara et al. (2018b). These observations raise the hope that the paramagnetic state in the intermediate field regime may be the sought-after chiral spin liquid with Majorana edge modes.

Theoretical models for $\alpha$-RuCl3 include substantial Kitaev and symmetric anisotropic $\Gamma$ interactions, both strongly dependent on the bond directions, with additional small exchanges such as the neareast neighbor Heisenberg $J$, the third nearest-neighbor Heisenberg $J_{3}$, and the anisotropic $\Gamma^{\prime}$Kim and Kee (2016); Winter et al. (2016); Wang et al. (2017), on the honeycomb lattice. Previous analyses are largely done on quantum models with small system sizes (typically a $24$-site cluster) via exact diagonalization (ED) Catuneanu et al. (2018); Rusnačko et al. (2019); Jiang et al. ; Gordon et al. ; Kaib et al. or in quasi-one dimensional limit via density matrix renormalization group (DMRG)Gohlke et al. (2018); Jiang et al. ; Gordon et al. , with varying degree of complexity. For example, a recent workGordon et al. on the $K\Gamma\Gamma^{\prime}$ model in an external magnetic field suggests that it allows an intermediate spin liquid phase continuously connected to the pure Kitaev model between the low field ZZ order and high field polarized state.

In this article, we investigate the possible competing phases in the classical $K\Gamma\Gamma^{\prime}$ model under a $[111]$ magnetic field for large system sizes. The purpose is to critically examine what kind of competing phases may be present and how these phases may be related to potential spin liquids in the quantum model. Rather surprisingly, we find a series of competing magnetic orders with large unit cells in the intermediate field regime. In particular, in the $K\Gamma\Gamma^{\prime}$ model with small $\Gamma^{\prime}$, the ground state in the zero field limit is the ZZ order, which is consistent with previous experiments and theoretical calculations. Upon increasing the field, the ZZ order is replaced by a series of magnetically ordered phases with $8$, $18$, $32$, $50$, $70$ and $98$-site unit cells before the system enters the polarized state (see Fig. 2). Hence the magnetic field reveals a series of competing orders, which form an intermediate region in the phase diagram. Most of these large unit cell orders had not been identified in previous works.

We compute the zero point quantum fluctuations for these magnetic orders and estimate the reduction of the size of the local moments. We find that quantum fluctuations are strong in the large unit cell orders so that the renormalized local moment is only about $50\%$ of the full magnitude on average. The flat and dense spin wave spectra in the large unit cell orders, in particular the $70$ and $98$-site orders, essentially look like continua of spin excitations. Furthermore, we calculate the thermal Hall conductivity due to magnons in some of the large unit cell orders and find that it is as large as that observed experimentally at low temperatures. While strong quantum fluctuations are present and hence it is likely for the series of competing phases to turn into spin liquids in the quantum limit, it is also evident that previous theoretical studies on quantum models with small system sizesGohlke et al. (2018); Catuneanu et al. (2018); Rusnačko et al. (2019); Jiang et al. ; Gordon et al. ; Kaib et al. cannot resolve many of these large unit cell orders. Therefore, in future analyses of such quantum models, it will be important to understand the role of quantum fluctuations in the large unit cell orders unveiled in the current work. Our findings demonstrate the possibility of novel and exotic ordering patterns in spin-orbital entangled Kitaev magnets, which provide an important basis for further investigations of the origin and the nature of quantum spin liquids that they may host.

Results

Model. We investigate the nearest neighbor $K\Gamma\Gamma^{\prime}$ model on the honeycomb lattice in a $[111]$ magnetic field $\mathbf{h}$,

[TABLE]

where $K$ is the Kitaev interaction, $\Gamma$ and $\Gamma^{\prime}$ are off diagonal spin exchanges, $(\lambda,\mu,\nu)$ is a cyclic permutation of $(x,y,z)$ and the field $\mathbf{h}=h(1,1,1)/\sqrt{3}$. We have also assumed an isotropic $g$ tensor. In (1), $h$ actually carries a factor of $S$ but for notational simplicity we will just write $h$ in units of the Kitaev interaction, for instance $h=0.1\lvert K\rvert$ instead of $h=0.1\lvert K\rvert S$, in the rest of this article.

In the experimentally relevant parameter regime, $K<0$ and $\Gamma>0$ are large while $\Gamma^{\prime}<0$ is small. In contrast to many of the previous studiesGohlke et al. (2018); Catuneanu et al. (2018); Rusnačko et al. (2019); Gordon et al. ; Jiang et al. ; Wang et al. ; Kaib et al. , we investigate the classical limit of this model, that is, by treating the spins $\mathbf{S}_{i}=(S_{i}^{x},S_{i}^{y},S_{i}^{z})$ in (1) as three dimensional vectors of fixed magnitude $\lvert\mathbf{S}_{i}\rvert=S$ for all $i$. We use simulated annealing to determine the ground state spin configuration of the system. Details of the simulated annealing calculation can be found in the Methods section.

Phase diagrams. We first consider the $K\Gamma$ model by setting $\Gamma^{\prime}=0$ in (1), with a ferromagnetic Kitaev interaction $K=-1$. We explore $\Gamma\in[0,0.5],h\in[0,1.2]$ and map out the phase diagram, as shown in Fig. 1. Apart from the extensively degenerate Kitaev limit $(\Gamma,h)=(0,0)$, we find that the vast majority of the parameter space favors particular magnetic orders. All these ordered phases, except the zigzag (ZZ) order and the ferromagnet (FM), are labeled by the number of sites contained in their respective magnetic unit cells.

In the zero field limit, the degeneracy of the Kitaev manifold is lifted as the ZZ order and a $12$-site order (a $6$-site order and an $18$-site order) are selected at small (intermediate) $\Gamma$. These two phases have exactly the same energy at $h=0$, but the ZZ order or the $6$-site order is prefered once $h\neq 0$. However, the 18-site order reemerges at higher fields and replaces the 6-site order as the ground state. Tracing back to the parameter region with small $\Gamma$ and $h$, we see that the ZZ, $6$-site, $12$-site and $18$-site orders are continuously connected to the Kitaev limit. The $6$-site order (the $18$-site order) was first reported in Ref. Lampen-Kelley et al., (Ref. Janssen et al., 2016) and termed the X phase (the diluted star phase). At sufficiently large values of $\Gamma$ and $h$, even larger cluster ordering patterns like the $30$-site and $50$-site orders are stabilized. There is also an 18-site order with $C_{3}$ symmetry, which we label by $18$-$C_{3}$ to distinguish it from the previous $18$-site order as they are described by different arrangements of spins on the honeycomb lattice.

Next, we set $\Gamma^{\prime}=-0.02$ and map out the phase diagram within the same ranges of $\Gamma$ and $h$, as shown in Fig. 2. The addition of such a small $\Gamma^{\prime}$ term to the $K\Gamma$ model alters the phase diagram quite significantly. The degenerate manifold in the Kitaev limit and its neighborhood are replaced by the FM phase. The ZZ order is stabilized over a large portion of the parameter space at zeroRau and Kee and low fields. Once again, we find at intermediate fields several large cluster ordering patterns, a 32-site order, a 70-site order, a 98-site order and a 50-site order with $C_{3}$ symmetry which we label by $50$-$C_{3}$. Finally, the strong $\Gamma$ high field regime of the phase diagram displays some similarities to the $\Gamma^{\prime}=0$ case, where the same $50$-site, $18$-site and ZZ order are the lowest energy spin configurations, before the system becomes a FM.

Details of the magnetic orders (the real space spin configurations and the static spin structure factors, etc.) that show up in the phase diagrams Figs. 1 and 2, from the four sublattice ZZ order to the 98-site order, can be found in the Supplementary Materials. We make some qualitative observations as follows. Firstly, stronger $\Gamma$ interaction stabilizes magnetic orders with larger unit cells. This is true for both zero and finite $\Gamma^{\prime}$. We expect that ordering patterns with even larger unit cells than those mentioned above may appear if $\Gamma$ is further increased beyond $0.5$. Secondly, the large unit cell orders, like the $70$-site and $98$-site orders, are closely competing in the parameter region where they are stabilized. The difference in energy is typically $10^{-3}$ to $10^{-4}$ of the energies of these orders. Thirdly, the magnetic orders can be classified into two categories, one with an inversion symmetry and the other with a three fold rotational symmetry. The ZZ order and the magnetic orders labeled by numbers fall into the former, while the magnetic orders labeled by numbers appended with -$C_{3}$ fall into the latter. More details can be found in the Supplementary Materials.

Magnetization. The proposed spin model for the material $\alpha$-RuCl3 is parametrized by dominant $K$ and $\Gamma$ exchanges, with $K<0$ and $\Gamma\approx-K/2$, plus some small additional interactions like $\Gamma^{\prime}$, $J$ and $J_{3}$, where $J$ ($J_{3}$) is the (third) nearest neighbor Heisenberg exchangeKim and Kee (2016); Winter et al. (2016); Wang et al. (2017). Therefore, in the phase diagram Fig. 2 of the $K\Gamma\Gamma^{\prime}$ model in a $[111]$ magnetic field, we take a cut along $\Gamma=0.5$ and plot the magnetization $m=\mathbf{S}\cdot\hat{\mathbf{h}}$ as a function of the field $h$, as shown in Fig. 3a. The magnetization increases monotonically with the field, and jumps at the phase transitions. The discontinuities are not very obvious at the transitions between the large unit cell orders, but are significant when the system enters to (exit from) a large unit cell order from (to) ZZ, and from ZZ to FM. This suggests the difficulty of detecting phase transitions at intermediate fields by inspecting the magnetization, if they exist at all in the quantum model.

Linear spin wave theory. As a first approach to study the effect of quantum fluctuations on the classical orders, we apply the linear spin wave theoryJones et al. (1987); Cookmeyer and Moore (2018) to calculate the reduction of ordered moments in the zero temperature limit. For simplicity, we assume the same underlying magnetic orders, and do not consider how the classical phase diagram may be changed due to quantum correction to the energy because there are too many competing phases. Details of the linear spin wave calculation can be found in the Methods section. In Fig. 3b, we plot the average fraction of spins achieved in the linear spin wave theory with $S=1/2$ as a function of the field,

[TABLE]

where $b_{i}^{\dagger}$ is the magnon creation operator at site $i$ and $N_{\mathrm{site}}$ is the total number of sites in the system. Blank regions indicate that the spin wave Hamiltonian is not positive definite at one or more momenta, i.e. the lowest magnon band becomes gapless. At low and intermediate fields $h\lesssim 1$, the average reduction of ordered moments is about $50\%$ of the full magnitude $S$, hinting at strong quantum fluctuations. At high fields $h\gtrsim 1$, $\langle S_{\mathrm{LSWT}}\rangle/S$ increase monotonically with $h$ in the ZZ phase, but the spin wave solution becomes unstable in the region $h\in[1.29,1.37]$, where the system is in the FM phase with the spins not completely aligning with the $[111]$ field (see Fig. 3a). Not only is this region likely to host quantum spin liquids, but it is also interesting from the classical aspect, which we will discuss in details later. Finally, for $h>1.37$, the system enters the fully polarized state and $\langle S_{\mathrm{LSWT}}\rangle/S=1$ achieves saturation.

The spin wave dispersion of a (very) large unit cell order typically appears flat and dense. As an example, we show the spin wave dispersion of the $98$-site order along certain high symmetry directions in the first Brillouin zone of the honeycomb lattice in Fig. 4.

Thermal Hall conductivity. We calculate the thermal Hall conductivity due to magnonsKatsura et al. (2010); Matsumoto et al. (2014); Murakami and Okamoto (2017),

[TABLE]

Details of the calculation can be found in the Methods section. Expressing the field in (1) as $\mathbf{h}=g\mu_{\mathrm{B}}\mu_{0}\mathbf{H}$, assuming the $g$ factor $g=2.3$Cookmeyer and Moore (2018); Wu et al. (2018) and the magnitude of the Kitaev interaction $\lvert K\rvert\approx 80\,\mathrm{K}$Kim and Kee (2016); Kasahara et al. (2018a), $\mu_{0}H=12\,\mathrm{T}$ roughly corresponds to $h=0.23$. At this field and with the parametrization $(K,\Gamma,\Gamma^{\prime})=(-1,0.5,-0.02)$, the system is in the $50$-$C_{3}$ order. We plot the thermal Hall conductivity $\kappa_{xy}$ as a function of temperature $T$, as shown in Fig. 5a. We show only data below $T_{\mathrm{c}}\approx 11\,\mathrm{K}$, defined as the temperature at which $\langle S_{\mathrm{LSWT}}\rangle/S$ drops to zero, i.e. the magnetic order is destroyed by thermal fluctuations (see Fig. 5c). We find that $\kappa_{xy}$ is close to zero but slightly negative at $10\,\mathrm{K}$. It gradually develops a positive value as $T$ decreases, and peaks at $5\,\mathrm{K}$ before diminishing again as $T\longrightarrow 0$. The maximum value of $\kappa_{xy}/T$ is about $0.3\times 10^{-3}\,\mathrm{W}/\mathrm{K}^{2}\mathrm{m}$, which is of the same order of magnitude as the half quantized value $0.826\times 10^{-3}\,\mathrm{W}/\mathrm{K}^{2}\mathrm{m}$ measured in Ref. Kasahara et al., 2018b.

We also calculate the thermal Hall conductivity for another large unit cell order, the $32$-site order, at the field $h=0.82$ (which would roughly correspond to $\mu_{0}H=43\,\mathrm{T}$) and with the same parametrization, as shown in Fig. 5b. This time $T_{\mathrm{c}}\approx 6\,\mathrm{K}$ (see Fig. 5c) and $\kappa_{xy}$ is negative. Starting from zero temperature, $\kappa_{xy}$ grows in magnitude as $T$ increases, and reaches $-0.9\times 10^{-3}\,\mathrm{W}/\mathrm{Km}$ at $6\,\mathrm{K}$. The trend and the magnitude of $\kappa_{xy}$ are similar to those reported in Ref. Kasahara et al., 2018a at lower fields ($\mu_{0}H=6$, $12$ and $15\,\mathrm{T}$). Hence the opposite signs of $\kappa_{xy}$ may indicate the presence of different magnetic orders.

Frustrated ferromagnet. We notice that there is a window of $h$ where the system is a FM but not fully polarized, i.e. the spins align uniformly but not in the direction of the $[111]$ field. Such a phase is also stabilized in the high field regime at other parametrizations $(K,\Gamma,\Gamma^{\prime})$ including the $K\Gamma$ model, and the width of the window is usually larger for stronger $\Gamma$. In the following, we attempt to derive some analytical understanding of why this situation occurs.

We start from the $K\Gamma$ model with $K<0$, $\Gamma>0$, and $\Gamma^{\prime}=0$ in (1). Assuming a FM state, that is, $\mathbf{S}_{i}=\mathbf{S}$ for all sites $i$, the Hamiltonian reduces to

[TABLE]

with the matrices

[TABLE]

and $N$ being the total number of unit cells. The Kitaev interaction becomes “isotropic” in the FM state, behaving like the Heisenberg interaction. The $\Gamma$ interaction still appears quite anisotropic at this stage, but a change of basis will bring it to a simpler and more illuminating form. Switching from the cubic $xyz$ coordinates to the crystallographic $abc$ coordinates, where the $a$, $b$ and $c$ axes point in the directions $[11\bar{2}]$, $[\bar{1}10]$ and $[111]$ respectivelyChaloupka and Khaliullin (2015); Rusnačko et al. (2019), the spin is given by

[TABLE]

In the $abc$ basis, the Kitaev interaction $\tilde{H}_{K}=RH_{K}R^{\mathrm{T}}=H_{K}$ remains the same, while the $\Gamma$ interaction assumes the form of an $\mathrm{XXZ}$ model

[TABLE]

We can then analyze (4) in the $abc$ basis term by term. It can be shown analytically that the energy of the classical Kitaev model is $K\lvert\mathbf{S}\rvert^{2}$ per unit cellBaskaran et al. (2008). Thus any FM phase will minimize the energy of the $K$ term. On the other hand, the $\Gamma$ term attains maximum (minimum) when the spin points along (lies on) the $c$-axis ($ab$-plane). The energy profile of the $\Gamma$ term is shown in Fig. 6.

Suppose that the field $\tilde{\mathbf{h}}=(0,0,h)$ is along the $c$-axis or the $[111]$ direction. The $h$ term wants to align the spin with the $c$-axis, but this will be costly in energy for the $\Gamma$ term. The competition between $\Gamma$ and $h$ tilts the spin away from the $c$-axis. Therefore, such a FM phase can be stabilized between the fully polarized state and some other orders, e.g. ZZ and $18$, in the high field regime. In contrast, if the field is along any of the in-plane directions, then all the $K$, $\Gamma$ and $h$ terms in (4) can be minimized simultaneously.

Now let us consider the $K\Gamma\Gamma^{\prime}$ model. A finite $\Gamma^{\prime}$ term acts similarly as $\Gamma$. One can easily show that, assuming a FM state, $H_{\Gamma^{\prime}}$ has the same structure as $H_{\Gamma}$ in (5). Thus a small $\Gamma^{\prime}<0$ ($\Gamma^{\prime}>0$) weakens (enhances) the effect of $\Gamma$. A similar FM but not fully polarized phase due to the presence of a large $\Gamma$ interation in the $JK\Gamma J_{3}$ model under a $[001]$ field was also found and discussed in Ref. Janssen et al., 2017.

Discussion

Classcially, the pure Kitaev model is extremely sensitive to an external magnetic field. It is polarized whenever the field $h\neq 0$. From the result of our simulated annealing calculation, a finite $\Gamma$ interaction on top of $K$ gives rise to a multitude of ordered phases, many of them possess fairly large magnetic unit cells, at finite fields. As $\Gamma$ increases, the window of these nontrivial magnetic orders becomes wider, and the system becomes polarized at greater value of $h$. Thus the combination of $\Gamma$ and $h$ effects like a prism that produce a rich and colorful phase diagram. Adding a small $\Gamma^{\prime}$ term stabilizes even larger cluster magnetic orders at intermediate fields. We successfully demonstrate that the $K\Gamma\Gamma^{\prime}$ honeycomb model is a playground for many exotic field induced magnetic orders, not simply the zigzag (ZZ) order and the polarized state as largely perceived in the past.

We discuss the implications of these large unit cell orders. First of all, the size of the system has to be sufficiently large to host them. If the system is smaller than or incommensurate with the magnetic order, the ground state spin configuration may appear like a disordered state. This calls for a serious reconsideration of the results from quantum calculations on small systems where finite size effect can be important, such as ED on the $24$-site clusterGordon et al. and iDMRG on the cylinder geometryGohlke et al. (2018), which report a quantum spin liquid ground state. The large unit cell magnetic orders found in this work cannot be captured by these and similar computationsCatuneanu et al. (2018); Rusnačko et al. (2019); Jiang et al. ; Kaib et al. on quantum models with small system sizes. Nevertheless, the possibility of a quantum spin liquid still exists, especially in the vicinity of the large unit cell orders where quantum fluctuations are strong. The large unit cell orders are very close in energy in the parameter region where they are stabilized. In addition, the average spin wave correction to the ordered moments in the large unit cell orders for a representative parametrization of $\alpha$-RuCl3 is found to be more than $50\%$. One can imagine that quantum fluctuations may melt these competing orders and promote a spin liquid state, but we will not know whether this is true until the magnetic orders are explicitly taken into account in the quantum model. If the large unit cell orders (partially) survive under quantum fluctuations, the magnon bands are typically flat and very close to each other such that they appear like the excitation continuum seen in inelastic neutron scattering experiments, which is often interpreted as fractionalized excitations in a quantum spin liquidBanerjee et al. (2016, 2017); Winter et al. (2018); Banerjee et al. (2018); Balz et al. . The resulting two magnon excitations will also form a very broad continuum at low energies. Moreover, we calculate the magnon thermal Hall conductivity for two of the large unit cell orders and show that it resembles the trend and/or the magnitude as that measured in experimentsKasahara et al. (2018a, b) below the ordering temperature. In contrast, as computed in Ref. Cookmeyer and Moore, 2018, the magnon thermal Hall conductivity in the ZZ order is in general quite small in magnitude.

We also discover the existence of a ferromagnetic (FM) but not fully polarized state at high fields in the $K\Gamma\Gamma^{\prime}$ model with zero or small $\Gamma^{\prime}$, which can be understood through the competition between the $\Gamma$ and $h$ terms. Here $K<0$ and $\Gamma>0$ are assumed. While the field always wants to orient the spins in its direction, the $\Gamma$ interaction is only minimized (maximized) when the spins are all lying on the $ab$-plane (pointing along the $c$-axis). This may explain why the system is more prone to polarization when the tilting angle of the field from the $[111]$ direction is larger. This also suggests that frustration is stronger (weaker) when the field is along (in) the $c$-axis ($ab$-plane). For instance, the simulated annealing calculation on the classical $JK\Gamma J_{3}$ honeycomb model in an in-plane fieldLampen-Kelley et al. with the parametrization $(J,K,\Gamma,J_{3})=(-0.035,-1,0.5,0.035)$ only yield the $6$-site order (termed the X phase) at intermediate fields, between the ZZ order at low fields and the polarized state at high fields, thus leading to a relatively simple phase diagram.

Methods

Simulated annealing. The simulated annealing calculation is performed on a honeycomb lattice with $L\times L$ unit cells (or $L\times L\times 2$ sites) with periodic boundary conditions. Most of the computations are done with $L=12,15,20$, but sometimes other $L$ is used when the ground state spin configuration is not obvious. The procedure of simulated annealing is outlined as follows. In the beginning, we generate a totally random spin configuration on the honeycomb lattice and define a “temperature” parameter $T$. We randomly select a site on the honeycomb lattice and propose a random orientation for the spin on that site. Next, we calculate the difference in energy and accept the change with the probability $\min\{1,\exp(-\Delta E/T)\}$. This step is repeated for $\sim 10^{7}$ times at a fixed $T$, which is then decreased gradually. Once $T<T_{\mathrm{c}}$ for some critical temperature $\sim 10^{-8}\lvert K\rvert$, we update the spin at site $i$ deterministically by aligning it in the direction of the local fieldJanssen et al. (2016) defined as

[TABLE]

where $H_{ij}$ is the three dimensional matrix that encodes the interaction between the spins at $i$ and $j$.

If the sublattice structure of a magnetic order is known, we can carry out the above procedure for a small number of spins and calculate the energy to very high precision. This allows us to better determine the phase boundary between competing orders.

Linear spin wave theory. The content in this section is mainly derived from Ref. Jones et al., 1987. For each sublattice $i$ in the magnetic unit cell, we first choose a local coordinates system in which the spin $\mathbf{S}_{i}$ aligns in the $z$-direction. The amounts to a change of basis characterized by the rotation matrix

[TABLE]

where $\theta_{i}$ and $\phi_{i}$ are the two angles parametrizing the orientation of $\mathbf{S}_{i}$ in the cubic $xyz$ coordinates, $(S_{i}^{x},S_{i}^{y},S_{i}^{z})=S(\sin\theta_{i}\cos\phi_{i},\sin\theta_{i}\sin\phi_{i},\cos\theta_{i})$. The third column of $R_{i}$ is precisely $\mathbf{S}_{i}$ up to the factor $S$, while the first and second columns are chosen such that the three columns are mutually orthonormal and satisfy the right hand rule. We define $\mathbf{S}_{i}=R_{i}\tilde{\mathbf{S}}_{i}$. Classically, we have $\tilde{\mathbf{S}}_{i}=(0,0,S)$. Quantum effects on the ordered moments are introduced through spin wave excitations (magnons),

[TABLE]

where we have used the linear spin wave approximation that neglects the third and higher order terms in $b_{i}$ in the series expansion of (10b) and (10c). Next, we rewrite the spin Hamiltonian as

[TABLE]

where $\tilde{H}_{ij}=R_{i}^{\mathrm{T}}H_{ij}R_{j}$ and $\tilde{\mathbf{h}}_{i}=R_{i}\mathbf{h}$. Representing $\tilde{\mathbf{S}}_{i}$ using (10a)-(10c), keeping only terms quadratic in $b_{i}$, and performing a Fourier transform

[TABLE]

where, from now on, $i$ denotes the position in the Bravais lattice defined by the translational symmetries of the magnetic order, $s$ denotes the sublattice in the magnetic unit cell, and $N$ is the total number of magnetic unit cells. We then obtain the spin wave Hamiltonian in momentum space

[TABLE]

where $\Psi_{\mathbf{k}}=(b_{\mathbf{k},1},\ldots,b_{\mathbf{k},\mathcal{N}},b_{-\mathbf{k},1}^{\dagger},\ldots,b_{-\mathbf{k},\mathcal{N}}^{\dagger})$ and $\mathcal{N}$ is the total number of sublattices in the magnetic unit cell. $\mathrm{D}_{\mathbf{k}}$ is a $2\mathcal{N}$ dimensional matrix of the form

[TABLE]

where $\mathrm{A}_{\mathbf{k}}$ and $\mathrm{B}_{\mathbf{k}}$ are $\mathcal{N}$ dimensional matrices. To obtain the spin wave dispersion, we diagonalize $\mathrm{D}_{\mathbf{k}}$ by a Bogoliubov transformation in order to to preserve the canonical commutation relation of the bosons,

[TABLE]

where $\mathcal{E}_{\mathbf{k}}=\mathrm{diag}(\varepsilon_{\mathbf{k},1},\ldots,\varepsilon_{\mathbf{k},\mathcal{N}},\varepsilon_{-\mathbf{k},1},\ldots,\varepsilon_{-\mathbf{k},\mathcal{N}})$ and $\sigma^{3}$ is a diagonal matrix with the first $\mathcal{N}$ entries equal to $1$ and the last $\mathcal{N}$ entries equal to $-1$. The average reduction of ordered moments (10a) at temperature $T$ can be calculated from the matrix elements of the Bogoliubov transformation,

[TABLE]

where $f_{\mathrm{BE}}$ is the Bose-Einstein distribution,

[TABLE]

Thermal Hall conductivity. We explain the various symbols that appear in the formula (3) for the calculation of the thermal Hall conductivityMatsumoto et al. (2014); Murakami and Okamoto (2017). $n$ is the magnon band index that runs from $1$ to $\mathcal{N}$. The function $c_{2}$ is given by

[TABLE]

where $\mathrm{Li}_{2}$ is the dilogarithm. $f_{\mathrm{BE}}$ is the Bose-Einstein distribution as defined in (17). $\Omega_{\mathbf{k},n}$ is the Berry curvature defined as

[TABLE]

where $\sigma^{3}$ and $\mathrm{T}_{\mathbf{k}}$ are defined as in (15). For the calculation of the total volume $V$ of the system, we use the inter-plane distance $5.72$ Å  between the honeycomb layers in $\alpha$-RuCl3Kasahara et al. (2018a, b). The exact value of the in-plane lattice constant does not enter the calculation explicitly because, while $1/V$ contribute two inverse factors of it, $\Omega_{\mathbf{k},n}$ contribute two factors, so they cancel out. When performing the summation over momenta in (3), we partition the first Brillouin zone (of the magnetic order) evenly such that it contains a total number of $L\times L$ $\mathbf{k}$ points. We check the convergence of $\kappa_{xy}$ with increasing $L$ up to $L=800$. We also ensure that the Chern number of each magnon band,

[TABLE]

where $A$ is the total area of the system, converges to an integer with increasing $L$.

Data availability. The data that support the findings of this study are available from the corresponding author upon reasonable request.

Details of the Magnetic Orders

For each of the nontrivial magnetic orders, we show in Figs. S1-S14 the sublattice structure in the magnetic unit cell, the real space spin configuration (except for the $70$-site and $98$-site orders) and the static spin structure factor

[TABLE]

where $\mathbf{r}_{i}=\mathbf{R}_{i}+sa/\sqrt{3}\hat{\mathbf{y}}$, $\mathbf{R}_{i}$ is the coordinates of the Bravais lattice, $a$ is the lattice constant, $s=0$ ($s=1$) if $\mathbf{r}_{i}$ belongs to the even (odd) sublattice of the honeycomb lattice and $N_{\mathrm{site}}$ is the total number of sites in the system. We choose the direct lattice vectors to be $\mathbf{a}_{1}=a\hat{\mathbf{x}},\mathbf{a}_{2}=a(\cos\frac{\pi}{3}\hat{\mathbf{x}}+\sin\frac{\pi}{3}\hat{\mathbf{y}})$, so that the reciprocal lattice vectors are $\mathbf{b}_{1}=\frac{2\pi}{a}\hat{\mathbf{x}}-\frac{2\pi}{\sqrt{3}a}\hat{\mathbf{y}},\mathbf{b}_{2}=\frac{4\pi}{\sqrt{3}a}\hat{\mathbf{y}}$. The structure factors are plotted within the second Brillouin zone, and the first Brillouin zone is enclosed by dashed lines. Whenever our discussion involves more than one of the three inequivalent $\mathrm{M}$ points in the momentum space, we use the notation $\mathrm{M}_{i},i=1,2,3$ to distinguish them.

We observe that the magnetic orders respect either an inversion $\mathcal{I}$ symmetry or a three fold rotational $C_{3}$ symmetry. Therefore, the magnetic orders can be divided into two classes. In Class $\mathcal{I}$, if two sites $i$ and $j$ are related by inversion, then the spins on these sites align exactly in the same direction. $i$ and $j$ belong to distinct sublattices of the honeycomb lattice, either $i$ even and $j$ odd, or $i$ odd and $j$ even. The inversion only acts in real space, i.e. it does not flip the direction of the spin. In Class $C_{3}$, the spin configuration remains invariant under a $C_{3}$ rotation about the $[111]$ direction, with the axis of rotation piercing a site. The rotation takes place in both the real space and the spin space,

[TABLE]

In the spin space, $C_{3}$ permutes the $x$, $y$ and $z$ components of the spin. We illustrate the $\mathcal{I}$ and $C_{3}$ symmetries using the $18$-site and $18$-$C_{3}$ orders as examples in Figs. S5b and S7b. In addition, each magnetic order in Class $\mathcal{I}$ has three different domains with the same energy, which are related by a $C_{3}$ rotation in both the real space and the spin space. Consequently, the profile of the structure factors corresponding to these three domains differs by a $C_{3}$ rotation. For clarity, however, there is no $C_{3}$ symmetry in each domain. Each magnetic order in Class $C_{3}$ has only one domain by definition.

Furthermore, we notice that there is perhaps some kind of number rule governing the size of the magnetic unit cells, which has yet to be understood. The $98$-site order appears like an augmented version of the $50$-site order, which in turn appears like an augmented version of the $18$-site order (see Figs. S5b and S10b, for instance). Their respective magnetic unit cells contain $7\times 7$, $5\times 5$ and $3\times 3$ unit cells of the honeycomb lattice. Besides, the $6$-site, $30$-site and $70$-site order, whose magnetic unit cells contain $3\times 1$, $5\times 3$ and $7\times 5$ physical unit cells, appear like parts of the $18$-site, $50$-site and $98$-site orders, respectively.