The Heart of Entanglement: Chiral, Nematic, and Incommensurate Phases in the Kitaev-Gamma Ladder in a Field

The bond-dependent Kitaev model on the honeycomb lattice with anyonic excitations has recently attracted considerable attention. However, in solid state materials other spin interactions are present, and among such additional interactions, the off-diagonal symmetric Gamma interaction, another type of bond-dependent term, has been particularly challenging to fully understand. A minimal Kitaev-Gamma (KG) model has been investigated by various numerical techniques under a magnetic field, but definite conclusions about field-induced spin liquids remain elusive and one reason may lie in the limited sizes of the two-dimensional geometry it is possible to access numerically. We therefore focus on the KG model defined on a two-leg ladder which is much more amenable to a complete study, and determine the entire phase diagram in the presence of a magnetic field along [111]-direction. Due to the competition between the interactions and the field, an extremely rich phase diagram emerges with fifteen distinct phases. Focusing near the antiferromagnetic Kitaev region, we identify nine different phases solely within this region: several incommensurate magnetically ordered phases, spin nematic, and two chiral phases with enhanced entanglement. Of particular interest is a highly entangled phase with staggered chirality with zero-net flux occurring at intermediate field, which along with its companion phases outline a heart-shaped region of high entanglement, the heart of entanglement. We compare our results for the ladder with a C3 symmetric cluster of the two-dimensional honeycomb lattice, and offer insights into possible spin liquids in the two-dimensional limit.


I. INTRODUCTION
The Kitaev model on a two-dimensional honeycomb lattice is a rare example of an exactly solvable model offering a quantum spin liquid with fractional excitations. 1 Under a time-reversal symmetry breaking field, it exhibits non-Abelian anyons with half-quantized thermal Hall conductivity originated from Majorana edge mode. Since the original proposal, finding a solid state material possessing such a quantum spin liquid has attracted great attention. A microscopic mechanism for realizing the Kitaev model in solid-state material was first suggested using the combined effects of strong spinorbit coupling and electron-electron interactions 2,3 . Later, the nearest neighbor generic spin model on an ideal honeycomb lattice was re-derived, and it was found that there are additional bond-dependent interactions present with the so-called Gamma (Γ) interaction among the most intriguing. 4 From the material perspective, α-RuCl 3 was proposed as a leading candidate with a weaker coupling between layers making the material close to two dimensional (2D). [5][6][7][8][9][10] Furthermore, the Γ interaction has been found to be as large as the Kitaev interaction in α-RuCl 3 . [11][12][13] Since then, RuCl 3 has been explored by several experimental and theoretical techniques. [14][15][16][17][18] In particular, early inelastic neutron scattering 10 and Raman spectroscopy 19 measurements have suggested a strong frustration well above the magnetic ordering temperature, indicating strong frustration which may originate from the bond-dependent Kitaev and Γ interactions. Remarkably, a half-quantized thermal Hall conductivity was recently reported in α-RuCl 3 in a certain range of the magnetic field 20 when a zig-zag magnetically ordered phase is destroyed by the magnetic field. 21 Other physical quantities accessed by several experimental techniques in RuCl 3 also suggested that there is a nontrivial intermediate phase under a magnetic field which is different from a trivially polarized spin state in high field regime. [22][23][24][25][26][27][28][29][30][31] However, the experimental evidence for a field-induced intermediate disordered phase in RuCl 3 is still under debate. [32][33][34] In parallel to the experimental progress, theoretical attempts to find nontrivial field-induced phases in extended Kitaev model have been pursued extensively. 12,[35][36][37][38][39][40][41][42][43][44][45][46][47][48] Most numerical studies are limited to either near the antiferromagnetic (AFM) Kitaev or near the ferromagnetic (FM) Kitaev region, as the exactly solvable Kitaev point offers a starting point. In particular, a minimal Kitaev-Gamma (KG) model under a [111] magnetic field has been studied near FM Kitaev and AFM Γ region relevant to RuCl 3 . 46-49 A 24-site exact diagonalization (ED) study showed a field-revealed Kitaev spin liquid near FM Kitaev region with a finite AFM Γ interaction, when the magnetic field is tilted away from the [111]-axis. 49 However, the infinite tensor product state (iTPS) found a small confined Kitaev spin liquid, and broken C 3 rotational phases are induced under the magnetic field. 47 Interestingly various large unit cell magnetic orderings have been reported in the classical KG model under the magnetic field along [111]axis. 46 , which are replaced by these broken rotational phases in quantum model. 48 Whether the C 3 broken phases are quantum spin liquids or not is at present not clear and will require further studies.
Numerical studies near AFM Kitaev limit under a magnetic field have found intriguing results. 37,[39][40][41][42]45,[50][51][52] It was suggested that a gapless U(1) spin liquid is induced by the magnetic field. 45 The energy spectra obtained by 24-site ED showed putative gapless excitations in the intermediate field, which then transition to a polarized state (PS) in high field regime. Several iDMRG studies also reported a change of central charge depending on the number of legs in the DMRG which indicates a finite spinon Fermi surface. 50 However one may question if the dense energy spectra are due to incommensurate order which is difficult to detect due to the finite size of ED and limited access to momentum points in iDMRG. Indeed, different iDMRG studies have found different gapless points in momentum space. 39,50,52 Despite intensive studies, definite conclusions on possible phases and the nature of numerically determined phases near the Kitaev regions remain indefinite. One reason for the controversial results among the previous studies may lie in the limited sizes of the two-dimensional honeycomb geometry that one can access numerically. Furthermore, the zerofield and field-induced phases of the entire phase space of KG model are yet to be determined. We therefore focus on the KG model defined on a two-leg ladder which is much more amenable to a detailed study and a complete phase diagram in the presence of a magnetic field along [111]-direction can be determined. Using high through-put iDMRG calculations we map out the entire phase diagram of the KG ladder which shows an extremely rich structure. After we determine the entire phase diagram, we focus on the region of the phase diagram where the Kitaev interaction is predominantly AFM. In this region, in the absence of a magnetic field, we identify a novel spin nematic (SN) phase with quadropolar order in addition to the AFM Kitaev phase (AK) and a phase connected to the isotropic FM ladder through a local 6-site spin rotation, FM U6 . [53][54][55][56] In zero field the KG ladder can be mapped to a ladder with four-spin exchange closely related to JQ models 57 extensively studied as models of deconfined criticality 58 .
When a magnetic field in the [111] direction is introduced, field and spin interactions compete, and a proliferation of phases is observed. We identify phases with scalar chiral ordering and several phases with magnetic ordering some of which might display incommensurate or very large unit cell ordering. Of particular interest is two chiral ordered phases characterized by a staggered chirality (SC) and uniform chirality (UC). The SC phase is a magnetically disordered and highly entangled phase occurring at intermediate magnetic field above the AK phase. It has the staggered chirality with zero-net flux despite it is under a rather large external field, and shows clear chain end excitations. Rather poetically, this phase along with its companion phases outline a heartshaped region of high entanglement, 'the chiral heart of entanglement'. On the other hand, the UC phase with the uniform chirality leading to a finite net flux appears between the SN and and the rung singlet (RS U6 ) phase connected to the isotropic AFM ladder through a local 6-site spin rotation. [53][54][55][56] It emerges at extremely low field, as if three phases, SN, UC, and RS U6 may meet at a critical point. All together, near the AFM Kitaev region alone, we identify 9 possible distinct phases in addition to the FM U6 phase and the PS occurring at high magnetic fields.
Below, in section II, we first review the KG model. In addition to the extensive 2D cluster studies, a one-dimensional (1D) KG chain model including only xand y-bonds was studied using non-Abelian bosonization and DMRG and reported SU(2) emergent phases. 59 . The two-leg ladder is made of two such KG chains by connecting them by z-bond. Technical details relevant fo the iDMRG and DMRG numerical methods are subsequently discussed in the following section, III. An overview and detailed discussion of the full phase diagram of the ladder is presented in section IV along with a discussion of the connection between the KG chain and the ladder. In section V we focus on the vicinity of the AFM Kitaev region under a magnetic field, where a rich phase diagram with various phases with enhanced entanglement is found. Finally, in section VI we compare our results for the ladder with 24-site ED results obtained in the honeycomb geometry, and discuss the implications of the ladder results to the 2D limit.

II. THE TWO-LEG KITAEV-GAMMA (KG) LADDER
The two-leg KG ladder is formed out of a strip of the honeycomb lattice. The two KG chains with only xand y-bonds are coupled by adding the z-bond as shown in Fig. 1. Periodic boundary conditions in the direction perpendicular to the ladder is then imposed by directly coupling the dangling z-bonds thereby forming a regular ladder. Sometimes, the additional z-bonds from imposing periodic boundary conditions are taken to be of opposite sign (and or strength) in which case the resulting model is usually referred to as a honeycomb ladder 60 . Here, all z-bonds are identical and a regular ladder is formed. In addition to the bond-dependent Kitaev interaction, the KG Hamiltonian incorporates another bond-dependent interaction, Γ. 4 For the KG ladder we orient the bonds so that the Kitaev z-bond connects the two legs of the ladder as shown in Fig. 1. The complete Hamiltonian is then given by where (α, β) takes on the values (y, z)/(x, z)/(x, y) for γ = x/y/z, and i, j refers to the nearest neighbor sites. We keep K = cos φ and Γ = sin φ and interpolate between the Kitaev ladder and Γ ladder by varying φ from 0 to 2π. We denote the total number of sites in the ladder (including both legs) by N . The pure Kitaev ladder at φ = 0, π is exactly solvable 61 and at both points it is in a disordered gapped phase 61 . However, recently it was shown that a non-local string order parameter (SOP) can be defined 62 at the Kitaev points. The SOP remains non-zero in the presence of a small Heisenberg coupling, J, at both the Kitaev points. Close to the FM Kitaev point, φ = π, the phase diagram of the KG ladder has recently been investigated in the presence of a magnetic field and additional interactions 49 . On the other hand, relatively little is known about the rest of the phase diagram of the KG ladder which is our focus here. In one dimension in the absence of a magnetic field, the closely related KG chain has been investigated in considerable detail 59,63 . An extended disordered phase close to the AFM Kitaev point, φ = 0 has been identified along with an adjacent spin nematic phase. For the KG chain it can rigorously be established that the phase diagram is symmetric with respect to Γ → −Γ, a symmetry that is clearly absent in the KG ladder.
It is of particular interest to also consider the effect of a magnetic field. Here we exclusively consider a field in the [111] direction, perpendicular to the honeycomb plane of the ladder. The magnetic field leads to a Zeeman coupling as where we choose the direction of h along the [111]-axis normalized as h = h √ 3 (1, 1, 1), g = 1, and S = σ 2 with σ is Pauli matrix. We use units with = 1 and µ B = 1.

III. NUMERICAL METHODS
As our main tools for investigating the KG ladder we use exact diagonalization (ED), finite size density matrix renormalization group (DMRG) and infinite DMRG (iDMRG) techniques. The iDMRG calculations are performed in two different ways. A high through-put mode with a small unit cell of size 24 or 60 and a maximal bond dimension of 500 and a high precision mode with a unit cell of 60 and a maximal bond dimension of 1000. Typical precisions for the two iDMRG modes are = 10 −8 and = 10 −10 , respectively, and = 10 −10 for the finite size DMRG calculations. Finite size DMRG calculations are performed both with open boundary conditions (OBC) and for smaller system sizes with periodic boundary conditions (PBC). In order to establish the phase diagram we focus on several different characteristics. With e 0 the ground-state energy per spin we define the energy susceptibilities where χ e h could equally well be called a magnetic susceptibility. For finite size systems of size N it has been established 64 that the energy susceptibility at a quantum critical point (QCP) diverges as Here ν and z are the correlation and dynamical critical exponents and d is the dimension. It follows that χ e may not necessarily detect the phase transition if the critical exponent ν is sufficiently large. We have therefore found it useful to supplement the analysis of χ e by a study of the entanglement spectrum 65 (ES) at different sections of the ladder. We have found it most useful to use a partition of the ladder where a cut is introduced at site N/2 − 1 thereby intersecting a rung. With λ α the eigenvalues of the resulting reduced density matrix, ρ N/2−1 the entanglement spectrum can be defined as − ln λ α and yields a characteristic signature of a phase. Close to the QCP the λ α rapidly change whereas they remain approximately constant inside a phase. We therefore focus on the largest of the λ α 's which we denote by λ 1 and study the lowest edge of the entanglement spectrum defined by − ln λ 1 .
In a product phase λ 1 = 1 (− ln λ 1 = 0) and such phases, with zero entanglement, are therefore easily detected by tracing out − ln λ 1 . Conversely, phases with high entanglement will have − ln λ 1 0 and of course, due to the ordering of the λ α one must have − ln λ α > − ln λ 1 for any α > 1. Sometimes changes in − ln λ 1 are imperceptible and we have therefore found it useful to define an entanglement spectrum susceptibility as follows Since the ES has to change at the QCP, χ λ1 should be able to detect any phase transition with the exception of unlikely scenario's where λ 1 only changes linearly at the QCP with all non-trivial changes in the higher λ α 's. We have found χ λ1 to be an extremely sensitive measure, often changing many orders of magnitude at a QCP, and we therefore typically focus on ln χ λ1 .

IV. FULL KG LADDER φ, h[111] PHASE DIAGRAM
We start with a discussion of the full phase diagram covering the entire range φ ∈ [0, 2π] and h[111] ∈ [0, 1.75]. Our results for the full phase diagram as obtained from iDMRG calculations are shown in Fig. 2. We first show ln χ λ1 φ in Fig. 2(a). The divergence of χ λ1 φ at a phase transition is so strong that it is most sensible to plot ln χ λ1 φ . A well defined phase, where λ 1 is close to constant is then visible in Fig. 2(a) as a dark blue coloring. On the other hand, a divergent χ λ1 φ , indicating a phase transition is visible as a dark red color. As is clearly evident from Fig. 2(a) the complexity of the phase diagram due to the many competing phases is truly remarkable. Fig. 2(b) and (c) show χ e φ and − ln λ 1 , and the spin excitation gap at zero field, respectively. We will discuss them in detail later.
A second view of the full phase diagram is shown in Fig. 2(d) where the bipartite von Neumann entanglement entropy S rung is shown. We define S rung = −Trρ N/2−1 ln ρ N/2−1 , where ρ N/2−1 is the reduced density matrix obtained from a partition of the ladder after site N/2 − 1, a partition that will cut a rung in the ladder. Highly entangled phases are here visible as bright yellow colors where as phases with negligible or no entanglement, where the ground-state wave function is well described by a simple product form, are shown as dark blue. We note that considerable scattering are clearly visible in certain regions at finite fields, for instance above φ = π/4 and φ = 0. The noise is due to poor convergence of the iDMRG due to the high frustration present. In these regions a more careful analysis with either exact diagonalizaiton or finite-size DMRG is necessary.

A. Zero field phase diagram
Let us now focus on the phase diagram in zero field. A detailed high precision calculation of χ e φ at zero field is shown in Fig. 2(b) along with the lower edge of the entanglement spectrum (ES), − ln λ 1 . A large number of well defined phase transitions are clearly visible which we now discuss.
Part of the phase diagram close to the FM Kitaev point φ = π, has previously been discussed 49 and in zero field a gapped spin liquid phase denoted KSL between φ/π = 1 and φ/π = 0.883 along with a second gapped phase KΓSL starting below φ/π = 0.883 have been identified. Since we here discuss the full phase diagram we shall refer to the KSL phase as FK and the KΓSL phase as AΓ to distinguish them from the phases ocurring at the AFM Kitaev point. The notation AΓ makes sense since this phase surrounds φ = π/2 where K = 0, Γ = 1.
A local unitary transformation, U 6 , is also known 59,63,66 . The U 6 transformation locally rotates the spins in a manner so that the Γ couplings are transformed into Heisenberg like (xx, yy, or zz) couplings with a changed sign. The transformation can be applied equally well to the chain, the ladder and the honeycomb plane. At the points φ = π/4 and φ = 5π/4 where the Kitaev and Γ couplings are of equal strength the KG ladder is therefore transformed into an isotropic FM and AFM Heisenberg ladder. These two points therefore have hidden SU(2) symmetry. It is well established that the isotropic AFM Heisenberg ladder is in a gapped disordered rung singlet (RS) phase 67,68 and we therefore denote the corresponding phase for the KG ladder as RS U6 . In the two-dimensional honeycomb lattice limit the RS U6 phase becomes a 120 degree ordered phase. 4 For the FM point, φ = π/4 we denote the magnetically ordered gapless phase by FM U6 . Since the FM U6 is well approximated by a product wave function with negligible entanglement it is distinctly visible in Fig. 2(d) with its almost black coloring. The FM and AFM points are shown as solid red circles along the φ-axis in Fig. 2. The same unitary U 6 transformation transforms the FM Kitaev point, φ = π, to the AFM Kitaev point, φ = 0. The energy spectrum at these two points must therefore be identical, a property that does not hold for any non-zero Γ.
Starting from right to left we observe the transition from RS U6 to the FK phase occurs at φ/π = 1 with the subsequent transition from the FK phase to the AΓ phase occurring at φ/π = 0.883. Between the gapless FM U6 phase and the gapped AΓ phase we observe a rapid sequence of several well defined phase transitions at φ/π = 0.440, φ/π = 0.428, φ/π = 0.396 and finally at φ/π = 0.385. See the zoomed inset in Fig. 2(b). While χ e φ only detects a single transition at φ/π = 0.385 the other 3 transitions are clearly identifiable in − ln λ 1 as shown in the zoomed inset. The precise nature of the intervening phase is at present unclear and left for future study. A clear transition out of the FM U6 phase to the gapped spin liquid phase AK is observed at φ/π = 0.086. As previously mentioned, precisely at φ = 0 a non-local string order has been found 62 and the AK phase is clearly identifiable as a spin liquid phase. In the 2D honeycomb lattice limit the AK phase becomes the AFM Kitaev spin liquid. We shall discuss the AK phase further below. We now turn to a discussion of the last phase observed in zero field, to the left of the AK phase.

Nematic Phase in zero Field, SN
In a recent study 63 the KG chain was investigated and a phase with spin-nematic (spin-quadropole) order adjacent to the spin-liquid phase at φ = 0 was identified. From a symmetry analysis of the following 4 order parameters were identified: The two order parameters Q c and Q d describe off-diagonal ordering between sites not coupled by the same terms in the Hamiltonian and Q e , Q f diagonal ordering again between sites not coupled the same way in the Hamiltonian. The above definitions therefore depend on a specific ordering of the couplings along the leg. Considering only the Kitaev coupling, site 1,2 would be xx coupled, site 2,3 yy coupled and so forth. (Note that along the legs of the KG ladder no S z S z coupling occurs so Q f can be defined between any two nearest neighbor sites). We can use the same order parameters to study nematic ordering along the legs of the KG ladder. Our results are shown in Fig. 3 as obtained from iDMRG. Due to the translational invariance the Q's are the same among all sites and are easily calculated. They all four become non-zero at φ/π = −0.155 which coincides with divergences in χ e φ and χ λ1 φ . The transition by varying φ is also clearly visible directly in − ln λ 1 as can be seen in Fig. 2 However, this is not a conventional spin-quadropole phase. The DRMG with OBC shows two different magnetic orderings depending on the size of the system. One has the AFM order along the leg and FM between the rung, and the other has a 6-site ordering mapping to AFM order after 6-site transformation, while the iDMRG finds the first one. The AFM order corresponds to a stripe order, if we continue the ordering pattern by increasing the number of ladders to the the 2D limit. Let us briefly consider what happens if an additional Heisenberg coupling J is introduced along side the K and Γ terms. In that case, the striped phase appears for an AF Heisenberg interaction J > 0 in 24-site ED calculations on the C 3 cluster, while the second ordering is likely a spiral order occurring for J < 0. 4 This suggests that this particular window of φ with J = 0 is in fact a line of first order transitions (in J) between these two orderings. To confirm such a possibility, we have studied the phase boundary by sweeping J ( parameterized as J ≡ K cos θ). Indeed we find a clear first order transition occuring at J = 0 in χ e θ , the second derivative of ground state energy per spin with respect to θ, as shown Appendix B.
While it is a line of first order transitions, we would refer to this as SN for spin-nematic, as the SN is a common feature of these orderings coexiting along the transition line. As the magnetic field becomes finite, it develops a magnetic order with almost zero entanglement which is shown later. The transition out of the SN phase into the RS U6 phase occurs at around φ/π = −0.265. As we shall discuss later, this quantum phase transition is actually a multi-critical point where the first order transition line ends, and two other phases occur. A high precision determination of the location of the critical point is therefore significantly more difficult than for the other QCP's. In particular so, since the entanglement for φ <φ/π = −0.265 is exceedingly high. It is therefore shown as a broader dashed line in Fig. 3.

Excitation gap in Zero field
In Fig. 2(c) we show the spin excitation gap to the first 4 lowest lying states, ∆ 1 , ∆ 2 , ∆ 3 and ∆ 4 as obtained from high precision finite size DMRG on ladders with N = 60. We have verified that finite size effects are relatively small if not in the proximity of a QCP. This is indicated by the solid blue circles in Fig. 2(c) which indicate extrapolations to N = ∞ of δ 1 by fitting data for N = 24, 36, 48, 60 and 72 to the form ∆(L) = ∆ ∞ + a exp(−L/ξ)/L.
Starting from the right, we find that the RS U6 phase, as expected, has a single ground-state with a well defined triplet excitation through-out most of the phase. The triplet excitation merges with higher lying excitations at φ = 1.56π. Precisely at the FM Kitaev point, φ = π, there is a level crossing leading to a clear first order transition. This is exact also for finite systems. The same holds true by symmetry at the AFM Kitaev point φ = 0. However, while the FK phase has a single ground-state below a well defined gap and φ = π is a transition point, in the AK phase the ground-state remains double degenerate below a gap, only split by finite-size effects. The AΓ phase also has a single ground state with a well defined triplet excitations above it. The FM U6 phase is gapless, but occasionally excited state DMRG calculations get trapped in higher lying states and the gapless nature of the phase does not appear clearly in Fig. 2(c). The final phase we discuss here, the SN phase occurring between φ/π = −0.265 and φ/π = −0.155 is clearly gapless as can be seen in Fig. 2(c) with all four gaps close to zero.

B. Relation to 1D KG Chain
As eluded to above, the phase diagram of the ladder is closely connected to that of the KG chain and the KG model defined on the full two-dimensional honeycomb lattice.
Comparing to the phase diagram of the KG chain 59,63 , a gapped FK phase appeared in the ladder, while it was gapless in the chain model. Similarly, the AK persists near AFM Kitaev region, and it is gapped in the ladder, while gapless in the chain. RS U6 is gapped in the ladder similar to AFM Heisenberg model, where as it was gapless in the chain. FM U6 remains magnetically ordered with gapless excitations as was the case for the chain. For the chain a gapless nematic phase was identified 63 on either side of the AK phase. For the ladder a similar gapless phase, SN, occurs but this time only on one side of the AK phase. AΓ remains disordered like the chain. There is therefore a close connection between the phases identified in the KG chain but we expect the KG ladder to be much closer to the two-dimensional honeycomb lattice and to represent most of the phases occurring in that limit although we in some cases expect phases to become gapless in the twodimensional limit. For instance, in the KG ladder AK and FK phases are gapped and we expect these phases to become the AFM and FM Kitaev spin liquid, respectively in the 2D limit, which are gapless.

C. Non-zero Field
When a magnetic field in the [111] direction is introduced, the full h[111], φ phase diagram is revealed as shown in Fig. 2(a) and (d). From the rung entanglement, S rung , shown in Fig. 2(d), it is clear that the introduction of a field in many cases tend to increase the entanglement. New highly entangled phases appear until the fully field polarized state is attained where all spins are polarized along h[111]. We denote this polarized state by PS. Trivially, it is a product state with zero entanglement. We start by discussing the fate of the RS U6 phase. This phase is adiabatically connected to the PS phase and no phase transition is observed at any nonzero field strength. This is contrary to the FM U6 phase state which cannot be adiabatically connected to the PS phase. This follows from the fact that the phase is not an ordinary FM state but only related to one through the local unitary rotation U 6 . An alignment of the spins along the [111] direction is therefore energetically costly. As can be seen in Fig. 2(a) and (d) a line of phase transitions around field strengths of h[111] = 0.7 − 1.0 occur, signalling the transition to the PS phase. A new phase at finite field above φ = π/2 is clearly visible. This is a gapless, magnetically ordered phase with the spins arranged in a zig-zag manner, in opposite directions on the two legs of the ladder, and we refer to this phase as the ZZ phase. The fate of the phases occurring between φ/π = 0.440 and φ/π = 0.385 is not known and left for future study. The FK and AΓ phases survive in the presence of a magnetic field and survive up to large field strengths. The nature of these two phases in the presence of h[111] has previously been discussed 49 . We therefore leave that part of the phase diagram aside an instead concentrate on the part of the phase diagram close to the AFM Kitaev point, φ = 0.
In the vicinity of φ = 0 it is clear from Fig. 2(a) and (d) that several new phase are induced by the magnetic field, in many cases with significantly increased entanglement. The nematic phase, SN, identified in zero field, survive in the presence of a non-zero field and is clearly visible in Fig. 2(a) and (d). However, the high through-put iDMRG calculations used in Fig. 2 are in certain regions having trouble achieving convergence as can be seen by the noise in the figure and a precise determination of the phase diagram in the vicinity of φ = 0 is difficult from the data presented in Fig. 2(a) and (d). We therefore focus on a high resolution study of that part of the phase diagram combined with finite-size DMRG calculations for the regions where it is not possible to achieve good convergence using iDMRG.
The most exotic part of the phase diagram of the KG model under the field is the AFM Kitaev region around φ = 0 where |Γ| K. In zero field we have identified the AK spin liquid phase and the nematic phase SN. However, it is clear that there is a very fine balance among the different couplings and the presence of a magnetic field will substantially increase this competition. We therefore investigate this part of the phase diagram with extreme care. Our results are presented below.
In order to get a detailed picture of the phase diagram we have performed high through-put iDMRG calculations in the region φ ∈ [−0. 35 By necessity, we have to use a relatively small unit cell and a maximal bond dimension of 500. Our results are shown in Fig. 4. In panel (a) is shown the lower edge of the entanglement spectrum, − ln λ 1 , determined from the reduced density matrix that cuts a rung. Phases with negligible entanglement (λ 1 ∼ 1) will show as dark blue where as highly entangled phases (λ 1 1) are colored yellow to dark green. The AK phase is clearly visible as the bright yellow phase centered around φ = 0, extending to non-zero fields. On the other hand, the SN phase which is close to a product state is clearly visible as dark blue. In addition to these two previously discussed phases there is a proliferation of phases occurring at higher fields. In order to more clearly identify phase transitions we show χ λ1 φ in Fig. 4(b) on a logarithmic scale. Deep blue coloring in Fig. 4(b) corresponds to a stable λ 1 , a well defined phase, while dark red coloring signals rapid change in λ 1 and a likely associated phase transition. Note the logarithmic scale, where the darkest red coloring is more than 10 orders of magnitude larger than the blue colors. A significant advantage analyzing the phase diagram in the way shown in Fig. 4(a) and (b), is that the entanglement spectrum, and therefore also λ 1 , has to change at a quantum phase transition.
In Fig. 4(b) many clear phase transitions are visible as dark red lines. However, there are also some extended regions with dark red coloring or noise appearing. These regions marked, η, τ and µ in Fig. 4(a), (b), are regions where the small unit cell idmrg calculations have difficulty reaching good convergence. Likely this is due to incommensurability effects and a further investigation using finite size DMRG or ED is warranted. From the data in Fig. 4(a), (b) we identify six well defined phases, the previously discussed AK and SN phases and 4 new phases that we name SC, β, γ and UC. Here SC and UC refer staggered chirality (SC) and uniform chirality (UC) because these phases exhibit a staggered and uniform pattern of chirality without any magnetic ordering, respectively. As we discuss in more detail below, it appears that the µ region is distinct from the SC phase and it seems quite plausible that the η, τ and µ regions are in fact well-defined phases and we therefore discuss them as such below. However, at present we cannot exclude the possibility that, for instance, what appears as a phase transition between the SN phase and the µ phase is instead a 'disorder' line marking the onset of incommensurate short range correlations thereby hindering the convergence of the iDMRG calculations. It is therefore possible that the η, τ and µ regions are not distinct phases but simply parts of the adjoining phases where short-range correlations are different. We shall return to this point below. Surprisingly, it is clear that some of the phases occurring at finite field have significantly increased entanglement, most notably the SC phase that together with its accompanying phases, γ and µ outline a heart shaped region of extraordinary high entanglement.
Complementary views of the same phase diagram as obtained from S rung and χ e h can be found in Appendix A.
A. Overview of phases, S α i Several of the phases visible in Fig. 4(a) and (b) are magnetically ordered phases. We have therefore performed high precision finite size DMRG calculations with N = 400 and OBC at the points marked with a red × in Fig.. 4(a). Results for the local magnetization S α i are shown in Fig. 5 and 6, here the green circles are S x i , the red circles S y i and the blue circles S z i along a single leg of the ladder. In some cases we find magnetic ordering with a surprisingly large unit cell, in other cases indications of incommensurate ordering. Note that a simple polarization of the spins along the field direction h[111] would have all S α i equal. Before a more detailed discussion of some of the phases we summarize some of the main findings in Fig. 5 and 6.

Magnetically ordered phases
We start with the six phases that show clear signs of ordering.
• η: This is a high field phase (region). As shown in Fig. 5(a) the local magnetization appear incommensurate. The ground-state is degenerate and the phase is likely gapless.
• β: As Γ is made slightly more negative the system transitions from the η phase to the β phase. The local ground-state magnetization along the leg of the ladder, shown in Fig. 5(b), displays a characteristic Möbius form, showing a single twist in S x i and S y i from one end of the open chain to the other while S z i remains constant. The ground-state is degenerate and the phase is possibly gapless.
• τ : This phase (region) is adjacent to the β phase but the beautiful intricate ordering along the leg shown in Fig. 5(c) is clearly distinct from that observed in the β phase with a large variation in S z i along the leg. The unit cell for the ordering appears at this value of Γ to be approximately 20 lattice spacings along the ladder leg. Again we find a degenerate ground-state and a possible gapless phase. At neighboring values of Γ we find similar large unit cell ordering.  • SN: In zero field the SN was clearly identified. There is no indication of a phase transition as the h[111] is introduced and we therefore assume that the well defined uniform phase visible in Fig. 4(a) and (b) is adiabatically connected to the SN phase in zero field. The phase is gapless and the presence of the non-zero magnetic field induces an incommensurate local magnetization as shown in Fig. 5(f).
As outlined above, the large unit cell ordering occurring in the µ and τ phases (regions) varies with Γ and h [111]. The same effect is observed for the incommensurate ordering in the η phase. In fact, the stripes occurring in these phases visible in Fig. 4(a) are likely caused by the variations in the local ordering best compatible with the 24 site unit cell. The lines could therefore represent lock in transitions. These phases are therefore not uniform in a conventional sense but could possibly be a series of phases with shifting sizes of unit cells for the magnetic ordering. We have not been able to resolve this.

Non magnetic phases
We now turn to a discussion of the three remaining disordered phases which do not show any conventional local magnetic ordering apart from that induced by the magnetic field.
• SC: In Fig. 6(a)  As φ is increased from -0.08 to -0.10 approaching the µ phase the size of the chain end excitations visibly grow. The SC phase is highly entangled, significanly more so than the AK phase. The ground-state does not appear degenerate but at h[111] = 0.38, φ = −0.05 we can limit the gap to the first excited state by ∆ 1 < 0.019. The gap is likely smaller in other parts of the phase. The bipartite entanglement entropy S(x) is close to constant throughout most of the ladder which is also consistent with a gapped phase (see Appendix C). The phase show signs of scalar chiral ordering as we discuss in more detail below.
• AK: This is the AFM Kitaev phase previously discussed. As shown in In Fig. 6(b) the local magnetization is aligned with the h[111] field and only faint signs of chain end excitations are visible. The phase has a gap and in zero field a two fold degenerate ground-state. It is possible to define a string order parameter (SOP) 62 as we discuss below.
• UC: On the left side of the SN phase a new phase appears in the presence of a h[111] field. As we shall discuss below this phase has scalar chiral ordering. The ground-state is degenerate and the phase is possibly gapless. However, the bipartite entanglement entropy S(x) is close to constant throughout most of the ladder which would suggest a gapped phase (see Appendix C).

Total magnetization, m ⊥ and m
It is useful to analyze the total magnetization of the open chain by separating the components parallel, m , and perpendicular, m ⊥ , to the [111] field direction. We define s α = i s α i /N and then It follows that To facilitate visualization it is most convenient to plot |m ⊥ | with a sign that we determine as sign(m ⊥ ·â) withâ = (1, 1, −2)/ √ 6. Surprisingly, m ⊥ is completely aligned or anti-aligned withâ for the two cases we shall now discuss and the angle m ⊥ forms withâ in theâ,b plane is either 0 or π. Our results are shown in Fig. 7(a),(b) for φ = 0 and φ = −0.08π, respectively. The results are from high precision iDMRG calculations with a unit cell of 60 and a maximal bond dimension of 1,000. An important point to notice in Fig. 7 is that m ⊥ is 2-3 orders of magnitude smaller than m .
Let us first discuss the results shown in Fig. 7(a) obtained from a field sweep from 0 to 0.9 at φ = 0. As the η phase is entered the iDMRG calculations fail to converge and that part of the plot is therefore colored lightly red. Starting from zero field m (blue circles) is found to approximately linearly increase with h[111] until the SC phase is reached where a kink in m is observed, consistent with a divergence of χ e h . Through the SC phase m increase more rapidly, this phase is therefore 'softer', consistent with the large chain end excitations shown in Fig. 6(a),(b). A second kink in m is observed as the γ phase is entered but the increase in m through out the γ phase is less pronounced. As the η phase is entered and excited kinks in m are again observed and in the PS phase the m tend toward the fully polarized value of 1/2. On the other hand m ⊥ (red circles) is non monotonic throughout the AK phase, approaches small values in the SC phase before jumping to larger values in the γ phase. In the polarized state (PS) phase m ⊥ changes sign and approach zero from below.
Our results for a similar field sweep at φ = −0.08π is shown in Fig. 7(b). The position of this field sweep is indicated as the vertical dotted line in Fig. 4(b). Again, m is seen to increase more rapidly through the SC and β phases as compared to the AK phase. The phase transitions between the AK SC and β phases are clearly visible as kinks in m and as the PS phase is approached m approach 1/2 in a characteristic cusp that was not observed between the η and PS phase. The nature of the β-PS transition and η-PS transition therefore clearly appear different. In this case m ⊥ is featureless at the AK to SC transition but the SC−β and β−PS transitions are clearly visible. m ⊥ changes sign in the SC and β phases. However, m ⊥ is in this case 3 orders of magnitude smaller than m .
For both field sweeps at φ = 0 and φ = −0.08π we emphasize that m ⊥ is either aligned or anti-aligned withâ.

Field and angle sweeps
To further investigate the sharpness of the phase transitions occurring in Fig. 4 Fig. 9. In panels (a), (b) and (c) we show results for χ e h , χ λ1 h and − ln λ 1 , respectively. Five transitions are clearly visible. While the PS-UC transition at φ = −0.289π is easy to miss in χ e h it is very well defined in χ λ1 h and − ln λ 1 . However, there is a precursor peak in χ e h not associate with a transition (see also Appendix A). The transition SC-AK is clearly defined at φ = −0.132π as the µ phase is approached from the SC side it is also clear that χ e h , χ λ1 h diverge at the µ-SC transition at φ = −0.146π. This transition is therefore well defined. On the other hand, the iDMRG fail to achieve good convergence in the light red colored region so the transition from SN to µ is not clear. As discussed above it is possible that the µ phase only marks the onset of short range incommensurate correlations in the SN phase and it is not a distinct phase. It is also possible that improved convergence of the iDMRG would show a divergence in χ e h and χ λ1 h inside the lightly red colored region.

B. Scalar Chirality in the UC and SC phases
The presence of a non-zero Γ term or the magnetic field raises the possibility of chiral ordering. The chirality without magnetic ordering is rare, unless there are three or four-spin interactions. In 1D system, it was shown that a four-spin interaction produces a long-range scalar chirality. 69 To check the presence of chiral ordering, we label the i'th spin on twolegs of the ladder as S i,1 and S i,2 where 1 and 2 refer to the bottom-leg and top-leg, respectively. We then define the scalar chiral order parameter with S = σ/2 as follows.
This clockwise definition is kept for all triangles made of three spins, for example κ = σ i,2 ·(σ i+1,2 × σ i+1,1 ) for upper triangles. If the κ is positive (negative), we assign blue (red) arrows i → j → k for κ = σ i · (σ j × σ k ) and all even permu- tations of i, j, k, which leads to the clockwise (anti-clockwise) circulation. It is also of interest to define the scalar chiral correlation function, such that C κ (r) → κ 2 . The scalar chiral order parameter breaks spatial symmetries and time reversal symmetry, but not SU(2). In Fig. 10(a) we show κ along a field sweep at constant φ = −0.27π through the UC phase. Results are from high precision iDMRG with a unit cell of 60 sites and a maximal bond dimension of 1,000. In this case we determine κ from C κ (r) shown in the inset Fig. 10(d) at h[111] = 0.125. C κ (r) reaches a constant value relatively quickly and κ is clearly non-zero throughout the phase, reaching significant values at the center of the UC phase. The scalar chirality, κ is negative throughout the chain and a sketch of the spatial modulation is shown in Fig. 12. Note the edge states appearing on the two opposing legs. Here weaker lines indicate a weaker κ. Since the h[111] field favors alignment of the spins, which would result in κ = 0, it is rather surprising to observe such a well defined scalar chirality at finite fields. In Fig. 10(b), (c) are shown χ e h and χ λ1 h , respectively. While the transitions delimiting the UC phase are almost absent in χ e h , and clearly do not coincide with the broad maximum of χ e h around h [11] ∼ 0.25, they are very well defined in χ λ1 h and in both cases they coincide with the results for κ in Fig. 10(a). As far as we can tell the UC phase does not intersect the zero field axis, instead, as can be seen in Fig. 4 the UC and SN phases meet at a triple point close to φ = −0.265π. Entanglement close to this triple point is therefore very elevated which is why − ln λ 1 (blue squares in Fig. 10(c)) is so high close to zero field.
In Fig. 11 we show iDMRG calculations for κ versus h[111] at a fixed φ = −0.08, crossing the AK, SC and β phases. χ e h and χ λ1 h along the same line in the phase diagram are shown in Fig. 8. As before, κ is obtained from calculations of C κ (r) in the large r limit and the sign of κ from local direct estimates of κ. In this case there is a spatial +-alternation of κ as shown in Fig. 13 but the magnitude of κ is the same for each triangle leading to the zero flux in the system. Although there is weak chirality in the AK and β phases, κ is an order of magnitude larger in the SC phase and jumps rather abruptly at the critical points indicated by the dotted lines in Fig. 11 which attain a constant value for modest values of r ∼ 20. Panel (c) in Fig. 11 show the bipartite entanglement obtained from a bi-partition of the system at site r. The resulting entanglement entropy S(r) = −Trρ r ln ρ r is shown versus r for a finite ladder with N = 400 and OBC. Clearly, S(r) is close to constant in the middle of the ladder which would be consistent with the existence of a non-zero gap in the SC phase.
(See also Appendix C). To make a comparison to the AK phase, we also compute the chirality in the AK phase, and its pattern is shown in Fig. 14. κ has distinctly different circulations. It has a staggering circulation between upper and lower triangles, but different staggering from the SC phase. For both the AK and SC phases, the magnitude of κ is the same for all pairs of triangles, implying that there is zero net-flux in the system. The phase transition from the AK to SC phases is accompanied by the sharp change of both magnitude and distinct pattern of κ.

C. String Order and Mapping to KQ Model
In Ref. 62 a non-local unitary transformation, V , was introduced of the following form: We then define the following unitary (disentangling) operator for an N -site chain with OBC: With the individual U (j, k) given as follows: At back to the original Kitaev ladder one arrives at a string correlation function of the following form: r odd (15) Note that, in Ref. 62 some of the indices in Eq. (15) in the expression for r even were incorrect. With this definition we find that the usual plaquette operator 1 W p ≡ O z (r = 4). W p is often used to characterize the AK phase. See Appendix D. Results for O z are shown in Fig. 15. In the presence of a non-zero magnetic field or a non-zero Γ the string order correlation function O z (r) is not long-range. Instead, as shown in the inset in Fig. 15(a) for h[111] = 0.25, it decays exponentially to zero. However, the length scale describing this exponential decay is extremely large, often exceeding hundreds or for small enough Γ, h[111], thousands of lattice spacings. The extent of the AK phase can therefore be determined by determining O z (r = 100) or O z (r = 900) which remain non-zero through out the AK phase. This is illustrated in Fig. 15(a)  Clearly, O z (r = 100) drops abruptly towards zero at the phase transition between AK and the SC phase. Fig. 15(b) shows O z (r = 100) versus φ at fixed h[111] = 0.15 again abrupt changes in O z (r = 100) are observed at the boundary of the AK phase.

Mapping to KQ model
It is of considerable interest to explore non-local unitary operators that will lead to string order correlation functions showing long-range order also for Γ = 0. We begin by considering how the ladder is transformed under the U 6 transformation previously described. If the spins on one leg of the ladder are numbered i = 1 . . . N/2 we assign them a second label k = (i − 1)mod6 + 1 on the second leg we assign the label k = (i + 2)mod6 + 1. With this labelling we introduce the following notation for the transformed bonds Here j is a nearest neighbor site. With this notation we can represent the transformed ladder using the following picture, We then consider the zero field case with OBC and change notation slightly compared to Ref. 62 by using an equivalent non-local unitary operator: With the individual w(j, k) given as follows: and W † = W . With this definition of W we see that on the vertical bonds of the ladder W leaves all interactions unchanged However, on horizontal bonds we find W S y 2 S y 4 W = −S y 1 S y 4 , W S y 3 S y 5 W = −S y 4 S y 5 , . . . (20) Note that W effectively moves the bond and changes the sign of the interaction. Likewise, we get for the horizontal zz bonds However, the horizontal xx bonds give rise to non-trivial 4spin interactions. Specifically: thereby coupling the four spins around a plaquette. We then introduce additional notation for transformed bonds With this notation in hand we can now apply the W transformation to the U 6 transformed ladder shown in Fig. 16. The resulting Hamiltonian can be drawn in the manner shown in Fig. 17. It is quite remarkable that the W transformation has generated four-spin exchange terms. We call this H KQ Hamiltonian the KQ ladder since the model is a close cousin of JQ models 57 extensively studied as models of deconfined criticality 58 . Unexpectedly, the ground-state for the KQ model in the AK phase has a significant overlap with a rung-triplet state. If we define: Then we can pictorially draw the rung-triplet state as shown in Fig. 18. With OBC there is another energetically equivalent rung-triplet state obtained by translation as shown in Fig. 18. The two states approximates the doublet ground-state of the KQ model. Surprisingly, the unit cell for these triplet states are 12 sites and not as one might have expected from the structure of H KG , 6 sites. Ordering of this type has previously been studied using SOP's inspired by the studies of s = 1 spin chains. If τ α i = S α i,1 + S α i+1,2 are the sum of two diagonally situated spins, one defines 70,71 : This order parameter has been used to distinguish topologically distinct phases in Heisenberg ladders and is non-zero in the rung-singlet phase where the topological number is even. In Fig. 19 we show results for O z even estimated from O z even (r) in the large r limit on systems with N = 120 and OBC. This SOP drops abruptly to zero at the limits of the AK phase and, as can be seen from the inset, there is no exponential decay observed in O z even (r) which instead quickly attains a constant value. As expected, O z even is also non-zero in the RS U6 phase, but clearly zero in the SN phase. Also shown in Fig. 19 are the overlaps with the rung-triplet states A|ψ and B|ψ as obtained from small N = 12 systems with OBC. There are significant finite-size effects at the boundaries of the AK phase for the overlaps, although they drop to zero outside the AK phase the critical points are not as well defined as for O z even . The rung-triplet states shown in Fig. 19 are approximate and for N > 12 different linear combinations enter such that with the simple definitions of the rung-triplet states above A|ψ and B|ψ tend to zero as N → ∞. However, O z even is clearly non-zero in the AK phase and in the KQ model this phase is therefore topologically equivalent to the Haldane like phases observed in Heisenberg ladders with even topological number. Since the KQ model is related to the KG model through unitary transformations the same must be true for the KG ladder throughout the AK phase. Calculations in the FK phase show that O z even is also clearly non-zero in zero field throughout that phase but drops to zero in the AΓ phase.

VI. COMPARISON: THE HONEYCOMB KG MODEL
It is important to note that the pure AFM Kitaev under the magnetic field studied by DMRG and 24-site ED exhibits an intermediate gapless phase before it polarizes. A U(1) spin liquid was suggested for this field-induced gapless phase 45 . However, in the two-leg ladder, we found five distinct phases including AK SC, γ, η, and PS at the pure AFM Kitaev point (white thin line at φ = 0 in Fig. 4), as the magnetic field increases. There are several factors including the obvious geometry difference that may result in the different results between the 24-site honeycomb cluster and current results. To understand possible origins of the difference, we investigate the following systems at AFM Kitaev point, φ = 0.
First we study the 24-site ladder using ED and compare the result with iDMRG with a unit cell of 60, to understand the finite size effects. The results of χ e h are shown in Fig.  19 (a) where the red and green dots are obtained by iDMRG and ED, respectively. Note that, due to problems with convergence there are no iDMRG results in the η phase. There are a couple of sharp features within η phase only found in ED, which we assign to finite size effects. Other than these additional peaks in η phase, the results are remarkably similar.
We also investigate the 24-site C 3 symmetric honeycomb cluster using the ED. χ e h is shown in Fig. 19 (b), where h sweeps from 0 to 0.8 by δh = 0.001, much smaller steps than previous studies. 45 45 only find 2 transitions performing ED on the same 24-site C 3 cluster, presumably due to a larger δh. However, based on the finite size effects found in 24-site ladder ED, we suspect that there are significant finite size effects in this system that makes it hard to determine whether these phases correspond to the SC, γ and η phases, or only one (SC) or two phases survive in the thermodynamic limit. Indeed when the field is tilted away from the c-axis, there are only two transitions found in the 24-site C 3 cluster 45 , while the three intermediate phases -SC, γ and η -found in the iDMRG ladder persist even in the titlted field (see Appendix E). Given that γ and η are incommensurate magnetic phases, and thus sensitive to the shape of cluster, we speculate that they likely turn into another type of incommensurate phase confined in the C 3 symmetric cluster. The chirality of the SC phase may survive in the honeycomb cluster defined at a triangle made of nearest neighbors or next nearest neighbors of honeycomb lattice. This is an excellent topic for future study, as it requires a bigger size system to check the chiral correlation C(r).
To shed further light on the occurrence of gapless excitations in 24-site C 3 cluster in ED, we choose an even smaller cluster of 2 × 6 ladder, and compute various quantities under the [111] field. The energy spectrum is shown in ably similar to 24-site C 3 cluster, suggesting that the intermediate gapless excitation feature is insensitive to the cluster size and shape, even though the critical fields where such phases arise change depending on the shape. The qualitative behaviour is similar to the 24-site ladder and C 3 ED presented above. However, it is not clear if there is one or two intermediate phases, as the incommensurability of γ and η phases would suffer from the change of cluster size and shape, as discussed above. The dynamical structure factor at a particular momentum Q 2 (shown in the inset) and the specific heat in the field are shown in Fig. 22 (a) and (b), respectively. The first intermediate phase is likely disordered, while the second intermediate phase likely exhibits incommensurate ordering. They do not exhibit well-defined excitation spectra in the specific heat similar to what is observed in the C 3 24-site ED. 45 While it is a 12-sites ladder geometry, the qualitative results are incredibly similar to the DMRG phase diagram at the φ = 0 AFM Kitaev point under the field and the 24-site C 3 symmetric honeycomb geometry results. 45 It shows a dense energy spectrum in the intermediate states of both disordered (SC) and field-induced incommensurate (η or γ) phase.
The above analysis with small ladder and C 3 cluster suggests that the AFM Kitaev 2D honeycomb model under the field may also display a richer phase diagram than what has been reported, and a high resolution numerical calculation is required to refine the phase diagram. Since the γ or η phases also show a dense energy spectrum, it is important to differentiate the U(1) spin liquid from the incommensurate phases in the honeycomb AFM Kitaev limit. Comparing the critical field above which gapless spin liquid occurs in 24-site ED 45 , it is likely that the disordered SC phase with enhanced entanglement and edge excitations extends to the gapless spin liquid in the honeycomb lattice, and the incommensurate phase is mixed with a polarized state, which was missed in 24-site ED of the honeycomb lattice. We conclude that the ladder model at the AFM Kitaev limit captures both disordered and incommensurate magnetically ordered phases under the magnetic field, and offers future directions in searching for a spin liquid in the honeycomb KG model.

VII. SUMMARY AND DISCUSSION
The KG model consists of two bond-dependent interactions, namely the Kitaev and Gamma interactions. The Kitaev interaction on the honeycomb lattice exhibits a spin liquid with fractionalized excitations. In particular, under a timereversal symmetry breaking term, the excitations obey non-Abelian statistics. The Gamma interaction is another highly frustrated interaction leading to a macroscopic degeneracy in the classical limit, and quantum fluctuations do not lift the degeneracy found in the AFM classical Gamma model. 72 Since these two frustrating interactions are dominant interactions in realistic descriptions of emerging Kitaev candidates such as RuCl 3 , the minimal KG model was initially proposed to understand RuCl 3 . 49,73 The magnetic field has been a crucial parameter, as the system may undergo a transition into a field-induced disordered phase before the trivial PS appears. Aside from its relevance to Kitaev materials, the minimal KG model may offer a playground to discover exotic spin liquids due to the combined frustration of the K and Gamma term, and thus has been extensively studied for the last few years. Given the huge phase space of AFM and FM Kitaev, AFM and FM Gamma, and the field, most studies are limited to a narrow phase space focusing on the FM Kitaev and AFM Kitaev regions. Many numerical methods have been used to identify phases of the extended Kitaev model under the field and intriguing results were reported near AFM Kitaev and FM Kitaev regions including the field-induced gapless U(1) spin liquid near the AFM Kitaev region. 45 However, it is not clear if the gapless excitations are due to more conventional physics such as incommensurate ordering.
Here we investigate the entire phase space of the KG ladder model under the magnetic field. While the geometry is limited to the ladder, it has the great advantage of allowing for high numerical precision such as accessing iDMRG with a high precision mode with a unit cell of 60 and a maximal bond dimension of 1000. Numerical calculations are therefore very well controlled. We found an extremely rich phase diagram of the KG model under the field. Among fifteen distinct phases identified, nine phases appear near the AFM Kitaev region alone. In the zero field, there is a quadrupole ordered phase named SN, two magnetically ordered phases (FM U6 and RS U6 ) straightforward to understand from the mapping of 6-site transformation, and the disordered AK phase. It is interesting that the SN phase found in the KG chain 63 survives in the ladder. Other than the AK phase (which becomes the Kitaev spin liquid in the 2D limit), the zero-field phases are ordered and the entanglement entropy is rather low. Under the field, highly entangled phases emerge. Apart from several incommensurate magnetic ordered phases, two highly entangled phases denoted by SC referring staggered chirality and UC uniform chirality are induced by the field. These phases exhibit distinct chirality orderings and high entanglement entropy with gapless edge excitations when the boundary is open.
The ladder results presented here offer several important insights in possible spin liquids and not-yet-identified phases in the 2D honeycomb lattice. We would like to recall that the pure Kitaev model in the ladder corresponding to AK and FK in this study is gapped, where the ladder can be viewed as a coupled chains. 74 As the number of chains grows, the ground state changes between gapped and gapless depending on the even and odd numbers of the chains, and eventually maps to the 2D Kitaev spin liquid in a true 2D limit. The AK and FK are gapped due to the geometry of the ladder, but its nature, magnetically disordered with high entanglement, is captured in the ladder model. Applying similar logic, we suggest that the disordered SC phase is related to the spin liquid in the 2D limit. The SC phase has a staggered chirality but different patterns from the AK phase in field, which differentiate the two phases. While it is gapped in the ladder, it may become gapless as the number of chains grows.
Interestingly the previous DMRG studies on the pure AFM Kitaev point (φ = 0) under the [111]-field with three to five number of legs reported different central charge in the intermediate phase. For three-leg chains, the central charge c = 1 50 was reported, while for four-and five-leg chains, c = 0 50 and c ∼ 4 39 respectively were found. Based on the central charge arguments, these studies indicate that there are gapless excitations associated with a spinon Fermi surface in the intermediate field region. It was suggested that the spinon Fermi surface pockets are around K/K -and Γ-point of the first Brillouin zone 50 , while the other DMRG study proposed the pockets around M -and Γ-points. 52 . The existence of a spinon Fermi surface in momentum space is yet to be determined.
While the most previous studies reported the field-induced intermediate phase as a single phase at the pure AFM Kitaev φ = 0 point in the [111]-field 37,40,41,45,50,52 , a separation of the intermediate region into three phases was noted in 39 , and the middle phase, which corresponds to the γ phase in the ladder, grows in extent with larger bond dimensions where five-legs were used. We find three different intermediate phases, SC, γ and η in the ladder at φ = 0. While the γ phase occupies a tiny phase space in the ladder, it is possible that this gapless incommensurate γ extends its phase space, as the number of chain grows. This implies that the chiral spin liquid candidate SC may require a finite FM Γ interaction, as it generates more frustration. The SC phase with long-range chirality and enhanced entanglement appearing at the intermediate field region of the ladder likely evolves to a field-induced spin liquid with a finite staggered chirality.
The UC phase is another candidate of spin liquid. It has uniform chirality pattern with high entanglement with a finite net flux, and it appears at very low field between SN and 6site transformed FM phase. This phase space has not been well explored in honeycomb clusters, and we suggest further studies in this region to look for a possible spin liquid. The nature of the SC and UC phases and statistics of excitations in these phases are excellent topics for future study.
Possible incommensurate orderings in the 2D limit also deserves some discussion. Incommensurate orderings are generally difficult to pin down, as they depend on the size of cluster, and the ordering wave-vector itself changes even inside a phase due to the nature of incommensuration. In the ladder, we found several incommensurate orderings. Some has high entanglement indicating a quantum order coexisting with an incommensurate ordering. The possibility of incommensurate ordering has been excluded in the C 3 symmetric cluster, mainly because of technical difficulty set by a limited size. Our ladder study suggests several incommensurate orderings may be present in the 2D honeycomb lattice, and future studies on such a possibility on larger honeycomb clusters are desirable. yellow, however, this feature is not associated with any real phase transition.
Another useful depiction of the AFM Kitaev region phase diagram can be obtained from S rung . With ρ N/2−1 the reduced density matrix of the first N/2 − 1 sites of the ladder, the bipartite entanglement entropy is obtained as S rung = −Trρ N/2−1 ln ρ N/2−1 , with the partition cutting the middle rung (and both legs) of the ladder. Our results for this quantity are shown in Fig. 23(b). In this case the UC phase along with all the other phases are clearly defined. Regions of increased values of S rung are visible as bright yellow colored bands in the µ and τ phases. These bands likely describe lock-ins to particular magnetic orderings compatible with the unit cell used in the calculations. The 'heart of entanglement', with its elevated entanglement, encompassing the µ, SC and γ phases is beautifully illuminated in yellow and orange colors. Even higher bipartite entanglement is present in the η phase with its bright yellow colors. χ e θ and (b) ln χ λ1 θ versus θ at a fixed φ = −0.21π. This is to check if there is a first order transition as a function of J occuring along the line of the SN. Indeed a clear first order transition is found in the energy suscepbility at θ = π/2, i.e. J = 0 as shown in Fig. 24. We conclude that the phase space denoted by the SN is a line of first order transitions separating two different magnetic ordering states. Since the nematic (spin-quadropole) order is finite in both ordered states, we keep the name the SN for this line of first order transitions, throughout the paper.

Appendix C: Bipartite entanglement in the AFM Kitaev region
If we instead of defining the reduced density matrix at x = N/2 − 1, as we did when considering S rung , but instead at a general x always taken to be odd, we obtain the bipartite rung entanglement entropy, S rung (x) = −Trρ x ln ρ x . Our results for S rung (x) at various points in the AFM Kitaev region are shown in Fig. 25 and 26 corresponding to the same points as the ones shown in Fig. 5, 6 in the main text. It is well established that the entanglement entropy is bounded 75,76 in systems with a gap. For a system with a gap S rung (x) should then attain a plateau in the middle of the chain for long enough ladders when that limit is attained. On the other hand, for a 1D critical gapless system with open boundary conditions we expect 77 S(x) = c 6 ln 2N π sin πx N + ln g + s 1 /2.
Here c is the central charge, g the universal ground-state degeneracy 78 and S 1 a non-universal constant. As is well known, the entanglement therefore grows logarithmically and is not bounded. The fact that a plateau is observed in the middle of the ladder for S rung is therefore consistent with a gap. From the results shown in Fig. 25 and 26 the behavior of S rung (x) in the SC, AK, UC and γ phases are therefore indicative of a gapped phase. The results for the remaining phases are more consistent with a gapless phase, although without The plaquette operator, W p , introduced by Kitaev 1 , is a crucial tool for characterizing the Kitaev spin liquid. On the honeycomb lattice W p is defined as a product of the six spins surrounding a hexagon in the lattice. In the setting of the KG ladder derived from a strip of the two-dimensional honeycomb lattice W p then corresponds to the 6 spins around 2 plaquettes of the ladder. Following the numbering of the sites on the ladder shown in Fig. 16 in the main text we then have: Our results for W p are shown in Fig. 27. At both the AFK and FM Kitaev points, φ = 0 and φ = π we have a fluxfree state with W p = 1 on all plaquettes. Away from these points W p can be different from 1. However, as can be seen in Fig. 27, W p remains high throughout the AK phase only decreasing when the phase is exited. However, the different phases in the AFK region are only partly visible in Fig. 27. Interestingly, W p is negative in the β and η phases shown as blue in Fig. 27.