SciPost Phys. 14, 152 (2023)

Topological and quantum critical properties of
the interacting Majorana chain model
Natalia Chepiga1 and Nicolas Laflorencie2
1 Kavli Institute of Nanoscience, Delft University of Technology,
Lorentzweg 1, 2628 CJ Delft, The Netherlands
2 Laboratoire de Physique Théorique, Université de Toulouse, CNRS, UPS, France

Abstract
We study Majorana chain with the shortest possible interaction term and in the presenceof hopping alternation. When formulated in terms of spins the model corresponds to thetransverse field Ising model with nearest-neighbor transverse and next-nearest-neighborlongitudinal repulsion. The phase diagram obtained with extensive DMRG simulations
is very rich and contains six phases. Four gapped phases include paramagnetic, period-2
with broken translation symmetry, Z2 with broken parity symmetry and the period-2-Z2phase with both symmetries broken. In addition there are two floating phases: gaplessand critical Luttinger liquid with incommensurate correlations, and with an additional
spontaneously broken Z2 symmetry in one of them. By analyzing an extended phase diagram we demonstrate that, in contrast with a common belief, the Luttinger liquid phasealong the self-dual critical line terminates at a weaker interaction strength than the end
point of the Ising critical line that we find to be in the tri-critical Ising universality class.
We also show that none of these two points is a Lifshitz point terminating the incommensurability. In addition, we analyzed topological properties through Majorana zero modes
emergent in the two topological phases, with and without incommensurability. In the
weak interaction regime, a self-consistent mean-field treatment provides a remarkableaccuracy for the description of the spectral pairing and the parity switches induced by
the interaction.
Copyright N. Chepiga and N. Laflorencie.
This work is licensed under the Creative Commons
Attribution 4.0 International License.
Published by the SciPost Foundation.

Received 06-12-2022
Accepted 11-04-2023
Check forupdatesPublished 13-06-2023
doi:10.21468/SciPostPhys.14.6.152

Contents
1 Introduction
1.1 Generalities
1.2 The interacting Kitaev-Majorana chain model
1.3 Paper outline

3
3
3
5

2 Main results and phase diagram
2.1 Different phases
2.1.1 Parity (Z2 ) broken antiferromagnetic (topological) order
2.1.2 Disordered paramagnetic phase
2.1.3 Floating phases

6
6
6
6
7

1

SciPost Phys. 14, 152 (2023)

2.2

2.1.4 Period-2 orders
Phase transitions and multicritical points

7
7

3 Floating phases
3.1 The boundaries of the floating phase
3.2 Floating-1 vs Floating-2
3.3 The eight-vertex point

8
8
9
11

4 h = 1 line
4.1 The boundaries of the Luttinger liquid
4.2 Continuation of the Ising transition and its end point
4.3 The persistency of incommensurability

12
12
14
14

5 Majorana zero modes
5.1 Definition and properties of strong Majorana zero modes
5.2 Non-interacting Kitaev chains
5.2.1 Transverse-field Ising chain
5.2.2 Incommensurability and parity switches for Kitaev wires
5.3 Interacting case: DMRG results
5.3.1 Z2 antiferromagnet at small g
5.3.2 Z2 period-2 antiferromagnet at large g
5.4 Self-consistent mean-field theory
5.4.1 Mean-field Hamiltonian
5.4.2 Numerical results
5.4.3 Mean-field phase diagram

15
15
16
16
16
18
18
19
20
20
20
21

6 Discussions and conclusions

22

A Useful transformations
A.1 Duality
A.2 Jordan-Wigner and Majorana fermions

24
24
24

B Strong zero modes for the non-interacting Kitaev model
B.1 SZM construction: iterative procedure
B.2 Surface magnetization
B.3 Finite size parity gap and parity switches

24
24
26
26

C Self-consistent mean-field

26

D Ising transition for small g

27

E Comparison to the previous study at h = 1

28

References

29

2

SciPost Phys. 14, 152 (2023)

1

Introduction

1.1

Generalities

Models of strongly-correlated low-dimensional systems attracted a lot of attention in the past
decades [1]. Frustrated spin chains, models of spinless fermions and hard-boson models
have been studied intensely over the years and have lead to a number of fascinating results. One of them is the emergence of Majorana edge states predicted by Kitaev [2] in achain of spinless fermions - an effective model of p-wave superconductors. This theory predictions stimulated impressive experimental activity [3–9] motivated by the potential usageof Majorana zero modes in quantum computing [10–12]. By construction the Kitaev chain
describes non-interacting Majorana fermions, therefore shortly after the model has been proposed it has also been extended to include interactions of various forms between the Majorana
fermions [13–19]. Interacting Majorana fermions have also been theoretically studied in thepresence of quenched disorder [20–26], in particular more recently in the context of manybody localization physics at high energy [27–31].

1.2

The interacting Kitaev-Majorana chain model

Here we focus on low-energy properties for the following disorder-free microscopic spin Hamiltonian on a one-dimensional (1D) chain:

H=

X


Jσ xj σ xj+1 − hσzj + gz σzj σzj+1 + g x σ xj σ xj+2 ,

(1)

j
x,z

where σ j are Pauli matrices on site j. By means of Jordan-Wigner transformation this model
can be rewritten in terms of interacting (Dirac) fermions (see Appendix A.2)
X 




t c †j c j+1 + h.c. + ∆ c †j c †j+1 + h.c. + 2hn j
H=
j

+ gz

X

+ gx

X



1 − 2n j+1



1 − 2n j+1

1 − 2n jj

c †j − c j †c j+2 + c j+2 ,(2)jwith equal hopping and pairing terms t = ∆ = J, and interaction strengths gz and g x . Wheng x,z = 0 we immediately recognize the celebrated transverse-field Ising (TFI) chain model [32,33] in Eq. (1), equivalent to the so-called Kitaev chain model as given by the first line of Eq. (2).
The two other terms encode the interactions: the second line a nearest-neighbor densitydensity repulsion for gz > 0 (attraction for gz < 0), while the last one can be understood as adensity-assisted hopping and pairing term at second neighbor distance.
This form of interaction takes a more symmetric expression in the Majorana language (seeAppendix A.2)X X
H = −iJ b j a j+1 − ha j b j −gz a j b j a j+1 b j+1 + g x b j a j+1 b j+1 a j+2 ,
(3)
j

j



where we have introduced two Majorana fermions at each site: a j = c †j +c j and b j = i c †j − c j .
A sketch of the model is given in Fig. 1 for both spin and Majorana representations. One caneasily see that tuning the transverse field of the Ising model in Eq. (1) introduces a bondalternation of the kinetic term in the Majorana chain, see Fig. 1.
3

SciPost Phys. 14, 152 (2023)

(a)

gx
h

σ1

J
gz

σ3

σ2

σL

gx
a1 h b1 J a2

b2

a3

gz

b3

aL

bL

(b)

Figure 1: Schematic picture for the spin Hamiltonian (a), and its Majorana fermion
representation in (b), both with open boundary conditions. In this work, we fix J = 1
and g x = gz = g.
Various types of interactions have been studied over the past decade, but clearly most of
the attention has been devoted to models with density-density interactions gz ̸= 0 [15, 18, 19,
34–41]. In particular, the emergence of a stable Luttinger liquid phase [19], the eight-vertex
criticality [41] and the exact disorder line [18] have been reported. In contrast, models withg x ̸= 0 have not retained much consideration, except a few cases involving correlated hoppingof fermions, as studied recently in Refs. [42–44].
In this paper we consider symmetric interactions gz = g x = g, which display a KramersWannier self-duality, see Appendix A.1. In the rest of the paper we will set J = 1 and focus on
the repulsive g > 0 regime, leaving the attractive case for further studies. In Refs. [16, 17],
using field-theory and density matrix renormalization group (DMRG) methods Rahmani et al.
have achieved a detailed study of symmetric interacting problem for the special case J = h.
Let us briefly summarize their results for the repulsive situation that we want to consider
in this work: For g ≲ 0.28 the system can be described by the Ising critical theory; when
0.28 ≲ g ≲ 2.86 it was suggested that the effective critical theory can be characterized by acentral charge c = 3/2 and corresponds to the combination of the Ising and the Luttinger liquidcriticalities; the latter is accompanied by an emergent U(1) symmetry in agreement with thelater analysis by Verresen et al. [19]. The point g ≈ 2.86 has been identified as a generalizedcommensurate-incommensurate transition beyond which the system is gaped with four-folddegenerate ground-state.
The goal of the present paper is to study an extended phase diagram based on which werediscuss some of the conclusions drawn in the previous works. Building on DMRG simulations [45–49], we map the global phase diagram of the extended interacting Kitaev-Majorana
chain model Eq. (1) with g = g x = gz as a function of a transverse field h and interaction
strength g, see Fig. 2. We perform simulations with two-site DMRG on systems with up to
N = 2000 sites with open boundary conditions (OBC) and up to N = 200 with periodic boundary conditions (PBC). We achieve the convergence by performing up to 8 sweeps keeping up
to 3000 states and discarding singular values below 10−8 . When using OBC we either fix
them by polarizing the edge spins along the chosen direction or keep them free. The former isused to compute Friedel oscillations profiles in the floating phases and to study quantum phasetransitions. The latter is chosen when we study the emergence of Majorana edge states and associated exact zero modes. We also complement the DMRG simulations using a self-consistent
Hartree-Fock treatment of the problem in the limit of weak interaction, see Appendix C forsome details.

4

SciPost Phys. 14, 152 (2023)

3.5

TCI

3

2

Ising

g

2.5

period-2- period-2

1.5

Kosterlitz-ThoulessS
Floating-2

Floating-1

c=3/2

1

0.5

0

IC

IC
C

0

Ising

8-vertex
0.5

1MParamagnetic

C
1.5

2

h

Figure 2: Phase diagram of the interacting Majorana chain model Eq. (4) as a function of coupling constants g and h, obtained from DMRG simulations. It contains
four gapped phases: Z2 , paramagnetic, period-2 and Z2 -period-2; and two criticalfloating phases in Luttinger liquid universality class. See main text, Section 2.1 fora brief description of each phase. The model is self-dual for h → 1/h and g → g/h.
The floating phases are separated from the gapped ones by the Kosterlitz-Thoulesstransitions (red lines) with indicated critical values of the Luttinger liquid parameters
K c . The multicritical point at h = 0 and g ≈ 0.4105 is in the eight-vertex universality class. Blue lines are Ising transition that terminates at g ≈ 3 with the tri-critical
Ising point (blue star). For 0.29 ≲ g ≲ 1.3 the Ising transition is superposed with
the Luttinger liquid phase resulting in a critical line with central charge c = 3/2.
Dotted black line states for the disorder line above which the dominant wave-vector
is incommensurate.

1.3

Paper outline

The rest of the paper is organized as follows. In section 2 we overview the phase diagram
briefly discussing the properties of each phase and the nature of the quantum phase transitions between them. In section 3 we discuss in more details the floating phases - Luttinger
liquid phases with incommensurate correlations and locate the boundaries of these phasesthat corresponds to two Kosterlitz-Thouless phase transitions. We then discuss the multicritical point along h = 0 line that belongs to the universality class of the eight-vertex model.
Equipped with the understanding of the extended phase diagram we revise the nature of thecritical lines along h = 1 in the section 4. In particular we will provide numerical evidence
that there is no generalize commensurate-incommensurate transition predicted in Ref. [17],5SciPost Phys. 14, 152 (2023)

the Ising transition persist beyond the Luttinger liquid phase and terminates at the tri-critical
Ising point. In section 5 we discuss the Majorana zero modes that appear in the Z2 regimes,
with or without incommensurability. DMRG results are successfully compared with a selfconsistent Hartree-Fock treatment in the limit of weak interaction. We summarize our results
and put them in perspective in section 6.

2

Main results and phase diagram

We investigate the phase diagram of the interacting Majorana chain model defined by Eq. (1)with g x = gz = g and J = 1 fixed, as a function of the transverse field h and the interaction
coupling constant g, as rewritten here
X


H=
σ xj σ xj+1 − hσzj + g σzj σzj+1 + σ xj σ xj+2 .
(4)j

The phase diagram obtained from extensive DMRG calculations is presented in Fig. 2 andconsists of six main phases. Below we provide a short summary of our results, describing thedifferent regimes in Sec. 2.1 and the associated phase transitions and (multi) critical points in
Sec. 2.2. More details and extended discussions are provided in the corresponding sections of
the paper.

2.1

Different phases

Before proceeding to a list of various phases let us point out a very important property ofthe phase diagram. The model defined by Eq. (1) up to boundary terms obeys the KramersWannier self-duality (see Appendix A.1) with h → 1/h and g → g/h. It immediately defines
h = 1 as a very special line in the phase diagram where the transitions between each pair ofdual phases take place. It turns out that this critical line is associated to a spontaneous parity
symmetry breaking for one of the side of the transition, depending on the interaction strength.
2.1.1

Parity (Z2 ) broken antiferromagnetic (topological) order

Realized for small g and h < 1 with spontaneously broken parity, this gapped phase is topologically non-trivial with a two-fold (parity) degeneracy with emergent Majorana edge states.
Note that this phase has no topological interest for the (Ising) spin degrees of freedom forwhich the so-called “topological phase” there boils down to a more conventional magnetic
order.
• Z2 -C: Commensurate region of the Z2 phase with zero modes showing a vanishing (parity) gap, exponentially suppressed with the system size.
• Z2 -IC: A region of the Z2 phase with incommensurate short-range order and with zeromodes that due to incommensurability are exact at some points even for the finite lengthof the chain.
2.1.2

Disordered paramagnetic phase

This gapped phase occurs at small g and h > 1. The ground-state is short-ranged correlated,and corresponds to all spins polarized along the field h.
• Paramagnetic-C: Paramagnet with commensurate short-range order• Paramagnetic-IC: Paramagnet with incommensurate short-range order showing the same
incommensurate vector as its dual phase Z2 -IC.
6

SciPost Phys. 14, 152 (2023)

2.1.3

Floating phases

These are gapless Luttinger liquid phases, with quasi-long-range incommensurate order. The
corresponding critical theory is characterized by a central charge c = 1. The floating phases
are stabilized by an emergent U(1) symmetry [17, 19]. In the Floating-2 phase realized forh > 1 the parity symmetry Z2 is spontaneously broken.
2.1.4

Period-2 orders• Period-2 (for h < 1): Gapped phase with spontaneously broken translation symmetry
and a two-fold degenerate ground-state. The phase is characterized by an antiferromagnetic order along z with |〈σzi − σzi+1 〉| ̸= 0, and on top of it there is an incommensurate
short-range order.
• Z2 -period-2 (for h > 1): Gapped phase with spontaneously broken translation and parity
symmetries. For chains with even number of sites there are Majorana zero modes.
Both phases, at least big portions of them, have incommensurate short-range order, however
the dominant wave-vector is very close to the commensurate value and, as we will see, approaches it with a zero slope. It means that for strong repulsion we cannot exclude the existence of the disorder line beyond which the correlations are commensurate. It also implies
that any finite-size result pointing towards commensurate wave-vector should be treated withcaution - it might take much larger system to detect a presence of incommensurability.

2.2

Phase transitions and multicritical points

The phase diagram contains a wide variety of quantum phase transitions:• Ising transition: Quantum phase transition between the Z2 and the paramagnetic
phases is in the Ising universality class along h = 1 and for g ≲ 0.29.
• Ising+LL: Inside the floating phase Ising transition is superposed with the Luttinger
liquid phase resulting in a critical line characterized by the central charge c = 3/2. Thistransition separates two critical floating phases one of which has broken Z2 symmetry.
According to our results Luttinger liquid terminates at g ≈ 1.3 marked in Fig. 2 with S,
for much smaller interaction strength than suggested in Ref. [17].
• Second Ising transition: The Ising critical line continues beyond the point S and terminates around g ≈ 3.
• Tricritical Ising point: The location of the end point of the Ising critical line
agrees within the error-bars with the location of the generalized commensurateincommensurate transition reported in Ref. [17]. However, we arrived to a different
conclusion regarding the nature of this end point: in Sec. 4 we provide numerical evidences that (i) the end point is in the tri-critical Ising (TCI) universality class, (ii) the
Luttinger liquid terminates at much weaker interaction strength and does not affect thenature of the end point, and (iii) the incommensurate short-range correlations persist
beyond this end point which is not a Lifshitz point separating commensurate and incommensurate regimes.
• Kosterlitz-Thouless transitions: The Luttinger liquid phase is stable against superconducting instability for K < 1/2 and against broken translation symmetry for K > 1/4.
When the Luttinger liquid exponent reaches these critical values, the system undergoes Kosterlitz-Thouless transition into the corresponding gapped phase. Note that in7

SciPost Phys. 14, 152 (2023)

the period-2 phase short-range order remains incommensurate and the transition belongs to the Kosterlitz-Thouless universality class by contrast to the Pokrovsky-Talapov
commensurate-incommensurate transition realized in a related model that takes placeat the same critical value of K c = 1/4 [19, 41].
• Multi-critical point M: The pairing operator becomes relevant and lead to a gapped Z2
phase for K > 1/2. The spin-flip operation always creates a pair of domain walls in aparamagentic phase and thus has the same critical value K c = 1/2. However, along theIsing critical line neither the pairing nor the spin-flip becomes a relevant perturbationand the Kosterlitz-Thouless transition to the floating phase takes place at K c = 1 that corresponds to free fermions. Two disorder lines (dotted black) separating commensurate
and incommensurate regions of paramagnetic and Z2 phases terminate at this point.
• Eight-vertex critical point: The multi-critical point located at g ≈ 0.41 and h = 0 is in
the universality class of the eight-vertex model [15, 41, 50, 51].

3

Floating phases

3.1

The boundaries of the floating phase

Let us first focus on the h < 1 side of the phase diagram. The stability of the Luttinger liquid
phase against superconducting pairing and spontaneous translation symmetry breaking due
to nearest-neighbor repulsion has been studied recently in Ref. [19]. It has been argued that
due to an emergent U(1) symmetry the Luttinger liquid phase is stable when the Luttinger
parameter lies within the interval 1/4 < K < 1/2. For K > 1/2 the pairing term becomesrelevant; for K < 1/4 the phase is unstable with respect to spontaneously broken translation
symmetry.
By choosing open boundary conditions we explicitly break translation symmetry and therefore can observe the floating phase with Friedel oscillations of local magnetization. In orderto reduce log-corrections we polarize first and last spins with strong boundary field. Then
according to the boundary conformal field theory, Friedel oscillations on a finite-size systembehave as [52]cos(q j)〈σzi 〉 ∝,(5)[(N /π) sin(πi/N )]K
where K is the corresponding scaling dimension that for σzi operator (equivalent to the local
density operator) coincides with the Luttinger liquid parameter, q is the incommensurate wavevector.1 By fitting Friedel oscillations one can get an extremely accurate estimate of both, theLuttinger parameter K and the incommensurate wave-vector q. Examples of such fits arepresented in Fig. 3 where the uniform part of the spin density 〈σzi 〉 has been subtracted, andthe fitting window is cted to avoid edge effects.
By keeping track of the Luttinger liquid exponent K we can now locate the boundaries of
the Floating-1 phase as shown in Fig. 4 (a). Note that due to exponentially slow opening of the
gap at the Kosterlitz-Thouless transition it is possible to extract an effective critical exponent
K on both sides of the transitions. We therefore associate the location of the transitions withthe crossing points of the obtained curve with the corresponding field theory predictions. An
important observation is that the wave-vector q is incommensurate on both sides of the floatingphase. Therefore the transition to the period-2 phase with spontaneously broken translation
1

Similar results can be obtained for the local pairing term 〈σ xj σ xj+1 〉 under appropriate boundary conditions.

8

SciPost Phys. 14, 152 (2023)

0.1

0.4

(a)

DMRG
ﬁt

0.05

0.2

0

0

-0.05

-0.2

-0.1

100

200

300

400

(b)

500

600

-0.4

700

100

200

300

400

500

600

700

Figure 3: Examples of the Friedel oscillations inside the Floating-1 phase obtainedon a finite-size system with N = 8012 sites with polarized boundary conditions.
Blue points are DMRG data, red points are the result of the fit with Eq. (5) (bluedots are completely hidden under red ones). Note that the uniform part of the spin
density 〈σzi 〉 has been subtracted, and the fitting window is restricted to the rangei ∈ [121, 680] in order to avoid edge effects.
Floating-1

0.8

period-2

Floating-1

1

(a)

0.7

period-2

(b)

0.8

0.6

0.6

0.5
0.4

0.4

0.2

0.3
0.2
0

0.5

1

1.5

2

0

2.5

0

0.5

1

1.5

2

2.5

Figure 4: (a) Luttinger liquid parameter K and (b) incommensurate wave-vector q
along the vertical cut at h = 0.8. The system is in the Luttinger liquid phase for
1/4 < K < 1/2 (white region). It is clear that the wave-vector remains incommensurate on both sides of the floating phase. Blue diamonds and red squares are results of
the Friedel oscillations similar to those presented in Fig. 3; blue crosses state for thewave-vector extracted from the density-density correlation inside the gapped phases.
symmetry is also in the Kosterlitz-Thouless universality class in contrast to the PokrovskyTalapov commensurate-incommensurate transition to the period-2 phase with commensurate
long- and short-range orders [41].
Relying on the Kramers-Wanier duality, the location of the Kosterlitz-Thouless transition
for h > 1 can be deduced from the results obtained for h < 1. These results are in excellent
agreement with independent calculations performed for h > 1. In the paramagnetic phaseelementary excitations - spin flips - create domain walls in pairs and thus lead to the samevalue of the critical Luttinger liquid parameter K = 1/2. For strong interaction g the same
symmetry is broken (the translation) and thus, very naturally, the critical parameter is thesame K = 1/4.

3.2

Floating-1 vs Floating-2

Let us now focus on the similarities and differences between the two floating phases. Because of the self-duality of the model one might expect similar scaling the z component of2

Odd number of sites is chosen to realize symmetric boundary conditions in the period-2 phase with antiferromagnetic order along σz . Inside the floating phases the even-odd effect is relaxed due to incommensurate
wave-vector q but whenever we polarize edges in the longitudinal direction we always take the odd number ofsites for consistency.

9

SciPost Phys. 14, 152 (2023)

1

(a)

0.2

(c)

0.5

0.1

0

0
-0.1

-0.5

-0.2
0

200

400

600

800

400

1

KT
0.5

1
0.9

10

-4

10 0

10 1

0

10 2

(e)

1.1

0.5

10

1

800

1.5

1.2

10 -3

600

2

1.5
1.4
1.3

1

10 -2

200

(d)

1.5

(b)

10 -1

-1
0

1.2

-2

10 -1

1.4

0

Floating-2 period-2-Z2
1

1.5

2

2.5

Figure 5: (a) Friedel oscillations of the local magnetization σzi (blue) and of thexpairing operator σix σi+1
(magenta) inside the floating-1 phase (h = 0.8, g = 0.7).
Both have been centered around their mean values, the former has been re-scaledwith a non-universal pre-factor const ≈ 0.84. (b) Decay of the magnetization correlations 〈σzi σzj 〉 − 〈σzi 〉〈σzj 〉 as a function of distance at the two dual points inside
the Floating-1 (blue) and Floating-2 (red). (c) Local magnetization along x direction inside the Floating-2 phase with boundaries x-polarized. The amplitude of theoscillations A(σix ) is extracted as a difference between the largest and smallest values that 〈σix 〉 takes over an interval of length N /4 in the middle of the chain. (d)The amplitude A(σix ), extracted as in panel (c) along an horizontal cut at g = 1, isplotted upon approaching the Ising transition at h = 1. Extracted critical exponent
β ≈ 0.116 agree within 8% with the theory prediction for Ising transition β = 1/8.
Inset: the same plot in a log-log scale. (e) The amplitude A(σix ) along vertical cut at
h = 1.5 across the Kosterlitz-Thouless transition between Floating-2 and period-2-Z2 .
the correlations 〈σzi σzj 〉 − 〈σzi 〉〈σzj 〉 for a given point in the Floating-1 (h, g) and the pairingxxcorrelations 〈σix σi+1
σ xj σ xj+1 〉 − 〈σix σi+1〉〈σ xj σ xj+1 〉 in its dual point in Floating-2 (1/h, g/h),and vice versa. But it turns out that the connection is even stronger since at the same point,x
these two operators: local magnetization σzi and pairing σix σi+1
have the same scaling dimension K. This can be observed by comparing the profiles of the Friedel oscillations presented
in Fig. 5 (a). This connection also manifests itself in the same decay of the correlations, forinstance 〈σzi σzj 〉 − 〈σzi 〉〈σzj 〉 at the two dual points as shown in Fig. 5 (b). In both cases there
is a non-universal pre-factor in the correlators, see Fig. 5 (a-b).
What distinguishes the two floating phases is the Z2 symmetry broken in the Floating-2
phase. In order to demonstrate it we look at the order parameter that we associate with anamplitude of the x-component of the magnetization A(σix ). To detect Ising transition at small gx
one traditionally define the order parameter as |〈σix −σi+1〉|(with x-polarized boundaries, see
Appendix D). Here we have to adjust it to make it suitable for incommensurate correlations.
This is done by computing the maximal and minimal values of the x component of the local
magnetization in the middle of the chain and over an interval equal to the quarter of thechain. The example is shown in Fig. 5 (c). Next, we keep track of the amplitude A(σix )upon approaching the transition at h = 1, by fitting the data to ∝ (h − 1)β we extract thecorresponding critical exponent β ≈ 0.116. Our results are in a reasonable agreement with
β = 1/8 for the Ising transition. As mentioned above the pairing operator is relevant for
10

SciPost Phys. 14, 152 (2023)

1.2

-0.8

(a)
-1

2

(b)

1

-1.2

1

0.8

-1.4
-1.6

0.6

0
0

1

2

-1.8
0.4

-2
-2.2

0.2

-2.4
5

6

0

7

0.4

0.45

0.5

Figure 6: Identification of the eight-vertex universality class. (a) Finite-size scaling of
the amplitude of oscillations in the local magnetization. Separatrix is associated withthe critical point, the slope corresponds to the scaling dimension d. (a) Upper andlower bounds are identified. (b) Scaling of the order parameter with the distanceto the critical point. Symbols are DMRG data, lines are fits ∝ (g − g c )β . Inset:
Comparison between the extracted values of β and d for the interacting Majorana
chain (black point) with Baxter’s eight-vertex model prediction Eq. (7) (red line).
K > 1/2. Thus for h > J the pairing that appears in the Floating-2 and period-2-Z2 phases
never constitute a relevant perturbations and remains insensitive to the Kosterlitz-Thoulesstransition at large-g as shown in Fig. 5 (e).

3.3

The eight-vertex point

At h = 0 the floating phase collapses into a single point, thus the transition between the period2 and the Z2 phase becomes direct. When written in terms of Dirac fermions, the model of
Eq. (2) obeys particle-hole symmetry along h = 0. Based on the previous study the transition
between the Z2 and the period-2 phases along the particle-hole symmetry line is expected to
be in the universality class of the eight-vertex model [15, 41, 50, 51].
According to Baxter [53] all critical exponents of the eight vertex model depend on a single
parameter ρ:3
ν = π/(2ρ) , β = (π − ρ)/(4ρ) .
(6)In special cases, e.g. for the integrable model, the parameter ρ can be expressed in therms ofcoupling constants [51, 53], however in general this is not possible. Nevertheless, the special
relation between ν and β should hold for any value of ρ. It means that we can express one
exponent in terms of another, say β = β(ν) excluding ρ from the equation. In practice,
however, numerical estimate of the scaling dimension d = β/ν is way more accurate thanthe estimate of the correlation length critical exponent ν. Note that by contrast to β andν that regulates the order parameter and the correlation length as a function of distance tothe transition, the scaling dimension d can be extracted at the critical point and reflect thedependence on the system size. After a simple algebra, we got the following prediction for theeight vertex model:
β = d/(2 − 4d) .
(7)
3

In Baxter’s original notations this control parameter was called µ.

11

SciPost Phys. 14, 152 (2023)

In order to verify this prediction numerically we extract independently β and d. We extract
the scaling dimension d of the local order parameter that in the present case is given by thealternation between local magnetization on even and odd sites |〈σzi 〉 − 〈σzi+1 〉| - vanishing in
the Z2 phase and finite in the period-2 one. The critical point is associated with the separatrix
in a log-log scale and according to Fig. 6 (a) is located in the interval 0.41 < g c < 0.411. The
slope corresponds to the scaling dimension and lies in the range 0.389 ≲ d ≲ 0.41.
We measure the second critical exponent β by looking how alternation of the local magnetization vanishes upon approaching the transition. The results are presented in Fig. 6 (b).
Note that the shift of the location of the critical point with system size is typical and has been
observed in the previous study of the eight-vertex point [41]. But the critical exponent βchanges very little with the system size and we take the maximal difference as an estimate forthe errorbar. The agreement between obtained numerical results and the theory prediction of
Eq. (7) is spectacular and is presented in the inset of Fig. 6 (b). Note that there is no fittingor adjustment parameter and here the comparison between the numerics and the theory is
direct.

4

h = 1 line

The h = 1 line is self-dual, and it hosts very rich critical behaviors. The exploration of theextended phases away from this line helps us to get deeper insights into the criticalities realized
at h = 1. Although our numerical results in general agree with the numerical results provided
in Ref. [17] ( see the Appendix E for more comparison), by exploring the surroundings of thecritical line we arrived at a different conclusion regarding the nature of the critical phases
along h = 1 for large g. Below we provide a detailed overview of this very challenging line.

4.1

The boundaries of the Luttinger liquid

Let us come back to the original arguments regarding the stability of the floating phase. For theoperator that simultaneously create p domain walls the scaling dimension is given by p2 /4K,
the operator becomes relevant when the corresponding scaling dimension is smaller than 2.
This implies that the floating phase is stable against superconducting instability of the form†ci† ci+1+ h.c. when 22 /4K > 2, in other words, for K < 1/2. Similar argument applies in theparamagnetic phase: an elementary excitation in this case is a spin flip that cannot createa single domain wall but only a pair of them. Therefore, the transition to the Z2 or to theparamagnetic phase takes place when the Luttinger liquid exponent reaches its critical value
K c = 1/2 [19]. However, these arguments are not valid along h = 1 critical line because noneof the two operators become relevant - the theory remains critical and the parity symmetry is
not broken along this line. We thus associate the transition along h = 1 with the point where
K reaches the value K = 1 of free-fermion theory.
We extract the Luttinger liquid parameter K by fitting the profile of the Friedel oscillations.
We polarize edge spins in z direction and fit the data as shown in Fig. 3. Note that the Isingtransition with the scaling dimension d = 1/8 is controlled by different operator that also
requires different boundary conditions (see for instance Appendix D). According to our data
presented in Fig. 7 (a) the Luttinger parameter reaches its critical value K c = 1 at g ≈ 0.29,
in good agreement with Ref. [17]. At the same point we observe that the central charge
jumps from c ≈ 1/2 typical for the Ising transition to c ≈ 3/2 as shown in Fig. 7 (c). Aroundthe same point, the incommensurability appears, therefore the multicritical point M is also a
Lifshitz point. In Fig.fig:h1(d) we show that the lowest excitation energies scale to zero as
1/N z with the dynamical critical exponent of the Lifshitz point z = 3. All these results agreewith the previous study by Rahmani et al. [17].
12SciPost Phys. 14, 152 (2023)

1

(a)

1.2
1

0.8

0.8

0.6

0.6

0.4

0.4

(b)

0.2

0.2
0

0.5

1

1.5

2

0.8

(c)

1.6

0
0

1.4

0.5

1

1.5

2

(d)

0.6

1.2
0.4

1
0.8

0.2

0.6
0

1

2

3

4

0

0

0.4

0.8

1.2

10 -3

Figure 7: (a) Luttinger liquid parameter K and (b) incommensurate wave-vector q extracted by fitting the Friedel oscillations along h = 1 line with open and fixed boundary conditions. (c) Central charge as a function of coupling constant g extracted
from the entanglement entropy for chains with periodic boundary conditions. (d)Finite-size scaling of the three lowest excitation energies at the multicritical point M
for systems with up to 19 sites with open and fixed boundary conditions. The scaling
is in excelent agreement with the theory prediction ∆ ∝ N −z for the Lifshitz pointwith z = 3.
When does the Luttinger liquid terminate ? For large values of g the operator that becomesrelevant is the one that breaks translation symmetry and lead to the period-2 phases. This
happens for K < 1/4 [19]. By contrast to the low-g case, however, the translation symmetryis broken in both surrounding gapped phases and at the critical line at h = 1 between them.
So the critical value of the Luttinger liquid exponent is always equal to K c = 1/4. According
to our data presented in Fig. 7 (a), along h = 1 the Luttinger liquid parameter reaches thecritical value K c = 1/4 around g ≈ 1.3.
Note that the central charge extracted from entanglement entropy in periodic system starts
to deviate from c = 3/2 for large g. As one can see in Fig. 7 (c) this deviation is very slow
(as expected for a Kosterlitz-Thouless transition) but it is also very systematic. According toour finite-size data for the entanglement entropy (restricted to small lengths due to periodic
boundary conditions) the deviation starts around g ≈ 1.8. This overestimates the end of
the Luttinger liquid compared to the more reliable numerical estimates of Luttinger parameter. In Ref. [17] such deviation starts in the interval 0.5 < g −1 < 0.75 (this corresponds to
1.33 < g < 2). In both cases the deviation from the c = 3/2 starts much earlier than the
expected end point at g ≈ 3 (g ≈ 2.86 in Ref. [17]) where the transition turns into 1st order.

13

SciPost Phys. 14, 152 (2023)

1
0.5
0
-0.5
-1
0

100

200

300

400

500

600

700

800

Figure 8: Local magnetization along x direction inside the period-2-Z2 phase. The
amplitude of the oscillations A(σix ) is extracted as a difference between the largest
and smallest values that 〈σix 〉 takes over an interval of length N /4 in the middle ofthe chain. Presented results are for N = 800 site with boundaries polarized in the x
direction.

4.2

Continuation of the Ising transition and its end point

What happens along h = 1 line after the Luttinger liquid terminates ? For larger values of g
we expect the transition to continue in the Ising universality class in agreement with spontaneously broken Z2 symmetry between the two gapped phases. In order to confirm thisprediction numerically we extract the critical exponent β by fitting the order parameter upon
approaching the transition. As discussed above in order to adjust the Ising order parameter tothe incommensurate case we consider an amplitude of the x-component of the magnetizationA(σix ) over the finite window in the middle of the chain. An example for the period-2-Z2 phaseis shown in Fig. 8. We keep track of the amplitude A(σix ) when approaching the transition fromthe period-2-Z2 phase, by fitting the data we extract the corresponding critical exponent β. In
Fig. 9 we show that along g = 2 and g = 2.5 the amplitude vanishes with the critical exponent
that is in excellent agreement with the prediction β = 1/8 for Ising critical line. According toour results the end point of the Ising critical line is located around g ≈ 3 (which is within theerrorbar of the end point predicted in Ref. [17]). At this point the extracted critical exponent
β ≈ 0.0413 agrees within 1% with the CFT prediction β = 1/24 for the tri-critical Ising point.
Beyond this end point the amplitude remains finite (see Fig. 9) indicating the first ordertransition between the phases. The possibility of the first order transition in an extended
version of the phase diagram has been already discussed in Rahmani et al. [17] but it wassuggested that the ground-state is four-fold degenerate along this line. We think that thedegeneracy is actually higher - six-fold - this corresponds to the level crossings of the two-folddegenerate state from the period-2 phase and four-fold degenerate one from the period-2-Z2phase. Six-fold degeneracy, or in other words three-fold degeneracy in each sector of broken
translation would also agree with the tri-critical end point detected numerically.

4.3

The persistency of incommensurability

So far we have shown that the Luttinger liquid phase terminates at the point which is different
from the end point of the Ising critical line. Let us now show that at none of these two points
the incommensurability terminates. In Fig. 10 we show the scaling of the spin-spin correlation
with distance, where on top of the main decay one can clearly distinguish the helices typicalfor incommensurate correlations. Note that on N = 801 the presence of incommensurability
is still clearly visible for g = 3.5. We expect that the correlations slowly approach but probably
never reach the commensurate value of q for any finite g.
14

SciPost Phys. 14, 152 (2023)

1.6
1.5

1.4

1.4

1.2

1.3

1

1.2

0.8
1.1

0.6
1

0.4

0.9

0.2
0

(b)

(a)
1

1.05

1.1

1.15

1.2

1.25

1.3

1.35

1.4

10

-2

10

-1

Figure 9: Scaling of the amplitude A(σix ) upon approaching the transition between
the period-2 and the period-2-Z2 phase. Symbols are DMRG data ; red, magenta andgreen lines are fits of the results for N = 800 to the power-law A(σix ) ∝ (h − 1)β ;blue and cyan lines are guide to the eyes. Along g = 2 and g = 2.5 the extracted
exponent β is in excellent agreement with the theory prediction for Ising criticalexponent β = 1/8, while for g = 3 the critical exponent agrees with the tri-critical
Ising end point with β = 1/24. For g = 4 and g = 6 the amplitude remains finite
up to h = 1 in agreement with a first order transition. Inset shows the same data in
log-log scale.

5

Majorana zero modes

5.1

Definition and properties of strong Majorana zero modes

The spontaneous breaking of Z2 symmetry, commonly signaled by a local (magnetic) order inthe thermodynamic limit, may take a non-local form in terms of fermions, as first observedby Kitaev [2] for non-interacting chains. Indeed, in the ordered phase of the TFI chain (seeEq. (1) with g x = gz = 0 and J > h), fermions display non-trivial topological properties, with
“unpaired” zero-energy Majorana edge states localized at the boundaries of an open chain.
This brought the notion of so-called strong zero-mode (SZM), popularized by Fendley [54].
This terminology “strong” refers to the fact that the entire many-body spectrum is concerned,in contrast with “weak” zero modes which only touch the low-energy part [55, 56]. We recall
that a SZM operator Ψzm must have the three following properties:†
1) Normalizable: Ψzm
Ψzm = 1.
2) Commuting with the Hamiltonian (at least in the thermodynamic limit) [H, Ψzm ] → 0.
Qz3) Anti-commuting with the discrete symmetry: {P, Ψzm } = 0, where P =i σi is theparity operator.
A striking consequence of these properties is that Ψzm provides a map between even andodd sectors for the entire many-body spectrum of H which therefore displays a pairing degeneracy that becomes exact in the thermodynamic limit. For a finite system, this pairing
degeneracy is usually lifted and a small residual finite-size parity gap exists between even andodd sectorsN∆p (N ) = Eodd − Eeven ∼ exp −
,
(8)
ξzm
exponentially vanishing with a correlation length ξzm which also controls the edge state localization [2, 54]. Below we discuss Majorana zero modes physics, first for non-interacting
systems where SZM operators can be exactly constructed. This will help us to understand15

SciPost Phys. 14, 152 (2023)

10 -1

10 -2

10 -3

10 -4

0

50

100

150

200

250

300

350

400

Figure 10: Correlations as a function of distance for h = 1 and and various values ofg. The incommensurability persists even beyond the tri-critical end point at g ≈ 3.
the behavior of the parity gap in the Z2 ground-state of the interacting problem at small g, a
regime well captured by a self-consistent mean-field decoupling.

5.2

Non-interacting Kitaev chains

5.2.1

Transverse-field Ising chain

There are only a few cases where exact expressions for SZM operators Ψzm are known, thesimplest and most famous example being the TFI chain model, defined for an open wire of N
sites by
!
N
−1
N
N
N
−1
X
X
X
X
σzj = −i Jb j a j+1 − h
aj bj .
(9)
H=J
σ xj σ xj+1 − h
j=1

j=1

j=1

j=1

With open boundary conditions, one can easily express the SZM operators, localized at the
edges
N
N
X
X
Left
Right
Ψzm =
Θ j a j and Ψzm =Θ j bN +1− j .
(10)
j=1

j=1

The amplitude Θ j decays exponentially away from each open endjΘ j ∝ exp −
,
ξzm
5.2.2

with

ξzm =

1
.
ln(J/h)(11)

Incommensurability and parity switches for Kitaev wires

One can also derive exact expressions in the more general case provided by the non-interacting
lattice model of p-wave superconducting wires [2], when hopping and pairing terms are notnecessarily equal t ̸= ∆. This model can be written using the three equivalent languages

16

SciPost Phys. 14, 152 (2023)(spins, Majorana and Dirac fermions: see Appendix A.2), as follows
X 




H=
t c †j c j+1 + h.c. + ∆ c †j c †j+1 + h.c. + 2hn j
j

=

X

y

y

X σ xj σ xj+1 + Y σ j σ j+1 − hσzj



j

= −i

X


X b j a j+1 − Ya j b j+1 − ha j b j ,

(12)

j

where t = X + Y and ∆ = X − Y . The condition of existence for SZM is quite simple (seeAppendix B) as it simply requires X + Y > h. However, as first discussed for spin correlationsin their seminal work by Barouch and McCoy [57], there is an interesting oscillatory regime
for the pair-wise correlations where precisely the SZM also display incommensurate (IC) modulation [58, 59]. This occurs under the simple condition h2 < 4X Y , with an IC wave-vector q
given byh.
(13)
cos q = p
2 XY
The amplitude Θ j of the left and right zero-modes operators Eq. (10) then obeysΘ j ∝ sin (q j) exp −

j
ξIC
zm


,

with

v
tX
1
.
=
ln
ξIC
Y
zm

(14)

In such a case, one can show that the finite-size parity gap also displays some IC modulation. It
is straightforward to obtain an exact analytical expression for the parity gap (see Appendix B.2
for details), which in the incommensurate regime readsNs 2 sin [q(N + 1)](IC)∆p (N ) = 2X M x
exp − IC ,
(15)
sin q
ξzm
where Msx is the surface magnetization (see Appendix B.2). Contrary to the TFI model, here
the finite-size parity gap may vanish exactly if sin [q(N + 1)] = 0 and q ̸= 0. On a finite chainof length N , this occurs N /2 times, i.e. when
 nπ 
h∗n
= cos
,
p
N +1
2 XY

n = 1,...,

N
,
2

(16)

corresponding to a level crossing within each doublet of opposite parity p = ±1, associatedwith a parity switch [59]. The commensurate-incommensurate
crossover inside the Z2 phasepoccurs for n = 0, i.e. when hIC = 2 X Y . On the other side (commensurate regime for
X + Y > h > hIC ) one recovers a monotonous exponential decay of the form Eq. (8) with


1
X
= ln
.
(17)
p
ξCzm
h + h2 − 4X Y
The SZM localization length ξzm and the incommensurate vector q, given by Eqs. (13), (14)and (17), are both shown in Fig. 11 (a) for the non-interacting Kitaev chain model Eq. (12)with hopping t = 1 and pairing ∆ = 0.8 (X = 0.9, Y = 0.1) as a function of the field h. The
surface magnetization Msx (which plays the role of the order parameter of the Z2 phase, see
Appendix B.2) is also plotted. In panel (b) we show for the same set of parameters how theparity gap Eodd − Eeven behaves as a function of h. For various system lengths, ∆p oscillates
in the incommensurate regime (with N /2 oscillations), as predicted.
17

SciPost Phys. 14, 152 (2023)
1
(a)

1
0.8

0.8

0.6

1/ξ
1/xi
zm
q/π
PHI2
s
Mx
(M
x)

1

0.8

0.6 0.6

0.4

1/xi
PHI
Mx

0.4 0.4

0.2
0

0 0

ℤ2

1
(b)
1

1

(Topological)

0.001

1e-06

L=12
Commensurate
L=16
L=20
L=24

Δp

1e-06

1e-06

1e-09

1e-09

1e-09

1e-12

1e-12

1e-12

0

0.2
0

hc

Incommensurate

0.001

0.001

PM
(Trivial)hIC

0.2 0.2

1/xi
PHI
Mx

0

0.2

0.4

0.4
0.2

0.6
0.4

0.6

h 0.8

0.6

0.8

L=12
L=16N L=12
= 12
L=20N = 16
L=16
L=24
N L=20= 20
N L=24
= 24
1

1
0.8

1

Figure 11: Free-fermion results for the non-interacting Kitaev chain model Eq. (12)with open boundary conditions, hopping t = 1 and pairing ∆ = 0.8 (X = 0.9,
Y = 0.1). (a) Field h-dependence of the SZM localization length ξzm and the incommensurate vector q, given by Eqs. (13), (14) and (17), plotted together with thesurface magnetization Msx (see Appendix B.2). (b) The finite-size parity gap ∆p (symbols: exact diagonalization results ; lines: analytical expressions, see Appendix
B.3)p
displays incommensurate oscillations, as predicted by Eq. (15) for h < 2 X Y = 0.6,
followed by a commensurate Z2 topological phase, and then the PM (trivial) regime
for h > 1.

5.3

Interacting case: DMRG results

In contrast with free-fermion systems, the situation is much more complicated in the interacting case [13, 14]. Despite recent efforts addressing the very existence of possible SZM
operators in the presence of interaction [60–63], no explicit and exact construction of SZM isknown, except the notable exemple of the integrable XYZ chain [64]. For the interacting andnon-integrable model Eq. (1) with g x ̸= gz ̸= 0, Kemp et al. [38] have shown the existenceof an “almost” SZM that almost commutes with the Hamiltonian. Nevertheless, evidences foredge modes have been provided for various interacting models, but only at low energy. Interestingly, incommensurability may help this detection with the possibility to continuously
tune the effective coupling (and therefore the parity gap) between the edges [7, 65, 66], andeventually make it vanishing. However, the strong character of the SZMs turns out to be much
more difficult to grasp because it addresses the whole many-body spectrum, as recently discussed in the context of many-body localization at infinite temperature [31, 67]. Here we will
focus on ground-state physics, and therefore will not address this issue.
5.3.1

Z2 antiferromagnet at small g

We first numerically investigate the effect of weak interaction strength g on the Majorana zeromodes of the TFI chain model in the Z2 regime h < 1 using DMRG simulations, addressing
the spectral pairing of the two lowest energy levels of the many-body spectrum. We detect thepresence of zero modes by computing finite-size energy difference between the two groundstates with open boundary conditions. This is done by targeting multiple states of the effective
Hamiltonian and keeping track of the energy as a function of DMRG iterations. The method
18

SciPost Phys. 14, 152 (2023)

10 -8

0.2

(a)

3
2

(b)

0.1

1
0

0
-1

-0.1

-2
-3
0

0.2

0.4

0.6

-0.2

0.8

2

3

4

5

6

Figure 12: DMRG results for the ground-state zero modes of Z2 phases of the interacting Majorana chain model Eq. (4). E0 states for the ground-state energy andE1 is the energy of the in-gap excitation. (a) Weak interaction Z2 phase: relative
energy of the in-gap states as a function of field h along an horizontal cut at g = 0.2
computed for N = 30 sites. Incommensurate correlations tuned by h lead to exact
level crossings between the two in-gap states with even and odd parity. (b) Stronglyinteracting regime Z2 -period-2 phase: same as (a) but as a function of g for h = 1.5.
Exact crossings (also induced by incommensurability) are clearly visible, together
with an overall decay with increasing system sizes N = 30, 40, 80.
is described in details in Ref. [49].
In Fig. 12 (a) we present results for the parity gap ∆p along an horizontal cut in the phase
diagram (Fig. 2) at g = 0.2 as a function of h. Results, shown for N = 30 sites, are displayedwith respect to the average energy of the two lowest levels, such that exact level crossingsare clearly visible. Note that similar oscillations (and level crossings) have been observed inRef. [39] with small chains for the interacting spin model Eq. (1) with g x = 0 and gz > 0. Thisis a nice example of field-induced parity switches, as already discussed above for the IC regimeof the non-interacting Kitaev chain, see Fig. 11 (b) and Eq. (15). In this exactly solvable case,the parity switch condition Eq. (16) yielded N /2 level crossings. In the present interacting case,
our DMRG results strongly suggest a similar condition with N /2 parity switches, as observed
numerically in Fig. 12 (a). This conjecture is also supported by the self-consistent MF approachat weak interaction, as we will explain below in Section 5.4.
This exact level crossing regime with parity switches is dictated by the incommensurability
vector q on top of the AF Z2 order. Upon increasing further h, this oscillatory regime disappears and is replaced by the more conventional commensurate Z2 topological order with amonotonous parity gap, but still exponentially vanishing with N as in Eq. (8) ; this occurs forh > hIC with hIC ≈ 0.87 for g = 0.2. Below in Fig. 14 we also discuss similar effects for othervalues of the parameters.
5.3.2 Z2 period-2 antiferromagnet at large g
Majorana zero-modes are also expected for the upper right part of the phase diagram (seeFig. 2) in the strong interaction regime for h > 1. Indeed, we expect two discrete symmetries
to be spontaneously broken there: the translation and the Z2 . A sketchy representation of the
4 associated ground-states would be for the two sectors of translation| →→←←→→← ···〉 ± | ←←→→←←→ ···〉

(18)

| ←→→←←→→ ···〉 ± | →←←→→←← ···〉.

(19)

Although the translation is explicitly broken in our DMRG simulations on open chains, one
can still probe the breaking of Z2 by monitoring the parity gap. Panel (b) of Fig. 12 provides
19

SciPost Phys. 14, 152 (2023)X̃
a1 h b1 X a2

b2

a3

b3

aL

bL

Y
Figure 13: Schematic picture for the MF Majorana Hamiltonian where the interaction
terms now correspond to 2-body Majorana couplings at distance 3. The MF model
becomes self-dual when Y = Xe , which precisely occurs if X = h.
DMRG results in this strong interaction regime for h = 1.5. One sees exact crossings inducedby the incommensurability. Due to proximity of the wave-vector to its commensurate value(at the Kosterlitz-Thouless transition the wave-vector is q ≳ 0.95π) we can detect only one
crossing for N = 30 and N = 40 and two crossings for N = 80. As expected the location of
the crossing changes with the system size. Of course, one might expect more crossings forlarger systems but the accuracy of our simulations are not sufficient to resolve in-gap states inthis part of the phase diagram for longer chains. Qualitatively our results are consistent withexponential decay of the amplitude of the energy difference with the system-size N .

5.4

Self-consistent mean-field theory

We now present a self-consistent mean-field (MF) treatment of the interacting terms which
gives surprisingly good results. This approach allows us to gain some physical insight and abetter understanding of the interacting Z2 topological regime, at least for weak interactions.
5.4.1

Mean-field Hamiltonian

We perform a MF decoupling of the g interaction terms (see Appendix C) in the weakly interacting limit. This leads to the following effective Hamiltonian in terms of Dirac and Majorana
fermions (up to irrelevant constant terms)
X 





HMF =
t j c †j c j+1 + h.c. + ∆ j c †j c †j+1 + h.c. + 2h j n j + Xe j c †j c j+2 + c †j c †j+2 + h.c.
j

= −i

X


X j b j a j+1 − Y j a j b j+1 − h j a j b j + Xe j b j a j+2 ,

(20)

j

where t j = X j + Y j and ∆ j = X j − Y j . Rewriting the same MF Hamiltonian in the spin language
is also useful:
X

y y
HMF =
X j σ xj σ xj+1 + Y j σ j σ j+1 − h j σzj + Xe j σ xj σzj+1 σ xj+2 .
(21)
jInterestingly, the non-interacting MF Hamiltonian is very close to a non-interacting Kitaev
chain model Eq. (12), with an additional second neighbor hopping term −iXe j b j a j+2 whichtakes the form of an α-chain model [68] with α = 2.4 A sketch of the MF model is given in
Fig. 13.
5.4.2

Numerical results¦©
The new couplings X j ; Y j ; h j ; Xe j of the MF Hamiltonian Eqs. (20) and (21) are obtained
from the self-consistent equations (see Eqs. (C.4 ) in Appendix C where some technical details
about the microscopic and convergence aspects are also discussed, OBC are used for the MF
4

Also equivalent to the so-called cluster model in spin language [69, 70].

20

SciPost Phys. 14, 152 (2023)
(a)

-6

10

h = 0.5
N=20
N
= 20
N
= 30
N=30
L=20
N
= 40
N=40
SCMF
MF

-2

-8

10

-10

10

10

10

-8

0

-10

10

g = 0.2

-8

30
20

-6
-10

10

-6

10

No MF

-8
-12

10

N = 302030
N = 20
30
20

10

-6

-8

10 10

-12

-10

MF
L=20

-1010

10

-14

0

0.2

1

SCMF
MF

10

10

10

= 20
L=20
0.8 N

-12

10
-8

10

0.6

0.4

-10

10

10

0.2

10

-6

-12

g = 0.1

10

Δp

10

-4

(b)

-6

10

0.4 10 -0.1

0.6

0

0.8

0.1

g1

0.2

0.3

0.4

-10

0

10

-8

0.2

0.4

h

0.6

0.8

1

10

Figure 14: Quantitative comparison between DMRG (symbols) and self-consistentmean-field (SCMF) (lines) for the parity gap ∆p . (a) Vertical scan in the phase diagram of Fig. 2 at h = 0.5: data are shown for N = 20, 30, 40. The agreement for g
not too large is quite remarkable, the green shaded area signals the region where no
MF convergence can be obtained. (b) Horizontal scan at g = 0.1 (upper panel) and
g = 0.2 (bottom panel): the agreement, quantitatively excellent for g = 0.1, slightly
deteriorates as g increases, but the periodicity of the oscillations remains perfectly
captured by the SCMF approach.
-6parameters). One can easily reach a converged MF solution at weak interaction g, which is
the target regime of such Hartree-Fock decouplings, but we also observe a breakdown of the
MF convergence when g exceeds a certain value ∼ 0.25, see below in Section 5.4.3 where theMF phase diagram is shown.
diagonalized iteratively till the convergence of the MF parameters¦ HMF is numerically©
X j ; Y j ; h j ; Xe j . In Fig. 14, we illustrate the exact diagonalization results of the converged
Hamiltonian HMF for the parity gap ∆p . A direct comparison with DMRG calculations is provided for the evolution of ∆p (a) vs. g for h = 0.5, and (b) against h for g = 0.1 and 0.2.
The agreement is quantitatively excellent for weak interaction strength, typically below ∼ 0.2,
with both the amplitudes and the incommensurate oscillations very well captured despite theapproximation. Interestingly the parity gap displays a similar behavior as compared to thestandard non-interacting Kitaev chain (Fig. 11), which undoubtedly captures the essential
features of this Z2 regime of the phase diagram.
5.4.3

Mean-field phase diagram

One can further gain insights about the qualitative effects produced by the interactions on
the non-interacting Z2 line of the TFI model, by looking more closely to the MF decoupling in
Eqs. (C.4 ). Indeed, keeping only the dominant corrections, we approximate the new couplings
by




xx
xx
X j ≈ J − g |C j−1| + |C j+1| , Y j ≈ g|C jx x | , h j ≈ h − g mzj−1 + mzj+1 , Xe j ≈ g mzj+1 , (22)where mzj = 〈σzj 〉 and C jx x = 〈σ xj σ xj+1 〉. These expressions can be even more simplified in
the limit of small field h, where mzj is small and C jx x ≈ −1 (AF order). Assuming uniform¦©
¦©
couplings, the toy-model with X j ; Y j ; h j ; Xe j
= J − 2g ; g ; h ; 0 , exactly solvable, is∀jexpected to capture some features of the interacting Majorana chain model in the Z2 regime
at small enough interaction g and field h.
In Fig. 15 (a) we verify that the disorder line separating commensurate and incommensurate regions of the Z2 phase, trivially given by h2IC = 4g for the toy-model, reproduces very well
21

SciPost Phys. 14, 152 (2023)
(a)

0.2
0.2

g
0.1
0.1

0.2

0.4

0.8

0.6
00

MF breakdownMF
breakdown
NN=8
=8
NN=16
= 16
NN=24
= 24
NN=32
= 32
DMRG
DMRG
h 2 /4

00

0.2
0.2

(b)

0.3

No MF

0.3
0.3

MF breakdown
N=8
N=16
N=24
N=32
DMRG

0.2

0.4
0.4

h

0.6
0.6

0.8
0.8

10−4

10−8

0.1

ℤ2 − C

1No MFℤ2 − ICgℤ2 − ICΔp10−12ℤ2 − C
0

11

0

0.2

0.4

h

0.6

0.8

1

10−16

Figure 15: SCMF phase diagram in the weak-interaction regime displayed in theplane field h – interaction g. (a) The C-IC disorder line, computed using SCMF forvarious chain lengths (as indicated on the plot), is found in excellent agreement with
DMRG data (the magenta line is a guide to the eyes).The full line g = h2 /4 is the
toy-model estimate for the C-IC line (see text) which gives an excellent description.
Green diamonds show the boundary of the SCMF convergence: in the top green area
SCMF does not converge anymore. (b) Color map of the parity gap ∆ p , numericallyobtained for1
N = 16. The parity switches lines are clearly visible where the parity
gap vanishes.

8

the more involved SCMF results. Quite remarkably, this remains true for the entire Z2 regime,
well beyond the h ≪ 1 limit. A direct comparison with DMRG results for the commensurateincommensurate line is also is excellent agreement.
The right panel of Fig. 15 (b) provides the MF phase diagram of HMF as a color map of theparity gap in the plane field h – interaction g. In addition to the disorder line, one sees N /2
dark lines showing the minima of ∆p which signal N /2 parity switches. For readability andclarity, this is shown for a small system size N = 16.

6

Discussions and conclusions

In this work we have presented an extended phase diagram of the interacting Majorana chain
in the case of symmetric interaction gz = g x , as overviewed in Fig. 2. The obtained phase
diagram is very rich and contains four gapped phases and two floating phases, i.e. Luttinger
liquid regimes with incommensurate correlations. In addition, our results motivated us torevisit the nature of the effective critical theories along the self-dual critical line at h = 1. Weargue that there is no generalized commensurate-incommensurate transition where Ising andLuttinger liquid criticalities both terminate. We demonstrate that the Luttinger liquid phasestops at g ≈ 1.3, while the Ising critical line persists beyond it up to g ≈ 3 where it ends witha tri-critical Ising point. Incommensurability persists beyond these two points and, if it ends,
it does so for a much stronger g interaction.
We have also provided numerical evidence of a multi-critical point in the eight-vertex universality class at h = 0 and g ≈ 0.41. Note that due to duality one can expect a second criticalpoint in the eight-vertex universality class at h → ∞ (or J → 0) and g/h ≈ 0.41 between
the paramagnetic and the period-2-Z2 phases. This implies that both, translation and parity,
symmetries will be broken on one side of the transition and, quite surprisingly this does notchange the nature of the critical point.
Interestingly, the floating phases in the present phase diagram are separated from thephases with broken translation symmetry by Kosterlitz-Thouless transitions, by contrast to
previous studies where a Pokrovsky-Talapov transition has been reported [19, 41]. This is be-

22

SciPost Phys. 14, 152 (2023)

cause the short-range order in the gapped phases is characterized by incommensurate wavevectors, and thus the transition between floating and period-2 phases are incommensurateincommensurate (and not commensurate-incommensurate as in the case of Pokrovsky-Talapov
transition). This short-range incommensurability is realized when gz = g x and absent wheng x = 0. Therefore by tuning g x /gz one might expect a crossover between the PokrovskyTalapov and the Kosterlitz-Thouless transitions similar to the one reported recently in a chainof spinless fermions with next-nearest-neighbor interactions [71]. It would be interesting tocheck this prediction numerically.
The tri-critical Ising conformal field theory is supersymmetric, thus we might expect thatthe end point of the Ising critical line will have the features of the supersymetric critical point,
for example, doubly degenerate excitation spectra. But the tri-critical Ising point is not the
only place in this phase diagram where one might expect supersymmetry to appear. The Ising
critical line that lies within the gapless region with U(1) and Z2 symmetry might have an
emergent N = (1, 1) supersymmetry [72, 73]. In the present case, the Z2 is satisfied by the
Hamiltonian, while the U(1) symmetry is an emergent symmetry stabilizing the Luttinger liquid phase [17,19]. Moreover, at the multicritical point where Ising transition exit the Luttinger
liquid phase and hits the Kosterlitz-Thouless transition one might expect the spontaneously
emergent N = (3, 3) supersymmetry. According to Ref. [74] this higher supersymmetry canbe realized under the following three conditions: i) the system has to be invariant with respect to U(1) and Z2 symmetry ; ii) it has to be tuned to the multicritical point where Isingand Kosterlitz-Thouless transitions coincide ; and iii) the velocity of the fermionic degree of
freedom should be smaller than or equal to the velocity of the bosonic degree of freedom.
The first two conditions are formally satisfied at the both multicritical points terminating thec = 3/2 line (see points M and S at the phase diagram of Fig. 2). Although these two points
appears very differently: point S appears as an intersection between the Ising the KosterlitzThouless lines; while due to continuity of the equal-K lines the Kosterlitz-Thouless transitionis expected to approach the point M at an infinite slope. The last condition on the velocities
requires further investigation. We leave this problem for future studies.
Finally, we have explored the topological properties of the systems by looking at the possible existence of Majorana zero-modes, signalled by the (two-fold) parity degeneracy. The Z2
topological order of the (Kitaev) free-fermion (for h < 1) remains robust against weak interactions, as unambiguously shown by large-scale DMRG simulations, and remarkably captured
by a self-consistent fermionic mean-field theory. Due to incommensurability, one observes a
succession of exact zero-mode level crossings already for finite chain lengths. A similar observation is also reported at large g in the Z2 period-2 regime. Despite such good evidences forthe existence of Majorana edge modes at low energy, their existence at higher energies, as itis the case in the absence of interaction, remains highly hypothetical.

Acknowledgements
NC acknowledges useful discussions with Naoki Kawashima and Kareljan Schoutens. NL isgrateful to Pasquale Marra for stimulating discussions. This work has been supported by
Delft Technology Fellowship (NC) and by the French National Research Agency (NL), project
GLADYS ANR-19-CE30-0013, and (NL) the EUR grant NanoX No. ANR-17-EURE-0009 in theframework of the “Programme des Investissements d’Avenir”. Numerical simulations havebeen performed at the Dutch national e-infrastructure with the support of the SURF Cooperative.

23

SciPost Phys. 14, 152 (2023)

A

Useful transformations

A.1

Duality

The model defined by the Hamiltonian Eq. (1) up to boundary terms transforms into itself by
Kramers-Wannier duality transformation:xx
σix σi+1→ τzi and σzi → τix τi+1
,

(A.1)

where σ and τ are Pauli matrices. Then the duality transformation is given by:
X
X
X

x
x
x
H=
σix σi+1+h
σzi + gσix (σi+1)2 σi+2+ σzi σzi+1
i

=h

i


X1
i

h

τzi +

X
i

i


X
g
x
x
τix τi+1+(τzi τzi+1 + τix τi+2
) ,
h i

(A.2)

x
where we used that (σi+1)2 = I and the fact that the model is symmetric with respect to thesign of the field h. Up to a pre-factor the Hamiltonian of Eq. (A.2 ) is equivalent to the original
Hamiltonian in Eq. (1) with h → h−1 and g → g/h. This duality transformation allows us to
study the nature of the quantum phase transitions only for h ≤ 1, and re-use these results forh > 1.

A.2

Jordan-Wigner and Majorana fermions

The Jordan-Wigner transformation maps Pauli operators to Dirac fermionsσzj = 1 − 2c †j c j ,
σ xj = K j c †j + c j ,yσ j = iK j c †j − c j ,
with

Kj =

j−1
Y

σzk .

(A.3)
(A.4)
(A.5)
(A.6)

k=1

Therefore the original spin model Eq. (1) can be rewritten as an interacting fermions problem,see Eq. (2). One can also introduce two Majorana fermion operators a j and b j
a j = c †j + c j = K j σ xj ,yb j = i(c †j − c j ) = K j σ j ,a j b j = i(1 − 2c †j c j ) = iσzj ,which satisfy the usual rules: (a/b)† = (a/b), (a/b)2 = 1, {ai , a j } = {bi , b j } = 2δi j , and{ai , b j } = 0. Then the interacting problem can be re-written in terms of Majorana variables asgiven by the Hamiltonian Eq. (3).

B

Strong zero modes for the non-interacting Kitaev model

B.1

SZM construction: iterative procedure

In the Majorana representation, the Kitaev model readsX
H = −iX b j a j+1 − Ya j b j+1 − ha j b j .
j

24

(B.1)

SciPost Phys. 14, 152 (2023)

Assuming simple linear combinations for the (normalized) SZM operators
N
1 X
Ψleft =
Θ a
NN j=1 j jN
1 Xand Ψright =
Θ b
,
(B.2)
NN j=1 j N +1− j






and using H, a j = 2i X b j−1 − hb j + Y b j+1 and H, b j = 2i −Ya j−1 + ha j − X a j+1 , weiteratively arrive at the simple recursion relation for Θ j :Θ j+1 =Yh
Θ j − Θ j−1 ,
X
X

(B.3)

such that[H, Ψleft ] = 2i

1
(Y ΘN −1 − hΘN ) bN
NN

and




1
H, Ψright = 2i
(Y ΘN −1 − hΘN ) a1 . (B.4)
NN

It is easy to solve the above recursion Eq. (B.3 ) with initial conditions Θ0 = 0 and Θ1 = 1. We
restrict the present discussion to positive couplings h, X , Y ≥ 0 and X ≥ Y (it is straightforward
to get results for the other cases).
The phase diagram of the Kitaev chain model can simply be infered from the existence of
normalizable SZM, which requires h < X + Y . Contrary to the TFIM case, here the topologicalregime is richer as one can distinguish two types of SZM decays.
(i) Commensurate regime if h2 ≥ 4XY.
¨



∞,
if h > X + Y ,
X
h j 
j
j
p
Θj =
× (1 + α) − (1 − α) −−→ X − j/ξC
zm ,
j≫1
αh 2X
if X + Y > h > 2 X Y ,αh epwhere α = 1 − 4X Y /h2 , and the edge mode localization length given by


1
2X
= ln
.
(B.5)
ξCzm
(1 + α)hpNote that in the limit α → 0 (h → 2 X Y from above)


h j−1
Θ j −−→ j
.
(B.6)
α→0
2X
The SZM normalization factor can be expressed in the large N limit:11−−−−→ 2α 
N N →∞(1 + α)−2 −

 +
h 2
2X

1
(1 − α)−2 −

 −
h 2
2X

−1/2

2
(1 − α2 )−1 −


h 2
2X



.

(B.7)

(ii) Incommensurate (IC) regime if h2 < 4XY.
IC
2X
Θj = psin (q j) e− j/ξzm
4X Y − h2(B.8)

displays oscillations and exponential decay, controlled byhcos q = p
2 XY

and

v
tX
1
=
ln
.
ξIC
Y
zm

The normalization factor is given by

−1/2
Y
p
1
1
X − cos (2q)−−−−→ 2 sin q 
+
.
 
Y
Y
Y 2
N
→∞
N
1− X
1 − 2 cos (2q) +
X

25

X

(B.9)

(B.10)

SciPost Phys. 14, 152 (2023)

B.2

Surface magnetization

An hallmark of the topological phase is the presence of edge states. The surface magnetization,a key-quantity in that respect, is defined by [75]xMsx = 〈 E − |σ1,L| E+ 〉 ,(B.11)

where | E ± 〉 are many-body eigenstates of parity p = ± associated to a partner | E ∓ 〉 of opposite
parity −p, obtained by acting the SZM operator, i.e. Ψleft/right | E ± 〉 ≈ | E ∓ 〉. From the definitionof the SZM operators Eq. (B.2 ), we easily arrive at
1
,
N

Msx =

(B.12)

for which we have analytical expressions in both C and IC regimes, see Eqs. (B.7 ), (B.10 ).

B.3

Finite size parity gap and parity switches

The (finite-size) parity gap is obtained from the simple expression∆parity = 〈 E − |H| E − 〉 − 〈 E + |H| E + 〉 ≈ 〈 E − | [H , Ψzm ] | E + 〉

 h
X −
Y
+
−
+
≈
〈 E |a1 | E 〉 − i〈 E |bN | E 〉ΘN − ΘN −1
N
X
X



h
Y
2
≈ 2X Msx
ΘN − ΘN −1 .
X
X

(B.13)We then obtain closed forms for both regimes. In the commensurate case



1+α
h N 
× (1 + α)N − (1 − α)N
≈ 2X
2α
2X
¨

 −N /ξCs 2 1+αzm ,2X M xe
if α ̸= 0 , 2α≈
s 2 −N /ξCzm
2X M x e
×N,
if α → 0 ,(C)∆paritywhere α =

p

2
Msx



(B.14)

1 − 4X Y /h2 and ξCzm is given by Eq. (B.5 ). For the incommensurate case we get(IC)∆parity ≈ 2X Msx

N2 sin [q(N + 1)] − ICe ξzm ,
sin q

(B.15)

where ξICzm and q are given by Eq. (B.9 ).

C

Self-consistent mean-field

Mean-field decoupling and self-consistent equations. We start from the interacting Majorana chain modelXX
H = −iJ b j a j+1 − ha j b j − g
a j b j a j+1 b j+1 + b j a j+1 b j+1 a j+2 ,
(C.1)
j

j

where the quartic term that cannot be treated exactly. We make a mean-field decoupling whichleads to a noninteracting system quadratic fermionic problem, with new coupling constants
which are self-consistently determined. The quartic terms are decoupled in all mean-fieldchannels which are consistent with the Wick’s theorem [39]:26

SciPost Phys. 14, 152 (2023)

a j b j a j+1 b j+1 ≈ a j b j 〈a j+1 b j+1 〉 + 〈a j b j 〉a j+1 b j+1+ a j b j+1 〈b j a j+1 〉 + 〈a j b j+1 〉b j a j+1 ,b j a j+1 b j+1 a j+2 ≈b j a j+1 〈b j+1 a j+2 〉 + 〈b j a j+1 〉b j+1 a j+2+ b j a j+2 〈a j+1 b j+1 〉 + 〈b j a j+2 〉a j+1 b j+1 .

This leads to the following mean-field Hamiltonian (up to irrelevant constant terms)
X
HMF = −i
X j b j a j+1 − Y j a j b j+1 − h j a j b j + Xe j b j a j+2 ,

(C.2)

(C.3)

j

where the parameters are computed from the following self-consistent equations



xz x
h j = h + ig 〈a j−1 b j−1 〉 + 〈a j+1 b j+1 〉 + 〈b j−1 a j+1 〉 = h − g mzj−1 + mzj+1 + C j−1,
j, j+1 ,



yy
xx
xx
X j = J − ig 〈b j−1 a j 〉 + 〈b j+1 a j+2 〉 + 〈a j b j+1 〉 = J + g C j−1,
j + C j+1, j+2 − C j, j+1 ,Y j = ig〈b j a j+1 〉 = −g C j,x xj+1 ,Xe j = −ig〈a j+1 b j+1 〉 = g mzj+1 .

(C.4)

The local magnetization mzj = 〈σzj 〉, the correlators at distance one (involving two sites)
ββ

β

β

xz x
x
z x
C j, j+1 = 〈σ j σ j+1 〉, and at distance two (involving three sites) C j−1,
j, j+1 = 〈σ j−1 σ j σ j+1 〉 arecomputed in the ground-state of HMF . We use OBC and therefore there is a dependence on theposition along the open chain, which dies off exponentially away from the boundaries. However, we have checked that such a dependence is a finite-size effect which does not influence
the MF phase diagram. The self-consistent loop is performed numerically, and convergence isobtained when a stationary solution is reached for the coupling parameters in Eq. (C.4 ).
In Fig. 16 we show the evolution of the effective MF parameters (evaluated in the middleof the open chain for j = N /2) as a function of the interaction g for a three values of thetransverse field h = 0.1, 0.5, 0.9. Several remarks are in order here. Both the transverse
magnetic field and the nearest neighbor coupling X decrease with g, while in the same time
the newly generated terms Y and Xe both grow. This can be naturally understood from the
SCMF equations Eq. (C.4 ), which can be approximated by h j ≈ h − 2g m, X j ≈ J + 2g C ,Y j ≈ −g C , and Xe j ≈ g m, where in the limit of small field h ≪ J, the on-site magnetizationm ≈ h and the nearest-neighbor antiferromagnetic correlator C ≈ −1.

Convergence. In panel (d) of Fig. 16 we also show for number of MF iterations τMF required
to globally converge all local MF parameters (to a relative precision of 10−5 ). τMF gradually
increases with g, with a notable enhancement at not too large interaction strength, due to tiny
parity gap in the system, but MF still converge (at least provided the parity gap does not reachmachine precision). However there is a clear breakdown of MF theory at larger g where τMF
diverges (vertical arrows).

D

Ising transition for small g

In the main text we extracted Luttinger liquid exponent K by fitting the local magnetization in
the z direction that in terms of bosons/fermions corresponds to the local density. For this we
fix the boundary condition to be polarized along z. However, in order to detect Ising transitionx
breaking Z2 symmetry we have to look at the different local operator - 〈σix − σi+1〉. It is very
27

SciPost Phys. 14, 152 (2023)

hN/2 − h0
h = 0.1-0.1

h = 0.5-0.2(a)h = 0.9

XN/2

1

(b)

0.8

YN/2 and X̃ N/2

0.6

(c)

Y

0.1

X̃

0.01

τMF

1000

N = 16N = 24

100
10

0

0.1

0.1

0.2

0.3

(d)
0

0.4 0.05

0.1
0.5

0.150.6 0.2

interaction g

0.25

0.3

Mean-field parameters
after numerically0.6
solving the SCMF equa0.2 Figure 16:0.3
0.4 obtained 0.5
tions Eq. (C.4 ) for h = 0.1, 0.5, 0.9 plotted against g.

intuitive to see why this term plays a role of an order parameter for the Ising transition: forsmall g and h the ground-state corresponds to the antiferromagnetic alignment of spins along
x-direction, thus the absolute difference between the two projection is large in the Z2 phase,
while it vanishes in the paramagnetic phase. Important to notice, that in order to use this local
operator one has to fix the boundary conditions by applying the boundary field along x. In
order to keep the profile symmetric we apply the same boundary field at each edge and takeodd number of sites N = 801. The results are presented in Fig. 17. From the DMRG data (blueand green) one can see that the difference between the profiles at g = 0 and upon approaching the floating phase g = 0.25 is negligible. In both cases fit to the theory prediction of the
profile gives the estimate of the scaling dimension d ≈ 0.126 which is in excellent agreementwith the CFT prediction for the Ising model d = 1/8.

E

Comparison to the previous study at h = 1

In the main text we mentioned that although we arrived to a different conclusions regarding critical regimes along h = 1 line our numerical results to a large extent agrees with theprevious study by Rahmani et al [17]. In Fig. 18 we compare our results for the Luttinger
liquid parameter K extracted in the present paper by fitting the Friedel oscillations and thosecomputed from the finite-size energy spectra by Rahmani et al. [17]. Overall we see a very
good agreement between the results. There is slight discrepancy on where the Luttinger liquid

28

SciPost Phys. 14, 152 (2023)

1.8
1.6
1.4
1.2
1
0.8
0

100

200

300

400

500

700

600

800

Figure 17: Fridel oscillation profile at the Ising transition at h = 1 and g = 0 (green)and g = 0.25 (blue). Red line is the result of the fit with d ≈ 0.126.
1.2
1from Rahmani et al.

0.8
0.6
0.4
0.2
0

0.5

1

1.5

2

2.5

3

3.5

4

Figure 18: Comparison between the Luttinger liquid parameter K extracted in thispaper with Friedel oscillations (blue and red) versus Luttinger liquid parameter extracted frpm the plots presented in Ref. [17] and obtained with energy of the lowlying excitations.
parameter crosses 1/4 line: according to our data this happens around g ≈ 1.3 (1/g ≈ 0.75),
according to Rahmani et al’s data, this happens at 1/g ≈ 0.5 or g = 2. The main source of
this discrepancy is probably associated the finite-size effects in the energy spectra computedon systems with up to 200 Majoranas (100 spins) and used in Ref. [17] to extract Luttinger
liquid exponent K. But in any case, according to all available data K reaches its critical value
1/4 well below the end point of the continuous Ising transition located at g ≈ 3.

References[1] T. Giamarchi, Quantum physics in one dimension, Oxford University Press, Oxford, UK,
ISBN 9780198525004 (2003), doi:10.1093/acprof:oso/9780198525004.001.0001.

29

SciPost Phys. 14, 152 (2023)[2] A. Y. Kitaev, Unpaired Majorana fermions in quantum wires, Phys.-Uspekhi 44, 131(2001), doi:10.1070/1063-7869/44/10s/s29.
[3] V. Mourik, K. Zuo, S. M. Frolov, S. R. Plissard, E. P. A. M. Bakkers and L. P. Kouwenhoven,
Signatures of Majorana fermions in hybrid superconductor-semiconductor nanowire devices,
Science 336, 1003 (2012), doi:10.1126/science.1222360.
[4] L. P. Rokhinson, X. Liu and J. K. Furdyna, The fractional a.c. Josephson effect in asemiconductor-superconductor nanowire as a signature of Majorana particles, Nat. Phys.
8, 795 (2012), doi:10.1038/nphys2429.
[5] M. T. Deng, C. L. Yu, G. Y. Huang, M. Larsson, P. Caroff and H. Q. Xu, Anomalous zero-bias
conductance peak in a Nb-InSb nanowire-Nb hybrid device, Nano Lett. 12, 6414 (2012),
doi:10.1021/nl303758w.
[6] A. Das, Y. Ronen, Y. Most, Y. Oreg, M. Heiblum and H. Shtrikman, Zero-bias peaks andsplitting in an Al-InAs nanowire topological superconductor as a signature of Majorana
fermions, Nat. Phys. 8, 887 (2012), doi:10.1038/nphys2479.
[7] R. Toskovic et al., Atomic spin-chain realization of a model for quantum criticality, Nat.
Phys. 12, 656 (2016), doi:10.1038/nphys3722.
[8] C. W. J. Beenakker, Search for Majorana fermions in superconductors, Annu. Rev. Condens.
Matter Phys. 4, 113 (2013), doi:10.1146/annurev-conmatphys-030212-184337.
[9] J. Alicea, New directions in the pursuit of Majorana fermions in solid state systems, Rep.
Prog. Phys. 75, 076501 (2012), doi:10.1088/0034-4885/75/7/076501.
[10] C. Nayak, S. H. Simon, A. Stern, M. Freedman and S. Das Sarma, Non-Abelian
anyons and topological quantum computation, Rev. Mod. Phys. 80, 1083 (2008),
doi:10.1103/RevModPhys.80.1083.
[11] S. Das Sarma, M. Freedman and C. Nayak, Majorana zero modes and topological quantum
computation, npj Quantum Inf. 1, 15001 (2015), doi:10.1038/npjqi.2015.1.
[12] P. Marra, Majorana nanowires for topological quantum computation, J. Appl. Phys. 132,
231101 (2022), doi:10.1063/5.0102999.
[13] L. Fidkowski and A. Kitaev, Effects of interactions on the topological classification of freefermion systems, Phys. Rev. B 81, 134509 (2010), doi:10.1103/PhysRevB.81.134509.
[14] L. Fidkowski and A. Kitaev, Topological phases of fermions in one dimension, Phys. Rev. B
83, 075103 (2011), doi:10.1103/PhysRevB.83.075103.
[15] E. Sela, A. Altland and A. Rosch, Majorana fermions in strongly interacting helical liquids,
Phys. Rev. B 84, 085114 (2011), doi:10.1103/PhysRevB.84.085114.
[16] A. Rahmani, X. Zhu, M. Franz and I. Affleck, Emergent supersymmetry fromstrongly interacting Majorana zero modes, Phys. Rev. Lett. 115, 166401 (2015),
doi:10.1103/PhysRevLett.115.166401.
[17] A. Rahmani, X. Zhu, M. Franz and I. Affleck, Phase diagram of the interacting Majorana
chain model, Phys. Rev. B 92, 235123 (2015), doi:10.1103/PhysRevB.92.235123.
[18] H. Katsura, D. Schuricht and M. Takahashi, Exact ground states and topological order in interacting Kitaev/Majorana chains, Phys. Rev. B 92, 115137 (2015),
doi:10.1103/PhysRevB.92.115137.
30

SciPost Phys. 14, 152 (2023)[19] R. Verresen, A. Vishwanath and F. Pollmann, Stable Luttinger liquids andemergent U(1) symmetry in constrained quantum chains, (arXiv preprint)
doi:10.48550/arXiv.1903.09179.
[20] A. M. Lobos, R. M. Lutchyn and S. Das Sarma, Interplay of disorder and interaction in Majorana quantum wires, Phys. Rev. Lett. 109, 146403 (2012),
doi:10.1103/PhysRevLett.109.146403.
[21] F. Crépin, G. Zaránd and P. Simon, Nonperturbative phase diagram of interacting disordered Majorana nanowires, Phys. Rev. B 90, 121407 (2014),
doi:10.1103/PhysRevB.90.121407.
[22] A. Milsted, L. Seabra, I. C. Fulga, C. W. J. Beenakker and E. Cobanera, Statistical translation invariance protects a topological insulator from interactions, Phys. Rev. B 92, 085139(2015), doi:10.1103/PhysRevB.92.085139.
[23] N. M. Gergs, L. Fritz and D. Schuricht, Topological order in the Kitaev/Majorana
chain in the presence of disorder and interactions, Phys. Rev. B 93, 075129 (2016),
doi:10.1103/PhysRevB.93.075129.
[24] M. McGinley, J. Knolle and A. Nunnenkamp, Robustness of Majorana edge modes andtopological order: Exact results for the symmetric interacting Kitaev chain with disorder,
Phys. Rev. B 96, 241113 (2017), doi:10.1103/PhysRevB.96.241113.
[25] J. F. Karcher, M. Sonner and A. D. Mirlin, Disorder and interaction in chiralchains: Majoranas versus complex fermions, Phys. Rev. B 100, 134207 (2019),
doi:10.1103/PhysRevB.100.134207.
[26] B. Roberts and O. I. Motrunich, Infinite randomness with continuously varying critical exponents in the random X Y Z spin chain, Phys. Rev. B 104, 214208 (2021),
doi:10.1103/PhysRevB.104.214208.
[27] J. A. Kjäll, J. H. Bardarson and F. Pollmann, Many-body localization ina disordered quantum Ising chain, Phys. Rev. Lett. 113, 107204 (2014),
doi:10.1103/PhysRevLett.113.107204.
[28] R. Sahay, F. Machado, B. Ye, C. R. Laumann and N. Y. Yao, Emergent ergodicity at thetransition between many-body localized phases, Phys. Rev. Lett. 126, 100604 (2021),
doi:10.1103/PhysRevLett.126.100604.
[29] S. Moudgalya, D. A. Huse and V. Khemani, Perturbative instability towards
delocalization at phase transitions between MBL phases, (arXiv preprint)
doi:10.48550/arXiv.2008.09113.
[30] T. B. Wahl, F. Venn and B. Béri, Local integrals of motion detection of localization-protected
topological order, Phys. Rev. B 105, 144205 (2022), doi:10.1103/PhysRevB.105.144205.
[31] N. Laflorencie, G. Lemarié and N. Macé, Topological order in random interacting IsingMajorana chains stabilized by many-body localization, Phys. Rev. Res. 4, L032016 (2022),
doi:10.1103/PhysRevResearch.4.L032016.
[32] E. Lieb, T. Schultz and D. Mattis, Two soluble models of an antiferromagnetic chain, Ann.
Phys. 16, 407 (1961), doi:10.1016/0003-4916(61)90115-4.
[33] P. Pfeuty, The one-dimensional Ising model with a transverse field, Ann. Phys. 57, 79 (1970),
doi:10.1016/0003-4916(70)90270-8.
31

SciPost Phys. 14, 152 (2023)[34] S. Gangadharaiah, B. Braunecker, P. Simon and D. Loss, Majorana edge statesin interacting one-dimensional systems, Phys. Rev. Lett. 107, 036801 (2011),
doi:10.1103/PhysRevLett.107.036801.
[35] E. M. Stoudenmire, J. Alicea, O. A. Starykh and M. P. A. Fisher, Interaction effects in
topological superconducting wires supporting Majorana fermions, Phys. Rev. B 84, 014503(2011), doi:10.1103/PhysRevB.84.014503.
[36] F. Hassler and D. Schuricht, Strongly interacting Majorana modes in an array of Josephson
junctions, New J. Phys. 14, 125018 (2012), doi:10.1088/1367-2630/14/12/125018.
[37] R. Thomale, S. Rachel and P. Schmitteckert, Tunneling spectra simulation of interacting
Majorana wires, Phys. Rev. B 88, 161103 (2013), doi:10.1103/PhysRevB.88.161103.
[38] J. Kemp, N. Y. Yao, C. R. Laumann and P. Fendley, Long coherence times for edge spins, J.
Stat. Mech.: Theor. Exp. 063105 (2017), doi:10.1088/1742-5468/aa73f0.
[39] G. Vionnet, B. Kumar and F. Mila, Level crossings induced by a longitudinal
coupling in the transverse field Ising chain, Phys. Rev. B 95, 174404 (2017),
doi:10.1103/PhysRevB.95.174404.
[40] I. Mahyaeh and E. Ardonne, Study of the phase diagram of the Kitaev-Hubbard chain,
Phys. Rev. B 101, 085125 (2020), doi:10.1103/PhysRevB.101.085125.
[41] N. Chepiga and F. Mila, Eight-vertex criticality in the interacting Kitaev chain, Phys. Rev.
B 107, L081106 (2023), doi:10.1103/PhysRevB.107.L081106.
[42] S. Jalal and B. Kumar, Edge modes in a frustrated quantum Ising chain, Phys. Rev. B 90,
184416 (2014), doi:10.1103/PhysRevB.90.184416.
[43] J. Ruhman and E. Altman, Topological degeneracy and pairing in a one-dimensional gas of
spinless fermions, Phys. Rev. B 96, 085133 (2017), doi:10.1103/PhysRevB.96.085133.
[44] L. Gotta, L. Mazza, P. Simon and G. Roux, Two-fluid coexistence in a spinless fermions chain with pair hopping, Phys. Rev. Lett. 126, 206805 (2021),
doi:10.1103/PhysRevLett.126.206805.
[45] S. R. White, Density matrix formulation for quantum renormalization groups, Phys. Rev.
Lett. 69, 2863 (1992), doi:10.1103/PhysRevLett.69.2863.
[46] U. Schollwöck, The density-matrix renormalization group, Rev. Mod. Phys. 77, 259(2005), doi:10.1103/RevModPhys.77.259.
[47] S. Östlund and S. Rommer, Thermodynamic limit of density matrix renormalization, Phys.
Rev. Lett. 75, 3537 (1995), doi:10.1103/PhysRevLett.75.3537.
[48] U. Schollwöck, The density-matrix renormalization group in the age of matrix productstates, Ann. Phys. 326, 96 (2011), doi:10.1016/j.aop.2010.09.012.
[49] N. Chepiga and F. Mila, Excitation spectrum and density matrix renormalization group
iterations, Phys. Rev. B 96, 054425 (2017), doi:10.1103/PhysRevB.96.054425.
[50] R. J. Baxter, Partition function of the eight-vertex lattice model, Ann. Phys. 70, 193 (1972),
doi:10.1016/0003-4916(72)90335-1.
[51] M. den Nijs, The domain wall theory of two-dimensional commensurate-incommensurate
phase transitions, Ph. Transit. Crit. Phenom. 12, 219 (1988).
32

SciPost Phys. 14, 152 (2023)[52] N. Chepiga and F. Mila, Lifshitz point at commensurate melting of chains of Rydberg atoms,
Phys. Rev. Res. 3, 023049 (2021), doi:10.1103/PhysRevResearch.3.023049.
[53] R. J. Baxter, One-dimensional anisotropic Heisenberg chain, Ann. Phys. 70, 323 (1972),
doi:10.1016/0003-4916(72)90270-9.
[54] P. Fendley, Parafermionic edge zero modes in Zn -invariant spin chains, J. Stat. Mech.: Theor.
Exp. P11020 (2012), doi:10.1088/1742-5468/2012/11/P11020.
[55] A. Alexandradinata, N. Regnault, C. Fang, M. J. Gilbert and B. A. Bernevig, Parafermionicphases with symmetry breaking and topological order, Phys. Rev. B 94, 125103 (2016),
doi:10.1103/PhysRevB.94.125103.
[56] K. Wada, T. Sugimoto and T. Tohyama, Coexistence of strong and weak Majorana zeromodes in an anisotropic X Y spin chain with second-neighbor interactions, Phys. Rev. B
104, 075119 (2021), doi:10.1103/PhysRevB.104.075119.
[57] E. Barouch and B. M. McCoy, Statistical mechanics of the X Y model. II. Spin-correlation
functions, Phys. Rev. A 3, 786 (1971), doi:10.1103/PhysRevA.3.786.
[58] H.-c. Kao, Chiral zero modes in superconducting nanowires with Dresselhaus spin-orbit
coupling, Phys. Rev. B 90, 245435 (2014), doi:10.1103/PhysRevB.90.245435.
[59] S. S. Hegde and S. Vishveshwara, Majorana wave-function oscillations, fermion parity switches, and disorder in Kitaev chains, Phys. Rev. B 94, 115166 (2016),
doi:10.1103/PhysRevB.94.115166.
[60] G. Kells, Many-body Majorana operators and the equivalence of parity sectors, Phys. Rev.
B 92, 081401 (2015), doi:10.1103/PhysRevB.92.081401.
[61] G. Kells, Multiparticle content of Majorana zero modes in the interacting p-wave wire, Phys.
Rev. B 92, 155434 (2015), doi:10.1103/PhysRevB.92.155434.
[62] A. Wi˛eckowski, M. M. Maśka and M. Mierzejewski, Identification of Majorana modes ininteracting systems by local integrals of motion, Phys. Rev. Lett. 120, 040504 (2018),
doi:10.1103/PhysRevLett.120.040504.
[63] T. Koma, Stability of Majorana edge zero modes against interactions, (arXiv preprint)
doi:10.48550/arXiv.2205.11222.
[64] P. Fendley, Strong zero modes and eigenstate phase transitions in the X Y Z/interacting
Majorana chain, J. Phys. A: Math. Theor. 49, 30LT01 (2016), doi:10.1088/17518113/49/30/30LT01.
[65] G. Vionnet, B. Kumar and F. Mila, Level crossings induced by a longitudinal
coupling in the transverse field Ising chain, Phys. Rev. B 95, 174404 (2017),
doi:10.1103/PhysRevB.95.174404.
[66] N. Chepiga and F. Mila, Exact zero modes in frustrated Haldane chains, Phys. Rev. B 96,
060409 (2017), doi:10.1103/PhysRevB.96.060409.
[67] G. Kells, N. Moran and D. Meidan, Localization enhanced and degraded topological order in interacting p-wave wires, Phys. Rev. B 97, 085425 (2018),
doi:10.1103/PhysRevB.97.085425.

33

SciPost Phys. 14, 152 (2023)[68] R. Verresen, R. Moessner and F. Pollmann, One-dimensional symmetry protected topological phases and their transitions, Phys. Rev. B 96, 165124 (2017),
doi:10.1103/PhysRevB.96.165124.
[69] M. Suzuki, Relationship among exactly soluble models of critical phenomena. I, Prog. Theor.
Phys. 46, 1337 (1971), doi:10.1143/PTP.46.1337.
[70] W. Son, L. Amico, R. Fazio, A. Hamma, S. Pascazio and V. Vedral, Quantum phase transition between cluster and antiferromagnetic states, Europhys. Lett. 95, 50001 (2011),
doi:10.1209/0295-5075/95/50001.
[71] N. Chepiga, From Kosterlitz-Thouless to Pokrovsky-Talapov transitions in spinless fermionsand spin chains with next-nearest neighbour interactions, Phys. Rev. Res. 4, 043225(2022), doi:10.1103/PhysRevResearch.4.043225.
[72] O. Foda, A supersymmetric phase transition in Josephson-tunnel-junction arrays, Nucl.
Phys. B 300, 611 (1988), doi:10.1016/0550-3213(88)90615-3.
[73] M. Sitte, A. Rosch, J. S. Meyer, K. A. Matveev and M. Garst, Emergent lorentz symmetry with vanishing velocity in a critical two-subband quantum wire, Phys. Rev. Lett. 102,
176404 (2009), doi:10.1103/PhysRevLett.102.176404.
[74] L. Huijse, B. Bauer and E. Berg, Emergent supersymmetry at the Ising-BerezinskiiKosterlitz-Thouless multicritical point, Phys. Rev. Lett. 114, 090404 (2015),
doi:10.1103/PhysRevLett.114.090404.
[75] I. Peschel, Surface magnetization in inhomogeneous two-dimensional Ising lattices, Phys.
Rev. B 30, 6783 (1984), doi:10.1103/PhysRevB.30.6783.

34