# Deciphering competing interactions of Kitaev-Heisenberg-$\Gamma$ system in
clusters

Sheikh Moonsun Pervez ID<EMAIL_ADDRESS>Institute of Physics, Sachivalaya
Marg, Bhubaneswar-751005, India  Homi Bhabha National Institute, Training
School Complex, Anushakti Nagar, Mumbai 400094, India  Saptarshi Mandal ID
<EMAIL_ADDRESS>Institute of Physics, Sachivalaya Marg,
Bhubaneswar-751005, India  Homi Bhabha National Institute, Training School
Complex, Anushakti Nagar, Mumbai 400094, India

###### Abstract

We study Kitaev-Heisenberg-$\Gamma$ Hamiltonian on finite size clusters and
investigate various zero and finite temperature properties in detail and
explore the key differences due to relative sign and magnitude of Kitaev($K$),
Heisenberg $(J)$) and $\Gamma$ interaction in the presence of magnetic field.
Nearest-neighbor spin-spin correlation function shows step like character when
$J$ is positive only irrespective of sign of $K$. Similarly the low energy
spectrum under magnetic field depends on the sign of $K$. Dynamics of Majorana
fermions depend greatly on external magnetic field which can even stabilize
the gauge fields to Kitaev limit in the presence of $J$. Majorana dynamics
also shows that $K-J-\Gamma$ model to be a testbed to manipulate quantum speed
limit. For small magnetic field, magnetization shows a dome like structure for
an intermediate temperature range for positive $K$. However a positive $J$ can
cause this for negative $K$ also. The zero temperature magnetisation and
susceptibility shows large dependence on external magnetic field and relative
sign and magnitude of $K$ and $J$. Either a negative $K$ or $J$ is shown to
give qualitatively similar susceptibility for different magnetic field only.
For the specific heat a negative $K$ or $J$ favors mostly single peak
structure while a two peak structure is favored for positive $K$ or $J$,
though the magnitude of magnetic field can causes re-entrant behavior also. We
briefly discuss effect of $\Gamma$ and compare with recent experiments.

##  I Introduction

Starting from its inception Kitaev modelkitaev-2006  has witnessed a
remarkable journey which is still continuing with overwhelming excitement.
Though initially it was proposed as a mere toy model with certain aspect of
quantum computation in perspective, it has now traversed almost all the aspect
a quantum many body system can offer. To mention a few of them, exact short
range spin-spin correlation smandal-2007  established its spin liquid nature,
effect of disorder and localization kao-2021 ; nasu-2020 has been investigated
for such spin liquid. Realization of Kondo effects and its competition with
Kitaev physics has been investigated tathagata-2020  and spin-spin correlation
function for perturbed Kitaev model tikhonov-2011 ; lunkin-2019  has found a
power low decaying spin-spin correlation function. The dynamics of perturbed
Kitaev model have explored several aspects such as confinement-deconfinement
transition, quantum phase transitions in perturbed Kitaev model saptarshi-2011
; animesh-2020 ; animesh-2021 ; knolle-2018 . Kitaev model for classical spins
baskaran-2008 ; samarth-2010  and in variety of lattices in two and three
dimensions kivelson-2007 ; vala-csl-2010 ; naveen-2008 ; eschmann-2020  showed
interesting feature such as order-disorder mechanism, existence of chiral spin
liquid, gapless fermi contour etc.
Apart from the above theoretically motivated works, akin to a dream realized
in a mysterious way, very soon a plethora of materials have been proposed to
posses Kitaev like interaction along with other kind of interactions
chalaupka-2010 ; simon-2022 ; takagi-2019 ; motome-2020 ; hermanns-2018 ;
abanerjee-2016 . The salient feature of these materials are that the presence
of non-Kitaev interactions make the coveted spin-liquid state non-existent at
very low temperature. However at higher temperature the Kitaev physics mainly
the fractionalization of spins into Majorana fermions and fluxes and the
intrinsic differences in the length and time scale of their dynamics enable
the spin liquid state to be revived at intermediate temperaturejanssen-2017 ;
chern-2021 ; czajka-2021 ; niravkumar-2019 ; wulferding-2020 ; berke-prb-2020
; haoxiang-nature-2021 ; balz-prb-2021 ; plessis-2020 . Further recent studies
show that proximity effect of $\alpha-RuCl_{3}$ in graphene heterostructure
sananda-2019 ; sananda-2021  increases the strength of Kitaev interaction and
moreover induces an insulator to metal transition which might lead to exotic
superconducting states. Such heterostructure study holds a great promise for
interesting Kitaev related physics to be observed in near future.

One of the fascinating aspect of condensed matter system is that physical
interactions at microscopic ranges govern the physical properties at
macroscopic length anderson ; fifty-years . Commonly all interactions are
extended over few atoms or molecules governed by mainly electrostatic
interactions combined with quantum mechanical laws and inherent symmetries. It
is also observed how the physical properties of the system changes as the size
of the system varies from microscopic to macroscopic scale. One such example
is quantum spin system. Various phenomena such as long range order, symmetry
breaking, fractionalization of spins, anyonic excitations, certain properties
of spin-spin correlation functions, topological properties are only realized
over large enough system sizessavary-review . However in this regards, there
are some exception in Kitaev Model where one of the defining property of
Kitaev system i.e bond-dependent short range correlation function and spin-
fractionalization is realized exactly even for the smallest size defined
properly. This renders some aspect of Kitaev magnetism or exactly for Kitaev
model in small clusters tathagata-2020 ; kataoka-jpsj-2020 . Further the
primary aspect of frustration or absence of it is mainly governed by local
considerations of exchange coupling of spin-spin interaction though globally
they drive the system completely differently. However it is an interesting
question to itself how these local dynamics of commonly prevalent local
exchange interaction competes with Kitaev type exchange interactions and
whether there is a definite pattern in it wang-2020 . To the best our
knowledge such detail exact comparative study is absent within a single work.
This motivates us to study Kitaev model in small clusters to explore several
aspects in detail. Further due to inherent difficulties to deal with an
interacting many body system it is not possible to to analysis in a single
study various competing mechanism exist in Kitaev systems. For example the
relative sign and strength of Kitaev interaction and Heisenberg interaction
and $\Gamma$ interaction would make the analysis very demanding for a single
study. To this end we consider exact analysis of Kitaev-Heisenberg-$\Gamma$
systems in small clusters and investigate some very essential quantities such
as eigenspectrum, nearest-neighbor spin-spin correlation function, flux
operators, detail dynamics of gauge and matter Majorana fermions, specific
heat, susceptibility and specific heat for all possible combinations of Kitaev
interaction, Heisenberg interaction, $\Gamma$ and external magnetic field in
$\hat{z}$ direction. We hope that our study would serve as a valuable
reference for Kitaev-Heisenberg-$\Gamma$ system.

Figure 1: Schematics of Kitaev clusters with (a) 4-site, (b) 6-site, (c)
8-site, and, (d) 12-site are shown in panel a, b, c and d respectively. The
red bonds, blue bonds and green dashed bonds are having $x-x$-type, $y-y$ type
and $z-z$ type Kitaev interaction respectively.

Our plan of presentation is as follows. We begin by investigating how the
spectrum of Kitaev-Heisenberg-$\Gamma$ system evolves as we vary the strength
of Heisenberg interaction for various values of $\Gamma$ for different cluster
of sizes in Sec. II. We also considered the effect of external magnetic field
on the low energy eigenvalues. Our analysis shows that the eigenenergies of
the Kitaev-Heisenberg-$\Gamma$ system van be grouped into two branches. In the
first branch, energies vary linearly for all ranges of Heisenberg coupling $J$
and for the second branch energies vary linearly only for large magnitude of
$J$ but for small values they shows a cusp or dip near $J=0$ whose exact
position depends on the value of $\Gamma$ and cluster sizes. Interestingly the
linear branches changes their energies from positive to negative when the sign
of $J$ is reversed but for second type it remains negative or positive as we
change the sign of $J$. The incidental crossing between the two kind of
energies can be controlled by tuning $\Gamma$.

Figure 2: Eigenvalues are plotted for various size of clusters as shown. The
eigenvalues are renormalized by dividing with $6\times N$ where $N$ is the
number of sites in a given cluster. First, second, third and fourth columns
are for four site, six site, eight site and twelve site clusters as denoted.
The first three rows are for antiferromagnetic Kitaev interaction with $K=1$
and the last three rows are for ferromagnetic Kitaev interaction in Kitaev
clusters for varying $J$ with fixed $K=-1$. The values of $\Gamma$ taken as
shown inside each panel.

After investigating the spectrum, we focus on the nearest-neighbor correlation
function in Sec. III and study how it depends on cluster sizes and on
different system parameters as before. We have considered both
$\sigma^{z}_{i}\sigma^{z}_{j}$ and $\sigma^{y}_{i}\sigma^{y}_{j}$. We find
that correlation function shows interesting variations depending on the
relative sign and magnitude of Kitaev interaction ($K$), Heisenberg
interaction $J$ and $\Gamma$. Very interestingly we find a clear distinction
between the ferromagnetic and antiferromagnetic Heisenberg interaction on the
nearest neighbor correlation function. For the antiferromagnetic Heisenberg
interaction the correlation function shows plateau like behavior where spin-
spin correlation does not change as we change $J$. This suggest a competing
mechanism between the quantum entanglement induced by Kitaev interaction and
antiferromagnetic Heisenberg interaction. On the other hand for ferromagnetic
Heisenberg interaction such plateaus are absent and shows various non-
monotonous behavior showing a different competition between the quantum
fluctuation by Kitaev interaction in response to ferromagnetic Heisenberg
interaction which tries to induce a collinear order in the system.

Next we study the fate of conserved plaquette operator in Sec. IV and dynamics
of Majorana fermions by exact numerical computation in Sec. V. We find that
the amplitude and time period of oscillations of Majorana fermions greatly
depend on the relative sign of the Kitaev interaction ($K$) and Heisenberg
interaction $(J)$. On an average a ferromagnetic Kitaev interaction is seen to
increase the time period of oscillations. We noticed that the external
magnetic field can modulate the nature of the dynamics of the Majorana fermion
greatly and depending on the directions of it, the external magnetic field can
stabilize the fluctuation of conserved gauge fermions. In certain cases it
shows when the probability of Majorana fermion is zero describing a scenario
to define quantum speed limit of Majorana fermions in these system.

Finally we present finite temperature magnetization, susceptibility and
specific heat in Sec. VI, VII. VIII respectively. We give a detail description
how these quantities depend on various magnetic exchange coupling
($K,J,\Gamma$) including external magnetic field $h_{z}$. The primary findings
on magnetization is that for antiferromagnetic Kitaev interaction and low
magnetic field, there is an intermediate range of temperature where
magnetization posses a smooth dome like behavior which is absent for
ferromagnetic Kitaev interaction. However depending on the Heisenberg
interaction this dome like structure can be re-entrant for ferromagnetic
Kitaev interaction. The external magnetic field is shown to control the zero
temperature magnetization which also depends differently for different
combinations of $J,K$ and $\Gamma$. Similar to magnetization, we find that
ferromagnetic and antiferromagnetic Kitaev interaction has distinct features
which again can change differently depending on the sign and magnitude of
Heisenberg interaction and applied external magnetic field.

Figure 3: Few lowest eigenvalues are plotted with external magnetic field
applied in $z$-direction. Each column denotes a particular cluster as
mentioned at the top. Each row denotes a combination of Kitaev ($K$) coupling
and Heisenberg coupling (J) as mentioned at the right side of each row. The
ground state energy ($E_{0}$) is shown by green dots and first excited state
$E_{1}$ is shown by blue bots and few more excited eigenstates $E_{n}$ with
$n>1$ are shown by red dots.

##  II Kitaev-cluster-Basics and spectrum analysis

The defining Hamiltonian we investigate in our study is $J-K-\Gamma$
modelwang-2021  which is given as,

$\displaystyle H=\sum_{\langle i,j\rangle}\left(J{\bf S}_{i}\cdot{\bf
S}_{j}+KS^{\gamma}_{i}S^{\gamma}_{j}+\Gamma(S^{\alpha}_{i}S^{\beta}_{j}+S^{\beta}_{i}S^{\alpha}_{j})\right)$
(1)

We begin our study of various sizes of Kitaev cluster by first analyzing the
eigenvalue spectrum and understand qualitatively the effect of Heisenberg and
$\Gamma$ interaction as plotted in Fig. 2. We start with four site cluster and
discuss $\Gamma=0$ at first. We notice that eigenvalues could be grouped into
two categories. Firstly there are eigenstates whose energies are linearly
varying with $J$ and secondly there are eigenstates whose energies vary non-
linearly for small values of $J$ but eventually vary linearly for large values
of $J$. This is true for all clusters. We find that a finite $\Gamma$ affects
the eigenspectra in two ways. Firstly it breaks the degenerates of eigenvalues
found for $\Gamma=0$. It is interesting to note that only those eigenvalues
which vary linearly become non-degenerate but the resultant eigenenergies
again vary linearly with $J$. Secondly the band width (the separation between
the highest and lowest eigenvalues) increases for non-zero $\Gamma$ and higher
for positive $\Gamma$. Also it is evident that $\Gamma$ determines the gap at
$J=0$. The physical reason why $\Gamma$ interaction breaks the degeneracy of
the eigenvalues is that it breaks the reflection symmetry under
$\sigma_{x}\rightarrow-\sigma_{x}$ but $\sigma_{y,z}\rightarrow\sigma_{y,z}$.

For large positive $J$, the four site cluster has large
singlet($\frac{1}{2}\ket{\uparrow_{1}\downarrow_{2}-\downarrow_{1}\uparrow_{2}}\times\ket{\uparrow_{3}\downarrow_{4}-\downarrow_{3}\uparrow_{4}}$)
contribution to the ground state. The other states which are varying linearly
can be approximated by classical state such as
$\left(\uparrow_{1}\downarrow_{2}\uparrow_{3}\downarrow_{4}\right)$ having
energy $E_{cl}=-2J+2K$ and
$\left(\uparrow_{1}\downarrow_{2}\downarrow_{3}\uparrow_{4}\right)$
($\left(\uparrow_{1}\uparrow_{2}\downarrow_{3}\downarrow_{4}\right)$) having
energy $E_{cl}=-2J-2K$. The highest (lowest) eigenvalues for positive
(negative) $J$ can be approximated as $E_{cl}=6J+2K$ ($E_{cl}=-6J+2K$) which
is determined by an approximated ferromagnetic state though quantum
fluctuation is more effective for antiferromagnetic $J$ to increase the
magnitude of the eigenvalue further. It is evident that as we decrease the
value of $J$ from a very large positive value to a large negative value, the
quantum ground state is continually adjusted from a large singlet component to
a large ferromagnetic component for $\Gamma=0$. However for non-zero $\Gamma$,
as we go from large antiferromagnetic $J$ to large ferromagnetic $J$ there is
a crossover between the states $E_{0}$ to $E^{n}$ at some finite $J$. Here
$E_{0}$ is the lowest energy state at $J\rightarrow-\infty$ and $E^{n}$ is
some positive high energy states at $J\rightarrow-\infty$. These two states
continually come closer to each other and for $\Gamma=0$, the asymptotically
merge but for finite $\Gamma$ there is a crossover between these two states at
some critical value of $J$.
In the appendix, we enlist all the eigenvalues and eigenstates for four site
cluster for finite $K,J$ and $\Gamma=0$ for future references. We also yield
exact expression for correlation function for comparison in the appendix. To
complete our discussion on eigenvalue spectra on finite size cluster we note
that the general features present in four site cluster prevail for larger site
clusters also. Mainly there exists two group of eigenstates one with linearly
varying eigenvalues and one with non-linearly varying with $J$. We always find
the number of linearly varying eigenvalues constitute the majority of
eigenstates and thus determine the finite temperature property.
In previous studies a lot of attention as given to the effect of magnetic
field on the Kitaev interaction in general and it is observed that it
stabilizes the Kitaev spin-liquid phase at intermediate temperature berke-
prb-2020 ; nasu-motome-2015 ; koga-nasu-2019 ; li-2021-ncom . Hence it is
important to examine what intricacies an external magnetic field brings in to
the eigenspectrum. We plotted the few low energy eigenstates for various
clusters as shown in Fig. 3.

Figure 4: Nearest neighbor correlation on a $z$-type bond in Kitaev clusters
in the presence of magnetic field $\vec{h}$ in the $\hat{z}$ direction is
plotted. A given column denotes a particular cluster as mentioned on the top.
Each row is plotted for a particular $(K,J)$ as written at the right side.
Only specific magnetic field with representative characteristics has been
plotted though magnetic field has been taken from 0.2 to 2 insteps of 0.2. The
color legends for magnetic field is shown at the left.

For the sake of simplicity we first start analysing the effect of external
magnetic field for $J=0$ and $J=-0.5$ as shown by the first, second and third
and fourth rows respectively. We observe that for these two choices of
Heisenberg interaction, when Kitaev interaction is ferromagnetic the
eigenvalues mainly increase with magnitude of $h_{z}$ as evident from the
second and fourth rows. However for the antiferromagnetic Kitaev interaction
the eigenvalues decrease or increase depending on the cluster size and
strength of external magnetic field. Notably for twelve site cluster for
$(J,K)=(-0.5,-1)$ (panel ‘$(\tilde{l})$’) we observe the gap between the
ground state and first excited states are gradually increasing with $h_{z}$
but for antiferromagnetic Kitaev interaction i.e for $(J,K)=(-0.5,1)$ (panel
‘$(l)$’) the ground state and the first excited states remain nearly
degenerate. Also in the former case the other excited states move away from
the ground state energy but in the later the low energy excited states move
closer toward the ground state energy. This highlights that the relative sign
and magnitudes of Kitaev and Heisenberg interaction compete very differently.
Similar analysis can also be done for antiferromagnetic Heisenberg
interactions. It is interesting to compare our result to that found in an
earlier study zhu-2018  where effect of an external magnetic field in the
direction of $(1,1,1)$ has been considered. Ferromagnetic and
antiferromagnetic Kitaev interaction have been shown to produce qualitatively
different results. In the case of ferromagnetic Kitaev interaction the first
excited state and other excited states shows monotonous increase of the gap
from the ground state energy which we also observe as evident from first row
Fig. 3. However for the antiferromagnetic Kitaev interaction a qualitatively
similar evolution of eigenenergies have been observed where higher excited
states come closer to the ground state energy and at some critical external
field it becomes nearly degenerate with the ground state zhu-2018 . From the
first row of Fig. 3 we observe a qualitatively different evolution of low
lying eigenstates with external magnetic field than found for ferromagnetic
Kitaev interactions. However the difference between the earlier result
zhu-2018  and in our case is attributed to the finite sizeholzmann-2016  and
the direction of magnetic field which has been taken in $(0,0,1)$ directions.

##  III Correlation function

In the pure Kitaev limit the two-spin correlation function vanishes beyond
nearest neighbor and also they are bond dependent i.e on $\alpha$-bond only
$\langle s^{\alpha}_{i}s^{\alpha}_{j}\rangle$ exists, a consequence of $Z_{2}$
conserved operators. On the other hand in the pure Heisenberg limit, the spin-
spin correlation function for nearest neighbor spins are isotropic i.e all
$\langle s^{\alpha}_{i}s^{\alpha}_{j}\rangle$ exist for any $\alpha=x,y,z$ and
they are equal, a consequences of $SU2$ symmetry. The characteristics bond-
dependent correlation function in the Kitaev limit do not survive when any
other interaction such as Heisenberg or $\Gamma$ or external magnetic field is
applied. As the nature of the bond-dependent short-range correlation function
at Kitaev limit does not depend on the system size and it is realized in exact
fashion even for small clusters, it is thus natural to ask what is the
manifestation of other non-Kitaev interaction such as Heisenberg interaction,
$\Gamma$ interaction and effect of external magnetic field. The previous
theoretical and many experimental findings have shown presence of a plethora
of phases for such generalized Kitaev-Heisenberg-$\Gamma$ modelkhwang-2022 ;
feng-2022 ; han-li-2023 ; yilmaz-2022 ; shuyi-2023 ; smit-2020 ; kimchi-2014 ;
yamaji-2016 ; young-baek-kim-2015 ; wang-2021 . There are long range ordered
various collinear phase, Kitaev-spin liquid phase and other quantum paramagnet
phase. It is generally understood that long-ranged classically ordered
collinear order arises mainly by local constraint to minimize the bond-
interactions which has some reciprocity with the Kitaev interaction that here
also the local $Z_{2}$ conserved operator induced constraint are the deciding
factor for short range correlation function. Thus we think it is indeed of
great relevance to investigate the correlation function even for small
clusters. However it is known that we will not be able to infer the complete
information regarding long-ranged correlation brought in by non-Kitaev
interaction such that the spin-liquid phase has long-ranged algebraically
decaying correlation function. Nevertheless the correlation function study for
the small cluster promises to show some signature of competition between
different kind of interactions. To begin with we first discuss the scenario
$\Gamma=0$.

###  III.1 $z-z$ correlation function

Without loss of generality we consider the bonds $(1,3)$, $(1,4)$, $(1,8)$
$(1,12)$ for the four site, six site, eight and twelve site cluster. We note
that they are all $z$-type bond as shown in Fig. 1.
Correlation function for $\Gamma=0$: Let us first analyze the correlation
function for $\Gamma=0$ and $h_{z}=0$. As we see from the red line plot that
for ferromagnetic or antiferromagnetic $J$, the $z-z$ correlation function
monotonically increases (decreases) for ferromagnetic(antiferromagnetic)
Kitaev interaction. This causes a point of inflexion at $J=0$. However we
observe that for four site cluster the correlation function changes
drastically with external magnetic field and become independent of whether
Kitaev interaction is ferromagnetic or antiferromagnetic. In both the cases
the correlation function saturates to $\pm 1$ for antiferromagnetic or
ferromagnetic Heisenberg interaction. However for $K=1$ as we increase the
value of the external magnetic field we find that after a critical magnetic
field $h_{z,c}$, the $z-z$ spin correlation function remains constant for
certain range of antiferromagnetic Heisenberg interaction. The ground state
does not changes much such that spin-spin correlation function remains
constant for a range of Heisenberg interaction. The above nature for $z-z$
correlation function changes as we increase the cluster sizes with two notable
differences. Firstly the zero field correlation function for antiferromagnetic
Heisenberg interaction yields a saturated $z-z$ correlation function. For
ferromagnetic interaction the $z-z$ correlation function monotonically
increases and reach a saturation whose upper bound is 1 for very large
external magnetic field. However, interestingly, for very large external
magnetic field, the spin-correlation-plateau is universally observed for all
cluster for certain ranges of antiferromagnetic Heisenberg interaction for
example consider $(N,K,h_{z})=(6,1,2),(6,-1,2),(12,1,2),(12,-1,2)$. In short
the main differences between the ferromagnetic and antiferromagnetic
Heisenberg interaction on the nearest-neighbor $z-z$ correlation function is
that for the former the $z-z$ correlation function behaves mostly monotonously
but for the later the $z-z$ correlation function shows plateau like
characteristics where $z-z$ correlation function remains constant for certain
range of antiferromagnetic Heisenberg interaction and the width of the
plateaus increase with increasing magnetic field.

Figure 5: Nearest neighbor $y-y$ correlation on a $z$-type bond in Kitaev
clusters in the presence of magnetic field $\vec{h}$ in the $\hat{z}$
direction is plotted. A given column denotes a particular cluster as mentioned
on the top. Each row is plotted for a particular $(K,J)$ as written at the
right side. Only specific magnetic field with representative characteristics
has been plotted though magnetic field has been taken from 0.2 to 2 insteps of
0.2. The color legends for magnetic field is shown at the left.

Correlation function for $\Gamma=\pm 0.5$: For the four size cluster, we
observe that for antiferromagnetic Kitaev interaction plateau-like property
continues to exist though for ferromagnetic Kitaev interaction it ceases to
exist. Otherwise the $z-z$ spin correlation function mostly increases
(decreases) as we increase the external magnetic field for antiferromagnetic
(ferromagnetic) Heisenberg interactions. For six and eight site cluster the
behavior of $z-z$ spin correlation function are qualitatively similar to that
of $\Gamma=0$, the plateau like behavior persists for multiple ranges of
antiferromagnetic Kitaev interactions and occasional oscillations appear near
$J=0$. Apart from this occasional oscillations which may arise due to
degenerate ground states, the $\Gamma$ interaction brings in two important
differences. First we note that for $\Gamma=0$, the asymptotic value of $z-z$
correlation function for large Heisenberg interaction is independent of the
strength of external magnetic field but in the presence of $\Gamma$ the
asymptotic value of $z-z$ correlation function for ferromagnetic Heisenberg
interaction depends on the external magnetic field and increases with the
strength of it. Secondly we find, for $K=1$, as we increase the cluster size,
the $\Gamma$ interaction brings the $z-z$ correlation function to lower value
for some small range of ferromagnetic Heisenberg interaction such that for
twelve site cluster the $z-z$ correlation function decreases for small values
of $h_{z}$. This initial decrease is more pronounced in the presence of
$\Gamma$ interaction. Thus we see multiple effects being in competition such
as quantum correlation due to increased size of the cluster as well as due to
the competition between Heisenberg and $\Gamma$ interaction. We note that for
$\Gamma=-0.5$, the $z-z$ spin-spin correlation function is qualitatively
similar to that of $\Gamma=0.5$, though for certain values of Heisenberg
interaction there is some differences. For example for small values of $h_{z}$
and for ferromagnetic Heisenberg interaction we observe a value of $J_{c}$ for
which the $z-z$ correlation function shows a maxima, this is pronounced for
six site cluster and eight and twelve site cluster for some intermediate value
of ferromagnetic $J$.

Figure 6:  Various plaquette operators are plotted against Heisenberg
interaction $J$. Each column denotes a particular cluster as mentioned at the
top. First row and second row correspond to antiferromagnetic and
ferromagnetic Kitaev interaction with $K=1$ and $K=-1$ respectively as
mentioned at the right side. The color legends show the value of magnetic
field taken. Inside each panel of the first row, the expression for particular
plaquette operator is given.

###  III.2 $y-y$ correlation function

In the Kitaev limit $y-y$ correlation function on a $z$-type bond has zero
value. However the Heisenberg and $\Gamma$ interaction brings in some non-zero
values to $y-y$ correlation functions. When $\Gamma=0,h_{z}=0$, $y-y$
correlation function saturates to some constant positive (negative) value for
very large ferromagnetic or antiferromagnetic Heisenberg interaction. For
intermediate ferromagnetic Heisenberg interaction there are some non-
monotonicity due to the competition of Kitaev and Heisenberg interaction. As
we apply external field $h_{z}$, the asymptotic value of $y-y$ correlation
function is zero for ferromagnetic Heisenberg interaction indicating that
ground state is polarized in the $z$-direction. But the asymptotic $y-y$
correlation for the large antiferromagnetic Heisenberg interaction is finite
even for large external magnetic field. This suggests that the ground state is
still a quantum ground state and this qualitatively marks the difference
between the effect of ferromagnetic and antiferromagnetic Heisenberg
interaction to the Kitaev interaction.
For intermediate strength of Heisenberg interaction interesting variations of
$y-y$ correlation appear which is counter intuitive. For ferromagnetic
Heisenberg interaction the $y-y$ correlation function first increase (opposite
to the naive expectation) with $h_{z}$ being non-zero which is remarkable.
Because naively we do not expect $y-y$ correlation to gain in strength due to
a field applied in $z$-direction. This suggests the quantum nature of the
ground state for intermediate strength of ferromagnetic Heisenberg interaction
and as it is departing from Kitaev limit the Heisenberg interaction is
manifesting itself in non-zero large $y-y$ correlation function. As we
increase the strength of ferromagnetic Heisenberg interaction, the $y-y$
correlation function starts to decrease. This results in a critical $J_{c,fm}$
for each $h_{z}$ where $y-y$ correlation function is maximum. For the
antiferromagnetic Heisenberg interaction we observe the plateau like behavior
for small values of $J$. The width of the plateaus increase as we increase the
strength of external magnetic field. It is important to note that this plateau
like behavior persists even for twelve site cluster.

Study of correlation function has attracted a lot of previous interest
tikhonov-2011 ; song-2016 ; lunkin-2019 ; yang-2020 ; berke-prb-2020  where
the spin-spin correlation function has been calculated in the presence of an
external magnetic field or other interactions. However either these studies
examined the long time and long distance behavior of correlation function
tikhonov-2011 ; lunkin-2019  or approximately for short distances yang-2020
perturbatively. On the other hand our exact analysis uncovered interesting
nature of correlations due to competition of Heisenberg and Kitaev
interaction. This reveals that the qualitative nature of nearest neighbor
correlation function for antiferromagnetic Heisenberg interaction is
completely different than ferromagnetic Heisenberg interaction in the presence
of external magnetic field which has not been explored before.

##  IV Expectation of plaquette conserved operator

The exact solvability of Kitaev model resides on the fact that it reduces to
non-interacting Majorana fermion hopping problem due to presence of conserved
plaquette operators which manifest itself as local $Z_{2}$ gauge field on a
given bond. These $Z_{2}$ gauge fields are bilinear in Majorana fermions. We
note that these plaquette conserved operators are product of multispin
operators and cover almost same size of the cluster for the four and six site
clusters. In the Kitaev limit they are conserved and static but in the
presence of non-Kitaev interactions such as Heisenberg and external magnetic
field they are no longer conserved and acquire dynamics. The $B_{p}$ operators
can always be expressed as
$B_{p}=\prod\sigma^{\alpha_{1}}_{1}...\sigma^{\alpha_{n}}_{n}$ where ‘$1$’ to
‘$n$’ denote the indices of sites on the plaquette in clockwise or
anticlockwise manner and $\alpha_{m}$ denotes the $m$-type outgoing Kitaev
bond at the site $m$. For the Kitaev cluster defined on a tetrahedra there are
three sites in an elementary plaquette while there are a four site conserved
operator as well which is given by
$B_{1-4}=\sigma^{z}_{1}\sigma^{z}_{2}\sigma^{z}_{3}\sigma^{z}_{4}$. For the
six and eight site cluster also the elementary plaquettes are squares. On the
other hand for the twelve site cluster the elementary plaquettes are hexagons.
The eigenvalue of the plaquette conserved operators are $\pm 1$ as evident
from the identity $B^{2}_{p}=1$. In the Kitaev limit when $J=0=h_{z}$, the
eigenstates of the Hamiltonian are simultaneous eigenstates of $B_{p}$’s.
However for the ground state, the eigenvalue of $B_{p}$ could be $+1$ or $-1$
depending on the total number of sites in the cluster. Various symmetries such
as time reversal symmetry, reflection symmetry also determine the ground state
eigenvalues of $B_{p}$ for various clusters lieb-1994 .

Figure 7: Dynamic correlation function $|g(t)|$ is plotted as a function of
time $t$ for $\Gamma=0$. Each column represents Majorana dynamics for various
cluster as mentioned at the top. In the right we denoted the Heisenberg
strength $J$ for which $|g(t)|$ has been plotted. Odd rows and even rows
represent plots for antiferromagnetic and ferromagnetic Kitaev interaction
respectively. Various color legends explain the magnitude of external magnetic
field $h_{z}$.

We begin our discussion with tetrahedral system as shown in Fig. 6, panel
$(a)$ and $(f)$ for $K=1$ and $K=-1$ respectively. For this system ground
state is degenerate with $B_{p}=\pm 1$ as being product of three spin operator
it is odd under time reversal symmetry. However we choose $B_{p}=1$ for
conveniences which has been used in Fig. 6 for $J=0,h_{z}=0$. However for non-
zero $J$ such a state selection is not possible and the exact ground state
eigenfunction determines its expectation value. The value of $B_{p}$ being
zero at finite $J$ for $K=\pm 1$ denotes that ground state eigenfunction is
equal superposition of $B_{p}=\pm 1$. This can also be explained on the basis
that Hamiltonian is time reversal invariant but the $B_{p}$ operators are not
and thus its expectation value is zero. In panel $(b)$ and $(g)$ of Fig. 6 we
plotted expectation of $B_{1-4}$ which can be obtained by the multiplication
of two plaquette operators associated with triangles. Thus it can be thought
of as a correlation between the two adjacent plaquette operators associated
with two triangles. We observe that at $J=0$, its value is $+1$ by our choice.
But for antiferromagnetic Heisenberg interaction it changes its value to -1
implying that the two $B_{p}$ associated with triangles are anti-correlated
and this constant value increases with increasing magnetic field. After a
critical $J_{afm}$, $S_{1-4}$ reduces to zero implying that the conserved
quantities associated with triangles are no longer correlated. It manifests an
unique competition between the quantum paramagnet states pertaining to Kitaev
limit and quantum paramagnet states in the Heisenberg limit.
For the six site clusters the $B_{p}$ shows $+1$ for $J=0,h_{z}=0$. There are
six four site plaquettes and one six site plaquette and here we have plotted
the $B_{p}$ associated with the four site plaquette. Each two four site
plaquette yields the six-site plaquette. Due to absence of reflection symmetry
for the elementary square we find the expectation value of plaqeutte operator
as 1 instead of -1 (which should be according to Lieb theorem lieb-1994  if
reflection symmetry is present). However for the six site plaquette operator
we obtain the expectation value $1$ at the Kitaev limit i.e $h_{z}=0=J$. For
nonzero $J$ and $K=1$ we observe that there is a window of $J$ for which
$B_{p}$ expectation takes a constant positive value and it goes to zero after
a critical ferromagnetic $J$ whose magnitude increases with external magnetic
field. For ferromagnetic Kitaev interaction this critical Heisenberg
interaction shift to antiferromagnetic $J$ suggesting an interesting
competition between the relative sign of $J$ and $K$.
Having described four site and six site clusters in detail we now explain the
ground state $B_{p}$ expectation value for eight site and twelve site cluster.
We note that the elementary plaquettes in eight site cluster is four and hence
Lieb theorem gives $B_{p}=-1$ for the ground state. We find that for
$J=0,h_{z}=0$ $B_{p}$ to be $-1$ as expected as seen from Fig. 6. For $K=1$
any large values of $\pm J$, $B_{p}$ vanishes as expected because of
eigenstates mixing with all eigenstates of $B_{p}$ in the presence of $J$.
However for small value of $\pm J$ and $K=1$, there are high values of $B_{p}$
for some values of $h_{z}$. This probably indicate different plaquette
operator transform differently under external magnetic field. We note that in
the eight site cluster there are plaquettes operators which contains only
$\sigma^{z}$ and hence commutes with the external magnetic field. Thus there
are preferences for some eigenstates to be selected and this causes finite
large $B_{p}$ for small $\pm J$. However there is a remarkable differences for
ferromagnetic Kitaev interaction in comparison to antiferromagnetic Kitav
interaction. For $K=-1$ for large $-J$ the $B_{p}$ saturates for $h_{z}=0$ and
also the half-width maxima of $B_{p}$ for other $h_{z}$ is broader in
comparison to antiferromagnetic $K$. For the 12 site cluster we observe that
for $J=0,h_{z}=0$ and $K=\pm 1$ $B_{p}=1$ satisfying Lieb’s theorem. For
finite $h_{z}$ and antiferromagnetic $K$, $B_{p}$ shows a rapid oscillations
(having positive and negative values of $B_{p}$) for small window of $|J|$ and
vanishes thereafter. The same is true for ferromagnetic $K$ though for small
value of $J_{z}$, the $B_{p}$ shows a saturation for large ferromagnetic $J$
zhu-2018 ; han-li-2023 ; nandini-2019 ; liang-2018 ; yilmaz-2022 .

##  V Dynamics of Majorana fermions

Much of the recent interest in explaining various experimental observation and
theoretical results was devoted to understand the gauge field dynamics in
Kitaev-Heisenberg-$\Gamma$ model. It is well known that a pure Kitaev model is
described by a non-interacting Majorana fermion hopping problem
($H=\sum_{ij}J_{ij}ic_{i}c_{j}$) in the background of conserved $Z_{2}$ gauge
fields $u_{ij}=ib_{i}b_{j}$. In appropriate representation these conserved
gauge field operators can be re-written as an equivalent number operator as
$u_{ij}=(2\chi^{\dagger}\chi-1)$, named bond-fermion mandaljpa . Also the
Majorana operator ($c_{i}$) can also be appropriately recasted interms of
complex fermion $\psi$ which on some bonds appear as
$2\psi^{\dagger}_{i}\psi_{i}-1$. Depending on the interest there were much
effort investigating the gauge field and matter field dynamics by studying the
expectation value $<ic(t)_{i}c(0)_{j}>$ and also $<ib_{i}(t)b_{j}>$ in the
presence of other non-Kitaev interaction such as Heisenberg interaction,
$\Gamma$ interaction and external magnetic field knolle-2018 ; nasu-
motome-2019 . These external perturbations make the gauge fields to acquire a
dynamic as they are no longer conserved and a quantitative analysis of these
matter field and gauge field dynamics enable us to understand the possible
spin-fractionalization leading to spin-liquid phase appearing at intermediate
temperature berke-prb-2020 ; nasu-motome-2015 ; koga-nasu-2019 ; li-2021-ncom
. This is also intimately connected to confinement-deconfinement transitions
associated with magnetically ordered state to quantum paramagnetic states
including Kitaev spin liquid state. To study the gauge field and matter field
dynamics here we use much simpler quantities ($g_{\phi}(t)~{}{\rm
and}~{}g_{\chi}(t)$) defined below which are equally capable to characterize
the dynamics of Majorana fermions brought in by non-Kitaev interaction.

$\displaystyle g_{\phi}(t)=\langle
N_{\phi,i}(t)N_{\phi,i}(0)\rangle,~{}~{}~{}N_{\phi,i}=\phi^{\dagger}_{i}\phi_{i}$
$\displaystyle g_{\chi}(t)=\langle
N_{\chi,i}(t)N_{\chi,i}(0)\rangle,~{}~{}~{}N_{\chi,i}=\chi^{\dagger}_{i}\chi_{i}$
(2)

The reason that we choose this operator is that $N_{\phi,i}$
($N_{\chi,i}$)acts as a projector on the ground state such that it projects
only states with ‘$\phi$’ ($\chi$) fermion on the bond ‘$i$’. Thus
$g_{\phi}(t)$ (or $g_{\chi}(t)$) yields the probability amplitude of those
states to return back to itself. In the Kitaev limit as $N_{\chi,i}$ is static
and thus $g_{\chi}(t)$ is constant in time but acquire time dependence once
non-Kitaev terms are considered. In the presence of non-Kitaev interactions
$\psi$ and $\chi$ interact with each other thus estimating those quantities in
Eq. V yield an understanding of the dynamics associated with non-Kitaev
interactions. In Fig. 7 and Fig. 8, we show how $g_{\chi}(t)$ and
$g_{\phi}(t)$ depends on time for various cluster and other system parameters
for $\Gamma=0$. Below we summarize some salient characteristics of gauge and
matter Majorana fermion dynamics.

Figure 8: Dynamics of matter Majorana fermion as defnied by the quantity
$g_{\phi}(t)$ is plotted for various cluster size and parameter values.
Various panels describe a particular cluster, a given combination of Kitaev
interaction strength $K$, Heisenberg interaction strength $J$ as described in
Fig. 7.

###  V.1 Dynamics of matter and gauge fermions at $\Gamma=0$

To describe the dynamics of majorana fermions we consider bond (1,3), (1,4),
(1,4), (1,6) for four site, six site, eight site and twelve site cluster. Thus
in the below text whenever we mention $\chi$ or $\phi$ fermion, it is to be
understood that they are defined on this bonds. Occasionally for eight site
and twelve site cluster we mention $\chi_{2}$ and $\phi_{2}$ and they are
defined on bonds (2,7) and (2,9) for eight and twelve site cluster
respectively. Qualitative differences between four site cluster and other
clusters: In all the parameter regime there are qualitative differences in the
dynamics of gauge fermions between the four site cluster and other larger
clusters. For $\Gamma=0,J=0$ the dynamics of gauge fermions are very regular,
periodic as we vary the external magnetic field. For intermediate
antiferromagnetic $J$, the oscillations are very regular and amplitude
decreases monotonously as we increase the strength of $h_{z}$ (see Fig. 7,
panel ‘$i$’), though for ferromagnetic $J$ appearances of additional frequency
could be noticed which points out that the ground state is modified
differently than the antiferromagnetic case.
The same qualitative difference as discussed above can be seen for matter
fermions also. However we notice that in general the average amplitude of
oscillations for gauge fermions are always larger than the matter fermions.
This can be attributed to non-local nature of the gauge fermions interms of
the spin-operators and also to the fact that each flux sectors are separated
by an energy gap which is much larger than the fermionic spectrum for each
flux sector.
Regular oscillations and irregular oscillations: It is observed that depending
on the strength of external magnetic field the oscillation of gauge fermions
could be regular to irregular type. For example consider six site cluster
$J=0$. Up to $h_{z}=0.6$, the oscillations are very small but after that it
suddenly increasesyamada-2020  and reaches maximum for $h_{z}=1.2$ and then
again it decreases, this has been shown in Fig. 7, panel $(b)$. In numerous
occasions it have been observed that $h_{z}$ can reverse the response of gauge
fermions amplitude or time period of oscillations. On the other hand for eight
site cluster a different scenario is observed comparing $J=-0.5,0,0.5$ i.e
Fig. 7, panel $(k)$. For $J=-0.5$ the amplitude of oscillation decreases with
increase of $h_{z}$, this behavior is opposite to that found for $J=0$ and for
$J=0.5$, the oscillations remain more or less same for all $h_{z}$ (see Fig.
7, panel $(g)$ ), no special feature is seen. However for twelve site cluster
such non-monotonous response is absent. very interestingly we see from Fig. 7
panel $(c),(d),\tilde{c},\tilde{d}$ that the gauge fermion does not show any
dependency on time for any magnetic field. This can be referred as
stabilization of gauge fields under external magnetic field on certain bonds.
However on a different bonds the gauge fermion dynamics is oscillatory as
shown by extra four panel at the right side of the main panel. The above
controlled manipulation of gauge field dynamics could be of practical use.
The identical trends can be seen for matter fermions as well where depending
on the relative strength of external magnetic field $h_{z}$ and Heisenberg
interaction $J$. For example one can compare
$(N,K,J,h_{z})=(6,-1,0,0.4)~{}{\rm and}~{}(6,-1,0,1.2)$. In the former case
there are many minimum found but for the later only a single minima found, see
Fig. 8, panel $(b)$, the nature of oscillation completely changed. Similarly
for $(N,K,J,h_{z})=(6,-1,0.5,0),(6,-1,0.5,0.4),(6,-1,0.5,0.6)$ panel
$(\tilde{f})$ in Fig. 8(similar trends are seen for $\phi_{2}$ fermion as
well). These all three set of parameters yield three different behavior for
$\phi_{2}$ fermion. For eight site cluster such differences are also found for
$(K,J)=(1,-0.5),(1,0),(1,0.5)$ (see Fig. 8, panel
$(\tilde{c}),(\tilde{g}),(\tilde{k})$ and also for $(K,J)=(-1,0),(-1,0.5)$
where depending on external magnetic field $\phi_{1}$ fermions dynamic changes
a lot (see Fig. 8, panel $(k)$). Similarly for twelve site clusters the
relevant set of parameters for such change in the nature of $\phi_{2}$
fermions can be seen for $(K,J)=(1,-0.5),(1,0)$.

Dependency of average amplitude and time period of oscillations on $h_{z}$ and
$J$: Interestingly it has been noticed that for $J=0$, the time period of
oscillations is much larger for $K=-1$ than $K=1$ case, this can be noticed
from the first and second rows of Fig. 8. In the first row, for finite
$h_{z}$, the oscillations has more amplitude and time period than the second
row. However as we turn on finite $J$ the time period of oscillations become
comparable though it greatly depends on the relative strength of $K,J$ and
$h_{z}$. It is also observed that the external magnetic field helps in
reviving the amplitude of matter fermion on that site. For example consider
twelve site cluster where the rapid fluctuations are absent in comparison to
other smaller cluster. If we look at $(K,J)=(-1,0.5)$, Fig. 8, panel
$(\tilde{h})$ where in between $t=1$ and $t=3$, there is a revival of
amplitude of matter fermions. Similarly we find for $(K,J)=(-1,0),(\pm
1,-0.5)$ in panel $(d)$ with $h_{z}=0.4$ and panel $(\tilde{l})$ and $(l)$.
Similarly we see that for gauge fermions, $K=-1,J=0$ corresponds to larger
time period in comparison to $K=1,J=0$ (see Fig. 7 first and second row ).

Nature of oscillations for large $\pm|J|$: Nature of oscillations in large $J$
(ferromagnetic of antiferromagnetic) limit is more regular. Yet depending on
the relative sign of $J$ and $K$ there could be presence of one single
frequency or more. This is reflected in the amplitude also. In some cases
amplitude is constant when there is a single frequency and amplitude varies
when multiple frequencies are present. For four site cluster large positive
$J$ the oscillation contains single frequency (independent of $h_{z}$) though
large $-J$ contains multiple frequency centered around two main frequencies
for non-zero $h_{z}$. For six site cluster the oscillations are still regular
though multiple frequencies appear but for eight site cluster it becomes
regular due to special symmetry of the cluster. For twelve site cluster there
are presence of many frequencies but they vary slowly with $h_{z}$ so that the
oscillations appear coherent. For matter fermion, the amplitude of
oscillations are always smaller in comparison to gauge fermions except special
cases where the external magnetic fields are seen to stabilize the gauge
fields on some particular bond.

Two different channels of dynamics of oscillations: It is also noted that
depending on initial position of the matter fermions its time dependency could
be very different. For comparison we note the results for twelve site cluster
with the parameter value $(K,J)=(-1,0.5)$. We observe that whether $\phi_{1}$
is decaying with rapidly, $\phi_{2}$ is having quite small oscillations and
vary smoothly. The average magnitude of $g_{\phi}(t)$ corresponding to
$\phi_{2}$ also gradually increases consistently with increasing $h_{z}$. This
different behavior depending on the initial position of the matter fermion
indicates the local competition between the different interaction such as
Heisenberg and Kitaev. The same is true for ferromagnetic Heisenberg
interaction with parameter $(K,J)=(-1,-0.5)$. There is a marked difference on
the dynamics of $phi_{1}$ and $\phi_{2}$. The same can be said for $J=0$. This
observation could be very much useful for the controlled manipulation of
Majorana fermions for useful quantum state operations.
Stabilization of $Z_{2}$ gauge field by external magnetic field: It has been
observed that as the applied magnetic field is in $z$-direction it can
stabilize the $Z_{2}$ gauge fieldkhwang-2022  as seen in for eight site
cluster ($K,J$)= ($-1,-2$) or ($-1,-0.5$), Fig. 7, panel $(\tilde{k})$ and
also for $\tilde{l}$ and other). For intermediate value of $J$ positive the
gauge field fluctuation decreases as we increase $h_{z}$ as seen for ($K,J$)=
($-1,0.5$), Fig. 7, panel $(\tilde{f})$. This can be thought of as stabilizing
the spin-liquid or fractionalized state by external magnetic field. This is
also true for $K=1$ though there are some fluctuations of the gauge field
around values closed to unity. Identical trends are seen for twelve site
cluster, a finite $h_{z}$ greatly stabilizes $Z_{2}$ gauge fields on certain
bonds. This can be seen for $(K,J)=(1,-2),(1,-0.5)$ (Fig. 7, panel $(l)$) and
$(K,J)=(-1,-2),(-1,-0.5)$(Fig. 7, panel $(\tilde{l})$). For antiferromagnetic
$J$, the $Z_{2}$ gauge fields show more oscillatory pattern in comparison to
ferromagnetic $J$, they become closer to unity.
Scenarios to define quantum speed limit: Recently various aspect of quantum
speed limit as defined by the evolution of a quantum state from a given state
to a distinct different state attracted wide interest for various reasons
deffner-2017 ; gal-ness-2022 ; pati-2022 ; arindam-2022 . In this context we
find that quantum-speed limit could be investigated as far as the evolution of
Majorana fermions are considered. Initially they begin with a given occupancy
at a given bond and if at any future time, the occupation of the fermion on
that bond becomes zero, we can define quantum speed limit for our system. In
our study of dynamics of gauge and matter fermion we witnessed numerous
occasions where such situations occur. For reference purpose we denote by
$(N,K,J,h_{z})$ the set of parameter where $N$ denotes the number of sites in
the cluster, $K$ and $J$ are the strength of Kitaev and Heisenberg interaction
respectively and $h_{z}$ being the strength of external magnetic field applied
in $z$-direction.
Most fascinating occurrence of quantum speed limit happens for four site
cluster with $(N,K,J,h_{z})=(4,-1,0,2),(4,-1,0.5,0.8)$. This happens in both
the matter and gauge fermion sector. For the second set of parameters they are
remarkably repeating with time. For $K=-1,J=-0.5$ the quantum speed limit
phenomena can also be observed depending on the strength of external magnetic
field. On the other hand for antiferromagnetic Kitaev interaction such
realization does not happen except for very large ferromagnetic Heisenberg
interaction. From this we note that when the sign of Kitaev interaction and
Heisenberg interaction is opposite then only it happens.
For the six site cluster and for the gauge fermions we find
$(N,K,J,h_{z})=(6,1,0,1.2)$ being such an example. For large antiferromagnetic
$J$ it also happens for matter fermions. For eight site cluster also there are
occasions for $\phi_{2}$ with
$(N,K,J,h_{z})=(8,1,0,1.6),(8,1,0,1.4),(8,1,0,1.8)$. There are occasions for
gauge fermions also for such phenomena as well as for ferromagnetic Kitaev
interactions for example $(N,K,J,h_{z})=(8,-1,0,0.6)$. For twelve site cluster
also there are special cases where such phenomena are observed. A careful
study and estimation of this quantum speed limit will be done in future study
as this is not the primary focus here.

Figure 9:  Zero field spontaneous magnetization has been plotted for various
cluster. Each row denotes a particular cluster as denoted at the top. At the
right side of each row we mention the value of $(K,\Gamma)$ taken.
Magnetisation has been plotted for various values of Heisenberg interaction
($J$) as denoted by the color legends at the left.

The main differences in the dynamics of majorana fermions between twelve site
cluster and other clusters: The most obvious difference between the twelve
site and other smaller cluster can be understood by analyzing $J=0$. As we see
that the gauge fermion remains constant (for cluster of all sizes) with time
but the matter fermion oscillates with a fixed frequency and returns back to
unity after characteristic time for four site and eight site cluster due to
the back flow from the boundary. On the other hand for six site cluster the
matter fermions return back to unity after few local maxima though for twelve
site cluster the matter fermions never returns back to unity. Secondly for
finite $J$, the oscillations in $\chi$ or $\phi$ fermions are having regular
pattern for different $h_{z}$ but for smaller cluster the pattern of dynamics
vary greatly depending on the magnitude of the external magnetic field.

The appearance of multiple frequencies in the dynamics of Majorana fermions
can be understood as follows. For a fixed gauge field configurations, the
effective Hamiltonian for $\psi$ fermions reduce to a superconducting
Hamiltonian with the presence of $\eta_{ij}\phi_{i}\phi_{j}+h.c$. The ground
state $|G\rangle$ is constructed as the vacuum of the quasi-particles where
the states with even number of $\psi$ fermions (including zero) appear with
different probability amplitude. Now the operator
$N_{\phi,i}=\phi^{\dagger}_{i}\phi_{i}$ selects (or projets) only those states
with the dimer ‘$i$’ being occupied. These daughter states could be re-
expressed as linear combinations of all the eigenstates (including excited and
ground states) with different amplitudes. Thus the final expression of
$g_{\phi}(t)$ contains the summation of $\sum_{n}A_{n}\exp^{-(E_{n}-E_{0})}$
where the detail analytic expression of $A_{n}$ is in general very complex.
With system size very large, a large number of contribution to $A_{n}$ coming
from numerous daughter states makes the probability amplitude oscillate with
no fixed time period.
Now for finite non-Kitaev interaction being present the gauge fermions are no
longer conserved and hence they also oscillate with time. For the four site
cluster due to special symmetry the matter fermion and gauge fermion returns
back to unity with identical frequency. As the cluster size is increased, the
oscillations of matter and gauge fermions become irregular due to
contributions from increasing number of excited states for larger clusters.
However few universal characteristic features are present. One notes that the
on an average the gauge fermion amplitude is always larger than the matter
fermion and closer to unity. Also matter fermions complete many cycles of
oscillations within a given time window in comparison to gauge fermions. This
points out that gauge fermions can be still approximated as slow varying
background conserved operators for the matter fermions and hence may support
survive of decondined phase.

##  VI Magnetization

Having discussed various zero temperature properties including the eigen-
spectrum analysis, nearest-neighbor correlation function, expectation of flux
operators and detailed dynamics of Majorana fermions we now focus on some
finite temperature properties of various clusters. Till now we have found that
various competing interactions such as Kitaev interaction ($K$), Heisenberg
interactions($J$) and $\Gamma$ interactions have effected the above mentioned
zero temperature and ground state properties in distinct ways. Thus it is
natural to investigate some finite temperature properties to see what are the
distinct ways each of those different interaction effect. Here we investigate
magnetization, susceptibility and specific heat as described below.

###  VI.1 Spontaneous magnetization

Figure 10: Finite temperature magnetization has been plotted for various
external magnetic field. Each column describes result for a specific cluster
as mentioned at the top. Odd rows are for antiferromagnetic Kitaev interaction
$K=1$ and the even rows are for ferromagnetic Kitaev interaction. The value of
Heisenberg interaction taken for each rows are mentioned at the right of each
row.

To calculate the spontaneous magnetization in the absence of external magnetic
field, we considered a very small external magnetic field in the Hamiltonian
and calculated the spontaneous magnetization along $\hat{z}$ direction. In the
Fig. 9 we have plotted the magnetization-vs-temperature for various
parameters. In the below we describe some of the salient features.

For ferromagnetic Kitaev interaction, $K=1$ : The zero field spontaneous
magnetization reveals that for $\Gamma=0,~{}J=0$ the magnetization linearly
increases as we decrease the temperature from a very high value. However there
is an intermediate temperature $T^{*}$ below which different scenarios may
happen depending on the sign and magnitude of Heisenberg coupling $J$ as well
as $\Gamma$ term. For $\Gamma=0$ and ferromagnetic $J$, the linear increase of
magnetization with decreasing temperature saturates to a constant value at
some critical temperature that we call $T_{c,fm}$. On the other hand for anti-
ferromagnetic Heisenberg interaction, there exist two critical temperatures
what we define as $T_{>,afm}$ and $T_{<,afm}$ with $T_{<,afm}<T_{>,afm}$ such
that below $T_{>,afm}$ the zero field magnetization starts to decrease until
$T_{<,afm}$ is reached and then it saturates to a constant value
$\mathcal{M}$. We note that the constant magnetization $\mathcal{M}$ generally
depends on the strength of magnitude of Heisenberg interaction and in general
increase or decrease with the increase of magnitude of ferromagnetic or
antiferromagnetic Heisenberg interaction.
Having explained the behavior of magnetization for $\Gamma=0$ we now describe
what happens for finite $\Gamma$. We note that generally at high very high
temperature, the magnetization is independent on the sign and magnitude of
$\Gamma$. However below the critical temperature $T^{*}$ the effect of
$\Gamma$ can be seen. For ferromagnetic $\Gamma$ and $J$, the magnetization
increases linearly with decrease in temperature and deviates from the
saturating trends below $T^{*}$ for $\Gamma=0$. However for antiferromagnetic
$\Gamma$ and ferromagnetic $J$ the the trends to reach saturating
magnetization remains intact though the value of saturation magnetization
decrease with antiferromagnetic $\Gamma$.
Very interestingly for antiferromagnetic $J$, the behaviors of magnetization
seems to be indifferent of the sign of $\Gamma$. The dome like behavior of
magnetization in between $T_{<,afm}$ and $T_{>,afm}$ exist though the
saturation magnetization decreases in the presence of $\Gamma$ irrespective of
its sign.
For antiferromagnetic Kitaev interaction, $K=1$ : The overall qualitative
nature of variation of magnetization for antiferromagnetic Kitaev interaction
in the presence of $J$ and $\Gamma$ with their possible signs and magnitudes
are similar to those of ferromagnetic Kitaev interaction. All the critical
temperature $T_{c,fm},T_{<,afm},T_{>,fm}$ are present characterizing the
particular behavior of magnetization for anti-ferromagnetic Kitaev
interaction. However it is found that for ferromagnetic $K$, the value of
$T_{c}$ decreases as evident from comparing Fig. 9 ($b_{8}$) and Fig. 9
($\tilde{b}_{8}$) for eight site cluster. Similar observations are made for
twelve site cluster by comparing Fig. 9 ($b_{12}$) and Fig. 9
($\tilde{b}_{12}$). Though $T_{>}$ seems to be remain identical to that of
ferromagnetic Kitaev interaction, the value of $T_{<}$ seems to decrease a
little for ferromagnetic Kitaev interaction. The most remarkable difference
between the ferromagnetic and antiferromagnetic Kitaev interaction is that the
spread of saturation magnetization (or the band width of magnetization) at
$T=0$ decrease for ferromagnetic interaction signifying the role of quantum
fluctuation playing a dominant role for antiferromagnetic Kitaev interaction.
Also we have noticed a decrease in bandwidth of magnetization at zero
temperature as we increase the size of the cluster which is expected as due to
quantum confinement quantum fluctuation plays a dominant role in smaller size
clusters.

Figure 11: Finite temperature susceptibility has been plotted for various
external magnetic field. Each column describes result for a specific cluster
as mentioned at the top. Odd rows are for antiferromagnetic Kitaev interaction
$K=1$ and the even rows are for ferromagnetic Kitaev interaction. The value of
Heisenberg interaction taken for each rows are mentioned at the right of each
row. In the inset we have plotted susceptibility for $J=\pm 2$ as mentioned at
the right of the figure.

###  VI.2 Magnetization in the presence of external field

Having described the spontaneous magnetization at zero field, we now proceed
to briefly describe the magnetization at finite external magnetic field. The
finite field magnetization has been shown in Fig. 10. Here we have taken
external magnetic field ranging from $0.2$ to $2.0$ insteps of $0.2$ but only
those magnetization with a given external field has been shown for which there
is a characteristic differences than the other. Below we describe some
essential features.
Ferromagnetic ($K=-1$) vs antiferromagnetic Kitaev ($K=1$) for $J=0$: In the
absence of Heisenberg interaction we find that there is a characteristic
difference between the ferromagnetic and antiferromagnetic Kitaev interaction
as evident from the first two rows of Fig. 10. The dome like structure
appearing in the intermediate temperature in the absence of external magnetic
field only appear for antiferromagnetic Kitaev interactions for small to
intermediate magnetic fields. However for ferromagnetic Kitaev interaction the
dome like behavior in the magnetization in the intermediate temperature are
completely absent even for small external magnetic field. The other important
difference is that the saturation magnetization for given external magnetic
field is higher for ferromagnetic Kitaev interaction than the
antiferromagnetic Kitaev interaction. The saturation magnetization is almost
independent of external magnetic field for ferromagnetic Kitaev interaction
for all cluster. But for antiferromagnetic Kitaev interaction the saturation
magnetization greatly depends on the cluster size and gradually increases with
cluster size.
Ferromagnetic ($K=-1$) vs antiferromagnetic Kitaev ($K=1$) for $J=0.5$: In the
presence of antiferromagnetic Kitaev interaction we find that when $K=1$ then
the dome like structure in the magnetization is more pronounced and even exist
for very high magnetic field. The most remarkable things for $K=1$ is that
there are two limits of saturation magnetization $\mathcal{M}_{1}$ and
$\mathcal{M}_{2}$ where $\mathcal{M}_{2}>\mathcal{M}_{1}\geq 0$. The low and
intermediate magnetic fields yield a constant saturation magnetization
$\mathcal{M}_{1}$ but the high magnetic field saturation magnetization yields
another constant magnetization $\mathcal{M}_{2}$. On the other hand for
ferromagnetic Kitaev interaction the saturation magnetization is different for
each magnetic field though there is a re-appearances of dome like structure
for small magnetic field. Thus we see remarkable characteristic differences in
the magnetization profile for ferromagnetic and antiferromagnetic Kitaev
interaction and the fact that dome like structure in the magnetization
reappear in the presence of antiferromagnetic Heisenberg interaction points
out that the quantum frustration brought in by the antiferromagnetic
Heisenberg interaction is responsible for that.
Ferromagnetic ($K=-1$) vs antiferromagnetic Kitaev ($K=1$) for $J=0.5$: For
the ferromagnetic Heisenberg interaction we observe that the zero temperature
magnetization limit for $K=1$ depends on the external magnetic field and
almost proportional to the external magnetic field. The dome like structure
exist though becomes weaker. For the ferromagnetic kitaev interaction the zero
temperature saturation magnetization is very high and becomes independent of
external magnetic field after a critical external magnetic field. The dome
like structure for the intermediate temperature range is completely absent.
Thus we see that effect of competing interactions manifested distinctly in the
magnetization profile.

##  VII Susceptibility

Having described the behavior of magnetization in great detail we now proceed
to discuss the susceptibility. As before we find characteristic signatures of
competition among different interactions. In Fig. 11 we have plotted
susceptibilities for various clusters and different combinations of
parameters. As before though we have investigated susceptibility for finite
magnetic field $h_{z}$ from $0.2$ to $2.0$, we have only plotted the
representative one to describe how susceptibilities vary with magnetic fields.
We fist begin with a comparison of susceptibility for $J=0$.
Ferromagnetic ($K=-1$) vs antiferromagnetic Kitaev ($K=1$) for $J=0$: Let us
first describe how the susceptibility varies when Heisenberg interaction is
absent and what is the differences between ferromagnetic and antiferromagnetic
Kitaev interactions. As we see from the second row of the Fig. 11, the
susceptibilities for different clusters for ferromagnetic Kitaev interaction
are similar. The susceptibility is zero at very high temperature and maximum
at some intermediate temperature and becomes constant as temperature tends to
zero. The peak height gradually decreases with increasing magnetic field and
also appears at higher temperature as expected. The saturation magnetization
also depend on magnetic field and decreases with increasing magnetic field.
However for antiferromagnetic Kitaev interaction the effect of external
magnetic field on susceptibilities are completely different. First of all,
there is no divergences in the susceptibilities for any magnetic field which
has been observed for ferromagnetic Kitaev interaction. Also the behavior of
susceptibility is not uniform for different magnetic field over the cluster
sizes. For example for the twelve site cluster, the susceptibilities can be
high or low at a given temperature depending on the external magnetic field.
Identical trends can bee seen for eight site cluster. However though there are
peaks in the susceptibilities at a given magnetic field for twelve site
cluster, for eight site clusters peaks are very mild and only appear for small
magnetic field. Owing to the small sizes and intricate nature of ground state,
the susceptibility for four and six sites cluster shows complex behavior.
Ferromagnetic ($K=-1$) vs antiferromagnetic Kitaev ($K=1$) for $J=0.5$: An
antiferromagnetic Heisenberg interaction brings in a very important difference
to the susceptibilities of ferromagnetic Kitaev interaction which is that the
divergences found for $J=0$ at low magnetic field now appear at some finite
magnetic field for twelve site cluster. On the other hand four and eight site
cluster there are no pronounced maxima at any given magnetic field but such
pronounced maxima appears for six site cluster. In a similar way we find that
though in the absence of Heisenberg interaction there is no pronounced maxima
for antiferromagnetic Kitaev interaction, an antiferromagnetic Heisenberg
interaction causes susceptibility to have pronounced maxima or diverge at some
large magnetic field. Also at low temperature the susceptibilities converge to
very low values for antiferromagnetic Kitaev interaction but for ferromagnetic
Kitaev interaction it saturates to quite large values. Thus we see clear
evidences of competing effect on the susceptibilities due to Heisenberg and
Kitaev interaction.
Ferromagnetic ($K=-1$) vs antiferromagnetic Kitaev ($K=1$) for $J=-0.5$: In
the presence of ferromagnetic Heisenberg interaction we observe that for
$K=-1$, the peak structure re-appears in susceptibilities and it is divergent
mostly for small magnetic field. As we increase the magnetic field, the peak
height gradually decreases across all clusters and at very low temperature
they tend to saturate at zero for intermediate to large magnetic field as
evident from the last raw of Fig. 11. For the antiferromagnetic Kitaev
interaction the pronounced peak appears for twelve site cluster at some
intermediate critical magnetic field. Except this critical magnetic field, the
susceptibilities are generally regular though at a given temperature the value
of susceptibilities is not monotonously increasing or decreasing on external
magnetic field. For four and six site cluster pronounced peaks re-appear at
very small magnetic field but for eight site cluster there is no peak like
structure rather there is a mild dome like structures. Thus we observe that
external magnetic field greatly influences the susceptibilities. These detail
observations could be helpful in understanding the experiments in real
material.

Figure 12: A cartoon picture of the electronic spectrum in the Kitaev limit.
Here $\rm{B_{1},B_{2},B_{3},B_{4}}$ denote the flux configurations and the
rectangular bar denote the eigenspectrum for a given flux configuration. The
minimum of energy for each flux configurations as well as the separation
between two successive energy levels depends on the flux configuration. Figure
13: Normalized specific heat $C_{v}$ is plotted in the absence of magnetic
field for various values of Heisenberg interaction $(J)$ as shown inside the
panel (a). Each column denotes a particular cluster as mentioned at the top.
Solid lines denote antiferromagnetic Kitaev interaction and dashed lines
denote the ferromagnetic Kitaev interaction with $K=1,-1$ respectively.

##  VIII Specific heat

We now proceed to analysis the specific heat for the Kitaev-
Heisenberg-$\Gamma$ system. Much of the recent experimental results invested a
great deal of analysis of some specific characteristics of specific heat for
this system as it can effectively signify the existence of spin-liquid phase
or magnetically ordered states. A theoretical scenario for eigenspectrum at
the Kitaev limit is shown in Fig. 12. The effective Hamiltonian is described
by various conserved flux configurations and each flux configurations yield a
eigenspectrum as shown in Fig. 12. Generally the lowest energy for each flux
configuration is different for different flux configurations. In the Fig. 12,
we have denoted by $T_{1}$ and $T_{2}$, two temperature such that
$T_{2}>T_{1}$. At low temperature as we see that the details of low energy
eigenstates contributions matter and hence the ground state energy of each
flux configuration dominate the contribution to specific heat yielding a first
peak. However at very high temperature the low energy details of the spectrum
for the each flux configurations are dominated by the high energy part of the
eigenspectrum yielding a second peak where the gap between the minimum
energies of each flux configuration does not matter. This two gap structure of
the specific heat drew a lot of attention as the deciding criteria for the
spin-liquid phase or deconfined phase to be there. As we will observe that
depending on the relative strength of Kitaev interaction $K$, Heisenberg
interaction $J$, $\Gamma$ and external magnetic field $h_{z}$ we find
existence of well established two peak structure, a single sharp peak or a
wide broad peak. Primarily we assign the existence of clearly established flux
sectors and the eigenspectrum above it when there are two well established
peaks. On the other hand a single sharp peak would imply a strongly
interacting confining phase with no emergent flux structure to be present.
However the broad wide peak would either mean that two peak structure survive
but they are close enough and hence can not be separated or the excitations
are dominated by interacting magnon as found in recent studies.
Now we discuss how the specific heat behaves for various values of $J$ and
$\Gamma$ at zero magnetic field as shown in Fig. 13. When the non-Kitaev
interaction are present, the fluxes no longer remains conserved. The direct
implications of this is that the fermionic excitations can no linger be
grouped together by flux quantum number and the low lying excitations over
each flux sectors may not be separated by a gap and they may ultimately form a
continuous band which can be continuously extended to the fermionic
excitations at high energy. Thus we may lose the concept of two separate
energy scales in the problem and find only a single peat in the specific heat.
The presence of $J$ and $\Gamma$ destroys the above structure of the spectrum
and fluxes are no longer conserved. This implies any eigenstates can not be
identified with a given flux quantum number. However one must note that a
given energy eigenstate can always be written in linear combinations of
appropriate eigenstates of various flux quantum numbers. As the probability
amplitude will in general be unequal for different flux quantum states, one
can ask the question whether the notion of fermionic excitations and flux
excitations still holds in the presence of non-Kitaev interactions. In this
connections we have done extensive numerical analysis for all possible
combinations of signs and magnitude of Kitaev interaction (K), Heisenberg
interaction ($J$) and $\Gamma$ interactions and we find a great
controllability of excitation spectra as revealed by specific heat as
described below.

Figure 14: Normalized specific heat $C_{v}$ in the presence of magnetic field
is plotted. Each column denotes a particular cluster as mentioned at the top.
At the right of each row, we mention the value of $(K,J)$ taken. In the inset
$C_{v}$ for very large values of $J=\pm 2$ is shown. In the right the color
legends show the value of magnetic field taken.

General trends of specific heat : In the Kitaev limit i.e when $K=\pm
1,J=\Gamma=0$, we find that specific heat has double peak structure for all
the cluster except the four site cluster which has a weak shoulder like hump.
Also the nature of specific heat is identical for $K=\pm 1$ for $J=\Gamma=0$
for a given cluster. Presence of Heisenberg interaction acts differently
depending on whether it is ferromagnetic or antiferromagnetic. For
antiferromagnetic $J$ with $\gamma=0$, we find that only single peak in
specific heat exist but the position of the peaks can shift right or left
depending on the magnitude of $J$. On the other hand for ferromagnetic $J$,
after a critical $J_{c}$ we see that distinct appearances of two peak
structure in specific heat. This mainly happens for $|J|>K$, though depending
on the cluster size the two peak structure may appear or disappear as seen
comparing the specific heat for six site and eight site cluster for $J=-0.5$ .
On the other hand for ferromagnetic Kitaev interaction, the two peak structure
is more frequent. For all ferromagnetic Heisenberg interaction it is present
and for antiferromagnetic Heisenberg interaction the two peak structure occurs
very occasionally only, for example for six site cluster with $J=1.0$.
Effect of $\Gamma$ on specific heat : For $J=0$, we find that a ferromagnetic
$\Gamma$ stabilizes the two peak structure for $K=1$ and for $K=-1$, the
effect of $\Gamma$(be it ferro or antiferro) destroys the two peak structure.
However when the magnitude of $J$ becomes larger than the magnitude of
$\Gamma$, it is the Heisenberg interaction which controls the nature of
specific heat. Interestingly at certain cases the effect of $\Gamma$ is seen
to yield an additional third peak such as
$(N,K,J,\Gamma)=(6,-1,1,0),(6,-1,1,0.5),(6,1,0,0.5),(12,1,-2,0)$ This implies
the flux sectors further divided into two distinct energy scale in the
systems. To elucidate the effect of $\Gamma$ we refer
$(N,K,J,\Gamma)=(12,1,1.0,0.5)$ and $(12,-1,1,0.5)$ and we observe that in the
former there is a distinct two peak structure which is absent on the later.
Another universal feature is that for antiferromagnetic $J$ the peaks in the
specific heat is more pronounced. Possibly the more distinct effect of
$\Gamma$ can be seen for eight site cluster with $K=-1$. We see that for
$\Gamma=0.5$, only $J=1$ yields a two peak structure. For $\Gamma=0$ all
ferromagnetic $J$ (including $J=0$) has two peak structure but with
$\Gamma=-0.5$ the two peak structure for negative $J$ becomes asymmetric and
two peak structure reappear for positive $J=0.5,1.0$. Thus we can say that
$\Gamma$ can control the re-entrance of two peak structures in the specific
heat.

###  VIII.1 Effect of magnetic field on specific heat

Ferromagnetic ($K=-1$) vs antiferromagnetic Kitaev ($K=1$) for $J=0$: The
primary difference between ferromagnetic and antiferromagnetic Kitaev
interaction in the absence of Heisenberg interaction is that the two peaks in
the case of antiferromagnetic Kitaev interaction are more frequent for many
values of external magnetic field where for ferromagnetic Kitaev interaction
the two peak structure only appears for small magnetic field, see Fig. 14
first row for antiferromagnetic Kitaev interaction and second row for
ferromagnetic Kitaev interaction. This imply that the exact quantum
eigenstates are renormalized for $K=1$ and $K=-1$ differently under the
external magnetic field. The frequent occurrences of two peak structure
implies that the eigenstates in the presence of external magnetic field are
mildly modified in the case of antiferromagnetic Kitaev interaction. On the
other hand for ferromagnetic Kitaev interactions, the effect of external
magnetic field is more stronger to destroy the two peak structure implying
confined phases. The broader wide peaks in the case of antiferromagnetic
Kitaev interaction strengthen this understanding further that two peak
structure mostly survives denoting robustness of deconfined phase in the
presence of external magnetic field. With this background understanding we now
proceed to explain what happens when we take into consideration a finite
antiferromagnetic Heisenberg interaction.
Ferromagnetic ($K=-1$) vs antiferromagnetic Kitaev ($K=1$) for $J=0.5$: In the
case of antiferromagnetic Kitaev interaction the effect of having $J=0.5$ is
that it makes the two peak structure remain intact but with a qualitative
difference. If we compare Fig. 14 first row (with $J=0$) and third row (with
$J=0.5$), we observe that for $J=0.5$, the second peak is much more stronger
and the first peak is much weaker. Also the distance between the peaks has
increased for $J=0.5$. This denotes that antiferromagnetic Heisenberg
interaction stabilizes the eigenspectrum determined by the conserved flux-
sector qualitatively in the presence of magnetic field. Remarkably for
ferromagnetic Kitaev interaction also we see re-entrance of two peak structure
in the presence of antiferromagnetic Heisenberg interaction, see Fig. 14,
fourth row. We thus conclude that a antiferromagnetic Heisenberg interaction
helps to stabilize the Kitaev-spin-liquid phase in comparison to absence of
Heisenberg interaction in the presence of external magnetic field.
Ferromagnetic ($K=-1$) vs antiferromagnetic Kitaev ($K=1$) for $J=-0.5$: For
the ferromagnetic Heisenberg interaction we observe that two peak structure
are rarely occurring for antiferromagnetic Kitaev interaction in comparison to
$J=0.5$, compare Fig. 14 third row and fifth row. Though the exact pattern of
specific heat depends on the cluster sizes also. For example for all the
cluster except twelve site, there are mostly a single peak. Similarly for the
case of ferromagnetic Kitaev interaction we observe that there exists mostly
single peak. This tells us that a ferromagnetic Heisenberg interaction works
against the Kitaev-spin-liquid phase. It is interesting to note that an
antiferromagnetic Heisenberg interaction ,in general, induces a tendency to
have a quantum paramagnet state which seems to favor a Kitaev-spin-liquid
phase in comparison to a ferromagnetic Heisenberg interaction which induces a
collinear order in the system which seems to destabilize the Kitaev-spin-
liquid state.

Having described the effect of various competing interactions of Kitaev-
Heisenberg-$\Gamma$ model in finite temperature properties like
magnetization,susceptibility and specific heat, we now cover a brief
comparison with the existing results. The two peak structure of specific heat
has been vindicated as the hall mark signature of Kitaev spin liquid phase
both theoretically and experimentally kataoka-jpsj-2020 ; andrade-2020 ;
tathagata-2020 ; bachus-2021 ; balz-prb-2021 ; feng-2020 ; yogesh-2017 ;
nasu-2015 ; yamaji-2016 . The effect of magnetic field on these quantities
have been very interesting some of which resonate with our results. For
example it has been shown that the magnetic field help to stabilize the
magnetically ordered state before melting into a quantum paramagnetic state in
$\alpha-RuCl_{3}$balz-prb-2021  where Kitaev interaction and Heisenberg
interaction both has been considered ferromagnetic. In accordance we find from
the fourth row of Fig. 14 that when the $K$ and $J$ are both ferromagnetic
only for zero or very small magnetic field we find two peak structure and an
increased magnetic field stabilizes the one peak structure across all
clusters. Similarly susceptibility shows qualitatively similar behavior found
in Ref. andrade-2020,  though we have only shown the total susceptibilities.
In fact in accordance with the previous results andrade-2020 , we find in some
cases, the susceptibility having maximum at some intermediate magnetic field
and not at zero field, see Fig. VII, third, fourth and fifth row. However the
combinations of $(J,K)$ is different than those taken in andrade-2020,  where
$J$ and $K$ both are taken to be ferromagnetic and contain other complex
interactions, it is interesting that the small cluster is able to capture the
qualitative essencekataoka-jpsj-2020 . Similarly in another previous study
yamaji-2016  which considered a strong ferromagnetic Kitaev interaction and a
week antiferromagnetic Kitaev interaction in the presence of next-nearest
neighbor interaction, the specific heat is shown to change the width and
height of the two peaks depending on the external magnetic field. This
scenarios can be compared with Fig. 14, fourth row, where across all the
cluster similar behavior is observed. We comment that due to the limited
cluster size, the effect of next-nearest neighbor interaction is simulated
easily.

##  IX Discussion

We now summarize our extensive investigation of some key feature of Kitaev-
Hesineberg-$\Gamma$ model on finite size clusters. We started with the
inspection of eigenvalues of the system. It has been found that after a
critical values of magnitude of Heisenberg interaction the spectrum varies
linearly and conforms to the classical values for ferromagnetic Heisenberg
interaction or appropriate singlet states for antiferromagnetic Heisenberg
interaction. In the large ferromagnetic limit the energy density varies as
$E=-6J+2K$ which does not depends on the cluster size. However, very
interestingly, in the large anti-ferromagnetic Heisenberg interaction the
asymptotic energy density depends on the cluster sizes which implies that the
ground state is still a quantum ground state with large contribution of
various possible singlets which forms a linear combination due to Kitaev and
$\Gamma$ interactions. Near Kitaev limit, there are some branches of
eigenvalues which bend and do not change sign across zero as J is changed from
negative to positive. Also there are large number of eigenenergies whose
values vary linearly and changes sign as the sign of $J$ is reversed. Finally
we have investigated how the eigenvalues vary with external magnetic field and
we have seen remarkable closeness of the results found for large system.
Though for ferromagnetic Kitaev interaction the effect of magnetic field is
similar for all cluster, for antiferromagnetic $J$, the effect of magnetic
field is shown to depend largely on the cluster sizes. This follows from the
fact that the antiferromagnetic Heisenberg interaction induces a significant
quantum fluctuations which depend greatly on cluster sizes.
Next we investigated the bond-dependent correlation function on a given $z$
bond and looked at $<\sigma^{z}_{i}\sigma^{z}_{j}>$ and
$<\sigma^{y}_{i}\sigma^{y}_{j}>$ correlations. In the Kitaev limit only
$<\sigma^{z}_{i}\sigma^{z}_{j}>$ is non-zero on $z$-bond while $x-x$ and $y-y$
type correlation is zero. In the presence of other non-Kitaev interaction
$x-x$ and $y-y$ type correlation function becomes non-zero. On the other hand
in the Heisenberg limit, the bond-correlation function
$<\sigma^{\alpha}_{i}\sigma^{\alpha}_{j}>$ is isotropic implying it is
identical for $\alpha=x=y=z$. Thus while in the Kitaev limit the bond-
dependent correlation function represents the extreme anisotropic character,
the Heisenberg limit represents an isotropic character and we would expect a
smooth transition between this as we increase the strength of Heisenberg
interaction. We observe that when $h_{z}$, the external magnetic field is
absent this transition is indeed continuous. However in the presence of some
intermediate values of external magnetic field the correlation function shows
non-analytic and anisotropic character which manifests as presence of
plateaus. This plateaus are only present for antiferromagnetic Heisenberg
interaction implying that the different quantum ground state in the Kitaev
limit and antiferromagnetic Heisenberg limit does not smoothly change from one
another and there is discontinuity in it.
As mentioned earlier Kitaev model is known to contain two types of Majorana
fermions named matter fermions and gauge fermions. In the Kitaev limit, two
gauge fermions combine to constitute a conserved $Z_{2}$ flux operator. On the
otherhand matter fermions constitute a non-interacting hopping problem coupled
to this conserved $Z_{2}$ field. However in the presence of other non-Kitaev
interaction gauge fermions acquire dynamics. Various interesting consequences
arise due to this such as confinement-deconfinement transition which is
related to revival of spin-liquid phase in intermediate temperature. As we see
that time evolution of Majorana fermions and gauge fermions follows
oscillatory yet decaying pattern depending on the different interaction. There
are various interesting pattern of this oscillations but must significant is
their dependency on external magnetic field. External magnetic field is seen
to greatly modulate the amplitude of oscillations and also the average
magnitude of it. Depending on the strength of external magnetic field, it can
freeze the gauge field without having any oscillations showing remarkable
dependence of magnetic field. Interesting realization of quantum speed
limit(QSL) phenomena also appear where probability of Majorana fermion
vanishes from the initial value of unity. A detail analysis of this QSL and
their dependency on Heisenberg coupling and external magnetic field will be
followed in future study.
Finally we have provided a detail analysis of finite temperature
magnetization, susceptibility and specific heat for various parameters. It has
been shown explicitly what is the salient difference between ferromagnetic and
antiferromagnetic Kitaev interaction and what is effect of magnetic field and
Heisenberg interaction (ferromagnetic and antiferromagnetic). Magnetization
shows clear distinction between ferro and antiferromagnetic Kitaev
interaction. For low field and antiferromagnetic Kitaev interaction,
magnetization shows a dome like pattern which is absent for ferromagnetic
Kitaev interaction. Also the zero temperature magnetization have large
dependency on external magnetic field. It can be evenly spaced from zero to a
maximum or close to unity depending on whether Kitaev interaction is
antiferromagnetic or antiferromagnetic. The Heisenberg interaction is shown to
modulate the qualitative nature of the magnetization further showing distinct
quantum effects depending on whether it is antiferromagnetic or ferromagnetic.
Similar intriguing dependencies susceptibility and specific heat on various
exchange coupling and external magnetic field has been observed.
With the recent intense effort to understand Kitaev-Heisenberg-$\Gamma$ system
we think that our extensive analysis would serve as an useful reference.
However we mention that all the aspect found here may not be readily
generalized to thermodynamic system. Rather it would be prudent to consider
this as an interesting platform to quantify the quantum and thermal
fluctuations pertaining to Kitaev-Heisenberg-$\Gamma$ system and use it for
more better understanding of the thermodynamic system. For example in all the
cases, universally, we have found that ferromagnetic and antiferromagnetic
Kitaev interaction has distinct responses to external magnetic field and
Heisenberg interactions. It appears that in the case of antiferromagnetic
Heisenberg interaction, the ground state has more quantum correlations
irrespective of sign of Kitaev interaction which manifests itself in the
plateaus like behavior in the correlation function, dome like shape in
magnetization and stable two peak structure in specific heat. On the other
hand a ferromagnetic Heisenberg interaction is seen induce different
characteristics. Thus a proper quantification of change of quantum fluctuation
from the Kitaev limit to Heisenberg limit for both $J$ positive and negative
would provide valuable insight to Kitaev physics.
Kitaev-Heisenberg-$\Gamma$ system is shown to posses various collinear
magnetic order at zero temperature which is destabilized at intermediate
temperature towards a spin-liquid phase. A study of classical ground state on
Kitaev cluster shown to exhibit interesting order by singularity mechanism
within classical ground state space sarvesh-2020 ; anjana-2018 . An extension
of this study to larger clusters and in the presence of Heisenberg interaction
could provide valuable insight to explore new mechanism of quantum fluctuation
khatua-2021 . Further within the slave-fermion language implementing a PSG
calculationwen-2002 ; yasir-2020 ; atanu-2023  for this finite size system
utilizing various symmetries available and their response to various non-
Kitaev interaction could be useful to compare the exact analysis done here and
will be taken as future study. On a purely speculative way, we think that
implementing a quench dynamics for the Kitaev system using external magnetic
field or pressure could reveal interesting signatures of various phases
obtained in Kitaev-Heisenberg-$\Gamma$ system.

## Acknowledgement

SMP acknowledges SAMKHYA (High-Performance Computing Facility provided by the
Institute of Physics, Bhubaneswar) for the numerical computation. SM
acknowledges support from the ICTP through the Associate’s Programme
(2020-2025). SM also thanks G. Baskaran and Nicola Seriani for interesting
discussions.

## References

  * (1) A. Kitaev, Annals of Physics 321, (2006) 2-111.
  * (2) G. Baskaran, S. Mandal, and R. Shankar, Phys. Rev. Lett. 98, 247201 (2007).
  * (3) Wen-Han Kao, Natalia B. Perkins, Annals of Physics 435 (2021) 168506.
  * (4) Joji Nasu and Yukitoshi Motome, Phys. Rev. B 102, 054437 (2020).
  * (5) Tathagata Chowdhury, Achim Rosch, and Ralf Bulla, Phys. Rev. B 101, 115133(2020).
  * (6) K. S. Tikhonov, M. V. Feigelman, and A. Yu. Kitaev, Phys. Rev. Lett 106, 067203(2011).
  * (7) A. V. Lunkin, K. S. Tikhonov, M. V. Feigelman, Journal of Physics and Chemistry Solids, Volume 128, May 2019, Pages 130-137.
  * (8) Saptarshi Mandal, Subhro Bhattacharjee, K Sengupta, R Shankar, G Baskaran, Phys. Rev. B 84,155121 (2011).
  * (9) Animesh Nanda, Kusum Dhochak, and Subhro Bhattacharjee, Phys. Rev. B 102, 235124 (2020).
  * (10) Animesh Nanda, Adhip Agarwala, and Subhro Bhattacharjee, Phys. Rev. B 104, 195115 (2021).
  * (11) Johannes Knolle, Subhro Bhattacharjee, and Roderich Moessner, Phys. Rev. B 97, 134432 (2018).
  * (12) G. Baskaran, Diptiman Sen, and R. Shankar, Phys. Rev. B 78, 115116 (2008).
  * (13) Samarth Chandra, Kabir Ramola, and Deepak Dhar, Phys. Rev E 82, 031113(2010).
  * (14) Saptarshi Mandal and Naveen Surendran, Phys. Rev. B 79 024426 (2008).
  * (15) H. Yao and S. A. Kivelson, Phys. Rev. Lett. 99, 247203 (2007).
  * (16) G. Kells, D. Mehta, J. K. Slingerland, and J. Vala, Phys. Rev. B 81, 104429 (2010).
  * (17) Tim Eschmann, Petr A. Mishchenko, Kevin O’Brien, Troels A. Bojesen, Yasuyuki Kato, Maria Hermanns, Yukitoshi Motome, and Simon Trebst, Phys. Rev. B 102, 075125 (2020).
  * (18) Jiri Chaloupka, George Jackeli, and Giniyat Khaliullin Phys. Rev. Lett. 105, 027204(2010).
  * (19) Simon Trebst, Ciaran Hickey, Physics Reports 950 (2022)1-37.
  * (20) Hidenori Takagi, Tomohiro Takayama, George Jackeli, Giniyat Khaliullin & Stephen E. Nagler, Nature Reviews Physics 1, 264-280 (2019).
  * (21) Yukitoshi Motome, Ryoya Sano, Seonghoon Jang, Yusuke Sugita1 and Yasuyuki Kato, J. Phys.: Condens. Matter 32 404001 (2020).
  * (22) M. Hermanns, I. Kimchi, and J. Knolle, Annu. Rev. Condens. Matter Phys. 9, 17–33 (2018).
  * (23) A. Banerjee, C. A. Bridges, J.-Q. Yan, A. A. Aczel, L. Li, M. B. Stone, G. E. Granroth, M. D. Lumsden, Y. Yiu, J. Knolle, S. Bhattacharjee, D. L. Kovrizhin, R. Moessner, D. A. Tennant, D. G. Mandrus, and S. E. Nagler, Nature Materials 15, 733–740 (2016).
  * (24) Li Ern Chern, Finn Lasse Buessen and Yong Baek Kim, npj Quantum Materials 6, 33 (2016).
  * (25) Lukas Janssen, Eric C. Andrade, and Matthias Vojta, Phys. Rev. B 96, 064430 (2017).
  * (26) Peter Czajka, Tong Gao, Max Hirschberger, Paula Lampen-Kelley, Arnab Banerjee, Jiaqiang Yan, David G. Mandrus, Stephen E. Nagler and N. P. Ong, Nature Physics 17, 915-919 (2021).
  * (27) Niravkumar D. Patel and Nandini Trivedi, PNAS 116 (25), 12199-12203 (2019).
  * (28) Haoxiang Li, T. T. Zhang, A. Said, G. Fabbris, D. G. Mazzone, J. Q. Mandrus, Gabor B. Halasz, S. Okamoto, S. Murakami, M. P. M. Dean, H. N. Lee & H Miao, Nature Communications 12, 3513 (2012).
  * (29) C. Balz, L. Janssen, P. Lampen-Kelley, A. Banerjee, Y. H. Liu, J.-Q. Yan, D. G. Mandrus, M. Vojta, and S. E. Nagler, Phys. Rev. B 103, 174417 (2021).
  * (30) Christoph Berke, Simon Trebst, and Ciarán Hickey, Phys. Rev. B 101, 214442 (2020).
  * (31) D. Reig-i-Plessis, T. A. Johnson, K. Lu, Q. Chen, J. P. C. Ruff, M. H. Upton, T. J. Williams, S. Calder, H. D. Zhou, J. P. Clancy, A. A. Aczel, and G. J. MacDougall, Phys. Rev. Mat. 4, 124407 (2020).
  * (32) Dirk Wulferding, Youngsu Choi, Seung-Hwan Do, Chan Hyeon Lee, Peter Lemmens, Clement Faugeras, Yann Gallais & Kwang-Yong Choi, Nature Communications 11 1603 (2020).
  * (33) Sananda Biswas, Ying Li, Stephen M. Winter, Johannes Knolle, and Roser Valenti, Phys. Rev. Lett 123, 237201 (2019).
  * (34) V. Leeb, K. Polyudov, S. Mashhadi, S. Biswas, Roser Valenti, M. Burghard, and J. Knolle, Phys. Rev. Lett, 126, 097201 (2021).
  * (35) P. W. Anderson, Science, 177, 393 (1972).
  * (36) Steven Strogatz, Sara Walker, Julia M. yeomans, Corina tarnita, Elsa Arcaute, Manlio De Domenico, Oriol Artime and Kwang-II Goh, Nature Reviews Physics 4, 508-510 (2022).
  * (37) Lucile Savary and Leon Balents, Rep. Prog. Phys. 80 016502 (2017).
  * (38) Kota Kataoka, Daigoro Hirai, Takeshi Yajima, Daisuke Nishio-Hamane, Rieko Ishii, Kwang-Yong Choi, Dirk Wulferding, Peter Lemmens, Shunichiro Kittaka,Toshiro Sakakibara, Hajime Ishikawa, Akira Matsuo, Koichi Kindo, and Zenji Hiroi, J. Phys. Soc. Jpn. 89, 114709 (2020).
  * (39) Yiping Wang, Gavin B. Osterhoudt, Yao Tian, Paige Lampen-Kelley, Arnab Banerjee, Thomas Goldstein, Jun Yan, Johannes Knolle, Huiwen Ji, Robert J. Cava, Joji Nasu, Yukitoshi Motome, Stephen E. Nagler, David Mandrus & Kenneth S. Burch, npj Quantum Materials 5, 14 (2020).
  * (40) Shi Wang, Zhongyuan Qi, Bin Xi, Wei Wang, Shun-Li Yu, and Jian-Xin Li, Phys. Rev. B 103, 054410 (2021).
  * (41) Joji Nasu, Masafumi Udagawa, and Yukitoshi Motome, Phys. Rev. B 92, 115122 (2015).
  * (42) Akihisa Koga and Joji Nasu, Phys. Rev. B 100, 100404(R) (2019).
  * (43) Han Li, Hao-Kai Zhang, Jiucai Wang, Han-Qing Wu, Yuan Gao, Dai-Wei Qu, Zheng-Xin Liu, Shou-Shu Gong $\&$ Wei Li, Nature Communications v12, 4007 (2021).
  * (44) Zheng Zhu, Itamar Kimchi, D. N. Sheng and Liang Fu, Phys. Rev. B 97, 241110(R) (2018).
  * (45) Markus Holzmann, Raymond C. Clay III, Miguel A. Morales, Norm M. Tubman, David M. Ceperley, and Carlo Pierleoni, Phys. Rev B 94, 035126 (2016).
  * (46) Han Li, Wei Li and Gang Su, Phys. Rev. B 107, 115124 (2023).
  * (47) F. Yilmaz, A. P. Kampf, and S. K. Yip, Phys. Rev. Research 4, 043024 (2022).
  * (48) Kyusung Hwang, Ara Go, Ji Heon Seong, Takasada Shibauchi and Eun-Gook Moon, Nature Communications 13, 323 (2022).
  * (49) Y. Yamaji, T. Suzuki, T. Yamada, S. I. Suga, N. Kawashima, and M. Imada, Phys. Rev. B 93, 174425 (2016).
  * (50) Shi Feng, Adhip Agarwala, Subhro Bhattacharjee, Nandini Trivedi, condmat, arXiv:2206.12990v2
  * (51) Shuyi Li, Han Yan, and Andriy H. Nevidomskyy, condmat arXiv:2301.07401v2
  * (52) R. L. Smit, S. Keupert, O. Tsyplyatyev, P. A. Maksimov, A. L. Chernyshev, and P. Kopietz, Phys. Rev. B 101, 054424 (2020).
  * (53) Itamar Kimchi, James G. Analytis, and Ashvin Vishwanath, Phys. Rev. B 90, 205126 (2014).
  * (54) Eric Kin-Ho Lee and Yong Baek Kim, Phys. Rev. B 91, 064407(2015).
  * (55) Xue-Yang Song, Yi-Zhuang You, and Leon Balents, Phys. Rev. Lett 117, 037209 (2016).
  * (56) Fan Yang, Kirill Plekhanov, and Karyn Le Hur, Phys. Rev. Research 2, 013005(2020).
  * (57) E. H. Lieb, Phys. Rev. Lett. 73, 2158 (1994).
  * (58) David C. Ronquillo, Adu Vengal and Nandini Trivedi, Phys. Rev. B 99, 140413(R) (2019).
  * (59) Shuang Liang, Bo-Sen He, Zhao-Yang Dong, Wei Chen, Jian-Xin Li, and Qiang-Hua Wang, Phys. Rev. B 98, 104410 (2018).
  * (60) Saptarshi Mandal, R Shankar, G Baskaran, Journal of Physics A: Math. and Theor. 45, 335304 (2012).
  * (61) Joji Nasu and Yukitoshi Motome, Phys. Rev. Research 1, 033007 (2019).
  * (62) Masahiko G. Yamada, npj Quantum Materials 5, 82 (2020).
  * (63) Sebastian Deffner and Steve Campbell, J. Phys. A: Math. Theor. 50, 453001 (2017).
  * (64) Gal Ness, Andrea Alberti and Yoav Sagi, Phys. Rev. Lett. 129, 140403 (2022).
  * (65) Brij Mohan and Arun Kumar Pati, Phys. Rev. A 106, 042436 (2022).
  * (66) Srishty Aggarwal, Subhashish Banerjee, Arindam Ghosh and Banibrata Mukhopadhyay, New J. Phys. 24, 085001 (2022).
  * (67) S. Bachus, D. A. S. Kaib, A. Jesche, V. Tsurkan, A. Loidl, S. M. Winter, A. A. Tsirlin, R. Valenti, and P. Gegenwart, Phys. Rev. B 103, 054440 (2021).
  * (68) Eric C. Andrade, Lukas Janssen, and Matthias Vojta, Phys. Rev. B 102, 115160 (2020).
  * (69) Kexin Feng, Natalia B. Perkins , and F. J. Burnell, Phys. Rev. B 102, 224402 (2020).
  * (70) Kavita Mehlawat, A. Thamizhavel, and Yogesh Singh, Phys. Rev. B 95, 144406 (2017).
  * (71) J. Nasu, M. Udagawa and Y. Motome, Phys. Rev. B 92, 115122(2015).
  * (72) Anjana M. Samarakoon, Gideon Wachtel, Youhei Yamaji, D. A. Tennant, Cristian D. Batista, and Yong Baek Kim, Phys. Rev. B 98, 045121 (2018).
  * (73) Sarvesh Srinivasan, Subhankar Khatua,G. Baskaran, and R. Ganesh, Phys. Rev. Research 2, 023212(2020).
  * (74) Subhankar Khatua, Sarvesh Srinivasan, and R. Ganesh, Phys. Rev. B 103, 174412 (2021).
  * (75) X.-G. Wen, Phys. Rev. B 65, 165113 (2002).
  * (76) J. Sonnenschein, A. Chauhan, Y. Iqbal, and J. Reuther, Phys. Rev. B 102, 125140 (2020).
  * (77) Atanu Maity, Francesco Ferrari, Ronny Thomale, Saptarshi Mandal, and Yasir Iqbal, Phys. Rev. B 107, 134438 (2023).

##  X Appendix

##  Appendix A Eigenstates of Kitaev-Heisenberg 4-site cluster

In 4-site cluster, with only Kitaev interaction the Hamiltonian is easily
diagonalized and the eigenvalues are $\\{\pm 2K_{x},\pm 2K_{y},\pm 2K_{z},\pm
2\sqrt{K_{x}^{2}+K_{y}^{2}+K_{z}^{2}}\\}$ with each energy level having two-
fold degeneracy. The Hamiltonian can be diagonalized into two sectors one with
even number of spin up and the other one with odd number of up. Two degenerate
ground states in Kitaev limit, are:

$\displaystyle\ket{{\rm GS}_{1}}$ $\displaystyle=$ $\displaystyle
K_{y}\left(\ket{\uparrow\uparrow\uparrow\uparrow}+\ket{\downarrow\downarrow\downarrow\downarrow}\right)-K_{x}\left(\ket{\uparrow\downarrow\uparrow\downarrow}+\ket{\downarrow\uparrow\downarrow\uparrow}\right)$
$\displaystyle+\left(K_{z}+\tilde{K}\right)\left(\ket{\uparrow\downarrow\downarrow\uparrow}+\ket{\downarrow\uparrow\uparrow\downarrow}\right),$
$\displaystyle\ket{{\rm GS}_{2}}$ $\displaystyle=$ $\displaystyle
K_{y}(\ket{\uparrow\uparrow\downarrow\downarrow}+\ket{\downarrow\downarrow\uparrow\uparrow})+K_{x}(\ket{\uparrow\downarrow\downarrow\uparrow}+\ket{\downarrow\uparrow\uparrow\downarrow})$
$\displaystyle+(K_{z}-\tilde{K})(\ket{\uparrow\downarrow\uparrow\downarrow}+\ket{\downarrow\uparrow\downarrow\uparrow}),$
$\displaystyle{\rm where}$
$\displaystyle\tilde{K}=\sqrt{K_{x}^{2}+K_{y}^{2}+K_{z}^{2}}.$ (3)

The spin-spin correlations (CR) can also be written in closed analytical form.
The non vanishing correlations are:

$\displaystyle{\rm CR}^{x-x}$ $\displaystyle=$
$\displaystyle-4\frac{K_{x}}{K_{y}^{2}}(K_{z}+\sqrt{K_{x}^{2}+K_{y}^{2}+K_{z}^{2}})/N_{0};$
$\displaystyle{\rm CR}^{y-y}$ $\displaystyle=$
$\displaystyle-4\frac{1}{K_{y}}(K_{z}+\sqrt{K_{x}^{2}+K_{y}^{2}+K_{z}^{2}})/N_{0};$
$\displaystyle{\rm CR}^{z-z}$ $\displaystyle=$
$\displaystyle-4\frac{K_{z}}{K_{y}^{2}}(K_{z}+\sqrt{K_{x}^{2}+K_{y}^{2}+K_{z}^{2}})/N_{0};$
$\displaystyle N_{0}$ $\displaystyle=$
$\displaystyle\sqrt{2+2\frac{K_{x}^{2}}{K_{y}^{2}}+\frac{2}{K_{y}^{2}}(K_{z}+\sqrt{K_{x}^{2}+K_{y}^{2}+K_{z}^{2}})^{2}}.$

With the application of Heisenberg interaction $J$ the Hamiltonian looses some
of the symmetries and the solution is more complex. The eigenvalues of the
odd-particle sector are $-2J\pm 2K_{y},-2J-2K_{y},-2J\pm
2K_{x},-2J-2K_{x},6J+2K_{y},6J+2K_{x}$ and the even-particle sector has
eigenvalues $6J+2K_{z},-2J\pm 2K_{z},-2J-2K_{z},\pm
2\left(6J^{2}+\sum_{\alpha}(J+K_{\alpha})^{2}\pm
J\sqrt{8\sum_{\alpha}(K_{\beta}-K_{\gamma})^{2}}\right)^{1/2}$, where each of
$\alpha,\beta,\gamma$ can take values $x,y,z$ and $\alpha\neq\beta\neq\gamma$.

In the presence of only Heisenberg interaction $J$ (i.e., $K$=0), the cluster
is SU(2) symmetric. In this limit, the ground state is a singlet state of the
form
$\frac{1}{2}\left(\ket{\uparrow_{1}\downarrow_{2}}-\ket{\downarrow_{1}\uparrow_{2}}\right)\times\left(\ket{\uparrow_{3}\downarrow_{4}}-\ket{\downarrow_{3}\uparrow_{4}}\right)$.
When we introduce external magnetic field $h_{z}$, we find presence of
$\ket{\uparrow\uparrow\downarrow\downarrow},\ket{\downarrow\downarrow\uparrow\uparrow}$
states in the ground state with some amplitude. When we increase the magnitude
of $h_{z}$, this amplitude changes but only these six states make the ground
state. For ferromagnetic $J$, even a small magnetic field yields a large
contribution of $\ket{\uparrow\uparrow\uparrow\uparrow}$ to the ground state.
The 4-site cluster in the presence of Kitaev and Heisenberg interaction, can
be solved analytically. The Hamiltonian becomes block-diagonalized of two
blocks for even and odd particle sector. The basis states used for these two
sectors are:

$\displaystyle\ket{\rm
even}=\begin{pmatrix}\ket{\uparrow\uparrow\uparrow\uparrow}\\\
\ket{\uparrow\downarrow\downarrow\uparrow}\\\
\ket{\downarrow\uparrow\uparrow\downarrow}\\\
\ket{\uparrow\uparrow\downarrow\downarrow}\\\
\ket{\downarrow\downarrow\uparrow\uparrow}\\\
\ket{\uparrow\downarrow\uparrow\downarrow}\\\
\ket{\downarrow\uparrow\downarrow\uparrow}\\\
\ket{\downarrow\downarrow\downarrow\downarrow}\end{pmatrix};$
$\displaystyle\ket{\rm
odd}=\begin{pmatrix}\ket{\downarrow\uparrow\uparrow\uparrow}\\\
\ket{\uparrow\downarrow\uparrow\uparrow}\\\
\ket{\uparrow\uparrow\downarrow\uparrow}\\\
\ket{\uparrow\uparrow\uparrow\downarrow}\\\
\ket{\uparrow\downarrow\downarrow\downarrow}\\\
\ket{\downarrow\uparrow\downarrow\downarrow}\\\
\ket{\downarrow\downarrow\uparrow\downarrow}\\\
\ket{\downarrow\downarrow\downarrow\uparrow}\end{pmatrix};$ (5)

and the odd particle sector is diagonalized in the rotated basis:

$\displaystyle\frac{1}{\sqrt{8}}\begin{pmatrix}1&-1&-1&1&-1&1&1&-1\\\
1&1&1&1&-1&-1&-1&-1\\\ 0&\sqrt{2}&-\sqrt{2}&0&0&-\sqrt{2}&\sqrt{2}&0\\\
\sqrt{2}&0&0&-\sqrt{2}&-\sqrt{2}&0&0&\sqrt{2}\\\
-\sqrt{2}&\sqrt{2}&0&0&-\sqrt{2}&\sqrt{2}&0&0\\\ 1&1&-1&-1&1&1&-1&-1\\\
1&1&1&1&1&1&1&1\\\ 0&0&\sqrt{2}&-\sqrt{2}&0&0&\sqrt{2}&-\sqrt{2}\\\
\end{pmatrix}\begin{pmatrix}\ket{\downarrow\uparrow\uparrow\uparrow}\\\
\ket{\uparrow\downarrow\uparrow\uparrow}\\\
\ket{\uparrow\uparrow\downarrow\uparrow}\\\
\ket{\uparrow\uparrow\uparrow\downarrow}\\\
\ket{\uparrow\downarrow\downarrow\downarrow}\\\
\ket{\downarrow\uparrow\downarrow\downarrow}\\\
\ket{\downarrow\downarrow\uparrow\downarrow}\\\
\ket{\downarrow\downarrow\downarrow\uparrow}\end{pmatrix}.$ (6)

The corresponding eigenvalues are
$-2J+2K_{y},6J+2K_{y},-2J-2K_{y},-2J-2K_{y},-2J-2K_{x},-2J+2K_{x},6J+2K_{x},-2J-2K_{x}$.
When the even-sector basis is rotated to a new one:

$\displaystyle\begin{pmatrix}\ket{\uparrow\uparrow\uparrow\uparrow}\\\
\ket{\uparrow\downarrow\downarrow\uparrow}\\\
\ket{\downarrow\uparrow\uparrow\downarrow}\\\
\ket{\uparrow\uparrow\downarrow\downarrow}\\\
\ket{\downarrow\downarrow\uparrow\uparrow}\\\
\ket{\uparrow\downarrow\uparrow\downarrow}\\\
\ket{\downarrow\uparrow\downarrow\uparrow}\\\
\ket{\downarrow\downarrow\downarrow\downarrow}\end{pmatrix}$
$\displaystyle\rightarrow$
$\displaystyle\frac{1}{\sqrt{2}}\begin{pmatrix}(1+\Delta)\ket{\uparrow\uparrow\uparrow\uparrow}\\\
(1+\Delta)\ket{\uparrow\downarrow\downarrow\uparrow}\\\
(1+\Delta)\ket{\uparrow\uparrow\downarrow\downarrow}\\\
(1+\Delta)\ket{\uparrow\downarrow\uparrow\downarrow}\\\
(1-\Delta)\ket{\uparrow\uparrow\uparrow\uparrow}\\\
(1-\Delta)\ket{\uparrow\downarrow\downarrow\uparrow}\\\
(1-\Delta)\ket{\uparrow\uparrow\downarrow\downarrow}\\\
(1-\Delta)\ket{\uparrow\downarrow\uparrow\downarrow}\\\ \end{pmatrix},$ (7)

where $\Delta=\sigma_{1}^{x}\sigma_{2}^{x}\sigma_{3}^{x}\sigma_{4}^{x}$, the
even-sector of the Hamiltonian is block diagonalized in two blocks. The lower
four are eigenstates (which are part of the second of the two blocks) with
eigenvalues $6J+2K_{z},-2J-2K_{z},-2J-2K_{z},-2J+2K_{z}$. The upper four basis
states reduce the Hamiltonian to a $4\times 4$ block:

$\begin{pmatrix}6J+2K_{z}&-2K_{y}&2K_{x}&0\\\
-2K_{y}&-2J-2K_{z}&4J&4J+2K_{x}\\\ 2K_{x}&4J&-2J-2K_{z}&4J+2K_{y}\\\
0&4J+2K_{x}&4J+2K_{y}&-2J+2K_{z}\end{pmatrix},$ (8)

and this gets diagonalized in the basis of,

$\displaystyle\frac{(1+\Delta)}{\sqrt{2}}\left(\upsilon_{1}^{(1)}\ket{\uparrow\uparrow\uparrow\uparrow}+\upsilon_{1}^{(2)}\ket{\uparrow\downarrow\downarrow\uparrow}+\upsilon_{1}^{(3)}\ket{\uparrow\uparrow\downarrow\downarrow}+\upsilon_{1}^{(4)}\ket{\uparrow\downarrow\uparrow\downarrow}\right)$
$\displaystyle\frac{(1+\Delta)}{\sqrt{2}}\left(\upsilon_{2}^{(1)}\ket{\uparrow\uparrow\uparrow\uparrow}+\upsilon_{2}^{(2)}\ket{\uparrow\downarrow\downarrow\uparrow}+\upsilon_{2}^{(3)}\ket{\uparrow\uparrow\downarrow\downarrow}+\upsilon_{2}^{(4)}\ket{\uparrow\downarrow\uparrow\downarrow}\right)$
$\displaystyle\frac{(1+\Delta)}{\sqrt{2}}\left(\upsilon_{3}^{(1)}\ket{\uparrow\uparrow\uparrow\uparrow}+\upsilon_{3}^{(2)}\ket{\uparrow\downarrow\downarrow\uparrow}+\upsilon_{3}^{(3)}\ket{\uparrow\uparrow\downarrow\downarrow}+\upsilon_{3}^{(4)}\ket{\uparrow\downarrow\uparrow\downarrow}\right)$
$\displaystyle\frac{(1+\Delta)}{\sqrt{2}}\left(\upsilon_{4}^{(1)}\ket{\uparrow\uparrow\uparrow\uparrow}+\upsilon_{4}^{(2)}\ket{\uparrow\downarrow\downarrow\uparrow}+\upsilon_{4}^{(3)}\ket{\uparrow\uparrow\downarrow\downarrow}+\upsilon_{4}^{(4)}\ket{\uparrow\downarrow\uparrow\downarrow}\right),$
$\displaystyle\upsilon_{1}$ $\displaystyle=$
$\displaystyle\bot_{1}\begin{pmatrix}\left(-3J^{2}(K_{\Sigma}-3K_{z})-2JK_{R}(K_{z}-\lambda_{-})-J\left(-K_{R}^{2}-K_{z}^{2}+K_{x}K_{y}+6JK_{R}\right)+J\lambda_{-}(K_{\Sigma}-3K_{z})\right)\\\
3J^{2}(K_{x}-K_{y})+2JK_{y}K_{R}+J\left(K_{R}^{2}-K_{z}^{2}+K_{x}K_{y}\right)-J\lambda_{-}(K_{x}-K_{y})\\\
3J^{2}(K_{x}-K_{y})-2JK_{x}K_{R}-J\left(K_{R}^{2}-K_{z}^{2}+K_{x}K_{y}\right)-J\lambda_{-}(K_{x}-K_{y})\\\
J(K_{x}-K_{y})(3J+K_{\Sigma}-\lambda_{-})\end{pmatrix}$ (10)
$\displaystyle\upsilon_{2}$ $\displaystyle=$
$\displaystyle\bot_{2}\begin{pmatrix}\left(-3J^{2}(K_{\Sigma}-3K_{z})-2JK_{R}(K_{z}+\lambda_{-})-J\left(-K_{R}^{2}-K_{z}^{2}+K_{x}K_{y}+6JK_{R}\right)-J\lambda_{-}(K_{\Sigma}-3K_{z})\right)\\\
3J^{2}(K_{x}-K_{y})+2JK_{y}K_{R}+J\left(K_{R}^{2}-K_{z}^{2}+K_{x}K_{y}\right)+J\lambda_{-}(K_{x}-K_{y})\\\
3J^{2}(K_{x}-K_{y})-2JK_{x}K_{R}-J\left(K_{R}^{2}-K_{z}^{2}+K_{x}K_{y}\right)+J\lambda_{-}(K_{x}-K_{y})\\\
J(K_{x}-K_{y})(3J+K_{\Sigma}+\lambda_{-})\\\ \end{pmatrix}$ (11)
$\displaystyle\upsilon_{3}$ $\displaystyle=$
$\displaystyle\bot_{3}\begin{pmatrix}\left(-3J^{2}(K_{\Sigma}-3K_{z})+2JK_{R}(K_{z}-\lambda_{+})+J\left(K_{R}^{2}+K_{z}^{2}-K_{x}K_{y}+6JK_{R}\right)+J\lambda_{+}(K_{\Sigma}-3K_{z})\right)\\\
3J^{2}(K_{x}-K_{y})-2JK_{y}K_{R}+J\left(K_{R}^{2}-K_{z}^{2}+K_{x}K_{y}\right)-J\lambda_{+}(K_{x}-K_{y})\\\
3J^{2}(K_{x}-K_{y})+2JK_{x}K_{R}-J\left(K_{R}^{2}-K_{z}^{2}+K_{x}K_{y}\right)-J\lambda_{+}(K_{x}-K_{y})\\\
J(K_{x}-K_{y})(3J+K_{\Sigma}-\lambda_{+})\\\ \end{pmatrix}$ (12)
$\displaystyle\upsilon_{4}$ $\displaystyle=$
$\displaystyle\bot_{4}\begin{pmatrix}\left(-3J^{2}(K_{\Sigma}-3K_{z})+2JK_{R}(K_{z}+\lambda_{+})+J\left(K_{R}^{2}+K_{z}^{2}-K_{x}K_{y}+6JK_{R}\right)-J\lambda_{+}(K_{\Sigma}-3K_{z})\right)\\\
3J^{2}(K_{x}-K_{y})-2JK_{y}K_{R}+J\left(K_{R}^{2}-K_{z}^{2}+K_{x}K_{y}\right)+J\lambda_{+}(K_{x}-K_{y})\\\
3J^{2}(K_{x}-K_{y})+2JK_{x}K_{R}-J\left(K_{R}^{2}-K_{z}^{2}+K_{x}K_{y}\right)+J\lambda_{+}(K_{x}-K_{y})\\\
J(K_{x}-K_{y})(3J+K_{\Sigma}+\lambda_{+})\\\ \end{pmatrix}$ (13)

with proper normalization $\bot_{1},\bot_{2},\bot_{3},\bot_{4}$.
$\upsilon_{i}^{(j)}$ is the $j^{\rm th}$ component of $\upsilon_{i}$. Here we
have used the following abbreviations:

$\displaystyle K_{\Sigma}$ $\displaystyle=$ $\displaystyle K_{x}+K_{y}+K_{z};$
(14) $\displaystyle\tilde{K}$ $\displaystyle=$
$\displaystyle\sqrt{K_{x}^{2}+K_{y}^{2}+K_{z}^{2}};$ (15) $\displaystyle
K_{R}$ $\displaystyle=$
$\displaystyle\sqrt{\tilde{K}^{2}-K_{x}K_{y}-K_{y}K_{z}-K_{z}K_{x}};$ (16)
$\displaystyle\lambda_{\pm}$ $\displaystyle=$
$\displaystyle\sqrt{9J^{2}+\tilde{K}^{2}+2JK_{\Sigma}\pm 4JK_{R}}.$ (17)

In this basis, the $4\times 4$ block is diagonalized with eigenvalues
$-2\lambda_{-},2\lambda_{-},-2\lambda_{+},2\lambda_{+}$.