This is the html version of the file https://arxiv.org/abs/2310.11243.
Google automatically generates html versions of documents as we crawl the web.
arXiv:2310.11243v1 [cond-mat.mes-hall] 17 Oct 2023
  Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                
Page 1
arXiv:2310.11243v1 [cond-mat.mes-hall] 17 Oct 2023
Coexistence of magnetic and topological phases in spin-1/2 anisotropic extended XY
chain
Rakesh Kumar Malakar1, and Asim Kumar Ghosh1,
1Department of Physics, Jadavpur University, 188 Raja Subodh Chandra Mallik Road, Kolkata 700032, India
In this study, a spin-1/2 anisotropic extended XY chain has been introduced in which both time
reversal and SU(2) symmetries are broken but Z2 symmetry is preserved. System exhibits the faithful
coexistence of magnetic and topological superconducting phases even in the presence of transverse
magnetic field. Location of those phases in the parameter space has been determined precisely.
Quantum phase transition is noted at zero magnetic field, as well as magnetic long range order is
found to withstand magnetic field of any strength. Exact analytic results for spin-spin correlation
functions have been obtained in terms of Jordan Wigner fermionization. Existence of long range
magnetic order has been investigated numerically by finding correlation functions as well as the
Binder cumulant in the ground state. Dispersion relation, ground state energy, and energy gap are
obtained analytically. In order to find the topologically nontrivial phase, sign of Pfaffian invariant
and value of winding number have been evaluated. Both magnetic and topological phases are robust
against the magnetic field and found to move coercively in the parameter space with the variation of
its strength. Long range orders along two orthogonal directions and two different topological phases
are found and their one-to-one correspondence has been established. Finally casting the spinless
fermions onto Majorana fermions, properties of zero energy edge states are studied. Three different
kinds of Majorana pairings are noted. In the trivial phase, next-nearest-neighbor Majorana pairing
is found, whereas two different types of nearest-neighbor Majorana pairings are identified in the
topological superconducting phase.
PACS numbers:
I. INTRODUCTION
Quantum fluctuation driven physical phenomena are so
widespread and exotic that they have constituted several
branches of study within the condensed matter physics.
Among them the most important branch is quantum
phase transition (QPT)13. QPT deals with the tran-
sition between phases of different types of quantum or-
ders driven solely by quantum fluctuations as they oc-
cur at zero-temperature (T =0). Nowadays, there have
been an upsurge of investigations on the topological or-
ders of quantum systems and transition between different
topological phases. So, it is pertinent to understand the
interplay between QPTs and topological transitions as
well as their possible interconnections. The most studied
quantum systems which exhibit QPTs are spin-1/2 one-
dimensional (1D) Ising and anisotropic XY models un-
der the transverse magnetic field49. Exact analytic form
of spin-spin correlation functions and other quantities
had been obtained by converting spin operators to spin-
less fermions by Jordan-Wigner (JW) transformations10.
Those systems exhibit magnetic long-range-order (LRO)
and undergoes transition to disorder phase at a specific
value of magnetic field, where spin fluctuations play cru-
cial role for the QPTs. Signature of QPT has been ob-
served in ferromagnet CoNb2O6, ferroelectric KH2PO4,
and other quantum matters1113.
Quantum fluctuation led another important phe-
nomenon which is known as entanglement exhibits un-
usual feature in the vicinity of QPT point. For examples,
derivatives of concurrences for anisotropic XY chain in
transverse field show logarithmic divergences at the tran-
sition points14,15. Entanglement in the ground state of
the anisotropic XY chain has been evaluated in terms of
von Neumann entropy. The constant entropy surfaces are
either ellipse or hyperbolas and they are found to meet
at the QPT points16. Entanglement has been recognized
as the key feature for the development of quantum com-
putation and communication17.
Investigation on topological matters have begun with
the discovery of Quantum Hall Effect by K Von
Klitzing18. In this phenomenon, Hall conductivity of two-
dimensional (2D) electron gas got quantized under strong
magnetic field, which was subsequently characterized by
the topological invariant known as Thouless-Kohmoto-
Nightingale-Nijs (TKNN) number19. This branch of
study has been enriched further with the discovery of
quantum anomalous Hall effect, where quantization of
Hall Conductivity was made possible only with the
time reversal symmetry (TRS) breaking complex hop-
ping terms, replacing the true magnetic field20. Quan-
tum matters exhibiting this particular quantized feature
are termed as topological insulator. This phenomenon
can be understood easily with the help of paradigmatic
Su-Schrieffer-Heeger (SSH) model which is nothing but a
1D tight-binding system composed of two-site unit cells
with staggered hopping parameters2123. The nontrivial
phase is characterized by a nonzero topological invariant
known as winding number which is determined by inte-
grating the Berry curvature over the 1D Brillouin zone
(BZ). Band gap must be nonzero for the nontrivial phase
which is accompanied additionally by the symmetry pro-
tected zero energy boundary states those are found lo-
calized on the edges of the open chain. Number of edge
states in a particular phase is determined by the bulk-
boundary correspondence rule24. However, at the topo-
logical transition point band gap must vanish.
The same topological order has been achieved sepa-
rately through another investigation in terms of 1D Ki-
taev model, which opens up another branch of study
known as topological superconductivity25. In this case,
band gap of the system can be identified specifically to the

Page 2
2
superconducting gap for Cooper pairing of parallel spin.
After expressing fermionic operators in terms of Majo-
rana fermions, different types of Majorana pairings are
found to emerge in the trivial and topological phases. Im-
portance of topological matter attributes to the fact that
topological robustness present in the nontrivial phase pro-
tects these states from any sort of imperfection present
in the materials. This robustness which is essentially led
by quantum fluctuations of the topological matter may
provide higher efficiency in electronic transport for the
development of quantum processing devices26,27. There-
fore, better understanding on the role of quantum fluctu-
ations is required where both long-range quantum orders
and topological phases are found to emerge simultane-
ously and coexist.
In order to study the effect of spin fluctuations for
topological matter, a 1D spin-1/2 model composed of
anisotropic XY chain and a three spin term has been
introduced in this work, and its property has been inves-
tigated extensively under the presence of transverse mag-
netic field. Faithful coexistence of magnetic and topolog-
ical superconducting phases has been established even in
the presence of transverse magnetic field. The system has
been solved exactly by expressing the spin operators in
terms of JW fermions where analytic forms of spin-spin
correlations functions have been obtained. In addition,
magnetic property of the system has been studied nu-
merically by estimating the ground state correlations as
well as the Binder cumulant of the magnetic LRO for a
finite chain. Topological phases of the system have been
characterized analytically by obtaining the eigenvectors
of quasiparticle dispersion spectra. Magnetic phase is
detected in the ground state, while the topological phase
is found in the single particle spinless fermionic excita-
tion. System exhibits two distinct magnetic LROs as
well as two different topologically nontrivial supercon-
ducting orders and they are found to coexist in such a
way that one-to-one correspondence between magnetic
and topological phases could be established. Comprehen-
sive phase diagrams for magnetic and topological phases
are drawn. Transition between different magnetic phases
is also associated to that between respective topological
phases. Therefore, this model serves as a prototype ex-
ample where quantum fluctuation driven magnetic and
topological phases are found to emerge and coexist. Pe-
culiarities of those phases attribute to the fact that they
could withstand the magnetic field of any strength while
transitions between them would occur even without the
magnetic field.
There is plenty of materials which exhibit the coex-
istence of magnetic and superconducting orders. The
phases of superconductivity and weak ferromagnetism
have been found to coexist in ErRh4B4
28. In addition,
experimental data for HoMo6S8 and HoMo6Se8 has given
strong evidence in favor of the coexistence of supercon-
ductivity and ferromagnetic ordering below transition
temperature29. The same coexistence has been observed
further in the materials UGe2, ZrZn2 and URhGe3032.
Theoretical investigation on the coexistence of p-wave su-
perconductivity and itinerant ferromagnetism for those
materials has been carried out33. It has been indi-
cated that strong magnetic fluctuation induces analogous
gap equations for the superconductivity. The coexis-
tence of ferromagnetism and superconductivity in the 2D
LaAlO3/SrTiO3 interface has been confirmed by magne-
tometry and susceptometry images34. Recently, coexis-
tence of superconducting state separately with ferro and
antiferromagnetic phases has been discovered in magneti-
cally anisotropic system (Eu, La)FeAs2
35. So in this con-
text, it is pertinent to investigate the interplay of mag-
netic and topological superconducting phases in great de-
tails.
The article has been arranged in the following order.
Hamiltonian is introduced in the Sec II, as well as sym-
metries of the model have been described. Various phys-
ical quantities are defined in order to study the character
of magnetic LRO. Properties of the Hamiltonian at sev-
eral limits have been studied in great details and com-
pared with the known results. Magnetic properties of
the system has been investigated in the Sec III, which
begins with the results for four-spin plaquette and fol-
lowed by numerical results obtained by exact diagonal-
ization and analytic results by JW fermionization. Dis-
persion relation, band gap, spin-spin correlation func-
tions, Binder cumulant and magnetic phase diagram have
been obtained. Topological properties of the system have
been presented in Sec IV. Values of Pfaffian invariant
and winding number for different topological phases have
been estimated here. Edge states are determined and
their properties in terms of Majorana pairing have been
studied. Topological phase diagram has been drawn. Fi-
nally, an extensive discussion based on all the results has
been made available in Sec V.
II. ANISOTROPIC EXTENDED XY CHAIN IN
TRANSVERSE MAGNETIC FIELD
Hamiltonian (Eq. 1) of the 1D anisotropic extended
XY model in the presence of transverse magnetic field is
written as
H =
N
j=1 [J((1 + γ)Sx
j Sx
j+1 + (1 − γ)Sy
j Sy
j+1)
+J(Sx
j Sx
j+2 + Sy
j Sy
j+2)Sz
j+1 + hzSz
j ],
(1)
where, N is the total number of sites and J is the near-
est neighbor (NN) exchange interaction strength. γ is the
anisotropic parameter while Jis the three-spin exchange
interaction strength. Sα
j , α = x, y, z, is the α-component
of spin-1/2 operator at the site j, and hz is the strength of
magnetic field acting along the z direction. Hamiltonian
breaks the rotational symmetry, U(1), about any direc-
tions, since [H, Sα
T] = 0, where Sα
T is the α-component of
the total spin, ST. The TRS is broken due to the presence
of odd-spin terms, which means the Hamiltonian retains
TRS when J= 0, and hz = 0. Symmetries of H in the
spin space under five different transformations and their
consequences are studied extensively which are discussed
below.

Page 3
3
A. Symmetries of H:
Hamiltonian obeys the symmetry relation,
UzH(J, J,γ,hz) Uz = H(−J, J,γ,hz),
(2)
where Uz = ∏
N
j=1 eiπjSz
j , and when N is assumed even.
The operator Uz rotates the spin vector at the j-th site
about the z-axis by the angle . This symmetry claims
that energy spectrum of H must remain unchanged upon
sign inversion of J. In the same way, it can be shown
that H satisfies the relation,
VβH(J, J,γ,hz)V
β = H(J, J,γ, hz),
(3)
where Vβ = ∏
N
j=1 eiπSβ
j , β = x, or y. Vβ performs ro-
tation of the every spin vector about the β-axis by the
angle π. It signifies that energy spectrum is unaltered
upon simultaneous sign inversion of both Jand hz. As
a consequence, Hamiltonian remains invariant under the
combined operations, W = VβUz, since
WH(J, J,γ,hz)W= −H(J, J,γ,hz).
(4)
The relation, WHW= −H, corresponds to the fact
that energy spectrum of H inherits the inversion symme-
try around the zero energy, or, it reflects the particle-hole
like symmetry of the system. It further implies that zero
energy states must appear in pair if the spectrum pos-
sesses them. This feature has a special importance in the
context of nontrivial topological phase where the emer-
gence of pair of zero-energy edge states in the 1D open
system corresponds to the nonzero value of bulk topolog-
ical invariant24.
Under the rotation of each spin vector about the z-axis
by the angle π/2, which is accomplished by the operator
Rz = ∏
N
j=1 ei π
2 Sz
j , H undergoes the transformation like
RzH(J, J,γ,hz)Rz = H(J, J, γ,hz).
(5)
It reveals that energy spectrum of H remains invariant
under sign reversal of γ. In other words, it implies that
as long as γ = 0, H lacks the U(1) symmetry. However,
H remains invariant finally under the transformation of
Vz = ∏
N
j=1 eiπSz
j , which means
VzH(J, J,γ,hz)V z = H(J, J,γ,hz).
(6)
As Vz performs rotation of the every spin vector about the
z-axis by the angle π, the Hamiltonian possesses the Z2
symmetry. Effect of these symmetry on the properties of
the system will be discussed in the proper context. Mag-
netic and topological properties of H will be described in
the subsequent sections. Antiferromagnetic (AFM) phase
will appear when J is assumed positive, and that will be
replaced by ferromagnetic (FM) phase when J is turned
negative. In the next section, several operators for identi-
fying magnetic orders employed in numerical or analytic
approaches will be discussed.
B. Long range correlations:
Existence of long range correlation of magnetic order
in a system can be studied by evaluating either the (i)
uniform correlation functions Cα
FM(n), α = x, y, z, for the
FM order or (ii) Neel (staggered spin-spin) correlation
functions Cα
Néel(n), for the AFM order. The expressions
of those functions are given as
Cα
FM(n) = 〈Oα
FM(n)〉, n = 0, 1, 2, ···(N −1), where,
Oα
FM(n) =
1
N
N
j=1
Sα
j Sα
j+n, and,
Cα
Néel(n) = 〈Oα
Néel(n)〉, where,
Oα
Néel(n) =
1
N
N
j=1
(−1)nSα
j Sα
j+n,
(7)
where Oα
FM/Néel(n) is the operator of FM/AFM order.
The number n denotes the separation between two spins
between which the correlation is to be measured. In order
to mark the phase transition points, Binder cumulant for
FM/AFM orders36,37:
Bα
FM/Néel = 1 −
〈(Mα
FM/Néel)4
3 〈(Mα
FM/Néel)22
,
(8)
can be evaluated numerically for the chains of finite
length, where
Mα
FM =
1
N
N
j=1
Sα
j , Mα
Néel =
1
N
N
j=1
(−1)jSα
j ,
(9)
are the operators for uniform and staggered magnetiza-
tions, respectively. In every case, expectation value 〈∗〉
has been evaluated with respect to the ground state. One
can further check that
(M
α
FM/Néel)2 = N2
N−1
n=0
Oα
FM/Néel(n).
(10)
Definition of Bα
FM/Néel (Eq. 8) indicates that the max-
imum possible value of it is 2/3. This maximum value
will be assumed when the system develops the perfect
LRO. Correlation functions for H obey the relations:
Cx
FM/Néel(n) = Cy
FM/Néel(n) = Cz
FM/Néel(n) as long as
γ = 0. But Cx
FM/Néel(n) = Cy
FM/Néel(n), when γ = 0, as
U(1) symmetry is being preserved by the system in this
case. Similar relation also holds for the Binder cumulants,
Bα
FM/Néel. However, in this study, Binder cumulant has
been evaluated to establish the existence of AFM LRO
only.
C. Limiting cases of H:
At several limits of the Hamiltonian, system is found
identical to well known spin models whose characteristics
are extensively studied long before. Here, properties of
such five different models are described briefly.

Page 4
4
1. For γ = ±1, J= 0, and hz = 0.
Hamiltonian reduces to Ising models, however the spin
operators align along the directions orthogonal to each
other, since
Hx = 2J
N
j=1
Sx
j Sx
j+1, and, Hy = 2J
N
j=1
Sy
j Sy
j+1,
respectively, when γ = ±1. Those Hamiltonians are con-
nected to each other via rotation of spin operators about
the z-axis by ±π/2. Apart from this those two Hamil-
tonians are equivalent to each other in every aspect for
obvious reasons like they possess LRO but exhibit no
QPT. The ground state is doubly degenerate for each
Hamiltonian in which Z2 symmetry is broken and the
system possesses a gap. At those points, Hamiltonians
retain their rotational symmetry about either x and y
directions, where AFM (FM) LRO is favoured along the
respective directions as long as J > 0 (J < 0). As a re-
sult, ground states correspond to a pair of Néel states in
terms of tensor product of eigenspinors of Sx and Sy op-
erators. So, in order to construct the exact ground state
wave functions at the degenerate point γ = 1, for the
thermodynamic limit, the normalized eigenspinors of Sx
operator are defined as Sx|χ±〉 = ±1
2 |χ±〉. Eigenspinors
can be further expressed as |χ±〉 = 1
√2 (|↑〉 ± |↓〉), where
Sz|↑〉 = +1
2 |↑〉, and Sz| ↓〉 = −1
2 | ↓〉. At this point of
the parameter space, Hamiltonian, Hx is found to com-
mute with the staggered spin-spin operator, Oα
Néel(n), for
α = x, but not for α = y and z. The exact form of dou-
bly degenerate normalized ground states at this point is
given by
±x〉 =
N/2
m=1,2,3,··· (|χ±2m−1 ⊗ |χ2m) ,
where total number of sites, N is even. Obviously, the
pair of ground states satisfies the orthonormality condi-
tion, 〈Ψµ
x ν
x〉 = δµν, where (µ, ν) = ±. Further, they are
connected to each other by lattice translation of unity.
Finally, it can be shown that ground state, |Ψ±x〉 satis-
fies the eigen value equations,
{ Hx±x〉 = −J
2 N±x,
Ox
Néel(n)|Ψ±x〉 = 1
4 ±x,
(11)
where the eigenvalues actually correspond to the exact
ground state energy per site, EG = −J
2
, and the exact
value of the x-component for staggered correlation func-
tions, Cx
Néel(n)=1/4, which is independent of n. At the
same time, they yield
Cα
Néel(n) = 〈Ψ±x|Oα
Néel(n)|Ψ±x〉 = 0,
when α = y and z. It means AFM LRO of the x-
component of spin-spin correlation exists and no LRO
persists for the y and z components of that. It is worth
mentioning at that point that no other types of correla-
tions like FM and chiral orders are found here. However,
AFM phase will be replaced by the FM one if J < 0,
with the same value of exact ground state energy per
site, EG = −
|J|
2 , and the exact value of x-component for
spin-spin correlation functions, Cx
FM(n)=1/4.
Similarly, at another degenerate point, say, γ = −1,
J= 0, and hz = 0, ground state of Hy, i. e., Ψ±y can
be constructed in terms of eigenspinors of Sy operator as
shown below:
±y〉 =
N/2
m=1,2,3,··· (|η±2m−1 ⊗ |η2m) ,
where Sy|η±〉 = ±1
2 |η±〉, and |η±〉 = 1
√2 (| ↑〉 ± i | ↓〉).
Obviously, |Ψ±y〉 satisfies the eigen value equations,
Hy±y〉 = −J
2 N±y,
Oy
Néel(n)|Ψ±y〉 = 1
4 ±y,
(12)
So, again EG = −J
2
, and the y-component for spin-spin
staggered correlation functions, Cy
Néel(n) = 1/4. At the
same time, they yield the relations,
Cα
Néel(n) = 〈Ψ±y|Oα
Néel(n)|Ψ±y〉 = 0,
when α = x and z for obvious reasons. So, both the
points γ = ±1, system exhibits AFM LRO at tempera-
ture T = 0, as long as J > 0. FM phase will replace the
AFM phase whenever J < 0.
Instead, interestingly enough, system at those points
(γ = ±1) hosts nontrivial topological characteristics39.
This feature can be shown when the single particle states
are studied in terms of spinless fermions by means of
JW fermionization. Under JW transformation, system
assumes nothing but the form of p-wave superconduc-
tor, which is known as the isotropic Kitaev chain with
zero chemical potential25. It is isotropic in a sense that
value of hopping parameter becomes equal to that of su-
perconducting parameter, and that equals to J. At this
parametrization it always exhibit the nontrivial topolog-
ical phase with ν = ±1, for γ = ±1, since the chemical
potential is absent.
2. For γ = ±1, J= 0, and hz = 0.
When hz = 0, Hamiltonian turns into the transverse
Ising model, where the system hosts LRO as long as
hz/J < 1, for T = 0. The system undergoes a phase
transition at hz/J = 1, to the disordered phase for
hz/J ≥ 1, in which ground state preserves the Z2 sym-
metry of the Hamiltonian13,9. Since hz plays the role of
chemical potential whenever the Hamiltonian expressed
in JW fermions, the resulting model becomes equal to
isotropic Kitaev chain with nonzero chemical potential.
As the values of hopping and superconducting parame-
ters are equal, system suffers a topological phase tran-
sition at the point hz/J = 1, separating the nontrivial
topological phase (ν = 1) for hz/J < 1, with the trivial
phase (ν = 0) for hz/J ≥ 1. So, the results clearly show
that topological and magnetic phases coexist for both the
cases, hz = 0, and hz = 0.

Page 5
5
3. For γ = 0, J= 0, and hz = 0.
At this point, Hamiltonian preserves the U(1) symme-
try. Correlation functions for hz = 0, behave as
Cβ
Néel(n) = {
(−1)
n+1
2
2
, if n ∈ odd,
0,
if n ∈ even,
where β = x, y. The system exhibits short range order
and so hosts no magnetic LRO, since
lim
n→∞
Cβ
Néel(n)=0.
Also there is no question for LRO for arbitrary hz. En-
ergy gap vanishes at this point and it separates two or-
dered phases around it. The system hosts two additional
multicritical points each one at |hz| = J3. Spin polarized
phase appears when |hz| ≥ J5. However, no topological
phase exists in this case for any value of hz.
4. For γ = 0, J= 0, and hz = 0.
Under this condition, the model is known as anisotropic
XY model in a transverse magnetic field. The resulting
spin model is exactly solvable, and there is an energy
gap. For γ > 0, the asymptotic value of the correlation
functions leads to7
lim
n→∞
Cx
Néel(n)=

1
2(1+γ) [γ2 {1−
(hz
J )2}]
1/4
, if |hz| < J,
0,
if |hz| ≥ J,
(13)
and
lim
n→∞
Cy
Néel(n)=0.
With the sign reversal of γ, correlation function, Cx
Néel(n),
will be replaced by Cy
Néel(n). The magnetic LRO ex-
ists within the range, −1 < hz/J < 1, irrespective of
the value of γ. System hosts the topological phase with
ν = +1 (ν = −1), for γ > 0 (γ < 0). Obviously, the
asymptotic behavior of correlations functions stated for
the previous cases 1, 2, 3 and 4 can be derived from the
more general expression in Eq. 13. Interplay of magnetic
correlation and entanglement in this case seems much in-
teresting by noting that von Neumann entropy has min-
imum on the circle16,
(hz
J )
2
+ γ2 = 1.
(14)
Remarkably, asymptotic behavior of the correlations,
Cβ
Néel(n → ∞), is found different around this circle. They
contain oscillatory terms within the circle, while outside
the circle they are monotonic7.
5. For γ = 0, and J= 0.
Disperson relation and magnetic properties of the sys-
tem under this condition have been studied before in
terms of spinless fermion38. The system exhibits no long-
range magnetic order in the absence of magnetic field,
however, it undergoes transition between two different
spin-liquid phases when J/J = 2. The ground state
phase diagram in the presence of both uniform and stag-
gered magnetic field is obtained, where spin-polarized FM
and AFM phases are found to appear, respectively, when
the strength of the field is very high. True magnetic
LRO has been developed in a system in order to minimize
the cooperative exchange energy, while the spin-polarized
phase appears when the Zeeman energy overcomes the
exchange energy. In this study FM version of the model
was considered since J < 038. Ground state phase dia-
gram with different types of spin-liquid phases have been
described. The system is no longer topological for any
values of J/J, and hz in this case.
Summarizing the results described in the last five cases
it indicates that no LRO exists when γ = 0, as shown in
the cases 3 and 5. Magnetic LRO exits when γ = 0 as
shown in the cases 1, 2, and 4, but at the same time
it indicates that AFM Cx
Néel survives for J > 0, when
γ > 0, while Cy
Néel survives for γ < 0. In the same way,
topological phases with ν = 1 exits when γ > 0, while
that with ν = −1 appears for γ < 0. It means LRO
with Cx
Néel (Cy
Néel) corresponds to topological phase with
ν =1(ν = −1). So one-to-one correspondence between
magnetic and topological phases has been shown. The
above argument holds but AFM correlations, Cx
Néel and
Cy
Néel would be replaced by FM correlations, Cx
FM and
Cy
FM, in the respective cases when J < 0. The same
correspondence is still valid in the presence of three-spin
terms which is being established in the present work.
6. For γ = 0, J= 0, and hz = 0.
The most general case has been studied extensively in
this work. It is still exactly solvable by means of spin-
less fermion. Both magnetic LRO and topological phases
are found present for this model and they coexist in the
parameter space. It means non-zero value of γ simul-
taneously induces both magnetic and topological orders.
Phase diagrams for both the orders have been made. The
effect of three-spin interacting term has indicated that
unlike the previous cases, magnetic field is not only able
to destroy the LRO but it is no more indispensable for
QPTs. As a result, parameter regime exhibiting LRO is
not bounded within a limited parameter space as long
as hz/Jis finite, since with the variation of hz, region
holding LRO is found to shift its location but without
altering its extend. Topological phases are detected as
usual by evaluating the Pfaffian invariant, winding num-
ber and symmetry protected zero-energy edge states. For
this purpose, the low energy single particle dispersion
relation is obtained in terms of JW fermionization. In
order to find magnetic LRO, spin-spin correlation func-
tions, Cβ
Néel(n) have been evaluated. Additionally, AFM
LRO has been detected numerically by finding the value
of them for the system with finite number of sites using
Lanczos exact diagonalization. Precise boundary of AFM
phase is found by evaluating the Binder cumulant in addi-

Page 6
6
tion. In the subsequent sections magnetic and topological
properties of this proposed model will be discussed and
it begins with the analytic results for four-spin plaquette
of H as described below.
III. MAGNETIC PROPERTIES OF H:
A. Four-spin plaquette for H:
The four-site Hamiltonian (N = 4) under the peri-
odic boundary condition (PBC) can be obtained analyt-
ically. The exact expression for all the 16 eigen values
(Em,m = 1, 2, ··· , 16) and eigen functions (ψm) of H
are shown in Appendix A. Energy spectrum consists of
two pairs of zero energy states. This feature is consistent
with the symmetry of the Hamiltonian, although they are
not related to topological edge modes by any means. ψ1
is the ground state as long as hz = 0 and γ = 1 for any
values of J. In addition, ψ5 is found degenerate with ψ1
only when J= 0. However, no ground state crossover is
there for hz = 0 and γ = 1, but ground state is doubly
degenerate when J= 0. This does not hold for γ = |1|.
Whereas, for hz = 0, ψ1 is the ground state in the region,
J′c < J< J′c, where
J′c = 2√(√η + √η2 − 8h2
zJ2 J)
2
J2γ2 − 2hz,
and η = J2(1 + γ2)+2h2
z. ψ5 is the ground state beyond
this region in this case. It means ground state crossover is
there, however, this occurrence cannot be related to the
phase transition anymore as it attributes to finite size
effect.
E
G
/
J
C
x N
e
e
l(1
)
(a)
(b)
γ =1
γ =1
hz =0
hz =0
J/J
0
1
2
−1
−3
−2
3
−0.5
−0.6
−0.7
0.15
0.20
0.25
FIG. 1: Variation of ground state energy per site, EG (a),
and correlation function, Cx
Néel(1) (b), for −3 ≤ J/J ≤ 3,
when γ = 1, and hz = 0. Red circles are the numerical data
while black dashed line for four-spin plaquette. Exact result,
Eq. 20, in (a) and Eq. 22, in (b) are plotted in cyan line for
comparison. Analytic and numerical results show an excellent
agreement.
The expressions of the correlation functions for four-
Cx
Néel(n)
n
(a)
γ = 1
J/J
hz =−2
0
0
246810
1214
−1
−3
−2
−4
−5
−6 −7
0.00
0.05
0.10
0.15
0.20
0.25
Cx
Néel(n)
n
(b)
γ = 1
J/J
hz =−1
0
0
1
2
2
4681012
14
−1
−3
−2
−4 −5
0.00
0.05
0.10
0.15
0.20
0.25
Cx
Néel(n)
n
(c)
γ = 1
J/J
hz =0
0
0
1
2
2
3
4
4
681012
14
−1
−3
−2
−4
0.00
0.05
0.10
0.15
0.20
0.25
Cx
Néel(n)
n
(d)
γ = 1
J/J
hz =1
0
0
1
2
2
3
4
4
5
681012
14
−1
−2
0.00
0.05
0.10
0.15
0.20
0.25
Cx
Néel(n)
n
(e)
γ = 1
J/J
hz =2
0
0
1
2
2
3
4
4
5
6
6
7
810
1214
0.00
0.05
0.10
0.15
0.20
0.25
FIG. 2: Correlation, Cx
Néel(n), when γ = 1, for hz = −2 (a),
hz = −1 (b), hz = 0 (c), hz = 1 (d), and hz = 2 (e).
spin plaquette are given by


Cz
Néel(n)=0,
Cx
Néel(n) = (ζ1+1)2
8(ζ2
1 +1)
,
Cy
Néel(n) = (ζ1−1)2
8(ζ2
1 +1)
,
(15)
with ζ1 =
γJ
hzJ+J/2−E1 . Those relations hold for n = 0,

Page 7
7
Cx
Néel(n)
n
(a)
γ = 2
J/J
hz =−2
0
0
2468101214
−1
−3
−2
−4
−5
−6
−7
0.00
0.05
0.10
0.15
0.20
0.25
Cx
Néel(n)
n
(b)
γ = 2
J/J
hz =−1
0
0
1
2
2
468101214
−1
−3
−2
−4
−5
0.00
0.05
0.10
0.15
0.20
0.25
Cx
Néel(n)
n
(c)
γ = 2
J/J
hz =0
0
0
1
2
2
3
4
4
6
810
12
14
−1
−3
−2
−4
0.00
0.05
0.10
0.15
0.20
0.25
Cx
Néel(n)
n
(d)
γ = 2
J/J
hz =1
0
0
1
2
2
3
4
4
5
68101214
−1
−2
0.00
0.05
0.10
0.15
0.20
0.25
Cx
Néel(n)
n
(e)
γ = 2
J/J
hz =2
0
0
1
2
2
3
4
4
5
6
6
7
8101214
0.00
0.05
0.10
0.15
0.20
0.25
FIG. 3: Correlation, Cx
Néel(n), when γ = 2, for hz = −2 (a),
hz = −1 (b), hz = 0 (c), hz = 1 (d), and hz = 2 (e).
since Cα
Néel(n = 0) = 1/4. Otherwise, they are in-
dependent of n, which is true only for N = 4. Be-
cause of the PBC, correlation functions satisfy the re-
lation, Cα
Néel(N n) = Cα
Néel(n). The result shows that
Cx
Néel(n) = Cy
Néel(n), when γ = 0. Variation of ground
state energy per site, EG, and Cx
Néel(n = 1) with respect
to J/J for this plaquette are shown by dashed line in
Fig. 1, when γ = 1 and hz = 0.
The results of four-spin plaquette has no importance
in general since they are drastically different from their
respective values for the system of thermodynamic limit,
N → ∞, specifically for the region where the strong quan-
tum fluctuation persists. But symmetries of the physical
quantities for the spin plaquette as obtained here remain
unaltered in the thermodynamic limit. Therefore, they
provide useful clue for obtaining results of larger system.
For example, Eq. 15 shows that Cz
Néel(n) = 0, for any
values of the parameters. It is found true in the asymp-
totic limit. Further, Cx
Néel(n) = Cy
Néel(n), when γ = 0,
and moreover, Cx
Néel(n, ±γ) = Cy
Néel(n, γ), for any val-
ues of Jand hz. It corresponds to the fact that values of
Cα
Néel(n), α = x, y, are interchangeable about the point,
γ = 0. Subsequent studies reveal that these properties
are independent of N by virtue of the symmetry. Finally
it leads to the fact that no LRO is there when γ = 0, irre-
spective of the values of Jand hz, as pointed out before
in cases 3 and 5.
In this context, EG, and Cx
Néel(n = 1) are compared
with the exact analytic and numerical results as displayed
in Fig. 1 (a) and (b), respectively. It reveals that value
of both EG, and Cx
Néel(1) for the plaquette is equal to the
respective exact results when γ = 1, but only for hz = 0,
and J= 0. It happens due to the fact at this point
(case 1), no quantum fluctuation persists. As a result,
both the values of EG, and Cx
Néel(1) attain their respective
maximum value. So, Cx
Néel(1) touches its saturated value,
i. e., Cx
Néel(1) = 1/4, at this point. It is another instance
where important results of the larger system could have
been captured in its four-spin replica.
B. Numerical results
In order to detect the existence of AFM LRO in the
ground state of H, x-component of staggered spin correla-
tions, Cx
Néel(n) have been computed. No LRO other than
Cx
Néel(n), (like FM and chiral phases) are found to appear
here. In order to investigate the properties of Cx
Néel(n),
ground states of spin chains with sites N = 20, 24, and
28 are obtained using the Lanczos exact diagonalization
techniques. As the Hamiltonian does not preserve the
U(1) symmetry, Hilbert space accommodates the states
of all Sz
T sectors. So the Hamiltonian matrix has been
spanned in the extended Hilbert space comprises with
Sz
T = ±N/2, ±(N − 1)/2, ··· , 0, sectors. Ultimately
Hilbert space is reduced manifold by taking into account
the translational symmetry of one lattice unit. Therefore
momentum wave vector k is introduced to associate this
symmetry by invoking PBC. The unique ground state
corresponds to either k = 0 or k = π, depending on the
values of N, and other parameters.
The variation of Cx
Néel(n) with respect to n, J/J,
and hz has been studied extensively, which is shown
in Fig 2 and 3, for γ = 1 and 2, respectively, when
N = 28. For each value of γ, Cx
Néel(n) is plotted for
hz = −2, −1, 0, 1, 2. In Figs 2 (c) and 3 (c), Cx
Néel(n)
are shown for γ = 1 and 2, when hz = 0. They indicate
that AFM LRO exists in the region, −2 < J/J < 2,
extended uniformly around its center J/J = 0. It con-
firms the fact that QPT occurs in this model at the

Page 8
8
points J/J = −2, +2, even in the absence of magnetic
field. The value of Cx
Néel(n) is also found symmetric about
the point, J/J = 0. In the regime, −2 < J/J < 2,
correlation function is very close to its saturated value,
Cx
Néel(n) ≈ 1/4, for both γ = 1 and 2. With the increase
(decrease) of hz, the region hosting LRO shifts toward
more positive (negative) side of J/J without changing
the width of the region. The diagrams 2 (d) and 3 (d),
indicate that LRO exists in the regime, 0 < J/J < 4,
when hz = 1. Similarly, diagrams 2 (e) and 3 (e), show
that LRO persists in the regime, 2 < J/J < 6, when
hz = 2.
On the other hand, AFM LRO is found to exist in the
regions, −4 < J/J < 0, and −6 < J/J < −2, when
hz = −1, and hz = −2, respectively, as depicted in two
pairs of diagrams, 2 (b), 3 (b), and 2 (a), 3 (a). The pro-
file of Cx
Néel(n) is absolutely symmetric around the center
of the region (J/J = 0) when hz = 0, whereas for hz = 0,
those shapes become more asymmetric around the re-
spective center of those region for nonzero LRO. An-
other interesting feature of those diagrams is that value
of Cx
Néel(n) decreases with the increase of |hz|. This is
due to the fact that contribution of Zeeman energy to
the ground state increases with the increase of |hz| at
the expense of exchange energy. Since the origin of mag-
netic LRO is associated with the cooperative exchange
energy, increase of |hz| leads to the lowering of the val-
ues of Cx
Néel(n), as the contribution of exchange energy is
becoming less. However, magnetic LRO is not at all di-
minished by the magnetic field as long as hz/Jis finite.
But the region of LRO is found to shift as an effect of the
field. Hence the region of nonzero magnetic LRO can be
identified by the relation, 2(hz J) < J< 2(J + hz), for
arbitrary hz. Cx
Néel(n) is found to diminish exponentially
with the separation n outside this region which confirms
the absence of LRO. Furthermore, in order to mark the
boundary of this region sharply, Binder cumulant for the
staggered correlation function has been evaluated numer-
ically as described below.
Evolution of Binder cumulant, Bx
Néel(N), (Eq. 8) with
respect to the number of sites, N helps to identify the
region where LRO exists. For this purpose, value of
Bx
Néel(N), for a system of finite number of spins, N =
20, 24, 28, have been estimated for different sets of pa-
rameters, hz = −2, −1, 0, 1, 2, and γ = 1 and 2. If there
exists LRO, Bx
Néel(N) grows with the increase of system
size, N, in contrast, Bx
Néel(N) decays with the increase
N, where short-ranged order persists36,37. As a result,
at the transition points, value of Bx
Néel(N) remains unal-
tered for any values of N, which are shown in Figs 4 for
γ = 1, and 5 for γ = 2, in great details.
According to the definition of Bx
Néel, (Eq. 8), the max-
imum possible value of Bx
Néel is 2/3, which is shown by
the horizontal dashed line in the Figs 4 and 5. It will
be observed when the corresponding magnetic order will
attain its saturated value. In this model, it happens for
hz = 0, and J/J = 0, as shown in Figs 4 and 5 (c),
where the peak of Bx
Néel touches the horizontal dashed
line (maximum value) for any value of N. It occurs due
to the fact that at this point, 〈Mx
Néel〉 = 1/2, which is
equal to its saturated value for spin-1/2 systems and cor-
responds to the maximum value of correlation function,
Cx
Néel(n)=1/4, for arbitrary n, as they are related by the
Eq. 10. Those results has been derived analytically be-
fore for γ = 1 in Eq. 11. As a result, 〈(Mx
Néel)4〉 = 1/16,
which finally leads to, Bx
Néel = 2/3.
In contrast, for hz = 0, value of Bx
Néel tends to increase
with the increase of N, indicating to touch its maximum
value for higher values of N, beyond N = 28. This phe-
nomenon indicates the existence of LRO in a specific re-
gion. Obviously, the reverse phenomenon confirms the
absence of LRO. However, value of the peak for Bx
Néel
decreases with the increase of |hz|, as clearly depicted in
the Figs. 4 and 5, (a, b, d, e). This phenomenon can
be explained by examining the relations among the rel-
evant quantities, Bx
Néel, Mx
Néel, and Cx
Néel(n) as stated in
the Eqs 8, 9 and 10, respectively. With the increase of
|hz|, contribution of Zeeman energy to the ground state
of the system increases. The competing Zeeman term re-
duces the effect of exchange term, which in turn reduces
the value of 〈Mx
Néel〉 below to its saturated value (1/2) in
such a fashion that Bx
Néel < 2/3.
Bx
Néel(N) for N = 20, 24 and 28 are drawn with red
dashed, green solid and blue dotted lines in each dia-
grams of Figs. 4 and 5. The shaded region bounded by
two vertical lines indicates the regime where LRO exists.
The region for nonzero LRO corresponds to the relation,
2(hz J) < J< 2(J + hz), for both γ = 1 and 2. Den-
sity plot for Bx
Néel when γ = 1 and 2 are shown in Fig.
6 (a) and (b), respectively. It indicates LRO is not af-
fected by the value of γ as long as it is nonzero and finite.
Analytic derivation of the correlation functions, Cx
Néel(n),
in terms of JW fermionization have been carried out in
the next section. The numerical results obtained in this
section shows an excellent agreement with the analytic
counterpart.
C. Jordan-Wigner fermionization
Several properties of the Hamiltonian, e. g., dispersion
relations, ground state energy, energy gap and spin-spin
correlations for H have been obtained analytically in this
section after expressing H in terms of spinless fermions.
Under JW transformations10:


S+
j = c
j
j−1
l=1
(1 − 2 nl),
S
j =
j−1
l=1
(1 − 2 nl) cj,
Sz
j = nj
1
2
,
(16)
where nl = c
l cl, is the number operator, the Hamiltonian
looks like
H =
1
2 ∑j [J (c
jcj+1 +c
j+1cj) −
J
2 (
c
jcj+2 +c
j+2cj)
+(c
jc
j+1 + cj+1cj) + 2hz (c
jcj
1
2)]
. (17)

Page 9
9
B
x N
e
e
l
(a)
γ =1
J/J
N=28
N=20
N=24
hz =−2
0
−1
−3
−2
−4
−5
−6
−7
0.1
0.2
0.3
0.4
0.5
0.6
0.7
B
x N
e
e
l
(b)
γ =1
J/J
N=28
N=20
N=24
hz =−1
0
1
2
−1
−3
−2
−4
−5
0.1
0.2
0.3
0.4
0.5
0.6
0.7
B
x N
e
e
l
(c)
γ =1
J/J
N=28
N=20
N=24
hz =0
0
1
2
−1
−3
−2
3
0.1
0.2
0.3
0.4
0.5
0.6
0.7
B
x N
e
e
l
(d)
γ =1
J/J
N=28
N=20
N=24
hz =1
0
1
2
−1
−2
3
4
5
0.1
0.2
0.3
0.4
0.5
0.6
0.7
B
x N
e
e
l
(e)
γ =1
J/J
N=28
N=20
N=24
hz =2
0
1
2
3
4
5
6
7
0.1
0.2
0.3
0.4
0.5
0.6
0.7
FIG. 4: Binder cumulant when γ =1 for hz = −2 (a), hz = −1
(b), hz = 0 (c), hz = 1 (d), and hz = 2 (e).
Comparing this to the form of Kitaev 1D model25, it
reveals that J and Jturns out to be equivalent to the
NN and next-nearest-neighbor (NNN) hopping integrals,
respectively, while γ and hz serve as the superconducting
and chemical potentials, respectively. Under the Fourier
transformation,
cj =
1
N
k∈BZ
ck eikj,
Hamiltonian acquires the Bogoliubov-de Gennes (BdG)
form, so,
H = ∑k [ψ
k H(k) ψk + ǫk
hz
2 ]
,
(18)
where ψ
k = [c
k ck], ǫk = J
2 (cos (k) − J
2J
cos (2k)) +
hz/2. Now introducing the vectors, g = (gx,gy,gz), and
B
x N
e
e
l
(a)
γ =2
J/J
N=28
N=20
N=24
hz =−2
0
−1
−3
−2
−4
−5
−6
−7
0.1
0.2
0.3
0.4
0.5
0.6
0.7
B
x N
e
e
l
(b)
γ =2
J/J
N=28
N=20
N=24
hz =−1
0
1
2
−1
−3
−2
−4
−5
0.1
0.2
0.3
0.4
0.5
0.6
0.7
B
x N
e
e
l
(c)
γ =2
J/J
N=28
N=20
N=24
hz =0
0
1
2
−1
−3
−2
3
0.1
0.2
0.3
0.4
0.5
0.6
0.7
B
x N
e
e
l
(d)
γ =2
J/J
N=28
N=20
N=24
hz =1
0
1
2
−1
−2
3
4
5
0.1
0.2
0.3
0.4
0.5
0.6
0.7
B
x N
e
e
l
(e)
γ =2
J/J
N=28
N=20
N=24
hz =2
0
1
2
3
4
5
6
7
0.1
0.2
0.3
0.4
0.5
0.6
0.7
FIG. 5: Binder cumulant, Bx
Néel, when γ = 2, for hz = −2 (a),
hz = −1 (b), hz = 0 (c), hz = 1 (d), and hz = 2 (e).
σ = (σxyz), where σα, α = x, y, z, are the Pauli
matrices, H(k) = g · σ. Now,
gx = 0,
gy = ∆k,
gz = ǫk,
where ∆k = 1
2 sin (k).
Under the Bogoliubov transformation, φk = Pkψk,
where Pk = (
uk vk
v
k
u
k ), and φ
k = [γ
k γk], the Hamil-
tonian becomes
H = ∑k [φ
k (P−1
k )
H(k) P−1
k φk + ǫk] ,
= NEG + ∑k
Ek [γ
k γk + γ
k γk] ,

Page 10
10
(a)
Bx
Néel
J/J
h
z
1.0
2.0
−2.0
−1.0
−1.5
−0.5
1.5
−2
−4
−6
2
4
0
6
0.1
0.2
0.3
0.4
0.0
0.0
0.5
0.5
0.6
0.7
(b)
Bx
Néel
J/J
h
z
1.0
2.0
−2.0
−1.0
−1.5
−0.5
1.5
−2
−4
−6
2
4
0
6
0.1
0.2
0.3
0.4
0.0
0.0
0.5
0.5
0.6
0.7
FIG. 6: Density plot of Binder cumulant, Bx
Néel, when γ = 1
(a), and γ = 2 (b), for −6 ≤ J/J ≤ 6.
when |uk| = √(1 + ǫk/Ek) /2, |vk| = √(1 − ǫk/Ek) /2,
where dispersion relation of the Bogoliubon (quasiparti-
cle) excitation is
Ek = √ǫ2
k + |∆k|2,
(19)
and the ground state energy per site is
EG = −
1
2π +π
π
Ekdk.
(20)
It indicates that superconducting phase survives as long
as γ = 0, where Cooper pairs are formed between parallel
spins. Value of ground state energy obtained numerically
is compared with the exact result, EG (Eq. 20) which is
shown in Fig. 1 (a), when γ = 1, and −3 ≤ J/J ≤ 3,
for hz = 0 . Eq. 20 is plotted in black line while red
circles are the numerical data. The excellent agreement
confirms the extreme accuracy of the numerical result.
However, for γ = 0, dispersion relation of spin excita-
tion is gapless as the Hamiltonian becomes
H(γ =0)=2 ∑k
ǫkc
kck,
since in this case |uk| = 1 and |vk| = 0. The ground state
energy per site for hz = 0, is estimated now by summing
the negative energy states38,
EG =
1
π k∈(ǫk<0)
ǫk dk.
(21)
Magnetic LRO, along with the topological superconduct-
ing phase cease to exist in this case. Dispersion relations
(Ek) and energy gap (EGap) of the system have been
evaluated for various cases as described below. EGap cor-
responds to the minimum value of Ek, which has been
obtained by solving the equation, dEk
dk
= 0, numerically,
for the fixed values of other parameters. Value of EGap is
important while determining the topological phase tran-
sition points. Because at the transition point EGap must
vanish.
D. Dispersion relations
For γ = 0, and hz = 0, energy dispersion is always
gapless as shown in Fig. 7. The positive and negative
portions of ǫk are drawn in red and blue, respectively.
It shows that, ǫk vanishes exactly at two distinct points,
k = ± arccos (J
J2 + 2J′2/2J), as long as J/J <
2. However for J/J > 2, ǫk vanishes at two additional
points marked by k = ± arccos (J + √J2 + 2J′2/2J)38.
In the presence of magnetic field, energy dispersion shifts
along the vertical direction depending on the sign of hz.
So the extent of positive and negative portions for ǫk shift
accordingly when hz = 0.
γ = 0
hz = 0
π
π
2
+π
2
+π
ǫk
/
J
J/J
k
0
0
0
1
1
2
2
−1
−1
−2
−2
−3
3
FIG. 7: Variation of ǫk/J with respect to J/J, within the BZ
for γ = 0, and −3 < J/J < 3, when hz = 0. Positive and
negative portions are drawn in red and blue, respectively.
In Fig 8 (a), (b) and (c), dispersion relations for hz =
−1, 0, and 1 are shown respectively, when γ = 1. They
are qualitatively different. For examples, when hz = 0,
three broad peaks and three valleys are there and they
are symmetric about the point J/J = 0. However, the
positions of valleys and peaks are interchanged around
J/J = 0. When J/J < 0, peaks appear at k = 0, ±π/2.
Height of the peaks increases with the increase of |J/J|.
Dispersion becomes gapless, Ek = 0, for J= −2J when
k = ±π, and for J= 2J when k = 0. There is always a
gap otherwise.
On the other hand, for hz = 1, there is one large peak at
k = 0, and two small peaks at k = ±π/2, when J/J < 0.
When J/J > 0, two large peaks appear at k = ±π/2,
while the deep valley appear at k = 0. Gap vanishes when
J/J = 0 and k = ±π. Even though the gap function is
nonzero since γ = 0, energy gap vanishes due to the effect
of Jand hz.
In order to identify the gapless region, minimum value
of E(k) has been estimated numerically. True band gap is

Page 11
11
(a)
γ =1
hz =−1
π
π
2
+π
2
+π
E
k
/
J
J/J
k
0
0
1
−1
−2
−3
−4
−5
0.0
1.0
0.5
1.5
(b)
γ =1
hz =0
π
π
2
+π
2
+π
E
k
/
J
J/J
k
0
0
1
2
−1
−2
−3
3
0.0
1.0
0.5
(c)
γ =1
hz =1
π
π
2
+π
2
+π
E
k
/
J
J/J
k
0
0
1
2
−1
3
4
5
0.0
1.0
0.5
1.5
FIG. 8: Dispersion relations for γ = 1, when hz = −1 (a),
hz = 0 (b), and hz = 1 (b).
equal to the twice of EGap. The variation of EGap in the
γ-J/J parameter space is shown in Fig. 9 (a), (b) and
(c), for hz = −1, 0, and 1. Obviously, EGap = 0, when
γ = 0 and hz = 0, since the superconducting phase does
not survive as shown in (b). It also serves as a line over
which phase transition occurs. Additionally EGap = 0,
along the lines J= ±2J, irrespective of the values of
γ. These lines serve as the boundaries of the trivial and
topological superconducting phases since the topological
phase persists in the annular region, −2 < J/J < 2,
which will be shown later.
For γ = 0, and hz = 0, there is a gap in the region,
−2hz/J < J/J < 0 when hz > 0 and in the region 0 <
J/J < 2hz/J, when hz < 0, as shown in 9 (c) and (a),
respectively, for hz = 1, and hz = −1. In contrast, when
γ = 0, there is a gap in the spectrum, except over two
parallel lines as described here. For example, band gap
vanishes at J= (−4J, 0), and J= (0, 4J), respectively,
when hz = −1, and hz = 1 as shown in (a) and (c), for any
value of γ. The results indicate that band gap vanishes
at the points, J= ±2J + 2hz, for arbitrary values of
hz. Further, it occurs exactly at k = π, and k = 0 for
(a)
E
G
a
p
/
J
γ
J/J
hz = −1
0
0
1
1
2
2
−1
−3
−2
−4
−5
−6
3
4
0.1
0.2
0.3
0.4
0.0
0.5
(b)
E
G
a
p
/
J
γ
J/J
hz = 0
0
0
1
1
2
2
−1
−3
−2
−4
3
3
4
4
0.1
0.2
0.3
0.4
0.0
0.5
(c)
E
G
a
p
/
J
γ
J/J
hz = 1
0
0
1
1
2
2
−1
−2
3
3
4
456
0.1
0.2
0.3
0.4
0.0
0.5
FIG. 9: Variation of band gap, EGap/J with respect to both
J/J and γ, when hz = −1 (a), hz = 0 (b), and hz = 1 (c).
J= ±2J +2hz, respectively. Otherwise, there is nonzero
band gap and the minimum of band gap appears when
k = 0, π. Although the value of EGap depends on γ,
location of those boundary lines is insensitive to the value
of γ.
E. Correlation functions
The analytic expression of spin-spin correlation func-
tion can be obtained easily by casting the spin operators
to spinless fermions. In terms of JW fermions, the cor-
relation function is given by the products of fermionic
operators4,6,7,
Cx
Néel(l m) =
1
4〈
BlAl+1Bl+1 ···Am−1
Bm−1
Am,
where Aj = c
j + cj, and Bj = c
j cj. The above expec-
tation value has been simplified by using Wick’s theorem
along with the fact that 〈AlAm〉 = 〈BlBm〉 = δlm. Fi-

Page 12
12
Cx
Néel(100)
(a)
γ =1
J/J
hz
0
0
1
2
2
4
6
−1
−2
−2
−4
−6
0.00
0.05
0.10
0.15
0.20
0.25
Cy
Néel(100)
(b)
γ =−1
J/J
hz
0
0
1
2
2
4
6
−1
−2
−2
−4
−6
0.00
0.05
0.10
0.15
0.20
0.25
FIG. 10: Correlation functions, Cx
Néel(n = 100) when γ = 1,
(a), and Cy
Néel(n = 100) when γ = −1, (b).
nally it assumes the from of Toeplitz determinant,
Cx
Néel(n) =
1
4
G−1
G−2
···
Gn
G0
G−1
··· Gn+1
...
...
...
...
Gn−2
Gn−3 ···
G−1
,
(22)
where the PBC is imposed and as a result, elements of
the determinant are nothing but the thermal expectation
value at temperature T, which is given by
Gn = 〈Bn+lAl,
=
1
π π
0
dktanh ( Ek
2kBT )
Ek
(ǫk cos (kn) − ∆k sin (kn)) .
Similarly, y- and z-component of correlation functions in
thermal equilibrium are given by
Cy
Néel(n) =
1
4
G1
G2
···
Gn
G0
G1
··· Gn−1
...
...
...
...
Gn+2 Gn+3 ···
G1
,
(23)
and
Cz
Néel(n) = 〈Sz2 GnGn/4,
(24)
where uniform magnetization along z-direction,
Sz〉 = −
1
π π
0
dk ǫktanh ( Ek
2kBT )
Ek
.
However, when T = 0, the ground state expectation val-
ues are given by
Sz〉 = −
1
π π
0
dk
ǫk
Ek
,
Gn =
1
π π
0
dk
ǫk cos (kn) − ∆k sin (kn)
Ek
.
Few limiting values can be derived easily. For exam-
ple, for γ = 0, Cx
Néel(n) = Cy
Néel(n). It happens due
to the fact that Gn = Gn, when γ = 0. In addition,
Cx
Néel(n, ±γ) = Cy
Néel(n, γ). When γ ≫ 1, Cx
Néel(n) = 0,
for any values of n, since the values of Gn along the
columns in the Eqs. 22, and 23 are becoming the same.
Which means Néel correlation does not survive when the
NN parallel spin cooper-pairing is very strong. Variation
of Correlation functions, Cx
Néel(n = 100) when γ = 1, and
Cy
Néel(n = 100) when γ = −1, are shown in Fig 10 (a)
and (b), respectively. They look identical reflecting the
fact that with the sign reversal of γ, correlations along x
and y directions are interchangeable. Additionally, cor-
relations preserve the symmetry, Cβ
Néel(n)(−J, hz) =
Cβ
Néel(n)(J,hz), which actually corresponds to the sym-
metry of Hamiltonian shown in Eq. 3. Cβ
Néel(100) has
the maximum value at hz = 0, which decreases with the
increase of |hz|. It means quantum fluctuations reduces
with the increase of |hz|, since the effect of competing
exchange integrals, J and Jis losing as a consequence.
IV. TOPOLOGICAL PROPERTIES OF H:
A. Pfaffian invariant and winding number
For γ = 0, H(k) satisfies the time-reversal, particle-
hole, and chiral symmetries as it obeys the following re-
lations,


PH(k)P−1 = −H(−k),
KH(k)K−1 = H(−k),
CH(k)C−1 = −H(k),
respectively. Here K, P = σxK, and C = σx are the
complex conjugation, particle-hole and chiral operators,
respectively. Conservation of any two symmetries corre-
sponds to that of the remaining one as they constitute
the BDI class. In order to characterize the topological
phase and location of the phase transition points, sign
of Pfaffian invariant and value of winding number have
been evaluated. Topological transition is marked by the
change of sign of the Pfaffian at the transition point. Fol-
lowing the general technique, Hamiltonian H(k) is being
converted to a skewsymmetric albeit hermitian matrix by
the transformation, ˜H(k) = DH(k)D, where
D =
1
√2[
1
i
1 −i ].
In order to locate the phase transition point, momentum-
space Pfaffian of the matrix i ˜H(k) is determined, which

Page 13
13
(a)
(b)
(c)
gz
gz
gz
gy
gy
gy
ν = 0
ν = 0
ν = 1
(d)
(e)
(f)
gz
gz
gz
gy
gy
gy
J< −2J
−2J<J< 2J
J> 2J
ν = 0
ν = 0
ν = 1
(g)
(h)
(i)
gz
gz
gz
gy
gy
gy
ν = 0
ν = 0
ν = 1
FIG. 11: Closed contours in the gz-gy plane for γ = 1. Figures
(a), (d) and (g) are plotted for J/J = −3, when hz = −1, 0, 1,
respectively. Similarly, (b), (e), (h) are for J/J = 0, and (c),
(f), (i) are for J/J = 3 when hz = −1, 0, 1, respectively.
In every case, direction of winding is counter clockwise. For
γ = −1, shape of the contours will be the same but with
winding of clockwise direction.
is defined as pf[i ˜H(k)] = √det{i ˜H(k)}40. As the EGap
vanishes exactly at k = π, and k = 0, for J= 2J + 2hz,
and J= −2J + 2hz, respectively, Pfaffians, pf[i ˜H(k =
π)], and pf[i ˜H(k = 0)] are evaluated. Those values are
given by
pf[i ˜H(π)] =
1
2 (
J
J
2
+ hz),
pf[i ˜H(0)] =
1
2 (
J
J
2
+ hz).
Hence the sign of pf[i ˜H(π)] and pf[i ˜H(0)] obey the rela-
tions:
sign(pf[i ˜H(π)]) = 

ve,
J> 2J + 2hz,
ve, −2J + 2hz < J< 2J + 2hz,
+ve,
J< −2J + 2hz,
sign(pf[i ˜H(0)]) = 

ve,
J> 2J + 2hz,
+ve, −2J + 2hz < J< 2J + 2hz,
+ve,
J< −2J + 2hz.
Now, the Pfaffian invariant, Q = sign(pf[i ˜H(π)]) ×
sign(pf[i ˜H(0)]), is given by40
Q = 

+ve,
J> 2J + 2hz,
ve, −2J + 2hz < J< 2J + 2hz,
+ve,
J< −2J + 2hz.
This result indicates that Q is negative in the annular
region enclosed by the boundary lines J= 2J +2hz, and
J= −2J + 2hz, while Q is positive elsewhere. In order
to characterize the topology, value of bulk topological
invariant, i. e., winding number (ν) has been determined,
which is defined as
ν =
1
2π C (ˆg(k) ×
d
dk
ˆg(k))dk,
where ˆg(k) = g(k)/|g(k)|, and C is a closed curve in the
gz-gy plane. Winding number enumerates the number of
winding around the origin, and at the same time, it will
be accounted as positive when the curve C is traversed
along the counter clockwise direction. For γ = 1, three
different sets of contours, {(a), (b), (c)}, {(d), (e), (f)},
and {(g), (h), (i)}, as shown in Fig 11 are drawn for
hz = −1, 0, 1, respectively. In each triplet set, three
different contours are consecutively drawn for J/J =
−3, 0, +3. However, in each diagram direction of the
contour is counter clockwise. So, the value of winding
number, ν = 1 for the diagrams (b), (e) and (h), since
in each case it encloses the origin once. Whereas, for the
remaining diagrams it is zero as they do not enclose the
origin. On the other hand, for γ = −1, the shapes of the
contours would be same but winding around clockwise
direction. Hence the nontrivial topological phase is define
by ν = −1, in this case. When γ > 0 (γ < 0), value of
winding number for H in the parameter space is given by
ν = {
1 (−1), −2J + 2hz < J< 2J + 2hz,
0,
J> 2J + 2hz and J< −2J + 2hz.
So, the points, γ = 0, and J= 2(hz ± J), are actually
the multicritical points.
In order to check the bulk-boundary correspondence
rule24, variation of eigen energies for H under open
boundary condition is shown in Fig. 12, with respect
to J/J, for hz = −1 (a), hz = 0 (b), and hz = 1
(c) when J = 1, and γ = 1. Pair of zero energy edge
states are present in the topologically nontrivial region,
2(hz J) < J< 2(J + hz). Finally a combined phase
diagrams for magnetic and topological phases is shown in

Page 14
14
-2.0
-1.5
-1.0
-0.5
0.0
0.5
1.0
1.5
2.0
-6
-5
-4
-3
-2
-1
0
1
2
E
n
erg
y
/
J
(a)
J/J
-1.5
-1.0
-0.5
0.0
0.5
1.0
1.5
-4
-3
-2
-1
0
1
2
3
4
E
n
erg
y
/
J
(b)
J/J
-2.0
-1.5
-1.0
-0.5
0.0
0.5
1.0
1.5
2.0
-2
-1
0
1
2
3
4
5
6
E
n
erg
y
/
J
(c)
J/J
FIG. 12: Variation of eigen energies with respect to J/J for
hz = −1 (a), hz = 0 (b), and hz = 1 (c) when J = 1, and
γ = 1. Zero energy edge states are present in the topological
region.
Fig. 13. Faithful coexistence of AFM LRO with nontriv-
ial topological order is noted in the Jhz space. The
parallel boundary lines enclosing this region are denoted
by the equations: J/2 = ±J + hz. Energy gap vanishes
over those line, as well as, the system undergoes simulta-
neous phase transition of magnetic and topological orders
on those lines.
B. Domino model
The spinless fermions have been converted to Majo-
rana fermions in order to study the difference of Majo-
rana pairings in trivial and topological phases. According
to the domino model a normal fermion at the j-th site is
made of two different Majoranas, Υa
j and Υb
j, as shown
in Fig. 14 (a), by blue and red dots, respectively4143.
Their positions over the chain is shown in Fig. 14 (b).
So, under the transformation,
Υa
j = cj + c
j, Υb
j = i (c
j cj),
(25)
-2
-1
0
1
2
-6
-4
-2
0
2
4
6
h
z
/
J
Antiferromagnet
Trivial
ν = ±1
ν = 0
J/J
FIG. 13: Magnetic and topological phase diagram in the
Jhz space. Topological phases of ν = ±1 is found to
coexist with AFM phase with ordering along x and y direc-
tions, respectively for J > 0, when γ > 0 and γ < 0. So, the
points, γ = 0, and J= 2(hz ±J), act like multicritical points.
AFM phase will be replaced by FM phase if J < 0.
Υa
j
Υb
j
(a)
(b)
(c)
(d)
(e)
j
j =1
j =1
j =1
j =1
2
2
2
2
3
3
3
3
4
4
4
4
N
N
N
N
···
···
···
···
N−1
N−1
N−1
N−1
FIG. 14: Domino model: (a) Spinless fermion at the j-th
site is composed of two Majoranas, Υa
j and Υb
j . (b) Lattice
composed of Majoranas. Majorana pairing in trivial phase
(c), and topological phases (d) and (e), for the extreme cases.
the Hamiltonian in Eq. 17, for hz = 0, is written as
H =
i
4
N−1
j=1 [J(1 − γa
j Υb
j+1 J(1 + γb
jΥa
j+1]
i
8
J
N−2
j=1 a
j Υb
j+2 + Υa
j+2Υb
j).
(26)
Three distinct phases can be understood in terms of three
different types of Majorana pairings for three extreme
cases. Among them trivial phase hosts one unique ma-
jorana pairing, while two distinct pairings are found in
the nontrivial phase. For example, the relations, J = 0,

Page 15
15
but J= 0, trivial corresponds to the trivial phase. Total
Hamiltonian now reads as
H = −i
J
8
N−2
j=1 a
j Υb
j+2 + Υa
j+2Υb
j).
In this case, an unique NNN intercell Majorana pairing
has been formed that is drawn by red and blue dotted
lines as shown in Fig. 14 (c). Obviously, all the Majo-
ranas participate in the pairing with no Majorana is left
unpaired. This picture is totally opposite to the trivial
phase found in Kitaev model, where intracell pairing has
been formed25. So unlike the Kitaev model, NNN Majo-
rana pairing is found here in the trivial phase.
Another extreme case leads to a pair of topological
phases which corresponds to the relations, J = 0, but
J= 0. Total Hamiltonian for γ = ±1 becomes
H = 
i J
2 N−1
j=1 Υb
jΥa
j+1, for γ = +1,
i J
2 N−1
j=1 Υa
j Υb
j+1,
for γ = −1.
In this case, two different NN intercell Majorana pairings
are permitted. The first (Υa
1) and the last (Υb
N ) Majo-
ranas are left unpaired leading to the appearance of zero-
energy edge states when γ = 1. This feature is identical
to the topological phase in Kitaev chain which is shown
in Fig. 14 (d)25. Additionally this picture corresponds
to the topological phase ν = 1. However, for γ = −1
the second (Υb
1) and the last but second (Υa
N ) Majoranas
are left unpaired which corresponds to the another phase
with ν = −1, as depicted in Fig. 14 (e). This particular
type of intercell Majorana pairing helps to construct new
fermionic quasiparticle basis in which the Hamiltonian,
H becomes diagonalized in the real space. In order to
accomplish this transformation, new fermionic operators,
¯cj are constructed by linear superposition of two adjacent
Majorana operators from NN sites42,
¯cj =
1
2 (
Υb
j + Υa
j+1) .
When γ = 1, in terms of the new fermionic operators,
H = −J N−1
j=1 ¯c
j ¯cj, indicating the cost of creating a new
fermion at any site j, is J in the topological phase. No
such fermionic basis can be constructed for the trivial
phase in which H is diagonalized, since two different types
of pairing are superposed in this case.
V. DISCUSSION
In this work we have introduced a spin-1/2 1D
anisotropic extended XY model in order to study the in-
terplay of magnetic and topological phases. In this sys-
tem a three-spin term has been added to the anisotropic
XY model which could be solved in terms of JW fermions.
Existence of long range magnetic correlations has been
confirmed numerically. Both magnetic and topological
properties have been studied extensively. System exhibits
AFM LRO along two orthogonal directions, say x and y
directions at the same location in the parameter space, if
the anisotropic parameter γ picks up +ve and −ve signs,
respectively. Topological superconducting phases have
been characterized by evaluating the Pfaffian invariant,
winding number, and zero-energy edge states. System
hosts topological phases with ν = ±1 in the same way
such that the coexistence of magnetic and topological
phases in the parameter space is found along with the
associated one-to-one correspondence between magnetic
and topological phases is obtained.
Magnetic and topological properties of this model at
various extreme limits are investigated. It shows that the
anisotropic parameter γ plays the crucial role behind the
origin of magnetic LRO as well as the nontrivial topologi-
cal phases. Finally, the trivial and topological phases are
explained in terms of different types of Majorana pairings.
Two different NN pairings have appeared in the topolog-
ical phases for the two extreme cases defined by J= 0,
but γ = ±1. On the other hand, an unique NNN Majo-
rana pairing is found in the other extreme case, J = 0,
which corresponds to the trivial phase. In contrast, only
NN Majorana pairing is found in the Kitaev model for
both trivial and topological phases25.
Those phases as well their one-to-one correspondence
could be destroyed by applying transverse magnetic field,
beyond hz/J ≥ 1, in the absence of three-spin term.
However, in the presence of three-spin term, coexistence
of those phases as well their one-to-one correspondence
would never be broken or eliminated by applying the mag-
netic field of finite strength. So, in contrast, this model
exhibits magnetic LRO in the magnetic field of any value.
Location of those phases in the parameter spaces is given
by 2(hz J) < J< 2(J + hz). Position of multicritical
points are given by γ = 0, and J= 2(hz ± J). On the
other hand, QPT could occur even in the absence of the
magnetic field. It occurs due to the fact that magnetic
field always tend to destroy the LRO resourced by the ex-
change integrals J and J. However, in this model J and
Jalso compete to each other. So, QPT might occur due
to the competing effect of J and Jalone. In this sense,
the three-spin terms induces exotic magnetic phases those
are not explored before. Magnetic field always oppose the
LRO in systems of any dimensions68,4446. On the other
hand, there are several instances where magnetic field
is indispensable for the emergence of topological phases
specially in the magnetic systems4751. However, in this
work, effect of magnetic field both on the magnetic and
topological phases are investigated and their interplay has
been explored.
VI. ACKNOWLEDGMENTS
RKM acknowledges the DST/INSPIRE Fellow-
ship/2019/IF190085.
Appendix A: ENERGY EIGENVALUES AND
EIGENSTATES OF THE FOUR-SITE
HAMILTONIAN, H
Expressions of 16 eigenvalues of the four-site Hamilto-
nian, Hare as follows.

Page 16
16
E1 = −J − √J2γ2 + (hz + J/2)2,
E2 = −J + √J2γ2 + (hz + J/2)2,
E3 = J − √J2γ2 + (hz + J/2)2,
E4 = J + √J2γ2 + (hz + J/2)2,
E5 = −√J2(1+γ2)+2h2
z +√(J2(1+γ2)+2h2
z)2 − 8J2h2
z,
E6 = √J2(1+γ2)+2h2
z +√(J2(1+γ2)+2h2
z)2 − 8J2h2
z,
E7 = −√J2(1+γ2)+2h2
z −√(J2(1+γ2)+2h2
z)2 − 8J2h2
z,
E8 = √J2(1+γ2)+2h2
z −√(J2(1+γ2)+2h2
z)2 − 8J2h2
z,
E9 = E11 = E12 = E13 = 0,
E13 = E14 = −hz + J/2,
E15 = E16 = hz J/2.
To express the eigenfunctions in a compact form, fol-
lowing notations are used.
χ2
n = Tn−1 |2〉 (n = 1),|2〉 = |↑↑↑↑〉,
χ1
n = Tn−1 |1〉 (n = 1, 2, 3, 4),|1〉 = |↓↑↑↑〉,
χ0
n,1 = Tn−1 |0〉(n = 1, 2, 3, 4),|0〉 = |↑↑↓↓〉,
χ0
n,2 = Tn−1 |0〉(n = 1, 2),|0〉 = |↑↓↑↓〉,
χ−1
n
= Tn−1 |−1〉 (n = 1, 2, 3, 4),|−1〉 = |↑↓↓↓〉,
χ−2
n
= Tn−1 |−2〉 (n = 1),|−2〉 = |↓↓↓↓〉.
Here the operator T behaves as T |pqrs〉 = |spqr〉. The
normalized eigenstates now can be expressed as follows.
ψi =
1
2ζ2
i +1 4
n=1(−1)n(χ−1
n
+ ζiχ1
n), i = 1, 2,
ψk =
1
2ζ2
k+1 4
n=1(−χ−1
n + ζkχ1
n), k = 3, 4,
ψl = 1
2Nl
[(E2
l − 4h2
z)(El 4
n=1 χ0
n,1 + 2J 2
n=1 χ0
n,2)
+2γJEl{(El − 2hz)χ−2
1
+ (El + 2hz)χ2
1}], l = 5, 6, 7, 8,
ψm = 1
√2 (χ0
n χ0
n+1), n = 1, 2, 3 and m = 9, 10, 11,
ψ12 = 1
√2 2
n=1(−1)n(n−1)/2χ0
n,2,
ψ13 = 1
2 4
n=1(−1)n(n−1)/2χ−1
n ,
ψ14 = 1
2 4
n=1(−1)n(n+1)/2−1χ−1
n ,
ψ15 = 1
2 4
n=1(−1)n(n+1)/2χ1
n,
ψ16 = 1
2 4
n=1(−1)n(n−1)/2χ1
n,
with ζi =
γJ
hz J+J/2−Ei
,
ζk =
γJ
hz+J+J/2−Ek
and
Nl = √(
√2γJEl)2 + (E2
l + 2J2)(E2
l − 4h2
z)2 + (2√2γJElhz)2.
Electronic address: rkmalakar75@gmail.com
Electronic address: asimk.ghosh@jadavpuruniversity.in
1 S Sachdev, Quantum Phase Transitions, Cambridge Uni-
versity Press, Cambridge, (1999).
2 B K Chakrabarti, A Dutta, P Sen, Quantum Ising Phases
and Transitions in Transverse Ising Models, Springer,
Berlin (1995).
3 A Dutta et. al., Quantum Phase Transitions in Transverse
Field Spin Models, Cambridge University Press, Cam-
bridge, (2015)
4 E. Lieb, T. Schultz, and D. Mattis, Ann. Phys. 16, 407
(1961)
5 S. Katsura, Phys. Rev. 127, 1508 (1962).
6 E. Barouch, B. M. McCoy and M Dresden, Phys. Rev. A
2, 1075 (1970).
7 E. Barouch and B. M. McCoy, Phys. Rev. A 3, 786 (1971).
8 E. Barouch and B. M. McCoy, Phys. Rev. A 3, 2137 (1971).
9 P. Pfeuty, Ann. Phys. 57, 79 (1970).
10 P. Jordan, E. Wigner, Z. Phys. 47, 631 (1928)
11 R Coldea et. al., Science, 327, 177 (2010).
12 R. Blinc, J. Phys. Chem. Solids 13, 204 (1960).
13 R B Stinchcombe, J . Phys. C: Solid State Phys., 6, 2459
(1973).
14 Osterloh A, Amico L, Falci G and Fazio R, Nature 416 608
(2002)
15 T. J. Osborne, and M. A. Nielsen, Phys. Rev. A 66, 032110
(2002).
16 F Franchini, A R Its, B-Q Jin and V E Korepin, J. Phys.
A: Math. Theor. 40 8467 (2007)
17 Nielsen M and Chuang I, Quantum Computation and
Quantum Communication (Cambridge Univ. Press, Cam-
bridge, 2000).
18 Klitzing K. V., Rev. Mod. Phys. 58, 519 (1986).
19 Thouless D. J., Kohomoto M., Nightingale P. and den Nijs
M., Phys. Rev. Lett. 49, 405 (1982).
20 Haldane F. D. M., Phys. Rev. Lett. 61, 2015 (1988).
21 W. Su, J. Schrieffer and A. J. Heeger, Phys. Rev. Lett. 42,
1698 (1979).
22 Heeger A. J., Kivelson S., Schrieffer J. R. and Su W. -P.,
Rev. Mod. Phys. 60, 781 (1988).
23 R K Malakar and A K Ghosh, J . Phys. Condens. Matter,
35, 335401 (2023).
24 Y. Hatsugai, Phys. Rev. Lett. 71, 3697 (1993).
25 A. Y. Kitaev, Phys. -Usp. 44, 131, (2001)
26 Y. -X. Chen and S.-W. Li, Phys. Rev. A 81, 032120 (2010).
27 J. D. Sau, R. M. Lutchyn, S. Tewari, S. Das Sarma, Phys.
Rev. Lett. 104, 040502 (2010).
28 L. N. Bulaevskii, A. I. Buzdin and S. S. Crotov, Solid State
Commun. 48, 719 (1983)
29 Bulaevskii L. N., Buzdin A. I., Kulic M. L. and Panjukov
S. V., Adv. Phys. 34, 175 (1985).
30 S. S. Saxena et. al., Phys. Rev. B 63, 144519 (2001).
31 C. Pfleiderer, M. Uhlarz, S. M. Hayden, R. Vollmer, H. von
Lohneysen, N.R. Bernhoeft, and G.G. Lonzarich, Nature
412, 58 (2001).
32 D. Aoki, A. Huxley, E. Ressouche, D. Braithwaite, J. Flo-
quet, J-P. Brison, E. Lhotel, and C. Paulsen, Nature 413,
613 (2001).
33 T. R. Kirkpatrick and D. Belitz, Phys. Rev. B. 67, 024515
(2003)
34 J. A. Bert et. al., Nature Phys. 7, 767 (2011).
35 Yu et. al., npj Quantum Mater. 6, 63 (2021).
36 K. Binder, Phys. Rev. Lett. 47, 693 (1981).
37 K. Binder, Z. Phys. B 43, 119 (1981).
38 I. Titvinidze and G. I. Japaridze, Eur. Phys. J. B 32, 383-
393 (2003)
39 K Chhajed, Resonance 26 1539 (2021)
40 J. D. Sau and S Tewari, Phys. Rev. B 88, 054503 (2013).
41 C.W.J. Beenakker, Annu. Rev. Condens. Matter Phys. 4,
113 (2013).
42 J Alicea, Rep. Prog. Phys. 75, 076501 (2012).
43 M Leijnse and K Flensberg, Semicond. Sci. Technol. 27,
124003 (2012).
44 H -J Mikeska, A Ghosh and A K Kolezhuk, Phys. Rev.
Lett. 93 217204 (2004).
45 A Ghosh, J. of Phys: Condens. Matter 13 5205 (2001).
46 A K Ghosh, Phys. Rev. B 80 214418 (2009).
47 Owerre S A, J. Appl. Phys. 120, 043903 (2016).

Page 17
17
48 D G Joshi, Phys. Rev. B 98, 060405(R) (2018).
49 M Deb and A K Ghosh, J. Phys.: Condens. Matter 31,
345601 (2019).
50 Sil A. and Ghosh A. K., J. Phys.: Condens. Matter 32,
205601 (2020).
51 M Deb and A K Ghosh, J. Magn. Magn. Mater. 533,
167968 (2021).