Twice hidden string order and competing phases in the spin-1/2 Kitaev–Gamma ladder

Abstract

Finding the Kitaev spin liquid in candidate materials involves understanding the entire phase diagram, including other allowed interactions. One of these interactions, called the Gamma (Γ) interaction, causes magnetic frustration and its interplay with the Kitaev (K) interaction is crucial to comprehend Kitaev materials. Due to the complexity of the combined KΓ model, quasi-one-dimensional models have been investigated. While several disordered phases are found in the 2-leg ladder, the nature of the phases are yet to be determined. Here we focus on the disordered phase near the antiferromagnetic Γ limit (denoted by AΓ phase) next to the ferromagnetic Kitaev phase. We report a distinct non-local string order parameter characterizing the AΓ phase, different from the string order parameter in the Kitaev phase. This string order parameter becomes evident only after two unitary transformation, referred to as a twice hidden string order parameter. The related entanglement spectrum, edge states, magnetic field responses, and the symmetry protecting the phase are presented, and its relevance to the two-dimensional Kitaev materials is discussed. Two newly identified disordered phases in the phase diagram of KΓ ladder is also reported.

Introduction

Over the course of investigating spin S = 1/2 two-dimensional (2D) honeycomb Kitaev materials^1,2,3,4,5, candidates of long-sought quantum spin liquids, an additional bond-dependent spin exchange term named Gamma (Γ) interaction⁶ was found along with the bond-dependent Kitaev (K) interaction^7,8. Unlike the standard Heisenberg interaction, both the Kitaev and the Γ interactions are highly nontrivial and extremely frustrated. While the same sign of K and Γ cancels the frustration leading to magnetically ordered phases, the regions with different signs of these interactions are highly frustrated⁵. Since these two interactions are known to dominate over other symmetry-allowed interactions in the emerging candidate material α-RuCl₃, the combined Kitaev–Gamma (KΓ) model has attracted considerable attention in the theoretical community and it is widely accepted^9,10,11 that a realistic description of α-RuCl₃ should be sought in the regime with antiferromagnetic (AFM) Gamma (Γ > 0) and ferromagnetic (FM) Kitaev (K < 0) interactions. Understanding the phases that arise in this frustrated regime of the 2D honeycomb KΓ-model is therefore of crucial importance.

This particular region of the KΓ-model has been extensively studied by a range of numerical methods⁵. Despite detailed studies, the nature of the phase next to the FM Kitaev spin liquid arising due to AFM Γ interaction remains controversial. Most studies have found that it is in a disordered phase^{12,13,14,15,16,17,18,19}, denoted by KΓSL for KΓ spin liquid, or nematic paramagnets, but functional renormalization approaches found magnetically ordered phases²⁰ while variational Monte Carlo calculations²¹ observed a narrow disordered phase next to the FM Kitaev spin liquid with most of the antiferromagnetic Γ region dominated by a zig-zag ordered phase.

Motivated by such discrepancy, another approach to investigate the 2D limit of the KΓ model was taken by starting from low-dimensional models with the hope of furthering the understanding of the honeycomb model by determining the phases of n-leg (brick-wall) models. Despite the obvious challenge in connecting the two limits, it is reasonable to expect that potential spin liquid phases arising in the honeycomb model should correspond to regions where the n-leg models display disordered phases. Such an approach was employed earlier for the pure Kitaev model²². Disordered phases in the anisotropic Kitaev 1-leg chain were found, and they were characterized by non-local string order parameters (SOPs)²². It has also been shown that the isotropic Kitaev 2-leg ladder model exhibits a disordered phase, characterized by an unconventional SOP different from that of the anisotropic chain Kitaev phase²³.

The one-dimensional (1D) chain and ladder version of the KΓ model were investigated numerically with very high precision, as the reduced dimensionality allows access to bigger system sizes^{24,25,26,27,28}. In the 1D chain model, it was found that the pure Gamma model belongs to a Luttinger liquid phase governed by the gapless hidden SU(2) Heisenberg chain, a fact revealed after a 6-site transformation, i.e, a duality mapping^24,29. A study of the same KΓ model on a quasi 1D ladder using DMRG and iDMRG techniques^18,30 found a magnetically disordered phase, possessing a small gap near the AFM pure Gamma limit. This phase surrounding the pure AFM Gamma point next to the FM Kitaev phase (denoted by FK) was referred to as the AΓ phase. Even though the AΓ phase in the ladder occurs in the same part of the phase diagram as the proposed KΓSL in the 2D limit, a distinct name was introduced, as it is yet to be determined how the AΓ phase is connected to the 2D limit KΓSL.

While it is clear that there is no magnetic order in this phase, the precise nature of the AΓ phase has not yet been settled due to its complex nature. The presence of a gap indicates that the AΓ-phase is likely a symmetry-protected topological (SPT) phase^31,32,33. If so, it is of important to identify a corresponding SOP, edge states, and the symmetry that protects this phase, which are characteristic of the SPT, and to determine how this phase responds to an external magnetic field. If strong evidence for a non-trivial SPT nature of the AΓ phase can be established, it would establish a next step to the proposed KΓSL in the 2D limit. In the following sections, we will systematically examine these inquiries and provide comprehensive responses. To perform a thorough analysis, we start by reviewing the full phase diagram prior to focusing on the nature of the AΓ phase. As shown below, we also report two other disordered phases.

Results

Model

The KΓ Hamiltonian is given by

$${H}_{K\Gamma }=\mathop{\sum}\limits_{{\langle i,j\rangle }_{\gamma \in (x,y,z)}}K{S}_{i}^{\gamma }{S}_{j}^{\gamma }+\Gamma \left({S}_{i}^{\alpha }{S}_{j}^{\beta }+{S}_{i}^{\beta }{S}_{j}^{\alpha }\right)$$

(1)

where (α, β) takes on the values (y, z)/(x, z)/(x, y) for γ = x/y/z, and 〈i, j〉 refers to the nearest neighbor sites. An alternative representation of the honeycomb lattice is as a brick-wall lattice²², and its two-leg limit with periodic boundary conditions is simply a ladder shown in Fig. 1a. The dotted bonds indicate Kitaev z-bonds arising from periodic boundary conditions, and we shall always take such bonds to be identical to the regular (solid) z-bonds, in which case the honeycomb strip can be viewed as a regular rectangular ladder.

Fig. 1: The KΓ ladder and corresponding phase diagram.

a A strip of the KΓ honeycomb lattice corresponding to a two-leg KΓ ladder with alternating x and y bonds along the leg and z-bond between the chains, with the numbering of the sites used throughout the paper. The dotted z-bonds arise from imposing periodic boundary conditions in the direction perpendicular to the ladder. However, all z-bonds are taken to have equal strength. The dashed red line indicates the partition used for ρ_B while the dashed blue line indicates the partition used for ρ_A. b Schematic phase diagram of the KΓ ladder, Eq. (1), for Γ > 0.

We parameterize the model by taking

$$K=\cos \phi ,\quad {{{\rm{and}}}}\quad \Gamma =\sin \phi ,$$

(2)

and interpolate between the Kitaev and Γ interactions by varying ϕ. Our main interest is in the region with ϕ/π ∈ [0, π] where Γ > 0 and the Kitaev term, K, changes from AFM to FM at ϕ = π/2, as this region is relevant to most two-dimensional (2D) Kitaev candidate materials. The total number of sites in the ladder (including both legs) is denoted by N.

Phase diagram

A full phase diagram is shown in Fig. 1b. It is obtained by various quantities presented in the next subsection. Moving from ϕ = 0 to π, the AFM Kitaev phase (denoted by AK), a FM phase denoted by FM${}_{{U}_{6}}$, AΓ, and FK phases are found consistent with the earlier works^18,30. At the special point ϕ = π/4, the FM${}_{{U}_{6}}$ phase can be mapped to the ferromagnetic Heisenberg ladder by applying a local unitary U₆ transformation (see Supplementary Note 1), and is therefore gapless. However, for ϕ ≠ π/4 a small gap appear, as we show in Supplementary Note 3. Surprisingly, two additional phases denoted by SPT_α and SPT_β can be identified between the FM${}_{{U}_{6}}$ and AΓ phases. As we will show below, they are magnetically disordered and display the characteristics of SPT phases, i.e., doubled entanglement spectrum. Beyond ϕ = π, the rung singlet phase denoted by RS${}_{{U}_{6}}$ phase delineates the FK phase.

In addition to the expected Kitaev phases AK and FK the appearance of the FM${}_{{U}_{6}}$ and RS${}_{{U}_{6}}$ phases are well established in the KΓ honeycomb and n-leg models. After the local U₆ transformation²⁹, corresponding to local spin rotations, at ϕ = π/4 and 5π/4, i.e., K = Γ, the KΓ 2D honeycomb model is equivalent to the FM and AFM Heisenberg model, respectively. The application of the U₆ transformation is specified in Supplementary Note 1. To understand the nature of the other three phases, AΓ, SPT_α and SPT_β, we first performed a detailed analysis of the entanglement spectrum.

Entanglement spectrum

Our results for the entanglement spectrum, as well as for the susceptibility, ${\chi }_{\phi }^{e}$ are shown in Figs. 2, 3. Due to the complexity of the phase diagram, we split it into two regions centered around the AΓ phase: Fig. 2 represents the left part of the AΓ and Fig. 3 the right part of the AΓ phase. We consider two different partitions shown in Fig. 1a corresponding to ρ_A and ρ_B. In both Figs. 2, 3 iDMRG results are shown for both partitions with ρ_A shown in panel (a) and ρ_B in panel (b) with ${\chi }_{\phi }^{e}$ in panel (c).

Fig. 2: Entanglement spectrum and susceptibility for ϕ/π < 0.5.

iDMRG results for: a Entanglement spectrum from ρ_A, red dashed line in Fig. 1a. b Entanglement spectrum from ρ_B, blue dashed line in Fig. 1a. The numbers refer to the degeneracy of the eigenvalue. c ${\chi }_{\phi }^{e}$. Two distinct phases, SPT_α and SPT_β, are visible between the FM${}_{{U}_{6}}$ and AΓ phase. The red shading between ϕ/π = 0.428-0.442 denotes a transitional region of limited convergence due to a field instability.

Fig. 3: Entanglement spectrum and susceptibility for ϕ/π > 0.5.

iDMRG results for: a Entanglement spectrum from ρ_A, red dashed line in Fig. 1a. b Entanglement spectrum from ρ_B, blue dashed line in Fig. 1a. The numbers refer to the degeneracy of the eigenvalue. c ${\chi }_{\phi }^{e}$. The AΓ FK and RS${}_{{U}_{6}}$ phases are clearly delineated.

To the left of the AΓ-phase, two other phases SPT_α and SPT_β are clearly separated from the FM${}_{{U}_{6}}$ and AΓ, as noted from their entanglement spectrum for ρ_A. Furthermore, the double degeneracy of the entanglement spectrum shown in Fig. 2b from ρ_B is a clear signal of a SPT phase^{33,34,35,36,37}. Note that there is only a weak signature in ${\chi }_{\phi }^{e}$ of the transition between the two phases, visible at ϕ = 0.3956π in Fig. 2c, corresponding to a small discontinuity in ${\chi }_{\phi }^{e}$. While the FM${}_{{U}_{6}}$ phase appears abruptly at ϕ = 0.3823π, the transitions between the AΓ phase and SPT_β phase at ϕ ~ 0.435π is not discernible in ${\chi }_{\phi }^{e}$. The blurriness of the transition is likely due to a field-induced phase that pinches off to a single point at zero field at the AΓ-SPT_β transition, thereby obscuring it. We note that, the quantum critical points (QCPs) are immediately noticeable in the entanglement spectra.

Moving to the right within the AΓ phase which encompass the point Γ = 1, the transition to the FK phase from the AΓ phase occurs at ϕ = 0.88271(5)π, which is clearly visible in the entanglement spectra as well as ${\chi }_{\phi }^{e}$, as shown in in Fig. 3. The transition between the RS${}_{{U}_{6}}$ and FK phases at ϕ = π is less apparent in the structure of the entanglement spectrum in Fig. 3a, b. However, from the results for ${\chi }_{\phi }^{e}$ in Fig. 3c this transition is immediately visible, and it was also noted that the second derivative of λ₁ clearly detects the transition³⁰.

Similar to the AK phase, none of the phases SPT_α, SPT_β, AΓ and FK has any long-range magnetic ordering, nor is there any indication of nematic (quadropolar) or chiral ordering. As discussed in Supplementary Note 2, all four phases are gapped with a finite sizeable correlation length. The difference between them is captured in the entanglement spectrum. In the SPT_α, SPT_β and FK (and AK) phases, the entanglement spectrum have all entries doubled when considering ρ_B. For the AΓ-phase, the same applies to the spectrum for ρ_A. Since an entanglement spectrum where all eigenvalues have degeneracy larger than one is a signature of a topological non-trivial phase, in the next section, we investigate the projective symmetry analysis to confirm their non-trivial topology prior to presenting associated non-local SOPs.

Projective symmetry analysis

With the SPT_α, SPT_β, AΓ and FK phases as potential SPT phases, it is of considerable interest to investigate the projective representations^{33,37,38,39,40,41,42}U, of a site symmetry ${{{\mathcal{R}}}}$. This has been done for S = 1 chains^{37,38,39,40,41,42} and ladders⁴³ and also for S = 1/2 ladders^{41,44,45,46,47}. The matrices U can be obtained from the generalized transfer matrices in iDMRG calculations as described in ref. ⁴⁸.

For the KΓ ladder, it is a significant simplification to consider the transformed model obtained after applying the U₆ transformation. This transformation maps the original Hamiltonian H_KΓ to the transformed ladder, denoted by ${H}_{{{{\rm{K}}}}\Gamma }^{{U}_{6}}$. Under the U₆ transformation the x, y and z-bonds of the KΓ-ladder are transformed into anisotropic Heisenberg bonds, ${x}^{{\prime} }$, ${y}^{{\prime} }$ and ${z}^{{\prime} }$ in the following manner:

$$\begin{array}{rcl}{x}^{{\prime} }:-K{S}_{i}^{x}{S}_{j}^{x}-\Gamma ({S}_{i}^{y}{S}_{j}^{y}+{S}_{i}^{z}{S}_{j}^{z})\\ {y}^{{\prime} }:-K{S}_{i}^{y}{S}_{j}^{y}-\Gamma ({S}_{i}^{x}{S}_{j}^{x}+{S}_{i}^{z}{S}_{j}^{z})\\ {z}^{{\prime} }:-K{S}_{i}^{z}{S}_{j}^{z}-\Gamma ({S}_{i}^{x}{S}_{j}^{x}+{S}_{i}^{y}{S}_{j}^{y})\end{array}$$

(3)

In the transformed ladder ${H}_{{{{\rm{K}}}}\Gamma }^{{U}_{6}}$ model, the definition of ${x}^{{\prime} }$, ${y}^{{\prime} }$ and ${z}^{{\prime} }$-bonds, can pictorially be represented as shown in Fig. 4a. See Supplementary Note 1. We will consider two different open boundary conditions (OBC) as shown in Fig. 4. The unit cells depicted in panel (a) are referred to as regular unit cells with regular OBC, while the ones in panel (b) are referred to as slanted unit cells with slanted OBC.

Fig. 4: Unit cells of the U₆ transformed ladder.

Two unit cells of ${H}_{{{{\rm{K}}}}\Gamma }^{{U}_{6}}$, the KΓ ladder after the U₆ local transformation with N = 6n sites. a Two regular unit cells with regular open boundary conditions (OBC). b Two slanted unit cells with slanted OBC.

In order to understand the symmetries of the matrix product state (MPS) wave-function, it is useful to write it in the canonical form^49,50,51,52:

$$\left\vert \Psi \right\rangle =\mathop{\sum}\limits_{{j}_{1},\ldots ,{j}_{N}}{M}_{{j}_{1}}^{[1]}{\Lambda }^{[2]}{M}_{{j}_{2}}^{[2]}\ldots {\Lambda }^{[N]}{M}_{{j}_{N}}^{[N]}\left\vert {j}_{1},\ldots ,{j}_{N}\right\rangle ,$$

(4)

where the ${M}_{{j}_{n}}^{[n]}$ are complex matrices and the M^[n], real, positive, square diagonal matrices. In the iDMRG formulation, the set of matrices on any unit cell becomes the same ${M}_{j}^{[n]}$ = M_j, M^[n] = M for all n, although they may vary within the unit cell. For the translationally invariant state, it can be shown^34,53 that for any (site) symmetry operation g, represented in the spin basis by the unitary matrix, ${\Sigma }_{j{j}^{{\prime} }}(g)$, the M_j matrices must transform as^34,48:

$$\mathop{\sum}\limits_{{j}^{{\prime} }}{\Sigma }_{j{j}^{{\prime} }}(g){M}_{{j}^{{\prime} }}={e}^{i\theta }{U}^{{\dagger} }(g){M}_{j}U(g),$$

(5)

where the unitary matrix U(g) commutes with the M matrices and e^iθ is a phase factor. With D denoting the bond dimension, the U matrices form a D-dimensional projective representation of the symmetry group of the wave-function, and they can be determined from the unique eigenvector of the generalized transfer matrix^34,48 with eigenvalue ∣λ∣ = 1, where the generalized transfer matrix is defined as

$${T}_{\alpha {\alpha }^{{\prime} };\beta {\beta }^{{\prime} }}^{\Sigma }=\mathop{\sum}\limits_{j}\left(\mathop{\sum}\limits_{{j}^{{\prime} }}{\Sigma }_{j{j}^{{\prime} }}{M}_{{j}^{{\prime} },\alpha \beta }\right){({M}_{j,{\alpha }^{{\prime} }{\beta }^{{\prime} }})}^{* }{\Lambda }_{\beta }{\Lambda }_{{\beta }^{{\prime} }}$$

(6)

If the largest eigenvalue is ∣λ∣ < 1, the symmetry is not a property of the state being considered. The projective representation is reflected in the fact that if Σ(g)Σ(h)=Σ(gh), then

$$U(g)U(h)={e}^{i\phi (g,h)}U(gh),$$

(7)

where the phases ϕ(g, h) are characteristic of the topological phase.

Let us now consider the site symmetries, ${{{{\mathcal{R}}}}}^{x}$ and ${{{{\mathcal{R}}}}}^{y}$ defined as

$${{{{\mathcal{R}}}}}_{l}^{x}={e}^{i\pi {S}_{l}^{x}},\quad {{{{\mathcal{R}}}}}_{l}^{y}={e}^{i\pi {S}_{l}^{y}},\quad {{{{\mathcal{R}}}}}_{l}^{z}={e}^{i\pi {S}_{l}^{z}}.$$

(8)

The Hamiltonian ${H}_{{{{\rm{K}}}}\Gamma }^{{U}_{6}}$ is invariant under the operators ${\prod }_{l}{{{{\mathcal{R}}}}}_{l}^{\gamma }$, with γ = x, y, z with distinct quantum numbers for the low-lying states. If these symmetries are respected, their representations can differ by a phase, $\phi ({{{{\mathcal{R}}}}}_{x},{{{{\mathcal{R}}}}}_{y})$ that must obey ${e}^{i\phi ({{{{\mathcal{R}}}}}_{x},{{{{\mathcal{R}}}}}_{y})}$= ± 1:

$$U({{{{\mathcal{R}}}}}^{x})U({{{{\mathcal{R}}}}}^{y})=\pm U({{{{\mathcal{R}}}}}^{y})U({{{{\mathcal{R}}}}}^{x}).$$

(9)

Furthermore, the non-trivial value ${e}^{i\phi ({{{{\mathcal{R}}}}}_{x},{{{{\mathcal{R}}}}}_{y})}$ = −1 can only occur if all eigenvalues of the entanglement spectrum are at least twice degenerate³⁴. The phase factor can then be isolated by defining⁴⁸:

$${{{{\mathcal{O}}}}}_{{{{{\rm{Z}}}}}_{2}\times {{{{\rm{Z}}}}}_{{{{\rm{2}}}}}}\equiv \frac{1}{D}{{{\rm{Tr}}}}\left(U({{{{\mathcal{R}}}}}^{x})U({{{{\mathcal{R}}}}}^{y}){U}^{{\dagger} }({{{{\mathcal{R}}}}}^{x}){U}^{{\dagger} }({{{{\mathcal{R}}}}}^{y})\right),$$

(10)

with D the bond dimension, with similar definitions for other pairs of operators in Eq. (8). In the above definitions it is understood that the transformations are applied throughout the lattice and in order to obtain the matrices U, generalized transfer matrices representing the relevant unit cell has to be considered.

The site symmetries, ${{{{\mathcal{R}}}}}^{x}$, ${{{{\mathcal{R}}}}}^{y}$ and ${{{{\mathcal{R}}}}}^{z}$, forming the dihedral group, D₂, are respected by both of the unit cells in Fig. 4a, b. Using generalized transfer matrices obtained from unit cells of the shape shown in Fig. 4a when studying the AΓ-phase and of the slanted shape shown in Fig. 4b when studying the FK, SPT_α, and SPT_β phases, we obtain

$${{{{\mathcal{O}}}}}_{{{{{\rm{Z}}}}}_{2}\times {{{{\rm{Z}}}}}_{{{{\rm{2}}}}}}=-1,$$

(11)

for the SPT_α, SPT_β, AΓ, and FK phases.

Similar analysis can be made for the time-reversal $({{{\rm{TR}}}})$ symmetry, defined by ${M}_{j}\to {\sum }_{{j}^{{\prime} }}{\left[{e}^{i\pi {S}^{y}}\right]}_{j{j}^{{\prime} }}{M}_{{j}^{{\prime} }}^{\star }$, with ⋆ denoting complex conjugation. In this case, it can be established that³⁴${U}_{{{{\rm{TR}}}}}{U}_{{{{\rm{TR}}}}}^{\star }$ = ${e}^{i\phi ({{{\rm{TR}}}},{{{\rm{TR}}}})}{\mathbb{1}}$ where the phase $\phi ({{{\rm{TR}}}},{{{\rm{TR}}}})$ cannot be absorbed into the definition of ${U}_{{{{\rm{TR}}}}}$. One should note that for most other symmetries, with the notable exception of inversion, similar considerations will lead to U² = ${e}^{i\phi }{\mathbb{1}}$ in which case the phase ϕ in fact can be absorbed into the definition of U. For instance, this is the case for $U\left({{{{\mathcal{R}}}}}^{\alpha }\right)$ discussed above. However, for time reversal, the phase factor ${e}^{i\phi ({{{\rm{TR}}}},{{{\rm{TR}}}})}$ can directly be extracted by defining⁴⁸:

$${{{{\mathcal{O}}}}}_{{{{\rm{TR}}}}}\equiv \frac{1}{D}{{{\rm{Tr}}}}\left({U}_{{{{\rm{TR}}}}}{U}_{{{{\rm{TR}}}}}^{\star }\right),$$

(12)

and again one finds that $\phi ({{{\rm{TR}}}},{{{\rm{TR}}}})$ = 0,π, so that ${{{{\mathcal{O}}}}}_{{{{\rm{TR}}}}}$ = ± 1. As an example, for the S = 1 Heisenberg spin chain in the Haldane phase it is known that ${{{{\mathcal{O}}}}}_{{{{{\rm{{Z}}}_{2}\times Z}}}_{{{{\rm{2}}}}}}$ = −1, ${{{{\mathcal{O}}}}}_{{{{\rm{TR}}}}}$= −1^34,35.

Using generalized transfer matrices obtained from unit cells of the shape shown in Fig. 4a for the AΓ-phase and of the slanted shape for the FK, SPT_α, and SPT_β phases, we obtain

$${{{{\mathcal{O}}}}}_{{{{\rm{TR}}}}}=-1,$$

(13)

consistent with the presence of a doubled entanglement spectrum in all cases. As was the case for ${{{{\mathcal{O}}}}}_{{{{{\rm{{Z}}}_{2}\times Z}}}_{{{{\rm{2}}}}}}$, if the unit cells are interchanged, one finds instead ${{{{\mathcal{O}}}}}_{{{{\rm{TR}}}}}$ = 1.

A summary of our results from the projective analysis are provided in Table 1, negative values indicate that the state transforms non-trivially. For all 4 phases, it is seen that a unit cell can be chosen for which the state transforms non-trivially under both the ${{{\rm{TR}}}}$ and ${{{{\mathcal{O}}}}}_{{{{{\rm{{Z}}}_{2}\times Z}}}_{{{{\rm{2}}}}}}$ symmetries.

Table 1 Summary of projective analysis.

Based on the above analysis of entanglement spectrum and projective symmetry, we conclude that AΓ phase is an SPT phase. It is then important to further identify its SOP that differentiates this phase from the other disordered phases.

Twice hidden string order

To establish a non-local string order parameter (SOP) characterizing the AΓ phase, we need to exploit a non-local unitary transformation that maps the original Hamiltonian with OBC to a new Hamiltonian that exhibits a local long-range order^54,55,56. We found that it is difficult to identify such non-local transformation starting from the original H_KΓ, but it can be achieved by first applying the U₆ transformation to arrive at ${H}_{{{{\rm{K}}}}\Gamma }^{{U}_{6}}$. It is then possible to define a non-local unitary transformation W, mapping ${H}_{{{{\rm{K}}}}\Gamma }^{{U}_{6}}$ to a new local Hamiltonian. We denote the resulting Hamiltonian, where four-spin terms appear, by H_KQ. For the parameters relevant for the AΓ-phase, H_KQ exhibits long-range order in the spin-spin correlation functions, corresponding to a local order parameter. Due to the application of two separate unitary transformations, one might consider the resulting order to be twice hidden.

The non-local unitary operator W for a N-site ladder with OBC that maps ${H}_{{{{\rm{K}}}}\Gamma }^{{U}_{6}}$ to H_KQ takes the following form

$$W=\mathop{\prod}\limits_{\begin{array}{c}j+1\, < \,k\\ j\,{{{\rm{odd}}}},\,k\,{{{\rm{odd}}}}\\ j=1,\ldots N-3\\ k=3,\ldots N-1\end{array}}w(j,k).$$

(14)

With the individual w(j, k) given as follows:

$$w(j,k)={e}^{i\pi ({S}_{j}^{y}+{S}_{j+1}^{y})\cdot ({S}_{k}^{z}+{S}_{k+1}^{z})},$$

(15)

and W^† = W. The OBC are here crucial for the mapping to be exact. Evidently, all w(j, k) are unitary and therefore also W, and $\left[w(j,k),w(l,m)\right]=0\,\,\forall \,j,k,l,m$. Note that this is a different labeling of the unitary operator introduced in Refs. ^23,30. It can also be shown that other combinations of the spin operators S^α appearing in Eq. (15) lead to equivalent unitary operators, for instance, ${e}^{i\pi ({S}_{j}^{y}+{S}_{j+1}^{y})\cdot ({S}_{k}^{x}+{S}_{k+1}^{x})}$ is another valid choice. However, the specific choice made in Eq. (15) will influence the type of ordering that is observed in H_KQ, as well as the specific form of H_KQ. Schematically, the transformations can be viewed as shown in Fig. 5. The detailed form of H_KQ after the W transformation on ${H}_{{{{\rm{K}}}}\Gamma }^{{U}_{6}}$ is presented in Supplementary Note 4.

Fig. 5: Pictorial view of KG transformations.

The two transformations from the original KΓ model to ${H}_{{{{\rm{K}}}}\Gamma }^{{U}_{6}}$, and the subsequent transformation to H_KQ are sketched. The type of order parameter for each Hamiltonian is indicated.

We first discuss the ordering in H_KQ where we denote the spins by S^″α, where the double prime represents the two transformations of spin from the original Hamiltonian. It is then convenient to study correlation functions of the form 4〈S^″α(0)S^″α(r)〉 along the legs of the ladder with r measured in lattice spacings along the leg. To avoid boundary effects, r = 0 is usually taken to correspond to a site in the bulk of the chain. In Fig. 6 we show results for 4〈S^″y(0)S^″y(r)〉 starting from site 47 with ϕ = 0.85π. Long-range order is clearly present. Similar results can be obtained for the other leg of the ladder as well as for 4〈S^″z(0)S^″z(r)〉. However, due to the choice of spin operators in the definition of w in Eq. (15), there is no ordering in 4〈S^″x(0)S^″x(r)〉.

Fig. 6: Spin correlations in the H_KQ model at ϕ = 0.85π.

DMRG results with N = 384 for the correlation function 4〈S^″y(0)S^″y(r)〉 versus distance, r along one leg of the ladder. r = 0 corresponds to site 47. Results are for H_KQ with ϕ = 0.85π and r is measured in lattice spacings along the leg.

Using the inverse of the non-local unitary operator W from Eq. (14) the above results for 4〈S^″y(0)S^″y(r)〉 is reproduced as a non-local string order correlation function in ${H}_{{{{\rm{K}}}}\Gamma }^{{U}_{6}}$ where we denote the spin variables by ${S}^{{\prime} \alpha }$ = ${\sigma }^{{\prime} \alpha }/2$. We find $\langle {S}_{1}^{{\prime\prime} y}{S}_{r+1}^{{\prime\prime} y}\rangle$ is given by

$$\begin{array}{ll}{{{{\mathcal{O}}}}}^{y}(r)\,=\,4\langle {S}_{1}^{{\prime\prime} y}{S}_{r+1}^{{\prime\prime} y}\rangle ={(-1)}^{\lfloor r/2\rfloor }\\ \qquad\qquad\times \,\left\{\begin{array}{ll}\left\langle {\sigma }_{2}^{{\prime} y}\left(\mathop{\prod }\limits_{k=3}^{r}{\sigma }_{k}^{{\prime} y}\right){\sigma }_{r+1}^{{\prime} y}\right\rangle \quad &r\,{{{\rm{even}}}}\\ \left\langle {\sigma }_{2}^{{\prime} y}\left(\mathop{\prod }\limits_{k=3}^{r-1}{\sigma }_{k}^{{\prime} y}\right){\sigma }_{r+1}^{{\prime} y}\right\rangle \quad &r\,{{{\rm{odd}}}}\end{array}\right.\end{array}$$

(16)

Note that, to fully reproduce the results for H_KQ shown in Fig. 6 with the string correlation function in (16) for ${H}_{{{{\rm{K}}}}\Gamma }^{{U}_{6}}$, a relabeling of the sites needs to be done that we have skipped for clarity.

We can now apply the inverse U₆ transformation (see Supplementary Note 1) to the above expressions for ${{{{\mathcal{O}}}}}^{y}(r)$ to determine the string order correlation functions that define the AΓ-phase in the original H_KΓ:

$${\Xi }^{y}(r)={U}_{6}^{-1}({{{{\mathcal{O}}}}}^{y}(r)),$$

(17)

with the ${U}_{6}^{-1}$ transformation detailed in Supplementary Note 1. The SOP in H_KΓ Hamiltonian is then given by

$${\Xi }^{y}=\mathop{\max }\limits_{r\to \infty }{\Xi }^{y}(r).$$

(18)

It is interesting to note that for a small r, for example r = 7, Ξ^y(r = 7) corresponds to the plaquette operator found in the pure Gamma model in the honeycomb lattice⁵⁷.

In Fig. 7 we show iDMRG results for Ξ^y (orange circles). The inset shows iDMRG results for Ξ^y(r) versus r at ϕ = 0.85π which at large r can be compared to the results for H_KQ shown in Fig. 6. Note that the results in Fig. 6 show the correlations along a leg, whereas the inset in Fig. 7 show results from both legs without including the sign and the relabeling of the sites. Due to the different methods used, there is not an exact equivalence for very small values of r.

Fig. 7: String order parameters in the AΓ and surrounding phases.

iDMRG results for the string order parameters in the H_KΓ and ${H}_{{{{\rm{K}}}}\Gamma }^{{U}_{6}}$ models, shown alongside DMRG results for H_KQ. Orange circles, Ξ^y for H_KΓ. The inset shows Ξ^y(r) versus r at ϕ = 0.85π for H_KΓ. Magenta circles, DMRG results for ${{{{\mathcal{O}}}}}_{{{{\rm{even}}}}}^{z}$ for H_KQ. Light blue circles, ${{{{\mathcal{P}}}}}_{{{{\rm{even}}}}}^{z}$ for ${H}_{{{{\rm{K}}}}\Gamma }^{{U}_{6}}$. Green circles, ${{{\mathcal{Z}}}}$ for ${H}_{{{{\rm{K}}}}\Gamma }^{{U}_{6}}$.

We emphasize that the appearance of a non-local SOP in the AΓ phase of H_KΓ is equivalent to the presence of long-range order in H_KQ. Hence, Ξ^y is non-zero throughout the AΓ-phase and goes to zero at the critical points delineating this phase. It is absent in the other disordered phases, SPT_α, SPT_β, and FK, and thus uniquely defines the AΓ phase.

With the identification of the AΓ-phase with regular long-range ordering in the H_KQ model, it is natural to ask if a regular (local) order parameter can also be identified for the SPT_α- and SPT_β-phases in the H_KQ. However, all local order parameters that we have investigated have not shown any ordering in the SPT_α- and SPT_β-phases for H_KQ. Extending the string-order correlation function defined for the S = 1 Haldane chain, a heuristic string-order correlation function has been proposed for S = 1/2 ladders by pairing two S = 1/2^58,59. Following the numbering of Fig. 1a, if ${\tau }_{i}^{\alpha }={S}_{2i}^{\alpha }+{S}_{2i+1}^{\alpha }$ are the sum of two diagonally situated spins, one defines^58,59:

$${{{{\mathcal{O}}}}}_{{{{\rm{even}}}}}^{\alpha }(r)=\left\langle {\tau }_{i}^{\alpha }\exp \left(i\pi \mathop{\sum }\limits_{l=i+1}^{i+r-1}{\tau }_{l}^{\alpha }\right){\tau }_{i+r}^{\alpha }\right\rangle .$$

(19)

The associated SOP is non-zero in the phase surrounding ϕ = 0 in H_KQ³⁰, corresponding to the AFM Kitaev (AK) phase in H_KΓ. The magenta points in Fig. 7 show our results for ${{{{\mathcal{O}}}}}_{{{{\rm{even}}}}}^{z}$ for the H_KQ model, which clearly is non-zero in the SPT_α- and SPT_β-phases. This is consistent with the nonexistence of a local order parameter in these two phases for H_KQ. We note that, due to the heuristic nature of ${{{{\mathcal{O}}}}}_{{{{\rm{even}}}}}^{z}$, it is not clear how to associate it with a local order in a related model. Since H_KΓ and ${H}_{{{{\rm{K}}}}\Gamma }^{{U}_{6}}$ are related by a local unitary transformation, any ordering in the SPT_α- and SPT_β-phases in either model would immediately be apparent in both, and we have not observed any for either model.

Building on the above results for ${{{{\mathcal{O}}}}}_{{{{\rm{even}}}}}^{z}$ for H_KQ we propose a closely related heuristic string order correlation function for ${H}_{{{{\rm{K}}}}\Gamma }^{{U}_{6}}$ in the following way: Define ${l}_{i}^{\alpha }={S}_{2i}^{\alpha }-{S}_{2i+1}^{\alpha }$ as the difference of two diagonally situated spins, following the numbering from Fig. 1a. We then have

$${{{{\mathcal{P}}}}}_{{{{\rm{even}}}}}^{\alpha }(r)=\left\langle {l}_{i}^{\alpha }\exp \left(i\pi \mathop{\sum }\limits_{l=i+1}^{i+r-1}{\tau }_{l}^{\alpha }\right){l}_{i+r}^{\alpha }\right\rangle ,$$

(20)

with ${\tau }_{i}^{\alpha }={S}_{i,1}^{\alpha }+{S}_{i+1,2}^{\alpha }$ as above. Results for ${{{{\mathcal{P}}}}}_{{{{\rm{even}}}}}^{z}$ versus ϕ obtained from iDMRG calculations with ${H}_{{{{\rm{K}}}}\Gamma }^{{U}_{6}}$ are shown in Fig. 7 as the light blue points. The FK, SPT_α and SPT_β-phases are clearly defined by a non-zero ${{{{\mathcal{P}}}}}_{{{{\rm{even}}}}}$. By applying ${U}_{6}^{-1}$ to the definition of ${{{{\mathcal{P}}}}}_{{{{\rm{even}}}}}$ it is straightforward to perform similar calculations using H_KΓ by evaluating ${U}_{6}^{-1}({{{{\mathcal{P}}}}}_{{{{\rm{even}}}}}^{\alpha }(r))$. Even though the definition ${{{{\mathcal{P}}}}}_{{{{\rm{even}}}}}$ is heuristic, we interpret this result as a verification of the SPT nature of the FK, SPT_α- and SPT_β-phases.

For the FK-phase, it is also instructive to consider an even simpler heuristic string order correlation function, ${{{\mathcal{Z}}}}$ defined as follows²²:

$${{{\mathcal{Z}}}}=\left\langle \mathop{\prod }\limits_{i}^{i+r}{\sigma }_{i}^{z}\right\rangle$$

(21)

We consider this string correlator for ${H}_{{{{\rm{K}}}}\Gamma }^{{U}_{6}}$ or equivalently to H_KΓ through application of ${U}_{6}^{-1}$. Results for ${{{\mathcal{Z}}}}$ are shown in Fig. 7 versus ϕ (green points) as obtained from iDMRG calculations. Throughout the FK-phase ${{{\mathcal{Z}}}}$ is almost identical to 1, dropping to zero at the transition to the AΓ-phase. However, we note that ${{{\mathcal{Z}}}}$ remain sizable throughout much of the RS${}_{{U}_{6}}$-phase, reflecting its heuristic nature.

Edge states and response to magnetic field

Another signature of SPT phases is the presence of edge states under OBC related to a ground state degeneracy. For the SPT phases in the KΓ ladder, it is clear from the degeneracy of the entanglement spectrum that we need to consider different shapes of clusters (regular vs. slanted OBCs) for the different SPT phases. In this section, we exclusively study the original Hamiltonian H_KΓ and do not consider ${H}_{{{{\rm{K}}}}\Gamma }^{{U}_{6}}$ nor H_KQ. For AΓ phase, we use N = 4n with the regular OBC and for the remaining SPT phases, we use the slanted OBC with N = 4n + 2 in order to have an equal number of the different bond types.

We first demonstrate the presence of edge states in the SPT phases. For the AΓ phase, results for the 16 lowest states with the regular OBC at ϕ = 0.85π are obtained using ED (see Supplementary Fig. 7). Four low-lying states below the gap are clearly present. With increasing N, these 4 states quickly become degenerate while the gap stabilizes at a finite value. Similar results can be obtained for the AK and FK phase using the slanted cluster, and a degeneracy of 4 is also observed for this case (see Supplementary Fig. 8). For the SPT_α- and SPT_β-phases, it has not been possible to produce reliable results in the same manner, likely due to the significantly larger correlation lengths.

To understand the nature of the edge-states, let us explore how the AK, FK, and AΓ phases respond to an external magnetic field. An external magnetic field introduces an additional term in the Hamiltonian of the form ${H}^{{\prime} }={g}_{L}{\mu }_{B}{{{\bf{B}}}}\cdot {{{{\bf{S}}}}}_{{{{\rm{tot}}}}}$, where S_tot = ∑_iS_i, g_L is the Landé factor and μ_B the Bohr magneton. Following ref. ³⁸ we denote the 4 states ψ₁, ψ₂, ψ₃ and ψ₄ and consider S_tot,α in this four-fold degenerate space by defining

$${S}_{{{{\rm{tot}}}},\alpha }^{\gamma \beta }=\langle {\psi }_{\gamma }| {S}_{{{{\rm{tot}}}},\alpha }| {\psi }_{\beta }\rangle ,\,\,\gamma ,\beta =1,2,3,4.$$

(22)

Here, the components of the total spin S_tot,α are usually taken to be identical to x, y, z but given the underlying honeycomb structure we shall find it useful to instead consider α = $\hat{{{{\bf{a}}}}}$, $\hat{{{{\bf{b}}}}}$, $\hat{{{{\bf{c}}}}}$ corresponding to the three directions [11 − 2], [1 − 10] and [111] that correspond to the perpendicular and parallel to the z-bond, and perpendicular to the plane of the honeycomb (or n-leg brick-wall), respectively. One then finds that the eigenvalues of the matrices $({S}_{{{{\rm{tot}}}},\alpha }^{\gamma \beta })$ for four corresponding degenerate states are simply given by (s_α, − s_α, 0, 0).

Our ED results for s_α are shown in Fig. 8. For the AK phase at ϕ = −0.08π, shown in Fig. 8a, we find s_a ~ 3/4, s_b ~ 1/2 and s_c ≃ 0.85. Similarly, for the FK-phase at ϕ = 0.95π, shown in Fig. 8b, we find s_a = 1, s_b ~ 2 and s_c = 3/4. In both cases, we expect some variation in the values of s_α as ϕ is tuned. We note that for both the FK and AΓ phase, the values of s_α quickly saturate at a small, finite value as N is increased. This is indicative of excitations localized at the edges as opposed to an actual magnetically ordered ground-state which should show s_α continually growing with N. A calculation of s_α for the SPT_α and SPT_β-phases do not yield clear results for the range of N available with ED, as discussed in Supplementary Note 6. In a realistic experimental setting, the presence of impurities will always create finite open segments of ladders, with a resulting Curie-law behavior. The response to an applied magnetic field is in that case highly anisotropic and the low temperature Curie-law response should show a strong directional dependence^38,41,43 with χ_a(T), χ_b(T) and χ_c(T) clearly distinguishable. The results for the AΓ-phase, shown in Fig. 8c, are even more intriguing. They are not only more anisotropic, but only s_b is non-zero, approaching a value close to 2 at ϕ = 0.85π. At a slightly different point in the AΓ-phase with ϕ = 0.8π we instead find s_b ≃ 4/3 but again only s_b is non-zero. This implies that the phase does not respond at all to a field applied along the $\hat{{{{\bf{a}}}}}$ and $\hat{{{{\bf{c}}}}}$ directions, effectively g_a, g_c ≃ 0.

Fig. 8: Eigenvalues of total spin in the AK, FK and AΓ phases.

The eigenvalues s_α = s_a,b,c of the total spin S_α = ${\sum }_{i}{S}_{i}^{\alpha }$ in the four degenerate ground-states, in zero field. Results are from ED. a AK-phase at ϕ = −0.08π with N = 14, 18…30 b FK-phase at ϕ = 0.95π with N = 14, 18…30, c AΓ-phase at ϕ = 0.85π with N = 12, 14, …30 (filled symbols).

To gain a clearer picture of how the ladder in the AΓ phase responds to a magnetic field applied along the $\hat{{{{\bf{b}}}}}$-direction, we have performed ED calculations in the presence of a small field in the $\hat{{{{\bf{b}}}}}$-direction. The resulting site dependent magnetization $\langle {S}_{i}^{a,b,c}\rangle$ can then easily be obtained for the $\hat{{{{\bf{a}}}}},\hat{{{{\bf{b}}}}}$ and $\hat{{{{\bf{c}}}}}$-directions. Our results are shown in Fig. 9 for N = 24 as obtained from the ground-state with a small field in the $\hat{{{{\bf{b}}}}}$-direction of 0.002 at ϕ/π = 0.85 in the AΓ-phase. The green, blue and cyan colors represent positive expectation values, whereas orange, red and pink colors indicate negative values, with the size of the points proportional to the expectation value. The response along the $\hat{{{{\bf{a}}}}}$, $\hat{{{{\bf{c}}}}}$ directions is completely symmetric in the positive and negative directions, yielding ${\sum }_{i}{S}_{i}^{a,c}$ = 0. However, along the $\hat{{{{\bf{b}}}}}$-direction the response is much larger, and we clearly find ${\sum }_{i}{S}_{i}^{b}\,\ne$ 0 consistent with the results shown in the SI (See Fig. S9) (c). With the field along the $\hat{{{{\bf{b}}}}}$-direction, sizeable excitations are visible at both ends of the ladder.

Fig. 9: Ground state magnetization in an infinitesimal field in the AΓ phase.

ED results with N = 24 for $\langle {S}_{i}^{a,b,c}\rangle$ in the lowest state in the AG phase at ϕ/π = 0.85 with a field in the $\hat{{{{\bf{b}}}}}$-direction of strength 0.002. The green, blue and cyan colors indicate positive values while orange, red, and pink indicate negative values. The size of the circles are proportional to the value of $\langle {S}_{i}^{a,b,c}\rangle$. a$\langle {S}_{i}^{a}\rangle$, b$\langle {S}_{i}^{b}\rangle$, c$\langle {S}_{i}^{c}\rangle$.

The response of the edge states to a magnetic field correlates with the lifting (or absence of lifting) of the degeneracy in the entanglement spectrum. As previously discussed, non-trivial indices in the projective analysis can only arise if the degeneracy of all states in the entanglement spectrum (ES) is larger than one. This implies that if a finite strength of the perturbation is needed to remove the degeneracy, then a phase transition does not occur until that strength is reached and the phase persists till that point. On the other hand, if the degeneracy is lifted for any non-zero strength of the perturbation, the symmetry protection is broken without an associated phase transition. We can then investigate the response of the degeneracy of the ES to a magnetic field in the $\hat{{{{\bf{a}}}}}$, $\hat{{{{\bf{b}}}}}$ and $\hat{{{{\bf{c}}}}}$ directions. This is shown in Fig. 10 for the AK, FK, and AΓ phases. For simplicity, we focus exclusively on the difference in the two largest eigenvalues Δ = λ₁ − λ₂, and in each case we employ the reduced density matrix that has a two-fold degeneracy at zero field.

Fig. 10: The Schmidt gap in the AK, FK and AΓ phases.

The splitting of the entanglement spectrum Δ = λ₁ − λ₂ from iDMRG calculations with a field in the $\hat{{{{\bf{a}}}}}$, $\hat{{{{\bf{b}}}}}$ and $\hat{{{{\bf{c}}}}}$ directions. a AK-phase at ϕ = −0.08π using ρ_B. b FK-phase at ϕ = 0.95π using ρ_B. c AΓ-phase at ϕ = 0.85π using ρ_A.

As can be seen in Fig. 10a the degeneracy in the ES for the AK phase is immediately lifted by a field in any of the three directions which correlates with the response of the edge-states (See Fig. S9a in the SI). However, the response is rather weak, and a relatively large field has to be applied to see a significant splitting. For the FK phase, we have a similar effect as shown in Fig. 10b, but in this case the response is much stronger. However, for the AΓ phase, where we show results in Fig. 10c, it is clear that a field in the $\hat{{{{\bf{a}}}}}$ direction of around h_a = 0.06 is needed to lift the degeneracy of the ES. For a field applied in the $\hat{{{{\bf{c}}}}}$ direction, a significantly larger field is needed. On the other hand, a field in the $\hat{{{{\bf{b}}}}}$ direction immediately induces a large splitting in the entanglement spectrum even for infinitesimal field strengths. Since the degeneracy remains intact in the $\hat{{{{\bf{a}}}}}$ and $\hat{{{{\bf{c}}}}}$ directions, we conclude that the SPT character of the AΓ phase persists with respect to a field applied in the $\hat{{{{\bf{a}}}}}$ and $\hat{{{{\bf{c}}}}}$ directions.

The AΓ phase is then protected by the product of time-reversal (${{{\rm{TR}}}}$) and π rotation around the $\hat{{{{\bf{b}}}}}$-axis (${{{{\mathcal{R}}}}}_{b}$), ${{{\rm{TR}}}}\times {{{{\mathcal{R}}}}}_{b}$, the only remaining symmetry⁶⁰ when the field is in the $\hat{{{{\bf{a}}}}}\hat{{{{\bf{c}}}}}$ plane, but broken when it is along the $\hat{b}$-axis. Hence, if a field is applied in the $\hat{{{{\bf{a}}}}}\hat{{{{\bf{c}}}}}$ plane, a transition to the trivial polarized state can only occur at finite field strengths with potentially other phases intervening before the polarized state is encountered. Several such transitions were observed for the AΓ phase (denoted by KΓSL) for a field in the $\hat{{{{\bf{c}}}}}$ direction^18,30. On the other hand, the FK phase is not protected by the ${{{\rm{TR}}}}\times {{{{\mathcal{R}}}}}_{b}$ symmetry, and if a field is applied in the $\hat{{{{\bf{c}}}}}$ direction the ES degeneracy is lost as shown in Fig. 10b. However, the field induced FK phase can still be distinguished from the polarized state.

Discussion

Our initial inquiry in this paper pertains to the nature of the AΓ phase and whether there exists a defining quantity for its characterization. For example, the Kitaev phases (AK and FK) in the ladder display the character of SPT phases. It is likely that AΓ is another SPT phase. If so, we expect all the signatures of the SPT such as the degeneracy of the entanglement spectrum, ground state degeneracy under OBC, and the presence of a SOP. Using the iDMRG, DMRG, and ED techniques, we indeed found that the entanglement spectrum is degenerate and there exists four-fold ground state degeneracy under the regular OBC in the AΓ phase. It is interesting to note that the same results were obtained for the Kitaev phases, AK and FK, but under the slanted OBC.

Despite such clear signatures of the SPT, determining the corresponding SOP in the AΓ phase has been challenging. We found that the string order correlation function is related to ordinary local order in a regular correlation function in a model H_KQ obtained only after two consecutive unitary transformations. Hence, we term this order as ‘twice’ hidden.

To understand the symmetry that protects the AΓ phase, we also investigated the effects of the external magnetic field. From the magnetic field response, we noted that the AΓ phase is completely inert to the magnetic field when the field is applied in the $\hat{a}$ and $\hat{c}$ directions, which correspond to perpendicular to the z-bond and the ladder plane, respectively. Accordingly, the entanglement spectrum degeneracy remains intact when the field is applied in the $\hat{a}$ and $\hat{c}$ directions. This is in contrast to the effect of applying the magnetic field along the $\hat{b}$ direction, i.e., parallel to the z-bond, which immediately lifts the degeneracy of the entanglement spectrum. We note that the product of ${{{\rm{TR}}}}$ and ${{{{\mathcal{R}}}}}_{b}$, ${{{\rm{TR}}}}\times {{{{\mathcal{R}}}}}_{b}$ symmetry, is preserved when the magnetic field is applied in the $\hat{{{{\bf{a}}}}}\hat{{{{\bf{c}}}}}$-plane, which is valid for the generic honeycomb Kitaev model beyond the ladder⁶⁰. Thus, we conclude that the AΓ is protected by the ${{{\rm{TR}}}}\times {{{{\mathcal{R}}}}}_{b}$ symmetry. Another intriguing implication from the magnetic field study is that the edge state in the AΓ phase is not a isotropic free spin-1/2 unlike the standard S = 1 Haldane SPT. They act like spinless modes under the field in $\hat{a}$ and $\hat{c}$ directions. Further studies are needed to fully understand the nature of the zero-energy modes at the boundary of the system with the regular OBC.

In the context of Kitaev materials, let us revisit our motivations for investigating the AΓ phase in the ladder model. As previously mentioned in the introduction, the majority of d⁵ Kitaev materials prominently feature FM Kitaev and AFM Γ interactions. However, ongoing debates persist regarding the specific phase that arises in this region. Several numerical studies have suggested the presence of magnetic disorder^{12,13,14,15,16,17,18,19}, while others have indicated a magnetically ordered phase, such as a zig-zag order²¹. Should a zig-zag order indeed be manifest in the 2D limit, we would anticipate observing the same ordering pattern in the 2-leg ladder, as the magnetic unit cell of the zig-zag can be captured in the ladder geometry. This is indeed the case for the Kitaev-Heisenberg ladder model, where the zig-zag, stripy, and FM ordered phases reported in the 2D honeycomb clusters are found in the ladder geometry²³.

Our findings have substantiated the presence of disordered state in the AΓ phase of the 2-leg ladder, categorizing it as a SPT phase characterized by a SOP. The 2D limit can be constructed by stacking the ladders, and one possibility of the resulting 2D phase is a stacked SPTs with edge modes known as a weak-SPT⁶¹. However, the coupling between the ladders may generate a new phase or critical point. It is interesting to note that the evolution from a stacked weak-SPT chains to a gapless critical point supporting edge modes that do not hybridize with bulk modes was reported in the extended anisotropic Kitaev model approaching from the dimer limit⁶². Our findings hint at the possibility that as the 2D limit is approached, the AΓ phase may become a 2D spin liquid, denoted as the KΓ spin liquid. However, we cannot rule out a possibility of large unit cell⁶³ or incommensurate²⁰ magnetic orders whose magnetic unit cells are beyond the ladder geometry, and a definitive resolution to this question remains a subject for future investigation.

Methods

Numerical Methods

We use a fully parallelized implementation⁶⁴ of the Lanczos algorithm to perform the exact diagonalizations (ED) of ladders with up to N = 30 using both open and periodic boundary conditions. In addition to exact diagonalizations, we use finite size density matrix renormalization group^{65,66,67,68,69,70} (DMRG) to study both the KΓ model, Eq. (1) and ${H}_{{{{\rm{K}}}}\Gamma }^{{U}_{6}}$ under both periodic (PBC) and open (OBC) boundary conditions, with the main part of our results obtained from the infinite DMRG^70,71 (iDMRG) variant of DMRG. The iDMRG calculations were performed with unit cells of 12, 24 or 60 sites. Typical precision for both DMRG and iDMRG are ϵ < 10⁻¹¹ with a bond dimension in excess of 1000.

Energy susceptibility

To determine the phase diagram we study the susceptibility derived from the ground-state energy per spin, e₀:

$${\chi }_{\phi }^{e}=-\frac{{\partial }^{2}{e}_{0}}{\partial {\phi }^{2}},$$

(23)

At a quantum critical point (QCP) it is known⁷² that, for a finite system of size N, the energy susceptibility diverges as

$${\chi }^{e} \sim {N}^{2/\nu -(d+z)}.$$

(24)

Here ν and z are the correlation and dynamical critical exponents and d is the dimension. Hence, χ^e only diverges at the phase transition if the critical exponent ν is smaller than 2/(d + z). For the present case we have d = 1 and if we assume z = 1, then a divergence is observed only if ν < 1.

Entanglement entropy and spectrum, Schmidt gap

When studying the ladder shown in Fig. 1a it is important to realize that there are different ways of partitioning the system in two partitions of size x and N − x. This is crucial when considering the bipartite von Neumann entanglement entropy, EE, as well as for the entanglement spectrum⁷³ of central importance for understanding topological properties^34,73,74,75. Both are obtained from the spectrum of the reduced density matrix, ρ_x, of either one of the two partitions. Here we focus on two specific partitions shown in Fig. 1a as the red and blue dashed lines. With the numbering in Fig. 1a, they correspond to either an odd (N/2-1, red) or even (N/2, blue) number of sites in the partitions. We refer to the density matrix derived from the former case with N/2-1 as ρ_B and to the latter case with N/2 as ρ_A. We mainly focus on the case where the number of sites in the partition is close to the mid-point, either N/2 − 1 (ρ_B) or N/2 (ρ_A). but when considering the bipartite entanglement entropy, we let the number in the partition vary but only consider an even number of sites in the sub system corresponding to moving the blue dashed line in Fig. 1a along the ladder. For a subsystem, A, of size x the entanglement entropy is defined by:

$$EE(x)=-{{{\rm{Tr}}}}{\rho }_{x}\ln {\rho }_{x}.$$

(25)

Our results for EE(x) can be found in Supplementary Note 2. The eigenvalues, l_α, of the reduced density matrix, ρ_x, correspond to the Schmidt decomposition, l_α = ${\lambda }_{\alpha }^{2}$ and thereby the entanglement spectrum⁷³, which then will depend on whether ρ_A or ρ_B from Fig. 1a is used. The Schmidt gap is then defined as λ₁ − λ₂.