Dynamics of the Kitaev-Heisenberg Model

TL;DR

This paper presents a matrix-product state method to compute dynamical responses in the Kitaev-Heisenberg model, revealing high-energy features related to spin liquids that align with experimental observations.

Contribution

Introduces a novel matrix-product state approach for dynamical functions in 2D Hamiltonians, applied to the Kitaev-Heisenberg model, uncovering high-energy spin liquid signatures.

Findings

High energy features beyond spin-wave theory in ordered phases

Observation of experimental-like features in zig-zag ordered phase

Proximity of spin liquid features at high energies

Abstract

We introduce a matrix-product state based method to efficiently obtain dynamical response functions for two-dimensional microscopic Hamiltonians, which we apply to different phases of the Kitaev-Heisenberg model. We find significant broad high energy features beyond spin-wave theory even in the ordered phases proximate to spin liquids. This includes the phase with zig-zag order of the type observed in $\alpha$-RuCl$_3$, where we find high energy features like those seen in inelastic neutron scattering experiments. Our results provide an example of a natural path for proximate spin liquid features to arise at high energies above a conventionally ordered state, as the diffuse remnants of spin-wave bands intersect to yield a broad peak at the Brillouin zone center.

Figures

Tables

Table 1. Table 1: Transition points α / π 𝛼 𝜋 \alpha/\pi for different circumferences sectors compared to exact diagonalization (ED) Gotfryd  et al. ( 2016 ) and infinite Projected Entangled Pair States (iPEPS) Osorio Iregui  et al. ( 2014 ) .

	ED	iPEPS	DMRG
			$L_2=6$		$L_2=12$
			gapped	gapless	gapped	gapless
AF/KSL	0.488	0.487	0.484	0.494	0.485	0.487
KSL/ZZ	0.510	0.513	0.523	0.513	0.514	0.512
FM/KSL	1.399	1.432	1.405	1.44	1.421	1.428
KSL/ST	1.577	1.557	1.573	1.548	1.562	1.558

Equations

$$H=\sum _{⟨i,j{⟩}_{\gamma }}{K}_{\gamma }{S}_i^{\gamma }{S}_j^{\gamma }+J\sum _{⟨i,j⟩}{\mathbf{S}}_i\cdot {\mathbf{S}}_j.$$

$${\mathcal{S}}^{\gamma \gamma }(\mathbf{k},\omega )=\frac{1}{2\pi }\sum _{\mathbf{r}}{\int }_{-\infty }^{\infty }{e}^{i(\omega t-\mathbf{k}\cdot \mathbf{r})}{C}^{\gamma \gamma }(\mathbf{r},t)\mathrm{d}t$$

$$W_l=\prod _{j\in loop}{\sigma }_j^{{\gamma }_j},$$

$${\tilde{W}}_l=\prod _{(j,k)\in loop}{u}_{jk}.$$

Full text

Introduction. The interplay of strong interactions and quantum fluctuations in spin systems can give rise to new and exciting physics. A prominent example are quantum spin liquids (QSL), as fascinating as they are hard to detect: they lack local order parameters and are instead characterized in terms of emergent gauge fields. On the experimental side, spectroscopic measurements provide particularly useful insights into such systems, in particular by probing the fractionalised excitations (e.g. deconfined spinons) accompanying the gauge field. Such measurements can be related to dynamical response functions, e.g. inelastic neutron scattering to the dynamical structure factor. On the theoretical side, determining the ground state properties of such quantum spin models is already a hard problem, and it is even more challenging to understand the dynamics of local excitations.

Here we present a combination of the density-matrix renormalization (DMRG) ground state method and a matrix-product states (MPS) based dynamical algorithm to obtain the response functions for generic two-dimensional spin systems. With this we are able to access the dynamics of exotic phases that can occur in frustrated systems. Moreover it is also very useful for regular ordered phases where one would conventionally use large-$S$ approximations, which in some cases cannot qualitatively explain certain high energy features Dalla Piazza et al. (2015); Banerjee et al. (2016a).

We demonstrate our method by applying it to the currently much-studied Kitaev-Heisenberg model (KHM) model on the honeycomb lattice

[TABLE]

The first term is the pure Kitaev model exhibiting strongly anisotropic spin exchange coupling Kitaev (2006). Neighboring spins couple depending on the direction of their bond $\gamma$ with $S^{x}S^{x}$, $S^{y}S^{y}$ or $S^{z}S^{z}$ (Fig. 1). The second is the $SU(2)$-symmetric Heisenberg term. The KHM serves as a putative minimal model for several materials including Na2IrO3, Li2IrO3 Chaloupka et al. (2013), and $\alpha$-RuCl3 Plumb et al. (2014). The pure model is an exactly solvable spin-$1/2$ model stabilizing two different Kitaev quantum spin liquids (KSL): a gapped $\mathbb{Z}_{2}$ one with abelian excitations (“A phase”) and one hosting gapless Majorana and gapped flux excitations (“B phase”) Kitaev (2006). If not stated otherwise, we use the parametrization $J=\cos\alpha$ and $K_{\gamma}=K=2\sin\alpha$. If $J=0$ and $K_{\gamma}$ bond-independent, the Kitaev model is in the B phase, which is stable under time-reversal symmetric perturbations as pointed out by Kitaev. Numerical studies of the ground state phase diagram of the KHM have shown an extended QSL phase for small $J$ and four symmetry broken phases for larger $J$ Chaloupka et al. (2013).

The dynamical response functions of the pure Kitaev model are known exactly and reveal characteristic features Knolle et al. (2014a, b), such as a spectral gap due to a spin flip not only creating gapless Majorana but also gapped flux excitations. This feature is perturbatively stable to small $J$ Song et al. (2016), but the influence of $J$ on high-energy features (or non-perturbatively at low energies) is unclear and of ongoing interestGotfryd et al. (2016). More pressingly, there appear to be proximate spin liquids Suzuki et al. (2015); Banerjee et al. (2016b), such as possibly the currently much-studied $\alpha$-RuCl3 Plumb et al. (2014); Johnson et al. (2015); Kim and Kee (2016); Sandilands et al. (2016a, b); Chaloupka and Khaliullin (2016); Lang et al. (2016); Halász et al. (2016); Koitzsch et al. (2016); Banerjee et al. (2016b, a), whose low-energy physics is consistent with spin waves on an ordered background, but whose broad high-energy features resemble those of a KSL. In particular, for intermediate energy scales there are star-like features Banerjee et al. (2016a) apparently arising from a combination of spin wave and QSL physics.

In this article, we first revisit the ground state phase diagram and confirm the previously found phases. The infinite cylinder geometry allows us to numerically confirm that the gaplessness of the KSL is robust throughout the entire phase. Secondly we use a recently introduced MPS based time evolution algorithm Zaletel et al. (2015) to obtain the dynamical spin structure factor. We benchmark our method by comparing to exact results for the Kitaev model and find a good agreement. We calculate the spectra of different (non-soluble) phases of the KHM. Most notably, we identify broad high energy continua even in ordered phases that are reminiscent of the broad features observed in recent experiments on $\alpha$-RuCl3 and which are moreover similar to the high energy features in the spin liquid phase, thus providing a concrete realisation of the concept of a proximate spin liquid.

Ground state phase diagram. We use the iDMRG algorithm on the KHM on infinite cylinders to map out the phase diagram. We choose cylinder geometries such that the corresponding momentum cuts contain the gapless Majorana modes of the Kitaev spin liquid. For the pure isotropic Kitaev model, there are gapless Majorana cones on the corners of the first Brilluoin zone, Fig. 1b. The full KHM has a $C_{6}$ symmetry which means that in the 2D limit these cones cannot shift. The iDMRG method determines the ground state of systems of size $L_{1}\times L_{2}$ where $L_{1}$ is in the thermodynamic limit and $L_{2}$ a finite circumference of up to $12$ sites beyond what is achievable in exact diagonalization. While traditionally iDMRG is used for finding the ground state of one-dimensional systems, it has become a fairly unbiased method for studying two-dimensional frustrated systems.

The resulting phase diagram for $L_{2}=12$ is shown in Fig. 2 (for the iDMRG simulations we keep $\chi=1200$ states), which agrees with previous studies Chaloupka et al. (2010); Jiang et al. (2011); Chaloupka et al. (2013); Sela et al. (2014); Osorio Iregui et al. (2014); Shinjo et al. (2015). For this $L_{2}$, the system is compatible with the sublattice transformation that maps zigzag to AF and stripy to FMChaloupka et al. (2010). Plotted are the ground state energy and the entanglement or von-Neumann entropy $S=-\mathrm{Tr}\rho^{\mathrm{red}}\log\rho^{\mathrm{red}}$ of the reduced density matrix $\rho^{\mathrm{red}}$ for a bipartitioning of the cylinder by cutting along a ring. Both the cusps in the energy density and the discontinuities of the entanglement entropy indicate first order transitions. A careful finite size scaling is difficult because of the large bond dimension needed and thus it is not possible to make definite statements about whether the transitions remain first order in the limit $L_{2}\rightarrow\infty$. The symmetry broken phases can be identified by measuring the local magnetization. We identify a Néel phase ($-0.185<\alpha/\pi<0.487$) that extends around the pure anti-ferromagnetic Heisenberg 111Note that due to the mapping of the 2D lattice onto a 1D chain, Mermin-Wagner-ColemanColeman (1973) applies at the pure AF-Heisenberg point and suppresses long range Néel order. In fact it is replaced by a gapped symmetry-preserving state, which extends over a finite region in parameter space. For further discussion, see appendix A.

point, the corresponding zigzag phase ($0.513<\alpha/\pi<0.894$), a ferromagnetic phase around the pure FM Heisenberg point ($0.894<\alpha/\pi<1.427$), and its stripy phase ($1.559<\alpha/\pi<1.815$). The two KSLs between Néel and zigzag as well as between FM and stripy are confirmed to be gapless. In particular, if $L_{2}$ is a multiple of six we use the finite entanglement scaling approach Calabrese and Cardy (2004); Tagliacozzo et al. (2008); Pollmann et al. (2009) and extract the expected chiral central charge $c=1$ for both KSLs, each of the two Majorana cones contributing $c=1/2$. See also appendix B. Note that when a gapless spin liquid is placed on a cylinder, the gauge field generically adjusts to open a gap He et al. (2016). In order to see gapless behaviour, we have to initiate the iDMRG simulations in the gapless sector to access a metastable state (see appendix C for additional details). The gapped ground state having a non-zero flux through the cylinder overestimates the stability of the QSL phases. It is notable how well the phase boundaries agree with those from the infinite projected entangled pair state (iPEPS) simulations Osorio Iregui et al. (2014).

Dynamical structure factor $\mathcal{S}(\boldsymbol{k},\omega)$. Starting from a ground state obtained using iDMRG, we calculate $\mathcal{S}(\boldsymbol{k},\omega)$ by Fourier transforming the dynamical correlation function $C^{\gamma\gamma}(\mathbf{r},t)=\langle S^{\gamma}_{\mathbf{r}}(t)S^{\gamma}_{\mathbf{0}}(0)\rangle.$ The real-time correlations can be efficiently obtained using a recently introduced matrix-product operator based time evolution method Zaletel et al. (2015). This allows for long range interactions resulting from unraveling the cylinder to a one-dimensional system which render standard methods like the time-evolving block decimation inefficient. Following the general strategy laid out in Refs. [Kjäll et al., 2011; Phien et al., 2012; Zauner et al., 2015], we perform the simulations for an infinite cylinder with a fixed circumference. Note that the entanglement growth and the resulting growth of the required number of states is generically slow as we only locally perturb the ground state and thus long times can be reached even in the cylinder geometry. We show results obtained for $0\leq t\leq T$ and to avoid Gibbs oscillations we multiply our real-time data with a Gaussian ($\sigma_{t}\approx 0.43T$). This corresponds to a broadening in $\omega$-space ($\sigma_{\omega}\approx\frac{2.3}{T}$). We use linear prediction to allow room for the tail of the Gaussian in real-time, but confirm that the final results do not depend on its details White and Affleck (2008). Thence,

[TABLE]

normalized as $\int\mathcal{S}^{\gamma\gamma}(\boldsymbol{k},\omega)\;\mathrm{d}\mathbf{k}\mathrm{d}\omega=\int\mathrm{d}\mathbf{k}$. If not stated otherwise, we present results for $\mathcal{S}(\boldsymbol{k},\omega)=\sum_{\gamma}\mathcal{S}^{\gamma\gamma}(\boldsymbol{k},\omega)$.

We benchmark the method by comparing our numerical approach to exact results for the pure Kitaev model. Figure 3a shows a comparison for the gapped Kitaev model in the A phase with $K_{x}/K_{y,z}=6$, the exact solution for $\mathcal{S}^{zz}(\boldsymbol{k}=0,\omega)$ shown in black. Our numerics (with resolution $\sigma_{\omega}\approx 0.06$ in units shown) for an infinite cylinder with $L_{2}=10$ (red) agrees well with such features as gap, bandwidth and total spectral weight. In the real-time data (inset), whilst the numerics agrees with the exact solution for the cylinder geometry, it overlaps with the 2D result only until a characteristic time scale corresponding to the perturbation traveling around the cylinder and then feeling the static fluxes inserted by the spin-flip. More generally we expect such timescales (after which 2D physics becomes 1D) to be particularly significant for systems with fractionalization. For Fig. 3b we take $K_{x}=K_{y}=K_{z}=-2$ being in the gapless KSL phase at $\alpha=\frac{3\pi}{2}$. Comparing the exact 2D result (black) to our numerics for a cylinder of circumference $L_{2}=6$ (red), we see qualitative similarities, such as a spectral gap (dashed lines; slightly obscured by our finite-time window), a dip where the fluxes suppress the van Hove singularity of the Majorana spectrumKnolle et al. (2014a), comparable bandwidth and strong low-energy weight. To better resolve the spectral gap, we rely slightly on linear prediction White and Affleck (2008) by using a real-time Gaussian envelope with $\sigma_{t}=0.56T$, corresponding to $\sigma_{\omega}\approx 0.045$. Two striking quantitative differences are (i) the spectral gap which for this circumference is approximately half that of the 2D limit, and (ii) the presence of a delta-peak on this gap ($\approx 4\%$ of total spectral weight). The latter, present for any cylinder, vanishes as $L_{2}\to\infty$. The inset compares exact real-time results on the cylinder Knolle (2016) with our numerics. Despite the true ground state on this cylinder being gapless and MPS only being able to capture gapped ground states exactly, we still find good agreement for appreciable times.

After this benchmarking, we explore $\mathcal{S}(\boldsymbol{k},\omega)$ in different phases of the KHM shown in Fig. 4, all with $\sigma_{\omega}\approx 0.06$. The pure Heisenberg FM ($\alpha=\pi$) can be solved in terms of linear spin wave theory (LSWT) and numerically captured with bond dimension $\chi=2$. Instead of this special point, in Fig. 4a we show results for $\alpha=1.1\pi$ (corresponding to $K=0.65J$) where we still find excellent agreement with LSWT. Note that there is an extremely small gap ($\approx 0.05|J|$) despite the presence of anisotropic couplings, as the entire KHM is $SU(2)$-symmetric in LSWT. We do not observe any strong cylinder effects on the dynamics, which is presumably related to the short correlation length and the absence of fractional excitations. The pure Heisenberg AFM (with small XXZ anisotropy) in Fig. 4b shows appreciable deviations from LSWT, with second order SWT Weihong et al. (1991) giving better agreement. Moreover, the weight in the spin waves is approximately halved, indicating the importance of higher order magnon contributions. Staying within the Néel phase but approaching the QSL, spin wave theory cannot even qualitatively describe Fig. 4c, with much weight in very broad high energy features unaccounted for.

Lastly we focus on a parameter regime producing zig-zag ordering like that found in $\alpha$-RuCl3 Johnson et al. (2015); Banerjee et al. (2016b, a). Fig. 5 shows $\mathcal{S}(\boldsymbol{k},\omega)$ for four different choices of $\alpha$: the first row contains the exact solution for the pure AFM Kitaev model, and the subsequent rows are all numerical results within the zigzag phase with increasing $\alpha$. For each $\alpha$ we show $\mathcal{S}(\boldsymbol{k},\omega)$ at fixed $\omega$: the columns display representative low-, mid- and high-energy features, with parameters $L_{2}=12$ and time cut-off $T=10$ corresponding to $\sigma_{\omega}\approx 0.23$. We average over the different symmetry broken directions. In appendix D, we show results for $L_{2}=6$ and $T=40$, revealing that even at this resolution the high-energy features stay very broad. The first column shows the low-energy physics of the Kitaev model being reconstructed into spin wave bands, with minima on the edges of the first Brillouin zone. For $\alpha=0.7\pi,0.8\pi$ these obey the $C_{6}$-symmetry, indicating that the cylinder geometry locally looks like 2D. Interestingly, the high-energy physics of the ordered phases is very similar to that of the pure Kitaev model: we have broad features centered around $\boldsymbol{k}=0$ which are diffuse w.r.t. $\omega$, with its characteristic energy and width simultaneously decreasing as $\alpha$ increases. The interplay between these low- and high- energy features then gives rise to different mid-energy shapes. In fact the six spin wave bands start on the edges of the first Brillouin zone. As the energy increases, these bands become increasingly diffuse, eventually overlapping in a very broad blob above the symmetric $\Gamma$ point $\boldsymbol{k}=0$. Both spin waves and blob sharpen as one moves away from the nearby QSL. Comparing with inelastic neutron data for $\alpha$-RuCl3 Banerjee et al. (2016a), we find the best qualitative agreement in Fig. 5 around $\alpha=0.7\pi$. In particular at intermediate energies there is a six-pointed-star whose arms point towards the edges of the first Brillouin zone. It is interesting to note that if we do not average over different symmetry broken directions, the low-energy physics strongly breaks the $C_{6}$ symmetry yet the six-pointed star at intermediate energies persists: thus even if we interpret these high energy features as the overlap of broad spin waves, at this point the effect of symmetry breaking has disappeared. Under what conditions such a symmetry restoration occurs more generally is an interesting question.

Conlusion. We have presented a new method for obtaining the dynamical properties of generic lattice spin models in (quasi-)two dimensions, which we expect to be useful for many future studies. In the KHM, our study reveals several features beyond spin-wave theory even in the ordered phases, providing a more detailed picture for the concept of a proximate spin liquid as potentially realised in $\alpha$-RuCl3.

Acknowledgements. We are grateful to Roser Valenti, Mike Zaletel and Johannes Knolle for stimulating discussions. In particular we thank Johannes for providing unpublished data for the dynamical correlations of the isotropic Kitaev model on the cylinder. This work was supported in part by DFG via SFB 1143 and Research Unit FOR 1807 through grants no. PO 1370/2-1.

Appendix A 1D vs 2D physics: symmetry breaking

From Monte-Carlo studies Reger et al. (1989) it is known that the ground state of the Heisenberg antiferromagnet (AFM) on the honeycomb lattice displays symmetry breaking Néel order. However, when we place the Heisenberg AFM on an infinitely long cylinder of finite circumference, it is in principle a 1D system and the presence of a continuous symmetry in fact forbids spontaneous symmetry breaking Coleman (1973). Instead we numerically find a gapped state which preserves both spin rotation and translation symmetry. This is analogous to the results for stacking an even number of coupled spin-$\frac{1}{2}$ Heisenberg chains White et al. (1994). The transition from 1D to 2D can be understood by noting that this symmetry-preserving state is effectively Néel-like within a correlation length $\xi$, the latter growing with circumference. Similarly to how one determines spontaneous symmetry breaking from finite size scaling in the context of exact diagonalization, one can conclude that the 2D limit achieves Néel order by scaling with respect to circumference.

The presence of a gap implies this symmetry-preserving state is stable under $SU(2)$-breaking perturbations. For example for $L_{2}=6$ it extends over $-0.2\pi\leq\alpha\leq 0.43\pi$, with a Néel order arising for larger $\alpha$ until we hit the spin liquid. The stability of this symmetry-preserving state under Kitaev perturbations is presumably related to the fact that the Néel order which arises in the 2D limit would have a very small spin gap. This is different for XXZ-type perturbations, which induce Néel order for relatively small anisotropies as shown in Fig. 6 (with $\Delta=1.1$), where our state is numerically converged (for large $\chi$) and the physics quickly becomes independent of circumference.

The DMRG simulations use a parameter $\chi$ which gives an upper bound on the entanglement. By limiting $\chi$ we can find a variational state with $\xi<L_{2}$. Locally this state then looks 2D and hence we can have symmetry breaking even for the $SU(2)$-symmetric Heisenberg model, as confirmed in Fig. 6. As we increase $\chi$, eventually $\xi$ becomes of the order of $L_{2}$, which signals the transition of 2D to 1D physics and the symmetry-preserving state arises. For $L=12$ the necessary $\xi$ is already out of reach, explaining the effective Néel order we see in Fig. 2. Similarly, in the zigzag phase there is an extended region with a gapped symmetry-restored ground state. This is in keeping with the sublattice transformation, which maps the zigzag to the Néel phase (in particular $\alpha=\frac{3}{4}\pi$ maps onto $\alpha=0$).

Appendix B Entanglement scaling of the gapless KSL

Matrix-product states (MPS) cannot capture algebraic ground state correlations. However, increasing the bond dimension gives an increasingly accurate estimate of the wave function. Calabrese and Cardy (2004) have shown that the entanglement entropy $S$ scales logarithmically with the correlation length $\xi$. In the MPS formalism, this is known as Finite-Entanglement Scaling with $S_{\chi}=c/6\log\xi_{\chi}$, where $\chi$ is the bond dimension of the MPS and $c$ is the chiral central charge Tagliacozzo et al. (2008); Pollmann et al. (2009).

Fig. 7 shows $S$ and $\log\xi$ for various MPS bond dimensions $\chi$ of up to $1024$. The lines serve as a guide to the eye corresponding to a slope with $c=1$. We observe a good match of the scaling for the pure Kitaev spin liquid at $\alpha=3/2\pi$. This reflects the fact, that the KSL can be mapped to a free fermion problem with two Majorana cones in the first Brillouin zone, each contributing $1/2$ to the central charge. The gapless nature persists within the whole KSL phase and the scaling suggests $c=1$.

Appendix C Ground sectors of the KSL on the cylinder

Similar to the plaquette operators $W_{p}=\prod_{j\in\hexagon}\sigma^{\gamma_{j}}_{j}$ we define a loop operator around the cylinder as

[TABLE]

where $\gamma_{i}=\{x,yz\}$ corresponds to the bond that is not part of the loop at site $i$. Following Kitaev (2006), $W_{l}$ can be expressed in terms of $\mathbb{Z}_{2}$ gauge field variables $u_{jk}$

[TABLE]

For our choice of lattice periodicity, both loop operators are related by a minus sign. Thus, $\tilde{W}_{l}\rightarrow+1$ (periodic boundary condition of the fermions) translates to $W_{l}\rightarrow-1$, which corresponds to the gapless sector if the cylinder is chosen such that cuts in reciprocal space go through the nodes of the Majorana cones. The second sector (antiperiodic boundary condition of the fermions) is always gapped and has a lower ground state energy than the gapless sector.

Regarding the computation of the ground state, we can now make use of the loop operator and initialize DMRG with a state $|\psi\rangle$ that has $\langle\psi|W_{l}|\psi\rangle=\pm 1$ depending on the desired sector. Table 1 contains the phase transitions for the gapped and the gapless sector and compares it to exact diagonalization (ED) and infinite Projected Entangled Pair States (iPEPS). As the gapped sector has a lower energy, its stability is enhanced and widens the KSL phase. This effect is more pronounced for a small circumference $L_{2}=6$.

Appendix D Dynamics of $L_{2}=6$ cylinder

In Fig. 8 we show $\mathcal{S}(\boldsymbol{k}=0,\omega)$ for the same choices of $\alpha$ as in Fig. 5, but now with a sharper $\omega$-resolution (corresponding to $T=40$) which is possible due to a smaller circumference ($L_{2}=6$). The finer features are most likely discretization effects due to the finite circumference, but the main points are that the broadness in $\omega$-space persists despite a finer resolution, and that the high-energy feature gets squeezed downward as we get further away from the nearby spin liquid. Note that the latter is a meaningful statement and not just due to an overall $\alpha$-dependent scaling of the Hamiltonian since the minima of the spin bands (as shown in the first column of Fig. 5) do not come down in energy (all at approximately $\omega=0.4$).