Comprehensive study of the phase diagram of the spin-1/2 Kitaev-Heisenberg-Gamma chain

A central question on Kitaev materials is the effects of additional couplings on the Kitaev model which is proposed to be a candidate for realizing topological quantum computations. However, two spatial dimension typically suffers the difficulty of lacking controllable approaches. In this work, using a combination of powerful analytical and numerical methods available in one dimension, we perform a comprehensive study on the phase diagram of a one-dimensional version of the spin-1/2 Kitaev-Heisenberg-Gamma model in its full parameter space. A strikingly rich phase diagram is found with eleven distinct phases. In addition to the established phases in the Kitaev-Gamma chain in previous works, we find eight additional phases when a nonzero Heisenberg term is added, including four Luttinger liquid phases, a ferromagnetic phase, a N\'eel ordered phase, an ordered phase of distorted-spiral spin alignments, and another ordered phase which breaks $D_3$ symmetry, where $D_3$ is the dihedral group of order six. Our work paves the way for studying one-dimensional Kitaev materials and may provide hints to the physics in higher dimensional situations.


I. INTRODUCTION
The Kitaev spin-1/2 model on the honeycomb lattice [1] with anisotropic bond-dependent Ising interactions is proposed to host exotic quasiparticle excitations, including Majorana fermions, and non-abelian anyons under applied magnetic fields. A remarkable feature of these excitations is that their braiding and fusion operations can be used to realize topological quantum computations [2]. For this reason, the model has stimulated intense research interest in the past decade [2][3][4].
It was first proposed that the Mott insulating A 2 IrO 4 (A=Li, Na) compounds on the honeycomb lattice provide a platform for realizing the Kitaev spin-1/2 model [5]. The d 5 -configuration of the magnetic Ir 4+ ion is subject to strong cubic crystal fields due to the surrounding octahedral environment formed by O 2− ions. The Ir 4+ ion is in its low spin configuration with all five electrons residing in the three t 2g orbitals, and the spin and effective orbital angular momenta are S = 1/2 and L = 1, respectively. Strong spin-orbit coupling plays a crucial role in transferring the oxygen-mediated orbital-dependent superexchange Hamiltonian to an anisotropic effective spin-1/2 model. Indeed, due to strong spin-orbit couplings, the ground state of a single ion is a Kramers doublet with a total angular momentum j = 1/2, and the projection of the superexchange model to the j = 1/2 subspace is the desired bond-dependent Kitaev spin-1/2 model [5].
Many theoretical efforts have been devoted to studying the phase diagram and fractional excitations of the Kitaev model on the honeycomb lattice augmented with Heisenberg and Gamma couplings [6,26,[31][32][33][34][35][36][37][38]. However, most of the studies are based on a classical analysis, mean field theories, or exact diagonalization (ED) on a small system, and a controllable understanding is usually lacking, which is a typical difficulty in two dimension (2D). On the other hand, in one dimension (1D), there are more reliable analytical and large-scale numerical methods to study the low energy properties, including bosonization [39,40], conformal field theory (CFT) [3,[41][42][43]45], and the density matrix renormalization group (DMRG) methods [46][47][48]. Hence, a detailed investigation of a 1D version of the Kitaev model may shed light on the physics in 2D. It also provides a starting point for an extrapolation to 2D by coupling the 1D chains together and tuning the interchain coupling strength from weak to strong. In addition, a study of the 1D generalized Kitaev model also has its own merit, since it can be realized in systems of Ruthenium stripes within a-or b-oriented superlattices of RuCl 3 [49].
The phase diagrams of the Kitaev-Heisenberg model in the 1D chain [50] and the two-leg ladder cases [51,52] have been analyzed before. In particular, the two-leg ladder model already has a phase diagram quite similar to the model on the 2D honeycomb lattice [31]. On the other hand, the Kitaev and Gamma couplings dominate over the Heisenberg coupling in α-RuCl 3 [26]. Thus a good approach is to first consider the Kitaev-Gamma model and then treat the Heisenberg term as a small perturbation. Along this line of logic, the phase diagram of the 1D spin-1/2 Kitaev-Gamma chain has been investigated in Refs. [53,54], which already shows a rich phase diagram. In addition to the emergent SU(2) 1 phase which covers the experimentally relevant parameter region of FM Kitaev and AFM Gamma couplings, the system also contains a rank-1 spin ordered phase with O h → D 4 symmetry breaking, a spin-nematic (i.e., spin-quadrupole) ordered FIG. 1: Sketch of the phase diagram for the spin-1/2 Kitaev-Heisenberg-Gamma chain with eleven phases in total, and cartoon plots for the long-range and quasi-long-ranger orders in the corresponding phases. The full parameter space is a two dimensional unit sphere parametrized by the polar and azimuthal angles θ and φ. Due to the equivalence (K, J, −Γ) (K, J, Γ), only the front half of the unit sphere corresponding to θ, φ ∈ [0, π] is shown. The (θ, φ) coordinates of the K, −K, Γ, AFM1, FM1, AFM2, FM2, AFM3, FM3 points are ( π 2 , 0), ( π 2 , π), ( π 2 , π 2 ), (0, #), (π, #), ( π 2 , π 4 ), ( π 2 , 3π 4 ), (π − arctan(2), 0) and (arctan (2), π), respectively, in which the symbol "#" is used when the value of φ can be arbitrarily chosen. We note that the "LL2" phase locates on the circular boundary of the half sphere; the "Kitaev", "Nematic" and "O h → D4" phases are restricted to the equator; and all the other phases are extended in a finite area in the phase diagram. In the cartoon plots, the z-direction in spin space is chosen to be pointing vertically up, except for the "Neel", "d-Spiral" and "FM" phases where the black dot in the upper-right corner indicates that the z-axis is chosen to be perpendicular to the plane. In the LLi (i = 1, 2, 3, 4) phases, the black arrows represent a site-dependent quantization axis for the direction of longitudinal fluctuations. In all cartoon plots, the spin directions refer to Eq. (1) without any sublattice rotation. phase with O h → D 3d symmetry breaking, and a Kitaev phase, where O h is the full octahedral group, D n is the dihedral group of order 2n, and D nd is D n × {e, i} where e and i are the identity element and the inversion operation, respectively. In the peculiar Kitaev phase, no numerical signatures of Luttinger liquid behaviors nor spin orderings are found [54]. Interestingly, it was conjectured in Ref. 54 that the phase transition between the rank-1 and spin-nematic ordered phases is possibly a deconfined quantum critical point.
In this work, we include a nonzero Heisenberg term, and make a comprehensive study of the phase diagram of the spin-1/2 Kitaev-Heisenberg-Gamma chain in its full parameter space, using a combination of extensive analytical, DMRG and ED calculations. Here we note that although the Kitaev coupling in the existing materials is FM, it is worth also to study the region with AFM Kitaev coupling, since there may be f -electron systems realizing AFM Kitaev materials [30]. Thus, our work provides a road-map to the exotic physics in both the existing and potential Kitaev materials. A strikingly rich phase dia-gram is found with a total of eleven distinct phases as shown in Fig. 1, where the Kitaev, Gamma, and Heisenberg couplings are parametrized as K = sin(θ) cos(φ), Γ = sin(θ) sin(φ), and J = cos(φ), respectively. The plot in Fig. 1 is only schematic, and the precise phase boundaries are shown in Fig. 3 as determined by DMRG numerics. Also, only half of the parameter space within the range θ, φ ∈ [0, π] is shown due to the equivalence (θ, φ) (θ, 2π − φ) (see Eq. (3)).
From an analytic point of view, there are six special points with explicit or hidden SU(2) symmetries, which provide starting points for a perturbative analysis in the regions nearby. More precisely, the AMF1 and FM1 points located at the north and south poles are explicitly SU(2) symmetric with AFM and FM couplings, respectively. On the other hand, the AFM2 and FM2 points have hidden SU(2) symmetry as revealed by a six-sublattice rotation [56], and the AFM3 and FM3 are SU(2) symmetric after a four-sublattice rotation [31,56].
For the regions close to the AFM points, we perform an RG analysis which is controllable as long as the devi-ations from the AFM points are not large. The low energy field theory is first obtained by directly projecting the perturbation Hamiltonian to the low energy SU(2) 1 degrees of freedom, and then subject to an RG analysis. We have also performed a careful symmetry analysis in each region to verify that the projected perturbation Hamiltonian contains all the symmetry allowed terms up to relevant and marginal couplings classified by the SU(2) 1 critical theory. In this way, five phases can be understood, including the Luttinger liquid phases "LL1", "LL2" and "LL3" ("LL" is "Luttinger liquid" for short), and the ordered phases "Néel" and "FM" with Néel and FM orders, respectively, which are all confirmed by our ED and DMRG numerical calculations. The Luttinger parameters in each Luttinger liquid phase are determined numerically, and a color plot of the results are shown in Fig. 3. We note that the above symmetry and RG analysis has a wider applicable range not limited to the spin-1/2 Kitaev-Heisenberg-Gamma model. The conclusions also hold when another off-diagonal symmetric Γ term [1] and terms involving non-nearest-neighbors are included as long as the lattice symmetries are respected and the perturbations are small; though of course, the coupling constants in the low energy theory will be renormalized.
There are three remaining phases which cannot be captured by the above RG analysis. We find a narrow "d-Spiral" phase close to the FM Kitaev point in the upper hemisphere as shown in Fig. 1. When Γ = 0, there is a duality transformation for the Kitaev-Heisenberg chain which maps the interval [−K, FM3] to [FM2, −K] on the circular boundary of Fig. 1. This establishes [FM3, −K] to have a spiral spin order as a consequence of the FM order in the arc [FM2, −K], and the symmetry breaking of the spiral order is determined to be (Z 2 × Z 2 ) D 4d → (Z 2 × Z 2 ) D 2 . Then the symmetry breaking pattern in the "d-Spiral" phase is inferred from the symmetry breaking of the spiral order in the arc [FM2, −K] based on the assumption that there is no phase transition between the "d-Spiral" and the spiral orders. The thus obtained symmetry breaking D 4d → D 2 in the "d-Spiral" phase predicts a "distorted-spiral" pattern of the spin alignments, which is supported by our DMRG numerics.
To understand the "LL4" and "D 3 -breaking" phases, we take a perturbative Luttinger liquid approach. The strategy is to separate the Hamiltonian in the sixsublattice rotated frame into two parts, such that one part is of the easy-plane XXZ type, and the other part is taken as a perturbation. Analysis shows that the first order effect of the perturbation Hamiltonian vanishes, hence the Luttinger liquid behavior can be stabilized in a window of intermediate values of J, which is identified as the "LL4" phase. DMRG numerics provide evidence for the existence of the "LL4" phase as shown in Fig. 3.
When J approaches zero, the Luttinger parameter diverges, hence higher order effects eventually become important and drive the system into a strong coupling limit. By a careful combination of strong coupling and symme-try analysis, two different types of symmetry breaking patterns are identified for small |J|, namely, D 3d → Z (I) 2 and D 3d → Z (II) 2 , where Z (I) 2 and Z (II) 2 are two different Z 2 groups. Since D 3d /Z 2 ∼ = D 3 , both phases break the D 3 symmetry, which is the origin of the name "D 3 -breaking" for the phase in Fig. 1. The type of D 3 -breaking order is selected by the sign of the coupling constant, the determination of which requires a third order perturbation calculation. We do not perform such a difficult high order perturbation calculation, but instead turn to a classical analysis. Based on the classical analysis, we find that the "D 3 -breaking" phase can be divided into two subregions corresponding to the "D 3 -breaking I" and the "D 3 -breaking II" phases separated by the dashed line in Fig. 1, which have Z (I) 2 and Z (II) 2 as the unbroken symmetry groups, respectively. However, our DMRG numerics provide evidence for the coexistence of the two orders in the entire region marked as "D 3 -breaking" in Fig. 1. Since the two D 3 -breaking orders have distinct symmetry breaking patterns, generically the two orders should not coexist, and any coexistence has to be accidental from a symmetry point of view. Whether the revealed coexistence is a truth or a finite size artifact remains unclear and is worth further studies.
The rest of the paper is organized as follows, where each section is made self-contained for the convenience of the readers who are interested in specific phases in the phase diagram. In Sec. II, the model Hamiltonian is introduced, and the sublattice rotations are discussed which reveal the hidden SU(2) symmetric points AFMi and FMi (i = 1, 2, 3) shown in Fig. 1. In Sec. III and Sec. IV, we combine RG calculations, symmetry analysis, and numerics together to study the "LL1" and the "Néel" phases, respectively. A brief description of the "LL2" phase is also included in Sec. IV. In Sec. V, the "d-Spiral" phase is studied. The symmetry breaking pattern is identified to be D 4d → D 2 , and the spin alignments are shown to exhibit a "distorted" spiral pattern. In Sec. VI, the "LL3" phase is investigated, again by a combination of RG, symmetry, and numerical analysis. Sec. VII is devoted to a discussion of the "LL4" and "D 3 -breaking I, II" phases. Numerics provide evidence for a coexistence of the proposed two types of D 3 -breaking orders. The origin for the coexistence remains unclear. In Sec. VIII, the "FM" phase is discussed. Finally in Sec. IX, we briefly summarize the main results and open questions of the paper.

A. The Hamiltonian
We consider a spin-1/2 Kitaev-Heisenberg-Gamma (KHΓ) chain [1] in zero magnetic field defined as in which i, j are two sites of nearest neighbors; γ = x, y is the spin direction associated with the γ bond shown in Fig. 2 (a); α = β are the two remaining spin directions other than γ; K, J and Γ, are the Kitaev, Heisenberg and Gamma couplings, respectively. The terms in H are spelled out explicitly in Supplementary Materials [55]. Throughout this work, we parametrize K, J, Γ as in which θ ∈ [0, π] and φ ∈ [0, 2π]. It is straightforward to observe that a global spin rotation R(ẑ, π) : ) leaves K and J invariant but changes the sign of Γ, in which R(n, θ) represents a rotation in spin space around then-direction by an angle θ. Hence, i.e., (θ, φ) (θ, 2π − φ). Due to this equivalence, the phase diagram will be studied within the parameter range θ ∈ (0, π), φ ∈ (0, π). It is apparent that the north and south poles where K, Γ vanish have explicit SU(2) symmetries: θ = 0 corresponds to the AFM Heisenberg model, and θ = π is the FM Heisenberg model. These two SU(2) symmetric points are denoted as AFM1 and FM1 in the phase diagram in Fig. 1. In this subsection, we briefly review the six-and foursublattice rotations [31,53,56,57], and show that they unveil several points in the parameter space which have hidden SU(2) symmetries. These hidden SU(2) symmetric points provide starting points for a perturbative study of the regions surrounding them.

The six-sublattice rotation
The six-sublattice rotation U 6 is defined as [53,57] Sublattice 1 : (x, y, z) → (x , y , z ), in which "Sublattice i" (1 ≤ i ≤ 6) represents all the sites i + 6n (n ∈ Z) in the chain, and we have abbreviated S α (S α ) as α (α ) for short (α = x, y, z). The Hamiltonian in the six-sublattice rotated frame is in which γ = x, z, y has a three-site periodicity shown in Fig. 2 (b), and S i is denoted as S i for simplicity. The terms in H are spelled out explicitly in Supplementary Materials [55]. It is clear from Eq. (5) that while the Kitaev and Gamma terms acquire a form similar to the Heisenberg model (but with unequal couplings along different spin directions) in the six-sublattice rotated frame, the Heisenberg J term loses its form. Indeed, H is SU(2) invariant when K = Γ, J = 0. Combining U 6 with Eq. (3), we see that (θ = π/2, φ = π/4) and (θ = π/2, φ = 3π/4) have hidden SU(2) symmetries with FM and AFM couplings, respectively, which are denoted as the FM2 and AFM2 points in the phase diagram shown in Fig. 1.
It is straightforward to observe that the Kitaev points are self-dual under U 4 . Setting Γ = 0 and normalizing the transformed parameters according to K 2 + J 2 + Γ 2 = 1 (where K = K + 2J, J = −J, Γ = 0), U 4 establishes the equivalences: in which [A, B] represents the arc between the points A and B on the circular boundary of Fig. 1.

C. Summary of the phase diagram
Here we make a summary of the phase diagram shown in Fig. 1. On the equator, there are four phases, namely, "Kitaev", "Nematic", "O h → D 4 ", and "Emergent and [φ c , π], respectively, in which φ c 0.034π, φ c 0.10π and φ c 0.33π as determined in Refs. 53, 54. We note that the "Emergent SU(2) 1 " line is part of the "LL1" phase. In Fig. 1, the spin orientations in the "O h → D 4 " phase is shown in the original frame, where the z-direction in spin space is chosen to be perpendicular with the plane. A cartoon plot for the spin-nematic order in the "Nematic" phase is also shown in Fig. 1, where the z-direction is taken to be vertical. As discussed in Ref. 54, the two largest spin-nematic order parameters in the original frame arê In Fig. 4, the orbital analogues of the adjacent-site spinquadrupoles are plotted forQ e . More precisely, the orbital analogue is z 2 − y 2 (z 2 − x 2 ) on the x-(y-) bond in Fig. 2 (a).
In this work, we focus on the other phases corresponding to J = 0. The "LL1" and "FM" phases are analyzed by combining RG and symmetry analysis together taking the emergent SU(2) 1 line on the equator as the unperturbed part of the Hamiltonian. These two phases are revealed to be a Luttinger liquid and an FM ordered phase, respectively. In a similar manner, the "Néel" and "LL3" phases can be understood through a combination of RG and symmetry analysis by perturbing the AFM1 and AFM3 points, respectively, and they correspond to a Néel ordered and a Luttinger liquid phase. We note that the "LL2" phase has been discussed in Ref. 50, which can also be obtained from an RG analysis.
The remaining phases cannot be approached by direct RG analysis and we will take different routes. For the "d-Spiral" phase, we identify the symmetry breaking pattern to be D 4d → D 2 based on the assumption that there is no phase transition between the "d-Spiral" phase and the known spiral phase on the arc [FM3, −K]. The corresponding spin ordering is shown to exhibit a distortedspiral pattern which is confirmed by our numerics. For the "LL4" and "D 3 -breaking I, II" phases, a unified perturbative Luttinger liquid analysis is performed to understand the two phases. We find that the Luttinger liquid behaviors can be stabilized in an intermediate range of J, which is confirmed by the DMRG numerics in the "LL4" phase. At small |J|, analytical analysis reveals two different D 3 -breaking orders, occupying different subregions in the "D 3 -breaking" phase. However, numerics provide evidence for the coexistence of these two different orders. Whether the coexistence is a truth or only a finite size artifact is worth further studies.
The Luttinger parameters have been calculated numerically based on the method described in Ref. 58. The results are shown as color plots in Fig. 3, in which the blue color is used to signify the situation where no reliable Luttinger parameter can be extracted. As can be seen from Fig. 3, it is a very nice result that nearly all the phase boundaries in Fig. 1 can be determined from this Luttinger parameter calculation. We note that the phase boundaries in Fig. 1 are only schematic, and their precise shapes should be referred to Fig. 3.
The cartoon plots for the phases with J = 0 are collected in Fig. 1, all referring to the original frame without any sublattice rotation. The patterns of the spin orientations in the ordered phases "Neel", "d-Spiral", The blue color is used to denote the situation where no reliable value of the Luttinger parameter can be extracted. The phase boundaries are plotted as solid white lines, except the line segment separating the "LL4" and the "FM" phases close to the equator which is represented by a dashed white line, since the precise shape of the phase boundary in that region cannot be clearly identified due to the percolation of the Luttinger liquid behaviors in the "LL1" phase into the lower hemisphere. The names of the phases are written in white letters, where "SN" is "Spin-Nematic" for short. DMRG numerics are performed on systems of L = 96 sites with open boundary conditions. "D 3 -breaking I, II" and "FM" are plotted, where the zdirection is chosen to be vertical for the "D 3 -breaking I, II" phases, and perpendicular to the plane for the other ordered phases. The quasi-long range orders in the Luttinger liquid phases "LLi" (i = 1, 2, 3, 4) are also shown, in which the z-axes are all fixed to be pointing upwards. The black arrows represent the site-dependent directions of the quantization axes for the longitudinal fluctuations, whereas the shaded blue ellipses represent the plane of the transverse fluctuations which dominate over the longitudinal fluctuations in all the four Luttinger liquid phases. We note that since the low energy Hamiltonian for the "LL4" phase in the six-sublattice rotated frame has an FM-type quasi-long range order, there is no oscillation accompanying the power decay in the cartoon plot of the "LL4" phase in Fig. 1. Finally, we make a comment about the ED and DMRG numerics that we have performed in this work. For systems with open boundary conditions, the DMRG method was used on chains with length up to L = 144 sites. For some of the calculations, such as the ground state energy computations determining the boundaries of the phases, we used ED on chains up to L = 24 sites long, while DMRG with periodic boundary conditions was used for chains of L = 36 sites. In all the cases, we have checked that our DMRG results are converged using up to m = 1000 states with a truncation error below 10 −7 .

III. THE "LL1" PHASE
In this section, we show that the region denoted by "LL1" in Fig. 1 is described by the gapless Luttinger liquid theory. The system exhibits a site-dependent quantization axis for the longitudinal fluctuations within the original frame as shown in Fig. 5. The strategy for analyzing the "LL1" phase is to take the "Emergent SU(2) 1 " phase of the Kitaev-Gamma chain [53] on the equator of Fig. 1 as the unperturbed Hamiltonian, and treat the Heisenberg term as a small perturbation using a perturbative RG analysis. To facilitate analysis, we work in the six-sublattice rotated frame defined by Eq. (4) throughout this section unless otherwise stated. The black arrows denote the directions of the quantization axes, and the solid blue ellipses represent the transverse fluctuations. The red line represents the AFM quasi-long range order for the longitudinal and transverse fluctuations defined in terms of the six-sublattice rotated frame. The z-direction in spin space is chosen to be pointing upwards, and the ydirection is along the chain pointing to the right. gent SU(2) 1 " phase when φ c ≤ φ ≤ π and θ = π/2, which provides an RG perturbative starting point for analyzing the "LL1" and "FM" phases on the two different sides of the equator. Due to the equivalence Γ −Γ, we will consider the equivalent region in the other half of the equator, i.e., θ = π/2, φ ∈ (π, 2π − φ c ). In particular, the point φ = 5π/4 has explicit SU(2) symmetry in the rotated frame.
The low energy degrees of freedom in the emergent SU(2) 1 phase are J L , J R , g, where J L , J R are the WZW left and right currents and g is the SU(2) 1 primary field which is a 2 × 2 matrix. At low energies, the lattice spin operators S i 's can be expressed in terms of J L , J R and g using the following modified nonabelian bosonization formula [53], in which E = C, D. The low energy Hamiltonian of the Kitaev-Gamma chain is [53] in which v is the spin velocity, and g c > 0 is the marginally irrelevant coupling. The values of both v and g c depend on the microscopic details, which in principle can be obtained from DMRG numerical calculations on the Kitaev-Gamma chain. We note that the spin-spin correlation functions can be calculated using Eq. (11) and Eq. (13).

B. Low energy Hamiltonian
In this subsection, we derive the low energy perturbation Hamiltonian by directly projecting the Heisenberg Hamiltonian to the low energy space. The low energy field theory is found to be the same as that of an XXZ chain with a quantization axis along the (1, 1, 1)direction. In practice, the spin operators within are replaced by J L , J R and g using the modified nonabelian bosonization formula in Eq. (11). This method is essentially a first order perturbation treatment of the Heisenberg term, and the coupling constants thus obtained are not accurate in the sense that they acquire renormalizations along the RG flow. On the other hand, by performing a careful symmetry analysis on the low energy theory, we have justified the use of this first order perturbation method in capturing the essential physics.
To obtain the perturbation Hamiltonian, the following operator product expansion (OPE) formula for the Néel order fields is needed, (15) in which x is a spatial coordinate, a is the lattice constant, λ, µ = x, y, z are spin directions, {...} and [...] in the superscripts denote symmetrization and antisymmetrization of the indices, respectively, and only terms up to quadratic order in the WZW currents are kept. Eq. (15) can be derived from the affine symmetry of the SU(2) 1 WZW model as discussed in detail in Supplementary Materials [55]. Plugging Eq. (11) into Eq. (14) and using Eq. (15), we obtain in which only the relevant and marginal terms are kept; the arrow "→" indicates that it is not an exact equality but only a projection; and the coefficient u 2 is The values of the other coefficients u i 's (i = 1, 3, 4, 5) are not important for our purpose and can be found in Supplementary Materials [55], in which a detailed derivation of Eq. (16) is also included. Furthermore, we have performed a careful symmetry analysis in the low energy Hamiltonian showing that Eq. (16) contains all the relevant and marginal terms allowed by symmetry within the SU(2) 1 WZW model (for details, see Supplementary Materials [55]). Therefore, the low energy Hamiltonian in Eq. (16) is enough and complete to capture the physics as long as |J| is small. In the SU(2) 1 WZW theory, J 0 λ J 0 λ (λ = L, R) is equal to 1 3 J λ · J λ (see Supplementary Materials [55]), hence it does not give an independent contribution. In addition, the inversion breaking term J 0 L −J 0 R can be eliminated by a chiral rotation [62][63][64]. As a result, the only nontrivial SU(2) breaking term in the low energy theory of the KHG chain is J 0 L J 0 R . This shows that the low energy physics is the same as that of an XXZ chain with a quantization axis along the (1, 1, 1)-direction. Whether the system remains gapless or develops an order is determined by the sign of the coupling − 1 3 Ju 2 . To gain a simple understanding as to why the (1, 1, 1)direction is special, here we give a brief description of the symmetry group in the six-sublattice rotated frame. With a nonzero Heisenberg term, the symmetry transformations of the Hamiltonian H in Eq. (5) are in which T is time reversal; T a is translation by one lattice site; I is the spatial inversion around the point C in Fig.  2 (b); and R a = R(n a , −2π/3), R I = R(n I , π) wherê In Eq. (19), the choice of the inversion center in the definition of I is well defined modulo three. If another inversion center 5 + 3n is chosen, then S α 10−i (α = x, y, z) in the second line in Eq. (19) should be correspondingly replaced by S α 10+6n−i . According to Eq. (19), the symmetry group G 1 of H is In the continuum limit, the lattice is coarse-grained and the effect of T a is smeared. Hence, it is not a surprise that the symmetry operation R a T a picks outn a (i.e., the (1, 1, 1)-direction) to be the quantization axis in the low energy theory (which is also justified by the symmetry analysis on the low energy Hamiltonian as discussed in Supplementary Materials [55]). Finally, we make a comment on the group structure of G 1 which will be used in Sec. VII. Since T 3a = (R a T a ) 3 belongs to G 1 , it is legitimate to consider the quotient group G 1 /<T 3a >. In our case, time reversal plays the role of "inversion" in the spin space, since T changes the signs of the spin operators. Neglecting the spatial components in R a T a and R I I, the spin operations R a and R I generate the symmetry group of a regular triangle as shown in Fig. 6. Since D 3d = D 3 × Z 2 in which D 3 is the symmetry group of a regular triangle and Z 2 is generated by the inversion operation [66], we see that <T, R a , R I > ∼ = D 3d . As proved in Supplementary Materials [55], <T, R a T a , R I I>/<T 3a > is isomorphic to D 3d even when the spatial components in the operations are included. As a result, G 1 /<T 3a > D 3d . This shows that the group structure of G 1 is in which 3Z = <T 3a >.
In a perturbative treatment of the Heisenberg term, we will assume that |J| va, g c a. With a nonzero Heisen-berg term, the low energy Hamiltonian is modified to (23) in which according to Eq. (16), In particular, since u 2 is always positive (as can be seen from Eq. (18)), we have y > y ⊥ when J > 0, and y < y ⊥ when J < 0. We note that in Eq. (23), the chiral term J 0 L − J 0 R term is dropped according to the discussions in Sec. III B. The RG flow equations of y ⊥ and y are of the Kosterlitz-Thouless (KT) type [65]: which are obtained by integrating out the modes with wavelengths between e l a and e l+dl a. It is well known that the system flows to a strong coupling limit when y < y ⊥ , and remains gapless when y > y ⊥ [65]. The phase diagram by tuning J is shown in Fig. 7. Notice that this phase diagram can be easily understood from intuitive physical arguments. The J > 0 case introduces an FM-like term into the Hamiltonian since the coupling is −J in the rotated frame. Thus the AFM fluctuations in the longitudinal direction are suppressed and the system has an easy-plane anisotropy, corresponding to the "LL" phase in Fig. 7, where "LL" is "Luttinger liquid" for short. Similarly, the J < 0 case has an easy-axis anisotropy which corresponds to the "Ordered" phase in Fig. 7. We note that Fig. 7 is for a fixed background Kitaev-Gamma chain, i.e., the value of φ in Fig. 1 is fixed. By varying φ, the gapless phase in Fig. 7 extends to the "LL1" phase in Fig. 1, whereas the "Ordered" phase in Fig. 7 is below the equator, which will be shortly shown to be the "FM" phase in Fig. 1. We also note that the 0-direction (i.e., the longitudinal direction) in Eq. (23) is the (111)-direction, which becomes staggered in the original frame according to Eq. (4). A sketch of the site-dependent quantization axis is shown in Fig. 5. However, the quantization axes in Fig. 5 are not precise since they can be distorted due to the renormalization effects of the spin operators along the RG flow similar as the origins of the C 1 , C 2 coefficients discussed in Eq. (11).
When J is negative, we have seen that the system develops a Néel long range order along the (111)-direction.
However, the alignments of the spins are distorted due to the difference in C 1 and C 2 . Using Eq. (11) and assuming a nonzero expectation value for N 0 = 1 √ 3 (N x +N y +N z ), the spin polarizations have a six-site periodicity, in which n ∈ Z and N = N 0 . Recall that these are the results within the six-sublattice rotated frame. To obtain the spin orientations in the original frame, it is enough to perform the inverse of the transformation defined in Eq. (4). It is straightforward to verify that where the superscript "(0)" is used to denote the spin operators in the original frame. Clearly, Eq. (27) represents an FM order along the (C 1 C 1 C 2 )-direction. This provides an understanding of the FM phase in Fig. 1 in the region close to the equator.

D. Numerical results
To determine the range of the "LL1" phase, we study the ground state energy E as a function of θ, φ. ED calculations are performed for periodic systems of L = 24 sites. In Fig. 8 (a), ∂ 2 E/∂φ 2 is scanned horizontally for several fixed values of θ, in which different curves are shifted vertically. The peaks correspond to divergent positions of ∂ 2 E/∂φ 2 , indicating first order phase transitions. We emphasize that only the two largest peaks in Fig. 8 (a) for each θ represent first order phase transitions, which determine the left and right boundary points of the black dome in Fig. 8 (c) at the corresponding cut of θ. The middle smaller peaks in Fig. 8 (a) do not scale with systems sizes as shown in Fig. 8 (b), therefore there is no phase transition at these positions in the thermodynamic limit. The absence of phase transition in the middle region is also confirmed by our numerical calculations of the Luttinger parameter which will be discussed in Sec. III D. Fig. 8 shows the phase boundary of the "LL1" phase determined from Fig. 8 (a), as represented by the curve connecting the small solid black squares. Starting from the "LL1" phase, the system transits into the "Néel" phase and the "d-spiral" phases by going upwards and rightwards, respectively, as shown in Fig. 8 (c). The "Néel" and "d-spiral" phases will be discussed in Sec. IV and Sec. V.
Next, we numerically determine the Luttinger parameter K in the "LL1" phase following the method described in Ref. 58. The low energy field theory is given by the following Luttinger liquid Hamiltonian in which u is the velocity and K is the Luttinger liquid parameter. While K = 1/2 corresponds to the SU(2) symmetric case, the symmetry is reduced to U(1) when K = 1/2. Define the Hamiltonian density h(x) as For a finite size system with an open boundary condition, the energy density h(x) contains a uniform part E U (x) and a staggered part in which x = ja (j 1) is the distance measured from the boundary of the system. Consider a general function f (j) (j ∈ Z) which contains a uniform part u(j) and a stagger part s(j), i.e., A three-point formula can be used to extract u(j) and s(j), as Hence E A can be obtained from the numerical result of h(x) by applying Eq. (32).
We have used this method to study the Luttinger parameters in the whole "LL1" phase, and the results are displayed in Fig. 3 where the magnitudes of K are represented by different colors. This provides direct numerical evidence for the entire region enclosed by the phase boundary determined in Fig. 8 to be a Luttinger liquid phase. On the other hand, as can be seen from Fig. 3, the Luttinger liquids percolate into the "FM" phase. However, we note that this is a finite size artifact. As discussed in Sec. VI A, the phase transition between the "LL1" and "FM" phases is second order. Since a gap E g ∼ e −const/ √ |J| opens exponentially slowly in the "FM" phase close to the transition line [3], the crossover system size L c (only above which an order can be observed) grows exponentially at small J, i.e., L c ∼ e const/ √ |J| . Thus the Luttinger liquid behavior can still be observed in an extended region in the "FM" phase in a finite size system.

IV. THE "NÉEL" PHASE
In this section, we study the "Néel" phase shown in Fig. 1. The spin alignments within the original frame are plotted in Fig. 9. The strategy that we take to analyze the "Néel" phase is to perform a perturbative RG analysis in the neighborhood of the AFM1 point located at the north pole in Fig. 1. Numerics provide evidence for the Néel ordering to hold in the entire region marked as "Néel" in Fig. 1. Throughout the section, we work in the original frame, and the discussions will be based on the Hamiltonian given in Eq. (1).
FIG. 9: Spin alignments in the Néel phase within the original frame. The z-direction in spin space is chosen to be perpendicular to the plane, and the x-direction is along the chain to the right.
FIG. 10: Site-dependent quantization axes for the longitudinal fluctuations in the "LL2" phase within the original frame. The black arrows denote the directions of the quantization axes, and the solid blue ellipses represent the transverse fluctuations. The red line represents the AFM quasi-long range order for the longitudinal and transverse fluctuations. The z-direction in spin space is chosen to be pointing upwards, and the y-direction is along the chain.

A. RG analysis
We can perform a first order perturbation treatment of the Kitaev and Gamma terms, by projecting H K and H Γ to the low energy degrees of freedom using the nonabelian bosonization formula [3], where c is a constant. Comparing with Eq. (11), the SU(2) symmetry is not broken in Eq. (G22). Explicit calculations show that (for details, see Supplementary Materials [55]), and We have performed a careful symmetry analysis of the low energy field theories for both the Γ = 0 and Γ = 0 cases as discussed in details in Supplementary Materials [55], where it is shown that the first order perturbation Hamiltonian is already enough and complete to capture the low energy physics. Notice that the low energy Hamiltonian is of an XXZ type in the absence of Γ.
when Γ = 0, the system has easy-plane and easy-axis anisotropies for K > 0 and K < 0, respectively. The Luttinger liquid phase when K > 0 (and Γ = 0) revealed by this analysis is the "LL2" phase shown in Fig. 1. The "LL2" phase has already been discussed in Ref. [50], and our DMRG numerics have confirmed the existence of the "LL2" phase as shown in Fig. 3. We note that the percolation of the Luttinger liquid behavior in the "LL2" phase to nonzero φ's is a finite size artifact. The quantization axis for the longitudinal fluctuation in the "LL2" phase is along the z-axis as shown in Fig.  10. However, similar to the case of the "LL1" phase, the quantization axes in Fig. 10 are not precise since they can be distorted due to the renormalization effects of the spin operators along the RG flow.
Next we perform an RG analysis of the low energy Hamiltonian in Eq. (35) for nonzero Γ. We will see that it belongs to the XYZ class as long as Γ = 0, and an Néel order develops. The low energy perturbation Hamiltonian can be written as The matrix kernel in Eq. (36) can be straightforwardly diagonalized. The eigenvalues are with the corresponding unnormalized eigenvectors (38) in which k = 2+4c 2 3c 2 K Γ . Notice that as long as Γ > 0, we always have E + > 0 and E − < 0 regardless of the value of K.
Next, let's perform a rotation of the coordinate system, such that Then in terms of the new coordinates, the marginal terms become (40) in which λ x = (g c + E − )/(2πv), λ y = (g c + E + )/(2πv), and λ z = g c /(2πv). In particular, when K = 0, we have in which Γ is assumed to be small, and only terms up to O(Γ) are kept. The superscript "(0)" in Eq. (41) is used to indicate that these are the initial values of the couplings before the RG flow. In what follows, we consider the case Γ > 0, K = 0 for simplicity. The case of a nonzero K can be discussed similarly.
The RG flow equations of λ x , λ y , λ z are [3] which can be obtained from the OPE formula for the WZW current operators [3]. There are three constants of motion for Eq. (42) which take the following forms when K = 0: Solving Eq. (42) with initial conditions given in Eq. (41) (see Supplementary Materials [55] for derivation), we obtain where Λ → ∞ as l → ∞. It is clear that the system flows to strong couplings at low energies.
To identify the phase corresponding to the strong coupling limit, we perform a chiral rotation R L (ŷ, π), which maps (J Then it is clear that the three couplings after the chiral rotation becomeλ x ,λ y ,λ z −Λ → −∞. On the other hand, as shown in Ref. 3, this strong coupling limit corresponds to a dimer order, i.e., whereg = U L g is the SU(2) WZW field after the chiral rotation, in which U L = e i 1 2 σ y π = iσ y is the chiral rotation matrix. Therefore, Eq. (45) implies that i.e., there is a Néel order along the y -direction. In terms of the original coordinates, the Néel ordering is alonĝ 2 ) T , as can be seen from Eq. (39). Finally, we discuss the symmetry breaking of the Néel phase. The only broken symmetry in the "Néel" phase is the time reversal symmetry. The unbroken symmetry group H N can be determined as It is straightforward to verify that the most general pattern of the spin orderings which is invariant under H N is given by The values of the parameters a, b in Eq. (48) cannot be determined from pure symmetry analysis, and in general depend on K and Γ. A plot of the spin alignments in Eq. (48) is shown in Fig. 9.

B. Numerical results
Numerics provide evidence for the "Néel" phase to occupy the entire region above the equator in Fig. 1 excluding the "LL1" and the "d-Spiral" phases. Fig. 11 shows s αα (r) (α = x, y, z) vs. sin(πr/L) extracted from the three-point formula in Eq. (32) at a representative point (θ = 0.36π, φ = 0.25π) in the Néel phase, where (−) r s αα (r) is the staggered component of the spin-spin correlation function S α 1 S α 1+r . The DMRG numerical results for three different system sizes L = 48, 96, 144 are displayed with black, red and green curves, respectively. As can be seen clearly from Fig. 11, all the three s αα (r) (α = x, y, z) approach constant values when r 1, indicating a Néel long range order. In addition, the patterns of the spin orderings are fully consistent with Eq. (48), where a 2 , b 2 can be determined from the asymptotic values of s αα (r) at large r.

V. THE "d-Spiral" PHASE
In this section, we study the "d-Spiral" phase in Fig.  1. The spin orientations within the original frame are shown in Fig. 12. The symmetry breaking pattern in the "d-Spiral" phase is inferred from the symmetry breaking in the spiral phase [50] of the Kitaev-Heisenberg chain based on the assumption that there is no phase transition between the "d-Spiral" and the spiral phases. The spin orientations in the "d-Spiral" phase are predicted to exhibit a "distorted-spiral" pattern, which is supported by our DMRG numerics. Throughout this section, we work in the four-sublattice rotated frame defined in Eq. (6) unless otherwise stated. When Γ = 0, the system reduces to a Kitaev-Heisenberg chain even in the four-sublattice rotated frame by virtue of the duality property as mentioned in Eq. (8). Hence, in this case, there is no difference between the original and four-lattice rotated frames if the symmetry group structure is concerned. T in which the inversion center for I is taken to be the site 0.
We note that all the other symmetry transformations can be generated by the operations in Eq. (49) as discussed in Supplementary Materials [55]. Thus, the symmetry group of the Kitaev-Heisenberg chain is Next we briefly describe the group structure of G 0 (for a detailed proof, see Supplementary Materials [55]).
T a ] 4 belongs to G 0 , it is legitimate to consider the quotient group G 0 /<T 4a >, whose group structure is (see Supplementary Materials [55]) in which from left to right, This shows that in which 4Z = <T 4a >.
FIG. 13: R(ŷ, π) and R(ẑ, −π/2) as symmetry operations of a square. Within the xy-plane, the effect of the rotation R(ŷ, π) on the square is the same as that of a reflection with respect to the x = 0 plane. In the above figure, the z-direction is perpendicular to the plane.
We make some comments on the geometrical meaning for the origin of D 4d . Notice that D 4d = D 4 × Z 2 in which D 4 is the symmetry group of a square containing four reflections and four rotations, and Z 2 = <i> where i is the inversion operation. In our case, time reversal acts as "inversion" in spin space, since it flips the signs of the spin operators. Neglecting the spatial components in R(ẑ, − π 2 )T a and R(ŷ, π)T a I, the spin operations R(ẑ, − π 2 ) and R(ŷ, π) generate the symmetry group of a square as shown in Fig. 13. As proved in Supplementary Materials [55], <T, R(ẑ, − π 2 )T a , R(ŷ, π)T a I>/<T 4a > is isomorphic to D 4d even when the spatial components in the operations are included.

Symmetry breaking pattern for Γ = 0
In this subsection, we discuss the spin ordering and the symmetry breaking pattern for the Γ = 0 case, i.e., the Kitaev-Heisenberg chain, and later extend the analysis to Γ = 0 in Sec. V B 2. The spin ordering for θ ∈ [−K, F M 3] located on the circular boundary in Fig.  1 has been worked out in Ref. 50 which exhibits a spiral pattern with a four-site periodicity. On the other hand, this region of parameter can be mapped to the FM-xy phase in Fig. 1 by the four-sublattice rotation. Therefore, the system is FM aligned along ±(110)-direction in the four-sublattice rotated frame. In what follows, we determine the symmetry breaking pattern corresponding to this FM spin order within the four-sublattice rotated frame.
Let's consider the little group of the FM-xy order, in which the spin alignments are S i ≡ (f, f, 0) T . Apparently, the spin alignments are invariant under both T a I and T 2a . Furthermore, although R(ẑ, − π 2 )T a and R(ŷ, π)T a I do not leave the FM-xy order invariant, the spin orientations are invariant under the combinations It is easy to verify that c and d satisfy c 2 = d 2 = (cd) 2 = e modulo T 4a . Comparing with the generator-relation representation for the group D n [66] we see that relations in Eq. (E8) are are satisfied for c, d with n = 2, hence the group generated by c, d (modulo T 4a ) is isomorphic to D 2 . Therefore, the little group of the FM-xy order is Combining with the expression for G 0 in Eq. (E13), we conclude that the symmetry breaking pattern is in which T 4a has been dropped for simplicity. In particular, since |D 4d /D 2 | = 4, the ground states are fourfold degenerate in which the spins align ferromagnetically along (±1, ±1, 0)-directions. On the other hand, by transforming back to the original frame, it can be seen that the spin orientations show a spiral pattern as plotted in Fig. 12 (a).
B. Symmetry breaking for Γ = 0 1. Symmetry group for Γ = 0 When Γ is nonzero, the symmetry operations of H in the four-sublattice rotated frame in Eq. (7) can be verified to be the following,

1.
T in which the inversion center is taken to be site 2 in Fig.  2 (c). Therefore, the symmetry group is As discussed in Sec. V A 1, the group structure of G 3 is in which D 4d is given by Eq. (52), and 4Z = <T 4a >. Comparing Eq. (60) with Eq. (E13), G 3 is explicitly a subgroup of G 0 , and there are two more generators in G 0 than in G 3 , i.e., T a I and T 2a .

Symmetry breaking pattern for Γ = 0
In this subsection, we discuss the symmetry breaking and the spin ordering for the Γ = 0 case. Let's assume that T 4a is not broken, and we'll consider the quotient group G 3 = G 3 /<T 4a > in what follows. If there is no phase transition when Γ is tuned from nonzero to zero, then the symmetry breaking for a nonzero Γ has to be as a natural extension of the pattern in Eq. (57) for the Γ = 0 case. Assuming the symmetry breaking to be D 4d → D 2 , next we solve the most general spin ordering which is invariant under D 2 . Requiring the invariance of the spin ordering under the two generators (R(ẑ, − π 2 )T a ) 2 · T and R(ẑ, − π 2 )T a · R(ŷ, π)T a I, the most general D 2 -invariant spin configuration can be determined as  Rotating back to the original frame, the spin ordering is Eq. (63) exhibits a "distorted-spiral' pattern, which is the origin of the name "d-Spiral" for the phase where "d" is "distorted" for short. A schematic plot of the spin alignments in Eq. (63) is shown in Fig. 12 (b). Taking f = k, h = 0, the spin orientations in Eq. (63) reduce back to the Γ = 0 case in Fig. 12 (a). We note that the other three degenerate ground states can be obtained by applying the operations in the equivalent classes in D 4d /D 2 .

C. Numerical results
We first determine the range of the "d-Spiral" phase by studying the ground state energy E(θ, φ) as a function of θ, φ. Fig. 14 (a) shows the second order derivative To further confirm the nature of the peaks as first order phase transitions, we have studied the finite size scaling behaviors of ∂ 2 E(θ, φ)/∂θ 2 . Fig. 14 (b) displays the numerical results for ∂ 2 E(θ, φ)/∂θ 2 at φ = 0.95π with three different system sizes L = 12, 24, 36. As is clear from Fig.  14 (b), ∂ 2 E(θ, φ)/∂θ 2 scales linearly with L, indicating a first order phase transition. On the other hand, we find no divergence in ∂ 2 E(θ, φ)/∂θ 2 when the lower boundary of the "d-Spiral" phase is traversed. Recall that there is a divergence in ∂ 2 E(θ, φ)/∂φ 2 as a function of φ as discussed in Fig. 8 (c), which determines the lower boundary between the "d-Spiral" phase and the "LL1" phase shown by the solid black squares in Fig. 14 (c) within the region φ ∈ [0.96π, π].
Next we discuss the numerical evidence for the spin ordering in Eq. (62) within the four-sublattice rotated frame. Before going to that, we mention a subtlety in numerical calculations. The four symmetry breaking ground states only become exactly degenerate in the thermodynamic limit. The ground state of a finite size system can be some arbitrary linear combination of the four states, and the coefficients depend on the system size and numerical details. Because of this, random cancellations occur if the expectation values of the spin operators S i are directly computed. To circumvent this problem, we measure S i in the presence of a small field along (110)direction, which is able to polarize the system to reside in one of the four nearly degenerate states so that the spins orient according to the pattern in Eq. (63). The value of this field should satisfy ∆E |h|L E g , where ∆E is the finite size energy splitting among the four states, E g is the gap between the ground state multiplet and the excitations, and |h| is the magnitude of the applied field. Such choice of field leads to a degenerate perturbation within the ground state quartet, but no mixing is induced between the ground state subspace and the excitations. Fig. 15 shows the measured values of S α j (α = x, y, z) at a representative point (θ = 0.4π, φ = 0.99π) in the "d-Spiral" phase. The data are obtained by performing DMRG numerics on a periodic system of L = 36 sites under an h = 10 −4 field along the (110)-direction. As can be seen from Fig. 15, the pattern of the spin orientations are fully consistent with Eq. (62), with extracted values f, k, h as f −0.0936, k −0.0928879, and h 0.0127. This provides evidence for the existence of the "d-Spiral" phase. We note that numerics have also been done for several other points in the "d-Spiral" phase which are all consistent with Eq. (62).

VI. THE "LL3" PHASE
In this section, we study the "LL3" phase in Fig. 1, and demonstrate that the low energy physics in this phase is described by the Luttinger liquid theory. The system exhibits a site-dependent quantization axis for the longitudinal fluctuations within the original frame as shown in Fig. 16. Recall that as discussed in Sec. II B 2, the four-sublattice rotation reveals a hidden SU(2) symmetric AFM point, i.e., the AFM3 point in Fig. 1. We will study the region close to the AFM3 point using again a combination of RG and symmetry analysis. Numerics provide evidence for Luttinger liquid behaviors in the entire region of the "LL3" phase as shown in Fig. 3. To facilitate analysis, we work in the four-sublattice rotated frame defined by Eq. (6) throughout this section unless otherwise stated.

A. Perturbative analysis
We first remark that in the absence of Γ, there is a duality transformation for the Kitaev-Heisenberg chain as   Fig. 1. Since the latter (i.e., the "LL2" phase) has been shown to be described by Luttinger liquid theory in Sec. IV A, we conclude that [K, AFM3] is also in a Luttinger liquid phase [50].
Next we include a nonzero Γ. The Kitaev and Gamma terms will be treated as perturbations to the SU(2) 1 WZW theory which describes the low energy physics of the AFM3 point in the four-sublattice rotated frame. We perform a first order projection of the perturbing Hamiltonian to obtain the low energy field theory of the system. According to Eq. (7), the perturbing Hamiltonian ∆H is Notice that all of the operators J L , J R , g are invariant under translation by two sites, whereas the Γ term in Eq. (64) is staggered every two sites. Therefore, the Γ term vanishes after projecting to the the low energy degrees of freedom and thus has no effect. This means that at least up to first order projection, the low energy Hamiltonian of the KHΓ chain is of an XXZ type, having the same coupling constants as an Kitaev-Heisenberg chain. Based on this analysis, the phase diagram for fixed values of φ but changing θ in the vicinity of the AFM3 point can be derived as shown in Fig. 17. When θ > θ c , the system flows to a strong coupling limit where an order develops, whereas when θ < θ c , the system remains gapless. Up to first order approximation, θ c is equal to π − arctan(2) independent of φ. However, a concern at this point is whether higher order effects will spoil the XXZ-type low energy Hamiltonian, and in particular, whether the gapless Luttinger liquid phase in Fig. 17 is stable with respect to high order perturbations. To resolve this question, we have performed a careful symmetry analysis. The symmetry group of the system in the four-sublattice rotated frame has been demonstrated to be G 3 ∼ = D 4d 4Z in Sec. V B. We are able to show that up to relevant and marginal couplings, the symmetry allowed terms other than the XXZtype low energy Hamiltonian only include the chiral term J z L − J z R . Since this chiral term can be eliminated by a chiral rotation [62][63][64], the low energy physics remains to be of the XXZ type. The detailed symmetry analysis is included in Supplementary Materials [55]. This justifies the validity of the first order perturbation analysis and the phase diagram in Fig. 17. However, we emphasize that the value of θ c can be shifted due to high order effects, and the actual value of θ c is a function of φ. In the "LL3" phase, based on the above analysis, the quantization axis for the longitudinal fluctuations in the four-sublattice rotated frame is along the zdirection. Notice that by performing the inverse of the four-sublattice rotation defined in Eq. (6), the zdirection becomes staggered in the original frame, hence the quantization axes in the "LL3" phase are drawn in a staggered manner in Fig. 16. On the other hand, we emphasize that the quantization axes in Fig. 16 are not precise, since they can acquire site-dependent distortions due to the renormalization effects of the spin operators along the RG flow similar to the situation in Sec. III C.
We also note that the ordered phase in Fig. 17 is an Néel order along z-direction, which is in the foursublattice rotation frame. As can be checked by performing the inverse of the four-sublattice rotation, the Neel-z order becomes an FM-z order in the original frame. But this is not accurate due to the possible distortions of the quantization axes. To identify the precise direction of the ordering, we first figure out the unbroken symmetries of the Néel-z order in the four-sublattice rotated frame, and then solve the most general form of the spin orientations which are invariant under the unbroken symmetry group.
As can be easily checked, the unbroken symmetry group H 3 corresponding to the Neel-z order is Since T 4a = (T R(ẑ, − π 2 )T a ) 4 is an unbroken symmetry operation, we will consider spins in a unit cell in what follows where all site indices should be understood as modulo four. The most general spin alignments which are invariant under the group H 3 in Eq. (65) can then be solved as in which a, b are two parameters depending on the values of K, Γ, J. Performing the inverse of the four-sublattice rotation defined in Eq. (6), the spin ordering in the original frame can be determined as which is an FM order along the (a, a, b)-direction, providing another understanding to the "FM" phase in Fig.  1.

B. Numerical results
The Luttinger parameter K can be extracted numerically using the same method as Sec. III D. DMRG numerics are performed to calculate the energy density on an open system of L = 96 sites. We have studied the Luttinger parameter in the entire "LL3" phase and the results are shown in Fig. 3. On the right hand side of the phase boundary between the "LL3" and "D 3 -breaking" phases, no Luttinger parameter can be extracted, and in fact, the transition line in Fig. 1 between these two phases is determined in this way. On the other hand, as discussed in Sec. VI A, the phase transition between the "LL3" and "FM" phases is second order. Since a gap E g ∼ e −const/ √ |θ−θc| opens exponentially slowly in the "FM" phase close to the transition line [3], the Luttinger liquid behaviors percolate into the "FM" phase in a finite size system, which smears the phase transition as can be seen from Fig. 3.

VII. THE "LL4" AND "D3-BREAKING I, II" PHASES
In this section, we study the "LL4" and the "D 3breaking I, II" phases in Fig. 1. In the "LL4" phase, the system exhibits a site-dependent quantization axis for the longitudinal fluctuations within the rotated frame as shown in Fig. 18. In the "'D 3 -breaking I, II" phases, the spin orientations have six-site periodicities within the original frame as shown in Fig. 19. To facilitate analysis, we work in the six-sublattice rotated frame throughout this section unless otherwise stated.
The strategy is to separate the Hamiltonian into two parts, where one part is of the easy-plane XXZ type, and the other part is taken as a perturbation. We demonstrate that the first order projection of the perturbation term to the low energy degrees of freedom remains to be of the XXZ type, indicating that there exists a region of Luttinger liquid which is stable against the perturbation.
However, the Luttinger parameter diverges as J → 0, and higher order terms eventually become relevant driving the system into an ordered phase. As discussed in Sec. III B, the symmetry group in the six-sublattice rotated frame is isomorphic to D 3d 3Z. We demonstrate that there are two types of orders both having six-fold degenerate ground states. The symmetry breaking patterns are D 3d → Z are two different Z 2 groups. Since D 3d /Z 2 ∼ = D 3 , both ordered phases break D 3 symmetry (albeit different D 3 groups), and the two phases are named as "D 3 -breaking I" and "D 3 -breaking II". The regions occupied by these two different D 3 -breaking orders are determined by a classical analysis discussed in Sec. VII B 2. On the other hand, when |J| becomes large, the system is driven into a Néel ordered phase along the (111)-direction in the sixsublattice rotated frame, which corresponds to the "FM" phase in the original frame in Fig. 1.
Finally in Sec. VII C, we present and discuss DMRG numerical results, which seem to suggest the coexistence of the two D 3 -breaking orders. "D3 breaking II" phases within the original frame. The zdirection in spin space is chosen to be pointing upwards, and the y-direction is along the chain to the right.
A. Perturbative Luttinger liquid analysis and the "LL4" phase We will perform a perturbative analysis in a neighborhood of the FM2 point in the lower hemisphere. Thus Γ, K > 0 and J < 0. Performing the following twosublattice rotation to the six-sublattice rotated spin operators S α j 's, the Hamiltonian in Eq. (5) becomes in which where Notice that the S α j 's in ∆H (2) can be written in terms of S α j 's via Eq. (68). We note that the purpose of this additional two-sublattice rotation is to make the transverse directions to have an AFM coupling in H XXZ .
Clearly, H XXZ corresponds to an easy-plane XXZ chain when J < 0, Γ > 0. The low energy theory is described by the Luttinger liquid Hamiltonian in Eq. (28), and the Luttinger parameter K is known to be [65] in which J = Γ − |J| and J ⊥ = Γ. Next, we add ∆H and ∆H (2) as perturbations. Let's first consider ∆H (2) . In the Luttinger liquid phase, the continuum limit can be taken at low energies. To avoid the staggered sign, we enlarge the unit cell to six sites, and consider the following sum in which The first order effect of ∆H (2) can be obtained by projecting ∆H (2) to the low energy degrees of freedom using the following abelian bosonization formula for the spin operators [3], in which a is the lattice constant, n is the lattice site number, and x = na. We note that the first order perturbation projection can be figured out without even carrying out any calculation. Since different sites are smeared out in the continuum limit, both ∆H (2),a n and ∆H (2),b n become SU(2) symmetric ∼ S(x) · S(x + a). As a result, the effect of ∆H (2) is to shift the values of J and J ⊥ , as Since we still have |J | < |J ⊥ | when ∆H (2) is included, the system remains to be gapless at least for |∆| 1 where ∆ = (K − Γ)/Γ. Of course, the Luttinger parameter is renormalized due to ∆H (2) .
Next we consider the effects of ∆H. Again, the first order effect of ∆H can be obtained by projecting ∆H to the low energy degrees of freedom using Eq. (76). In fact, it can be observed that the projection vanishes without even carrying out any calculation. From Eq. (71), we obtain ∆H 1+6n,2+6n + ∆H 3+6n,4+6n + ∆H 5+6n,6+6n = in which In the continuum limit, Eq. (78) consists of total derivatives as can be seen by Taylor expanding F and G, which vanish in the Hamiltonian after the integration dx. Similar analysis can be performed on the other three terms ∆H 2+6n,3+6n +∆H 4+6n,5+6n +∆H 6+6n,7+6n . Hence we conclude that the projection of ∆H onto the Luttinger liquid degrees of freedom vanishes, and the system should remain in the Luttinger liquid phase up to first order in J.
However, the above analysis of first order projection cannot be trusted if higher order effects are included. Since K diverges as J → 0, higher order terms may become relevant for small enough J. Denote J α and N α as the uniform and staggered parts of the spin operators S α i , respectively, which, according to Eq. (76), are defined as and Since the scaling dimension of N ± (= N x ± iN y ) is 1/(4K), the operators involving powers of N ± become relevant for sufficiently small J. On the other hand, since N z has scaling dimension equal to K, the operators involving N z become relevant when K is small enough.
Hence, in what follows, we will consider terms only involving N α (α = x, y, z), which should be able to capture the spin orderings for both small and large values of J.
Next we perform a symmetry analysis to figure out what terms are allowed in the low energy theory. In the rotated basis defined via Eq. (68), the generators of the symmetry group G 1 in Eq. (21) become in which U 2 = Π n R 2n (ẑ , π)R which gives the transformation in Eq. (68). It is straightforward to verify that R a = R(ẑ , 2π 3 ) and R I = R(ŷ , π). Under the symmetries in Eq. (82), the transformation properties of N α are , and in which τ α (α = 1, 2, 3) are the three Pauli matrices, and τ 0 is the 2 × 2 identity matrix. Therefore, we see that the following terms are allowed in the low energy Hamiltonian, where n > 0, n ∈ Z. At large K, the symmetry allowed operator with the smallest scaling dimension (= 9/K) is which is ∼ J 3 cos(6 √ πθ) after bosonization. Notice that this term can only appear at the level of third order perturbations, hence the coupling should be proportional to J 3 . When K > 4.5, cos(6 √ πθ) becomes relevant. According to Eq. (73), the value of J which has the critical K value (= 4.5) is When ∆ = 0 (i.e., φ = 0.25π), the corresponding θ c1 determined from Eq. (87) is 0.514π. Of course, the value of J c1 (∆) in Eq. (87) is not accurate due to high order effects. Based on the above analysis, we see that the Luttinger liquid is stable only when |J| > |J c1 (∆)| for a fixed value of ∆. On the other hand, the above analysis applies only to the region of a small |J|. By increasing |J|, we note that the term (N z ) 2 ∼ cos(4 √ πφ) eventually has a smaller dimension than cos(6 √ πθ). When the Luttinger parameter K is smaller than 1/2, (N z ) 2 becomes relevant since its scaling dimension is 4K. Hence there exists another critical value |J c2 (∆)| above which the Luttinger liquid become unstable again. In conclusion, a Luttinger liquid can be stablized in an intermediate range B. The "D3-breaking" phase 1. The strong coupling limit for |J| < |Jc1| We will figure out the spin orders corresponding to the above mentioned two strong coupling limits when |J| < |J c1 | and |J| > |J c2 |. Let's first consider |J| < |J c1 (∆)|. The term λ cos(6 √ πθ) is minimized at the following values of the angle θ, in which λ is the coupling constant in the term λ cos(6 √ πθ), and 0 ≤ n ≤ 5, n ∈ Z. Notice that the system is six-fold degenerate regardless of the sign of λ. Naively, according to Eq. (76), the spins exhibit a Néel ordering in the x y -plane given by in which N ⊥ is the magnitude of the ordering. However, we note that the spin orders in Eq. (89) are not precise. The is because there can be U(1) breaking coefficients in the abelian bosonization formula Eq. (76) due to the renormalization effects in the high energy region along the RG flow, which is similar to the case encountered in Eq. (11) as discussed in details in Ref. 53. As a result, the abelian bosonization formula Eq. (76) should only respect the discrete lattice symmetry, not the emergent U(1) symmetry at low energies. Taking this into account, the true spin orientations will be distorted with respect to those in Eq. (89) due to the U(1) breaking coefficients. To figure out the correct pattern of the spin orientations, we first determine the symmetry breaking pattern corresponding to the spin ordering in Eq. (89), and then solve the most general form of the spin orderings which is consistent with the identified symmetry breaking. This approach is similar as the one used in Sec. VI A to discuss the "FM" phase based on the RG analysis in the vicinity of the AFM3 point.
When λ > 0, we choose a representative spin configuration corresponding to n = 1 in θ (I) n . As can be seen from Eq. (89), this is a Néel order along y -direction. As can be easily checked, the Néel-y configuration is invariant under T T 3a and R I I. Hence the little group of the Néel-y order is which can be rewritten as where Z (I) 2 = <R I I>, and 3Z = <T T 3a >. Recall that the symmetry group of the Hamiltonian H in Eq. (69) is G 1 = <T, R a T a , R I I>, which can be rewritten as where D 3d = <T, R a T a , R I I> mod T T 3a , and 3Z = <T T 3a >. Therefore, by taking the quotient of <T T 3a >, we conclude that the symmetry breaking pattern is Notice that since |D 3d /Z (I) 2 | = 6, the ground states are six-fold degenerate. This order is named as "D 3 -breaking I" since the broken symmetry group D 3d /Z (I) 2 is isomorphic to D 3 .
Next, we work out the most general form of the spin ordering which is invariant under H (I) . The invariance under T T 3a clearly requires FIG. 20: "Center of mass" directions of the three spins within a unit cell for the six degenerate ground states as represented by the six red (blue) solid circles in the "D3-breaking I (II)" phases in the six-sublattice rotated frame. In the "D3breaking I" phase, the red circles are connected with thin dashed red lines to indicate that the hexagon is coplanar. In the "D3-breaking II" phase, the plots are for J → 0 according to the classical analysis.
As can be checked, further requiring the invariance under R I I, the spin alignments are constrained to be in which x 2 + y 2 + z 2 = 1. The other five degenerate spin configurations can be obtained by performing the operations in the equivalent classes in D 3d /Z (I) 2 to the spin ordering in Eq. (95).
It is interesting to work out the spin alignments in the six-sublattice rotated frame by performing the inverse of the transformation defined in Eq. (68). The result is in which x, y, z can be expressed through x , y , z using Eq.
in which the superscript (0) is used to denote the original frame. A plot of the spin orientations in Eq. (97) is shown in Fig. 19 (a) where the unnormalized parameters are chosen as x = −1, y = 0, z = 1. A similar analysis can be performed to the case λ < 0. Choosing n = 0 in θ (I) n , we see that the system has a Néel ordering along the x -direction. It can be shown that the little group of the Néel-x configuration is This time, the symmetry breaking pattern is in which Z in which a 2 + b 2 + c 2 = 1, m 2 + n 2 = 1. We note that the ground states are again six-fold degenerate, and the spin orientations in the other five degenerate states can be obtained by performing the operations in the equivalent classes in D 3d /Z (II) 2 to the above spin configuration. This order is named as "D 3 -breaking II" since the broken symmetry group D 3d /Z (II) 2 is isomorphic to D 3 . We emphasize that although the broken symmetry is again isomorphic to D 3 , it is a different D 3 group compared with the "D 3 -breaking I" case.
Similar with the "D 3 -breaking I" case, the spin alignments in the "D 3 -breaking II" phase in the six-sublattice rotated frame can be determined as in which a, b, c, m, n can be expressed through a , b , c , m , n using Eq. (68). By performing the classical analysis discussed in Sec.VII B 2, we find that a, b, c, m, n approach 1, −1, 1, −1, 1, respectively, in the limit J → 0 for fixed ∆ = 0. Therefore, when J → 0, the "center of mass" direction of the three spins in a unit cell defined as S c 's in the six degenerate ground states when J → 0 are represented as the six blue circles as shown in Fig. 20. We also note that the corresponding spin orientations in the original frame are A plot of the spin orientations in Eq. (102) is shown in Fig. 19 (b) where the unnormalized parameters are chosen as a = −b = c = m = −n = 1. Bases on the above analysis, we see that when θ > θ c1 (φ), D 3 -breaking orders are developed, where there are two types of possible orders denoted as "D 3 -breaking I" and "D 3 -breaking II" depending on the sign of the coupling constants. The determination of the sign of the coupling constant requires a third order perturbation. We will not perform such a difficult calculation of third order perturbation, but turn to a classical analysis in the D 3 -breaking phase which will be discussed in Sec. VII B 2, where it is verified that the region for θ > θ c1 (φ) can be divided into two subregions numbered by "I" and "II", which have "D 3 -breaking I" and "D 3 -breaking II" orders correspondingly. The dashed line separating the two D 3 -breaking phases in Fig. 1 is determined from such classical analysis. We note that the line is plotted as dashed since it is not numerically verified which will be discussed in Sec. VII C.
We also note that the coupling constant ∼ J 3 is very small when θ > θ c1 (φ). Hence, from an RG point of view, the system has to flow a very long time to develop such order. This means that a very large system size is required to observe the D 3 -breaking orders, making the detection of the orders very difficult numerically. This might be the reason for why numerics detect a coexistence of the two D 3 -breaking orders on a finite size system as will be discussed in Sec. VII C.

Classical analysis
The classical analysis is the saddle point approximation in the spin path integral formalism which is valid in the large-S limit, where S is the spin value. In what follows, we neglect the quantum fluctuations of the spins and approximate them as classical three-vectors, in which S is the spin magnitude,n i = (x i , y i , z i ) T is a unit vector. Absorbing S 2 into a redefinition of Γ, K, J by defining and introducing the Lagrange multipliers {λ i } 1≤i≤3 to impose the constraints x 2 i + y 2 i + z 2 i = 1, the energy per unit cell in the six-sublattice rotated frame becomes in which the free energy corresponding to a general spin-S KHΓ chain is considered. We have numerically studied the minimization of the classical free energy, and verified that the solution of Eq. (95) (Eq. (100)) corresponds to the global minimum of the free energy in the range denoted by the "D 3 -breaking I" ("D 3 -breaking II") phase in Fig. 21.
3. The strong coupling limit for |J| > |Jc2| The above analysis completes the discussion for the strong coupling limit in the region |J| < |J c1 (∆)|. Next we identify the order in the parameter region |J| > |J c2 (∆)|. When cos(4 √ πφ) becomes relevant, φ n orders at (2n + 1) √ π/4 or n √ π/2 depending on the sign of the coupling constant. According to the bosonization formula in Eq. (76), the system develops a Néel ordering along z -direction. Performing the inverse of the transformation in Eq. (68), it corresponds to a Néel-n a order in the six-sublattice rotated frame, wheren a is the unit vector along the (111)-direction. As discussed in Sec. III C, by taking into account the distortions of the quantization axes in the bosonization formula and performing the inverse of the six-sublattice rotation, the Néel-n a ordering corresponds to an FM spin order in the original frame. Thus we conclude that the system should transit from the "LL4" phase to the "FM" phase by increasing the magnitude of J where J < 0. We can also make a tree-level estimation on the value of J c2 . As can be seen from Eq. (70), the anisotropy of the H XXZ + ∆H (2) Hamiltonian in Eq. (70) becomes easy-axis when J < −2Γ. Neglecting the effects of ∆H in Eq. (69), the critical value J c2 (∆) is determined to be J c2 (∆) ≡ −2Γ. At ∆ = 0, this gives the point (θ c2 = 0.804π, φ = 0.25π) in the phase diagram. Of course, the value of J c2 (∆) must be shifted due to the effects of ∆H and higher order effects of ∆H (2) .

Phase diagram around the FM2 point
Based on the above analytic analysis, we propose the following phase diagram which applies at least in a neighborhood of the FM2 point: "D 3 -breaking I, II", θ > θ c1 (φ); "LL4", θ c2 (φ) < θ < θ c1 (φ); "FM", θ < θ c2 (φ), where we have expressed J and ∆ in terms of θ and φ, and the ranges of θ all refer to the corresponding value of φ. In the "LL4" phase, the quantization axis for the longitudinal fluctuation is along the (111)-direction in the six-sublattice rotated frame which becomes staggered in the original frame with possible site-dependent distortions as discussed in Sec. III C. However, the staggered sign in the definition of the coordinates in Eq. (68) indicates an FM-type quasi-long range order in the sixsublattice rotated frame. Hence, no oscillation is drawn in the cartoon plot of the "LL4" phase in Fig. 18.

C. Numerical results
The numerical studies on the Luttinger parameter reveal that the "LL4" phase is rather narrow as shown in Fig. 3. While the phase transition line between the "LL4" and "FM" phases can be clearly identified at the far end of the "LL4 peninsula", it cannot be accurately determined close to the equator since the Luttinger liquid behaviors in the "LL1" phase percolate into the "FM" phase in finite size systems as discussed in Sec. III D. Thus the segment of the phase transition line between the "LL4" and "FM" phases in the vicinity of the equator is plotted as a dashed rather than solid line in Fig. 3. Here we note an interesting observation. According to the discussion in Sec. VII A, a rough estimation of the range of the "LL4" phase at φ = 0.25π is 0.51π < θ < 0.80π. However, as can be seen from Fig. 3, the actual range greatly shrinks compared with the above estimation. An explanation of why high order terms have such a huge effect is worth further considerations.
Next we numerically study the spin ordering in the "D 3 -breaking" phase in Fig. 1. We emphasize that the six-sublattice rotated frame is taken, not the frame after the further transformation defined in Eq. (68). As discussed in Sec. V C, a small field has to be applied to test the spin orders. We consider two types of fields: h I along (−1, 0, 1)-direction, and h II along (1, −1, 1)-direction. If the system is in the "D 3 -breaking I" phase, then the field h I is able to polarize the system such that the spins are aligned according to the pattern in Eq. (96); on the other hand, if the system is in the "D 3 -breaking II" phase, then an h II -field will polarize the spins into the pattern in Eq. (101). Fig. 22 shows the numerically measured expectation values of S α j (α = x, y, z) under h I -, h II -fields (both equal to 10 −4 ) at a representative point (θ = 0.51π, φ = 0.15π) within the "D 3 -breaking" phase in Fig. 1. ED numerics are performed on L=24 system with periodic boundary conditions. As can be seen from Fig. 22, the patterns of S α j are consistent with Eq. (96) (Eq. (101)) under the h I -(h II -) field shown as the black (red) data points. The magnitudes of both spin orders are huge (about 10 3 times larger than the applied fields), indicating a coexistence of the two types of classical orders. We have tested other points in the "D 3 -breaking" phase and found that the coexistence occurs in the entire phase. However, since the two types of classical orders correspond to different symmetry breaking patterns, generically the two orders should not coexist. Whether this is a finite size artifact or there is indeed a coexistence remains unclear and is worth further studies.

VIII. THE "FM" PHASE
In Sec. III C, Sec. VI A and Sec. VII B 3, we have inferred the "FM" phase based on RG analysis in three different regions, i.e. in the region close to the "Emergent SU(2) 1 " line on the equator of Fig. 1, the region close to the "LL3" phase, and the region close to the "LL4" phase. The spin alignments in the "FM" phase are shown to be

IX. CONCLUSIONS
In summary, we have studied the phase diagram of the spin-1/2 Kitaev-Heisenberg-Gamma chain. In addition to the already established phases for the Kitaev-Gamma chain, there are eight phases when a nonzero Heisenberg term is added, including four Luttinger liquid phases, an FM phase, a Néel phase, a narrow "d-Spiral" phase in which the spins align in a distorted-spiral pattern, and a "D 3 -breaking" phase. While good agreements are reached between analytic and numerical calculations for all other phases, numerics provide evidence for the coexistence of the proposed two D 3 -breaking orders in the "D 3 -breaking" phase, though the two orders are expected to occupy different subregions based on a combination of perturbative Luttinger liquid, symmetry and classical analysis. Whether this coexistence is a truth or a finite size artefact is worth further studies. Our comprehensive study of the phase diagram of the 1D generalized Kitaev model provides a road-map to the exotic physics in Kitaev materials.

Appendix A: The Hamiltonians in different frames
In this section, we spell out the terms in the Hamiltonians in different frames. In general, we write the Hamiltonian H as H = L j=1 H j,j+1 where H j,j+1 is the term on the bond between the sites j and j + 1. The forms of H j,j+1 will be written explicitly.
In the unrotated frame, the form of H j,j+1 has a two-site periodicity. We have In the six-sublattice rotated frame, the form of H j,j+1 has a three-site periodicity. We have In the four-sublattice rotated frame, the form of H j,j+1 has a four-site periodicity. We have In the literature, another Γ -term sometimes is also considered [1], which is defined as the following in the original frame Since Γ is much smaller than the other three couplings [1], the Γ term is neglected in this work.
Appendix B: Transformation rules of the SU(2)1 WZW fields In this section, we work out the symmetry transformations of the WZW primary field and its descendent fields. The strategy is to first work within the fermion representation, then use the nonabelian bosonization formula to translate them into the WZW language.
We first summarize the transformation properties of the primary and descendent fields in the SU(2) 1 WZW theory. The derivations are left for the subsequent subsections. The scaling fields in the SU(2) 1 WZW theory are known to be , in which (AB)(w) in the holomorphic sector is the O((z − w) 0 ) term in the operator product expansion (OPE) of A(z)B(w) as a Laurent series in z − w around the point w, and definitions are similar for the anti-holomorphic sector. The transformation laws of g and J L , J R under time reversal T , spatial translation T a , inversion I and spin rotation R ∈ SO(3) are summarized as in which x is the spatial coordinate; R α β (α, β = x, y, z) is the matrix element of the 3 × 3 rotation matrix R; (x) = trg(x) is the dimer order parameter; and N (x) = itr(g(x) σ) is the Néel order parameter [3].

Spatial inversion
In fermion language, Hence From this, we obtain, Using the bosonization formula, we have When the charge boson is gapped out, e −i √ 2πφ(x) is purely imaginary, hence e i √ 2πφ(x) = −e −i √ 2πφ(x) . This yields

Time reversal
In fermion language, Then This gives By canceling the charge boson, we obtain (B20)

Spin rotation
The transformation law under spin rotation R is in which the 2 × 2 matrix U is the SU(2) representation of the rotation R.
In summary, the transformation rules are Sometimes, it is more convenient to use the following combinations of g = tr(g), N = itr(g σ).
Appendix C: Redundancies in the descendent fields in the SU(2)1 WZW model In this section, we show that the existence of null fields in the SU(2) 1 WZW model leads to redundancies in the descendent fields. We focus on several cases that are relevant for our purposes, including the quadratic WZW current terms, and the dimension-3/2 fields J α L N β + J β L N α where α = β. The holomorphic sector is taken as an example, and the antiholomorphic sector can be treated in a similar manner.
First let's consider the quadratic WZW current terms. We compute (J a L J b L ) by working out the O(1) terms in the OPE J a L (x + a)J b L (x). In what follows, all higher order terms in a will be neglected. The free fermion contraction rule is [2] The OPE J a L (x + a)J b L (x) can be calculated as where n Lα =: ψ † Lα ψ Lα : (α =↑, ↓). Hence all the three (J a L J a L ) (a = x, y, z) are equal given by Eq. (C3). As a result, (J a L J a L ) = 1 3 ( J L · J L ). If a = b, the first term in the last line of Eq. (C2) vanishes. Taking the O(1) part of Eq. (C2), we obtain Using z = τ + ix, we also have If we only keep the O(1) terms, then it becomes Note that due to Fermi statistics, we must have β = δ. For our purpose, let's take c = x, d = y in accordance with Eq. (F18). Then the coefficient in Eq. (C8) becomes since both σ x and σ y are off-diagonal. Take β =↑ as an example. Then it is clear that the coefficient vanishes. The same conclusion holds for β =↓, and for the antiholomorphic term J c Appendix D: Operator algebra of the SU(2)1 WZW model

The SU(2)1 affine algebra
In this subsection, we fix the normalization conventions for the SU(2) 1 affine algebra. The commutation relations of the affine generators can be obtained from the Ward identities [2]. Alternatively, they can be determined from the 1D Dirac fermions, which we briefly review in this subsection.
We work with radial ordering [2]. The OPE of the free fermion fields are The WZW current operators are in which a = x, y, z. The anomalous constant term of the OPE J a (z)J b (w) (similar forJ) can be obtained from the current-current correlation function [4]. The result of the OPE is in which abc is the rank-3 totally antisymmetric tensor. Define the affine generator as J a n = 1 2πi dzz n J a (z).
It is straightforward to obtain the following commutation relations from Eq. (D3) [2], Note the 1/2 factor in the anomalous term in Eq. (D5), which is a consequence of the choice of normalization in Eq. (D2).

Operator algebra of the SU(2)1 WZW model
In this part, we use the affine symmetry to determine the OPE g αβ (z,z)g γδ (w,w) up to quadratic order of the current operators. Throughout this section, repeated indices indicate summation unless otherwise stated.
First, let's figure out what kinds of primary fields appear in the OPE. Notice that 1 and g are the only two primary fields in the SU(2) 1 WZW model. Thus the components of g αβ (z,z)g γδ (w,w) on these two primary fields are given by the vacuum expectation values g αβ (z,z)g γδ (w,w) and g αβ (z,z)g γδ (w,w)g −1 λµ (0, 0) , respectively. However, the three-point function g αβ g γδ g −1 λµ vanishes due to the chiral rotation symmetry. In fact, consider the holomorphic sector, then g αβ g γδ g −1 λµ is spin-1/2 ⊗ spin-1/2 ⊗ spin-1/2, which does not contain an SU(2) singlet. On the other hand, the vacuum is SU(2) invariant, so the expectation value is zero. As a result, the OPE g αβ (z,z)g γδ (w,w) only contains the primary field 1.
In the case of the Virasoro algebra, once the three-point functions of the primary fields are known, the full OPE can be obtained by using the conformal symmetry and the Virasoro algebra [2]. This is also true for the affine algebra. Thus we conclude that the OPE g αβ (z,z)g γδ (w,w) has no components on all affine descendent fields of the primary field g. Hence it only contains descendent fields of 1, which means that g αβ (z,z)g γδ (w,w) can be written in terms of affine current operators.
We write in which and for some field φ, (J a −n φ)(z) is defined as In Eq. (D7), φ(z) is taken to be the identity field 1.
Acting with the OPE on the vacuum |Ω , we have in which |g γδ >= g γδ (0, 0)|Ω . In what follows, we only keep the holomorphic part. The full expression is a product of the holomorphic part and the anti-holomorphic part. Define then we have We are going to show that all the coefficients of β K αγ can be determined from β 0 αγ = β αγ . We now act J a n on g αβ (z,z)|g γδ with n > 0. Then J a n g αβ (z,z)|g γδ = 1 On the other hand, since J a n (n > 0) annihilates |g γδ [2], we obtain J a n g αβ (z,z)|g γδ = [J a n , g αβ (z,z)]|g γδ .
We also note that similar treatment can be performed for the anti-holomorphic sector, and the result is

First order expansion of the OPE
Let's first consider the holomorphic sector. Take N = 0, n = 1, then J a n |1 αγ = ( Using (in which the superscript 1a is 1 a for short) and the affine algebra Eq. (D5), we have Thus it is clear that Then to lowest order, the holomorphic part of the OPE g αβ (z,z)g γδ (w,z) is given by Notice that in the above expression is just the WZW current at the point w, hence g αβ (z,z)g γδ (w,z) = 1 z 1/2 αγ − a=x,y,z (σ a iσ 2 ) αγ J a (w) + ... .
Also note that we only know g αβ (z,z)g γδ (w,w) ∝ αγ βδ 1 |z − w| . (D29) In the above calculation of the OPE, we have fixed the coefficient of this vacuum expectation value to be 1, which amounts to a multiplicative rescaling of the primary field g. Proceeding exactly similarly for the anti-holomorphic part, we obtain Combining the holomorphic and anti-holomorphic parts together, the result up to first order in the WZW currents is

Second order expansion of the OPE
We first consider the holomorphic sector. We write in which γ a and β ab are . We also note that in Eq. (D32), and in which (J a J b )(w) means the zeroth order term in the expansion of the OPE J a (z)J b (w) in powers of z − w.
First we take N = 0, n = 2 in Eq. (D20). Then On the other hand, we can use the affine commutation relations to calculate J c 2 |2 αγ as This gives or as a matrix identity Next we take N = 1, n = 1. Then On the other hand, using the affine algebra we get Thus we get the equation In conclusion, the equations determining γ a and β ab are The first equation can be simplified by multiplying −i cde on both sides. We then arrive at Plugging the first equation in Eq. (D43) into the second, we obtain We have three comments here. 1. We can only solve the sum d β dd as 1 2 iσ 2 , and the individual β xx , β yy , β zz cannot be determined individually. However, according to Sec. C, (J a J a ) are all equal for a = x, y, z. Hence we can well take 2. The information about the off-diagonal symmetric part of the tensor β ab is missing, i.e., there is no information about 1 2 (β ab + β ba ) (a = b). However, according to Sec. C, (J a J b ) + (J b J a ) = 0 when a = b. Thus this is not a problem.
3.We only know a combination of the anti-symmetric part of β ab and γ c as given by the first equation in Eqs. (D43), and we are not able to solve them individually. However, according to Sec. C, Thus we can drop (J c −2 1) and absorb it into the antisymmetric combination i abc (J a In summary, Eq. (D43) can be simplified as The solution is clearly, By doing exactly similar treatment to the anti-holomorphic sector, we obtain In summary, we are able to write down the OPE of g up to quadratic order in the WZW currents, as i.e., .
Using Eq. (D51), we obtain the OPE between the Neel order fields, in which λ, µ = x, y, z. We note that the normalization in Eq. (D52) is chosen as Eq. (D29), which is different from the one used in the main text which is On the other hand, the 1D spin-1/2 free Dirac fermion can be decomposed into a U(1) charge boson and an SU(2) WZW boson according to the nonabelian bosonization theory. Then the OPE calculated in the bosonic language must coincide with the fermionic calcuation. By properly taking into account the U(1) charge sector, we have verified that the fermionic approach gives the same result as Eq. (D52). they can all be generated by the above five operations: Therefore, the symmetry group of the Kitaev-Heisenberg chain is The group structure of G 0 will be shown in Sec. E 2 a to be G 0 D 4d (Z 2 2Z).
b. Γ = 0 Next we turn on a nonzero Γ and study the symmetry group of a general Kitaev-Heisenberg-Gamma chain. The symmetry transformations are now 1. T where in particular, with the presence of the cross terms S α i S β i+1 (α = β) in the Hamiltonian, the operations R(α, π) (α = x, y, z) are no longer symmetries. We thus conclude that the symmetry group G N is By a similar analysis with the G 0 case discussed in the main text, the group structure of G N when Γ = 0 is in which from left to right, Z 2 × Z 2 = <T > × <R(n N , π)T a >/<T 2a >, Z 2 = <T a I>, and 2Z = <T 2a >.
On the other hand, in the main text we have discussed the same model but in the six-sublattice rotated frame. The group structure in the six-subalttice rotated frame has been worked out to be G 1 ∼ = D 3d 3Z. Notice that G 1 and G N must be isomorphic since they are essentially the same group up to a unitary transformation. An explicit verification of this isomorphism is included in Appendix E 3. To determine the group structure of G 0 , first notice that N 0s = <T 2a , T a I> = <T a I> <T 2a > is a normal subgroup of G 0 , hence G 0 can be written as in which G 0s = <T, R(ẑ, π 2 )T a , R(ŷ, π)>/<T 2a >.
Then by defining a = R(ẑ, π 2 )T a , b = R(ŷ, π), it is straightforward to verify that a 4 = b 2 = (ab) 2 = e mod T 2a . Since the generator-relation representation for the group D n is D n = <α, β|α n = β 2 = (αβ) 2 = e>, (E8) and the relations in Eq. (E8) are satisfied for the two generators of G 0s /<T >, we see that G 0s /<T > must be a subgroup of D 4 . On the other hand, since the actions of {1, a, a 2 , a 3 , b, ab, a 2 b, a 3 b} (E9) restricted within the spin space are which are all distinct operations, there must be at least eight distinct group elements in G 0s /<T >. But the order of D 4 is eight, hence we conclude that G 0s /<T > D 4 , i.e., G 0s D 4d . In conclusion, the group structure of G 0 is of a regular triangle. There are three rotations and three reflections, in which the rotation group C 3 is a normal subgroup of D 3 . Thus D 3 can be written as in which Z 2 is the group generated by a reflection. For example, ab is a refelection and a 2 = a −1 is a rotation. Hence, This gives On the other hand, <ab> <a 2 > is also a normal subgroup of G 1 . This is because T 3a R I IT −1 3a = (R a T a ) −6 R I I. Therefore, G 1 can also be written as or alternatively, Using we obtain which is exactly G N .
Appendix F: Symmetry analysis of the low energy field theory 1. Symmetry analysis in the "LL1" phase In this section, we perform a symmetry analysis of the low energy field theory in the "LL1" phase. The symmetry group of the system in the six-sublattice rotated frame has been worked out to be G 1 = <T, R a T a , R I I> ∼ = D 3d 3Z. We will exhaust all the symmetry allowed relevant and marginal terms in the low energy field theory. Throughout this section, we work in the six-sublattice rotated frame unless otherwise stated.
All operators with scaling dimensions not greater than 2 in the SU(2) 1 WZW model are listed as follows, Next we inspect which of these terms are compatible with the D 3d 3Z symmetry.
1) The dimension 1/2 operators are forbidden since they change sign under odd powers of T a (see Ref. [3]) and in the our case T 3a is a symmetry of the system.
2) The dimension 3/2 operators also change sign under T 3a .
3) To analyze the dimension 1 operators, we first classify the current operators J L , J R according to the group <R a T a , R I I>/<T 3a > ∼ = D 3 . The three dimensional vector representation J s (s = 1, −1 corresponding to L, R) of D 3 can be decomposed as V 1s ⊕ V 2s , in which where "Span{...}" denotes the R-linear space spanned by the vectors in the bracket, and The actions of R a T a , R I I in V 1s , V 2s can be worked out as and respectively, in which R a T a keeps the chiral sector while R I I flips the chiral index, and Notice that in V 2s , both R a T a and R I I correspond to orthogonal matrices in the basis of {J 1 s , J 2 s }. According to Eq. (B1), J L − J R is invariant under the time reversal operation. On the other hand, because of Eqs. (F4,F5), the only combination compatible with the D 3 symmetry is J 0 is allowed in the low energy theory. 4) Finally, we analyze the dimension 2 operators. We consider the actions of the time reversal operation T and the group D 3 = <R a T a , R I I> separately. Since time reversal switches the chiral indices L and R, the combinations J α L J β L ± J α R J β R , J α L J β R ± J β L J α R are even (for +) and odd (for −) under the time reversal operation. Next consider the constraints imposed by the D 3 group. Neglecting the chiral index, since J 0 is an A 2 representation of D 3 , J 0 J 0 is apparently D 3 -invariant. Hence, the terms are allowed by D 3d = D 3 × <T >. On the other hand, {J 1 , J 2 } form an E representation of the D 3 group. Using the multiplication rule we know that J 1 J 1 + J 2 J 2 and J 1 J 2 − J 2 J 1 are within the A 1 and A 2 representations, respectively. Out of the A 1 representation, the following two terms do not change sign under the time reversal operation, hence allowed in the low energy theory. As to the A 2 representation, an odd combination of the two chiral sectors must be formed to ensure the invariance under R I I, but then it spoils the time reversal invariance. In summary, all the independent symmetry allowed relevant and marginal terms can be taken as in which we have used the relation J s · J s = J 0 s , J s · J s } as the two linearly independent terms.
We make a final remark that all the conclusions in this subsection are based on symmetry analysis, hence not dependent on the detailed forms of the microscopic Hamiltonian. Thus the analysis applies also to the KHΓΓ chain where a small Γ term is included, and it holds even when the terms beyond nearest neighbor couplings are taken into account. Of course, the coupling constants in the low energy theory will be renormalized by these additional terms.

Symmetry analysis in the Néel phase
In this section, we perform a symmetry analysis of the low energy field theory in the "Néel" phase. The symmetry group of the system in the original frame when Γ = 0 has been worked out to be G N = <T, T a I, R(n N , π)T a >. We will exhaust all the symmetry allowed relevant and marginal terms in the low energy field theory. Throughout this section, we work in the original frame unless otherwise stated.
We start from the AFM Heisenberg model whose low energy properties are described by the SU(2) 1 WZW model, and treat the Kitaev and Gamma terms as small perturbations, which applies to a neighborhood of the AFM1 point. Here we will focus on the case of a nonzero Γ.
The symmetry allowed relevant and marginal terms can be identified in a similar manner as Sec. F 1. The list of terms with dimensions not greater than two is the same as Eq. (F1).
1) For dimension 1/2 operators, is forbidden by R(ẑ, π 2 )T a , and N α is forbidden by T . 2) Time reversal symmetry restricts the dimension 1 operators to be the combinations J α L − J α R . However, this combination changes sign under T a I.
3) The dimension 3/2 operators are J α L , J α R , J α L N β , J α R N β . Time reversal restricts them to be the combinations (J α L − J α R ) and (J α L + J α R )N β . However, according to the transformation laws in Eqs. (B2,B3), both combinations change sign under T a I, hence forbidden. 4) To discuss the dimension 2 operators, let's first study the linear space span{J x s , J y s , J z s } (s = L, R) as a representation of the group <R(n N , π)> ∼ = Z 2 , where span{· · ·} represents the R-linear space spanned by the vectors within the bracket. The three operators J x s , J y s , J z s can be recombined into which have R(n N , π) eigenvalues −1, 1, −1, respectively. Therefore, the quadratic terms invariant under R(n N , π)T a are given by Time reversal and T a I further require L and R to appear symmetrically. In summary, in addition to the terms that are already in the SU(2) invariant AFM Heisenberg model, the linearly independent nontrivial additional terms in the KHG chain are We note that the term J x L J y R + J y L J x R does not appear in the first order perturbation low energy Hamiltonian derived in the main text. However, it can be generated upon the RG flow. In fact, this term appears in the low energy Hamiltonian of the KHΓΓ chain. At low energies, generically, all symmetry allowed terms will be generated, and the KHΓ and KHΓΓ chains contain the same set of terms though with different coefficients.

Symmetry analysis in the "LL2" phase
In this section, we analyze the phase diagram of the Kitaev-Heisenberg chain based on a combination of symmetry and RG analysis. We will perform a symmetry analysis of the low energy field theory in the "LL2" phase. The symmetry group of the system in the original frame when Γ = 0 has been worked out to be G 0 = <T, T 2a , T a I, R(ẑ, π 2 )T a , R(ŷ, π)>. We will exhaust all the symmetry allowed relevant and marginal terms in the low energy field theory. Throughout this section, we work in the original frame unless otherwise stated.
All the relevant and marginal terms are given in Eq. (F1) In summary, the symmetry allowed terms are Clearly, the low energy Hamiltonian is still of the XXZ type.
On the other hand, we can do a first order perturbation treatment to the Kitaev term. Explicit calculations show that (for details, see Sec. G 2), Hence, we see that all the symmetry allowed terms indeed appear within the low energy Hamiltonian. It is then straightforward to determine the phase diagram. The RG flow is again of the KT type. The phase diagram is shown in Fig. 25. When K > 0, i.e., the "LL2" phase, the system is described by the Luttinger liquid theory with an emergent U(1) symmetry at low energies. On the other hand, when K < 0, the system develops a Néel order along z-direction.

Symmetry analysis in the "LL3" phase
In this section, we perform a symmetry analysis of the low energy field theory in the "LL3" phase. The symmetry group of the system in the four-sublattice rotated frame when Γ = 0 has been worked out to be G 3 = <T, R(ŷ, π)T a I, R(ẑ, − π 2 )T a >. We will exhaust all the symmetry allowed relevant and marginal terms in the low energy field theory. Throughout this section, we work in the four-sublattice rotated frame unless otherwise stated.
The list of terms with dimensions not greater than two is the same as Eq. (F1).
1) The dimension 1/2 operators and N change sign under R(ẑ, − π 2 )T a and T , respectively. Hence both are forbidden.
2) The dimension 1 operators are restricted to the form J α L − J α R (α = x, y, z) due to time reversal symmetry, in which only J z L − J z R is allowed by R(ẑ, − π 2 )T a . It can be seen that J z L − J z R is invariant under R(ŷ, π)T a I. Therefore, the term J z L − J z R is allowed in the low energy theory. 2) The dimension 3/2 operators are restricted to the forms (J α L − J α R ) and (J α L + J α R )N β due to time reversal symmetry. The symmetry operation R(ẑ, − π 2 )T a forbids (J α L − J α R ) . This is because changes sign under T a , and it is impossible for J α L − J α R to change another sign under the rotation R(ẑ, − π 2 ) since the π/2-rotation does not have the eigenvalue of −1.
However, it is possible for suitable combinations of (J α L + J α R )N β to change sign under R(ẑ, − π 2 ). In a vector representation {v α } α=x,y,z of the SO(3) group, the eigenvalues of R(ẑ, − π 2 ) are 1, i, −i with eigenvectors v z , v x + iv y , v x − iv y , respectively. Hence, among the rank-2 tensors (J α L + J α R )N β , the following combinations change sign They can be alternatively rewritten into real combinations as which are allowed by R(ẑ, − π 2 )T a . However, the term (J x L + J x R )N x − (J y L + J y R )N y changes sign under R(ŷ, π)T a I, and (J x L + J x R )N y + (J y L + J y R )N x vanishes in the SU(2) 1 WZW model (see Sec. C). Hence, we conclude that there is no dimension-3/2 term in the low energy theory.
3) For the dimension 2 operators, let's focus on the left-right cross terms J α L J β R . The terms invariant under R(ẑ, − π 2 )T a are J z L J z R , J x L J x R + J y L J y R , and J x L J y R − J y L J x R . All these three terms are allowed by R(ŷ, π)T a I. However, the third term is forbidden by time reversal symmetry.
In summary, the nontrivial symmetry allowed terms with dimensions not greater than two are in which J z L − J z R can be removed by a chiral rotation as discussed in the main text, and the system is again of an XXZ type at low energies.
and the OPE formula for the Néel order fields, we obtain u.c.
Next consider the cross term. Using 1 a S y 1 = D 1 (J y L (x) + J y R (x)) + C 1 (−) x a 1 a itr(g(x)σ y ), 1 a S z 2 = D 1 (J z L (x + a) + J z R (x + a)) + C 1 (−) x a +1 1 a itr(g(x + a)σ z ), we obtain 1 a 2 S y 1 S z We note a subtlety here. Notice that the spin operators commute when the two operators are on different sites. However, Eq. (G7,G8) do not satisfy this. In fact, there is no reason to expect that such commutativity still holds when the operators are restricted to a low energy subspace. On the other hand, we do have an ambiguity in writing S α i S β j in terms of the low energy degrees of freedom, since S α i S β j and S β j S α i are different. To regularize such discrepancy, we use the symmetrized expressions 1 2 (S α i S β j + S β j S α i ), which removes the antisymmetric terms within Eq. (G7,G8). An alternative way is to keep using Eq. (G7), but consider the sum of pairs 1 2 (S y 1 S z 2 + S y 3 S x 2 ), which is ensured by the R I I symmetry operation. Then apparently, combining S y 3 S x 2 with S x 2 S y 3 together lead to an effective symmetrization. As a result, we obtain 1 a 2 (S y 1 S z 2 + S z 1 S y 2 ) Summing the three bonds up, we obtain in which J {z J y} + J {y J x} + J {x J z} = 1 2 (3J 0 J 0 − J · J) is used. Hence, these terms contribute to H J the following expression, Combining Eqs. (G5,G11) together, we arrive at the low energy Hamiltonian close to the emergent SU(2) 1 line as discussed in the main text, and the coefficients u's are given by In this section, we work in the original frame. We will project the Kitaev term to the low energy subspace, and confirm that indeed all the symmetry allowed terms will emerge. The Kitaev term is Using the nonabelian bosonization formula, we have itrg(x + a) itrg(x)(J L (x + a) + J R (x + a)) + c 2 (2πa) 2 tr(g(x)σ x )tr(g(x + a)σ x ). (G15) Using the following OPE, J L (w)g(z,z) = − 1 2 σg(z,z) w − z + ( J L g)(z,z), J R (w)g(z,z) = 1 2 g(z,z) σ w − z + ( J R g)(z,z), we obtain, J x L (x)tr(g(x + a)σ x ) = 1 4πia tr(g(x + a)) + tr((J x L g)σ x ), J x R (x)tr(g(x + a)σ x ) = 1 4πia tr(g(x + a)) + tr((J x R g)σ x ), tr(g(x)) + tr((J x L g)σ x ), tr(g(x)) + tr((J x R g)σ x ).
Therefore, the second and the third terms in Eq. (G15) sum up to i c 2πa (−) 1 4πia trg × 4. The fourth term in Eq. (G15) can be calculated using the OPEs Eqs. (D52) as In summary, Similarly, S y 2i S y 2i+1 = 1 2 a dx − 2c 2 (2πa) 2 + c 2πa 2 tr(g) + J y L J y L + J y R J y R + 1 3 c 2 ( J L · J L + J R · J R ) + (2 + 2c 2 )J y L J y Summing the two terms up, we arrive at the low energy Hamiltonian close to the AFM1 point along the circular boundary of Fig. 1 in the main text.

Low energy Hamiltonian in the "Néel" phase
In this section, we derive the low energy Hamiltonian in the "N 'eel" phase close to the AFM1 point. Throughout this section, we work in the original frame.
Trying the ansatz Eqs. (I4,I5) are reduced to Since there are three variables a, b, λ and an equal number of equations, a solution exists generically.
In what follows, we only discuss the special case K = 0, Γ > 0 for illustration. Other cases can be solved exactly similarly. The solution of Eq. (I7) is To check if this is a minimum of the free energy, the eigenvalues of the Hessian matrix can be calculated in a perturbative expansion over Γ. We have determined the two lowest eigenvalues to be 1 √ 2 Γ and √ 2Γ , which are both positive when Γ > 0.

Appendix J: DMRG numerical results for the Luttinger parameters
Recall that as discussed in the main text, for a finite size system with an open boundary condition, the energy density h(x) contains a uniform part E U (x) and a staggered part E A (x), where in which x = ja (j 1) is the distance measured from the boundary of the system. The provides a method to accurately determine the Luttinger parameter K. We apply this method to the "LLi" (i = 1, 3, 4) phases.
1. The "LL1" phase Fig. 27 shows E A thus obtained vs sin(πx/L) on a log-log plot at a representative point (θ = 0.42π, φ = 0.75π) in the "LL1" phase, where DMRG numerics are performed on a system of L = 96 sites with an open boundary condition. As can be seen from Fig. 27, an excellent linear fit can be obtained from which the Luttinger parameter is determined to be K 0.66. Fig. 28 shows E A (x) vs sin(πx/L) on a log-log scale at a representative point in the "LL3" phase (θ = 0.64π, φ = 0.11π), where a good linear fit is obtained giving a Luttinger parameter equal to 0.591.

The "LL4" phase
The plots of log E A (x) vs log sin(πx/L) at a representative point in the "LL4" phase close to the "peninsular end" and a nearby point in the "FM" phase are shown in Fig. 29, where DMRG numerics are performed on open systems of L = 96 sites and x is the distance measured from the boundary. As can be seen from Fig. 29, while a good linear fit can be obtained for the point (θ = 0.57π, φ = 0.30π) within the "LL4" phase, no linear relation can be fitted for the point (θ = 0.6π, φ = 0.30π) which is in the "FM" phase.