Unusual excitations and double-peak specific heat in a bond-alternating spin-$1$ $K$-$\Gamma$ chain

One-dimensional gapped phases that avoid any symmetry breaking have drawn enduring attention. In this paper, we study such phases in a bond-alternating spin-1 $K$-$\Gamma$ chain built of a Kitaev ($K$) interaction and an off-diagonal $\Gamma$ term. In the case of isotropic bond strength, a Haldane phase, which resembles the ground state of a spin-$1$ Heisenberg chain, is identified in a wide region. A gapped Kitaev phase situated at dominant ferromagnetic and antiferromagnetic Kitaev limits is also found. The Kitaev phase has extremely short-range spin correlations and is characterized by finite $\mathbb{Z}_2$-valued quantities on bonds. Its lowest entanglement spectrum is unique, in contrast to the Haldane phase whose entanglement spectrum is doubly degenerate. In addition, the Kitaev phase shows a double-peak structure in the specific heat at two different temperatures. In the pure Kitaev limit, the two peaks are representative of the development of short-range spin correlation at $T_h \simeq 0.5680$ and the freezing of $\mathbb{Z}_2$ quantities at $T_l \simeq 0.0562$, respectively. By considering bond anisotropy, regions of Haldane phase and Kitaev phase are enlarged, accompanied by the emergence of dimerized phases and three distinct magnetically ordered states.


I. INTRODUCTION
The Kitaev honeycomb model [1], consisting of bonddependent Ising couplings of spin-1/2 degrees of freedom, is a rare example which not only is exactly solvable but also hosts a quantum spin liquid (QSL) ground state with fractionalized excitations, e.g., itinerant Majorana fermions and localized fluxes [2]. These excitations account for the double-peak specific heat anomaly at two different energy scales [3]. During the last decade, a large family of rare-earth magnets which could realize bond-directional interactions have garnered huge interest, providing avenues for the exploration of exotic phases of matter and emergent phenomena (for a review, see Refs. [4,5]). According to the Jackeli-Khaliullin mechanism [6], it was suggested that the Kitaev interaction with J eff = 1/2 moment could be realized in the 4d/5delectron honeycomb compounds by the interplay of spinorbit coupling and electron correlations. Recently, new scenarios for the Kitaev interaction in f -electron systems have also been proposed (see Ref. [7] and references therein). This theoretical progress as well as relevant material realizations, paves the way for hunting Kitaev QSLs; yet this is hindered by several essential non-Kitaev terms, such as Heisenberg coupling and symmetric offdiagonal Γ interaction [8,9].
Aligning with the efforts to study the ground state and thermodynamical properties in the spin-1/2 Kitaev honeycomb model [1-3, 10, 11], high-spin analogs have also generated much interest [12][13][14], inciting the materialization of Kitaev interaction in magnetic compounds * shijiehu@csrc.ac.cn † hykee@physics.utoronto.ca with S > 1/2 [15]. Recently, it has been suggested that strong spin-orbit coupling between anion sites together with strong Hund's coupling in the e g orbital might be a practical way to achieve the spin-1 Kitaev interaction in real materials [16]. In spite of there being no exact solution for higher-spin counterparts, local conserved Z 2 quantities could still be constructed [12], and several interesting phenomena of the spin-1/2 model, including the double-peak specific heat [17] and field-induced intermediate gapless QSL [18][19][20], could also be retained at least for the spin-1 case. Notwithstanding the bidimensionality of real materials, quantum spin chains also play vital roles in understanding peculiar quantum phenomena in two dimensions as they promote strong quantum fluctuations [21][22][23][24][25][26][27][28].
Over the past few decades, quantum spin chains have attracted broad attention for their ability to host unconventional quantum criticality [29,30] and topological phases [31,32]. The Haldane phase in the antiferromagnetic (AFM) spin-1 Heisenberg chain is a paramount example which falls beyond Landau's paradigm of symmetry breaking, and is now recognized as a symmetry-protected topological (SPT) phase protected by time-reversal symmetry, bond-centered inversion symmetry, and/or a dihedral group of π rotations about cubic axes [33,34]. Its ground state is unique under the periodic boundary condition (PBC), while it is four-fold degenerate under the open boundary condition (OBC) as a result of two spin-1/2 edge states [35]. Nevertheless, a novel "Kitaev" phase, the ground state of the spin-1 Kitaev spin chain [36,37], emerges as another interesting phase which is also gapped and has the same ground-state degeneracy pattern as that of the Haldane phase.
In this respect, a couple of attractive questions are naturally raised. First of all, although the Kitaev phase owns all the three symmetries that protect the Haldane phase [37], the origin of the edge states that contribute to the ground-state degeneracy is still unclear. To understand whether it is a SPT phase or not, the entanglement spectrum is a useful quantity to clarify the puzzle. Next, in the Kitaev phase, degeneracy of the lowest-lying excited state goes with the increase in system size, exhibiting a large density of states just above the excitation gap. This phenomenon has been demonstrated to trigger double-peak specific heat anomalies in several frustrated systems [38,39]. Hence, it is necessary to check the lowtemperature behavior of the specific heat. Finally, in the spin-1/2 analogy, the ferromagnetic (FM) Kitaev limit is known to be a multicritical point as a confluence of several topological quantum phase transitions (QPTs) [27]. By contrast, a spin-1 FM Kitaev point should survive against competing interactions, giving rise to a possible region of the Kitaev phase.
In this paper, we study the quantum phase diagram of a bond-alternating spin-1 K-Γ chain using the densitymatrix renormalization group (DMRG) method [40][41][42]. For this model, it is composed of two bond-directional frustrated interactions, allowing us to explore the rich phase diagram by tuning both the interaction intensity and the bond strength relatively. Throughout the phase diagram, we identify the Kitaev phase and the Haldane phase in the vicinity of the Kitaev and Γ limits, respectively. The natures of phases are revealed by excitation gaps, spin-spin correlations, the nonlocal string order parameter (SOP) [43], and Z 2 quantities on bonds. Near the FM Kitaev limit, there is a first-order Kitaev-Haldane QPT at nonzero Γ/|K|, while a magnetically ordered state intervenes in the AFM Kitaev region. We also calculate thermodynamic quantities (e.g., specific heat and thermal entropy) of the Kitaev phase using the transfer-matrix renormalization group (TMRG) method [44,45].
The remainder of the paper is organized as follows. In Sec. II we introduce the theoretical model and present the quantum phase diagram. In Sec. III we study the excitations of the Haldane phase and Kitaev phase in the isotropic K-Γ chain. Effects of anisotropic bond strength are studied in Sec. IV, with an emphasis on the Haldanedimer transition and three magnetically ordered states. In Sec. V, we confirm the double-peak specific heat in the Kitaev phase, and explain the physical origins of the two peaks. Finally, a brief conclusion is presented in Sec. VI. Further information about the specific heat in the spin-1/2 and spin-1 Kitaev chains is given in Appendices A and B.

II. MODEL AND METHOD
The K-Γ spin model is defined on the zigzag chain as illustrated in Fig. 1(a), where the spins sit on the edges of each bond. The full Hamiltonian is composed of two analogical terms, where L is the chain length and g x (g y ) is the strength of the odd (even) bond. The exchange part contains the Kitaev (K) interaction and the off-diagonal Γ interaction, which is given by Here, γ could be either x or y and it refers to the type of bond that connects spins i and j, see Fig. 1 where the bonds γ =x,ỹ, andz circularly, as depicted in Fig. 1(b). In light of Eq. (3), it can instantly be found that SU (2) symmetry recovers along two lines K = ±|Γ|. Here, x (red bonds) and y (green bonds) stand for the γ-index and bond widths indicate their strengths relatively. (b) Pictorial bond structure of the rotated Hamiltonian. The shaded region represents a period of six sites in the U6 rotation.
As demonstrated in Ref. [27], the Hamiltonian in Eq. (1) possesses two peculiar properties concerning symmetries in the parameter space. One is a self-dual relation, which implies that the eigenvalue E of H satisfies E(g) = gE(1/g) where g ≡ g y /g x is the relative bond strength. The other is a mirror symmetry with respect to the Γ axis, i.e., E(K, Γ) = E(K, −Γ). These relations cause us to consider the phase diagram mainly in the right half circle of Fig. 2(a), where Γ ≥ 0 and g ∈ [0, 1].
The interactions are parameterized as K = sin θ and Γ = cos θ with θ ∈ (−π, π], and the ground-state phase diagram is mapped out by the DMRG method [40][41][42]. In the DMRG calculation, we mainly adopt the PBC to remove the edge effect, while the OBC is also used to study edge excitations. Given the structure of the unit cell in the U 6 rotated basis, we choose the chain length L to be strictly a multiple of six. In addition, up to 2000 block states are kept so as to maintain a small truncation error of ∼ 10 −7 at most. Figure 2(a) depicts the full phase diagram within the region g ∈ [0, 2] and θ ∈ (−π, π] for the bond-alternating ○ denotes the isotropic case of g = 1. Black solid circles at θ = π/4 and −π/4 represent the hidden SU (2) FM and AFM Heisenberg chains, respectively. In the right half circle, there is a Haldane phase, a gapped Kitaev phase whose two regions are distinguished as Kitaev-I (FM Kitaev region) and Kitaev-II (AFM Kitaev region) domains, two dimerized phases located at g < 1 and g > 1, respectively, and three magnetically ordered states termed FMU 6 , MI and MO (see text for details). (b) Zoom-in region of 0.10π ≤ θ ≤ 0.18π. spin-1 K-Γ chain. In either the FM or AFM Kitaev limit, there is a gapped Kitaev phase, whose regions shall be distinguished as the Kitaev-I phase and the Kitaev-II phase, respectively, for the sake of clarity. The Kitaev phases in the two regions share the same ground-state and thermodynamic properties, although the former is more fragile against Γ interaction. Remarkably, such an asymmetric stability of the Kitaev phases in the FM and AFM Kitaev limits also persists in the spin-1/2 analogy and may be a general feature of the K-Γ model. For example, in the spin-1/2 K-Γ chain, it is found that the FM Kitaev limit is merely a multicritical point while there is a finite region near the AFM Kitaev side [22,24,27]. In the spin-1/2 honeycomb-lattice K-Γ model, the region of the FM Kitaev QSL is considerably shrunken when compared with its AFM counterpart [8]. Of note is that the asymmetry is proposed to come from the interplay of two flux-pair hopping processes, whose magnitudes depend crucially on the sign of the Kitaev interaction [46]. Going back to the spin-1 K-Γ chain, by increasing Γ interaction from the isotropic FM Kitaev limit, the Kitaev-I phase survives up to |Γ|/K = −0.10(1), followed by a Haldane phase which holds a nonlocal SOP, a finite excitation gap, and also two gapless edge modes. The Haldane phase exists in a wide anisotropic region of 0.60 g 1.70, and two partially dimerized phases are then induced via continuous QPTs. The difference between the dimerized phases lies in that there is a stronger x (y)-type bond in the inner (outer) circle of g = 1. Oppositely, the Kitaev-II phase near the AFM Kitaev limit is more robust and occupies a larger territory. When K and |Γ| are comparable, a magnetically ordered state with eight-fold ground-state degeneracy is favored. In addition, two narrow regions of extra magnetic orderings are induced at modest anisotropy [see Fig. 2

A. Characters of Haldane phase
First of all, let us consider the isotropic spin-1 K-Γ chain with g = 1, which is amenable to the revealment of crucial features of the entire phase diagram. The (K = −1, Γ = 1) point at θ = −π/4 is a hidden SU (2) AFM Heisenberg point whose ground state is the Haldane phase [31]. For this phase, it owns a (bulk) excitation gap ∆ e 0.410479 [47,48] and possesses a nonlocal SOP defined as [43] where spins are situated in the U 6 rotated basis. The gapped Haldane phase could survive against competing interactions with θ = −π/4 since it is a SPT phase which usually undergoes a QPT with a closure of its excitation gap. To estimate the region of the Haldane phase, we calculate the SOP O z H in the range of θ ∈ [−π/2, π/2], as presented in Fig. 3(a). Except for the accidental FM point at θ = π/4 where O z H is equal to the saturated value, O z H is robust and suffers from a tiny finite-size effect in a wide region of −0.4685(5) < θ/π < 0.1270 (5). Here, the transition points and the corresponding error bars are estimated from the vanishing points of SOP. These values are consistent with the independent estimate which shall be shown later. With increasing θ, O z H slowly grows deep in the Haldane phase. Specially, the value of O z H is 0.4935(2) in the pure Γ limit of θ = 0, which is larger than the value of 0.3743(1) for the AFM Heisenberg chain [47]. Near phase boundaries, O z H has a jump at θ t,1 = −0.4685(5) and varies smoothly around θ t,2 = 0.1270(5), indicative of a first-order and a continuous QPT, respectively. In the vicinity of the FM SU (2) point, there is a magnetically ordered phase with an eight-fold ground-state degeneracy. These degenerate states could be unified by the so-called η-notation [49], from which the spins within the six-site unit cell are ) are the Ising variables. The three η's are free to choose either 1 or −1 without altering the energy, giving rise to the degenerate manifold. Figure 3(b) shows the values of (a, b, c) in the same parameter region as Fig. 3(a). a and b are equal as a reminiscence of g = 1, and they compete with c when changing the relative value of K and Γ. Since c is always different from a and b, the phase exhibits an out-of-plane spin structure (cf. Eqs. (5) and (6)) [27] and thus is termed the M O phase. The M O phase is found to exist in the parameter range of θ t,2 < θ < θ t,3 where θ t,3 = 0.3845 (5).
The Haldane phase is known to carry a finite excitation gap, which is the key point of Haldane's conjecture [31]. We find that finite-size effects of the energy of the ground state and low-lying excited states are less pronounced in the PBC, and the numerical result suggests that the Haldane phase indeed possesses a nonzero excitation gap in the whole region (not shown). Taking the isotropic Γchain (θ = 0) as an example, our result indicates that the ground-state energy e g = −0.96226390(3) and the excitation gap ∆ e = 0.229135 (7). Noteworthily, the entanglement spectrum [50] is two-fold (four-fold) degenerate for the lowest-lying levels under the OBC (PBC), which is a hallmark of the SPT phase [51]. In addition, the Haldane phase also acquires two free edge spin-1/2s which account for the four-fold degenerate ground state in the thermodynamical limit under the OBC [47]. In the isotropic spin-1 Heisenberg chain with an even number of sites, the spin-1/2 edge modes are coupled with an effective AFM interaction. Accordingly, the four-fold degeneracy is split into a single state and a Kennedy triplet state for finite-size systems [35]. However, this picture partially breaks down if |K| = |Γ|. To explain this, we have calculated the first four energy levels E υ (υ = 0-3) under an open chain with L = 24. For the sake of comparison, we use the average energȳ E = (E 0 + E 1 + E 2 + E 3 )/4 as a reference scale and introduce a relative energy as ε υ = E υ −Ē. Figure 4 presents ε υ in the window of θ/π ∈ [−0.5, 0.2]. Throughout the Haldane phase, the lowest energy level (black asterisk) is always unique, while higher ones are partially degenerate. When θ < −π/4 (i.e., K < −|Γ|), the quasi-degenerate states are of "1+2+1" structure, However, they are of "1+1+2" structure when θ > −π/4 (i.e., K > −|Γ|). Higher levels recover as a Kennedy triplet state once θ = −π/4 where the model is reduced to the spin-1 SU (2) Heisenberg chain.
For the spin-1 Heisenberg chain, the existence of two edge spin-1/2s could be measured numerically by calculating the on-site magnetization S z l in theS z tot = 1 subspace, which is found to decay exponentially towards the middle of the chain [47,52]. Notably, the edge states can also be probed by a couple of experimental methods (see Ref. [53] and references therein). Unfortunately, detecting the edge modes is more intractable in our case for the lack of U (1) symmetry. In this regard, we adopt a cumulant correlation function are the accumulated magnetization on the left and right edges, respectively [54]. The C z n represents the correlation between two edge spin-1/2s regardless of the singlet or triplet sector. We calculate the cumulant correlation function C z n on an open chain with 128 sites, the size of which is far larger than the correlation length of the Haldane phase. The behaviors of C z n in the Heisenberg chain (θ = −π/4) and Γ-chain (θ = 0) are shown in Fig. 5. While both curves have strong oscillations when n is small, they saturate to −1/4 and −1/6 or so, respectively, as n 20. We note that, for the lack of Kitaev interaction, the correlator C z n in the Γ-chain is thus only two-thirds of that in the Heisenberg chain. Given the discrepancy between the saturated values of the correlator, we confirm that there are two spin-1/2 edge modes which interact antiferromagnetically with each other as revealed by the minus sign of C z n .

B. Unusual excitations of Kitaev phase
We now turn to studying the nature of the Kitaev phase in the vicinity of two Kitaev limits. In the spin-1 Kitaev spin chain, the bond-directional Ising interactions allow for the existence of many Z 2 quantities that commute with the Hamiltonian. To this end, it is natural to introduce an on-site operator Σ α l = e iπS α l [36], which mutually commutes with each other and happens to be I−2(S α l ) 2 for the spin-1 case. The bond-parity operatorŝ W l on the odd/even bonds [37], commute with the Hamiltonian of the Kitaev spin chain with eigenvalues being ±1 (which explains the Z 2 nature). In the Kitaev limit, the ground state is unique under the PBC and lies in the sector with allŴ l = +1 [36]. It is gapped with an extremely short correlation length ξ 1 since the spin-spin correlation beyond the nearest-neighbor bonds is extremely small. The first excited state corresponds to flipping the eigenvalue of any of the bond parity operators to −1 while leaving the rest unchanged, generating an L-fold degenerate first excited state [37]. We find that the ground-state energy of the isotropic Kitaev chain is -0.603560592(3) 1 , which matches perfectly with a previous estimate by exact diagonalization method [36]. Also of interest is that the energy is very close to that of the two-dimensional spin-1 Kitaev honeycomb model [55]. To capture the QPT driven by Γ interaction, we define the averaged bond density which should still be +1 for the ground state in the Kitaev limit. Figure 6(a) shows W b as a function of θ for different length L in both FM and AFM Kitaev limits.
In the left panel, W b experiences a dramatic jump at θ t,1 /π = −0.4685(5) (i.e., |Γ|/K −0.10), favoring the first-order QPT as revealed by the SOP of the neighboring Haldane phase (cf. Fig. 3(a)). However, W b varies smoothly as θ decreases from π/2. To locate the transition point, we take the first-order derivative of W b with respect to θ (see inset). It can be observed that there is a peak whose position is θ t,3 /π = 0.3845(5). In Fig. 6(b) we show the excitation gap ∆ e = E 1 − E g , defined as the energy difference between the ground state (E g ) and the first excited state (E 1 ), around the transition points. In the range of −0.50π ≤ θ ≤ −0.40π, both Kitaev-I phase and Haldane phase are gapped, and transition at θ t,1 is of first order due to the level crossing. In addition, excitation gap of the Kitaev-II phase at θ = 0.50π is equal to that of the Kitaev-I phase at θ = −0.50π, with a value of ∆ e = 0.17634970 (2). Within the Kitaev-II phase, ∆ e decreases gradually and is vanishingly small when approaching the phase boundary of the Kitaev-II to M O transition.
So far, we have shown that the Kitaev-Haldane transition at θ t,1 /π = −0.4685(5) is a topological first-order transition which can be recognized by the jumps in nonlocal SOP and Z 2 quantities. Both of the two are gapped with unique ground states. The Haldane phase, in particular, is a preeminent example of SPT phase with shortrange entanglement. By contrast, the Kitaev phase is not a SPT phase since its lowest entanglement spectrum is unique rather than doubly degenerate. The entanglement spectrum is a quantity which encodes the spectral information of the subsystem A with a reduced density matrix ρ A . It is defined as − ln λ υ where λ υ is the eigenvalue of ρ A with υ λ υ = 1 [50]. The entanglement spectrum of the Haldane phase is universally known to be doubly degenerate; it is two-fold degenerate under the OBC, while it is four-fold degenerate under the PBC [51]. As shown in the inset of Fig. 6(b), the Haldane phase at θ = −0.40π indeed has a four-fold degeneracy entanglement spectrum under PBC. However, spectra of the Kitaev phase at θ = −0.50π and −0.48π are unique, in contradiction to the character of a SPT phase.
In the case of the OBC, the Haldane phase has a fourfold degenerate ground state owing to the two spin-1/2 edge states. Interestingly, the ground state of the Kitaev phase under the OBC is also four-fold degenerate, relating to the breakdown of two marginal Z 2 quantitieŝ W 1 andŴ L−1 , see the inset of Fig. 7(a) nontrivial observation is that summations of both Ŵ 1 and Ŵ L−1 within the degenerate ground state are extremely close to zero, which resemble those of zero edge modes. Above the ground state, the lowest excitations come from flipping further Z 2 quantitiesŴ 3 andŴ L−3 near the boundaries. To detect the trails of excitations, we define the local excitation energy ω l as where · g and · e represent the expectation values with respect to the first and fifth energy levels, respectively. Figure 7(a) shows the distribution of excitation energy ω l for the isotropic Kitaev spin chain (g = 1) under OBC. It can be found that the excitation energy is indeed localized at the very boundary parts, with ω l being zero in the central region. The total excitation gap ∆ e = l ω l is 0.11594909(2), which is roughly 2/3 of the bulk gap 0.17634970(2) reserved in the PBC case. We propose that the loss of excitation gap is attributed to Ising-type couplings and short correlation length ξ 1 of the Kitaev spin chain. Because of the OBC imposed, correlations within the marginal x-type bonds of (1, 2) and (L − 1, L) are enhanced, leading to effective x-type couplings among these four sites. Hence, the excitations are entangled and the marginal Z 2 quantities are no longer conserved. Because of the extremely short correlation length, the excitations are gathered at edges of the chain. Unfortunately, the emergent S x 1 S x L correlation does not contribute to the energy. Instead, it trammels the excitations and thus reduce the excitation gap. Our analysis implies that, if excitations are not bounded at the edges, then excitation gap ∆ e would not depend on the boundary condition severely. To check the conjecture, we consider a stronger y-type bond, say g = 1.2. We show the distribution of excitation energy ω l in Fig. 7(b), and it is clearly found that the excitation comes from the very middle of the chain. More importantly, discrepancy of the excitation gap between the OBC and PBC is insignificant. We have shown the excitation gap ∆ e as a function of anisotropy g under both OBC (red circle) and PBC (blue square) in Fig. 7(c). When g 1.06, the excitations are localized at edges and the excitation gap under OBC is somewhat smaller.

A. Haldane-Dimer transition
In this section we study the effect of imbalanced bond strength with g = 1. We recall that our model (1) is reduced to the bond-alternating spin-1 Heisenberg chain when θ = −π/4, which undergoes a Haldane-dimer transition at g = 0.587(2) that belongs to the Gaussian universality class with a central charge c = 1 [56][57][58][59]. The dimerized phase is characterized by a stable alternation of nearest-neighbor spin-spin correlations, which is equivalent to the difference of S z iS z j between the odd bonds and even bonds within the six-site unit cell, namely Due to the inherent bond alternation, the transitional symmetry of the neighboring sites is broken innately. When g < 1, the x-type bond is stronger than the y-type, and vice versa. Hence, the ground state of the dimerized phase is unique with a finite energy gap. Figure 8(a) shows the energy gap ∆ e for three different length L = 48, 96, and 144. As g goes from 0 to 2, ∆ e exhibits two remarkable valleys where the value becomes smaller and smaller as L is increased. We make a linear extrapolation of the series of minimal values around each valley marked by a pink diamond (g < 1) or a cyan pentagram (g > 1), and it is shown in the inset that the excitation gap ∆ e indeed closes at the transition point in the thermodynamic limit. Noteworthy, we find that the transition points are g t,1 = 0.613(2) and g t,2 = 1.633(2), satisfying the self-dual relation as g t,1 · g t, 2 1. In   (2) and is vanishingly small otherwise. Although O z D is nonzero as long as g = 1, it changes abruptly on the brink of the Haldane phase. We have taken the first-order derivative of O z D (not shown), and the peak position coincides with the critical point g c = 0.613 (2) where the SOP of the Haldane phase nicely vanishes.

B. Magnetically ordered phases
As demonstrated, the dimerized phase is favored by sufficient bond alternation. In the case of strong anisotropy where g 0.4 or g 2.5, there is a direct transition between the dimerized phase and the M O phase. Otherwise an intermediate region is stabilized as a consequence of the interplay between competing interactions and modest bond alternation. A more careful inspection shows that it harbors two magnetically ordered states (see Fig. 2(b)) with distinct spin patterns. On the side near the dimerized phase, it is a FM U6 phase whose spins align along one specific direction of ±x, ±ŷ, or ±ẑ.
For example, it could be S 1 , S 2 , S 3 ; S 4 , S 5 , S 6 = c, b, b; c, a, a ẑ (11) where a, b and c are intensities of the magnetization along theẑ direction and a, b, c ≤ 1. The other is an in-plane magnetic phase (termed M I ) where one of the spin components is zero. We find that spins within the magnetic unit cell reads and For this phase, the value of a, b, c should be bounded by √ a 2 + b 2 ≤ 1 and c ≤ 1/ √ 2. The series of QPTs between the magnetically ordered states could be signified by the magnetization of (a, b, c). In this regard, we focus on the line of g = 2.0 (which is equivalent to g = 0.5 because of the self-dual relation) and calculate these values from the spin-spin correlation functions on a chain of length L = 72. Figure 9 shows the values of (a, b, c) as a function of θ in the range 0.10π ≤ θ ≤ 0.18π. Due to the finite-size effect, a, b, and c are nonzero but very tiny when θ/π < 0.1345 (5). We have checked that they will eventually go to zero in the thermodynamic limit. As θ exceeds 0.1345(5), the system enters into the FM U6 phase whose values of (a, b, c) are accidentally bigger than 1/2. Afterwards, the system undergoes another two first-order QPTs at which the ground state evolves from the M I phase to M O phase, respectively.

V. DOUBLE-PEAK SPECIFIC HEAT IN THE KITAEV PHASE
In this section we go beyond the ground-state study by calculating thermodynamic quantities in the spin-1 Kitaev chain to elucidate the unusual excitations of the Kitaev phase. We recall that specific heat C v of the spin-1/2 Kitaev honeycomb model is well-recognized to exhibit a double-peak structure at two different energy scales, signifying two kinds of Majorana fermions resulting from fractionalization of spin degrees of freedom [2]. The high-temperature peak relates to the enhancement of short-range spin-spin correlations while the lowtemperature peak comes from freezing of fluxes [3]. In addition, the thermal entropy displays an approximately half plateau with the value 1 2 ln 2k B per site in the intermediate crossover region, in accordance with a half release of entropy around each peak. Moreover, existence of double-peak specific heat is also reported in the spin-1 Kitaev honeycomb model, yet precise position of the low-temperature peak is blurry because of the strong finite-size effect [17]. We note in passing that there has been a renascent interest in the thermodynamics of the Kitaev QSL quite recently [60][61][62][63].
In contrast to the spin-1/2 Kitaev chain which only possesses a sole peak in the specific heat (see Appendix A for detail), our study unambiguously suggests that the spin-1 Kitaev chain displays a double-peak specific heat. To illustrate this, we begin by introducing the partition function Ξ = Tre −βH with β = 1/k B T (herein, the Boltzmann constant k B = 1), which is generally the starting point to calculate thermodynamic quantities. Consequently, the free energy is F = −β −1 ln Ξ and the internal energy is U = − ∂ ln Ξ ∂β . The specific heat is thus calculated by and the thermal entropy is given by with S 0 being the residual entropy at zero temperature. The direct way to calculate the partition function Ξ is by diagonalizing the Hamiltonian from which the entire energy spectrum {E υ } is readily available. Needless to say, this route is strongly limited by the system size and we consider a six-site closed chain for simplicity. Surprisingly, physical quantities such as C v and S suffer from a weak finite-size effect and the results are fairly close to those in the thermodynamic limit. The first few energy levels are shown in Tab. I, and their degeneracies ρ υ are {1, 6, 6, 2, · · · }.
A more reliable way to obtain Ξ is by virtue of finite-temperature computational techniques such as the TMRG method [44,45], which is an extension of the DMRG method to finite temperature (T = 0). The · · · · · · · · · · · · TMRG method relies on the quantum-classical correspondence by way of the Trotter-Suzuki decomposition, and represents the partition function as a trace of a series of transfer matrices. It deals directly with an infinite spin chain and the errors come from the Trotter-Suzuki step τ = β/M (M is the Trotter-Suzuki number) and truncated number of states m [44,45]. We note that an additional reorthogonalization procedure is applied so that more block states can be kept [64]. The TMRG method has been successfully applied to various quantum spin chains where several thermodynamic quantities such as specific heat and magnetic susceptibility could be calculated with high precision [65][66][67][68]. For the sake of accuracy, we fix τ = 0.01 and m = 1024, which is enough to give a satisfactory precision in our simulation down to the lowest temperature of 0.0033. Figure. 10(a) shows the specific heat C v as a function of temperature T up to 10. Here, exact result on an extremely small system size of L = 6 (black dotted line) and the TMRG calculation on an infinite size system (red solid line) are shown. It can be found that C v is almost system-size independent as the difference between the two limit cases is quite small. The specific heat C v acquires two peaks at a low temperature T l 0.0582 and a high temperature T h 0.5860, respectively. The low-temperature peak is more pronounced and we propose that it relates to the large degeneracy of the low-lying excited states. To demonstrate this, we firstly consider a two-level system which consists of the ground state and the first excited state out of the sixsite closed chain (see Tab. I). The degeneracies of the two states are ρ 0 = 1 and ρ 1 = 6, respectively, with an energy gap ∆ κ ≈ 0.18018574 between them. For this system, the partition function Ξ = ρ 0 + ρ 1 e −β∆κ , and the specific heat is which depends on the energy gap ∆ κ and relative degeneracy ρ 1 /ρ 0 = 6. The specific heat in Eq. (16) shows a peak at the extreme temperature T p = ∆ κ /x p , where x p = 3.23565205 · · · is determined by the transcenden- tal equation x−2 x+2 e x = ρ 1 /ρ 0 = 6. This yields the extreme temperature T p ≈ 0.0557, which is fairly close to T l 0.0582. The drawback of this oversimplified approximation is that intensity of the specific heat is smaller than the actual value. However, the intensity could be significantly improved by considering a four-level system with a partition function The corresponding specific heat is also shown in Fig. 10(a) (blue dash-dotted line), and position and extremum are both close to the TMRG result on the infinite system. With increasing system size, the dimension of the Hilbert space will enlarge exponentially, and a further multi-level system should be considered to reproduce the low-T peak (see Appendix B for a discussion on a 12-site closed chain). Figure. 10(b) shows the behavior of entropy S defined in Eq. (15). The entropy decreases from its saturated value of ln 3 rapidly with the lowering of the temperature in the neighborhood of T h and T l . The 1/2-plateau of the entropy is smeared out in the intermediate region between the two temperature scales. Instead, a shoulder in entropy is observed and the entropy becomes ∼ 1 2 ln 3 at a temperature of T m 0.20. Interestingly, it is exactly the same temperature at which specific heat shows its local minimum. To understand physical mechanisms of the double-peak structure in the specific heat, we calculate the averaged bond density W b (see Eq. (8)) and the averaged nearest-neighbor correlator The results of W b and S b are shown in Fig. 10(c). Below the low-T peak at T l , W b is almost unchanged with a value of ∼ 1. It then decreases to 1/9 successively with the increase in temperature in the intermediate region 2 . Therefore, the low-T peak originates from the local conserved quantity. On the other hand, the high-T peak is closely related to the growth of short-range spin-spin correlations. Above the high-T peak at T h , S b is very small, signifying a paramagnetic phase. It is then enhanced dramatically by decreasing the temperature and 2 In the high-T paramagnetic phase, we have (S x l ) 2 = (S y l ) 2 = (S z l ) 2 = 2/3 for ∀ l in the spin-1 system. The bond-parity oper-atorŴ l is also uniformly distributed. Taking the x-type operator for instance,Ŵ l = (Σ x l ) 2 = [1 − 2(S x l ) 2 ] 2 = 1/9.
is stabilized around 0.60, which is merely the absolute value of the ground-state energy of the isotropic spin-1 Kitaev chain. We would emphasize that the double-peak specific heat is a universal behavior of the Kitaev phase in spite of anisotropy g. For this purpose, we have calculated the specific heat and entropy for several different g up to 2.0, see Fig. 11. It can be found that the high-T peak is very pronounced while the low-T peak is shifted to a lower temperature as g is increased. As can be seen from Fig. 11(b), there is a 1 2 ln 2-plateau of entropy in the middle region, and the plateau will last for a larger temperature scale with increasing anisotropy. Thus, in the large anisotropy limit where g is infinity (or equivalently, g → 0), the low-T peak vanishes and there is a residual entropy at the zero temperature due to the 2 L/2 -fold ground-state degeneracy. In this circumstance, a spin-1 Ising bond of either x-type or y-type is capable of revealing the residual entropy. The energy spectrum of this two-site model is {−1, 0, 1}, with a degeneracy of {2, 5, 2}, respectively. Thus, the partition function is known to be Ξ = 4 cosh(βK) + 5. According to Eq. (15), we arrive at the entropy S = 1 2 ln 4 cosh(βK) + 5 − 2(βK) cosh(βK) 4 cosh(βK) + 5 .
For large enough temperature, the entropy in Eq. (19) is expected to yield ln 3. As T becomes infinitely small, we find that S(T → 0) = 1 2 ln 2. In this circumstance, the low-temperature peak in the specific heat disappears, and the sole peak locates at T 0.3896 with the maximal valueC v 0.4974.

VI. CONCLUSION
In this paper, we studied a bond-alternating spin-1 K-Γ chain, focusing on the nonmagnetic Haldane phase and Kitaev phase that exhibit unusual excitations. The Haldane phase is an outstanding example of a SPT phase which is gapped with short-range entanglement, while the Kitaev phase is a one-dimensional incarnation of the Kitaev QSL and is not a SPT phase since its lowest entanglement spectrum is unique. Whereas both of the phases have unique ground states under the PBC, the degeneracies of the first-excited states are different. It is triplet degenerate for the former, while it is L-fold degenerate with L being the chain length for the latter as a result of the bond-resolved Z 2 quantities. Interestingly, they both possess four-fold degenerate ground states in the case of the OBC. The degeneracy in the Haldane phase comes from the spin-1/2 edge states, while it originates from two marginal Z 2 quantities in the Kitaev phase. On top of the degenerate ground state in the Kitaev phase, the spatial profile of the excitations highly relies on the relative bond strength g ≡ g y /g x . The excitations are confined at the boundaries of the chain when g 1.06 and locate at the very middle otherwise. In the former case, the excitations at edges are entangled, weakening the excitation gap as compared with its PBC counterpart. The quantum phase diagram also contains two dimerized phases which undergo continuous QPTs to the Haldane phase when tuning the anisotropy g. In addition, three magnetically ordered states are identified, of which the M O phase is the most energetically favored and is situated alongside the AFM Kitaev phase.
We also investigated the thermodynamic behaviors of the Kitaev phase, which is found to exhibit a fascinating double-peak structure in the specific heat. During the low-temperature and high-temperature crossover region, the entropy is released gradually without generating a plateau. Pertaining to the origin of the double peaks, we propose that the high-temperature peak relates to the enhancement of nearest-neighbor spin correlation while the low-temperature peak comes from freezing of Z 2 quantities and is relevant to the highly degenerate low-lying excited states. We also find that the double-peak specific heat is robust against anisotropy g, although the lowtemperature peak is shifted to lower temperature steadily as g deviates from 1. Notably, a 1 2 ln 2-plateau of entropy appears in the crossover region.
In closing, we would like to make some remarks on the connection between the Kitaev phase and the soughtafter QSL in the spin-1 Kitaev honeycomb model [17][18][19][20]. The spin-1 Kitaev QSL has been gaining much attention over the years because of the growing interest in high-spin Kitaev materials [16]. For both phases, there is no spontaneously symmetry breaking in the ground states, and they are gapped with extremely short-range spin-spin correlations. Noteworthily, they both exhibit anomalous double peaks in their specific heat. However, the intrinsic nature of the one-dimensional gapped Kitaev phase remains to be explored in future works. For example, it would be intriguing to know to which class the Kitaev phase belongs from the perspective of symmetry fractionalization [69]. Also, calculation of the low-energy excitation spectrum from the dynamic structure factor is capable of probing elementary excitations [70]. On all counts, by using a high-precision numerical study of the Kitaev phase, our work would offer some insights into the spin-1 Kitaev   The Hamiltonian of the Kitaev spin chain reads [71,72] where S l = (S x l , S y l , S z l ) is the spin-1/2 operator at site l, and L is the total number of sites. By using a spin duality transformation, it could be rewritten as [10] H K = L l=1 g xS x 2lS x 2l+2 + g y 2S y 2l , which is a diluted transverse field Ising model. For this model, the specific heat is exactly known as [73] C v = 1 2 π 0 dk π β k 2 2 sech 2 β k 2 (A3) where β = 1/(k B T ), and the dispersion energy is k = g x 1 + g 2 + 2g cos k/2 with g = g y /g x . We note that the prefactor 1/2 in Eq. (A3) comes from the fact that only one half of the spins (i.e., spins at even sites) are involved in the Hamiltonian of Eq. (A2). Figure 12(a) shows the specific heat C v of the isotropic (g = 1) spin-1/2 Kitaev chain at finite-size system of L = 6 (red), L = 8 (green), L = 10 (blue), and L = 12 (pink). It can be found that there is a pronounced peak at T 0.3162, below which there is a subleading peak at a lower temperature which becomes smaller and smaller as L is increased. We emphasize that the low-temperature peak is a finite-size effect and it will disappear as L goes to infinity (see the thick black line). As shown in Fig. 12(b), a linear extrapolation of the subleading peaks indeed gives a zero value when L → ∞. In the low tem-perature region where T 0.1, the specific heat is subject to the asymptotic behavior C v (T ) πT /6 [74], signifying a gapless system. This, in turn, demonstrates the dramatic difference between the spin-1/2 and spin-1 Kitaev chains as the latter is gapped and presents a stable double-peak structure in the specific heat.
Appendix B: Low-temperature peak in a 12-site spin-1 Kitaev chain In Fig. 10(a), we have demonstrated that the low-T peak of the specific heat in the spin-1 Kitaev chain relates to the large degeneracy of the low-lying excited states. To further strengthen such a conclusion, we now show the evolution of the low-T peak by increasing the number of energy levels in a 12-site closed chain. For this system, the ground state is unique while the firstexcited state is 12-fold degenerate, separated by an excitation gap ∆ κ ≈ 0.17630343. Furthermore, the degeneracies of the first nine energy levels are successively {1, 12, 12, 12, 6, 12, 12, 24, 12}. Figure 13 shows the low-T peak based on a two-level ansatz (with the first 13 energy levels, red dot-dashed line) and a nine-level ansatz (with the first 103 energy levels, blue dotted line). The TMRG result (thick black line) is also shown for comparison. We find that the two-level system could roughly recover the low-T peak, although the position and height of the peak deviate from the TMRG result. However, the nine-level system with 103 energy levels could significantly improve the result. We note that the number of energy levels is only a rather small portion (∼ 2 × 10 −4 ) of the whole energy spectrum whose dimension is 3 12 = 531441. With the increase in the system size, we believe that an even smaller portion of the whole energy levels will nicely produce the low-T peak. In this sense, we highlight the importance of the degenerate low-lying excited states in generating the low-temperature peak of the specific heat.