Interlayer Coupling Effect on a Bilayer Kitaev Model
Hiroyuki Tomishige Joji Nasu Akihisa Koga
cond-mat.str-el cond-mat.mtrl-sci








Department of Physics, Tokyo Institute of Technology, Meguro, Tokyo 152-8551, Japan

Department of Physics, Tokyo Institute of Technology, Meguro, Tokyo 152-8551, Japan

Department of Physics, Tokyo Institute of Technology, Meguro, Tokyo 152-8551, Japan


 

We investigate a bilayer Kitaev model, where two honeycomb layers are coupled by the Heisenberg interactions, to discuss effects of an interlayer coupling against the Kitaev quantum spin liquids (QSLs).
In this model, there exists a local conserved quantity, which results in no long-range spin correlations in the system.
Using the exact diagonalization, bond operator mean-field theory, and cluster expansion techniques, we study ground state properties in the system.
The obtained results suggest the existence of a first-order quantum phase transition between the Kitaev QSL and singlet-dimer states.
We find that one-triplet excitation from the singlet-dimer ground state is localized owing to the existence of the local conserved quantity.
To examine finite-temperature properties, we make use of the thermal pure quantum state approach.
We clarify that double-peak structure in the specific heat inherent in the Kitaev QSL
is maintained even above the quantum phase transition.
The present results suggest that the Kitaev QSL is stable against the interlayer interference.
Magnetic properties of multilayer Kitaev models are also addressed.





Interlayer Coupling Effect on a Bilayer Kitaev Model

    Akihisa Koga
    December 30, 2023
======================================================




§ INTRODUCTION



Exploring quantum spin liquids (QSLs) is one of the central subjects in condensed matter physics since the Anderson's suggestion <cit.>.
A lot of theoretical and experimental studies have been devoted to clarify the nature of QSLs but the thermodynamic properties and excitation spectra remain elusive.
Theoretically, it has been still difficult to analyze the properties of frustrated Heisenberg models, which are considered to be archetypal models of QSLs, without approximations.
On the other hand, one of the promising models to discuss QSLs is a quantum spin model on a honeycomb lattice with bond-dependent Ising interactions, which is known as the Kitaev model <cit.>.
This model is exactly solvable and its ground state is a QSL
with short-range spin correlations <cit.>.
Furthermore, quantum spins are fractionalized into
itinerant Majorana fermions and localized Z_2 fluxes.
This leads to gapless elementary excitations in the ground state and
double peak structure in the specific heat <cit.>.

Moreover, Jackeli and Khaliullin have revealed that
the Kitaev model should be realizable in Mott insulators
with the strong spin-orbit coupling and specific lattice structure,
where localized  j_ eff=1/2 spins are coupled
by the superexchange interactions <cit.>.
As the candidate materials, A_2IrO_3(A=Na, Li) <cit.>,
α-RuCl_3 <cit.> and H_3LiIr_2O_6 <cit.> have been intensively studied,
which stimulate further experimental and theoretical investigations
on the nature inherent in the Kitaev physics <cit.>, and the generalization of the Kitaev model <cit.>.


Nevertheless, it is also recognized that
the magnetic properties at low temperatures
in the candidate materials cannot be fully reproduced by
the two-dimensional Kitaev model
while this model should capture the magnetism at higher temperatures.
For example, the materials exhibit a long-range magnetic order
at low temperatures <cit.> and star-shape low-energy structure has been observed
by the inelastic neutron scattering experiments <cit.>.
To account for these features, a lot of theoretical investigations have been made for additional effects beyond the pure Kitaev model, such as Heisenberg and/or Γ terms with/without long-range
interactions <cit.>
and intermediate spin-orbit coupling <cit.>.
The previous works suggest that the Kitaev QSL survives
even in the small additional interactions owing to the existence of the spin gap
in the Kitaev model <cit.>, and the magnetic order observed in the real materials
is stable by the additional interactions within a two-dimensional honeycomb layer
in a plausible parameter range <cit.>.
On the other hand, the candidate materials
are composed of honeycomb layers,
and the role of the stacking structure
in α-RuCl_3 and H_3LiIr_2O_6
has been examined recently <cit.>.
Therefore, it is highly desired to study the stability of the Kitaev QSL
in the presence of an interlayer coupling.

In this paper, we investigate the bilayer Kitaev model,
where two honeycomb layers are coupled by the Heisenberg interaction,
as one of simple models [see Fig. <ref>(a)].
We show that this model possesses a local conserved quantity
on each pair of stacked two hexagons (bi-hexagon),
which leads to the absence of the long-range magnetic order.
To study the stability of the Kitaev QSL against the interlayer coupling,
we employ the bond-operator mean-field (MF)
approximation <cit.> and
cluster expansion technique <cit.>.
We also use the exact diagonalization (ED) on finite clusters
to discuss ground state properties and excitation spectrum.
The numerical results suggest the existence of
a first-order phase transition between the Kitaev QSL and
spin-singlet dimer states at zero temperature.
Moreover, using the thermal pure quantum (TPQ) state
approach <cit.>,
we clarity that the double-peak structure emerges
in the small interlayer coupling and retains
even above the quantum transition point expected from the ED calculations.
We also extend the argument based on the presence of
the local conserved quantities to the multilayer Kitaev model
and show the absence of three-dimensional long-range correlations in the multilayer Kitaev model with arbitrary stacking numbers as well.

This paper is organized as follows.
In Sec. <ref>, we introduce
the model Hamiltonian on the bilayer honeycomb lattice
and discuss the local conserved quantity and
the parity symmetry in the system.
Our methods are briefly summarized in Sec. <ref>.
In Sec. <ref>,
we show the numerical results for the ground state properties in the bilayer system, which suggest the existence of a first-order quantum phase transition
between the QSL and singlet dimer states.
Thermodynamic properties are discussed in Sec. <ref>.
In Sec. <ref>,
we discuss magnetic properties in the multilayer Kitaev model.
The summary is provided in the last section.





§ MODEL AND ITS SYMMETRY









To address the effect of the interlayer coupling between Kitaev models
on honeycomb lattices, we introduce a following simple model,
where two Kitaev models are coupled by the Heisenberg interaction
[see Fig. <ref>(a)]:

    H = -  J_K∑_ij_α,n S_i,n^α S_j,n^α
     + J_H ∑_i S_i,1· S_i,2,

where S_i,n^α=1/2σ_i,n^α (α=x,y,z) and
σ_i,n^α is Pauli matrix at site i of the n(=1,2)th layer.
J_K (>0) is the ferromagnetic Kitaev coupling in each layer
and J_H (>0) is the antiferromagnetic Heisenberg coupling between two layers.
We assume that each site on the layer 1 is located just above that
on the other.
In each layer, the anisotropy of the Ising-type interactions
depend on the bonds;
there are three kinds of nearest neighbor (NN) bonds,
⟨ ij ⟩_α (α=x,y,z),
which we refer to as the α bond, on the honeycomb lattice
[see Figs. <ref>(a) and <ref>(b)].

When J_H=0, the system is reduced to two single-layer Kitaev models.
In the model, there exists the local Z_2 conserved quantity
W_p,n=σ_p_1,n^x σ_p_2,n^y σ_p_3,n^z
σ_p_4,n^x σ_p_5,n^y σ_p_6,n^z on the nth layer ,
which results in the QSL ground state with long-range spin entanglement and the fractionalization of the quantum spins.
This allows us to map the single layer model onto free Majorana fermion system
with gapless elementary excitations
although there exists a gap in the spin excitation.
At finite temperatures, the spin fractionalization emerges as a peculiar temperature dependence of observables such as double-peak structure in the specific heat <cit.>.
On the other hand, in the case of J_K=0,
the system is composed of independent dimers.
The ground state is represented by the direct product of
interlayer dimer singlets with the spin gap, where inter-dimer wave functions are disentangled.
Although these two nonmagnetic ground states possess a spin gap,
their low-energy properties are different from each other.
Then, one naively expects a phase transition(s) between these two states
by changing the parameter λ = J_H/J_K.

The striking feature of this model is that
there exists a local conserved quantity.
In the presence of the interlayer coupling J_H, W_p no longer commutes with the Hamiltonian, but the product X_p=W_p,1W_p,2 remains a local Z_2 conserved quantity [see Fig. <ref>(b)].
This is because two spin operators with dimer site i in X_p
on plaquette p have the same spin component,
and X_p commutes with S_i,1· S_i,2.
Therefore, the eigenstates of the Hamiltonian given in Eq. (<ref>) is characterized by the eigenvalue of X_p, ±1, in each bi-hexagon p.
The existence of the local conserved quantity leads to the absence of spin correlations except for NN sites in Kitaev layers and interlayer dimers.
Then, we can say that the ground state is always nonmagnetic
even in the presence of interlayer coupling J_H.


To clarify another symmetry of the present bilayer Kitaev system given in Eq. (<ref>), we adopt the bond-operator representation <cit.>,
which is useful for spin-dimer
systems <cit.>.
In this representation, the four local bases on each dimer i are taken as

    s_i   =s_i^†0=1/√(2)(↑↓_i-↓↑_i),
    t_x_i   =t_xi^†0=-1/√(2)(↑↑_i-↓↓_i),
    t_y_i   =t_yi^†0=i/√(2)(↑↑_i+↓↓_i),
    t_z_i   =t_zi^†0=1/√(2)(↑↓_i+↓↑_i),

where s_i^† and t_α i^† (α=x,y,z) are the creation operators of the singlet and triplets on the dimer site i, respectively, and 0 is their vacuum.
We assume that these bond operators
behave as bosons and impose the local constraint s_i^† s_i +∑_α t_α i^† t_α i=1 on each dimer i so as to reproduce the commutation relation of an S=1/2 spin.
By means of the bond operators, the Hamiltonian given in Eq. (<ref>) is rewritten as

    H=   -J_K/2∑_ij_α(s_i s_j^† t_α i^† t_α j + s_i s_j t_α i^† t_α j^† +  H.c.
               
    -∑_ββ'γγ'ϵ_αβγϵ_αβ'γ't_β i^† t_γ it_β' j^† t_γ' j)
       +J_H∑_i(-3/4s_i^† s_i+1/4∑_α t_α i^† t_α i)
       -∑_iμ_i(s_i^† s_i +∑_α t_α i^† t_α i-1).

where ϵ_αβγ stands for the Levi-Civita symbol and μ_i is a Lagrange multiplier on each dimer to impose the local constraint.

We find that the number of each boson is not conserved, but
each boson are created or annihilated as a pair.
For example, the second line in Eq. (<ref>)
is expanded as
t_iy^† t_iz t_jz^† t_jy-t_iy^† t_iz t_jy^† t_jz+ H.c.
when α=x.
Therefore, the parity of the particle number is conserved
in each component of bosons,
In other words, the parity operators,
P_s=exp[iπ∑_is_i^† s_i] and
P_t_α=exp[iπ∑_ it_iα^† t_iα]
(α=x,y,z) commute with the Hamiltonian.






We note that a direct product of the four local states on a dimer,
s_i, t_x_i, t_y_i, and t_z_i is
the eigenstate of X_p
since

    σ_i1^ασ_i2^αs_i=-s_i,  σ_i1^ασ_i2^αt_β_i=
    (1-2δ_αβ)t_β_i.

In addition, the number operators of the bond operators are given as projectors onto the singlet and triplet states:

    s_i^† s_i=1/4- S_i1· S_i2,
       t_α i^†  t_α i =1/4+ S_i1· S_i2-2S_i1^α S_i2^α,

leading to the fact that each parity operator commutes with X_p.
Therefore, the Hilbert space of the Hamiltonian can be classified into
each subspace specified by [P_t_x, P_t_y, P_t_z, {X_p}], where
P_s is determined by P_t_α and the number of dimers N
due to the local constraint.
For example, for the state Φ_s=∏_i s_i^†0,
P_t_x=P_t_y=P_t_z=+1 and X_p=+1 for all plaquettes,
which is schematically shown in Fig. <ref>(a).
Here, we can consider a state Φ_t_xi=t_ix^†∏_j≠ is_j^†0=t̃_ix^†Φ_s,
where we have introduced the triplet excitation (triplon) creation operator
from the singlet state,
t̃_iα^† = t_iα^† s_i.
This state is also an eigenstate of X_p with the eigenvalues of -1
for the adjacent two plaquettes with the shared x bond,
where the triplon is sit on one of the edge sites.
[see Fig. <ref>(b)].
In the case, the eigenvalue of X_p' is +1 for the other plaquettes p'.
From the above discussion,
one can determine the spatial configuration of eigenvalues of X_p,
as shown in Fig. <ref>(b).
The eigenvalues of X_p for the multiple triplet-excited state can be also discussed in the same manner, as shown in Fig. <ref>(c).
Note that this is similar to the configuration of W_p
for the state where a spin is operated to the ground state
in the single-layer Kitaev model.
The details of the excitation spectrum of triplons will be given
in Sec. <ref>.




§ METHODS




To analyze the properties of the ground state in Eq. (<ref>),
we mainly employ the ED method in finite size clusters.
The clusters used in the present calculations are
shown in Fig. <ref>(c).
In the ED calculations, we utilize the presence of the conserved quantities,
X_p and the parities P_s, P_t_x, P_t_x, and P_t_z,
to reduce the matrix dimensions.
In addition to the ED calculations,
we also make use of the bond-operator MF theory <cit.> and
cluster expansion technique <cit.>
to perform the comprehensive analysis.
In the bond-operator MF theory, the local constraint is changed to the global one by introducing the uniform chemical potential μ=μ_i and the Bose condensation of singlets are assumed as s_i=s_i^†=s̅ <cit.>.
Owing to these assumptions, the Hamiltonian is reduced to the free boson system of the triplet excitations as

    H_ MF=   (-3/4J_Hs̅^2-μs̅^2+μ)N+(J_H/4-μ)∑_i,α t_α i^† t_α i
       -J_Ks̅^2/2∑_ij_α( t_α i^† t_α j +  t_α i^† t_α j^† + H.c.),

where we neglect the scattering terms between triplet excitations given in the second line of Eq. (<ref>), and μ and s̅ are determined self-consistently.

We also use the cluster expansion technique <cit.>,
which is one of the powerful methods to discuss the quantum phase transitions
in the frustrated quantum spin systems such as
J_1-J_2 <cit.>,
orthogonal-dimer <cit.>
and Kitaev-Heisenberg models <cit.>.
In the bilayer model,
we begin with interlayer dimer singlets <cit.>.
As discussed before, there exists the parity symmetry
for the number of singlets and triplets in our model, and thereby
the odd-order coefficients never appear
in the ground state energy E_λ/N
and spin correlations ⟨ S_i,1· S_i,2⟩.
We compute the power series up to the 30th order for the above quantities.
Furthermore,
exploiting the first-order inhomogeneous differential method <cit.>,
we deduce the quantities far from the dimer limit λ→∞.

To examine thermodynamic quantities at finite temperatures,
we employ the thermal pure quantum (TPQ) state
approach <cit.>.
In the paper, we treat the cluster with N=12 [see Fig. <ref>(c)].
We prepare more than 20 random vectors for initial states, and the physical quantities are calculated by averaging the values generated by these initial states.










§ NUMERICAL RESULTS





 §.§ Ground State Properties






In this section, we discuss ground state properties
of the bilayer Kitaev model given in Eq. (<ref>).
Figure <ref>(a) shows the λ dependence of
the ground-state energy E_λ.
We obtain the smooth curve by the ED calculations.
We find that the energy obtained by the bond-operator MF approximation
is lower than that of the ED.
This should be due to its artifacts originating from the fact that
triplet-triplet correlations are not taken into account correctly.
With increasing λ, the difference between the ED and bond-operator
MF results becomes negligible,
meaning that the bond-operator MF approximation is justified
in the large λ region.
We also show the energy obtained by the cluster expansion method
as the solid line in Fig. <ref>(a).
Surprisingly, this almost coincides with the ED results
except for λ≲ 0.2 though the method is also an approach
from the large λ limit.
This indicates that the dimer singlet state,
which is adiabatically connected to the direct product state Φ_s,
is realized in the region.

In the small λ region, a large size dependence appears
in the ED results, as shown in Fig. <ref>(b),
in contrast to the large λ region.
With increasing the system size, the energy in the limit of λ=0 should
approach the exact value although it depends on the cluster shape.
An important point is that
the energy curve tends to have the bend structure around λ∼ 0.06.

The peculiar λ dependence of the ground-state energy
can be clearly seen by taking its derivative.
Figure <ref>(c) shows the second derivative of E_λ
as a function of λ.
It is found that the peak structure in each cluster
develops with increasing the system size
while the dependence of the cluster shape is also observed.
This indicates the existence of a first-order quantum phase transition
at λ=λ_c(∼ 0.06) in the thermodynamic limit.





To clarify the nature of two distinct phases around λ_c,
we also calculate the interlayer spin correlations
C_s=1/N∑_i  S_i,1· S_i,2.
The results are shown in Fig. <ref>.
When λ→∞, the singlet dimer is realized in each site and
C_s=-3/4.
In the large λ region, the strong interlayer coupling stabilizes
the dimer singlet state with a short correlation length.
Therefore, the ED results little depend on the cluster size and are
in a good agreement with those obtained by the cluster expansion,
implying that the dimer singlet state is realized.
With decreasing λ, the Ising coupling J_K suppresses
interlayer spin correlations and thereby
the absolute value of C_s decreases monotonically.
Around λ=λ_c, large cluster-size and shape dependence appears
in the ED results and it is hard to extrapolate the quantity
in the thermodynamic limit.
Nevertheless, the results for the larger clusters are still consistent
with the cluster expansion, which implies that the dimer-singlet state
is realized
when λ>λ_c.

On the other hand, in the smaller λ region,
the spin correlation is suddenly changed and take a smaller value,
suggesting the realization of the Kitaev QSL expected
in the two independent Kitaev models.
It is worth noting that λ_c=(J_H/J_K)_c∼ 0.06 is close
to the spin gap of the Kitaev model, Δ_K/J_K∼ 0.065.
This implies that the Kitaev QSL against the interlayer Heisenberg coupling is maintained by the existence of the spin gap inherent in the Kitaev model.





 §.§ Excitation spectrum








Next, we discuss the excitation structure of the bilayer Kitaev model.
It is known that in the Kitaev model (λ=0),
there is a Majorana continuum in the excited states.
On the other hand, in the dimer limit λ→∞,
due to the spin gap with the excitation energy J_H,
discrete excited levels appear, corresponding to the number of triplons.
Here, we first analyze the excitation spectrum by means of
the bond operator method.
The MF Hamiltonian Eq. (<ref>) indicates that
the α component of the local triplet excitation is hybridized
with that only on the NN α bond.
Thus, the triplet excitation t_α is localized on each α bond
and forms the bonding and antibonding states.
Figure <ref> shows the one-body spectrum of the triplet excitation obtained by the bond operator MF approximation.
The introduction of the Kitaev interaction splits excitation energy of the spin gap into two levels, which corresponds to the formation of the bonding and antibonding states.

The fact that triplet excitations are localized on the corresponding NN bonds can be confirmed exactly by considering the spatial configuration of
the local conserved quantity X_p.
This is done in the similar manner as the proof for the absence of
the long-range spin correlations
in the single layer Kitaev model <cit.>.
As shown in Fig. <ref>(b), in the one-triplon state Φ_t_xi on dimer site i, the two eigenvalues of X_p are flipped from the ground state.
On the other hand, one-triplon state Φ_t_xj
on another dimer site j possesses
a different configuration of eigenvalues of X_p, and therefore
Φ_t_xi|Φ_t_xj=Φ_st̃_xit̃_xj^†Φ_s=0 expect for the case where the pair ij is on the NN x bond.
As for one-triplon state with a different component, Φ_t_α j  (α≠ x),
its parity is different from that of the state Φ_t_xj
and these two states are never hybridized.
The above consideration leads to localized bonding and antibonding states for one triplon excitation,
which are 3N/2-fold degenerate each.
The corresponding energies are
Δ_±/J_H=1± 1/(2λ) + 3/(8λ^2) + ⋯.
In addition, one can considers the multiple excitations,
which are shown as the dashed-dotted lines in Fig. <ref>.
Note that a pair of two triplons with the α component
on an α bond can be mixed with the singlet state Φ_s,
and hence, the bound state may have a dispersion
in the presence of the Kitaev interaction with O(1/λ^2).
Furthermore, two triplon states with different components should form
the bound state.
Therefore, two triplon states should be split into
some dispersionless bound states and dispersive bands
on the introduction of the Kitaev coupling J_K.

To confirm the above consideration, we calculate the excitation spectrum using the ED method.
The results are presented in Fig. <ref>.
The low-energy spectra obtained by the bond-operator MF approximation
is well reproduced by the finite-size calculation
in the large λ case because of
the existence of the localized excitations.
We have also confirmed that the first and second excited states
in λ≳ 1.5 are 3N/2-fold degenerate,
which is adiabatically connected to 3N-fold degenerate states
coming from local triplet excitations at λ→∞.
These degeneracies are consistent with those of the bonding and antibonding states of the triplet excitations discussed above.
At λ∼ 1.5, excitation energies originating
from one-triplon and two-triplon excited states intersect without mixing.
This is understood from the fact that the parity of the triplon number
is conserved in the present system,
which prohibits the mixing between one-triplon and two-triplon excited states.
We wish to note that the lowest excitation energy is always finite.
Furthermore, we can not find any tendency for closing the gap.
This is consistent with the fact that the singlet-dimer state is realized in the region λ>λ_c.

In the Kitaev limit, the excitation spectrum shows a continuum
coming from the Majorana fermions, as mentioned before.
On the other hand, low-energy one-triplon states are localized, whose excitation energy may be proportional to the interlayer spin correlation C_s.
We expect that, with decreasing λ,
the discrete levels disappear at
the first-order quantum phase transition point λ_c
and the Majorana continuum might appear below λ_c.
However, larger cluster calculations are needed
to clarify the spectral change around λ_c,
and this remains as a future work.






§ FINITE TEMPERATURE PROPERTIES








Finally, we discuss finite temperature properties of the bilayer Kitaev model.
Figure <ref> shows the temperature dependence of the specific heat
for several λ by using the TPQ state for the N=12 cluster.
When λ=0, the system is reduced to the decoupled Kitaev models and
the double peak structure appears
at T=T_L∼ 0.044J_K and T_H∼ 0.38J_K in the N=12 cluster.
This is consistent with
the fact that the spin degrees of freedom is split into itinerant Majorana
fermions and localized gauge fluxes in each single-layer Kitaev model.
By comparing with the curves with λ=0.2,
the lower characteristic temperature T_L
slightly shifts to higher temperatures, while the other little shifts.
Further increasing λ,
these two peaks merge into a single peak,
which is a Schottky-type peak characteristic of the gapped systems.






To clarify how the double peak structure in the single-layer Kitaev model is gradually changed,
we show in Fig. <ref> the density plot of the specific heat
on the plane of λ and temperature.
The double peaks exhibiting at λ=0 clearly appear
below λ∼ 0.2.
With increasing λ, the temperature of the higher-temperature peak hardly changes whereas the lower one depends on λ.
We find that these peaks merge into the single peak around λ∼ 0.6.
This suggests that the proximity effect of the Kitaev QSL emerges
at higher temperatures even above λ_c.
Similar effect has been discussed in the Kitaev-Heisenberg model,
where the double-peak structure is seen even in the magnetic ordered phase close to the Kitaev QSL ground state <cit.>.






§ DISCUSSION










In this section, we briefly discuss magnetic properties of
the multilayer Kitaev systems,
where the Kitaev layers are connected by the Heisenberg couplings,
as shown in Fig. <ref>.
In the system, we can define the local Z_2 conserved quantity
X̃_p=∏_n W_p,n (see Fig. <ref>),
which is similar to the bilayer system discussed in the previous sections.
When the number of layers is large enough,
this system is regarded as one-dimensional Heisenberg chains
coupled with Kitaev interactions.
In the case, its spin correlation length shows power-law decay
along each chain direction.
However, the existence of the local conserved quantity X̃_p
indicates that there are no spin correlations
except for NN sites in Kitaev layers and sites
belonging to the same Heisenberg chain.
This suggests that the ground state of the stacked Kitaev model
with the interlayer Heisenberg couplings shown Fig. <ref>
is nonmagnetic as three-dimensional longer range spin correlations are absent.






§ SUMMARY





In summary, we have investigated the ground-state and the finite temperature
properties of the bilayer Kitaev model.
The results obtained by the ED, bond-operator MF,
and cluster expansion methods suggest
the existence of a first-order quantum phase transition
between the Kitaev QSL and singlet-dimer states in the thermodynamic limit.
In the excitation spectrum, the one-triplon excitation is localized
in the singlet-dimer phase,
which is proved by the existence of a local conserved quantity.
We have also discussed finite-temperature properties, where
the double-peak structure intrinsic to the Kitaev QSL appears
even in the singlet-dimer phase near the phase boundary.
Furthermore, we have shown that three-dimensional long-range
spin correlations are also absent in the stacked multilayer Kitaev systems
with arbitrary layer numbers.




The authors would like to thank S. Suga for valuable discussions.
This work is supported by Grant-in-Aid for Scientific Research from
JSPS, KAKENHI Grant Number JP17K05536, JP16H01066 (A.K.) and JP16K17747, JP16H02206, JP16H00987 (J.N.)
Parts of the numerical calculations were performed
in the supercomputing systems in ISSP, the University of Tokyo.






