Temperature evolution of spin dynamics in two- and three-dimensional Kitaev models: Influence of fluctuating gauge fluxes

TL;DR

This study investigates how spin dynamics evolve with temperature in 2D and 3D Kitaev models, revealing distinct behaviors related to gauge flux fluctuations and phase transitions, using an advanced unbiased numerical method.

Contribution

Develops a novel QMC+CTQMC numerical approach to study low-temperature spin dynamics in Kitaev models, overcoming previous limitations and enabling analysis of gauge flux effects.

Findings

Smooth crossover in 2D honeycomb case

Singular behavior at phase transition in 3D hyperhoneycomb case

Low-temperature spin dynamics sensitive to gauge flux fluctuations

Abstract

The long-sought quantum spin liquid is a quantum-entangled magnetic state leading to the fractionalization of spin degrees of freedom. Quasiparticles emergent from the fractionalization affect not only the ground state properties but also thermodynamic behavior in a peculiar manner. We here investigate how the spin dynamics evolves from the high-temperature paramagnet to the quantum spin liquid ground state, for the Kitaev spin model describing the fractionalization into itinerant matter fermions and localized $Z_2$ gauge fluxes. Beyond the previous study [J. Yoshitake, J. Nasu, and Y. Motome, Phys. Rev. Lett. $\textbf{117}$, 157203 (2016)], in which the mean-field nature of the cluster dynamical mean-field theory prevented us from studying low-temperature properties, we develop a numerical technique by applying the continuous-time quantum Monte Carlo (CTQMC) method to statistical samples generated by the quantum Monte Carlo (QMC) method in a Majorana fermion representation. This QMC+CTQMC method is fully unbiased and enables us to investigate the low-temperature spin dynamics dominated by thermally excited gauge fluxes, including the unconventional phase transition caused by gauge flux loops in three dimensions, which was unreachable by the previous methods. We apply this technique to the Kitaev model in both two and three dimensions. Our results clearly distinguish two cases: while the dynamics changes smoothly through the crossover in the two-dimensional honeycomb case, it exhibits singular behaviors at the phase transition in the three-dimensional hyperhoneycomb case. We show that the low-temperature spin dynamics is a sensitive probe for thermally fluctuating gauge fluxes that behave very differently between two and three dimensions.

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^{\prime }}{c}_j-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 }},$$

$$\chi =\frac{1}{N}\sum _{j,j^{\prime }}{\int }_0^{\beta }d\tau ⟨S_j^z(\tau )S_{j^{\prime }}^z⟩,$$

$$1/T_1=S_{j,j}^x(\omega =0)+S_{j,j}^y(\omega =0),$$

$$1/T_1=S_{j,j^{\prime }}^x(\omega =0)+S_{j,j^{\prime \prime }}^y(\omega =0),$$

$$\displaystyle\Delta W_{p}=\frac{1}{N_{p}T^{2}}\Big{(}\big{\langle}\big{(}\sum_{p}W_{p}\big{)}^{2}\big{\rangle}-\big{\langle}\sum_{p}W_{p}\big{\rangle}^{2}\Big{)},$$

$$S(\text{q},\omega )=\frac{1}{N}\sum _{j,j^{\prime }}{e}^{i\text{q}\cdot ({\text{r}}_j-{\text{r}}_{j^{\prime }})}{S}_{j,j^{\prime }}^z(\omega ),$$

Full text

††preprint: APS/123-QED

I Introduction

The quantum spin liquid (QSL) is an exotic state of matter in insulating magnets showing no magnetic order down to zero temperature ($T$) Anderson1973 ; Balents2010 . It is not characterized by any conventional order parameter, but known to exhibit topological quantum entanglement resulting in fractionalization of the fundamental spin degrees of freedom Wen1991 ; Misguich2011 . This is purely quantum mechanical nature arising in strongly correlated many-body systems, as seen in fractional charges by the fractional quantum Hall effect Tsui1982 ; Stormer1999 . Although the spin fractionalization has attracted great attention for identifying the QSL in candidate materials, the unambiguous detection remains largely elusive Balents2010 ; Yamashita2008 ; Yamashita2009 ; Yamashita2010 .

The Kitaev spin model, originally introduced on a two-dimensional (2D) honeycomb lattice Kitaev2006 , has generated a new trend in the study of QSLs. This is because of the following virtues of this model. First of all, the model is exactly soluble in the ground state, and the exact ground state is a QSL. The exact solution is obtained by representing the spin operators by Majorana fermion operators, which simultaneously provides canonical formulation of the fractionalization: the elementary spin excitations are described by itinerant matter fermions and localized $Z_{2}$ gauge fluxes, both of which are composed of the Majorana fermions. Furthermore, the model can be extended to any tri-coordinate lattices with preserving the solubility, even in three dimensions (3D) Mandal2009 ; Hermanns2015 ; O'Brien2016 . Last but not least, the bond-dependent anisotropic interaction in this model has a realization in some magnetic materials with strong spin-orbit coupling Jackeli2009 . All these features have accelerated the combined studies between theory and experiment for realization and identification of Kitaev QSLs Nussinov2015 ; Trebst_preprint .

Among many consequences of the spin fractionalization unveiled by the recent studies of the Kitaev model is thermal fractionalization, i.e., thermodynamic signatures originating from different energy scales of the fractionalized quasiparticles Nasu2015 . The thermal fractionalization manifests itself in, for instance, two peaks in the specific heat at $T=T_{L}$ and $T_{H}$ ($T_{L}\ll T_{H}\sim J$, where $J$ is the dominant Kitaev coupling) and successive entropy release by a half of $\log 2$ around these temperatures. Besides, the spin dynamics is also of importance for experimental identification of the fractionalization. In the previous studies Yoshitake2016 ; Yoshitake_preprint , the authors calculated dynamical quantities for the 2D Kitaev model, developing the cluster extension of the dynamical mean-field theory (CDMFT) in a Majorana fermion representation and combining it with the continuous-time quantum Monte Carlo (CTQMC) method. The CDMFT+CTQMC study revealed an interesting aspect of the fractionalization: dichotomy between static and dynamical spin correlations. This was shown by the significant $T$ evolution of the magnetic susceptibility $\chi$, the NMR relaxation rate $1/T_{1}$, and the dynamical spin structure factor $S(\mathbf{q},\omega)$ in the $T$ regime below $T_{H}$ where the static spin correlations saturate and almost $T$ independent.

Despite the successful calculations of dynamical properties, the applicable $T$ range of the CDMFT+CTQMC method is limited: the method does not give reasonable results at very low $T$. This is due to the occurrence of phase transition at $T\sim T_{L}$ as an artifact of the mean-field approximation in the CDMFT. Moreover, the CDMFT+CTQMC method is not suitable for the Kitaev model on 3D lattices by the following reasons. One is that a larger cluster is necessary in the CDMFT, as the unit cell, or more strictly speaking, the smallest loop of lattice sites, for which the conserved $Z_{2}$ gauge flux is defined, becomes larger for 3D than 2D in general. Another reason is that the 3D extensions of the Kitaev model may cause a phase transition, which might be hard to capture by the CDMFT. For instance, the Kitaev model on a 3D hyperhoneycomb lattice exhibits an unconventional phase transition triggered by proliferation of loops composed of thermally excited gauge fluxes Nasu2014b ; Nasu2014 . The cluster approximation in the CDMFT is not suitable to describe such a topological transition characterized by global quantities beyond the cluster. An alternative method is desired to study the spin dynamics, including the low-$T$ behavior.

Besides such a theoretical demand, it is crucial to clarify the spin dynamics of the Kitaev model in the whole $T$ range also from the experimental point of view. Recently, many candidates have been explored in both quasi-2D and 3D materials Singh2010 ; Singh2012 ; Plumb2014 ; Takayama2015 ; Modic2014 . Some indications of the fractionalization were observed, for instance, in the specific heat Mehlawat2017 , magnetic Raman scattering Sandilands2015 ; Glamazda2016 , inelastic neutron scattering Banerjee2016 ; Banerjee_preprint ; Do_preprint , and thermal transport Hirobe_preprint ; Leahy2017 . However, such indications are for rather high-$T$ features, corresponding to the theoretical predictions around and below $T_{H}$ associated with itinerant matter fermions Nasu2015 ; Nasu2016 ; Yoshitake2016 ; Yoshitake_preprint ; Nasu_preprint . It is highly desired to experimentally capture another indications dominated by thermally excited gauge fluxes at lower $T$. Although all the candidate materials exhibit a magnetic order at low $T$, several efforts have been made for suppressing the order, e.g., by external pressure Takayama2015 ; Breznay_preprint , magnetic field Ruiz_preprint ; Hentrich_preprint ; Wolter_preprint , and chemical substitution Lampen-Kelley_preprint . Given such an upsurge of interest, it is highly important to clarify the dynamical behavior of the 2D and 3D Kitaev models down to the lowest $T$.

In this paper, we propose a new numerical method which overcomes the problems in the previous CDMFT+CTQMC method. We here adopt the quantum Monte Carlo (QMC) method, instead of the CDMFT, for generating statistical samples used in the CTQMC calculations. The QMC method is also formulated on the basis of a Majorana fermion representation, which has been used to compute thermodynamic properties in a series of previous studies for the Kitaev models on several tri-coordinate lattices Nasu2015 ; Nasu2014b ; Nasu2015b ; Nasu2016 ; Nasu_preprint . Thus, the new combined method, which we call the QMC+CTQMC method, provides a versatile technique, free from biased approximation. We demonstrate that the method is applicable in a wider $T$ range, including the low-$T$ region below $T\sim T_{L}$, which was not accessible by the previous CDMFT+CTQMC method. In the 2D honeycomb case, comparing the data of $\chi$ and $1/T_{1}$ by the CDMFT+CTQMC and QMC+CTQMC methods, we show that although the former works quite well above the fictitious critical temperature, only the latter can give reasonable results at lower $T$. In the 3D hyperhoneycomb case, we present the QMC+CTQMC results for $\chi$, $1/T_{1}$, and $S(\mathbf{q},\omega)$. From the comparison between the 2D and 3D results, we clarify the signatures arising from the difference of the system dimension. While everything changes smoothly through the crossover at $T=T_{L}$ in the 2D honeycomb case, the dynamical quantities exhibit singular behaviors in the 3D hyperhoneycomb case at the phase transition caused by the topological nature of excited gauge flux loops. Thus, the QMC+CTQMC is applicable to the unconventional phase transition in 3D, which is not accessible by the CDMFT+CTQMC method. Our results show that the dynamical properties at low $T$ depend substantially on the system dimension, despite almost dimension-independent behavior of the static spin correlations. This is the low-$T$ aspect of the dichotomy between static and dynamical spin correlations, which was found in the intermediate $T$ region in the previous study.

The structure of this paper is as follows. In Sec. II, we introduce the Kitaev model and its Majorana fermion representation. We also present the details of the QMC+CTQMC method. In Sec. III, we present the QMC+CTQMC results for the 2D and 3D cases in Sec. III.1 and III.2, respectively. Finally, Sec. IV is devoted to the summary.

II Model and method

In this study, we consider the Kitaev model on a 2D honeycomb lattice [Fig. 1(a)] and 3D hyperhoneycomb lattice [Fig. 1(b)], whose Hamiltonian is given in the common form Kitaev2006 ; Mandal2009

[TABLE]

Here, $p=x,y,z$ represents one of the three different types of bonds on the tri-coordinate lattices, and $\langle j,j^{\prime}\rangle_{p}$ denotes a set of neighboring sites $j,j^{\prime}$ on the $p$ bonds; see Figs. 1(a) and 1(b). $S_{j}^{p}$ represents the $p$ component of quantum spin $S=1/2$ at site $j$, and $J_{p}$ is the coupling constant for the $p$ bond.

A mathematically faithful representation of the Hamiltonian in Eq. (1) is obtained by applying the Jordan-Wigner transformation along the chains composed of the $x$ and $y$ bonds Chen2007 ; Feng2007 ; Chen2008 :

[TABLE]

where $c_{j}$ and $\bar{c}_{j}$ are two types of Majorana fermion operators at site $j$; $\eta_{r}=i\bar{c}_{j}\bar{c}_{j^{\prime}}$ is defined on each $z$ bond connecting sites $j$ and $j^{\prime}$. The sum over $(j,j^{\prime})_{p}$ is taken for the neighboring sites $j$ and $j^{\prime}$ colored by black and white, respectively, in Figs. 1(a) and 1(b). The bond variable $\eta_{r}$ commutes with the Hamiltonian as well as other $\eta_{r^{\prime}}$, and $\eta_{r}^{2}=1$; hence, $\eta_{r}$ is a conserved $Z_{2}$ variable taking $\pm 1$. The ground state is exactly obtained as the state with all $\eta_{r}=+1$ for both honeycomb and hyperhoneycomb cases. The exact ground state is shown to be a QSL, both gapless and gapped depending on the ratios between the coupling constants $J_{p}$ Kitaev2006 . The elementary excitations are also exactly described by the operators $\{c_{j}\}$ and $\{\eta_{r}\}$. In this Majorana fermion representation, therefore, the original spin operators $\{\mathbf{S}_{j}\}$ are fractionalized into $\{c_{j}\}$, which describe itinerant Majorana fermions called matter fermions, and $\{\eta_{r}\}$, which are the localized $Z_{2}$ variables.

The $Z_{2}$ variables $\{\eta_{r}\}$ are related with the $Z_{2}$ gauge fluxes $W_{p}$ discussed in the original paper by Kitaev Kitaev2006 . The gauge flux is also a conserved $Z_{2}$ quantity defined for each elementary plaquette [a hexagon in the 2D honeycomb case and a ten-site plaquette in the 3D hyperhoneycomb case; see Figs. 1(a) and 1(b)]: it is defined by the product of $\eta_{r}$ belonging to the plaquette $p$, as $W_{p}=\prod_{r\in p}\eta_{r}$. The ground state with all $\eta_{r}=+1$ corresponds to the state with all $W_{p}=+1$, which is called the flux-free state. At nonzero $T$, the gauge fluxes are thermally excited from the flux-free state by flipping $W_{p}$.

In the previous study, the authors have developed the CDMFT+CTQMC method for calculating the finite-$T$ spin dynamics of the Kitaev model in Eq. (1), by using the Majorana representation in Eq. (2) Yoshitake2016 ; Yoshitake_preprint . In this method, we generate the configurations of the $Z_{2}$ variables $\{\eta_{r}\}$ by the CDMFT, and compute the imaginary-time spin correlations by applying the CTQMC calculations to each configuration. The combined method successfully delivers precise data for the dynamical properties in a wide $T$ range. A problem in the CDMFT+CTQMC method is that the cluster approximation in the CDMFT part leads to a fictitious phase transition at $T=\tilde{T}_{c}$ by ordering of $\{\eta_{r}\}$. In the 2D Kitaev model on the honeycomb lattice, there is no phase transition at a nonzero $T$ and only two crossovers occur at very different $T$ scales, $T_{L}$ and $T_{H}$ Nasu2015 . In the isotropic case with $J_{x}=J_{y}=J_{z}=J=\pm 1$, $T_{L}\simeq 0.012$ and $T_{H}\simeq 0.375$; the fictitious $\tilde{T}_{c}\simeq 0.014$ is slightly higher than $T_{L}$. Meanwhile, in the 3D case on the hyperhoneycomb lattice, the model exhibits a phase transition at $T=T_{c}$ ($T_{c}\simeq 0.0039$ for the isotropic case), but it is not due to the ordering of $\{\eta_{r}\}$: the transition is caused by global objects, i.e., closed loops composed of thermally excited gauge fluxes $W_{p}$ Nasu2014b . Thus, the phase transition by ordering of $\{\eta_{r}\}$ in the CDMFT is an artifact arising from the mean-field nature. Because of this problem, the CDMFT+CTQMC method is not applicable to the very low-$T$ region around and below $T_{L}$ in 2D and $T_{c}$ in 3D 111At sufficient low $T$, where almost all $\eta_{r}=+1$, the CDMFT+CTQMC method reproduces well the quantum spin liquid nature..

In order to solve this problem, instead of the CDMFT, we here adopt the real-space QMC simulation, which has been used to calculate static quantities in the previous studies Nasu2015 ; Nasu2014b ; Nasu2015b ; Nasu2016 ; Nasu_preprint . Using the QMC simulation, we generate statistical samples of the configuration of localized $Z_{2}$ variables $\{\eta_{r}\}$, for which the dynamical spin correlations are computed by the CTQMC simulation. In this case, we can study much larger system sizes than the clusters used in the CDMFT, which enables us to systematically investigate the low-$T$ dynamical properties including the unconventional phase transition in 3D without biased approximation. We call this new combined technique the QMC+CTQMC method.

In Sec. III, we compute the dynamical properties for the isotropic case with $J_{x}=J_{y}=J_{z}=J=\pm 1$ by the QMC+CTQMC method; $J=+1$ corresponds to the ferromagnetic (FM) case, while $J=-1$ the antiferromagnetic (AFM) case. All the static quantities, such as the specific heat, behave in the same manner for the FM and AFM cases, and hence, the crossover and phase transition temperatures are common to the two cases. The configurations of $\{\eta_{r}\}$ are generated by the QMC calculations under the same conditions with the previous studies Nasu2014b ; Nasu2015 . Note that the QMC simulation is done for finite-size clusters with the open boundary condition, at least, in one direction. For each configuration, we perform the CTQMC calculations for the $z$ bonds, sufficiently far from the open boundaries. Typically, we select $60$ ($16$-$18$) bonds in the 2D (3D) case near the central region of each cluster (away from the open boundary), and average the results over the bonds. In each CTQMC calculation, we typically perform $5\times 10^{3}$ measurements at every $20$ MC steps, after $10^{5}$ MC steps for initial relaxation. To obtain the dynamical quantities as functions of the real frequency from the imaginary-time spin correlations, we perform the maximum entropy method (MEM) under the same conditions with the previous CDMFT+CTQMC study Yoshitake_preprint ; we use the Legendre polynomial up to $100$th order for $T\leq 0.006$, while we expand up to $50$th order for higher $T$ as well as for the 2D case.

III Results

III.1 2D honeycomb

First, we show the results for the 2D case on the honeycomb lattice. Figure 2 displays the QMC+CTQMC results for the magnetic susceptibility $\chi$ and the NMR relaxation rate $1/T_{1}$. $\chi$ is calculated from the imaginary-time spin correlations, without using the MEM, as

[TABLE]

where $N$ is the system size and $\beta=1/T$ is the inverse temperature (we set the Boltzmann constant $k_{B}=1$ and the reduced Planck constant $\hbar=1$). On the other hand, we compute $1/T_{1}$ by Yoshitake_preprint

[TABLE]

for the onsite component and

[TABLE]

for the nearest-neighbor(NN)-site component separately, where $S^{p}_{j,j^{\prime}}(\omega)$ is the spin correlations as a function of the real frequency $\omega$ obtained by the MEM from $\langle S^{p}_{j}(\tau)S^{p}_{j^{\prime}}\rangle$. Here, $j^{\prime}$ and $j^{\prime\prime}$ are the sites neighboring to site $j$ on the $x$ and $y$ bonds, respectively. Note that both $\chi$ and $1/T_{1}$ are isotropic in spin space for the current isotropic case with $J_{x}=J_{y}=J_{z}$ on the honeycomb lattice.

As shown in Figs. 2(a) and 2(b), the results for different system sizes $L=12$ and $20$ agree with each other ($N=2L^{2}$), indicating that the QMC+CTQMC results well converge with respect to the system size. In the figures, the previous CDMFT+CTQMC results are also plotted by gray symbols for comparison Yoshitake2016 . In the CDMFT+CTQMC method, as mentioned above, the cluster mean-field approximation leads to a fictitious phase transition at $\tilde{T}_{c}\simeq 0.014$, and hence, we plot the data above $\tilde{T}_{c}$. We find that the QMC+CTQMC results well agree with the CDMFT+CTQMC ones for $T\gtrsim\tilde{T}_{c}$, which supports the validity of the latter for $T\gtrsim\tilde{T}_{c}$. While such validity was claimed for the static quantities in the previous studies Yoshitake2016 ; Yoshitake_preprint , the present results demonstrate it explicitly for the dynamical quantities.

The present QMC+CTQMC method enables us to study the low-$T$ region around and below the low-$T$ crossover temperature $T_{L}\simeq 0.012$, beyond $\tilde{T}_{c}$ in the CDMFT+CTQMC result. $T_{L}$ is the temperature where the localized $Z_{2}$ gauge fluxes $W_{p}$ begin to be frozen into the flux-free state while decreasing $T$ Nasu2015 . Thus, our results show how the dynamical properties are affected by thermally excited gauge fluxes. Figure 2(a) indicates that, while decreasing $T$ around $T_{L}$, $\chi$ decreases slightly and changes the curvature from upward to downward convex, for both the FM and AFM cases. While further decreasing $T$, $\chi$ appears to converge to a nonzero value, as expected for the system which does not conserve the $z$ component of total spin. The asymptotic value is almost one order of magnitude larger for the FM case than the AFM case. On the other hand, as shown in Fig. 2(b), $1/T_{1}$ decreases below the peak slightly above $T_{L}$ as partly seen in the CDMFT+CTQMC results Yoshitake2016 , and continues to decrease around $T_{L}$ reaching to almost zero below $T\sim 0.005$. The low-$T$ suppression is due to a nonzero flux gap required to excite the $Z_{2}$ gauge fluxes from the flux-free ground state Kitaev2006 .

We also compute the $T$ derivatives of $\chi$ and $1/T_{1}$, as shown in Figs. 3(a) and 3(b), respectively. Both derivatives show a peak around $T_{L}$, but change smoothly without showing any singularity. For comparison, we also compute the thermal fluctuation of gauge fluxes $W_{p}$ by the QMC method, defined by

[TABLE]

where $N_{p}$ is the number of plaquettes in the system. Note that $\Delta W_{p}$ corresponds to the specific heat in the anisotropic limit (toric code), where the effective Hamiltonian is given in the form $\mathcal{H}\propto\sum_{p}W_{p}$ Kitaev2006 ; hence, $\Delta W_{p}$ measures the energy fluctuation related to the gauge fluxes. As shown in Fig. 3(c), $\Delta W_{p}$ also shows a broad peak around $T_{L}$, similar to the $T$ derivatives of $\chi$ and $1/T_{1}$. All these smooth changes with broad peaks are consistent with the fact that $T_{L}$ is not a phase transition but just a crossover in the 2D case Nasu2015 . Furthermore, the similar behavior between three quantities in Fig. 3 suggests that the $T$ derivatives of $\chi$ and $1/T_{1}$ provide good probes for the fluctuations of gauge fluxes.

Interestingly, $d\chi/dT$ behaves differently between the FM and AFM cases, as shown in Fig. 3(a): it is negative for $T_{L}\lesssim T\lesssim T_{H}$ and changes the sign to positive just above $T_{L}$ for the FM case, while mostly positive in the same $T$ range for the AFM case. The qualitative difference will be useful for identifying the sign of the dominant Kitaev interactions in candidate materials. The details of the difference between the FM and AFM cases, including the nonlinear components of the magnetic susceptibility, will be reported elsewhere.

III.2 3D hyperhoneycomb

Next, we turn to the 3D case on the hyperhoneycomb lattice. Figure 4 shows the QMC+CTQMC results for $\chi$ and $1/T_{1}$. The system size is given by $N=4L^{3}$: $256$, $500$, and $864$ sites for $L=4$, $5$, and $6$, respectively. Note that in the hyperhoneycomb lattice the $z$ bond is not equivalent to the $x$ and $y$ bonds from the lattice symmetry; we compute $\chi$ by Eq. (3) and $1/T_{1}$ by Eqs. (4) and (5) with replacing $\langle S^{p}_{j}(\tau)S^{p}_{j^{\prime}}\rangle$ ($p=x,y$) by $\langle S^{z}_{j}(\tau)S^{z}_{j^{\prime}}\rangle$ for simplicity. The overall $T$ dependence is similar to the 2D results as follows. The high-$T$ behaviors above $T_{H}$ are almost unchanged from the 2D cases, presumably because the bandwidth of matter fermions is independent of the dimensionality. With a decrease of $T$, $\chi$ begins to deviate from the Curie-Weiss behavior below $T\sim T_{H}$ and converges to a nonzero value after showing a peak, while $1/T_{1}$ increases below $T_{H}$ and strongly suppressed due to the flux gap after showing a peak at a low $T$. Nonetheless, there are quantitative differences. For instance, the peak of $\chi$ for the FM case is more than twice larger that that for the 2D case. Simultaneously, the change at low $T$ is much steeper in 3D than 2D. Similar behaviors are also seen in $1/T_{1}$. We will briefly comment on the quantitative differences in the end of this section.

However, we also find a qualitative difference between 3D and 2D in the low-$T$ behavior. The 3D hyperhoneycomb model exhibits a phase transition at $T_{c}\simeq 0.0039$ Nasu2014b . The phase transition takes place between the high-$T$ paramagnet and the low-$T$ QSL, driven by the proliferation of loops composed of the localized $Z_{2}$ gauge fluxes $W_{p}$. Thus, the transition is of topological nature, not characterized by local spin operators contrary to conventional magnetic ordering Nasu2014b ; Nasu2014 . Nevertheless, we find singular behaviors in both $\chi$ and $1/T_{1}$, as more clearly seen in the $T$ derivatives shown in Figs. 5(a) and 5(b). Both $T$ derivatives show a sharp peak at $T\simeq T_{c}$, which becomes sharper for larger system sizes. We also plot the thermal fluctuation of gauge fluxes $\Delta W_{p}$ in Eq. (6) in Fig. 5(c). In this 3D case, $\Delta W_{p}$ shows a similar sharp peak to $d\chi/dT$ and $d(1/T_{1})/dT$. All these behaviors are in stark contrast to the 2D case, where the crossover at $T_{L}$ leads to smooth $T$ dependence as shown in Fig. 3.

The low-$T$ behaviors of the dynamical quantities are substantially different from those in 2D, not only in the critical behavior associated with the phase transition but also the larger $T$ dependence. This clear difference depending on the spatial dimension is rather surprising when considering that the static spin correlations are not much different between 2D and 3D in the whole $T$ range Nasu2015 ; Nasu2014b . In the previous CDMFT+CTQMC study, the authors unveiled a prominent feature of the Kitaev QSL, the dichotomy between static and dynamical spin correlations, from the $T$ dependence of the static spin correlations for NN sites and $1/T_{1}$ Yoshitake2016 ; Yoshitake_preprint . The significant dimensional dependence at low $T$ found here is another aspect of the dichotomy.

We note that the difference of the sign of $d\chi/dT$ between the FM and AFM cases for $T_{c}\lesssim T\lesssim T_{H}$ is also seen in the 3D case, as shown in Fig. 3(c). We also note that the behavior of $d\chi/dT$ is similar to that found in the effective model in the anisotropic limit $J_{z}\gg J_{x},J_{y}$ Nasu2014 .

Finally, we show the QMC+CTQMC results for the dynamical spin structure factor $S(\mathbf{q},\omega)$ for the 3D case in Fig. 6. $S(\mathbf{q},\omega)$ is defined as

[TABLE]

where $\mathbf{r}_{j}$ represents the position vector for site $j$. The results are plotted along the symmetric lines in the first Brillouin zone shown in Fig. 1(c). The overall $T$ and $\omega$ dependence is similar to the 2D case reported in the previous study Yoshitake2016 : almost $\mathbf{q}$-independent incoherent response around $\omega=0$ for $T\gtrsim T_{H}$, growth of the incoherent spectra around $\omega=|J|$ below $T_{H}$, and a rapid increase of the quasi-elastic response while approaching $T_{c}$. Also, as in the 2D case, the difference in the sign of $J$ appears in the $\mathbf{q}$ dependence of the spectral intensity. We note that the lowest-$T$ data below $T_{c}$ in Fig. 4(a) agree well with the previous $T=0$ result Smith2015 .

Figures 7(a) and 7(b) display the low-$\omega$ part of $S(\Gamma,\omega)=S(\mathbf{q}=0,\omega)$ around $T_{c}$ for the FM case. Qualitatively similar behaviors are also seen for $S(\mathbf{q},\omega)$ near the zone boundary for the AFM case. With a decrease of $T$ across $T_{c}$, the quasi-elastic peak near $\omega=0$ shifts to a slightly higher $\omega$, leading to the opening of the flux gap below $T_{c}$. The peak height is almost unchanged across $T_{c}$. For comparison, we plot the corresponding data for the 2D honeycomb case around $T_{L}$ in Figs. 7(c) and 7(d). In the 2D case, the peak above $T_{L}$ is much broader with a lower peak height compared to the 3D case. When lowering $T$ across $T_{L}$, the peak becomes sharper with a shift of the peak position to a higher $\omega$. These differences between 2D and 3D are closely related with the quantitatively different behaviors of $\chi$ and $1/T_{1}$ observed in Figs. 2 and 4 as follows. The sharper peak of $S(\Gamma,\omega)$ near $\omega=0$ already existing above $T_{c}$ in 3D corresponds to much larger values of $\chi$ and $1/T_{1}$ just above $T_{c}$ compared to the 2D results above $T_{L}$. Furthermore, the shift of the peak across $T_{c}$ in Fig. 7(b) is related with the steep changes of $\chi$ and $1/T_{1}$ around $T_{c}$.

IV Summary

We have developed the numerical method for studying the spin dynamics of the Kitaev models by combining the QMC and CTQMC methods on the basis of a Majorana fermion representation. The QMC+CTQMC method overcomes the shortcoming in the previous CDMFT+CTQMC method, and enables us to investigate the very low-$T$ region where the $Z_{2}$ gauge fluxes play a role. The experimental observation of the gauge fluxes is one of the open issues in the Kitaev-type QSLs, and hence, the theoretical results obtained by our method provide the references for the experiments in candidate materials.

We have applied the QMC+CTQMC method to the 2D and 3D Kitaev models. Calculating the magnetic susceptibility, the NMR relaxation rate, and the dynamical spin structure factor, we discussed the influences of thermally fluctuating gauge fluxes, with focusing on the differences arising from the spatial dimensions. In the 2D honeycomb case, everything changes smoothly while lowering $T$, reflecting the crossover associated with particlelike gauge flux excitations. In contrast, in the 3D hyperhoneycomb case, the system exhibits a phase transition by the proliferation of looplike gauge flux excitations, which leads to singular behaviors in the dynamical properties. We found that the dichotomy between static and dynamical spin correlations, which begins below the high-$T$ crossover associated with itinerant matter fermions, persists down to the low-$T$ region, in a more peculiar form reflecting thermally excited gauge fluxes; while the dichotomy in the higher-$T$ region is rather universal independent of the spatial dimension, the low-$T$ one appears differently between 2D and 3D, reflecting the different nature of the localized $Z_{2}$ gauge flux excitations. We showed that the $T$ derivatives of the magnetic susceptibility and the NMR relaxation rate provide good probes for fluctuating gauge fluxes in both 2D and 3D. Our results will be useful for identifying the contributions from the $Z_{2}$ gauge fluxes in the experimental candidates and also the nature of their excitations.

Acknowledgements.

This research was supported by Grants-in-Aid for Scientific Research under Grants No. JP15K13533, No. JP16K17747, No. JP16H02206, and No. JP16H00987. Parts of the numerical calculations were performed in the supercomputing systems in ISSP, the University of Tokyo.