Level crossings induced by a longitudinal coupling in the transverse field Ising chain
Grégoire Vionnet Brijesh Kumar Frédéric Mila




Institute of Physics, École Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland

School of Physical Sciences, Jawaharlal Nehru University (JNU), New Delhi 110067, India.

Institute of Physics, École Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland



We study the effect of 
antiferromagnetic longitudinal coupling on the one-dimensional transverse 
field Ising model with nearest-neighbour couplings. In the topological phase 
where, in the thermodynamic limit, the ground state is twofold degenerate, 
we show that, for a finite system of N sites, the longitudinal coupling 
induces N level crossings between the two lowest lying states as a function 
of the field. We also provide strong arguments suggesting that these N 
level crossings all appear simultaneously as soon as the longitudinal coupling 
is switched on. This conclusion is based on perturbation theory, and a mapping 
of the problem onto the open Kitaev chain, for which we write down the 
complete solution in terms of Majorana fermions.


Level crossings induced by a longitudinal coupling in the transverse field Ising chain
    Frédéric Mila
    December 30, 2023
======================================================================================



The topological properties of matter are currently attracting a considerable 
attention <cit.>. One of the hallmarks of a topologically non 
trivial phase is the presence of surface states. In one dimension, the first 
example was the spin-1 chain that was shown a long time ago to have a
gapped phase <cit.> with two quasi-degenerate low-lying states 
(a singlet and a triplet) on open chains <cit.>. These low-lying 
states are due to the emergent spin-1/2 degrees of freedom at the edges of 
the chains which combine to make a singlet ground state with an almost 
degenerate low-lying triplet for an even number of sites, 
and a triplet ground state with an almost degenerate low-lying singlet when 
the number of sites is odd. In that system, the emergent degrees of freedom 
are magnetic since they carry a spin 1/2, and they can be detected by standard 
probes sensitive to local magnetisation such as NMR <cit.>. 

In fermionic systems, a topological phase is present if the model includes a 
pairing term (as in the mean-field treatment of a p-wave superconductor), 
and the emergent degrees of freedom are two Majorana fermions localised at 
the opposite edges of the chain <cit.>. Their detection is much less 
easy than that of magnetic edge states, and it relies on indirect
consequences such as their impact on the local tunneling density of 
states <cit.>, or the presence of two quasi-degenerate 
low-lying states in open systems. In that respect, it has been suggested to 
look for situations where the low-lying states cross as a function of an 
external parameter, for instance the chemical potential, to prove that there 
are indeed two low-lying states <cit.>. 

In a recent experiment with chains of Cobalt atoms evaporated onto a 
Cu_2N/Cu(100) substrate <cit.>, the presence of level crossings 
as a function of the external magnetic field has been revealed by scanning 
tunneling microscopy, which exhibits a specific signature whenever the 
ground state is degenerate. The relevant effective model for that system is 
a spin-1/2 XY model in an in-plane magnetic field. The exact diagonalisation 
of finite XY chains has indeed revealed the presence of 
quasi-degeneracy between the two lowest energy states, that are well separated 
from the rest of the spectrum, and a series of level crossings between them
as a function of the magnetic field <cit.>.  Furthermore, the 
position of these level crossings is in good agreement with the experimental 
data. It has been proposed that
these level crossings are analogous to those predicted in topological 
fermionic spin chains, and that they can be interpreted as a consequence 
of the Majorana edge modes <cit.>. 

The topological phase of the XY model in an in-plane magnetic field is 
adiabatically connected to that of the transverse field Ising model, in 
which the longitudinal spin-spin coupling (along the field) is switched off. 
However, in the transverse field Ising model, the two low-lying states never 
cross as a function of the field, as can be seen from the magnetisation curve 
calculated by Pfeuty a long time ago <cit.>, and which does not show 
any anomaly. The very different behaviour of the XY model in an in-plane field 
in that respect calls for an explanation. The goal of the present paper is to 
provide such an explanation, and to show that the presence of N level 
crossings, on a chain of N sites, is generic as soon as an antiferromagnetic 
longitudinal coupling is switched on. To achieve this goal, we have studied 
a Hamiltonian which interpolates between the exactly solvable transverse field 
Ising (TFI) and the longitudinal field Ising (LFI) chains.  
The approach that best accounts for these level crossings turns out to  be an 
approximate mapping onto the exactly solvable Kitaev chain, which contains all 
the relevant physics. In the Majorana representation, the level crossings are 
due to the interaction between Majorana fermions localised at each end of the 
chain.

The paper is organized as follows. In section <ref>, we present the 
model and give some exact diagonalisation results on small chains to get an 
intuition of the qualitative behaviour of the spectrum. We show in section 
<ref> that  perturbation theory works in principle but is rather 
limited because of the difficulty to go to high order. We then turn to an 
approximate mapping onto the open Kitaev chain via a mean-field decoupling in 
section <ref>. The main result of this paper is presented in section 
<ref>, namely the explanation of the level crossings in a Majorana 
representation. Finally, we conclude with a a quick discussion of some 
possible experimental realisations in section <ref>.





§ MODEL
 
We consider the transverse field spin-1/2 Ising model with an additional 
antiferromagnetic longitudinal spin-spin coupling along the field, i.e. the 
Hamiltonian

    H=J_x ∑_i=1^N-1S_i^x S^x_i+1 + J_z ∑_i=1^N-1S^z_i S^z_i+1 
    - h ∑_i=1^N S^z_i

with J_z ≥ 0 [This Hamiltonian is equivalent to an XY model in 
an in-plane magnetic field, but we chose to rotate the spins around the 
x-axis so that we recover the usual formulations of the TFI and LFI 
models as special cases.]. 
This model can be seen as an interpolation between the TFI model (J_z=0) 
and the LFI model (J_x=0). The case J_z=J_x corresponds to the effective 
model describing the experiment in Ref. toskovic, up to small 
irrelevant terms [In Ref. toskovic it is explained that 
the ± 3/2 doublet of the spin 3/2 Cobalt adatoms can be 
projected out by a Schrieffer-Wolff transformation due to the strong 
magnetic anisotropy. The resulting effective spin 1/2 model is the one of 
equation (<ref>) with J_x=J_z and additional nearest-neighbour 
out-of-plane and next-nearest-neighbour in-plane Ising couplings. 
These additional terms do not lead to qualitative changes because the model is still
symmetric under a π-rotation of the spins around the 
z-axis, and since their coupling constants are small (∼ 0.1J_x) 
they have only a small quantitative effect in exact diagonalisation results.]. 
Since we will be mostly interested in the parameter range 
0≤ J_z≤ J_x, we will measure energies in units of J_x by setting 
J_x=1 henceforth. The spectrum of the Hamiltonian in Eq. (<ref>) is 
invariant under h→ -h since the Hamiltonian is invariant if we 
simultaneously rotate the spins around the x-axis so that 
S^z_i → -S^z_i ∀ i. Hence, we will in most cases quote the results 
only for h ≥ 0.

The TFI limit of H can be solved exactly by Jordan-Wigner mapping onto a 
chain of spinless fermions <cit.>. In the thermodynamic limit, it is 
gapped with a twofold degenerate ground state for h < h_c=1/2, and undergoes
a quantum phase transition at h=h_c to a non-degenerate gapped ground state 
for h>h_c. The twofold degeneracy when h<h_c can be described by two 
zero-energy Majorana edge modes <cit.>. As a small positive J_z is 
turned on, there is no qualitative change in the thermodynamic limit, except 
that h_c increases with J_z. Indeed, the model is then equivalent to the 
ANNNI model in a transverse field which has been extensively studied before, 
see for example <cit.>. A second order perturbation 
calculation in 1/h yields h_c =1/2+(3/4)J_z+O(J_z^2) for small J_z and 
h_c=1/2+J_z +O(1/J_z) for large J_z <cit.>. Since, for 
J_z ≳ 1, there are other phases arising <cit.>, we shall 
mostly consider J_z ≲ 1 in the following in order to stay in the 
phase with a degenerate ground state.


For a finite size chain, the twofold degeneracy of the TFI model at 0<h<h_c 
is lifted and there is a small non-vanishing energy splitting 
ϵ = E_1-E_0 between the two lowest energy states,  where the 
E_k are the eigenenergies and E_k ≤ E_k+1 ∀ k. This splitting 
is exponentially suppressed with the system length, 
ϵ∼exp(-N/ξ) <cit.>. These two quasi-degenerate states 
form a low energy sector separated from the higher energy states. The spectrum 
for J_z=0 and N=3 is shown in Fig. <ref>a. For J_z > 0, the 
splitting ϵ has an oscillatory behaviour and vanishes for some values 
of h. For N=3, it vanishes once for h>0. See the spectrum for J_z=0.5 
and J_z=1 in Figs <ref>b-c. As J_z becomes large, there is no low 
energy sector separated from higher energy states any more. In the LFI limit, 
J_z→∞, the eigenstates have a well defined magnetisation in the 
z-direction and the energies are linear as a function of h, see 
Fig. <ref>d. In this limit, the level crossings are obvious.
As the field is increased, the more polarised states become favoured, which 
leads to level crossings.







The plots in Fig. <ref> are instructive for very small N but become 
messy for larger chains. In Figs <ref>a-b, we show the spectrum relative 
to the ground state energy, i.e. E_k-E_0, of a chain of N=6 sites for 
J_z=0 and J_z=0.75. The energies E_0 and E_1 are plotted in 
Figs <ref>c-d for the same parameters. The structure of the spectrum is 
similar to the N=3 case, except that now ϵ vanishes at three points 
for h>0. In general, there are N points of exact degeneracy where the 
splitting ϵ vanishes since the spectrum is symmetric under h→ -h. 
This is shown in Fig. <ref> for 2≤ N≤ 8. For N even, there are 
N/2 level crossings for h>0, and for N odd, there are (N-1)/2 level 
crossings for h>0 and one at h=0. 

As shown in Figs <ref>e-f, the level crossings lead to jumps in the 
magnetisation M(h)=- E_0/ h. The number of magnetisation 
jumps turns out to be independent of J_z for 0<J_z<∞, as illustrated 
in Fig. <ref>. In the LFI limit, most of the jumps merge together at 
h=J_z, with an additional jump persisting for even N at 
h=J_z/2 [In the LFI model the lowest energy with a given 
magnetisation  is E_0,M= 0 =-J_z(N-1)/4 and 
E_0,M≠ 0 = E_0,0+ J_z(|M|-1/2) - Mh. Thus for even N, the ground 
state has M=0 for 0<h<J_z/2, M=1 for J_z/2<h<J_z and M=N/2 for 
h>J_z, whereas for odd N the ground state has M=1/2 for 0<h<J_z 
and M=N/2 for h>J_z.]. In this large J_z region, however, there is no 
quasi-degeneracy and the magnetisation jumps indicate level crossings but no 
oscillation in contrast to the small J_z region. Since there are no level 
crossings in the TFI limit, one might expect the number of crossings to 
decrease as J_z decreases. However, the exact diagonalisation results do 
not support this scenario, and hint to all level crossings appearing at the 
same time as soon as J_z≠ 0. This is a remarkable feature that we shall 
explain in the following. 


A useful equivalent representation of the Hamiltonian in Eq. (<ref>) 
in terms of spinless fermions is obtained by applying the Jordan-Wigner 
transformation used to solve exactly the TFI model <cit.>,
 
    S_i^x = 1/2(c_i†  + c_i)exp(ıπ∑_j<ic_j† c_j)
    
    S_i^y = 1/2ı(c_i† - c_i)exp(ıπ∑_j<ic_j† c_j)
    
    S_i^z=c_i† c_i - 1/2,

which yields
      
    H    = 1/4∑_i=1^N-1(c_i†-c_i)(c_i+1†+c_i+1) 
    		- h ∑_i=1^N (c_i† c_i - 1/2)
          
    		+ J_z ∑_i=1^N-1(c_i† c_i 
    		- 1/2)(c_i+1† c_i+1 - 1/2)

where the c_i,c_i† are fermionic annihilation and creation operators. 
This is the Hamiltonian of a spinless p-wave superconductor with 
nearest-neighbour density-density interaction. As for the simpler TFI model, 
the Hamiltonian is symmetric under a π-rotation of the spins around the 
z-axis, S^x_i→ -S^x_i and S^y_i→ -S^y_i in the spin 
language. This leads to two parity sectors given by the parity operator

    P=e^ıπ∑_j=1^N c_j† c_j = (-2)^N S_1^z⋯ S_N^z.

In other words, the Hamiltonian does not mix states with even and odd number 
of up spins, or equivalently with even and odd number of fermions. The ground 
state parity changes at each point of exact degeneracy, and thus alternates as 
a function of the magnetic field for J_z > 0. This can be understood 
qualitatively by looking at Fig. <ref>f. The magnetisation plateaus are 
roughly at M=0,1,2,3. Hence to jump from one plateau to the next, one spin 
has to flip, thus changing the sign of the parity P.









§ PERTURBATION THEORY
 

As a first attempt to understand if the N level crossings develop 
immediately upon switching on J_z, we treat the  
V=J_z ∑_i=1^N-1S^z_i S^z_i+1 term as a perturbation to the 
exactly solvable transverse field Ising model. One may naively expect that 
degenerate perturbation theory is required since the TFI chain has a 
quasi-twofold degeneracy at low field. Fortunately, the two low-energy states 
live in different parity sectors <cit.> that are not mixed by the 
perturbation V. We can therefore apply the simple Rayleigh-Schrödinger 
perturbation theory in the range of parameters we are interested in, i.e. 
J_z ≲ 1.







Writing A_i=c_i† + c_i and B_i=c_i† - c_i, the perturbation can be 
rewritten as V=(J_z/4) ∑_i=1^N-1B_iA_iB_i+1A_i+1.
The unperturbed eigenstates are |m⟩ = Υ†_m|0⟩ where 
|0⟩ is the ground state and the Υ†_m are a product of the 
creation operators corresponding to the Bogoliubov fermions. 
The matrix elements are then

    ⟨n|V|m|=⟩J_z/4∑_i=1^N-1⟨0|Υ_n B_iA_iB_i+1A_i+1Υ†_m|0|⟩

which can be computed by applying Wick's theorem, similarly to how correlation 
functions are found in <cit.>. We computed the effect of V up to third 
order, with the basis of virtual states slightly truncated, namely by keeping 
states with at most three Bogoliubov fermions. Since the more fermions there 
are in a state, the larger its energy, we expect this approximation to be 
excellent.

As shown in Fig. <ref>, the number of crossings increases with the order 
of perturbation, and to third order in perturbation, the results for N=6 
sites are in qualitative agreement with exact diagonalisations. From the way 
level crossings appear upon increasing the order of perturbation theory, one 
can expect to induce up to 2m+1 level crossings if perturbation theory
is pushed to order m, see Fig. <ref>. So these results suggest that the 
appearance of level crossings is a perturbative effect, and that, for a given 
size N, pushing perturbation theory to high enough order will indeed lead to 
N level crossings for small J_z. However, in practice, it is impossible to 
push perturbation theory to very high order. Indeed, the results at order 3 
are already very demanding. So, these pertubative results are encouraging, but 
they call for an alternative approach to actually prove that the number of 
level crossings is indeed equal to N, and that these level crossings appear 
as soon as J_z is switched on.



§ FERMIONIC MEAN-FIELD APPROXIMATION
 




In the fermionic representation, Eq. (<ref>), there is a quartic 
term that cannot be treated exactly. Here, we approximate it by mean-field 
decoupling. In such an approximation, one assumes the system can be well 
approximated by a non-interacting system (quadratic in fermions) with 
self-consistently determined parameters. For generality, we decouple the 
quartic term in all three mean-field channels consistent with Wick's theorem,

    c†_i   c_ic†_i+1c_i+1≈
       ⟨c†_ic_i|c⟩†_i+1c_i+1 + ⟨c†_i+1c_i+1|c⟩†_ic_i - ⟨c†_ic_i|⟨%s|%s⟩⟩c†_i+1c_i+1
    
    -    ⟨c†_ic†_i+1|c⟩_ic_i+1 - ⟨c_ic_i+1|c⟩†_ic†_i+1  +⟨c†_ic†_i+1|⟨%s|%s⟩⟩c_ic_i+1
     
    +   ⟨c†_ic_i+1|c⟩_ic†_i+1 +⟨c_ic†_i+1|c⟩†_ic_i+1  - ⟨c†_ic_i+1|⟨%s|%s⟩⟩c_ic†_i+1.

Here, ⟨.|$⟩ denotes the ground state expectation value.
The3N-2self-consistent parameters⟨c†_ic_i|$⟩, 
⟨c†_ic†_i+1|$⟩ and⟨c†_ic_i+1|$⟩ can be found 
straightforwardly by iteratively solving the quadratic mean-field Hamiltonian. 



As it turns out, it is more instructive to consider only three self-consistent 
parameters. To do so, we solve the mean-field approximation of the 
translationally invariant Hamiltonian (c_N+1=c_1), 

    H'    = ∑_i=1^N{1/4  
    	(c_i†-c_i)(c_i+1†+c_i+1) - h(c_i† c_i - 1/2) . 
    	
       . + J_z (c_i† c_i - 1/2)(c_i+1† c_i+1 - 
    	1/2) }
       ≈∑_i=1^N{ 
    	-μ c_i† c_i +( t c_i+1† c_i +  h.c. )
    	 - (Δ c_i+1† c†_i +  h.c. ) }
          +  const,

where μ=h+J_z(1-2⟨c†_ic_i|)⟩, 
t=1/4 - J_z⟨c†_ic_i+1|$⟩ andΔ=1/4-J_z⟨c_ic_i+1|$⟩ are determined self-consistently. 
These parameters are found to be real, and are shown in Fig. <ref> 
as a function h, J_z and N.

Using these self-consistent parameters, the Hamiltonian in 
Eq. (<ref>) is then approximated by the following mean-field 
problem on an open chain:

    H_ MF   = -∑_i=1^N μ(c†_ic_i-1/2) 
    
          + ∑_i=1^N-1[(t c†_i+1c_i +  h.c.) 
    - ( Δ c†_i+1c†_i +  h.c. )   ],

up to an irrelevant additive constant [In the periodic chain 
used to get the mean-field parameters, there is a level crossing when 
μ=2t. To get good agreement with exact diagonalisation results and avoid a small 
discontinuity, we need to compute the expectation values in the state 
adiabatically connected to the ground state at μ<2t. Thus for μ>2t, 
the ⟨.|$⟩ are not computed in the ground state, but in the first 
excited state. Since all the level crossings arise forμ< 2t, this has 
no influence on the following discussion..


    < g r a p h i c s >
Critical fields, h_ crit, where the degeneracy is exact, 
as a function of J_z in the self-consistent mean-field approximation 
(<ref>) (blue crosses) compared to the exact diagonalisation result 
(red squares) for (a) N=6 and (b) N=7.


    < g r a p h i c s >
Energy splitting ϵ=E_1-E_0 as a function of h in the 
self-consistent mean-field approximation (<ref>) (blue solid line) 
compared to the exact diagonalisation result (red dashed line) for J_z=0.5
 and (a) N=6 and (b) N=7.



Since the self-consistent parameters are almost independent of the system size 
(see Fig. <ref>c), the boundaries are not very important and the bulk 
contribution is determinant. This partly justifies the approximation of 
playing with the boundary conditions to get the approximate model 
(<ref>) with just three self-consistent parameters. This approximation 
is also justified by the great quantitative agreement with the exact 
diagonalisation results for the critical fields forJ_z ≲0.8(see Fig. <ref>), and to a lesser extent for the energy splittingϵ=E_1-E_0between the two lowest energy states, see Fig. <ref>. 

ForNodd, the degeneracy ath=0is protected by symmetry for anyJ_zin the Hamiltonian (<ref>). Indeed, under the transformationS_i^z →-S_i^z ∀i, the parity operator transforms asP→(-1)^N P. Hence, forNodd andh=0, the ground state has to be 
twofold degenerate. As can be seen in Fig. <ref>b, the critical fieldh=0at lowJ_zevolves to a non-zero value for largeJ_z, thus showing 
that this symmetry is broken by the mean-field approximation (<ref>). 
The discrepancy is, however, small forJ_z ≲0.8as can also be seen 
in Fig. <ref>b. 

We observe from Fig. <ref>a that as a function of magnetic field, the 
parameterstandΔare almost constant, whereasμis almost 
proportional toh. Thus, we can understand the physics of the level 
oscillations by forgetting about the self-consistency and consideringμ,tandΔas free parameters, i.e. by studying the open Kitaev 
chain <cit.>, where the level crossings happen asμis tuned.
Compared to the TFI model for whichΔ=t, the main effect ofJ_z>0is to make0 < Δ< t, which, as we shall see in the next section, is 
the condition to see level oscillations.

Such a mapping between the two lowest lying energy states of the interacting 
Kitaev chain and of the non-interacting Kitaev chain can be made rigorous for a 
special value ofh>0, provided the boundary terms in 
equation (<ref>) are slighty modified <cit.>. But this 
particular exact case misses out on level-crossing oscillations.



§ LEVEL OSCILLATIONS AND MAJORANA FERMIONS

We define2NMajorana operatorsγ'_i, γ_i”as: 

    γ'_i = c_i+ c†_i
    γ”_i = -ı (c_i -  c†_i)

which satisfyγ'_i†=γ'_i,γ”_i†=γ”_i,{γ'_i,γ”_j}=0and{γ'_i,γ'_j}= 
{γ”_i,γ”_j}=2δ_ij. Since theμ, t, Δare real, theH_MFof Eq. (<ref>) reads

    H_ MF=   ı/2∑_i=1^N-1[ -(t+Δ)γ”_iγ'_i+1 
    + (t-Δ )γ'_iγ”_i+1] 
        -ıμ/2∑_i=1^N γ'_iγ”_i= 
    ı/2∑_i,j=1^Nγ'_iM_ijγ”_j.

From the singular value decomposition ofM, we writeM=UΣV^T, 
whereUandVare orthogonal matrices andΣ=diag(ϵ_1,…, ϵ_N)with realϵ_iand|ϵ_i|≤|ϵ_i+1| ∀i. Thus, the Hamiltonian reads

    H_ MF   = ı/2∑_i,j,k=1^Nγ'_i U_ikϵ_k V^T_kjγ”_j
    =ı/2∑_k=1^Nϵ_k γ̃'_k γ̃”_k 
       =∑_k ϵ_k (η†_kη_k-1/2)

where

    γ̃'_k = ∑_i=1^N γ'_i U_ik,     γ̃”_k = ∑_i=1^N γ”_i V_ik

are the rotated Majorana operators, and theη_k = 1/2 (γ̃'_k+ıγ̃”_k) are 
fermionic annihilation operators corresponding to the Bogoliubov 
quasiparticles.


 As derived in Appendix, in general the Majorana operators,γ̃'_kandγ̃”_k, are of the form
 
    γ̃'_k = ∑_j (a_+ x_+^j + b_+ x_+^N+1-j 
    		+ a_- x_-^j + b_- x_-^N+1-j ) γ'_j
    γ̃”_k = ∑_j ( a_+ x_+^N+1-j + b_+ x_+^j 
    		+a_- x_-^N+1-j + b_- x_-^j) γ”_j

 where thex_±,a_±andb_±are functions of the energyϵ_kwhich is quantised in order to satisfy the boundary 
 conditions. On can easily solve numerically the nonlinear equation for 
 theϵ_k. Here, we will instead focus on a simple analytical 
 approximation forγ̃'_1, γ̃”_1andϵ_1which works well to discuss the level crossings,
 and is equivalent to the Ansatz given in <cit.>.

From Eqs. (<ref>) and (<ref>), we see that forϵ=0, 
we have eithera_±=0orb_±=0. Without loss of generality, we can 
chooseb_±(ϵ=0) = 0. Since we expectϵ_1 ≪1, we 
approximate 
 
    b_± (ϵ_1) ≈ b_± (0)=0

 and 
 
    x_± (ϵ_1) ≈ x_±(0) = 
     	μ±√(μ^2-4t^2+4Δ^2)/2(t+Δ),

which yields
  
    γ̃'_1    ≈∑_j (a_+ x_+^j  + a_- x_-^j ) γ'_j
    γ̃”_1    ≈∑_j ( a_+ x_+^N+1-j +a_- x_-^N+1-j )γ”_j

 with∑_j (a_+ x_+^j  + a_- x_-^j )^2=1.
 The boundary conditions (<ref>) now read 
 
    a_+ + a_-    =    0 
    
    	 a_+ x_+^N+1 + a_- x_-^N+1   =    0

and in general cannot be both satisfied unlessϵ_1=0exactly. 

If|x_±|<1,γ̃'_1is localised on the left side of the chain 
with its amplitude∼e^-j/ξasj≫1withξ=-1/ln(max(|x_+|,|x_-|)). Furthermore,γ̃”_1is related 
toγ̃'_1by the reflection symmetryj→N+1-j.
Thus, in the thermodynamic limit, the boundary condition (<ref>) is 
irrelevant andϵ_1 →0asN→∞. 
Similarly, if|x_±| >1the boundary condition (<ref>) becomes 
irrelevant in the thermodynamic limit.
However, if|x_+|>1and|x_-|<1,  or|x_+|<1and|x_-|>1, thenγ̃'_1,γ̃”_1have significant weight on both 
sides of the chain and both boundary conditions (<ref>) and 
(<ref>) remain important in the thermodynamic limit. Hence, the 
approximationϵ_1 ≈0is bad, indicating a gapped system.
As discussed in <cit.>, for|μ|<2|t|we have either|x_±|<1or|x_±|>1which yieldsϵ_1=0in the thermodynamic limit. This 
is the topological phase with a twofold degenerate ground state. For a finite 
system, however, the boundary conditions (<ref>) and (<ref>) are 
in general not exactly satisfied and the system is only quasi-degenerate with 
a gapϵ∼e^-N/ξ.
For|μ|>2|t|, either|x_+|>1and|x_-|<1, or|x_+|<1and|x_-|>1, and the system is gapped. 

In the topological phase,|μ|<2|t|, there are parameters for which the 
boundary conditions (<ref>) can be exactly satisfied even forN<∞and thusϵ_1=0exactly. In such a case, there is an 
exact zero mode even for a finite chain. This was previously discussed in 
Ref. <cit.>, as well as in <cit.> where a more general method that
applies to disordered systems is described. Ifx_±∈ℝ, it 
is never possible to satisfy the boundary conditions (<ref>) and 
therefore the quasi-gap is always finite,ϵ_1 ≠0. However, ifx_+=re^ıϕ ∉ℝ, Eq. (<ref>) yieldsx_-=x_+^*and(x_+^N+1 - x_-^N+1) ∝r^N+1sin[(N+1)ϕ]. 
Thus it may happen for specific parameters thatϵ_1=0exactly. This 
degeneracy indicates a level crossing. The phaseϕ, defined for|μ|<μ_c=2√(t^2-Δ^2), is given by

    tanϕ  = √((μ_c/μ)^2-1).

It thus goes continuously fromϕ(μ=0^+)=π/2toϕ(μ→μ_c) →0. Hence, there are critical chemical potentials,0 ≤μ_⌈N/2⌉ <  …< μ_m <…<μ_1<μ_c, 
such thatϕ(μ=μ_m)=πm/N+1(see Fig. <ref>a). 
For these criticalμ_m, the system is exactly degenerate, i.e.ϵ_1=0. In the TFI limit, we haveΔ=tandμ_c=0, 
thus there are no level crossings.


    < g r a p h i c s >
(a) Phase ϕ(μ) of x_+=re^ıϕ within the approximation 
(<ref>) for several Δ with t=1 and N=6. The horizontal black 
dotted lines indicate the values ϕ=π m/N+1. (b) Splitting 
E_1-E_0=|ϵ_1| in the Kitaev chain calculated exactly solving 
numerically the full self-consistent equations described in Appendix 
(blue solid line) and with the analytical approximate result 
in Eq. (<ref>) (red dashed line) for N=6, t=1 and Δ=0.3.


For|μ|<2|t|, writingx_+=re^ıϕwithr>0, we have

    ϵ_1    =Σ_11=(U^T M V)_11
       ≈ 4(t+Δ) a_+^2 r^N+2sin(ϕ) sin[(N+1)ϕ],

where we used the approximations (<ref>), (<ref>) and the 
boundary condition (<ref>) [respectively (<ref>)] whentΔ>0(respectivelytΔ<0), since in this case|x_±|<1(respectively|x_±|>1). Note thatϕ(-μ) = ϕ(μ)-π, 
and thusϵ_1is an odd function ofμfor oddNand an 
even function ofμfor evenN. Sinceϵ_1changes sign 
wheneversin((N+1)ϕ)=0, the degeneracy points indicate level crossings. 
This approximate description works extremely well, as shown in 
Fig. <ref>b forΔ=0.3t. Becauseϕtakes all the values 
in]0,π/2]for0 < μ<μ_c, and in]-π,-π/2]for-μ_c < μ< 0, there are either exactlyNlevel crossings as a function 
ofμif0<μ_c∈ℝ, i.e. if|Δ|<|t|, and no zero 
level crossing otherwise. At the points of exact degeneracy,b_±(ϵ=0)=0, the zero-mode Majorana fermions are localised on 
opposite sides of the chain. When the degeneracy is not exact, however,b_±(ϵ≠0) ≠0and the zero-mode Majorana fermions mix 
together to form Majoranas localised mostly on one side but also a little 
bit on the opposite side.

In the XY model in an out-of-plane magnetic field, which is equivalent to 
the non-interacting Kitaev chain <cit.>, these level crossings lead 
to an oscillatory behaviour of the spin correlation functions <cit.>. 
In the context of p-wave 
superconductors, the level oscillations described above also 
arise in more realistic models and are considered one of the hallmarks of 
the presence of topological Majorana fermions <cit.>. Although 
it is still debated whether Majorana fermions have already been observed, 
strong experimental evidence for the level oscillations was reported in 
<cit.>.

Coming back to the mean-field Hamiltonian of Eq. (<ref>), we can get 
the phaseϕwithin the approximation (<ref>), i.e. the phase ofx_+(ϵ=0), as a function of the physical parametersh, J_zsince 
we know how the self-consistent parametersμ, t, Δdepend on them.
We plot in Fig. <ref> the phaseϕas a function ofhfor 
severalJ_zwhich yields a good qualitative understanding of the sudden 
appearance ofNlevel crossings as soon asJ_z > 0. As previously 
discussed, the self-consistent parameters are almost independent ofNand 
therefore the curvesϕ(h)are almost independent ofNas well. The 
main effect ofNis to change the conditionϕ(μ=μ_m)=πm/N+1for the boundary condition in 
Eq. (<ref>) to be satisfied and thus for the system to be exactly 
degenerate.



    < g r a p h i c s >
Phase ϕ(h) of x_+(ϵ=0) based on the self-consistent 
parameters μ, t, Δ of the mean-field decoupling for several J_z 
and (a) N=6, (b) N=9, (c) N=12. The horizontal black dotted lines
indicate the values ϕ=π m/N+1.


§ SUMMARY


The main result of this paper is that the level crossings between the two 
lowest energy eigenstates of the XY chain in an in-plane magnetic field 
are more generally a fundamental feature of the transverse field Ising chain 
with an antiferromagnetic longitudinal coupling  
howsoever small. These points of level crossings (twofold degeneracy) 
correspond to having Majorana edge modes in a Kitaev chain onto which the 
problem can be approximately mapped. The level crossings of the XY chains 
have been observed experimentally in <cit.> by scanning tunneling 
microscopy on Cobalt atoms evaporated onto a Cu_2N/Cu(100) substrate. 
By varying the adsorbed atoms and the substrate, it should be possible to 
vary the easy-plane and easy-axis anisotropies, and thus to explore the 
exact degeneracy points for various values of the longitudinal coupling.
The possibility to probe the two-fold degeneracy of this family of spin 
chains is important in view of their potential use for universal quantum
computation <cit.>.
Besides, one could also realise the spinless fermionic Hamiltonian 
(<ref>) in an array of Josephson junctions as described in 
<cit.>. The advantage of this realisation is that it allows a 
great flexibility to tune all the parameters of the model. We hope that 
the results of the present paper will stimulate experimental investigations 
along these lines.


We acknowledge Somenath Jalal for useful discussions and the Swiss National 
Science Foundation for financial support. B.K. acknowledges the financial 
support under UPE-II and DST-PURSE programs of JNU.
*



§ MAJORANA SOLUTIONS OF THE KITAEV CHAIN


To solve the Kitaev chain (<ref>), we need to find the singular 
value decomposition of

    M=[  -μ τ_-   0   ⋯; τ_+  -μ τ_-   0   ⋯;   0 τ_+  -μ τ_-   0   ⋯;           ⋱   ⋱   ⋱    ;       ⋯   0 τ_+  -μ τ_-;           ⋯   0 τ_+  -μ ]

withτ_±= t±Δ, i.e. find orthogonal matricesU,Vand a real diagonal matrixΣsuch thatM=UΣV^T. Writingu⃗_kandv⃗_kthek^thcolumns ofUandVrespectively, they satisfy

    Mv⃗_k    = ϵ_k u⃗_k 
    u⃗_k^T M    = ϵ_k v⃗_k^T.


Let's find two unit-norm column vectorsu⃗, v⃗andϵsuch thatMv⃗=ϵu⃗andu⃗^TM=ϵv⃗^T. 
First we forget about the normalisation and boundary conditions and focus 
on the secular equation. 
Setting the components ofu⃗, v⃗asu_j=a x^jandv_j = b x^j, 
we have

    M v⃗   = b/aτ_+ -μ x +τ_- x^2/xu⃗+  b.t.
    u⃗^T M    =a/bτ_- -μ x +τ_+ x^2/xv⃗^T +  b.t.

where b.t. stands for boundary terms. Hence,uandvsatisfy the secular 
equation provided 

    b/a=√(τ_- -μ x +τ_+ x^2/τ_+ -μ x +τ_- x^2)

and

    ϵ = 1/x√((τ_- -μ x +τ_+ x^2)
    	(τ_+ -μ x +τ_- x^2)).

Because of the reflection symmetryj →N+1-j, ifxis a solution of 
equation (<ref>) for someϵ, then1/xis also a solution. 
Assumingϵknown, the solutions arex_±,1/x_±and satisfy

    0    = ϵ^2 x^2 - (τ_- -μ x +τ_+ x^2)
    		(τ_+ -μ x +τ_- x^2) 
       ∝ (x-x_+)(x-1/x_+)(x-x_-)(x-1/x_-)

which by identification yields, writingρ_±= x_±+ 1/x_±,

    x_±   = 1/2(ρ_± +√(ρ_±^2 - 4)),
    ρ_±   = μ t ±√((t^2-Δ^2)ϵ^2 + 
    		Δ^2(μ^2-4t^2+4Δ^2))/t^2-Δ^2.

Taking into account the reflection symmetry, the general form of the 
components ofu⃗, v⃗is thus
 
    u_j = a_+ x_+^j + b_+ x_+^N+1-j + a_- x_-^j + b_- x_-^N+1-j
    
    		v_j =  a_+ x_+^N+1-j + b_+ x_+^j +a_- x_-^N+1-j + b_- x_-^j

with the ratiosb_+/a_+andb_-/a_-given by equation (<ref>) 
withx=x_+andx=x_-respectively.

Furthermore, we have the boundary conditions

    a_+ + b_+x_+^N+1 + a_- + b_-x_-^N+1   = 0 
    
    	a_+ x_+^N+1 + b_+ + a_-x_-^N+1 + b_-   = 0

which set the ratioa_-/a_+and give the quantisation condition on the 
energiesϵ_k.
The last degree of freedom, saya_+, is then set by normalisingu⃗(from equation (<ref>), u⃗ =  v⃗ ).

Note that for the special casest= Δandμ=0, we haveγ̃'_1 = γ'_1andγ̃”_1=γ”_Nwithϵ_1=0. We have a similar result fort=-Δandμ=0. 
For these two cases, the general formalism described above does not apply 
since it yieldsx_±=0,±∞.

99kane
M. Z. Hasan and C. L. Kane, Rev. Mod. Phys.  82, 3045 (2010).
zhang
X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys.  83, 1057 (2011).
haldane
F. D. M. Haldane, Phys. Lett. A  93, 464 (1983).
kennedy
T. Kennedy, J. Phys. Cond. Mat.  2, 5737 (1990).
tedoldi
F. Tedoldi, R. Santachiara, and M. Horvatić, Phys. Rev.
Lett.  83, 412 (1999).
kitaev
A. Y. Kitaev. Phys.-Usp.  44 131 (2001).
mourik 
V. Mourik, K. Zuo, S. M. Frolov, S. R. Plissard, 
E. P. A. M. Bakkers, and L. P. Kouwenhoven, 
Science  336, 1003 (2012).
nadj-perge 
S. Nadj-Perge, I. K. Drozdov, J. Li, H. Chen, S. Jeon, J. Seo, 
A. H. MacDonald, B. A. Bernevig, and A. Yazdani, Science  346, 602 (2014).
sarma
S. Das Sarma, J. D. Sau, and T. D. Stanescu,
Phys. Rev. B  86, 220506(R) (2012).
toskovic
R. Toskovic, R. van den Berg, A. Spinelli, I. S. Eliens,
B. van den Toorn, B. Bryant, J.-S. Caux, and A. F. Otte,
Nat. Phys.  12, 656 (2016).
dmitriev
D. V. Dmitriev, V. Y. Krivnov, A. A. Ovchinnikov, and A. Langari, 
J. Exp. Theor. Phys. 95, 538 (2002).
mila
F. Mila, Nat. Phys.  12, 633 (2016).
pfeuty
 P. Pfeuty, Ann. Phys.  57, 79 (1970).
chakrabarti
S. Suzuki, J.-i. Inoue, and B. K. Chakrabarti,  Quantum Ising Phases and 
Transitions in Transverse Ising Models (Springer, Lecture Notes in Physics, 
Vol. 862 (2013)).
jalal
S. Jalal and B. Kumar, Phys. Rev. B  90, 184416 (2014).
rujan
P. Ruján, Phys. Rev. B  24, 6620 (1981).
hassler
F. Hassler and D. Schuricht, New J. Phys.  14, 125018 (2012).
lieb
E. Lieb, T. Schultz, and D. Mattis, Ann. Phys.  16, 407 (1961).
katsura
H. Katsura, D. Schuricht, and M. Takahashi,
Phys. Rev. B  92, 115137 (2015).
kao
H.-C. Kao, Phys. Rev. B  90, 245435 (2014).
hedge
S. S. Hegde, and S. Vishveshwara, Phys. Rev. B  94,
115166 (2016).
barouch
E. Barouch and B. M . McCoy, Phys. Rev.  A3, 786 (1971).
loss
D. Rainis, L. Trifunovic, J. Klinovaja, and D. Loss, 
Phys. Rev. B  87, 024515 (2013).
markus
 	S. M. Albrecht, A. P. Higginbotham, M. Madsen, F. Kuemmeth, 
 	T. S. Jespersen, J. Nygård, P. Krogstrup, and  	C. M. Marcus, 	
	Nature  531, 206 (2016).
loss2 Y. Tserkovnyak and D. Loss, Phys. Rev. A  84, 032333 (2011).]