Study of the phase diagram of the Kitaev-Hubbard chain

We present a detailed study of the phase diagram of the Kitaev-Hubbard chain, that is the Kitaev chain in the presence of a nearest-neighbor density-density interaction, using both analytical techniques as well as DMRG. In the case of a moderate attractive interaction, the model has the same phases as the noninteracting chain, a trivial and a topological phase. For repulsive interactions, the phase diagram is more interesting. Apart from the previously observed topological, incommensurate, and charge density wave phases, we identify the “excited state charge density wave” phase. In this phase, the ground state resembles an excited state of an ordinary charge density phase, but is lower in energy due to the frustrated nature of the model. We find that the dynamical critical exponent takes the value z 1.8. Interestingly, this phase only appears for even system sizes, and is sensitive to the chemical potential on the edges of the chain. For the topological phase, we present an argument that excludes the presence of a strong zero mode for a large part of the topological phase. For the remaining region, we study the time dependence of the edge magnetization (using the bosonic incarnation of the model). These results further expand the region where a strong zero mode does not occur.


I. INTRODUCTION
Topological phases of noninteracting fermions have been studied extensively and the full classification of the possible phases has been found [1][2][3][4]. The classification scheme is based on time-reversal symmetry, particle-hole symmetry, and chiral symmetry and which phases are possible depends on the spatial dimension. To a given class and spatial dimension one can attribute a group, i.e., Z 2 or Z, that describes the possible gapped topological phases in that class. As an example the class BDI has all three symmetries (squaring to one) and in one dimension its gapped phases are labeled by the group Z. The prototypical model belonging to this class is the Kitaev chain (with real couplings) [5]. In the topological phase of the Kitaev chain, the model hosts Majorana zero modes on the edges of the chain. This zero mode results in a fully doubly degenerate many-body spectrum.
For interacting fermions there is no universal classification scheme. Nonetheless for the class BDI it was shown that in the presence of interactions the different topological classes correspond to the elements of the group Z 8 [6,7]. One should note that this classification only concerns gapped phases and the phase diagram of a general model has gapless points or regions as well. Therefore it is interesting not only to wonder about how the interaction affects a topological phase [8] but also how gapless phases emerge due to interactions.
The main aim of this paper is studying the phase diagram of the Kitaev chain in the presence of a density-density interaction, see Refs. [9][10][11][12][13][14] for earlier results for this model. We refer to this model as the Kitaev-Hubbard chain. We study both the attractive and repulsive regime. In addition to the trivial and topological phases which are inherited from the Kitaev chain, the model has an incommensurate phase, a charge density wave phase, and a new phase which is sensitive to the system size and the chemical potential at the boundaries of the chain. Moreover, we study the possibility of having a full doubly degenerate many-body spectrum in the topological phase of this interacting model.
The paper is organized as follows. In Sec. II we briefly review both the classical and quantum axial next-nearestneighbor Ising (ANNNI) model and present the model of interest for this paper, namely the Kitaev-Hubbard chain. We show that the quantum ANNNI model and the Kitaev-Hubbard chain are dual to each other and present the main result, that is the phase diagram of this model. We continue by presenting our analytical results based on bosonization and the numerical results based on density matrix renormalization group (DMRG) [15,16] for an attractive interaction in Sec. III. In Sec. IV we move on to the repulsive regime, for which there is an incommensurate phase and a new phase which was missed previously. Moreover, we will discuss the properties of the highest excited states as well as the finite temperature features of the model in Sec. V. We conclude the paper in Sec. VI.

II. THE MODEL AND THE PHASE DIAGRAM
To introduce the model we set the scene by reminding the reader of the classical Ising model, which is the cornerstone of our understanding of phase transitions [17][18][19][20].
To define this model consider a square lattice where on each site, with coordinate (i, j), one has a classical spin, s i, j = ±1. For a configuration of spins {s} the energy is We assume thatJ 1 ,J ⊥ > 0. The classical Ising model in two dimensions was initially solved by Onsager [17], later other solutions were found, see [18]. The model has two phases. For temperatures T less than the critical temperature T c , it is in the ferromagnetic phase and for T > T c it is in the disordered phase. As a generalization one can consider the ANNNI model where one adds an extra interaction along one axis, which could favor ferromagnetic or antiferromagnetic alignment [9,21,22]. The energy for a given configuration of spins of the ANNNI model is We follow the usual notation and defineκ = −J 2 J 1 . Forκ < 0 there is no frustration, and the axial term only modifies the critical temperature. Forκ > 0, however, the situation is more subtle. For small and positiveκ the model has two phases, a ferromagnetic and paramagnetic phase. The pointκ = 1/2 is a multicritical point, where four phases meet, namely the ferromagnet, paramagnet, a floating phase, and the so-called "antiphase" [9]. In the antiphase, the axial interaction dominates and there are four spin configurations that minimize the energy. For the floating phase however, all interactions are important. Whether the floating phase continues to infiniteκ is an open question, which we study in this paper by analyzing the corresponding quantum model.
The quantum counterpart of the classical Ising model is the transverse field Ising model (TFIM) [23][24][25]. This model is defined on a chain. On each site there is a spin-1/2 degree of freedom and the Hamiltonian for a chain of size L can be written as where τ α are Pauli matrices. One can relate the parameters J 1 and B in terms of coupling constants of the classical model and the temperature [25]. The TFIM shows a quantum phase transition at T = 0. For |B| < |J 1 | the model is in the ferromagnetic phase while for |B| > |J 1 | it is in the disordered phase. The classical Ising model in two dimensions and the one-dimensional TFIM are in the same universality class. The quantum ANNNI model was also studied [10][11][12]14]. The Hamiltonian reads [9] In analogy with the classical ANNNI model, we define κ = − J 2 J 1 . For κ < 1 2 the model has only ferromagnetic and paramagnetic phases. In this regime there is a second order phase transition between the ordered and the disordered phase at B c (κ ). The point κ = 1 2 is a multicritical point from which the floating phase emerges. For very large κ the dominant term is the next-nearest-neighbor term in Eq. (4) which gives rise to four ground states. Previous studies found that the floating phase survives till κ 5 [11,12]. We will present our results for infinite κ by studying the dual model.
To study large frustration, κ 1, we employ the duality map as follows: We further perform an on-site rotation, Using these two transformations one can transform the Hamiltonian for the quantum ANNNI model, Eq. (4), to an interacting quantum Ising model, Note that the sign of B and J 1 are not important and we will assume that both are positive. The sign of J 2 , however, is crucial. Hence we use B as our energy scale and study the following Hamiltonian: By performing a Jordan-Wigner (JW) transformation [26], one arrives at the Kitaev-Hubbard chain, For U = 0, the model reduces to the TFIM (or the Kitaev chain [5]) which is exactly solvable [24,27]. In the bosonic representation the two phases are the ordered (h < 1) and the disordered (h > 1) phases which are converted to the topological and the trivial phases, respectively, in the fermionic incarnation. The phase transition between these phases is described by the c = 1 2 Ising CFT. In the topological phase the model hosts Majorana zero modes which live on the edges and has zero energy up to exponentially small corrections in the system size. An interesting question is whether the topological phase is stable in the presence of interactions [8,28]. The U term in Eq. (8) gives rise to a density-density interaction, but renders the model nonintegrable.
Another special case is h = 0 for which the model is called the XY model [27]. This model is gapped, except when U = ±1, which are critical points with central charge c = 1. The gapped phases are ordered phases in the bosonic representation and topological in the fermionic incarnation. It is worth to mention that the model in Eq. (8) is actually relevant for an array of superconducting islands [29]. Moreover, for a specific set of couplings known as the Peschel-Emery (PE) line, the ground state is exactly doubly degenerate [30]. This line lies in the topological phase and host weak zero modes [31,32], see Sec. V for more detail on this.
We obtained the phase diagram of the model in Eq. (8) as presented in Fig. 1. Parts of this phase diagram were schematically drawn previously based on the known results for the quantum ANNNI model [10,11,13,28,29]. We believe, however, that in the previous studies one phase (the green region in the plot) has been missed. This phase will be described in Sec. IV C.

III. THE ATTRACTIVE INTERACTION CASE
For attractive interactions, U < 0, there is no frustration, κ < 0, and one expects a direct transition from the topological phase (ordered phase) to the trivial phase (disordered phase), as is the case for U = 0. To study the U < 0 case, one can start from the Ising critical point, namely h = 1 and U = 0, and perform perturbative calculations to find the location of the phase transition [33]. We will, however, start from another critical point, namely the XY critical point with h = 0 and U = −1. Close to this point we can use bosonization. To do so it is convenient to transform the Hamiltonian in Eq. (8). We first do an on-site rotation, which changes the Hamiltonian to in which we defined δU = U + 1 and dropped the lower and upper bounds of the sums, since we will be working in the continuum limit. We assume that δU, h 1. We can use the JW transformation, write the Hamiltonian in terms of spinless fermions ψ j , and perform a Fourier transformation to work in momentum space. Doing so, the first two terms, namely the XY model, give rise to a gapless cos(k) band, which can be linearized around k = ± π 2 , to get the continuum model, In the last equation a is the lattice spacing and x = ja will be the continuous spatial coordinate. We first consider the δU term in Eq. (12). Plugging Eq. (13) into this term we get a term proportional to To further simplify the result, we use the dual fields φ(x) and θ (x), which obey the commutator [φ(x), θ (y)] = i (y − x), where (x) is the Heaviside step function. Therefore, ∂ x θ (x) is the conjugate momentum of the field φ(x). We employ the bosonization dictionary [34], in which α is a cutoff in momentum space. Hence we can write the δU term in terms of the dual fields, which is proportional to The h term in Eq. (12) can be treated in the same way. Using the bosonization dictionary and by dropping the fast oscillatory terms which are proportional to (−1) j , one can show that Therefore the h term reduces to The full Bosonized Hamiltonian then reads in which H 0 is the free bosonic Hamiltonian and c 1 and c 2 are some constants which depend on a and α. Since H 0 is quadratic, we can calculate the renormalization group flow of the couplings, δ and h, up to first order We should mention that we assumed that we can drop the renormalization of K, one of the Luttinger parameters. These two equations tell us that close to the XY point, h = 0 and U = −1, the curve separating the topological and the trivial phases is given by To check this analytical result we performed a numerical study and compared to the results shown in Fig the phase transition we ran DMRG using the ALPS libraries [35][36][37]. This DMRG study also serves to get acquainted with the performance of the algorithm for the model at hand.
We should note that the second term in Eq. (10) [the magnetic field in the bosonic incarnation in Eq. (8)] has been implemented as follows: Therefore in the bosonic incarnation the strength of the magnetic field in the bulk is h while it is half of this strength on the first and the last sites, h 1 = h L = h 2 . We used finite size scaling (FSS) to pinpoint the transition. Within this framework, one finds the ground state and the first excited state and use the following scaling ansatz for the energy difference between them, δ, close to the critical point [38], where L is the system size, z is the dynamical critical exponent, ν is the correlation length exponent, and F is a scaling function. For the TFIM it is known that z = ν = 1 [39]. We present our data for U = −0.4 as an example. For a given U at h = h c (U ), the quantity L z δ does not depend on the system size. In Fig. 3 it is clear that Lδ for different system sizes cross at h c = 0.406. One can use this critical value and take ν = 1 to check the validity of the ansatz Eq. (24), see Fig. 3. Indeed, the data for different system sizes follow the same functional form if one scales the h axis appropriately.
We used FSS for −0.8 U 0.25 and calculated the critical field h c (U ). The data are presented in Fig. 2. We fitted a power law with the form to the data. Using the data in the range of −0. good agreement with the perturbative approach of Ref. [33] around the point h = 1, U = 0.
To determine the central charge c along the curve we used the scaling of the excited states' energy as a function of system size L [40] and the Calabrese-Cardy (CC) formula for the entanglement entropy (EE) of a subsystem of size l, in the case of an open chain at the critical point S(l, L) [41][42][43], Here c is the central charge and S 0 is a constant. Both methods result in the central charge c = 1 2 . Therefore we have the central charge c = 1 2 all along the transition line except for h = 0 and U = −1, where we have the XY critical point with central charge c = 1 [40].

IV. THE REPULSIVE INTERACTION CASE
In this section we study the Hamiltonian Eq. (8) with repulsive interaction U > 0, which is more involved. We can qualitatively understand the physics in the weak and strong interaction limits. For weak interaction, U 1, we still expect to have two phases, i.e., a topological (ordered) phase and a trivial (disordered) phase, as in the case of the TFIM. On the other hand for very strong interaction, the model is in the Néel phase which corresponds to the charge density wave (CDW) in the fermionic representation. Apart from these two easy limits it is hard to say something about the phase(s) for intermediate interaction strength. For the earlier studies of the model with repulsive interaction using DMRG and bosonization we refer to Refs. [10][11][12]14].
To find the phase diagram of the model we performed DMRG using the ALPS libraries [35][36][37]. We used the fermionic incarnation of the model, i.e., Eq. (10), to perform the DMRG. This is useful since the total parity P is a conserved quantity, Above the parity is given in both bosonic and fermionic incarnations. Since P 2 = 1, its eigenvalues are P = ±1. Because parity is a good quantum number, the eigenstates with P = +1 are called the even sector and the ones with P = −1 the odd sector. One can restrict DMRG to be done in one sector [37].
In a default setup we performed DMRG for a chain of size L = 240, used the bond dimension χ = 500 for the Schmidt decomposition, and swept eight times. By the convention in ALPS each sweep means an optimization from the first site to the last one and again from the last site to the first one. We have, however, checked our results to be sure that the DMRG is converged. Specifically, we increased the number of sweeps up to 12, and considered different system sizes. We already note that in the green region of the phase diagram ( Fig. 1), the system size modulo four plays an important role. We discuss this issue extensively below.
To begin our tour of the repulsive regime, we present our results for the lowest two eigenstates in each parity sector, the ground state (n = 0) and the first excited state (n = 1). We label the states with P and n and their energy will be denoted by E P n . The result for h = 0.7 along the dashed line in Fig. 1 is presented in Fig. 4. One can see features of different phases in this plot. For U 1.05 the model is in a gapped topological phase in which the ground state is doubly degenerate, with opposite parities. For 1.05 U 1.29 we have a gapless incommensurate phase in which observables oscillate and the parity of the ground state depends highly on details such as the system size and the couplings. Increasing U further, for 1.29 U 1.41 we have a new phase which we believed has been missed previously. The intriguing feature of this region is that the nature of the ground state differs for even and odd system sizes. Finally for U 1.41 the system has a CDW ground state. We will later discuss the degeneracy and the parity of the ground state in this phase for different system sizes.

A. The topological phase
In Fig. 4 we start in the topological phase (the white region in Fig. 1), where the model is gapped and the ground state is doubly degenerate. The two ground states have different parity. In Fig. 5 we present the EE S(l ) as a function of subsystem size l and the (expectation value of the fermion) occupation number at each site N ( j) for the ground state in each parity sector. The difference of the energy of the two states is of the order of 10 −7 in units of hopping strength, note that we set t = 1 in Eq. (10). As is evident, S(l ) saturates to a high value in the bulk for both ground states, which is a signature of the topological phase. The occupation number is constant in the bulk for both parity sectors.

B. The incommensurate phase
By increasing the interaction strength we enter the incommensurate (IC) phase, the orange region in  phase. In this phase the ground state is unique but its parity does strongly depend on the parameters and the system size.
In Fig. 6 we present the EE and occupation number N ( j) for a chain of size L = 240 at a point in the IC phase, h = 0.7 and U = 1.2. For this set of parameters the ground state happens to belong to the odd sector.
One signature of this phase is the presence of oscillations in the EE and occupation number. Another feature of this phase is that one can use the CC formula, Eq. (26), for the EE to fit the data. Previously this has been done to distinguish the IC phase from the topological and the trivial phases of a Z 3 parafermionic chain [44].
For fitting we usually drop the first and the last 20 sites because there are clear edge and finite size effects in the both EE and occupation number, see Fig. 6. Using the CC formula, Eq. (26), to fit the data, we get c = 1.03. The fitted curve is also plotted in Fig. 6 as a solid line. Hence the essential features of the EE in the IC phase can be captured by the CC formula with c = 1. One would, of course, get a better fitting by considering the fluctuations on top of the CC formula, but we did not attempt this. We checked that the gap within each parity sector scales as 1/L, consistent with the fact that the main features of the EE are well described by the CC formula.
Pinpointing the transition from the topological phase to the IC phase is rather a hard task and this transition was conjectured to be of Kosterlitz-Thouless type [9,11]. We looked at the energy and its derivatives and we did not find any signature of a first or second order transition. We set the transition point by the energy difference of the order of 10 −3 between the ground state in the two parity sectors and by requiring c = 1 for fitting the CC formula to the DMRG data. Note that for a point in the topological phase close to the transition the EE can also be fitted well to the CC formula with central charge c < 1. One can distinguish the topological phase from the IC phase by looking at larger system sizes. In the topological phase, for large enough system size we do see the saturation of EE, however, in the IC phase one would obtain c = 1 also for larger system sizes.
Another qualitative check for distinguishing the IC phase is the presence of oscillations in the occupation number. This can be verified by looking at the bulk of the chain, comparing Figs. 5(b) and 6(b).
From Fig. 4 it is evident that the parity of the ground state in the IC phase changes. From this figure one can see that there are several crossings between the ground states of the two parity sectors, i.e., E +1 0 and E −1 0 . One should keep in mind that both the precise point where the these crossings happen, as well the parity of the ground state in this region, depends on the system size. Nevertheless, there are two important crossings, namely the last crossing at U 1.41 and the nextto-last crossing at U 1.28. We will later discuss the region between these two crossings in detail in Sec. IV C. Note that for h = 0.7 the IC phase ends at the next-to-last crossing at U 1.28 in Fig. 4.
To verify that level crossings occur we looked at the ground state energy and its derivatives with the steps of 10 −2 along the U axis. The first and the second derivatives of the ground state energy with respect to U are shown in Fig. 7. There are discontinuities in ∂E G ∂U which give rise to the peaks in − ∂ 2 E G ∂U 2 .

085125-6
(a) Energy with respect to the ground state energy, E P n − E G , of the ground state (n = 0) and the first excited state (n = 1) in both parity sectors are plotted as function of U for h = 0.05. We have looked at the similar quantities for different system sizes and the discontinuities are always present in the first order derivative. These are clear signatures of level crossings. We note that both the number of crossings, as well as their location, depends rather strongly on the system size.
To address the question about the presence of the floating phase in the ANNNI model [Eq. (4)] in the highly frustrated limit, i.e., κ 1, we can use our knowledge of the Kitaev-Hubbard chain. Because the Kitaev-Hubbard chain is dual to the ANNNI model, the presence of the IC phase in the fermionic incarnation corresponds to the presence of the floating phase in the ANNNI model. From Eq. (9) we can see that for U 1 and small h, we will be in the regime of B J 1 , κ 1, which is exactly the regime where the presence of the floating phase in the ANNNI model was under debate. In Fig. 8(a) we plot the energy of the ground state (n = 0) and the first excited state (n = 1) in both parity sectors for a chain of length L = 240, at h = 0.05. This corresponds to κ 20. The situation is similar to h = 0.7, as plotted in Fig. 4. For low enough U the ground state is doubly degenerate and the states have different parities, i.e., a topological phase. For large enough U , the ground state is also degenerate and the states have the same parity, i.e., P = +1 in this case.
In Fig. 8(b) one observes a discontinuity in ∂E G ∂U and the corresponding peak in − ∂ 2 E G ∂U 2 for h = 0.1 (which corresponds to κ = 10) at U = 1.05, as well as a broad feature at U = 0.99. In addition, it is also possible to fit the EE with the CC formula, resulting in c = 1 for a small region in U , even for h = 0.05. Finally, one observes an onset to the oscillations in the occupation number, that is characteristic for the IC phase. Given these features, we conclude that the IC phase in the Kitaev-Hubbard chain continues down to the XY critical point. Therefore, we also conclude that the floating phase in the ANNNI model is present for arbitrarily large frustration.

C. The excited-state CDW phase
In our numerical studies we found that the ground state for h = 0.7 and U 1.41 has CDW order and a low EE, which is constant in the bulk, that is, independent of the subsystem size. By looking at the energy levels in Fig. 4 we see that there is a region, 1.28 U 1.41, just before the CDW phase in which the the ground state is in the odd sector for our default system size, L = 240. We believe that this part of the phase diagram, colored light green in Fig. 1, is a phase which had been missed before and has been considered to be part of the CDW. We will discuss the properties of this phase thoroughly and show how the behavior of the model in this part of the phase diagram depends on the system size. Moreover, we show that the transition point from the new phase to the CDW phase can be actually controlled by tuning the chemical potential on the first and the last sites.
We present the EE and occupation number for a generic point (h = 0.7 and U = 1.325) in the new phase in Fig. 9. It is clear the EE of the ground state grows as a function of subsystem size. The particle number has a long wavelength oscillation of twice the system size, accompanied by a π phase shift for neighboring sites. These observations show that the ground state does not have a conventional CDW ordering for our system size. In fact, the properties of the ground state resemble the properties of low-lying excited states in the CDW phase. We therefore refer to this phase as the "excited-state charge density wave" (esCDW) phase. In addition, we found that the model is gapless by looking at the first (n = 1) and the second (n = 2) excited states in the odd sector. The energy difference between them, E −1 n = E −1 n − E −1 0 , approaches zero as L −1.80 for n = 1, 2. The data for n = 1 and the fit are shown in Fig. 10. This shows this phase is a gapless phase with the dynamical critical exponent z 1.80 and hence cannot be described by a conformal field theory. We also found that there is a finite gap to the lowest energy state in the even sector. In Fig. 11 we plot E +1 0 − E −1 0 as a function 1 L which shows the existence of the gap.
An intriguing feature of this region of the phase diagram is that its features depend crucially on the system size. The ground state belongs to the odd sector for system sizes that are a multiple of four, L = 4n. For L = 4n + 2 however, the 085125-7  ground state shows the same behavior as described above, but it belongs to the even sector. For an odd number of sites, L = 4n + 1 (or L = 4n + 3), however, the model is gapped and has a unique ground state with the CDW ordering in the even (odd) sector. In Fig. 12 we present the two lowest states in each parity sector for L = 243. One should compare this plot with Fig. 4; one observes that after the IC phase there is only one phase. We checked that the ground state has a low EE and the occupation number showed CDW ordering. However, there is one important difference with the CDW phase for even system sizes. Namely, in Fig. 4 we see that as soon as one enters the CDW phase the ground state is doubly degenerate and both of them have the same parity, i.e., P = +1. For odd number of sites, however, the ground state is unique, with a finite gap to the lowest excited state, which has opposite parity. This excited state also has a low EE and CDW ordering. For large U , the gap between these states goes to zero, which we verified by considering values for U up to U ∼ 10. We conclude that for odd system sizes, upon increasing U , the IC phase gives way to a CDW phase (with a unique ground state) immediately.
We should comment on one more feature in Fig. 12. The energy of the ground state in the even sector (the green dots) makes a rather large jump between U = 1.5 and U = 1.525. We thoroughly checked that this jump is not an artifact of the DMRG algorithm. First of all, the states are clearly converged, because the observables all have smooth profiles (as a function of the site number), and the energy does not change if the number of sweeps is increased. In addition, this feature is present for all system sizes and parameters in this region, it does not go away upon slightly varying these parameters.
In the region where the esCDW phase is observed, the behavior of the model does not only depend on the (parity of the) system size, but also on how one implements the magnetic field term. As we mentioned before we have implemented the magnetic field term on bonds, see Eq. (23). We also considered the model with the same magnetic field on all sites, namely by using the following term in the fermionic incarnation in DMRG: As expected, we did not seen any difference between the two ways of implementing the magnetic field in the topological and the IC phases. The IC phase for h = 0.7 still ends around U 1.28, however, the esCDW phase extends up to a larger value of U , namely U 1.68, as can be observed by comparing Figs. 4 and 13, in which we plot the two lowest energies for both parity sectors as a function of U . Therefore the esCDW phase is sensitive to both the system size as well as the chemical potential on the first and the last sites. To be sure that these results are not due to convergence issues and unluckily chosen system sizes, we checked the results for various smaller system sizes as well, leading to identical results.
We should make two important remarks, before closing our discussion about the esCDW phase. First of all we should mention that if one uses infinite DMRG to study this phase, one would conclude that it belongs to the CDW phase [12]. This is due to the fact that the ground state energy per unit site (or per bond) for the odd number of sites is lower than the one for even number of sites by 10 −3 in the units of hopping energy.
Another issue concerns the presence of the esCDW phase for small h. We clearly see the presence of the esCDW phase for h = 0.1. For smaller h, however, it is hard to draw firm conclusions. As we discussed in Sec. IV B we see the signatures for the crossing between the energy levels and hence discontinuities in the first order derivative of the energy with respect to U . However, to conclude that both the IC and the new phase are present for small h we need to see two crossings, one which corresponds to the transition from the IC phase to the new phase and another one for the transition from the new phase to the CDW phase. With current accuracy and system sizes we only see one crossing. Because the ground state EE is consistent with the CC formula for c = 1 we conclude that at least the IC phase is present. However, the the width of these two phases shrinks as we get closer to the XY critical point and distinguishing them becomes harder and needs larger system sizes and more accurate treatments. We expect, however, that the new phase also extends to h = 0, i.e., down to the XY critical point, as indicated in the phase diagram, Fig. 1.

D. The charge density wave phase
By further increasing U , as we concluded above, we enter the regime where the ground state has CDW ordering. For L = 240 and implementing the h term as in Eq. (23), the EE and occupation number for h = 0.7 and U = 1.5, are presented in Fig. 14. In this case the two ground states belong to the even parity sector. In general, if we implement the h term on bonds as we have in Eq. (23), the ground state is doubly degenerate. For L = 4n (or L = 4n + 2), where n is a large integer, the two ground states belongs to the even (odd) sector. For L odd, however, one ground state has even parity, the other odd. This can easily be understood in the limit where U is infinite, so that one ignore the hopping and pairing terms in the Hamiltonian (which is then diagonal). As is shown in Fig. 14 the EE for both ground states is saturated to a low value and the occupation number shows the pattern which one expects in the CDW phase, i.e., (1010 · · · ) and (0101 · · · ).
If one instead implements the magnetic field on the sites, as in Eq. (28), the ground state for an even number of sites is still twofold degenerate. However, for an odd number of sites, the ground state is unique, and the first excited state has energy 2|h|. Both of these states have the usual characteristics of CDW states.

V. ON THE EXCITED STATES AND FINITE TEMPERATURE PROPERTIES
In this section we study some properties of the model which are related to the full many-body spectrum of the model and its properties at finite temperature.

085125-9
(a) Entanglement entropy as a function of subsystem size.  First of all we address the nature of the zero mode in the topological phase of the model. Usually in the studies on the topological phases the emphasis is on the degeneracy of the ground state. There are, however, a few cases where one can say more. For example in the Kitaev chain [5] which can be written in terms of free fermions, the presence of the zero mode guarantees that the full many-body spectrum is doubly degenerate [45]. Given an eigenstate in the many-body spectrum of the Kitaev chain with a given fermionic parity, one can construct another eigenstate with opposite parity. This can be done by looking at the occupancy of the zero mode. If it is empty, one can consider a state with an occupied zero mode and vice versa. Since filling the zero mode does not cost any energy (or an energy that is exponentially small in system size), the states are degenerate but with opposite parity. Therefore full many-body spectrum is doubly degenerate and the zero mode is called a strong zero mode. Another, more nontrivial, example is the XY Z chain for which Fendley found an operator which commutes with the Hamiltonian up to exponentially small corrections in the system size [46]. The operator, however, is only normalizable in the topological phase, and hence the model hosts a strong zero mode and full many-body spectrum degeneracy in this phase.
There is also a conjecture for general interacting models [47] which states that if one starts from a noninteracting model within a topological phase, where we know there is a strong zero mode, one can construct a strong zero mode for the general interacting model providing that the system does not go through a phase transition or has level crossings. Therefore the full many-body spectrum would be degenerate all over the topological phase and there is equivalence between different parity sectors. In another study on nonintegrable interacting models such as the Kitaev-Hubbard chain [48], it was argued that a strong zero mode is not present. Below we present an analytic argument that excludes the presence of a strong zero mode for a large part of the topological phase. We also study the edge magnetization in the region where our argument does not apply.

A. Absence of strong zero mode in some regions within the topological phase
As we mentioned previously, the Hamiltonian in Eq. (8) is solvable for U = 0 [24] and h = 0 [27], because it can be reduced to a free fermionic model in those cases. Based on the exact results we know that there is a strong zero mode all along the U axis, except for the critical points U = ±1 and along the h axis for |h| < 1. Thus it is natural to ask whether a strong zero mode survives throughout the topological phase.
Here we present an argument to exclude the possibility of having a strong zero mode for some regions of the topological phase.
Consider the solid blue line [h(s), U (s)] in Fig. 15, which lies in the topological phase. Hence the Hamiltonian H (s) = H[h(s), U (s)] has a doubly degenerate ground state. Actually, the line we chose is the Peschel-Emery line [30], for which the ground state degeneracy is exact, the splitting is not exponentially small in system size. We want to prove that the full many-body spectrum of H (s) is not doubly degenerate 085125-10 all along this line. To do so we examine the behavior of the highest excited state(s). This can be done by noting that for any Hamiltonian H, the ground state(s) of −H are the highest excited state(s) of H itself. Let us consider −H (s), Now we perform an on-site rotation, with which we get Therefore we can track the behavior of the highest excited state(s) of H (s), by looking at the ground state(s) of H[h(s), −U (s)]. The dashed blue line in Fig. 15 is the mirror of the solid blue line with respect to the h axis. Following this dashed line from h = 0 and U = 0 till P * , the ground state of H is doubly degenerate. By passing P * , however, we enter the trivial phase and hence the ground state becomes unique. We can use this information and interpret as the behavior of the highest excited state(s) of H (s), namely the Hamiltonian along the solid blue line. Doing so, we learn that by passing the point P, the highest excited state is not doubly degenerate anymore. As a result there cannot be a strong zero mode and hence no equivalence of the parity sectors beyond the point P along the blue solid line. Having this argument, we can draw the red dashed line which gives a bound on the region where the system potentially hosts a strong zero mode.
Below we present numerical results which show that in part of this region, the edge magnetization survives for a long time, both for low as well as high temperature, which is consistent with the presence of a strong zero mode. We cannot rigorously prove its existence, and it might well be that a strong zero mode is only present when the model reduces to a free-fermionic model.

B. Time dependence of the edge magnetization
In the following we consider the time dependence of the edge magnetization [48], to see if the system has a strong zero mode. Long time coherence of edge magnetization at infinite temperature is a signature of a strong zero mode [48]. We will also consider finite temperature, in which case the time-dependent edge magnetization at temperature T is given by where j are the energies of the system. We will refer to the time in units of 1/J, and temperature in units of J, where J = 1. We considered system sizes up to L = 16.
We first look at small h = 0.1, and a rather strong interaction |U | = 0.5, but still in the region in Fig. 15 where the analytic argument cannot exclude a strong zero mode. In these two cases, we find that the edge magnetization at very large T = 10 3 , so effectively infinite temperature, survives up to t ∼ 100, before dropping to zero, with only little variation upon increasing the system size from L = 10 to L = 14. with an amplitude of 0.1. For these parameters, it is clear that the system does not have a strong zero mode.
For h = 0.1 and |U | = 0.1, the situation is rather different, as we show in Fig. 16 for h = 0.1 and U = −0.1, where we plot |A T (t )| as a function of t on a log scale for temperatures T = 10 −3 , 1, and 10 3 , with blue, green, and red symbols, respectively. The magnetic field was implemented on the bonds, see Eq. (23).
We observe that the edge magnetization survives longer upon increasing the system size, for all temperatures. Interestingly, the edge magnetization survives longer for higher temperatures. However, at high temperature, it drops to zero, while for low temperature, the edge magnetization rapidly oscillates at long times. The high temperature behavior is consistent with the behavior seen for the XY Z chain [48], for which it was proven that the system has a strong zero mode [46]. In Fig. 17 we plot the T = 1000 edge magnetization for U = 0.1, h = 0.1 and even system sizes L = 8-16. We observe that the time that the edge magnetization survives grows exponentially with system size, at least for the system sizes that we could consider. So based on this data, we cannot rule out the presence of a strong zero mode for these parameters.
We note that for U < 0, the different ways of implementing the magnetic field give very similar results. However, for U > 0 at low temperature, the time at which the edge magnetization drops to zero for the first time can differ by as much as two orders of magnitude. For which implementation of the magnetic field this happens first, depends on the system size. That the time dependence of the edge magnetization depends on the value of the magnetic field on the edges is of course not so strange, but from the numerical results we obtained, a clear picture does not emerge.
For larger magnetic fields and moderate interaction, h = 0.7, and |U | = 0.1, the edge magnetization for infinite temperature survives to t ∼ 100, again with little dependence on the system size. Hence, for these parameters, there is clearly no strong zero mode. For h = 0.7 and U = 0.5, our analytical argument excludes a strong zero mode, because the point h = 0.7, U = −0.5 is in the trivial phase. Indeed, we found that at high temperature the edge magnetization drops sharply at t ∼ 1, while the edge magnetization for low temperature, probing the ground state, survives up to long times.

VI. DISCUSSION
We presented a detailed study the phase diagram of the Kitaev-Hubbard chain, that is the Kitaev chain in the presence of a nearest-neighbor density-density interaction. In the case of attractive interactions, the model exhibits a topological phase and a trivial phase, the same phases which also appear in the Kitaev chain. For repulsive interactions, however, the phase diagram is more interesting, with an IC phase, an esCDW and a CDW phase, in addition to the trivial and topological phases.
The bosonic incarnation of the interacting Kitaev chain is an interacting transverse field Ising chain, which is dual to the quantum ANNNI model. There was an open question in the quantum ANNNI model about the presence of the floating phase for high value of frustration. Using fermionic terminology, this translates to presence of the IC phase in the Kitaev-Hubbard chain for small values of the transverse field h that is close to the XY critical point. Using DMRG we see the presence of the IC down to h = 0.05 which corresponds to κ 20, although it becomes quite narrow. Therefore we strongly believe that the IC phase continues to h = 0, which is the limit of infinite frustration.
The IC phase is a gapless phase, for which the EE can be fitted to the CC formula with a value of the central charge c 1. However, the EE as a function of subsystem size shows additional oscillations, on top of the CC behavior. It would be interesting to find a good ansatz for the oscillations, because a quantitative knowledge of these oscillations can shed additional light on the nature of the IC phase. It would also be interesting to see how well they correlate with the oscillations observed in the site occupation number.
The esCDW phase, we believe, was missed in previous studies. This phase only appears for an even number of sites and its ground state looks like an excited state in the CDW phase. In this region, the model is gapless, with a dynamical critical exponent z 1.8. The EE in this phase grows as a function of subsystem size, but cannot be fitted to the CC formula, which is consistent with the dynamical critical exponent z > 1. For odd system sizes, the model shows a CDW phase, with a unique ground state. It is quite intriguing that there is finite region in the phase diagram which responds to the (parity of) the number of sites. In infinite system size studies, one observes a normal CDW phase [12]. This is consistent with our findings, because the ground state energy per site is slightly lower, of order 10 −3 , for odd system sizes. Because the energy difference is so small, the nature of this phase will in practice be determined by other perturbations, for instance additional interaction terms.
For repulsive interactions, we found that there are four phases which emerge from the XY critical point, namely the topological phase, the IC phase, the esCDW phase, and the CDW phase. Capturing these phases and finding the relevant operator(s) for each phase using for instance bosonization will be left for future research.

085125-12
Finally, we considered the nature of the zero mode in the topological phase. We presented an analytical argument that excludes a strong zero mode for a large part of the topological phase. To address the remaining region, we calculated the time dependence of the edge magnetization [48]. In the remaining region these results rule out a strong zero mode in the case that either h or |U | is large. For small h and |U |, we cannot exclude a strong zero mode. For U > 0, the time dependence of the edge magnetization does not show consistent behavior as a function of temperature and the implementation of the magnetic field. This warrants further study of the edge magnetization in this model.