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

TL;DR

This study uses dynamical mean-field theory and quantum Monte Carlo methods to analyze spin dynamics and fractionalization in the honeycomb Kitaev model across various temperatures, revealing unusual static-dynamic correlation dichotomies.

Contribution

It extends previous work by analyzing anisotropic couplings and provides detailed dynamical properties of the Kitaev model at finite temperatures using advanced computational techniques.

Findings

Identification of two temperature crossovers related to fractionalization

Unusual static-dynamic spin correlation dichotomy

Behavior reflecting gauge flux configurations

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 ($T$) 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 $T$ 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 $T$ scales. Extending the previous study for the isotropic coupling case [J. Yoshitake, J. Nasu, and Y. Motome, Phys. Rev. Lett. ${\textbf 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 and the continuous-time quantum Monte Carlo method based on a Majorana fermion representation. We find that the system exhibits peculiar behaviors below the high-$T$ crossover whose temperature is comparable to the average of the exchange constants, reflecting 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.

Equations

$$\mathcal{H}=-\sum _p{J}_p\sum _{⟨j,j^{\prime }{⟩}_p}{S}_j^p{S}_{j^{\prime }}^p,$$

$$\mathcal{H}=i\frac{J_x}{4}\sum _{(j,j^{\prime }{)}_x}{c}_j{c}_{j^{\prime }}-i\frac{J_y}{4}\sum _{(j,j^{\prime }{)}_y}{c}_j{c}_{j^{\prime }}-i\frac{J_z}{4}\sum _{(j,j^{\prime }{)}_z}{\eta }_r{c}_j{c}_{j^{\prime }},$$

$$\mathcal{H}=\sum _{\gamma ,{\gamma }^{\prime },j,j^{\prime }}\frac{1}{2}{\mathcal{H}}_{\gamma ,j;{\gamma }^{\prime },j^{\prime }}^0{c}_{\gamma ,j}{c}_{{\gamma }^{\prime },j^{\prime }}+\sum _{\gamma ,j,j^{\prime }}\frac{1}{2}{\mathcal{H}}_{j,j^{\prime }}^{\{\eta \}}{c}_{\gamma ,j}{c}_{\gamma ,j^{\prime }},$$

$$G(\mathbf{k},i{\omega }_n)=(i{\omega }_n-2{\mathcal{H}}^0(\mathbf{k})-\Sigma (\mathbf{k},i{\omega }_n){)}^{-1},$$

$${\mathcal{H}}_{j,j^{\prime }}^0(\mathbf{k})=\sum _{\gamma }{\mathcal{H}}_{\gamma ,j;0,j^{\prime }}^0{e}^{-i\mathbf{k}\cdot {\mathbf{r}}_{\gamma }},$$

$$G_{j,j^{\prime }}(i{\omega }_n)=\frac{1}{N^{\prime }}\sum _{\mathbf{k}}{\left[(i{\omega }_n-2{\mathcal{H}}^0(\mathbf{k})-\Sigma (i{\omega }_n){)}^{-1}\right]}_{j,j^{\prime }},$$

$${\mathcal{G}}_{j,j^{\prime }}^0(i{\omega }_n{)}^{-1}=G_{j,j^{\prime }}(i{\omega }_n{)}^{-1}+{\Sigma }_{j,j^{\prime }}(i{\omega }_n).$$

$$Z=\sum _{\{\eta \}}{Z}^{\{\eta \}},$$

$$Z^{\{\eta \}}=\int \mathcal{D}\chi \exp (-{\mathcal{S}}_{\text{eff}}^{\{\eta \}}).$$

$${\mathcal{S}}_{\text{eff}}^{\{\eta \}}=$$

$$+2T\sum _{j,j^{\prime },n\ge 0}{\chi }_{j,-{\omega }_n}{\mathcal{H}}_{j,j^{\prime }}^{\{\eta \}}{\chi }_{j^{\prime },{\omega }_n}.$$

$${\left[(G^{\{\eta \}}(i{\omega }_n){)}^{-1}\right]}_{j,j^{\prime }}={\left[({\mathcal{G}}^0(i{\omega }_n){)}^{-1}\right]}_{j,j^{\prime }}-2{\mathcal{H}}_{j,j^{\prime }}^{\{\eta \}}.$$

$$G_{j,j^{\prime }}^{\mathrm{imp}}(i{\omega }_n)=\sum _{\{\eta \}}P(\{\eta \})G_{j,j^{\prime }}^{\{\eta \}}(i{\omega }_n),$$

$$P(\{\eta \})=Z^{\{\eta \}}/\sum _{\{\eta \}}{Z}^{\{\eta \}}.$$

$$Z^{\{\eta \}}=\prod _{n\ge 0}\text{det}[-G^{\{\eta \}}(i{\omega }_n)].$$

$${\Sigma }_{j,j^{\prime }}(i{\omega }_n)={\left[({\mathcal{G}}^0(i{\omega }_n){)}^{-1}\right]}_{j,j^{\prime }}-{\left[(G^{\mathrm{imp}}(i{\omega }_n){)}^{-1}\right]}_{j,j^{\prime }}.$$

$$G(i{\omega }_n)=G^{\mathrm{imp}}(i{\omega }_n),$$

$$⟨S_j^z(\tau )S_{j^{\prime }}^z⟩=\sum _{\{\eta {\}}^{\prime },{\eta }_{r_0}=\pm 1}P(\{\eta {\}}^{\prime },{\eta }_{r_0})⟨S_j^z(\tau )S_{j^{\prime }}^z{⟩}^{\{\eta {\}}^{\prime }},$$

$${\mathcal{S}}_{\text{eff}}^{\{\eta {\}}^{\prime }}={\mathcal{S}}_{\mathrm{hyb}}^{\{\eta {\}}^{\prime }}+{\mathcal{S}}_{\mathrm{local}},$$

$${\mathcal{S}}_{\mathrm{hyb}}^{\{\eta {\}}^{\prime }}=$$

$${\mathcal{S}}_{\mathrm{local}}=$$

$${\tilde{G}}_{j,j^{\prime }}^{\{\eta \}}(i{\omega }_n)=G_{j,j^{\prime }}^{\{\eta \}}(i{\omega }_n).$$

$${\Delta }_{j,j^{\prime }}^{\{\eta {\}}^{\prime }}(i{\omega }_n)=[{\tilde{G}}^{\{\eta \}}(i{\omega }_n){]}_{j,j^{\prime }}^{-1}-(i{\omega }_n-2{\mathcal{H}}_{j,j^{\prime }}^{\{\eta \}}).$$

$${\Delta }_{j,j^{\prime }}^{\{\eta {\}}^{\prime }}(\tau )=\frac{T}{2}\sum _n{e}^{-i{\omega }_n\tau }{\Delta }_{j,j^{\prime }}^{\{\eta {\}}^{\prime }}(i{\omega }_n).$$

$$\frac{Z}{Z_{\mathrm{local}}}$$

$$=\sum _{d,i_0,...,i_{2d-1}}{\int }_0^{\beta }d{\tau }_0...{\int }_0^{\beta }d{\tau }_{2d-1}\frac{1}{d!}⟨{\chi }_{i_0}({\tau }_0)...{\chi }_{i_{2d-1}}({\tau }_{2d-1}){⟩}_{\mathrm{local}}\mathrm{Pf}({\hat{\Delta }}^{\{\eta {\}}^{\prime }}(d,i_0,{\tau }_0,...,i_{2d-1},{\tau }_{2d-1})),$$

$$⟨\mathcal{A}{⟩}_{\mathrm{local}}=\frac{\int \mathcal{D}\chi \mathcal{A}{e}^{-{\mathcal{S}}_{\mathrm{local}}}}{\int \mathcal{D}\chi e^{-{\mathcal{S}}_{\mathrm{local}}}}.$$

$${\hat{\Delta }}^{\{\eta {\}}^{\prime }}(d,i_0,{\tau }_0...,i_{2d-1},{\tau }_{2d-1}{)}_{m,n}={\Delta }_{i_m,i_n}^{\{\eta {\}}^{\prime }}({\tau }_m-{\tau }_n).$$

$$\frac{\mathrm{Pf}({\hat{\Delta }}^{\{\eta {\}}^{\prime }}(d,i_0,{\tau }_0,...,i_{2d-1},{\tau }_{2d-1}))}{\mathrm{Pf}({\hat{\Delta }}^{\{\eta {\}}^{\prime }}(d+1,i_0,{\tau }_0,...,i_{2d+1},{\tau }_{2d+1}))},$$

$$⟨S_j^z(\tau )S_{j^{\prime }}^z{⟩}^{\{\eta {\}}^{\prime }}=\frac{Z_{\mathrm{local}}}{Z}\sum _{d,i_0,...,i_{2d-1}}{\int }_0^{\beta }d{\tau }_0...{\int }_0^{\beta }d{\tau }_{2d-1}\frac{1}{d!}$$

$$\times \mathrm{Pf}({\hat{\Delta }}^{\{\eta {\}}^{\prime }}(d,i_0,...,i_{2d-1},{\tau }_0,...,{\tau }_{2d-1})).$$

$$\frac{⟨T_{\tau }{\chi }_{i_0}({\tau }_0)...{\chi }_{i_{2d-1}}({\tau }_{2d-1})S_j^z(\tau )S_{j^{\prime }}^z{⟩}_{\mathrm{local}}}{⟨T_{\tau }{\chi }_{i_0}({\tau }_0)...{\chi }_{i_{2d-1}}({\tau }_{2d-1}){⟩}_{\mathrm{local}}}.$$

Full text

††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 $e^{2}/h$ ($e$ is the elementary charge and $h$ 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($T$) 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 $S=1/2$. 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 $S=1/2$ 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 $S=1/2$ 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 $T$ scales; in each crossover, itinerant Majorana fermions and localized gauge fluxes release their entropy, a half of $\log 2$ per site Nasu2014 ; Nasu2015 . Also in the ground state, the dynamical spin structure factor $S(\mathbf{q},\omega)$ 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, $\alpha$-RuCl3 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 $\alpha$-RuCl3 Sandilands2015 . Furthermore, the $T$ 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 $T$ by a newly developed numerical technique, the Majorana dynamical mean-field method Yoshitake2016 . Indications of the spin fractionalization were identified in the $T$ dependences of $S(\mathbf{q},\omega)$, the relaxation rate in the nuclear magnetic resonance (NMR), $1/T_{1}$, and the magnetic susceptibility $\chi$. 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 $S(\mathbf{q},\omega)$, $1/T_{1}$, and $\chi$ 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 $T$ 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 $T$ and $\omega$ dependence of spin correlations and the $T$ 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

[TABLE]

where $p=x$, $y$, and $z$, and the sum of ${\langle j,j^{\prime}\rangle}_{p}$ is taken for the nearest-neighbor (NN) sites on three inequivalent bonds of the honeycomb lattice, as indicated in Fig. 1(a); $S_{j}^{p}$ is the $p$ component of the $S=1/2$ spin at site $j$. Hereafter, we denote the average of $J_{p}$ as $J$ and set the energy scale as $\sum_{p}|J_{p}|=3$, i.e., $J=1$, and parametrize the anisotropy of the exchange coupling constants as $J_{x}=J_{y}=\pm\alpha$ and $J_{z}=\pm(3-2\alpha)$, 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: $\langle S_{j}^{p}S_{j^{\prime}}^{p}\rangle$ are nonzero only for the NN sites $j,j^{\prime}$ on the $p$ bonds as well as the same site $j=j^{\prime}$ Baskaran2007 . Hereafter, we denote the NN correlations as $\langle S_{j}^{p}S_{j^{\prime}}^{p}\rangle_{\rm{NN}}$. There are two types of QSL phases depending on the anisotropy in $J_{p}$: one is a gapless QSL realized in the region with $0.75\leq\alpha\leq 1.5$ including the isotropic point $\alpha=1$ ($J_{x}=J_{y}=J_{z}=\pm 1$), while the other is gapful for $0\leq\alpha<0.75$. 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 $S=1/2$ spin Kitaev2006 . In this method, the Hilbert space in the original spin representation, $2^{N}$, is extended to $4^{N}$ in the Majorana fermion representation ($N$ 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 $J_{x}$ and $J_{y}$ bonds. Then, by introducing two Majorana fermions $c_{j}$ and $\bar{c}_{j}$ for the spinless fermions, the Hamiltonian in Eq. (1) is rewritten as

[TABLE]

where the sum over $(j,j^{\prime})_{p}$ is taken for the NN sites on a $p$ bond with $j<j^{\prime}$. $\eta_{r}=i\bar{c}_{j}\bar{c}_{j^{\prime}}$ is defined on each $z$ bond connecting $j$ and $j^{\prime}$ sites ($r$ is the index of the $z$ bond). Here, $\eta_{r}$ is considered as a $Z_{2}$ variable taking $\pm 1$, as $\eta_{r}$ commutes with the total Hamiltonian as well as with other $\eta_{r^{\prime}}$ and as $\eta_{r}^{2}=1$. Thus, the model in Eq. (2) describes itinerant Majorana fermions $\{c_{j}\}$ (called matter fermions) coupled to the $Z_{2}$ variables $\{\eta_{r}\}$ (called gauge fluxes). The ground state is given by all $\eta_{r}=1$, giving QSLs with gapless or gapful excitations depending on $\alpha$, as in the original Kitaev’s solution.

In the present numerical study at finite $T$, 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 $z$ 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 $\eta_{r}$ develop at low $T$. 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

[TABLE]

where $\gamma$ and $\gamma^{\prime}$ are the indices for the clusters, and $j$ and $j^{\prime}$ denote the sites in each $N_{c}$-site cluster. The coefficient $1/2$ 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

[TABLE]

where $\omega_{n}=(2n+1)\pi T$ is the Matsubara frequency ($n$ is an integer, and the Boltzmann constant $k_{\rm B}$ and the reduced Planck constant $\hbar$ are set to unity), $\Sigma(\mathbf{k},i\omega_{n})$ is the self-energy, and $\mathcal{H}^{0}(\mathbf{k})$ is the Fourier transform of $\mathcal{H}^{0}_{\gamma,j;\gamma^{\prime},j^{\prime}}$ in Eq. (3) given by the matrix:

[TABLE]

where $\mathbf{r}_{\gamma}$ is the coordinate of the cluster $\gamma$.

Following the spirit of the DMFT Metzner1989 ; Georges1996 , we omit the $\mathbf{k}$ dependence of the self-energy: $\Sigma(\mathbf{k},i\omega_{n})=\Sigma(i\omega_{n})$. In this approximation, local Green’s function is defined within a cluster as

[TABLE]

where $N^{\prime}$ is the number of clusters in the whole lattice ($N=N_{c}N^{\prime}$), and $j$ and $j^{\prime}$ denotes the sites in the cluster. The Weiss function is introduced to take into account the correlation effects in other clusters as

[TABLE]

In order to take into account the interaction $\mathcal{H}^{\{\eta\}}$ 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

[TABLE]

where

[TABLE]

Here, the sum of $\{\eta\}$ in Eq. (8) runs over all possible configurations of $\{\eta\}$, and $\mathcal{D}\chi=\prod_{j,n}d\chi_{j,\omega_{n}}$ in Eq. (9); $\chi_{j,\omega_{n}}$ is the Grassmann number corresponding to the Majorana operator $c_{j}$ (more precisely, $c_{j}/\sqrt{2}$ following the notation in Ref. Nilsson2013 ). The effective action is given by

[TABLE]

For a given configuration of $\{\eta\}$, 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

[TABLE]

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

[TABLE]

where $P(\{\eta\})$ is the statistical weight for the configuration $\{\eta\}$ given by

[TABLE]

$Z^{\{\eta\}}$ is obtained from Green’s functions as

[TABLE]

We note that $G^{\rm{imp}}(i\omega_{n})$ is obtained exactly by computing $G^{\{\eta\}}(i\omega_{n})$ and $P(\{\eta\})$ for all $2^{N_{c}/2}$ configurations of $\{\eta\}$ in the $N_{c}$-site cluster Udagawa2012 . The self-energy for the impurity problem is obtained as

[TABLE]

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

[TABLE]

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 $T$ dependences of static quantities of the Kitaev model, such as the specific heat and the equal-time spin correlations $\langle S_{j}^{p}S_{j^{\prime}}^{p}\rangle$. 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 $2^{N_{c}/2}$ configurations in Eq. (12) enables the exact calculations for the given cluster. Furthermore, the cluster-size dependence is sufficiently small at all the $T$ 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 $\alpha=1.0$ 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 $\langle S^{p}_{j}(\tau)S^{p}_{j^{\prime}}\rangle$ ($\tau$ 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 $\langle S^{p}_{j}(\tau)S^{p}_{j^{\prime}}\rangle$, we need to take into account the imaginary-time evolution of $\{\bar{c}\}$ that compose the conserved quantities $\{\eta\}$, e.g., $S^{z}_{j}(\tau)=\pm i\chi_{j}(\tau)\bar{c}_{j}(\tau)/\sqrt{2}$; 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, $\langle S^{z}_{j}(\tau)S^{z}_{j^{\prime}}\rangle$ on an $r_{0}$ bond is calculated as

[TABLE]

where $\{\eta\}^{\prime}$ represents the configurations of $\eta_{r}$ except for $\eta_{r_{0}}$ on the $r_{0}$ bond. $P(\{\eta\}^{\prime},\eta_{r_{0}})$ is obtained from the converged solution of the Majorana CDMFT in Sec. II.2. $\langle S^{z}_{j}(\tau)S^{z}_{j^{\prime}}\rangle^{\{\eta\}^{\prime}}$ is the dynamical spin correlation on the $r_{0}$ bond calculated by the CTQMC method for each configuration $\{\eta\}^{\prime}$. The sum of $\{\eta\}^{\prime}$ runs over all possible configurations of $\{\eta\}^{\prime}$ within the cluster. Note that Eq. (17) is derived from the fact that $S^{z}_{j}$ commutes with $\eta_{r}$ in $\{\eta\}^{\prime}$, whereas it does not commute with $\eta_{r_{0}}$. Thus, for a given $\{\eta\}^{\prime}$, the interaction lies only on the $r_{0}$ 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 $\chi_{j,\omega_{n}}$ whose $j$ does not belong to the $r_{0}$ bond. Then, we obtain

[TABLE]

where

[TABLE]

and $j,j^{\prime}$ in Eqs. (19) and (20) are the sites on the $r_{0}$ bond; $\beta=1/T$ is the inverse temperature. In Eq. (19), the hybridization function $\Delta_{j,j^{\prime}}^{\{\eta\}^{\prime}}(\tau)$ is calculated from $G_{j,j^{\prime}}^{\{\eta\}}(i\omega_{n})$ in the converged solution of CDMFT as follows. Let us define the matrix $\tilde{G}^{\{\eta\}}(i\omega_{n})$ as a $2\times 2$ submatrix of $G^{\{\eta\}}(i\omega_{n})$, as

[TABLE]

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

[TABLE]

Note that $\Delta_{j,j^{\prime}}^{\{\eta\}^{\prime}}(i\omega_{n})$ does not depend on $\{\eta\}$, which is straightforwardly shown by the matrix operations in the right hand side. Converting Eq. (22) to the imaginary-time representation, we obtain

[TABLE]

Given Eqs. (18)-(20), the partition function of the system is expanded in terms of $\mathcal{S}_{\rm{hyb}}^{\{\eta\}^{\prime}}$ as

[TABLE]

where $Z_{\rm{local}}=\int\mathcal{D}\chi e^{-\mathcal{S}_{\rm{local}}}$ is the partition function for the two sites described by $\mathcal{S}_{\rm{local}}$, and $\langle\mathcal{A}\rangle_{\rm{local}}$ represents the expectation value in the two-site problem as

[TABLE]

In the second line of Eq. (24), $d$ is the order of $\mathcal{S}_{\rm{hyb}}^{\{\eta\}^{\prime}}$ in the expansion of $e^{-\mathcal{S}_{\rm{hyb}}^{\{\eta\}^{\prime}}}$, Pf($M$) is the Pfaffian of skew-symmetric matrix $M$, and $\hat{\Delta}^{\{\eta\}^{\prime}}(d,i_{0},\tau_{0},...,i_{2d-1},\tau_{2d-1})$ is a $2d\times 2d$ matrix, whose $(m,n)$ element is given by

[TABLE]

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 $(d,i_{0},\tau_{0},...,i_{2d-1},\tau_{2d-1})$ 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 $d$ as $(d,i_{0},\tau_{0},...,i_{2d-1},\tau_{2d-1})$ to $(d+1,i_{0},\tau_{0},...,i_{2d-1},\tau_{2d-1},i_{2d},\tau_{2d},i_{2d+1},\tau_{2d+1})$ by adding $(i_{2d},\tau_{2d}),(i_{2d+1},\tau_{2d+1})$. 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 $d$, the ratio is calculated by adding two rows and columns in the matrix $\hat{\Delta}^{\{\eta\}^{\prime}}$ as

[TABLE]

whose calculation cost is in the order of $d^{2}$ by using the fast update algorithm. On the other hand, in Eq. (24), $\langle{T_{\tau}\chi_{i_{0}}(\tau_{0})...\chi_{i_{2d-1}}(\tau_{2d-1})}\rangle_{\rm{local}}$ 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 $\{\eta\}^{\prime}$, $\langle S^{z}_{j}(\tau)S^{z}_{j^{\prime}}\rangle^{\{\eta\}^{\prime}}$ in Eq. (17), is calculated as

[TABLE]

For the MC sampling, we need to evaluate

[TABLE]

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 $\alpha=1.0$, $\langle S^{p}_{j}(\tau)S^{p}_{j^{\prime}}\rangle$ for $p=x,y$ are equivalent to $\langle S^{z}_{j}(\tau)S^{z}_{j^{\prime}}\rangle$. Meanwhile, for the anisotropic case, we compute $\langle S^{p}_{j}(\tau)S^{p}_{j^{\prime}}\rangle$ for $p=x,y$ by the same technique described above with using the spin rotations $\{S^{x},S^{y},S^{z}\}\to\{S^{y},S^{z},S^{x}\}$ or $\{S^{x},S^{y},S^{z}\}\to\{S^{z},S^{x},S^{y}\}$.

In the CTQMC calculations in Sec. III, for each configuration $\{\eta\}^{\prime}$, we typically run $10^{7}$ MC steps and perform the measurements at every 20 steps, after $10^{5}$ 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, $\langle S^{p}_{j}(\tau)S^{p}_{j^{\prime}}\rangle$. 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 $\omega$, we need to inversely solve the equation given by the generic form $g(\tau)=\int d\omega\rho(\omega)e^{-\omega\tau}$. In our problem, $g(\tau)$ and $\rho(\omega)$ correspond to the dynamical spin correlations as functions of imaginary time $\tau$ and real frequency $\omega$: $g(\tau)=\langle S^{p}_{j}(\tau)S^{p}_{j^{\prime}}\rangle$ and $\rho(\omega)=S^{p}_{j,j^{\prime}}(\omega)$. In the following calculations, we utilize the Legendre polynomial expansion following Refs. Boehnke2011 ; Levy2017 :

[TABLE]

where $P_{m}(x)$ is the $m$th Legendre polynomials and $x(\tau)=2\tau/\beta-1$. Then, the inverse problem is given by

[TABLE]

where

[TABLE]

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 $\rho(\omega)$ to $\rho_{l}=\rho(\omega_{l})$, and determine $\rho_{l}$ to minimize the function

[TABLE]

where $\zeta$ and $\delta$ are the coefficients described below, and $C$ is a variance-covariance matrix of $g_{m}$; $\tilde{g}_{m}=\sum_{l}\Delta\omega\rho_{l}K_{m}(\omega_{l})$. We take the Legendre expansion up to $50$th order and $\Delta\omega=0.01125$ in the following calculations. In Eq. (33), $\rho^{(0)}_{l}$ is the advance estimate of $\rho_{l}$, 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 $\mathcal{F}$ is equivalent to the least squares method. The least squares method is unstable, as $g_{m}$ is rather insensitive to a change of $\rho_{l}$. The second term, called the entropy term, stabilizes the minimization process. In the following calculations, we set $\zeta=625$ to sufficiently take into account the effect of the entropy term, where the value of $\delta$ is determined self-consistently in each MEM calculation based on the maximum likelihood estimation, called the classical MEM Jarrell1996 (typically, $\delta\simeq 1$-$10$). We note that the deviations of $\tilde{g}_{m}$ from $g_{m}$ are typically comparable to the statistical errors in the CTQMC calculations. In the following results, we estimate the errors of $\rho(\omega)$ by the standard deviation between the data for $\zeta=100$, $625$ and $10000$ in the range where the MEM retains the precision.

In the MEM, $\rho(\omega)$ should be positive for all $\omega$. In our problem, the onsite correlation $S^{p}_{j,j}(\omega)$ satisfies this condition automatically, whereas $S^{p}_{j,j^{\prime}}(\omega)$ for the NN sites $j,j^{\prime}$ on the $p$ bond, which is denoted by $S^{p}_{\rm{NN}}(\omega)$ hereafter, can be negative. (Note that all the further-neighbor correlations beyond the NN sites vanish in the Kitaev model Baskaran2007 .) To obtain $S^{p}_{\rm{NN}}(\omega)$ properly, we calculate $S^{p}_{j,j}(\omega)+2S^{p}_{j,j^{\prime}}(\omega)+S^{p}_{j^{\prime},j^{\prime}}(\omega)$, which is positive definite for all $\omega$, and subtract the onsite contributions note . The accuracy of $S^{p}_{j,j^{\prime}}(\omega)$ obtained by the MEM are examined in Appendix B in the one-dimensional limit with $\alpha=1.5$, where $S^{p}_{j,j^{\prime}}(\omega)$ 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 $J_{x}$, $J_{y}$, and $J_{z}$. By comparing the results with those by the QMC method Nasu2015 , we confirm that the CDMFT is valid in the $T$ range above the artificial critical temperature close to the low-$T$ 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 $T$ 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 $C_{v}$ and equal-time spin correlations for NN pairs on the $p$ bonds, $\langle S_{j}^{p}S_{j^{\prime}}^{p}\rangle_{\rm{NN}}$, obtained by the Majorana CDMFT, with those by QMC in Ref. Nasu2015 . The data are calculated for the FM case with bond asymmetry: $\alpha=0.8$ ($J_{x}=J_{y}=0.8$ and $J_{z}=1.4$) and $\alpha=1.2$ ($J_{x}=J_{y}=1.2$ and $J_{z}=0.6$). While the data of $C_{v}$ are common to the FM and AFM cases, the sign of $\langle S_{j}^{p}S_{j^{\prime}}^{p}\rangle_{\rm{NN}}$ is reversed for the AFM case. Note that similar comparison was made for the isotropic case $\alpha=1.0$ ($J_{x}=J_{y}=J_{z}=1$) 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 $T_{\rm{L}}\simeq 0.012$ and $T_{\rm{H}}\simeq 0.375$ in the isotropic case. In the anisotropic cases, the low-$T$ crossover takes place at a lower $T$, i.e., $T_{\rm{L}}\simeq 0.0052$ for $\alpha=0.8$ and $T_{\rm L}\simeq 0.0075$ for $\alpha=1.2$, while the high-$T$ one is almost unchanged, i.e., $T_{\rm{H}}\simeq 0.375$. These behaviors are excellently reproduced by the Majorana CDMFT, except for the low-$T$ peak; the CDMFT results show a sharp anomaly at $\tilde{T}_{\rm c}\simeq 0.0063$ for $\alpha=0.8$ and $\tilde{T}_{\rm c}\simeq 0.013$ for $\alpha=1.2$. This is due to a phase transition by ordering of $\eta$, 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 $T$ range above the artificial phase transition temperature $\tilde{T}_{\rm c}$. Although they appear to be reproduced even below $\tilde{T}_{\rm c}$, there is a small anomaly at $\tilde{T}_{\rm c}$ associated with the artificial transition, while the QMC data smoothly change around $T_{\rm L}$. (Note that the appropriate sum of the NN spin correlations is nothing but the internal energy, and hence, the $T$ derivative corresponds to the specific heat.)

Thus, the comparison indicates that the Majorana CDMFT gives quantitatively precise results in the wide $T$ range above the artificial transition temperature $\tilde{T}_{\rm c}$: in the present cases with $\alpha=0.8$ and $1.2$, the CDMFT is reliable for $T\gtrsim 0.007$ and $T\gtrsim 0.014$, respectively. As discussed in the previous study Nasu2015 , the thermal fractionalization of quantum spins sets in below $T\simeq T_{\rm H}$, which is well above $\tilde{T}_{\rm c}$. Thus, the $T$ 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 $T$ 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 $S(\mathbf{q},\omega)$ at several $T$ for the FM case with $\alpha=1.0$, $0.8$, and $1.2$. $S(\mathbf{q},\omega)$ is calculated as

[TABLE]

where $S^{p}_{j,j^{\prime}}(\omega)$ is obtained by the MEM described in Sec. II.4 from the imaginary-time correlations $\langle S^{p}_{j}(\tau)S^{p}_{j^{\prime}}\rangle$ by CDMFT+CTQMC. As mentioned above, nonzero contributions in Eq. (35) come from only the onsite and NN-site components of $S^{p}_{j,j^{\prime}}(\omega)$; we present their $T$ dependences in Appendix C. The Brillouin zone and symmetric lines on which $S(\mathbf{q},\omega)$ is plotted are presented in Fig. 1(b). Although the results at $\alpha=1.0$ were shown in the previous study Yoshitake2016 , we present them (for a sightly different $T$ set) for comparison. We show the data at four temperatures: $T\simeq 2T_{\rm L}$, $\sqrt{T_{\rm L}T_{\rm H}}$, $0.64T_{\rm H}$, and $6.4T_{\rm H}$. Note that $T_{\rm L}\simeq 0.012$, $0.0052$, and $0.0075$ for $\alpha=1.0$, $0.8$, and $1.2$, respectively, while $T_{\rm H}\simeq 0.375$ for all the cases.

As shown in Fig. 3, at sufficiently high $T$ than $T_{\rm{H}}$, $S(\mathbf{q},\omega)$ does not show any significant $\mathbf{q}$ dependence for all $\alpha$ studied here; $S(\mathbf{q},\omega)$ shows only a diffusive response centered at $\omega\sim 0$, as shown in Figs. 3(j)-3(l). When lowering $T$ below $T_{\rm{H}}$, the diffusive weight is shifted to the positive $\omega$ region ranging up to above $\omega\sim J$ for all the cases, as shown in Figs. 3(g)-3(i). Simultaneously, a quasi-elastic component grows gradually at $\omega\sim 0$. Both the inelastic and the quasi-elastic components show a discernible $\mathbf{q}$ dependence; in particular, the latter increases the intensity around the $\Gamma$ point reflecting the FM interactions. While $S({\rm{K}}_{1},\omega)=S({\rm{K}}_{2},\omega)$ and $S({\rm{M}}_{1},\omega)=S({\rm{M}}_{2},\omega)$ for $\alpha=1.0$ from the symmetry, the quasi-elastic response is small (large) around the M1-K1 line compared to that around the M2-K2 line for $\alpha=0.8$ ($1.2$) because of the anisotropy.

When further lowering $T$ and approaching $T_{\rm L}$, the quasi-elastic component increases its intensity, while the inelastic response at $\omega\sim J$ does not change substantially. In particular, in the case of $\alpha=0.8$, the quasi-elastic component is sharpened and develops to a $\delta$-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 $T=0$, where the $\delta$-function peak appears due to the change of the parity between the ground state and the flux-excited state Knolle2014 (for the $\delta$-function peak, see also Fig. 16 in Appendix B). On the other hand, $S(\mathbf{q},\omega)$ at $\alpha=1.2$ 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 $\alpha=1.2$ are qualitatively similar to those for $\alpha=1.0$ in Figs. 3(d) and 3(a), except for different $\mathbf{q}$ dependence mentioned above.

Figure 4 shows the results for the AFM case. The overall $\omega$ dependence of $S(\mathbf{q},\omega)$ is similar to that for the FM case at all $T$: the diffusive response centered at $\omega\sim 0$ for $T\gtrsim T_{\rm{H}}$ [Figs. 4(j)-4(l)], the shift of the diffusive weight to the region of $\omega\sim J$ and the growth of a quasi-elastic component at $\omega\sim 0$ below $T_{\rm{H}}$ [Figs. 4(g)-4(i)], and the $\delta$-function like peak for $\alpha=0.8$ while approaching to $T_{\rm{L}}$ [Figs. 4(e) and 4(b)]. The similarity of the $\omega$ dependences of $S(\mathbf{q},\omega)$ between FM and AFM cases is partly understood by the relation $2S({\rm{K}}_{1},\omega)_{\rm{FM}}+S({\rm{K}}_{2},\omega)_{\rm{FM}}=2S({\rm{K}}_{1},\omega)_{\rm{AFM}}+S({\rm{K}}_{2},\omega)_{\rm{AFM}}$, which holds for $J_{x}=J_{y}$ [$S(\mathbf{q},\omega)_{\rm{FM}}$ and $S(\mathbf{q},\omega)_{\rm{AFM}}$ are $S(\mathbf{q},\omega)$ for the FM and AFM cases, respectively]. On the other hand, the $\mathbf{q}$ dependence is in contrast to the FM case: while the weight of the quasi-elastic response almost vanishes around the $\Gamma$ point, those on the zone boundary are enhanced in an almost opposite manner to the FM cases. In addition, the incoherent weight at $\omega\sim J$ also shows the opposite $\mathbf{q}$ dependence to the FM case: the weight is stronger around the $\Gamma$ point than that on the zone boundary. The opposite $\mathbf{q}$ dependences between the FM and AFM cases directly follow from the relation $S(\textbf{q},\omega)_{\rm{AFM}}=-S(\textbf{q},\omega)_{\rm{FM}}+(2/3)\sum_{p}S^{p}_{j,j}(\omega)$.

In order to show the $T$ dependences of $S(\mathbf{q},\omega)$ more explicitly, we present in Figs. 5-8 the $T$-$\omega$ plot of $S(\mathbf{q},\omega)$ at $\mathbf{q}=\Gamma$, ${\rm{K}}_{1}$, and K2 with the intensity profiles for the same set of $T$ used in Figs. 3 and 4. Figure 5 shows the result for the FM case at $\alpha=0.8$. The overall weight of $S(\mathbf{q},\omega)$ shifts from $\omega\sim 0$ to a large-$\omega$ region when the system is cooled down below $T\sim T_{\rm{H}}$. Below $T_{\rm{H}}$, quasi-elastic response gradually grows and develops to the $\delta$-function like peak. The peak intensity in $S(\Gamma,\omega)$ and $S({\rm{K}_{2}},\omega)$ is larger than that for $S({\rm{K}_{1}},\omega)$, reflecting the anisotropy of the interaction.

Figure 6 shows the corresponding plot for the AFM case at $\alpha=0.8$. In contrast to the FM case, the strong quasi-elastic response is seen for $\mathbf{q}={\rm{K}}_{1}$, which develops to the $\delta$-function like peak at low $T$. We note that the dip and shoulder like structures around $\omega=0$ in the intermediate $T$ for the result at $\mathbf{q}={\rm{K}}_{2}$ 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 $S^{p}_{j,j^{\prime}}(\omega)$ for the NN bonds by subtracting the onsite component $S^{p}_{j,j}(\omega)$ from $S^{p}_{j,j}(\omega)+S^{p}_{j,j^{\prime}}(\omega)$, both of which are obtained by the MEM. In the present case, as both of $S^{z}_{j,j}(\omega)$ and $S^{z}_{j,j}(\omega)+S^{z}_{j,j^{\prime}}(\omega)$ become large around $\omega=0$ due to the development of the $\delta$-function like peak, the relative error becomes large for $S({\rm{K}_{2}},\omega\sim 0)$, which may lead to artificial structures.

Figures 7 and 8 show the results at $\alpha=1.2$. As observed in Figs. 3 and 4, $S(\mathbf{q},\omega)$ for both the FM and AFM cases behave similarly to those at $\alpha=1.0$ Yoshitake2016 . In the anisotropic cases, however, the difference between $S({\rm{K}}_{1},\omega)$ and $S({\rm{K}}_{2},\omega)$ is obvious: the quasi-elastic peak for $S({\rm{K}}_{1},\omega)$ is larger (smaller) than that for $S({\rm{K}}_{2},\omega)$ in the FM (AFM) case.

As discussed in the previous study Yoshitake2016 , there is a relation between the static spin correlation and the average frequency of $S(\Gamma,\omega)$, $\bar{\omega}\equiv\int\omega S(\Gamma,\omega)d\omega/\int S(\Gamma,\omega)d\omega$, originating from the sum rule for $S(\mathbf{q},\omega)$. $T$ dependences of $\bar{\omega}$ are shown by white dashed curves in Figs. 5(b), 6(b), 7(b), and 8(b). In all cases, $\bar{\omega}$ is nearly zero for sufficiently high $T$, but it grows at $T\sim T_{\rm{H}}$ and becomes almost independent of $T$ for $T\lesssim T_{\rm{H}}$. These $T$ dependences are similar to those of the static spin correlation between the NN sites shown in Fig. 2.

III.3 NMR relaxation rate

Figure 9 shows the NMR relaxation rate $1/T_{1}$ obtained by the CDMFT+CTQMC method. While the results at $\alpha=1.0$ were presented in the previous study Yoshitake2016 , we present them for comparison in Fig. 9(a). $1/T_{1}$ in the magnetic field applied to the $z$ direction, which is denoted by $1/T^{z}_{1}$, is given by

[TABLE]

where $A_{\textbf{q}}$ is the hyperfine coupling constant, $\chi^{\perp}(\textbf{q},\omega)$ is the dynamical susceptibility for the spin component perpendicular to the magnetic field direction, and $\omega_{0}$ is the resonance frequency in the NMR measurement. The dynamical susceptibility $\chi(\textbf{q},\omega)$ is related with the dynamical spin structure factor through the fluctuation-dissipation theorem, as

[TABLE]

In the NMR experiments, $\omega_{0}$ is in general negligibly small compared to the typical energy scale of the system, $J$. Thus, by taking the limit of $\omega_{0}\rightarrow 0$ in Eq. (36) and using Eq. (37), we obtain

[TABLE]

where the coefficients $a_{0}$, $a_{1}$, $a_{2}$, and $a_{3}$ are determined by $A_{\textbf{q}}$. The similar equations are obtained for $1/T^{x}_{1}$ and $1/T^{y}_{1}$ by the cyclic permutation of $x,y,z$ ($1/T_{1}^{x}=1/T_{1}^{y}$ for the present cases from the symmetry). Because $A_{\textbf{q}}$ depends on the details of the system, we here compute the onsite and NN-site components of $1/T_{1}$ separately with omitting the coefficients: the onsite components are calculated as

[TABLE]

while the NN-site ones are

[TABLE]

where the sign is $+(-)$ for the FM (AFM) case. We note that, in the anisotropic cases $\alpha\neq 1.0$, the NN-site $1/T_{1}^{x}$ is not simply given by the sum in Eq. (42): it will be given by a linear combination of $S^{y}_{\rm{NN}}(\omega=0)$ and $S^{z}_{\rm{NN}}(\omega=0)$ with appropriate coefficients determined by $A_{\mathbf{q}}$. Such a linear combination, however, can be constructed from our data for Eqs. (41) and (42) by noting that $1/T_{1}^{z}=2S^{y}_{\rm{NN}}(\omega=0)$ for $J_{x}=J_{y}$. Hence, we present the results by Eqs. (41) and (42) in Fig. 9 for simplicity.

As shown in Fig. 9, for all cases, the onsite component of $1/T_{1}^{p}$ is nonzero and almost $T$ independent above $T\sim T_{\rm H}$, as expected for the conventional paramagnets Moriya1956 . On the other hand, the NN-site component is zero in the high-$T$ limit and increases as decreasing $T$. This behavior corresponds to the development of NN-site static spin correlations shown in Sec. III.1, as they have a relation through the sum rule, $\int S_{j,j^{\prime}}^{p}(\omega)d\omega=\langle S_{j}^{p}S_{j^{\prime}}^{p}\rangle$.

When lowering $T$ below $T_{\rm H}$, $1/T_{1}^{x}$ for $\alpha=0.8$ substantially increases, as shown in Fig. 9(b). The enhancement is much larger than the case of $\alpha=1.0$ in Fig. 9(a). This is due to the evolution of the $\delta$-function like peak in $S^{z}(\mathbf{q},\omega)$ discussed in Sec. III.2. In contrast, $S^{x}(\mathbf{q},\omega)$ and $S^{y}(\mathbf{q},\omega)$ do not develop such $\delta$-function like peaks, and hence, $1/T_{1}^{z}$ does not show enhancement unlike $1/T_{1}^{x}$. While further decreasing $T$, $1/T_{1}^{x}$ shows a peak slightly above $T_{\rm{L}}$. The decrease at low $T$ reflects a spin gap originating from the nonzero flux gap in the ground state Kitaev2006 . On the other hand, the onsite and NN-site components of $1/T_{1}^{z}$ are both suppressed below $T\sim T_{\rm H}$, after showing a plateau and broad peak, respectively. The suppression of $1/T_{1}^{z}$ is due to an increase of energy cost for a spin flip on the strong $z$ bond under the well-developed static spin correlations between NN sites in this $T$ range. Actually, the energy cost is represented by the average frequency of $S_{j,j}^{x}(\omega)$ as there is a relation

[TABLE]

On the other hand, $\bar{\omega}^{x}_{\rm{onsite}}$ is also written as $\bar{\omega}^{x}_{\rm{onsite}}=(J_{y}\langle S_{j}^{y}S_{j^{\prime}}^{y}\rangle_{\rm{NN}}+J_{z}\langle S_{j}^{z}S_{j^{\prime}}^{z}\rangle_{\rm{NN}})/2$ by the sum rule Yoshitake2016 . Thus, the energy cost becomes large below $T_{\rm{H}}$ according to the growth of $\langle S_{j}^{z}S_{j^{\prime}}^{z}\rangle_{\rm{NN}}$.

In contrast, as shown in Fig. 9(c), $T$ dependence of $1/T_{1}^{p}$ at $\alpha=1.2$ is similar to that at $\alpha=1.0$ in Fig. 9(a). Both $1/T_{1}^{x}$ and $1/T_{1}^{z}$ increase below $T_{\rm{H}}$ while decreasing $T$, in contrast to the case with $\alpha=0.8$. For $\alpha=1.2$, however, $1/T_{1}^{z}$ is larger than $1/T_{1}^{x}$, reflecting the stronger interactions on the $x$ and $y$ bonds than the $z$ bond. On further decreasing $T$, $1/T_{1}^{p}$ at $\alpha=1.2$ also shows the peak structure slightly above $T_{\rm{L}}$ and then decreases, as expected from the finite flux gap in the ground state.

Although the system is described by free Majorana fermions coupled to localized gauge fluxes, the NMR relaxation rate does not obey the Korringa law, $1/(T_{1}T)\sim{\rm constant}$, which is expected for free fermion systems. This is natural because the spin-flip excitation in the NMR process is a composite of both itinerant matter fermions and localized gauge fluxes. Nonetheless, for comparison to forth-coming experiments, we plot the Korringa ratio as a function of $T$ in Appendix D.

III.4 Magnetic susceptibility

Figures 10 and 11 show the $T$ dependences of the magnetic susceptibility $\chi^{p}$ for the FM and AFM cases, respectively. $\chi^{p}$ at $\alpha=1.0$ presented in the previous study Yoshitake2016 , are also presented in Fig. 10(a) and Fig. 11(a) for comparison. $\chi^{p}$ is calculated from the imaginary-time spin correlations as

[TABLE]

Note that this is obtained without the MEM. In all the cases, at sufficiently high $T$ compared to the dominant $J_{p}$, $J_{p}^{\rm max}={\rm max}\{J_{p}\}$, $\chi^{p}$ obeys the Curie-Weiss law,

[TABLE]

which is obtained by the standard mean-field approximation in the original spin representation. While decreasing $T$, $\chi^{p}$ shows a deviation from $\chi^{p}_{\rm{CW}}$ below $T\sim J_{p}^{\rm max}$.

Among the results, $\chi^{x}$ for the FM case and $\chi^{z}$ for the AFM case at $\alpha=0.8$ show peculiar $T$ dependences at low $T$. The former largely deviates from the Curie-Weiss behavior and saturates to a small nonzero value, as shown in Fig. 10(b) note1 . Meanwhile, the latter shows a broad hump at $T\sim T_{\rm H}$ and decreases as lowering $T$, as shown in Fig. 11(b). These $T$ dependences are qualitatively understood by considering a two-site dimer model on the $z$ bond obtained by setting $J_{x}=J_{y}=0$. The dimer model gives the analytical forms for the magnetic susceptibility as

[TABLE]

The results are plotted by the dashed-dotted curves in Figs. 10(b) and 11(b). $\chi_{\rm dimer}^{x}$ for the FM case almost saturates around $T\sim J_{z}/4$, as the dominant $J_{z}$ interaction suppresses the magnetization in the $x$ direction. This accounts for the behavior of $\chi^{x}$ in Fig. 10(b) qualitatively. Meanwhile, $\chi_{\rm dimer}^{z}$ also well reproduces a hump at $T\sim 0.5$ in $\chi^{z}$ for the AFM case in Fig. 11(b); $\chi^{z}$ remains nonzero down to low $T$ as nonzero $J_{x}$ and $J_{y}$ smear out the dimer gap.

In the case of $\alpha=1.2$, $T$ dependences of $\chi^{p}$ shown in Figs. 10(c) and 11(c) are similar to those for $\alpha=1.0$ in the previous study Yoshitake2016 replotted in Figs. 10(a) and 11(a), respectively; on decreasing $T$, $\chi^{p}$ continues to increase down to $T\sim T_{\rm{L}}$ in the FM case, whereas $\chi^{p}$ shows broad peak at a higher $T$ in the AFM case. The effect of anisotropic $J_{p}$, however, is clearly observed: the stronger interactions on the $x,y$ bonds than the $z$ bond result in larger (smaller) $\chi^{x}$ than $\chi^{z}$ in the FM (AFM) case. In addition, the temperature of the broad peak of $\chi^{z}$ ($\chi^{x}$) in the AFM case shifts to a lower (higher) $T$ than that for $\alpha=1.0$.

IV Discussion

As pointed out in the previous study for the isotropic case by the authors Yoshitake2016 and confirmed also for the anisotropic cases in the present study, a remarkable feature in the Kitaev model is the dichotomy between the dynamical and static spin correlations; namely, the NMR relaxation rate $1/T_{1}^{p}$ and the magnetic susceptibility, both of which reflect the dynamical spin correlations, show substantial $T$ dependences below $T_{\rm H}$ (Figs. 9-11), even though the static spin correlations $\langle S_{j}^{p}S_{j^{\prime}}^{p}\rangle_{\rm{NN}}$ almost saturate to the $T=0$ values (Fig. 2). The dichotomy is unconventional behavior hardly seen in conventional insulating magnets. This might be a signature of the fractionalization of quantum spins, as $T_{\rm H}$ is the temperature where the fractionalization sets in as indicated in the specific heat and entropy Nasu2015 .

To examine the dichotomy in more detail, we calculate the $T$ dependences of $1/T_{1}^{p}$ and $\langle S_{j}^{p}S_{j^{\prime}}^{p}\rangle_{\rm{NN}}$ for two extreme cases by assuming the configuration of $\{\eta\}$ by hand. One is the flux-free state with all $\eta_{r}=+1$, which is realized in the ground state. The other is the state with completely random $\{\eta\}$, corresponding to the high-$T$ limit. For this purpose, we regard a single $z$ bond $r_{0}$ as the cluster in CDMFT, and take $P(\eta_{r_{0}}=1)=1$ and $P(\eta_{r_{0}}=-1)=0$ for the former uniform state, while $P(\eta_{r_{0}}=1)=P(\eta_{r_{0}}=-1)=1/2$ for the latter random state, in Eq. (12) of the self-consistent equation of CDMFT note2 .

Figure 12 shows the results. In all cases, $\langle S_{j}^{p}S_{j^{\prime}}^{p}\rangle_{\rm{NN}}$ for both uniform and random $\{\eta\}$ shows almost similar $T$ dependence to the CDMFT results, as shown in Figs. 12(a)-12(c). However, $1/T_{1}^{p}$ exhibits considerably different $T$ dependence. For instance, in the isotropic case with $\alpha=1.0$, although $1/T_{1}$ is almost $T$ independent for $T>T_{\rm H}$ for both uniform and random $\{\eta\}$ similar to the result by the CDMFT+CTQMC method in Ref. Yoshitake2016 , it shows different behavior below $T_{\rm H}$ between the two cases, as shown in Fig. 12(d). For the case with uniform $\{\eta\}$, $1/T_{1}$ decreases to zero after showing a small hump. The suppression at low $T$ reflects the flux gap $\Delta\simeq 0.065J$ in the flux-free state Kitaev2006 ; Knolle2014 . On the other hand, for the case with random $\{\eta\}$, $1/T_{1}$ monotonically increases while decreasing $T$ in the calculated $T$ range. Similar $T$ dependences of $1/T_{1}^{p}$ are obtained for $1/T_{1}^{x}$ at $\alpha=0.8$ and $1/T_{1}^{x,z}$ at $\alpha=1.2$, as shown in Figs. 12(e) and 12(f), respectively. We note that $1/T^{z}$ for $\alpha=0.8$ behaves differently; we will comment on this point in the end of this section.

The results clearly indicate that the peculiar $T$ dependences of $1/T_{1}$ found in the CDMFT+CTQMC results are closely related with fluctuations of the gauge fluxes $\{\eta\}$ composed of localized Majorana fermions $\{\bar{c}\}$ emergent from the spin fractionalization. As seen in the equal-time spin correlations shown in Figs. 12(a)-12(c), itinerant matter fermions develop their kinetic energy to the saturation at $T\sim T_{\rm H}$ (the equal-time spin correlations correspond to the kinetic energy of matter fermions). Due to the fractionalization, however, the localized gauge fluxes are still disordered even below $T_{\rm H}$ Nasu2015 , which results in the enhancement of $1/T_{1}$, as indicated in Figs. 12(d)-12(f). When approaching $T_{\rm L}$, $\{\eta\}$ are aligned in a coherent manner Nasu2015 , and hence, $1/T_{1}$ is rapidly suppressed at $T\sim T_{\rm L}$. Thus, the $T$ dependence of $1/T_{1}$ is qualitatively explained by the crossover from that for the random $\{\eta\}$ to the fully-aligned $\{\eta\}$ while decreasing $T$. The crossover occurs well below $T_{\rm H}$ and close to $T_{\rm L}$. Of course, as the original quantum spin is a composite of itinerant matter fermions and localized gauge fluxes, the spin-flip dynamics is a composite excitation. Nevertheless, our results indicate that the peculiar $T$ dependence of the NMR relaxation rate as well as the magnetic susceptibility is dominated by the emergent gauge fluxes from the fractionalization.

As noted above, $1/T_{1}^{z}$ for $\alpha=0.8$ behaves differently from others: $1/T_{1}^{z}$ for the random $\{\eta\}$ is smaller than that for the uniform $\{\eta\}$ at low $T$, as shown in Fig. 12(e). This is presumably because of the peculiar $T$ dependence of the density of states (DOS) for the itinerant matter fermions at $\alpha=0.8$. In the gapless QSL region for $0.75\leq\alpha\leq 1.5$ but close to the gapless-gapful boundary at $\alpha=0.75$, the DOS opens a gap as $\{\eta\}$ are thermally disordered by raising $T$ Nasu2015 . Thus, the DOS for matter fermions is gapless for the uniform $\{\eta\}$, while gapped for the random $\{\eta\}$. As spin excitations by $S_{j}^{x}$ and $S_{j}^{y}$ are composite excitations of both itinerant matter fermions and localized gauge fluxes, the gap in the DOS for matter fermions suppresses $1/T_{1}^{z}$ for the random case compared to the uniform one. Since $\{\eta\}$ are aligned uniformly below $T_{\rm L}$, we expect that $1/T_{1}^{z}$ shows an abrupt increase while decreasing $T$ through $T_{\rm L}$. This indicates that while a rapid change of $1/T_{1}$ when approaching $T_{\rm L}$ is yielded by the coherent alignment of $\{\eta\}$, either increase or decrease of $1/T_{1}$ at $T_{\rm L}$ may be affected by the itinerant matter fermions.

V Summary

To summarize, we have presented numerical results for spin dynamics of the Kitaev model with the anisotropy in the bond-dependent coupling constants. We calculated the experimentally-measurable quantities, the dynamical spin structure factor $S(\bf{q},\omega)$, the NMR relaxation rate $1/T_{1}$, and the magnetic susceptibility $\chi$, in the wide $T$ range including the peculiar paramagnetic region where quantum spins are fractionalized. The results have been obtained by the Majorana CDMFT+CTQMC method, which were developed by the authors previously Yoshitake2016 ; we gave detailed descriptions of the method, including the MEM for analytical continuation. We also confirmed the Majorana CDMFT is precise enough in the range of $T$ and anisotropy that we investigated in the present study.

We found that the Kitaev model exhibits unconventional behaviors in spin dynamics in the finite-$T$ paramagnetic state in proximity to the QSL ground state. The prominent feature is the dichotomy between static and dynamical spin correlations as a consequence of the spin fractionalization. The dichotomy appears clearly in the increase of $1/T_{1}$ below $T_{\rm H}$ where the fractionalization sets in, despite the saturation of static correlations. Similar behavior was also seen in the isotropic case in the previous study Yoshitake2016 . Our results suggest that the dichotomy is found universally in the fractionalized paramagnetic region irrespective of the anisotropy in the system.

On the other hand, we also clarified interesting behaviors that depend on the anisotropy at low $T$. When one of the three bond-dependent interactions is stronger than the other two, the spin dynamics shows peculiar $T$ and energy dependences distinct from those in the isotropic coupling case as follows. As lowering $T$, $S(\bf{q},\omega)$ develops a $\delta$-function like peak, which is well separated from the incoherent continuum. $1/T_{1}$ monotonically decreases in the spin component for the stronger bond. $\chi$ increases and saturates to a nonzero value for the spin component for the weaker bonds, while it shows hump and then decreases for the stronger-bond component in the antiferromagnetic case. We also showed that the peculiar $T$ dependences of $\chi$ are qualitatively explained by the two-site dimer model. In contrast, when the anisotropy is opposite, i.e., when the two types of bonds become stronger, the results are qualitatively unchanged from those for the isotropic case, while the effect of anisotropy is obvious in the $\bf{q}$ dependence in $S(\bf{q},\omega)$ and the different components in $1/T_{1}$ and $\chi$.

Our results will stimulate further experimental and theoretical analyses of candidate materials for the Kitaev QSLs. As most of the materials are assumed to be anisotropic in the exchange constants Choi2012 ; Johnson2015 ; Yamaji2014 ; Winter2016 , our results will be helpful for understanding of unusual behaviors in the real compounds. We emphasize that our numerical data obtained by the Majorana CDMFT+CTQMC method are quantitatively reliable in the calculated paramagnetic regime, as the cluster approximation and the analytic continuation are both well controlled. Although there are residual interactions in addition to the Kitaev-type ones in real materials, our results provide good references in the limit of the pure Kitaev model for interpreting the role of the additional interactions.

While we have calculated dynamical quantities of the Kitaev model in the wide $T$ range, the calculations were limited above $T_{\rm{L}}$ due to the phase transition which is artifact of the mean-field nature of CDMFT. It is necessary to develop more sophisticated method to study the dynamical properties below $T_{\rm L}$. The low-$T$ spin dynamics will be interesting, in particular, for extensions of the Kitaev model to three-dimensional lattices, such as hyperhoneycomb and hyperoctagon lattices Hermanns2014 . In the three-dimensional cases, in general, the Kitaev models may cause a finite-$T$ phase transition between the paramagnetic and QSL phases. Indeed, such an exotic transition was found for the hyperhoneycomb Kitaev model Nasu2014 . The phase transition is triggered by the confinement and deconfinement of emergent loops composed of excited fluxes Nasu2014 . This is a topological phase transition that cannot be described by a local order parameter. Although it is expected that dynamical quantities exhibit peculiar behavior associated with the topological phase transition, the CDMFT is not able to describe such a transition. Thus, with bearing the fact in mind that there are some candidates for the three-dimensional Kitaev QSLs Okamoto2007 ; Kuriyama2010 ; Modic2014 ; Takayama2015 ; Kimchi2014 the calculation of dynamical quantities in all $T$ range beyond the CDMFT is an interesting challenge left for future works.

Appendix A Cluster size dependence

In the CDMFT, we replace the lattice model to the impurity model with a finite-size cluster. The CDMFT becomes exact in the limit of infinite size cluster. Although the cluster size dependence was examined for the isotropic case with $\alpha=1.0$ in Supplemental Material for the previous study Yoshitake2016 , here we present the cluster size dependences of $\chi^{p}$ and $1/T_{1}^{p}$ for $\alpha=0.8$ and $1.2$ in comparison with the $\alpha=1.0$ case. As the onsite and NN-site components of $1/T_{1}^{p}$ shows almost the same $T$ dependences below $T_{\rm{H}}$ (see Fig. 9), we present only the onsite one.

Figure 14 shows the cluster size dependence of $\chi^{p}$ and $1/T_{1}^{p}$ obtained by the CDMFT+CTQMC calculations for three different types of clusters shown in Figs. 1(a), 13(a), and 13(b). In each type, we change the cluster sizes in the width in the $xy$-chain direction while keeping that in the $z$-bond direction. This is because the width in the $xy$-chain direction is rather relevant compared to that in the $z$-bond direction in the present CDMFT, presumably due to the Majorana representation based on the Jordan-Wigner transformation along the $xy$ chains. Hereafter, we define the size of the cluster by the average width in the $xy$-chain direction: for instance, $4.3$ for the cluster in Fig. 1(a), while $4$ and $5$ for Figs. 13(a) and 13(b), respectively.

As shown in Figs. 14(a)-14(j), the CDMFT+CTQMC results for $\chi^{p}$ show quick convergence with respect to the cluster width for all the cluster types. Even close to the artificial critical temperature $\tilde{T}_{\rm c}$, the results for the width larger than $4$ are almost convergent to the large width limit for all types of the clusters: the remnant relative errors are $\lesssim 5$%. Note that $\tilde{T}_{\rm c}\sim 0.014$ for $\alpha=1.0$, $\tilde{T}_{\rm c}\sim 0.0063$ for $\alpha=0.8$, and $\tilde{T}_{\rm c}\sim 0.013$ for $\alpha=1.2$ (for the rotated lattice coordinate used to calculate $\langle S^{p}_{j}(\tau)S^{p}_{j^{\prime}}\rangle$ for $p=x,y$, $\tilde{T}_{\rm c}$ becomes slightly lower: $\tilde{T}_{\rm c}\sim 0.0052$ for $\alpha=0.8$ and $\tilde{T}_{\rm c}\sim 0.0094$ for $\alpha=1.2$).

On the other hand, as shown in Figs. 14(k)-14(o), the cluster-size dependences of $1/T_{1}$ remains up to relatively higher $T$ than $\chi^{p}$. But the remnant relative errors are $\lesssim 10$% for the cluster width larger than $4$, which are sufficiently small to observe the characteristic $T$ dependences of $1/T_{1}$ as shown in Figs. 9.

Appendix B Accuracy of the maximum entropy method

In the CDMFT+CTQMC calculations, we calculate $S^{p}_{j,j^{\prime}}(\omega)$ from $\langle S^{p}_{j}(\tau)S^{p}_{j^{\prime}}\rangle$ by the MEM as described in Sec. II.4. In this Appendix, we examine the accuracy of the MEM in the limit of decoupled one-dimensional chains, i.e., $\alpha=1.5$ ($J_{z}=0$), where $S^{p}_{j,j^{\prime}}(\omega)$ can be obtained directly without the MEM. We also examine the accuracy by comparing $S^{p}_{j,j^{\prime}}(\omega)$ at sufficiently low-$T$ with the analytical solution in the ground state.

First, we show the comparison in the limit of decoupled one-dimensional chains, i.e., $\alpha=1.5$ ($J_{z}=0$). In this limit, the Kitaev Hamiltonian in Eq. (1) is written only by itinerant matter fermions $\{c\}$, in the form of Eq. (2) with $J_{z}=0$. In this noninteracting problem, following Ref. Derzhko2000 , we can calculate $S^{x}_{j,j^{\prime}}(\omega)$ by considering the real-time evolution (RTE) of $\langle S^{x}_{j}(t)S^{x}_{j^{\prime}}\rangle$, instead of the imaginary-time correlation $\langle S^{x}_{j}(\tau)S^{x}_{j^{\prime}}\rangle$, as

[TABLE]

We call this method as the RTE in the following. In the RTE calculations, we consider an $xy$ chain with 600 sites under the open boundary condition and take a sufficiently small $\epsilon=0.04$ in Eq. (48).

On the other hand, $S^{z}_{j,j^{\prime}}(\omega)$ has a nonzero value only for the onsite component, which is given by $4\langle S^{z}_{j}(\tau)S^{z}_{j}\rangle=\langle c_{j}(\tau)c_{j}\rangle$. Hence, $S^{z}_{j,j}(\omega)$ is obtained as

[TABLE]

where $D(\omega)$ is the DOS for itinerant matter fermions in the one-dimensional limit:

[TABLE]

We call this method to estimate $S^{z}_{j,j}(\omega)$ the exact-DOS in the following.

Figure 15 shows the results of $S^{p}_{j,j^{\prime}}(\omega)$ obtained by the MEM, RTE, and exact-DOS methods for the FM case with $\alpha=1.5$ ($J_{x}=J_{y}=1.5$ and $J_{z}=0$). We present both onsite and NN-site components for $S^{x}_{j,j^{\prime}}(\omega)$, while only the onsite one for $S^{z}_{j,j^{\prime}}(\omega)$. We find that overall $\omega$ dependence of $S^{p}_{j,j^{\prime}}(\omega)$ is well reproduced by the MEM. In particular, the agreement is excellent in the low $\omega$ region; the growth of $S^{x}_{j,j^{\prime}}(\omega=0)$ on decreasing $T$, which contributes to $1/T_{1}$, is well reproduced by the MEM. On the other hand, the relatively sharp structures at $\omega\sim 1.5$ are blurred in the MEM results for both $p=x$ and $z$, presumably because $\langle S^{p}_{j}(\tau)S^{p}_{j^{\prime}}\rangle$ is more insensitive to $S^{p}_{j,j^{\prime}}(\omega)$ in the larger $\omega$ region. Nevertheless, as shown in Figs. 15(e) and 15(f), the MEM results reproduce the broad incoherent peak of $S^{x}_{j,j}(\omega)-S^{x}_{\rm{NN}}(\omega)$.

Next, we examine the accuracy of the MEM for the data at sufficiently low $T$ with the analytical solution in the ground state Knolle2014 . Figure 16 shows $S^{z}(\Gamma,\omega)$ obtained by the Majorana CDMFT+CTQMC method for the FM case with $\alpha=0.8$ at $T=0.003$. In the ground state, the energy required to flip a single $\eta_{r}$ is $\Delta\simeq 0.042$ at $\alpha=0.8$. Reflecting the flux gap, $S^{z}(\Gamma,\omega)$ at low $T$ has a $\delta$-function like peak at $\Delta\simeq 0.042$ Knolle2014 . As shown in Fig. 16, our CDMFT+CTQMC result shows a peak at this energy, which is considered to precisely reproduce the low-energy structure of the dynamical spin structure factor.

From these observations, we consider that the MEM results for $S({\bf q},\omega)$ and $1/T_{1}$ in Sec. III.2 and III.3 are accurate enough to discuss the $T$ and $\omega$ dependences.

Appendix C Spin correlations as functions of $T$ and $\omega$

In this Appendix, we present the spin correlations as functions of $T$ and $\omega$, which are obtained by the MEM. Figure 17 shows the results for onsite and NN-site components for $\alpha=1.0$, $0.8$, and $1.2$. The data are used to obtain the dynamical quantities in Sec. III.2 and III.3.

Appendix D $T$ dependence of the Korringa ratio

Figures 18 and 19 display the $T$ dependences of the Korringa ratio defined as

[TABLE]

which is computed by using the NMR relaxation rate $1/T_{1}^{p}$ and the magnetic susceptibility $\chi^{p}$ obtained in Sec. III.3 and III.4. Interestingly, as shown in Fig. 18(a), $\mathcal{K}^{p}$ for the isotropic FM case is almost constant close to $1$ for $T_{\rm L}\lesssim T\lesssim T_{\rm H}$, which is apparently consistent with the behavior expected for free electron systems. This is also the case for the $x$ component for the FM case with $\alpha=1.2$, as shown in Fig. 18(c). However, the suggestive behavior is presumably superficial, as the results for the AFM cases as well as for $\alpha=0.8$ behave differently with substantial $T$ dependence.

Acknowledgements.

The authors thank M. Imada, Y. Kamiya, K. Ohgushi, S. Takagi, M. Udagawa, and Y. Yamaji for fruitful discussions. Y. M. thanks A. Banerjee, C. D. Batista, K.-Y. Choi, S. Ji, S. Naglar, and J.-H. Park for constructive suggestions. This research was supported by Grants-in-Aid for Scientific Research under Grants No. JP15K13533, No. JP16K17747, and No. JP16H02206. Parts of the numerical calculations were performed in the supercomputing systems in ISSP, the University of Tokyo.