# Valence bond fluctuations in the Kitaev spin model

###### Abstract

We introduce valence bond fluctuations, or bipartite fluctuations associated to bond-bond correlation functions, to characterize quantum spin liquids and their entanglement properties. Using analytical and numerical approaches, we find an identical scaling law between valence bond fluctuations and entanglement entropy in the two-dimensional Kitaev spin model and in one-dimensional chain analogues. We also show how these valence bond fluctuations can locate quantum phase transitions between the three gapped and the gapless Majorana semi-metal phases in the honeycomb model. We then study the effect of a uniform magnetic field along the direction opening a gap in the intermediate phase which becomes topological. We still obtain a robust signal to characterize the transitions towards the three gapped phases. The linear scaling behavior of such bipartite fluctuations in two dimensions is also distinguishable from the one in the Néel magnetic state.

Introduction.— The quest of quantum spin liquids in the Mott regime anderson1987; balents2010 has been a great challenge these last decades in relation with the discovery of quantum materials coldea2001; shimizu2003; mendels2010. Quantum spin liquids show interesting topological and entanglement properties kalmeyer1987; kivelson1987; wen1989; read1991 which can be used for applications in quantum information kitaev2009. The Kitaev spin model on the honeycomb lattice kitaev2006 represents an important class of models, since it can be solved exactly in a Majorana fermion representation and demonstrates the significance of gauge fields on the low-energy properties. The model shows three gapped spin liquid phases carrying Abelian anyon excitations and an intermediate gapless phase which can be identified as a semi-metal of Majorana fermions. Applying a magnetic field along the three spatial directions kitaev2006, referring to a field in the direction, induces a gap in the intermediate phase, then producing a topological phase with non-Abelian anyons MooreRead closely related to a superconductor ReadGreen with chiral edge modes. It is important to emphasize that static spin-spin correlation functions are not sufficient to characterize the phase diagram of this model baskaran2007. Theoretical efforts have been performed to compute dynamical correlation functions kitaev2011; knolle2014; song2016; pollmann; kourtis as well as entanglement properties yao2010; pollmann. By bipartitioning a system spatially, the entanglement entropy measures how entangled the two subsystems are eisert2010. Related to these theoretical developments, quasi-two-dimensional quantum materials have been synthesized Kee; Trebst; Takagi, with recent measurements from neutron neutron and Raman scattering Nasu, nuclear magnetic resonance janvsa2018 and thermal transport nasu2017; kasahara20181; kasahara20182. One and two-dimensional spin systems, could also be engineered in quantum circuits and ultra-cold atoms circuit; Oreg; Duan.

In this Letter, we propose valence bond fluctuations as a probe of entanglement properties and correlations in the ground state of the Kitaev spin model.

A valence bond anderson1987 here corresponds to the spin-spin pairing between two nearest neighbor electrons. A relation between bipartite fluctuations of some observables such as charge or spin and entanglement properties of many-body Hamiltonians has been previously shown song2012; alex2014; loic2017; pollmann2014. This measure can also accurately detect quantum phase transitions between Néel and resonating valence bond spin phases stephan. Here, we show the relation between entanglement properties and valence bond fluctuations in the Kitaev spin model, through bipartite fluctuations associated to bond-bond correlation functions. In the three gapped phases, the valence bonds between nearest neighbors form a crystalline or dimer order kitaev2006. Approaching the transition(s) to the gapless intermediate phase, these bonds now resonate giving rise to gapless critical fluctuations. The valence bond fluctuations are distinct from the ones of a collection of resonating Einstein-Podolsky-Rosen pairs or Bell pairs alet2007, since the elementary block of the theory is a Majorana fermion. We present calculations in one and two dimensions on the honeycomb lattice, and check our mathematical findings with numerical calculations, e.g., through the Density Matrix Renormalization Group (DMRG). In one dimension, the gapless phase is reduced to a quantum critical point feng2007, which then develops into a plane for ladder systems PRBlong. In two dimensions, in the absence of a magnetic field, valence bond correlations in space shuo2008 show a similar scaling as the dynamical spin structure factor kitaev2011.

Model on the chain.— First, we address the quantum chain or wire model feng2007; PRBlong, shown in Fig. 1 (top). The Hamiltonian takes the form

(1) |

The sum acts on odd sites only such that is an integer with being the total number of unit cells. Applying the Jordan-Wigner transformation, one can then map quantum spins-1/2 onto spinless fermionic operators with occupancy and at each site: . We further consider the following representation with two Majorana fermions per site: for ; for . The key point is that in this basis all the Majorana fermions decouple from the chain and encode the double-degeneracy of the (spin) ground state on a given bond of nearest neighbors. A Majorana fermion then contributes to a entropy by analogy to the two-channel Kondo model affleck1991; 2CKME. Below, we address long-range quantum properties of the chain which are characterized the Majorana fermion Hamiltonian PRBlong

(2) |

To calculate valence bond correlation functions, it is then useful to introduce a bond operator acting in the middle of two sites, then forming a dual lattice.

Therefore, we introduce complex bond fermion operators on the dual lattice as . In momentum space, in the basis , then the Hamiltonian becomes

(3) |

Here , with the lattice spacing set to . The matrix has the two eigenvalues , reflecting a gap in the spectrum if . We check that the gap closes at when and that the free fermion chain on the dual lattice results in a power-law decay of correlation functions. We find its counterpart in the original spin basis through the observable “valence bond correlator”. We introduce the irreducible valence bond correlation function: where . The site index represents the -th unit cell of the sublattice . In the dual lattice, relates to the density of (free) fermions . From Wick’s theorem, therefore we get for , .

As a comparison, take the usual two-spin correlator . Decoupled Majorana fermions lead to . Applying the Wick’s theorem, vanishes in all phases, and similarly for and which involve Jordan-Wigner strings formed by the pairing of and Majorana fermions. Deviating from the gapless point, from the bond-fermion model and from the Ising symmetry of the spin chain, we predict that the correlation length is proportional to the inverse of the gap , with , and that for , otherwise. We then perform numerical calculations based on DMRG which verify these predictions with the associated critical exponent in the inset of Fig. 1 (bottom right).

To calculate the valence bond fluctuations, it is useful to define the two functions:

(4) |

While measures the fluctuations in subsystem only, also includes the correlations between and . There is an equality between the two quantities: . Below, we focus on , which then defines the valence bond fluctuations related to bipartite bond-bond correlations. We choose a bipartition of the spin chain into two parts and depicted in Fig. 1 (top). A direct lattice summation leads to SM

(5) |

where we take subsystem lengths and the Euler constant. On the contrary, always contains a higher order scaling in . Our findings are consistent with the DMRG results. At the gapless point shown in Fig. 1 (bottom left), we observe a logarithmic scaling of with a pre-factor , , which is in agreement with the free fermion representation in the dual lattice FS. Fig. 1 (bottom right) probes the gapped region. When , we check that goes to zero reflecting the crystallization of the dimers. Slowly closing the gap , as near the phase transition point, we check that exhibits a logarithmic behavior .

It is interesting to reveal now the relation between and the entanglement entropy eisert2010. Here, represents the reduced density matrix of the subsystem . In the limit , eigenstates are formed on strong -links and vanishes when the boundary is set on the weak -link (see Fig. 1, top). By increasing , long-range entanglement emerges with . The entropy reaches its maximum at the gapless point , on a finite critical chain under open boundary conditions and obeys the universal behaviour cardy2004; affleck1991: , with the central charge in Conformal Field Theory (CFT), the boundary entropy and a non-universal constant. In the inset of Fig. 1 (bottom left), from DMRG, we then check that the central charge . From Eq. (2), half of the spin degrees of freedom are disentangled by decoupling all -Majorana fermions.

Both and then share a logarithmic growth with subsystem size typical of (critical) conformal-field theories in one dimension. To evaluate the valence bond fluctuations we diagonalize the spectrum in momentum space in the basis, whereas the entanglement entropy reflects the real space degrees of freedom on the original lattice. While the gapless spectrum for complex bond fermions defined on the dual lattice has symmetry with the central charge , the original spin chain shares the central charge for -Majorana particles, reminiscent of the Kitaev honeycomb model kitaev2006; kasahara20182.

Model on the honeycomb lattice.— On the two-dimensional honeycomb lattice, due to the protection of -flux configurations in the ground state, the two-spin correlator also vanishes beyond nearest neighbors in all phases baskaran2007. However, the valence bond correlator preserves local flux pairs in neighboring plaquettes and supports gapless fermion excitations. It exhibits a power-law decay in the gapless phase and an exponential decay in gapped phases shuo2008. Here, we study the spatial dependence and anisotropy of the valence bond correlations to evaluate its fluctuations on a bipartite lattice. We also include the effect of a uniform magnetic field and will then compare the behavior of valence bond fluctuations with the one of the entanglement entropy.

We start from the Hamiltonian , under the perturbation , where represents two nearest-neighbor sites and denotes three different directions for the Ising couplings (see Fig. 2, top left). When , the leading order term in perturbation theory breaks time-reversal symmetry, and the effective Hamiltonian is simplified to kitaev2006

(6) |

Here and describes next-nearest neighbors and connected by site . On the honeycomb lattice, one can then map spin operators into four Majorana fermions per site: , with the constraint for a physical state . Fixing the gauge for in sublattice on Fig. 2 (top left), the effective Hamiltonian then presents a quadratic form in terms of as in one dimension kitaev2006:

(7) |

To obtain the eigenvectors, again it is more convenient to use the basis of complex bond fermions: with the unit cell coordinate defined on the -links and the sublattice index. In momentum space, one obtains a similar form as Eq. (3) with , . Here, and are functions defined on unit vectors and (shown in Fig. 2, top left): , . The gapless intermediate semi-metal phase satisfies with all permutations of .

On the honeycomb lattice, the valence bond correlation on the -links is defined through the operator: . From Wick’s theorem, we find

(8) |

with the total number of lattice unit cells. In the absence of magnetic field , has no singularities in the three gapped Abelian phases, therefore this results in an exponential decay of . In the intermediate gapless semi-metal phase, singularities appear at two Dirac points . In the supplementary material SM, we present a detailed analysis of the behavior of at long distances. Performing an expansion around the two Dirac points similar to the superconductor LoicPRB, we recover a power-law decay shuo2008, and establish that the coefficient depends on the anisotropic function :

(9) |

Here, is the angle between the vectors and . The space variable refers to .

Once the gapless intermediate phase is subject to a magnetic field along the direction, an identical power-law behavior (including the same angular dependence) emerges in the dynamical correlation function kitaev2011: . We check that at long distances and at a finite time , . Both observables are proportional to the density-density correlation function of the bond fermions . When , a gap opens in the intermediate phase and the valence bond correlator now reveals an exponential decay, with a sign change observed when increasing the strength of the magnetic field SM.

To gain some intuition on the behavior of bipartite fluctuations, we perform analytically the lattice summation assuming an isotropic form of SM. We then find that for a bipartition represented in Fig. 2 (top left) with the subsystem size , the system shows a scaling in and an area law in

(10) |

In a gapless phase, we obtain where is a constant for a given set of ’s in the -isotropic assumption. For a gapped phase, we obtain where is the correlation length. This approach then shows that must reveal a maximum when reaching a quantum phase transition from a gapped phase into the gapless intermediate regime. We then check these results numerically using Eq. (8) for the function . In Fig. 2 (top right), we recover the linear scaling of in the gapless phase (). The inset shows the scaling of , where the leading-order term () is dominated by the on-site bond fluctuations () consistent with the -isotropic approximation. It is further confirmed that the anisotropy effects in the -function are responsible for the decrease of when the system goes deeper into the gapless phase SM. Fig. 2 (bottom left) also includes the evolution of the linear scaling factor with different magnetic strengths. The signature of the divergence in across the phase transition line is robust against small fields.

Now, we address the behavior of the entanglement entropy in the Kitaev model. As pointed out in Ref. yao2010, the total entanglement entropy of the Kitaev honeycomb model consists of two pieces: the gauge field part and the fermionic contribution . Since the measurement of valence bond fluctuations preserves the gauge field structure, probes entanglement properties of the fermion sector. Fortunately, is responsible for all the essential differences between the Abelian and non-Abelian phases. We then extract the linear factor from following the methods of Refs. peschel2003; yao2010. We establish that shares the same response as across the phase transition lines shown in Fig. 2 (bottom right). Therefore, in two dimensions, we also find the relation . In Ref. SM, for completness, we provide some details on the numerical approach to evaluate the fermionic contribution to the entanglement entropy.

It is perhaps useful to compare the obtained behavior of bond-bond correlation functions from the ones of the two-dimensional Heisenberg model, i.e., of a Néel ordered phase subject to spin-wave excitations. When antiferromagnetic Heisenberg interactions are dominant, the modified spin-wave theory predicts that a staggered magnetic field is required to stabilize the Néel state at zero temperature for finite lattices takahashi1989; song2011. Performing a spin-wave analysis, then we find that the same valence bond correlation shows: with and . As a result, the bipartite fluctuations now follow a volume square law: , arising from the non-vanishing long-range correlation of . Measuring the precise leading order scalings then allows to probe the phase, Kitaev spin liquid versus Néel state, of a two-dimensional quantum material. We emphasize here that the entanglement entropy of the Néel state still reveals an area law song2011, as in the Kitaev spin model.

To summarize, we have found a general relation between the valence bond fluctuations and the entanglement entropy of the Kitaev spin model in one and two dimensions. Bond-bond correlation functions and valence bond fluctuations also appear as a relevant tool to identify phases and phase transitions of Majorana magnetic quantum systems. Generalisations to three-dimensional systems could be analyzed Maria; Mariareview.

Acknowledgements.— This work has benefitted from discussions with L. Henriet, L. Herviou, N. Laflorencie, C. Mora, F. Pollmann, S. Rachel, G. Roux, H.-F. Song, A. Soret, at the DFG meetings FOR 2414 in Frankfurt and Göttingen, at CIFAR meetings in Canada and at the conference in Montreal related to the workshop on entanglement, integrability, topology in many-body quantum systems. Support by the Deutsche Forschungsgemeinschaft via DFG FOR 2414 is acknowledged as well as from the ANR BOCA. KP acknowledges the support by the Georg H. Endress foundation. Numerical calculations performed using the ITensor C++ library, http://itensor.org/.

## References

## I Supplementary material

This additional material is designated to help the reader to re-derive the mathematical formulae and to check some numerical steps presented in the Letter. First, we provide more information on the derivation of Eq. (5) in the Letter regarding the one-dimensional case. Then, we give more detail on the two-dimensional situation for the bond-bond correlation function, and on the bipartite fluctuation scaling relations including a numerical analysis. We also justify the numerical approach for the entanglement entropy in two dimensions related to Fig. 2 in the Letter.

## I. Kitaev spin chain

Here, we give a proof of Eq. (5) in the Letter.

### A. Quantum critical point

For the critical Kitaev spin chain at , the bipartite fluctuation within subregion takes the form . The sums depend only on the difference of two variables and can be converted into a single sum through

(11) |

From the expression of : for ; for , we get

(12) |

where for and for . A finite sum can be written as the difference of two infinite sums

(13) |

For the second equality, we have used the relation: . We can apply the properties of the digamma function, which shares the series representation related to the Euler’s constant and the asymptotic expansion

(14) |

Therefore,

(15) |

Then for a bipartition , we get the scaling of and at the gapless point

(16) |

### B. Gapped regime

Now, we consider that the Kitaev spin chain is in the gapped phase with where and are negative such that the strong bonds occur on the -links. We assume that the bond correlators behave as for and for , where the correlation length satisfies . For ,

(17) |

Approximating the single summation by an integral, we obtain

(18) |

As , the third and fourth terms can be ignored and this results in

(19) |

## II. Kitaev honeycomb model: Bond-bond correlation functions

Here, we present a more detailed analysis of bond correlation functions for the Kitaev honeycomb model in different phases, with and without the uniform magnetic field .

### A.

For simplicity, we start from the gapless limit where three Ising couplings share the same strength: . When , the valence bond correlator can be simplified into the product of two sums

(20) |

Here, represents the number of unit cells in the system. The main contribution comes from the two Dirac points where . We can approximate the summation by an expansion around each Dirac point with a small radius : . For the first sum, around one Dirac point we get . Here is the angle between the relative vector around the Dirac cone and the axis. It is clear to see that is anisotropic. We denote the direction of the two unit cells as with the angle between vectors and . Then , with the relative angle between and . Now we can evaluate the summation by taking the continuum limit

(21) |

The factor comes from the contribution of two Dirac points, which can be verified as equal. The factor originates from a change of basis from in the Brillouin zone (with unit vectors and ) to . We have also used the relation . Via a change of variables , we reach

(22) |

where with a cutoff and inside the integral denotes the Bessel function of the first kind. The same expansion around the other Dirac point would give an additional phase factor and the total contribution reads

(23) |

For the second sum in , we change to and the relative angle from to :

(24) |

In the end, we recover the scaling of bond correlators in the gapless phase and the amplitude contains an anisotropic factor

(25) |

When the direction is perpendicular to , for instance , the spatial oscillations of the bond correlators disappear . As verified by Fig. 3 (left), it supports a smooth curve of revealing the scaling. In Eq. (25), the cutoff can be determined by setting the radius of the momentum integration and approximating to be the total system size . Analytically, we obtain with , . (Throughout the supplementary material, is equivalent to the natural logarithm with the base ). It is consistent with the numerical fitting results in Fig. 3 (left): for . Shifting to other directions, As a result, in Fig. 3 (right) the sampling points of along the direction oscillate. When the gapless phase is perturbed by a magnetic field along the direction, the irreducible dynamical correlation function in the long range and at a finite time kitaev2011 shares the identical power-law behavior (including the same angular dependence) as the static valence bond correlator:

(26) |

where is the strength of the magnetic field and the parameter can be estimated as . We find

(27) |

Both observables are proportional to the density-density correlation function of bond fermions.

Importantly, it can be checked that the forms of and the anisotropic -function in Eq. (25) are valid for the whole gapless region. Later, we will study these anisotropic effects on the bipartite fluctuations in relation with Fig. 5 (right). For the gapped phase, numerically follows an exponential decay with a fast decreasing correlation length shown in Fig. 3 (left). Also in Fig. 3 (middle), one observes less anisotropy effects in the gapped region.

### B.

As shown in Fig. 4, we further study the effects of a uniform magnetic field on the phase (). With a gap opening, the valence bond correlation function in Fig. 4 (left) shows an exponential decay, similar to three gapped phases, but it changes the sign from positive to negative when varying the strength of the magnetic field (see Fig. 4, right). Consequently, in Fig. 4 (middle) we see an enhancement in the amplitude of bond correlation functions once the magnetic field is sufficiently large. We find that this sign change originates from the competition between the Ising interactions and the external magnetic field. When , the valence bond correlator can be expressed in an alternative form as:

(28) |

While the Ising interactions give a positive contribution to the bond correlators, the external magnetic field gives a negative one. Changing the strength of the external magnetic field , it is verified that

(29) |

as and