Majorana dynamical mean-field study of spin dynamics at finite temperatures in the honeycomb Kitaev model

# Majorana dynamical mean-field study of spin dynamics at finite temperatures in the honeycomb Kitaev model

Junki Yoshitake Department of Applied Physics, University of Tokyo, Bunkyo, Tokyo 113-8656, Japan    Joji Nasu Department of Physics, Tokyo Institute of Technology, Meguro, Tokyo 152-8551, Japan    Yasuyuki Kato Department of Applied Physics, University of Tokyo, Bunkyo, Tokyo 113-8656, Japan    Yukitoshi Motome Department of Applied Physics, University of Tokyo, Bunkyo, Tokyo 113-8656, Japan
August 5, 2019
###### Abstract

A prominent feature of quantum spin liquids is fractionalization of the spin degree of freedom. Fractionalized excitations have their own dynamics in different energy scales, and hence, affect finite-temperature () properties in a peculiar manner even in the paramagnetic state harboring the quantum spin liquid state. We here present a comprehensive theoretical study of the spin dynamics in a wide range for the Kitaev model on a honeycomb lattice, whose ground state is such a quantum spin liquid. In this model, the fractionalization occurs to break up quantum spins into itinerant matter fermions and localized gauge fluxes, which results in two crossovers at very different scales. Extending the previous study for the isotropic coupling case [J. Yoshitake, J. Nasu, and Y. Motome, Phys. Rev. Lett. 117, 157203 (2016)], we calculate the dynamical spin structure factor , the NMR relaxation rate , and the magnetic susceptibility while changing the anisotropy in the exchange coupling constants, by using the dynamical mean-field theory based on a Majorana fermion representation. We describe the details of the methodology including the continuous-time quantum Monte Carlo method for computing dynamical spin correlations and the maximum entropy method for analytic continuation. We confirm that the combined method provides accurate results in a wide range including the region where the spins are fractionalized. We find that also in the anisotropic cases the system exhibits peculiar behaviors below the high- crossover whose temperature is comparable to the average of the exchange constants: shows an inelastic response at the energy scale of the averaged exchange constant, continues to grow even though the equal-time spin correlations are saturated and almost independent, and deviates from the Curie-Weiss behavior. In particular, when the exchange interaction in one direction is stronger than the other two, the dynamical quantities exhibit qualitatively different dependences from the isotropic case at low , reflecting the opposite parity between the flux-free ground state and the flux-excited state, and a larger energy cost for flipping a spin in the strong interaction direction. On the other hand, when the exchange anisotropy is in the opposite way, the results are qualitatively similar to those in the isotropic case. All these behaviors manifest the spin fractionalization in the paramagnetic region. Among them, the dichotomy between the static and dynamical spin correlations is unusual behavior hardly seen in conventional magnets. We discuss the relation between the dichotomy and the spatial configuration of gauge fluxes. Our results could stimulate further experimental and theoretical analyses of candidate materials for the Kitaev quantum spin liquids.

preprint: APS/123-QED

## I Introduction

Quantum many-body systems show various intriguing phenomena which cannot be understood as an assembly of independent particles. One of such phenomena is fractionalization, in which the fundamental degree of freedom in the system is fractionalized into several quasiparticles. A well-known example of such fractionalization is the fractional quantum Hall effect, in which the Hall conductance shows plateaus at fractional values of ( is the elementary charge and is the Planck constant) Tsui1982 (); Stormer1999 (). In this case, the quasiparticles carry a fractional value of the elementary charge, as a collective excitation of the elementary particle, electron. This is fractionalization of charge degree of freedom. On the other hand, another degree of freedom of electrons, spin, can also be fractionalized. Such a peculiar phenomenon has been argued for quantum many-body states in insulating magnets, e.g., a quantum spin liquid (QSL) state.

QSLs are the magnetic states which preserve all the symmetries in the high-temperature() paramagnet even in the ground state and evade a description by conventional local order parameters. A typical example of QSLs is the resonating valence bond (RVB) state, proposed by P. W. Anderson Anderson1973 (). The RVB state is a superposition of valence bond states (direct products of spin singlet dimers), which does not break either time reversal or translational symmetry. In the RVB state, the spin degree of freedom is fractionalized: the system exhibits two different types of elementary excitations called spinon and vison Kivelson1987 (); Senthil2000 (). Spinon is a particlelike excitation carrying no charge but spin . Meanwhile, vison is a topological excitation characterized by the parity of crossing singlet pairs with its trace. Another example of QSLs is found in quantum spin ice systems, in which peculiar excitations are assumed to be magnetic monopoles, electric gauge charges, and artificial photons resulting from fractionalization of the spin degree of freedom Hermele2004 (); Castelnovo2008 ().

Among theoretical models for QSLs, the Kitaev model has attracted growing interest, as it realizes the fractionalization of quantum spins in a canonical form Kitaev2006 (). The Kitaev model is a localized spin model defined on a two-dimensional honeycomb lattice with bond-dependent anisotropic interactions (see Sec. II.1). In this model, the ground state is exactly obtained as a QSL, in which quantum spins are fractionalized into itinerant Majorana fermions and localized gauge fluxes. The fractionalization affects the thermal and dynamical properties in this model. For instance, the different energy scales between the fractionalized excitations appear as two crossovers at largely different scales; in each crossover, itinerant Majorana fermions and localized gauge fluxes release their entropy, a half of per site Nasu2014 (); Nasu2015 (). Also in the ground state, the dynamical spin structure factor shows a gap due to the flux excitation and strong incoherent spectra from the composite excitations between itinerant Majorana fermions and localized gauge fluxes Knolle2014 (). Such incoherent spectra were indeed observed in recent inelestic neutron scattering experiments for a candidate for the Kitaev QSL, -RuCl Banerjee2016 (); Do_preprint (). The magnetic Raman scattering spectra also shows a broad continuum dominated by the itinerant Majorana fermions, in marked contrast to conventional insulating magnets Knolle2014b (). Such a broad continuum was experimentally observed also in -RuCl Sandilands2015 (). Furthermore, the dependence of the incoherent response was theoretically analyzed and identified as the fermionic excitations emergent from the spin fractionalization Nasu2016 (); Glamazda2016 ().

In the previous study, the authors have studied dynamical properties of the Kitaev model at finite by a newly developed numerical technique, the Majorana dynamical mean-field method Yoshitake2016 (). Indications of the spin fractionalization were identified in the dependences of , the relaxation rate in the nuclear magnetic resonance (NMR), , and the magnetic susceptibility . In the previous study, however, the results were limited to the case with the isotropic exchange constants, despite the anisotropy existing in the Kitaev candidate materials Johnson2015 (); Yamaji2014 (). In the present paper, to complete the analysis, we present the numerical results of the dynamical quantities for anisotropic cases. We also provide the comprehensive description of the theoretical method, including the details of the cluster dynamical mean-field theory (CDMFT), the continuous-time quantum Monte Carlo (CTQMC) as a solver of the impurity problem to calculate the dynamical spin correlations, and the maximum entropy method (MEM) for the analytic continuation. We discuss a prominent feature proximate to the QSL, i.e., the dichotomy between static and dynamical spin correlations, from the viewpoint of the fractionalization of spins.

The paper is organized as follows. In Sec. II, after introducing the Kitaev model and its Majorana fermion representation, we describe the details of the CDMFT, CTQMC, and MEM. In Sec. III, we show the numerical results for , , and while changing the anisotropy in the exchange constants. In Sec. IV, we discuss the dichotomy between the static and dynamical spin correlations by comparing the dependences for the uniform and random flux configurations. The cluster-size dependence in the CDMFT is examined in Appendix A. The accuracy of MEM is also examined in Appendix B in the one-dimensional limit where the dynamical properties can be calculated without analytic continuation. We also show the and dependence of spin correlations and the dependence of the Korringa ratio in Appendix C and D, respectively.

## Ii Model and method

In this section, we describe the details of the methods used in the present study, the Majorana CDMFT and CTQMC methods. After introducing the Majorana fermion representation of the Kitaev model in Sec. II.1, we describe the framework of the Majorana CDMFT in Sec. II.2, in which the impurity problem is solved exactly. In Sec. II.3, we introduce the CTQMC method which is applied to the converged solutions obtained by the Majorana CDMFT for calculating dynamical spin correlations. We also touch on the MEM used for obtaining the dynamical spin correlations as functions of real frequency from those of imaginary time in Sec. II.4.

### ii.1 Kitaev model and the Majorana fermion representation

We consider the Kitaev model on a honeycomb lattice, whose Hamiltonian is given by Kitaev2006 ()

 H=−∑pJp∑⟨j,j′⟩pSpjSpj′, (1)

where , , and , and the sum of is taken for the nearest-neighbor (NN) sites on three inequivalent bonds of the honeycomb lattice, as indicated in Fig. 1(a); is the component of the spin at site . Hereafter, we denote the average of as and set the energy scale as , i.e., , and parametrize the anisotropy of the exchange coupling constants as and , where and correspond to the ferromagnetic (FM) and antiferromagnetic (AFM) cases, respectively. We note that the FM and AFM cases are connected through unitary transformations Kitaev2006 ().

As shown by Kitaev Kitaev2006 (), the model is soluble and the exact ground state is obtained as a QSL. The spin correlations are extremely short-ranged: are nonzero only for the NN sites on the bonds as well as the same site  Baskaran2007 (). Hereafter, we denote the NN correlations as . There are two types of QSL phases depending on the anisotropy in : one is a gapless QSL realized in the region with including the isotropic point (), while the other is gapful for . The ground state has nontrivial fourfold degeneracy in the thermodynamic limit Mandal2012 ().

The exact solution for the ground state was originally obtained by introducing four types of Majorana fermions for each spin Kitaev2006 (). In this method, the Hilbert space in the original spin representation, , is extended to in the Majorana fermion representation ( is the number of spins). Thus, to calculate physical quantities, such as spin correlations, it is necessary to make a projection from the extended Hilbert space to the original one.

Soon later, however, another way of solving the model was introduced by using only two types of Majorana fermions Chen2007 (); Feng2007 (); Chen2008 (), in which the projection is avoided as the Hilbert space is not extended. In this method, the spin operators are written by spinless fermions by applying the Jordan-Wigner transformation to the one-dimensional chains composed of two types of bonds, say, the and bonds. Then, by introducing two Majorana fermions and for the spinless fermions, the Hamiltonian in Eq. (1) is rewritten as

 H=iJx4∑(j,j′)xcjcj′−iJy4∑(j,j′)ycjcj′−iJz4∑(j,j′)zηrcjcj′, (2)

where the sum over is taken for the NN sites on a bond with . is defined on each bond connecting and sites ( is the index of the bond). Here, is considered as a variable taking , as commutes with the total Hamiltonian as well as with other and as . Thus, the model in Eq. (2) describes itinerant Majorana fermions (called matter fermions) coupled to the variables (called gauge fluxes). The ground state is given by all , giving QSLs with gapless or gapful excitations depending on , as in the original Kitaev’s solution.

In the present numerical study at finite , we adopt the Majorana representation used in Eq. (2). This is because the form of Eq. (2) is suitable for the CDMFT calculations (see Sec. II.2), as the interaction term, the third term in Eq. (2), only lies on bonds. In this study, we apply the CDMFT to deal with thermal fluctuations and compute the static quantities. For calculating dynamical quantities, we apply the CTQMC method to the converged solutions obtained the CDMFT. While the framework was briefly introduced in Ref. Yoshitake2016 (), we describe further details in the following sections.

### ii.2 Cluster dynamical mean-field theory in the Majorana fermion representation

As presented in the previous study by the real-space QMC simulation Nasu2015 (), spacial correlations between develop at low . To take into account such spacial correlations, we adopt a cluster extension of DMFT Kotliar2001 (). As the Majorana Hamiltonian in Eq. (2) is formally similar to the Falicov-Kimball model or the double-exchange model with Ising localized moments, we follow the DMFT framework for the double-exchange model Furukawa1994 ().

In the CDMFT, we regard the whole lattice as a periodic array of clusters. The Hamiltonian in Eq. (2) is rewritten into the matrix form of

 H=∑γ,γ′,j,j′12H0γ,j;γ′,j′cγ,jcγ′,j′+∑γ,j,j′12H{η}j,j′cγ,jcγ,j′, (3)

where and are the indices for the clusters, and and denote the sites in each -site cluster. The coefficient in Eq. (3) is introduced to follow the notation in Ref. Nilsson2013 (). In Eq. (3), the first term corresponds to the first and second terms in Eq. (2), while the second term is for the third term. Green’s function for Eq. (3) is formally written as

 G(k,iωn)=(iωn−2H0(k)−Σ(k,iωn))−1, (4)

where is the Matsubara frequency ( is an integer, and the Boltzmann constant and the reduced Planck constant are set to unity), is the self-energy, and is the Fourier transform of in Eq. (3) given by the matrix:

 H0j,j′(k)=∑γH0γ,j;0,j′e−ik⋅rγ, (5)

where is the coordinate of the cluster .

Following the spirit of the DMFT Metzner1989 (); Georges1996 (), we omit the dependence of the self-energy: . In this approximation, local Green’s function is defined within a cluster as

 Gj,j′(iωn)=1N′∑k[(iωn−2H0(k)−Σ(iωn))−1]j,j′, (6)

where is the number of clusters in the whole lattice (), and and denotes the sites in the cluster. The Weiss function is introduced to take into account the correlation effects in other clusters as

 G0j,j′(iωn)−1=Gj,j′(iωn)−1+Σj,j′(iωn). (7)

In order to take into account the interaction in Eq. (3) within the cluster that we focus on, we consider the impurity problem for the cluster described by the effective action in the path-integral representation for Majorana fermions Nilsson2013 (). The partition function is given by

 Z=∑{η}Z{η}, (8)

where

 Z{η}=∫Dχexp(−S{η}% eff). (9)

Here, the sum of in Eq. (8) runs over all possible configurations of , and in Eq. (9); is the Grassmann number corresponding to the Majorana operator (more precisely, following the notation in Ref. Nilsson2013 ()). The effective action is given by

 S{η}eff= −T∑j,j′,n≥0χj,−ωn(G0(iωn))−1j,j′χj′,ωn +2T∑j,j′,n≥0χj,−ωnH{η}j,j′χj′,ωn. (10)

For a given configuration of , the impurity problem defined by Eq. (9) is exactly solvable because it is nothing but a free fermion problem. Green’s function is obtained as

 [(G{η}(iωn))−1]j,j′=[(G0(iωn))−1]j,j′−2H{η}j,j′. (11)

Note that we slightly modified the notation from the previous study in Ref. Yoshitake2016 (). Then, local Green’s function for the impurity problem is calculated by

 Gimpj,j′(iωn)=∑{η}P({η})G{η}j,j′(iωn), (12)

where is the statistical weight for the configuration given by

 P({η})=Z{η}/∑{η}Z{η}. (13)

is obtained from Green’s functions as

 Z{η}=∏n≥0det[−G{η}(iωn)]. (14)

We note that is obtained exactly by computing and for all configurations of in the -site cluster Udagawa2012 (). The self-energy for the impurity problem is obtained as

 (15)

In the CDMFT, the above equations, Eqs. (6), (7), (12), and (15), are solved in a self-consistent way. The self-consistent condition is given by

 G(iωn)=Gimp(iωn), (16)

namely, the calculation is repeated until local Green’s function in Eq. (6) agrees with Green’s function calculated for the impurity problem in Eq. (12).

The Majorana CDMFT framework provides a concise calculation method for dependences of static quantities of the Kitaev model, such as the specific heat and the equal-time spin correlations . It is worth noting that the CDMFT calculations can be performed without any biased approximation except for the cluster approximation: the exact enumeration for all the configurations in Eq. (12) enables the exact calculations for the given cluster. Furthermore, the cluster-size dependence is sufficiently small at all the range above the critical temperature for the artificial phase transition due to the mean-field nature of the CDMFT, as demonstrated for the isotropic case with in the previous study Yoshitake2016 () (see also Sec. III.1 and Appendix A for the anisotropic cases). On the other hand, for obtaining dynamical quantities, such as the dynamical spin correlations ( is the imaginary time), we need to make an additional effort beyond the exact enumeration in the CDMFT, as discussed in the next subsection.

In the CDMFT+CTQMC calculations in Sec. III, we use the 26-site cluster shown in Fig. 1(a). In Appendix A, we examine the dependence on the cluster size as well as shape.

### ii.3 Continuous-time quantum Monte Carlo method

In order to calculate the dynamical spin correlations , we need to take into account the imaginary-time evolution of that compose the conserved quantities , e.g., ; the sign depends on the sublattice on the honeycomb structure. For this purpose, we adopt the CTQMC method based on the strong coupling expansion Werner2006 (). In this method, on an bond is calculated as

 ⟨Szj(τ)Szj′⟩=∑{η}′,ηr0=±1P({η}′,ηr0)⟨Szj(τ)Szj′⟩{η}′, (17)

where represents the configurations of except for on the bond. is obtained from the converged solution of the Majorana CDMFT in Sec. II.2. is the dynamical spin correlation on the bond calculated by the CTQMC method for each configuration . The sum of runs over all possible configurations of within the cluster. Note that Eq. (17) is derived from the fact that commutes with in , whereas it does not commute with . Thus, for a given , the interaction lies only on the bond, and hence, it is sufficient to solve the two-site impurity problem in the CTQMC calculations. The two-site impurity problem is defined by the integration in Eq. (9) on whose does not belong to the bond. Then, we obtain

 S{η}′eff=S{η}′hyb+Slocal, (18)

where

 S{η}′hyb= −∑j,j′∫β0dτ∫β0dτ′χj(τ)Δ{η}′j,j′(τ−τ′)χj′(τ′), (19) Slocal= ∑j,j′∫β0dτχj(τ)(δj,j′2∂∂τ+H{η}j,j′)χj′(τ), (20)

and in Eqs. (19) and (20) are the sites on the bond; is the inverse temperature. In Eq. (19), the hybridization function is calculated from in the converged solution of CDMFT as follows. Let us define the matrix as a submatrix of , as

 ~G{η}j,j′(iωn)=G{η}j,j′(iωn). (21)

Then, the hybridization function is given as a function of the Matsubara frequency in the form

 Δ{η}′j,j′(iωn)=[~G{η}(iωn)]−1j,j′−(iωn−2H{η}j,j′). (22)

Note that does not depend on , which is straightforwardly shown by the matrix operations in the right hand side. Converting Eq. (22) to the imaginary-time representation, we obtain

 Δ{η}′j,j′(τ)=T2∑ne−iωnτΔ{η}′j,j′(iωn). (23)

Given Eqs. (18)-(20), the partition function of the system is expanded in terms of as

 ZZlocal =∫Dχe−S{η}′hybe−Slocal∫Dχe−Slocal=⟨e−S{η}′hyb⟩local =∑d,i0,...,i2d−1∫β0dτ0...∫β0dτ2d−11d!⟨χi0(τ0)...χi2d−1(τ2d−1)⟩localPf(^Δ{η}′(d,i0,τ0,...,i2d−1,τ2d−1)), (24)

where is the partition function for the two sites described by , and represents the expectation value in the two-site problem as

 ⟨A⟩local=∫DχAe−Slocal∫Dχe−Slocal. (25)

In the second line of Eq. (24), is the order of in the expansion of , Pf() is the Pfaffian of skew-symmetric matrix , and is a matrix, whose element is given by

 ^Δ{η}′(d,i0,τ0...,i2d−1,τ2d−1)m,n=Δ{η}′im,in(τm−τn). (26)

We note that this is the first formulation of the CTQMC method with using the Pfaffian in the weight function to our knowledge, whereas a QMC simulation in the Majorana representation has been introduced for itinerant fermion models Li2015 ().

In the CTQMC calculation, we perform MC sampling over the configurations by using the integrand in Eq. (24) as the statistical weight for each configuration. In each MC step, we perform an update from one configuration to another; for instance, an increase of the order of expansion as to by adding . To judge the acceptance of such an update, we need to calculate the ratio of the Pfaffian. This is efficiently done by using the fast update algorithm, as in the hybridization expansion scheme for usual fermion problems (for example, see Ref. Rubtsov2005 ()). For the above example of increasing , the ratio is calculated by adding two rows and columns in the matrix as

 Pf(^Δ{η}′(d,i0,τ0,...,i2d−1,τ2d−1))Pf(^Δ{η}′(d+1,i0,τ0,...,i2d+1,τ2d+1)), (27)

whose calculation cost is in the order of by using the fast update algorithm. On the other hand, in Eq. (24), is obtained as the average in the two-site problem, which can be calculated by considering the imaginary-time evolution of all the four states in the two-site problem.

Then, the dynamical spin correlation for the configuration , in Eq. (17), is calculated as

 ⟨Szj(τ)Szj′⟩{η}′=ZlocalZ∑d,i0,...,i2d−1∫β0dτ0...∫β0dτ2d−11d! ⟨χi0(τ0)...χi2d−1(τ2d−1)Szj(τ)Szj′⟩local ×Pf(^Δ{η}′(d,i0,...,i2d−1,τ0,...,τ2d−1)). (28)

For the MC sampling, we need to evaluate

 ⟨Tτχi0(τ0)...χi2d−1(τ2d−1)Szj(τ)Szj′⟩local⟨Tτχi0(τ0)...χi2d−1(τ2d−1)⟩local. (29)

This is again calculated by considering the imaginary-time evolution of all the four states in the two-site problem. In the isotropic case with , for are equivalent to . Meanwhile, for the anisotropic case, we compute for by the same technique described above with using the spin rotations or .

In the CTQMC calculations in Sec. III, for each configuration , we typically run MC steps and perform the measurements at every 20 steps, after MC steps for the initial relaxation.

### ii.4 Maximum entropy method

By using the CTQMC method as the impurity solver in the CDMFT, which we call the CDMFT+CTQMC method, we can numerically estimate the dynamical spin correlation as a function of the imaginary time, . To obtain the physical observables, such as the dynamical spin structure factor and the NMR relaxation rate, which are given by the dynamical spin correlations as functions of frequency , we need to inversely solve the equation given by the generic form . In our problem, and correspond to the dynamical spin correlations as functions of imaginary time and real frequency : and . In the following calculations, we utilize the Legendre polynomial expansion following Refs. Boehnke2011 (); Levy2017 ():

 gm=√2m+1∫β0dτPm(x(τ))g(τ), (30)

where is the th Legendre polynomials and . Then, the inverse problem is given by

 gm=∫dωρ(ω)Km(ω), (31)

where

 Km(ω)=√2m+1∫β0dτPm(x(τ))e−ωτ. (32)

For solving the inverse problem, we adopt the maximum entropy method (MEM) Jarrell1996 (). The following procedure is the standard one, but we briefly introduce it to make the paper self-contained. In the MEM, we discretize to , and determine to minimize the function

 F =12∑m,n(gm−~gm)ζ−1C−1m,n(gn−~gn) −δ∑l⎡⎣ρl−ρ(0)l−ρlln⎛⎝ρlρ(0)l⎞⎠⎤⎦, (33)

where and are the coefficients described below, and is a variance-covariance matrix of ; . We take the Legendre expansion up to th order and in the following calculations. In Eq. (33), is the advance estimate of , which we set to be a constant in this study.

Once neglecting the second term in the right hand side of Eq. (33), the minimization of is equivalent to the least squares method. The least squares method is unstable, as is rather insensitive to a change of . The second term, called the entropy term, stabilizes the minimization process. In the following calculations, we set to sufficiently take into account the effect of the entropy term, where the value of is determined self-consistently in each MEM calculation based on the maximum likelihood estimation, called the classical MEM Jarrell1996 () (typically, -). We note that the deviations of from are typically comparable to the statistical errors in the CTQMC calculations. In the following results, we estimate the errors of by the standard deviation between the data for , and in the range where the MEM retains the precision.

In the MEM, should be positive for all . In our problem, the onsite correlation satisfies this condition automatically, whereas for the NN sites on the bond, which is denoted by hereafter, can be negative. (Note that all the further-neighbor correlations beyond the NN sites vanish in the Kitaev model Baskaran2007 ().) To obtain properly, we calculate , which is positive definite for all , and subtract the onsite contributions note (). The accuracy of obtained by the MEM are examined in Appendix B in the one-dimensional limit with , where can be calculated without using the MEM.

## Iii Result

In this section, we present the results obtained by the CDMFT and the CDMFT+CTQMC methods. In Sec. III.1, we present the specific heat and equal-time spin correlations for the NN sites obtained by the CDMFT for the cases with anisotropic , , and . By comparing the results with those by the QMC method Nasu2015 (), we confirm that the CDMFT is valid in the range above the artificial critical temperature close to the low- crossover. In Sec. III.2, III.3, and III.4, we present the CDMFT+CTQMC results for dynamical quantities, i.e., the dynamical spin structure factor, the NMR relaxation rate, and the magnetic susceptibility, respectively, in the qualified range. We discuss the results in comparison with the isotropic case reported previously in Ref. Yoshitake2016 ().

### iii.1 Static quantities: comparison to the previous QMC results

Figure 2 shows the benchmark of the Majorana CDMFT. We compare the specific heat and equal-time spin correlations for NN pairs on the bonds, , obtained by the Majorana CDMFT, with those by QMC in Ref. Nasu2015 (). The data are calculated for the FM case with bond asymmetry: ( and ) and ( and ). While the data of are common to the FM and AFM cases, the sign of is reversed for the AFM case. Note that similar comparison was made for the isotropic case () in Ref. Yoshitake2016 ().

As indicated by two broad peaks in the specific heat in the QMC results, the system exhibits two crossovers owing to thermal fractionalization of quantum spins Nasu2015 (); the crossover temperatures were estimated as and in the isotropic case. In the anisotropic cases, the low- crossover takes place at a lower , i.e., for and for , while the high- one is almost unchanged, i.e., . These behaviors are excellently reproduced by the Majorana CDMFT, except for the low- peak; the CDMFT results show a sharp anomaly at for and for . This is due to a phase transition by ordering of , which is an artifact of the mean-field nature of CDMFT.

On the other hand, the QMC results for the NN spin correlations are also precisely reproduced by the Majorana CDMFT in the wide range above the artificial phase transition temperature . Although they appear to be reproduced even below , there is a small anomaly at associated with the artificial transition, while the QMC data smoothly change around . (Note that the appropriate sum of the NN spin correlations is nothing but the internal energy, and hence, the derivative corresponds to the specific heat.)

Thus, the comparison indicates that the Majorana CDMFT gives quantitatively precise results in the wide range above the artificial transition temperature : in the present cases with and , the CDMFT is reliable for and , respectively. As discussed in the previous study Nasu2015 (), the thermal fractionalization of quantum spins sets in below , which is well above . Thus, the ranges qualified for the CDMFT include the peculiar paramagnetic state showing the thermal fractionalization. In the following sections, we apply the CDMFT+CTQMC method in these qualified ranges to the study of spin dynamics, which was not obtained by the previous QMC method Nasu2015 ().

### iii.2 Dynamical spin structure factor

Figure 3 shows the CDMFT+CTQMC results for the dynamical spin structure factor at several for the FM case with , , and . is calculated as

 S(q,ω) =∑pSp(q,ω), (34) Sp(q,ω) =13N∑j,j′eiq⋅(rj−rj′)Spj,j′(ω), (35)

where is obtained by the MEM described in Sec. II.4 from the imaginary-time correlations by CDMFT+CTQMC. As mentioned above, nonzero contributions in Eq. (35) come from only the onsite and NN-site components of ; we present their dependences in Appendix C. The Brillouin zone and symmetric lines on which is plotted are presented in Fig. 1(b). Although the results at were shown in the previous study Yoshitake2016 (), we present them (for a sightly different set) for comparison. We show the data at four temperatures: , , , and . Note that , , and for , , and , respectively, while for all the cases.

As shown in Fig. 3, at sufficiently high than , does not show any significant dependence for all studied here; shows only a diffusive response centered at , as shown in Figs. 3(j)-3(l). When lowering below , the diffusive weight is shifted to the positive region ranging up to above for all the cases, as shown in Figs. 3(g)-3(i). Simultaneously, a quasi-elastic component grows gradually at . Both the inelastic and the quasi-elastic components show a discernible dependence; in particular, the latter increases the intensity around the point reflecting the FM interactions. While and for from the symmetry, the quasi-elastic response is small (large) around the M-K line compared to that around the M-K line for () because of the anisotropy.

When further lowering and approaching , the quasi-elastic component increases its intensity, while the inelastic response at does not change substantially. In particular, in the case of , the quasi-elastic component is sharpened and develops to a -function like peak as shown in Figs. 3(e) and 3(b). In addition, the broad incoherent weight splits from the coherent peak. These behaviors appear to asymptotically converge onto the result at , where the -function peak appears due to the change of the parity between the ground state and the flux-excited state Knolle2014 () (for the -function peak, see also Fig. 16 in Appendix B). On the other hand, at does not show such a drastic change, and the quasi-elastic component grows continuously, as shown in Figs. 3(f) and 3(c). We note that the results for are qualitatively similar to those for in Figs. 3(d) and 3(a), except for different dependence mentioned above.

Figure 4 shows the results for the AFM case. The overall dependence of is similar to that for the FM case at all : the diffusive response centered at for [Figs. 4(j)-4(l)], the shift of the diffusive weight to the region of and the growth of a quasi-elastic component at below [Figs. 4(g)-4(i)], and the -function like peak for while approaching to [Figs. 4(e) and 4(b)]. The similarity of the dependences of between FM and AFM cases is partly understood by the relation , which holds for [ and are for the FM and AFM cases, respectively]. On the other hand, the dependence is in contrast to the FM case: while the weight of the quasi-elastic response almost vanishes around the point, those on the zone boundary are enhanced in an almost opposite manner to the FM cases. In addition, the incoherent weight at also shows the opposite dependence to the FM case: the weight is stronger around the point than that on the zone boundary. The opposite dependences between the FM and AFM cases directly follow from the relation .

In order to show the dependences of more explicitly, we present in Figs. 5-8 the - plot of at , , and K with the intensity profiles for the same set of used in Figs. 3 and 4. Figure 5 shows the result for the FM case at . The overall weight of shifts from to a large- region when the system is cooled down below . Below , quasi-elastic response gradually grows and develops to the -function like peak. The peak intensity in and is larger than that for , reflecting the anisotropy of the interaction.

Figure 6 shows the corresponding plot for the AFM case at . In contrast to the FM case, the strong quasi-elastic response is seen for , which develops to the -function like peak at low . We note that the dip and shoulder like structures around in the intermediate for the result at may be an artifact originating from low precision in the MEM for this AFM case because of the following reason. As described in Sec. II.4, we calculate for the NN bonds by subtracting the onsite component from , both of which are obtained by the MEM. In the present case, as both of and become large around due to the development of the -function like peak, the relative error becomes large for , which may lead to artificial structures.

Figures 7 and 8 show the results at . As observed in Figs. 3 and 4, for both the FM and AFM cases behave similarly to those at  Yoshitake2016 (). In the anisotropic cases, however, the difference between and is obvious: the quasi-elastic peak for