Valence bond fluctuations in the Kitaev spin model

TL;DR

This paper introduces valence bond fluctuations as a new way to characterize quantum spin liquids and their phase transitions in the Kitaev model, revealing a universal scaling law and robustness under magnetic fields.

Contribution

It presents valence bond fluctuations as a novel tool to analyze quantum spin liquids, establishing their relation to entanglement entropy and phase transitions in the Kitaev model.

Findings

Valence bond fluctuations scale identically to entanglement entropy.

They can locate phase transitions between gapped and gapless phases.

The fluctuations remain a robust indicator under magnetic field effects.

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 111 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 Neel magnetic state.

Equations

$${\mathcal{S}}_A=-\text{Tr}{\rho }_A\ln {\rho }_A,$$

$$\mathcal{S}(\rho ,{\rho }^{\prime })=\text{Tr}\left[\rho (\log \rho -\log {\rho }^{\prime })\right],$$

$$\mathcal{S}(\rho ,{\rho }^{\prime })\ge \frac{1}{2}{\left||\rho -{\rho }^{\prime }|\right|}^2.$$

$${\mathcal{F}}_A={⟨{\left(\sum _{i\in A}{Q}_i\right)}^2⟩}_c,$$

$${\mathcal{F}}_{AB}=\left|\sum _{i\in A}\sum _{j\in B}{⟨Q_i{Q}_j⟩}_c\right|.$$

$${\mathcal{F}}_{AB}=\frac{1}{2}\left|{\mathcal{F}}_A+{\mathcal{F}}_B-{\mathcal{F}}_{A\cup B}\right|.$$

$$I(A,B)={\mathcal{S}}_A+{\mathcal{S}}_B-{\mathcal{S}}_{A\cup B}.$$

$$I(A,B)=\mathcal{S}({\rho }_{A\cup B},{\rho }_A\otimes {\rho }_B).$$

$$I(A,B)\ge \frac{{\left(⟨Q_A{Q}_B⟩-⟨Q_A⟩⟨Q_B⟩\right)}^2}{2{\left||Q_A|\right|}^2{\left||Q_B|\right|}^2}.$$

$$\frac{I(A,B)}{{\mathcal{F}}_{AB}}\ge \text{cst}.$$

$${\mathcal{F}}_A={\mathcal{F}}_B,{\mathcal{F}}_A+{\mathcal{F}}_B=2{\mathcal{F}}_A\ge 0.$$

$${\mathcal{S}}_A=\underset{K\to \infty }{\lim }\sum _{n=1}^{\lfloor (K+1)/2\rfloor }{\alpha }_{2n}(K)C_{2n}.$$

$$C_n={(-i{\partial }_{\lambda }{)}^n\ln \chi (\lambda )|}_{\lambda =0}.$$

$$\frac{{\mathcal{S}}_A}{{\mathcal{F}}_A}\simeq \text{cst}.$$

$$\text{cst}\cdot \frac{3}{{\pi }^2}=\{\begin{array}{ll}K,& \text{Luttinger liquids};\\ c/g,& \text{U(1) CFTs},\end{array}$$

$${\mathcal{F}}^{\text{TS}}$$

$${\mathcal{F}}^{\text{VB}}$$

$${\mathcal{S}}_{\text{2D HAF}}^{\text{VB}}\sim \{\begin{array}{ll}ax\ln x+bx& (\text{Neˊel})\\ b^{\prime }x& (\text{VB})\end{array}.$$

$$|⌢{⟩}_{\overline{(i,1)(j,2)}}=\frac{1}{\sqrt{2}}(|{\uparrow }_{i,1}⟩|{\downarrow }_{j,2}⟩-|{\downarrow }_{i,1}⟩|{\uparrow }_{j,2}⟩),$$

$$|{\Phi }_0⟩=\sum _p{\lambda }_p|{\phi }_p⟩,$$

$$|{\phi }_p⟩=\prod _{\overline{(i,1)(j,2)}\in p}|⌢{⟩}_{\overline{(i,1)(j,2)}}.$$

$${\mathcal{S}}^{\text{VB}}({\Phi }_0)=\frac{{\sum }_p{\lambda }_p{\mathcal{S}}^{\text{VB}}({\phi }_p)}{{\sum }_p{\lambda }_p}.$$

$${\mathcal{F}}^{\text{TS/VB}}({\Phi }_0)$$

$${\mathcal{S}}^{\text{VB}}({\Phi }_0)$$

$${\mathcal{F}}^{\text{TS/VB}}({\Phi }_0):=\sum _p|{\lambda }_p{|}^2{\mathcal{F}}^{\text{TS/VB}}({\phi }_p,{\phi }_p).$$

$${\mathcal{S}}^{\text{VB}}\propto \ln 2\cdot n,$$

$$\frac{{\mathcal{S}}^{\text{VB}}}{{\mathcal{F}}^{\text{TS}}+{\mathcal{F}}^{\text{VB}}}\ge \ln 2.$$

$${\mathcal{F}}_{\text{2D HAF}}^{\text{TS/VB}}=b^{\text{TS/VB}}\cdot x+\mathcal{O}(x).$$

$${\mathcal{F}}_{\text{Kitaev},{\mathbb{Z}}_2}^{\text{TS}}=0.$$

$$\frac{{\mathcal{S}}_F}{{\mathcal{F}}_{\text{Kitaev},{\mathbb{Z}}_2}^{\text{VB}}}\ge \text{cst}.$$

$${\mathcal{F}}_{\text{Kitaev, SU(2)}}^{\text{TS/VB}}=b^{\text{TS/VB}}\cdot x+\mathcal{O}(x).$$

$$\mathcal{H}=\sum _{j=2m-1}{J}_1{\sigma }_j^x{\sigma }_{j+1}^x+J_2{\sigma }_{j+1}^y{\sigma }_{j+2}^y,$$

$$c_j$$

Full text

I Introduction

The quest of quantum spin liquids in the Mott regime [1, 2, 3, 4, 5, 6] has been a great challenge these last decades in relation with the discovery of quantum materials [7, 8, 9, 10, 11, 12, 13, 14, 15]. Quantum spin liquids show interesting topological and entanglement properties [16, 17, 18, 19] which can be used for applications in quantum information [20]. The Kitaev spin model on the honeycomb lattice [21] represents an important class of models, since it can be solved exactly in a Majorana fermion representation and demonstrates the significance of $\mathbb{Z}_{2}$ 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 [21], referring to a field in the $[111]$ direction, induces a gap in the intermediate phase, then producing a topological $\mathbb{Z}_{2}$ phase supporting non-Abelian anyonic excitations [22] and closely related to a $p_{x}+ip_{y}$ superconductor [23] with chiral edge modes. It is important to emphasize that static spin-spin correlation functions are exactly zero beyond nearest neighbors in this model [24]. Theoretical efforts have been performed to compute dynamical correlation functions [25, 26, 27, 28, 29] as well as the entanglement entropy [30, 29]. By bipartitioning a system spatially, the entanglement entropy measures how entangled the two subsystems are [31]. Related to these theoretical developments, quasi-two-dimensional quantum materials have been synthesized [32, 33, 34], with recent measurements from neutron [35] and Raman [36] scatterings, nuclear magnetic resonance [37] and thermal transport [38, 39, 40]. One and two-dimensional Kitaev spin liquids, could also be engineered in ultra-cold atoms [41] and quantum circuits  [42, 43].

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

A valence bond (VB) [3] here corresponds to the spin-spin pairing between two nearest neighbor electrons. Our first insight comes from the system of SU(2)-symmetric quantum spins with resonating valence bonds (RVB) where we find the bond fluctuations can be related to valence bond entropy of Einstein-Podolsky-Rosen (EPR) pairs or Bell pairs [44]. Extending to the Kitaev spin liquids, in the three gapped phases, the valence bonds between nearest neighbors form a crystalline or dimer order [21]. Approaching the transition(s) to the gapless intermediate phase, these bonds now resonate giving rise to gapless critical fluctuations, which in principle encode information on quantum phase transitions and entanglement properties. Our calculations indeed reveal an identical scaling between valence bond fluctuations and entanglement entropy in one-dimensional chain and two-dimensional honeycomb lattice, and we 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 [45], which then develops into a plane for ladder systems [46]. In two dimensions, in the absence of a magnetic field, the long-range valence bond correlations in space [47] share a similar scaling as the dynamical spin structure factor [25]. We include also the effects of a uniform magnetic field in the perturbative regime, and discuss relevant consequences from the excitations of flux pairs [27] and from the formation of U(1) gapless spin liquids once three Ising couplings become anti-ferromagnetic (AFM) [48]. In the end, to make a closer link with quantum materials, we give a comparison of valence bond fluctuations in the Néel state favored by strong AFM Heisenberg exchanges.

II Review of Fluctuations as An Entanglement Probe

First, we begin with a brief review on the relation between bipartite fluctuations and entanglement entropy in many-body Hamiltonians characterized by different symmetries [49, 50, 51, 52]. In Sec. II.1, we define the general fluctuations on a bipartite lattice, which from the information theory, provide a lower bound for the mutual information related to entropy. We then remind in Sec. II.2, an exact expression for the entropy as a series expansion of even cumulants (with particle number or spin fluctuations the leading order) for the U(1) charge conserved systems [49] and an inequality between the two quantities emerging among the SU(2) quantum spins described by resonating valence bonds.

Generalizing these works to Kitaev spin liquids coupled to a gapped $\mathbb{Z}_{2}$ gauge field represents our central motivation in this work. The difficulty lies in finding the right observable encoding the long-range correlation of matter Majorana fermions in the gapless phase, hence the entanglement properties. Fortunately, by analogy to the RVB states in the three gapped phases, in Sec. II.3 we verify that the valence bond fluctuations represent non-vanishing lower bound for the entropy.

II.1 Generalities

We decompose a generic quantum system into two parts $A\cup B$. For the subsystem $A$, the entanglement is measured by the von Neumann entropy [31]

[TABLE]

where ${\rho}_{A}=\text{Tr}_{B}\rho$ represents the reduced density matrix of sub-system $A$. Once given two density matrices $\rho$ and $\rho^{\prime}$, the distance between two states can be probed by the relative entropy

[TABLE]

with a norm bound [53]

[TABLE]

Here the norm stands for $\left||\rho|\right|=\text{Tr}\sqrt{\rho^{\dagger}\rho}$ and we have assumed that $\hbox{Tr}\rho=\hbox{Tr}\rho^{\prime}=1$. Making an analogy with vectors, one may also write $S(\rho,\rho^{\prime})=S(\rho\parallel\rho^{\prime})$. For instance, for diagonal (density) matrices, each eigenvalue may refer to a coordinate along one direction. On the other hand, to evaluate the fluctuations, we introduce two measurements

[TABLE]

Here, $Q$ is a chosen operator for targeted systems: charge, particle number, one spin or two spins on a valence bond and $\left<Q_{i}Q_{j}\right>_{c}=\langle Q_{i}Q_{j}\rangle-\langle Q_{i}\rangle\langle Q_{j}\rangle$ denotes the reduced correlation function. It is easy to notice while $\mathcal{F}_{A}$ measures the fluctuations in subsystem $A$, $\mathcal{F}_{AB}$ covers the correlations between $A$ and $B$. There is an equality between the two quantities:

[TABLE]

An important finding so far established is to relate $\mathcal{F}_{AB}$ with $\mathcal{S}_{A}$ by mutual information [54]

[TABLE]

From the definition of relative entropy (2.2), the mutual information has an alternative expression

[TABLE]

Choosing any operator of the matrix form $Q=Q_{A}\otimes Q_{B}$ with $Q_{A}$ the bounded operator in region $A$ and applying the Schwarz inequality $\left||\rho|\right|\geq\text{Tr}(\rho Q)/\left||Q|\right|$ to the norm bound (2.3), one obtains [55]

[TABLE]

The numerator recovers $\mathcal{F}_{AB}$. Correspondingly, for $\mathcal{F}_{AB}\neq 0$, we arrive at

[TABLE]

Although the lower bound between the bipartite fluctuations and the mutual information is universal, it remains ambiguous what is the form of operator $Q$ one should choose for a given many-body system such that the fluctuations measured are non-vanishing. A second inquiry would be: under which circumstances we could reach the equality of (2.10) such that the fluctuations share the same scaling as the original entropy.

II.2 Exact relations and inequalities

In this subsection, we give two known examples where one can relate entanglement entropy directly to charge or spin fluctuations. We further show in Sec. II.2.2, for the SU(2)-symmetric RVB state, that the bond correlator is also a good option for operator $Q$.

II.2.1 Noninteracting fermions with conserved U(1) charge

Let us consider a system of non-interacting fermions with conserved total charge, or more precisely total number of particles: $\hat{N}_{\text{tot}}|\psi\rangle=N_{\text{tot}}|\psi\rangle$. At zero temperature, the ground state $|\psi\rangle$ of the total system becomes pure. Basic properties of the entanglement entropy then follow

• •

symmetric: $\mathcal{S}_{A}=\mathcal{S}_{B}$;

• •

subadditive: $\mathcal{S}_{A}+\mathcal{S}_{B}\geq\mathcal{S}_{A\cup B}=0$.

Without any calculation, one can already see the similarities between charge fluctuations and the entropy. If we take $\hat{Q}_{A}=\hat{N}_{A}$ in $\mathcal{F}_{A}$ (2.4), where $\hat{N}_{A}$ represents the number of particles in sub-region $A$, as a result of total charge conservation $\hat{N}_{A}-N_{A}=-(\hat{N}_{B}-N_{B})$, the fluctuations also inherit the symmetric and sub-additive characteristics

[TABLE]

In fact, one can relate the two quantities more rigorously by the cumulant expansion of entropy [49]

[TABLE]

Here, the coefficients are all positive and related to the unsigned Stirling number of the first kind: $\alpha_{n}(K)=2\sum_{k=n-1}^{K}S_{1}(k,n-1)/(k!k)$. The cumulants $C_{n}$ are given by the generating function $\chi(\lambda)=\langle e^{i\lambda\hat{N}_{A}}\rangle$ according to

[TABLE]

By definition, one verifies $\mathcal{F}_{A}=C_{2}$. For a Gaussian process, one can truncate the serie with $K=1$, but for non-Gaussian models one needs to check carefully the convergence of the series with the appropriate number of cumulants [49]. As a comparison, if we consider a Bell pair or an EPR pair, then this generally requires around 10 cumulants to reproduce the $\ln 2$ entropy. Although the equivalence of entropy and a complete set of even cumulants $(\ref{eq:sc})$ is unique to the systems with a mapping to non-interacting fermions, the general relation between entropy and fluctuations (2.14) can be further extended to the interacting one-dimensional (1D) critical systems that conserve total charge, and can be described by a Gaussian model through conformal field theory (CFT) or bosonization. In those cases, $\mathcal{S}_{A}$ can also be truncated by a $K=1$ upper-bound. Thus one gets

[TABLE]

The constant proves to encode rich information, for instance [49]

[TABLE]

where similarly to Refs. [49, 56], we also introduce the letter $K$ for the Luttinger parameter and $c$ represents the central charge in conformal field theory (CFT). Parameter $g=\pi v\kappa$ consists of the velocity $v$ and compressibility $\kappa=\partial n/\partial\mu$.

II.2.2 SU(2) quantum spins with EPR pairs

Intuitively, one may wonder what will happen when the system breaks U(1) charge conservation and when the system becomes higher dimensional, such as two-dimensional. Next, we give an example of the SU(2)-symmetric valence bond state [49]. Indeed, a direct correspondence in terms of inequality similar to relation $(\ref{eq:sf})$, subsists if we replace $\mathcal{F}_{A}$ with $\mathcal{F}_{AB}$, by analogy to mutual information (2.10). Furthermore, we would like to extend the result of Ref. [49] and show how the two-spin fluctuations and valence bond fluctuations capture different features of the entanglement entropy. Here, by “two-spin fluctuations”, we refer to fluctuations associated to spin-spin correlation functions. From Eq. (2.5), two-spin (TS) fluctuations and valence bond fluctuations between two subregions read:

[TABLE]

In this work, our aim is to show that the valence bond fluctuations are essential since they provide relevant information both for SU(2) and $\mathbb{Z}_{2}$ quantum spin liquids.

We consider the two-dimensional Heisenberg antiferromagnetic (HAF) model where arise two competing phases: a Néel state and a gapped VB state. In either configuration, the valence-bond entropy has proven to exhibit distinct behaviors [44]:

[TABLE]

Here $x$ denotes the length of the boundary between two subsystems. This measure can also accurately detect quantum phase transitions between Néel and RVB spin phases in quasi-one dimensional ladder systems [57]. For the moment, we focus on the fluctuations of the gapped VB state and will address the comparison with a Néel state later in Sec. V.2.

Suppose our system comprises $N$ sites on even and odd sublattices. Dimer coverings between different sublattices sharing the form

[TABLE]

minimize the energy from the antiferromagnetic Heisenberg interactions. The subscript $(i,\alpha)$ describes the site on $i$-th unit cell of the $\{\alpha=1,2\}$ sublattice. A singlet state $|\Phi_{0}\rangle$ can then be represented as a complex superposition of all possible dimer or pairing configurations

[TABLE]

For a given global pairing distribution $p$, the product state goes over all local internal dimers $\overline{(i,1)(j,2)}$.

The corresponding VB entanglement entropy (2.17) is defined as [44]

[TABLE]

Formally, acting on the singlet state $|\Phi_{0}\rangle$, the fluctuations and the VB entropy can be obtained from the decomposition (2.19)

[TABLE]

Eq. (2.21) indicates in general there is no simple correspondence between $\mathcal{F}^{\text{TS/VB}}$ and $\mathcal{S}^{\text{VB}}$. Yet, we may try to simplify (2.21) as

[TABLE]

For $\mathcal{F}^{\text{TS}}$, the relation (2.22) is exact owing to the sublattice symmetry of two-spin correlation functions. For $\mathcal{F}^{\text{VB}}$, it is a redefinition in the sense one counts the bond fluctuations inside each pairing pattern with probability $|\lambda_{p}|^{2}$ and at the same time, ignores the contributions from the overlaps of different pairing patterns. This redefinition is crucial for $\mathcal{F}^{\text{VB}}$ to resemble the behavior of the VB entanglement entropy.

We present in Appendix A a detailed analysis of the fluctuations in SU(2)-symmetric quantum spin systems.

On one hand, if the gapped VB state is composed of $N$-site singlets carrying equal weights ($\lambda_{p}=\lambda$), we have the following inequality reminiscent of the lower bound for mutual information (2.10)

[TABLE]

Here $n$ denotes the number of singlets that the boundary crosses. Both relations above take the equality $``="$ if the maximum resonating range $N\leq 4$. As $N\to\infty$, the system approaches the gapless critical point.

On the other hand, as soon as the singlet bonds decay exponentially with distance ($\lambda_{p}\sim e^{-r/\xi}$),

[TABLE]

A similar area law scaling is revealed in two types of fluctuations alongside the VB entanglement entropy (2.17).

II.3 Generalization to Kitaev $\mathbb{Z}_{2}$ spin liquids

Valence bond states in SU(2) spin systems and Kitaev $\mathbb{Z}_{2}$ spin liquids may be distinguishable from the form of correlation functions. For the former, the two-spin correlator follows an exponential decay in the gapped phase; for the latter, however, the static correlation between two spins becomes exactly zero beyond nearest neighbors [24].

In fact, the Kitaev honeycomb model is solved in the Majorana representation with one spin operator mapped onto the product of one matter and one gauge Majorana fermions [21]. Once acting on the ground state embedded with a static $\mathbb{Z}_{2}$ gauge field, the gauge Majorana fermion creates a pair of fluxes in two adjacent hexagons. It renders two-spin fluctuations irrelevant, if given an arbitrary boundary (not assigned on the same Ising links)

[TABLE]

One should resort to the bond-bond operator [45]. Since the excitation of a flux pair is annihilated simultaneously by the other spin on the same bond, valence bond fluctuations always give a relevant lower bound regardless of the boundary position

[TABLE]

In the equation above, we have rescaled the mutual information by Fermi entropy according to the area law of entanglement entropy [30].

Another observation comes consistently from the SU(2)-invariant Kitaev spin liquids [58]: there, the two-spin operator can be expressed solely in terms of matter Majorana fermions (preserving the gauge structure) and its correlation becomes non-vanishing. Similar to the Heisenberg anti-ferromagnet (2.24), the spin and bond fluctuations obey a linear growth in the gapped region

[TABLE]

We then conclude that both two-spin and valence bond fluctuations are appropriate as relevant probes of the entanglement entropy for the SU(2) spin systems, whereas for the Kitaev $\mathbb{Z}_{2}$ spin liquids only the valence bond fluctuations play a substantial role.

Below, we address valence bond fluctuations both for one-dimensional and two-dimensional Kitaev spin models. In Secs. III and IV, we prove how valence bond fluctuations (2.26) develop the same scaling as the entanglement entropy both for the one-dimensional chain and the honeycomb lattice.

III Model on the chain

Here, we address the quantum chain or wire model [45, 46], shown in Fig. 1 (a). The Hamiltonian takes the form

[TABLE]

The sum acts on odd sites only such that $1\leq m\leq M$ is an integer with $M$ being the total number of unit cells.

III.1 Valence bond correlator

Applying the Jordan-Wigner transformation, one can then map quantum spins-1/2 onto spinless fermionic operators with occupancy $1$ and [math] at each site: $\sigma_{j}^{+}=a_{j}^{\dagger}\cdot{\prod_{i<j}(-\sigma_{i}^{z})},\sigma_{j}^{-}=a_{j}\cdot{\prod_{i<j}(-\sigma_{i}^{z})},\sigma_{j}^{z}=2a_{j}^{\dagger}a_{j}-1$. We further consider the following representation with two Majorana fermions per site:

[TABLE]

The key point is that in this basis all the $d$ Majorana fermions decouple from the chain and encode the double-degeneracy of the (spin) ground state on a given bond of nearest neighbors. A $d$ Majorana fermion contributes to a $(\ln 2)/2$ entropy by analogy to the two-channel Kondo model [59, 60]. Below, we address long-range quantum properties of the chain captured by the $c_{j}$ Majorana fermion Hamiltonian [46]

[TABLE]

To calculate valence bond correlation functions, it is useful to introduce a complex bond fermion operator acting in the middle of two sites $(2m-1,2m)$

[TABLE]

thus forming a dual lattice with site index $m$. In momentum space, choosing the basis $\Psi^{\dagger}=(\psi_{k}^{\dagger},\psi_{-k})$, the Hamiltonian becomes

[TABLE]

The matrix elements read

[TABLE]

From now on, we set the lattice spacing $l$ of the original chain to unity $1$. Matrix $\widetilde{M}$ has two eigenvalues

[TABLE]

reflecting a gap in the spectrum if $J_{1}\neq J_{2}$. We check that the gap closes at $k_{F}=\pi/2$ when $J_{1}=J_{2}$ and that the chain on the dual lattice (3.5) results in a critical gapless theory of free fermions with central charge $c=1$. We find its counterpart in the original spin basis through the observable “valence bond correlator”.

Using definitions in Sec. II.1, we introduce the bond-bond correlation functions $I(i,j)=\langle Q_{i}Q_{j}\rangle_{c}$, with here:

[TABLE]

It is important to underline that in Fig. 1 (a), we have chosen the strong bonds associated to the $J_{1}$ coupling, referring to the $x$ spin Pauli operator in $Q_{j}$.

As before, the site index $(j,\alpha)$ represents the $j$-th unit cell of the sublattice $\alpha=\{1,2\}$. In the dual lattice, $Q_{j}$ relates to the density of bond fermions $\psi^{\dagger}_{j}\psi_{j}$. At the gapless point, we get from Wick’s theorem

[TABLE]

for $i\neq j$. As a comparison, take the usual two-spin correlator $\langle\sigma_{i}^{z}\sigma_{j}^{z}\rangle=(-1)^{i+j+1}\left<c_{i}d_{i}c_{j}d_{j}\right>$. Decoupled $d$ Majorana fermions lead to $\langle c_{i}d_{j}\rangle=\langle d_{i}d_{j}\rangle=0$. Applying Wick’s theorem, $\langle\sigma_{i}^{z}\sigma_{j}^{z}\rangle$ vanishes in all phases beyond nearest neighbors, and the same for $\langle\sigma_{i}^{x}\sigma_{j}^{x}\rangle$ and $\langle\sigma_{i}^{y}\sigma_{j}^{y}\rangle$ which involve Jordan-Wigner strings formed by the pairings of $c$ and $d$ Majorana fermions. Once the distance of two sites goes beyond the nearest neighbour $|i-j|>1$, one verifies

[TABLE]

Like the Kitaev honeycomb model, in its one-dimensional chain analogue, again the two-spin operator does not encode the long-range correlation of the gapless Majorana fermions.

Now 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 $\xi$ of the bond operator is proportional to the inverse of the gap $\Delta=|J_{2}-J_{1}|$,

[TABLE]

The valence bond correlations share the behavior

[TABLE]

We then perform numerical calculations based on DMRG which verify these predictions with the associated critical exponent $\nu=0.94\sim 1$ in the inset of Fig. 1 (c).

III.2 Results on fluctuations

III.2.1 Logarithmic growth

Below, we focus on $\mathcal{F}_{AB}$, which defines the fluctuations between two subsystems (2.5) associated to bond-bond correlations (3.9). Here the valence bond fluctuations $\mathcal{F}_{AB}$ also appear as the effective non-vanishing lower bound for mutual information (2.10).

Fig. 1 (a) depicts the bipartition we choose for the spin chain with subsystem lengths $l_{A}=l_{B}={L_{\text{total}}}/2=L$. A direct lattice summation in Appendix B leads to

[TABLE]

with $\gamma\simeq 0.57721$ the Euler constant. On the contrary, $\mathcal{F}_{A}$ always contains a higher order scaling linear in $l_{A}$.

Our findings (3.14) are confirmed by DMRG simulations. At the gapless point $J_{1}=J_{2}$ shown in Fig. 1 (b), we observe in $\mathcal{F}_{AB}$ a logarithmic scaling with respect to the length of subregion $A$: $\mathcal{F}_{AB}\propto\ln l_{A}$. Roughly, one can identify this term by taking the 1D integral of the bond correlation $I(i,j)$ in Eq. (3.10). Moreover, the pre-factor $\alpha/\pi^{2}$ is recovered with $\alpha=0.95$. Here, the fact that $\alpha$ reproduces central charge $c=1$ of the dual lattice, is in agreement with the free bond fermion representation [61].

Fig. 1 (c) further probes the gapped region. When $|J_{1}|\gg|J_{2}|$, we check that $\mathcal{F}_{AB}$ goes to zero reflecting the crystallization of the dimers. Slowly closing the gap $\Delta$, near the phase transition point, $|c_{1}|\gg|c_{2}|$. we check that the logarithmic behavior $\propto\ln\xi$ dominates in $\mathcal{F}_{AB}$.

III.2.2 Relation to entanglement entropy

It is interesting to go beyond the lower bound and reveal the relation between $\mathcal{F}_{AB}$ and the original entanglement entropy. Deep in the gapped phase driven by $|J_{1}|\gg|J_{2}|$, eigenstates are formed on strong $x$-links. $\mathcal{S}_{A}$ vanishes accordingly when the boundary is set on the weak $y$-link (see Fig. 1, a). By increasing $|J_{2}|$, long-range entanglement emerges among the dimers, which is accompanied by a logarithmic growth in entropy associated to the correlation length

[TABLE]

The same response is observed in valence bond fluctuations $\mathcal{F}_{AB}$ of the gapped region.

Meanwhile, the entropy reaches its maximum when the gap closes at $J_{1}=J_{2}$. Suppose the critical chain is finite with open boundaries, the entropy is proven to show the universal behaviour [62, 59]:

[TABLE]

where $g$ counts the boundary entropy and $s_{1}$ stands for a non-universal constant. In inset of Fig. 1 (b), from DMRG, we check that the central charge extracted from the entropy (3.16) reads: $c=0.49\sim 1/2$. It can be understood from the fact that after the Jordan-Wigner transformation, half of the spin degrees of freedom are disentangled from the Hamiltonian (3.3) by decoupling all $d$-Majorana fermions.

At the gapless point, both $\mathcal{F}_{AB}$ and $\mathcal{S}_{A}$ then share a logarithmic growth with subsystem size typical of (critical) conformal field theories in one dimension. Related to this finding, we would like to address the following comment: to evaluate the valence bond fluctuations we diagonalize the spectrum in momentum space in the $\psi$ basis, whereas the entanglement entropy reflects the real space degrees of freedom on the original lattice. This justifies why in our calculations the central charge $c=1$ is revealed in the valence bond fluctuations, whereas the central charge $c=1/2$ is observed in the entanglement entropy.

We establish that in both the gapped and gapless phases of one-dimensional Kitaev spin liquids, there is an identical scaling rule between the valence bond fluctuations and the entanglement entropy

[TABLE]

IV Model on the honeycomb lattice

As discussed briefly in Sec. II.3 and as will be described below, on the two-dimensional honeycomb lattice, due to the protection of [math]-flux configurations in the ground state, the two-spin correlator also vanishes beyond nearest neighbors in all phases [24]. The valence bond correlator, however, preserves flux pairs in neighboring plaquettes and supports gapless fermion excitations. Although Ref. [47] finds numerically that it exhibits a power-law decay in the gapless phase and an exponential decay in gapped phases, the valence bond correlator itself does not locate precisely the phase transition (see Fig. 2, a). It motivates us to develop an approach of evaluating its global fluctuations on a bipartite lattice. The enhanced features in fluctuations come intrinsically from the spatial dependence and anisotropy of the local bond correlator.

We demonstrate below that the valence bond fluctuations and the entanglement entropy allow us to locate quite accurately the phases and quantum phase transitions in the two-dimensional Kitaev honeycomb model [21].

IV.1 Majorana representation

We start from the Hamiltonian $\mathcal{H}=-\sum_{\langle ij\rangle_{a}}J_{a}\sigma_{i}^{a}\sigma_{j}^{a}$, under the perturbation $\mathcal{V}=-\sum_{j}h_{a}\sigma_{j}^{a}$, where $\langle ij\rangle$ represents two nearest-neighbor sites, forming the bonds and $a=x,y,z$ denotes one of three different Ising couplings assigned onto them (see Fig. 3, a). When $|h_{a}|\ll|J_{a}|$, the cubic term in perturbation theory breaks time-reversal symmetry, and the effective Hamiltonian is simplified to [21]

[TABLE]

Here $\kappa\sim h_{x}h_{y}h_{z}/(J_{a}J_{b})$ and $\langle\langle ik\rangle\rangle$ describes next-nearest neighbors $i$ and $k$ connected by site $j$. Below, we study the effect of $\kappa$ in the perturbative regime of Eq. (4.1) where $\kappa\ll J_{x}+J_{y}+J_{z}$.

On the honeycomb lattice, one can write spin operators in terms of four Majorana fermions per site [21]:

[TABLE]

with the constraint for a physical state

[TABLE]

Here, $c_{j}$ are matter Majorana fermion operators determining the fermionic spectrum and $c_{j}^{a}$ belong to gauge Majorana fermions which in the ground state, can be fixed artificially. Noticing that $[ic_{i}^{a}c_{j}^{a},\mathcal{H}]=0$ and $(ic_{i}^{a}c_{j}^{a})^{2}=1$, we choose the gauge $u_{\langle ij\rangle_{a}}=ic_{i}^{a}c_{j}^{a}:=+1$ for $i$ in sublattice $\{1\}$, as depicted in Fig. 3 (a), such that the total flux in the plaquette becomes zero consistent with Lieb’s theorem [63]. The definition of two sub-lattices, named $1$ and $2$ in Fig. 3 (a), follow the one-dimensional situation and the presence of two inequivalent bonds in a (zig-zag) given direction.

The effective Hamiltonian [21] then presents a quadratic form in terms of $c_{j}$ as in one dimension (3.3):

[TABLE]

To obtain the eigenvectors, again it is more convenient to use the basis of complex bond fermions:

[TABLE]

with $r$ the unit cell coordinate defined on the $z$-links and $\{1,2\}$ the sublattice index. In momentum space, one arrives at a similar form of Hamiltonian as Eq. (3.5) with matrix elements

[TABLE]

Here, $f(\vec{k})$ and $g(\vec{k})$ are functions of the form

[TABLE]

with two unit vectors $\vec{n}_{1}=(1/2,\sqrt{3}/2)$ and $\vec{n}_{2}=(-1/2,\sqrt{3}/2)$ shown in Fig. 3 (a). We check that the energy spectrum obtained from Eqs. (3.7), (4.6) and (4.7) reveals a gapless intermediate semi-metal phase  [21]

[TABLE]

where $(a,b,c)$ involve all permutations of $(x,y,z)$.

IV.2 Valence bond correlator

For the honeycomb lattice, we choose to measure the valence bond correlator onto the $z$-links, as compared to the $x$-links used for the one-dimensional chain in Sec. III:

[TABLE]

Since any physical observable is independent of a specific gauge [24], in the last equality, we can adapt our gauge choice for the gauge Majorana fermions. One thus sees the bond correlator probes the matter Majorana fermions without disturbing the gauge part.

On the contrary, a single spin operator $\sigma_{j,1}^{z}=ic_{j,1}c_{j,1}^{z}$ influences the gauge structure [24]. Defining the bond gauge fermion $\chi=(c_{j,1}^{z}+ic_{j,2}^{z})/2$, we find that the number operator $N_{\chi}=\chi^{\dagger}\chi=(u_{\langle j_{1}j_{2}\rangle_{z}}+1)/2$ takes the value $1$ and [math] depending on the gauge choice of $u_{\langle j_{1}j_{2}\rangle_{z}}=\pm 1$. In either gauge, $\sigma_{j,1}^{z}=ic_{j,1}(\chi^{\dagger}+\chi)$ changes the occupation number of the bond gauge fermion $\chi$. It is equivalent to flip the linking number $u_{\langle j_{1}j_{2}\rangle_{z}}$ to $-u_{\langle j_{1}j_{2}\rangle_{z}}$ and excite one $\pi$-flux pair in two neighbouring plaquettes. Therefore, the two-spin correlation is totally suppressed by the static $\mathbb{Z}_{2}$ gauge background beyond nearest neighbours (the latter case goes back to the measurement on a single bond).

Next, we focus on the local structure of the non-vanishing valence bond correlation. From Wick’s theorem, one obtains a compact form

[TABLE]

with $N$ the total number of lattice unit cells.

IV.2.1 Zero field

In the absence of magnetic field $\kappa=0$, $I(i,j)$ has no singularities in the three gapped Abelian phases, therefore this results in an exponential decay of $I(i,j)$. In the intermediate gapless semi-metal phase, singularities appear at two Dirac points $\pm\vec{k}^{*}$.

We first look at the behavior of bond correlations in the gapless region. A detailed analysis of the asymptotic behavior of $I(i,j)$ at long distances can be found in Appendix C. Performing an expansion around the two Dirac points similar to the $p_{x}+ip_{y}$ superconductor [64], we recover a power-law decay [47]

[TABLE]

and establish that the $\widetilde{c}_{1}$ coefficient depends on the cutoff function $t(\Lambda)$ and on the anisotropic function $Y(\vec{r})$. More precisely,

[TABLE]

Here, $\theta^{*}$ is the angle between the vectors $\vec{r}=\vec{r}_{i}-\vec{r}_{j}$ and $\vec{k}^{*}$. The space variable $r$ refers to $|\vec{r}_{i}-\vec{r}_{j}|$.

We can start from the simplest case by making the directions of $\vec{r}$ and $\vec{k}^{*}$ perpendicular to each other: $\vec{r}_{\perp}=\vec{r}_{j}-\vec{r}_{i}=(j-i)(\vec{n}_{1}+\vec{n}_{2})$. The spatial oscillations disappear in the bond correlator with $Y(\vec{r}_{\perp})=1$. Our analytic expression becomes consistent with the numerical fitting results of Ref. [47]. Shown in Fig. 2 (a), it supports a smooth curve of $I(i,j)$ revealing the $r^{-4}$ scaling in the gapless region ($J_{z}<0.5$).

We can derive a more precise analysis. At the gapless point $J_{x}=J_{y}=J_{z}=1/3$, in Appendix. C we derive the expression for the cutoff function

[TABLE]

where inside the integral $J_{1}(k)$ denotes the Bessel function of the first kind. Here, the cutoff $\Lambda=\xi r$ can be further approximated by setting the radius of the momentum integration $\xi=1$ and taking $r\simeq r_{\text{max}}$ to be the total system size $L=100$. Analytically, we obtain $\log|I(i,j)|=c_{1}-\alpha\log|i-j|$ with $c_{1}=\left.\log(\widetilde{c}_{1}/9)\right|_{\Lambda=100}=-4.63$, $\alpha=4$. Here, $``\log"$ is equivalent to the natural logarithm with the base ${e}$. It recovers well the numerical fitting result (see Fig. 2, a): $c_{1}=-4.60,\alpha=4.06$.

It is important to stress that Ref. [47] has not pointed out the role of the anisotropic $Y$-function. Once shifted to other directions $Y(\vec{r})\neq$ cst, $c_{1}=\log(|\cos^{2}(\vec{k}^{*}\cdot\vec{r})-\cos^{2}(\theta^{*})|)-4.63.$ Accordingly, as verified by Fig. 2 (b) and (c), the sampling points of $\log|I(i,j)|$ along the non-perpenticular direction oscillate rapidly. We also emphasize here the forms of $\tilde{c}_{1}$ and the anisotropic $Y$-function in Eq. (4.12) remain true for the whole gapless region. Later, we will study these anisotropic effects on the bipartite fluctuations in relation with Fig. 4.

For the gapped phase, on the other hand, from numerics $I(i,j)$ follows an exponential decay with a fast decreasing correlation length shown in Fig. 2 (a). Meanwhile, in Fig. 2 (b), one observes less anisotropy effects in the gapped region.

It may be relevant to mention that once the gapless intermediate phase is subject to a magnetic field along the $\hat{z}$ direction, an identical power-law behavior (including the same angular dependence) emerges in the dynamical correlation function [25]:

[TABLE]

where $h_{z}$ is the strength of the magnetic field and the parameter $h_{0}$ can be estimated as $h_{0}\sim J$. We find at long distances and at a finite time,

[TABLE]

Both observables are proportional to the density-density correlation function of the bond fermions $\psi_{r}$ in (4.5).

IV.2.2 Small finite field

We further study the effects of a small uniform magnetic field on the bond correlation in the intermediate phase. For simplicity, we take $J_{x}=J_{y}=J_{z}=J=1/3$.

When $0<\kappa\ll J$, a gap opens and the valence bond correlator in Fig. 2 (e) now reveals an exponential decay, similar to three gapped spin liquid phases. Yet its sign changes from positive to negative when increasing the strength of the magnetic field (see Fig. 2, d). Consequently, in Fig. 2 (f) we observe 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. For $\kappa\neq 0$, the valence bond correlator (4.10) can be expressed in an alternative form

[TABLE]

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 $\kappa$, it is verified that

[TABLE]

as ${\partial_{\kappa}I(i,j)}=2[{\partial_{\kappa}F(\vec{r})}\cdot F(\vec{r})-{\partial_{\kappa}G(\vec{r})}\cdot G(\vec{r})]$ and

[TABLE]

When $\kappa>0$, the derivative of $I(i,j)$ is always negative.

The monotonically decreasing bond correlation function is expected to cross zero around the point where the strengths of the Ising interactions and the magnetic field are comparable. We can roughly estimate the crossover point by starting from a relatively small $\kappa$ parameter. In this circumstance, $I(i,j)$ is still governed by an expansion $|\delta\vec{k}|\in\Omega(0,1)$ around two original Dirac points $\pm\vec{k}^{*}$. The denominator in Eq. (4.16) turns out to be

[TABLE]

When the parameter $\lambda=J/(6\kappa)>1$, $\kappa<\kappa_{c}=0.055$, the $F(\vec{r})^{2}$ term arising from the Ising interactions is dominant and $I(i,j)$ keeps a positive sign. Otherwise $\lambda\ll 1$, $\kappa\gg\kappa_{c}$, then the $-G(\vec{r})^{2}$ term from the external magnetic field grows steadily and has the tendency to drive $I(i,j)$ negative. In accordance with the numerical calculations shown in Fig. 2 (d) and (f), for different distances, all crossover points where $I(i,j)$ changes sign are located at $\kappa>\kappa_{c}$.

IV.3 Results on fluctuations

IV.3.1 Area law

Next, to gain some intuition on the behavior of bipartite fluctuations, in Appendix D we perform analytically the lattice summation by assuming an isotropic form of $I(i,j)$, namely with $Y(\vec{r})=1$.

Given a bipartition on the honeycomb lattice represented in Fig. 3 (a), we first derive a general scaling form for fluctuations within an arbitrary region $\Omega=l_{x}\times l_{y}$

[TABLE]

with $l_{x}$ and $l_{y}$ of the same order as $L$.

For the Kitaev honeycomb model, from Sec. IV.2 we see the valence bond correlator reveals a power law decay ($\alpha=4$) in the gapless phase and an exponential decay ($\alpha\to\infty$) in the gapped phases. Therefore, in all phases $\mathcal{F}_{A}$ shows the volume law: $\mathcal{F}_{A}\propto L^{2}=\mathcal{V}$. As usual, we can extract $\mathcal{F}_{AB}$ from the equality (2.6): $\left|\mathcal{F}_{A\cup B}-\mathcal{F}_{A}-\mathcal{F}_{B}\right|/2$. With a subsystem size $A=B={(L/2)}\times L$, it is noticeable that the volume term vanishes after the subtraction, leading to an area law in $\mathcal{F}_{AB}\propto L=\mathcal{A}$, where $\mathcal{A}$ refers then to an area.

Meanwhile, under the $Y$-isotropic form assumption, we establish in Appendix D the linear scaling factor of valence bond fluctuations in different phases:

[TABLE]

In the gapless phase, we obtain ${\alpha}_{\mathcal{F}}=3.84\tilde{c}_{1}$ where $\tilde{c}_{1}$ denotes the constant coefficient in the bond correlator (4.11) for a given set of $J_{a}$’s. In a gapped phase, we obtain ${\alpha}_{\mathcal{F}}\propto\xi^{3}$ with $\xi$ the correlation length. This approach then implies that with a rapidly growing correlation length, ${\alpha}_{\mathcal{F}}$ must reach a maximum when undergoing a quantum phase transition from a gapped phase into the gapless intermediate regime (see Fig. 4, a).

Numerically, we check these results by the method of finite-size scaling which starts from the anisotropic form (4.10) of the function $I(i,j)$. The exact anisotropic numerical calculations agree well with our previous $Y$-isotropic form approximation. In Fig. 3 (b), we recover the linear scaling of $\mathcal{F}_{AB}$ in the gapless phase ($J_{x}=J_{y}=J_{z}$). The inset shows the scaling of $\mathcal{F}_{A}$, where the leading-order $L^{2}$ term ($0.391$) is dominated by the on-site bond fluctuations ($0.362$, from Eq. (D13)). Since these on-site contributions are later subtracted, $\mathcal{F}_{AB}$ contains more information about the entanglement properties.

IV.3.2 Peak structure in linear scaling factor

Shown in Fig. 4 (a), we continue our study by extracting numerically the linear scaling factor of valence bond fluctuations from different regimes of the phase diagram. A peak structure centered at the quantum phase transition line is observed. While the gapped region can be understood from the simultaneous evolution with the correlation length ($\alpha_{\mathcal{F}}\propto\xi^{3}$), we check that the anisotropy effects in the $Y$-function are responsible for the decrease of $\alpha_{\mathcal{F}}$ when the system goes deeper into the gapless phase.

It is noted that after the double summation in $\mathcal{F}_{AB}$ (2.5), the bond fluctuations around the boundary (or domain wall) between two subsystems become the major contribution. Therefore, we can focus on the short-range behavior of the anisotropic factor $Y(\vec{r})$ in the bond correlator (4.12) along the direction perpendicular to the boundary. Illustrated by Fig. 4 (b), at short distances, $Y(\vec{r})$ reaches the maximum value when $J_{z}$ evolves to the phase transition line. Consequently, the amplitude of the linear scaling ${\alpha}_{\mathcal{F}}$ in $\mathcal{F}_{AB}$ would drop when we decrease $J_{z}$ in the gapless phase.

Fig. 4 (c) also includes the development of the linear scaling factor ${\alpha}_{\mathcal{F}}$ with different magnetic strengths. The signature of the peak structure in ${\alpha}_{\mathcal{F}}$ across the phase transition line is robust against small fields ($\kappa\leq 0.10$). By increasing $\kappa$, the gap is enlarged. The anisotropic effects of the $Y(\vec{r})$ function originally dominant in the gapless region become reduced significantly, thus making the cusp of $\alpha_{F}$ more smooth. Stronger magnetic field effects are discussed qualitatively in Sec. V.1.

IV.3.3 Relation to Fermion entropy

Now, we address the behavior of the entanglement entropy in the Kitaev model. As pointed out in Ref. [30], the total entanglement entropy of the Kitaev honeycomb model consists of two pieces: the gauge field part $\mathcal{S}_{G}=(L-1)\ln 2$ and the fermionic contribution

[TABLE]

Since the measurement of valence bond fluctuations preserves the $\mathbb{Z}_{2}$ gauge field structure, $\mathcal{F}_{AB}$ probes the entanglement properties of the fermion sector. Fortunately, $\mathcal{S}_{F}$ is responsible for all the essential differences between the Abelian and non-Abelian phases.

We then extract the linear factor ${\alpha}_{\mathcal{S}}$ from $\mathcal{S}_{F}$ following the methods of Refs. [65, 30]. For completeness, we refer the readers to Appendix E where some details on the numerical approach to evaluate the Fermi entropy are presented. It is found that in a small magnetic field $\kappa=0.01$, the linear scaling factor ${\alpha}_{\mathcal{S}}$ shares the same response as ${\alpha}_{\mathcal{F}}$ across the phase transition line shown in Fig. 4 (d). On top of that, once the magnetic field strength is increased, the peak structure of $\alpha_{\mathcal{S}}$ in Fermi entropy disappears slightly more quickly than the one in bond fluctuations. In addition, it is relevant to observe that $\alpha_{\mathcal{S}}$ remains zero for $J_{z}>0.5$. Another interesting observation is that as a function of $J_{z}$, the magnitudes of $\alpha_{\mathcal{F}}$ and $\alpha_{\mathcal{S}}$ vary approximately in the same range $[0,0.10]$.

To summarize, in two dimensions, we also find that valence bond fluctuations and the entanglement entropy of the Fermi sector show the same area law scaling in all phases:

[TABLE]

Moreover, their linear scaling factors act as signatures to characterize quantum phase transitions between the Abelian and non-Abelian phases in the Kitaev honeycomb model.

V Concluding remarks

In this final section, we would like to make brief remarks about the effects of stronger magnetic fields on the gapless phase of the Kitaev honeycomb model. Two scenarios are presently discussed in the literature: one with excitations of flux pairs [37, 27], and the other with a transition to another type of gapless spin liquid with U(1) symmetry [48]. We will propose possible responses from the valence bond fluctuations respectively.

Afterwards, a comparison with the Néel state supported by antiferromagnetic Heisenberg interactions will give us additional insights with regard to the application in quantum materials.

V.1 Influence of Perturbations

V.1.1 Perturbed Kitaev QSLs with gauge-flux pairs

For the Kitaev materials in a gapless spin liquid state, two types of gaps have been observed in the presence of a small tilted magnetic field [37]

[TABLE]

Here, $\Delta_{v}$ denotes the gap from the creation of a pair of fluxes (or visons) and $\Delta_{f}\propto h_{x}h_{y}h_{z}$ refers to the one induced by one matter Majorana fermion excitation. Our previous analysis of Sec. IV remains valid as long as $\Delta_{v}\gg\Delta_{f}$.

In Ref. [27], it was suggested that one can construct an exact perturbed ground state with even number of virtual fluxes by a unitary mapping $U$ from the unperturbed state $|\varphi_{0}\rangle$: $|\varphi\rangle=U|\varphi_{0}\rangle$. The transformed spin operator takes the form [27]

[TABLE]

with $Z=1,f=0$ for the pure Kitaev model and $Z<1,f\neq 0$ in the presence of perturbations. The nonzero $f$ parameters open a Majorana fermion gap instantaneously.

One immediately notices for the valence bond operator $Q_{i}=\sigma_{i,1}^{z}\sigma_{i,2}^{z}$, that the leading order contribution turns into

[TABLE]

Concerning the gauge Majorana fermions $c_{i}^{a}$, we have used the gauge convention form $u_{\langle i_{1}i_{2}\rangle_{z}}=+1$ and all others being zero, acting on the unperturbed state. The first term in the transformed bond operator (5.3) has a rescaling factor $Z^{2}$. It leads to a $r^{-4}$ decay in valence bond correlation $I(i,j)=\langle Q_{i}Q_{j}\rangle_{c}$ within the distance shorter than the correlation length ($r<\xi$). The second part contains products of four matter Majorana fermions, which then result in much faster decays in bond correlations, for instance $r^{-6}$ and $r^{-8}$ over short distances.

Neglecting these higher order corrections arising from the $f$-decomposition and taking into account the general scaling rule on honeycomb lattice (4.20), we establish that the valence bond fluctuations still show an area law

[TABLE]

The linear scaling factor is rescaled according to

[TABLE]

where $\alpha_{\mathcal{F},0}$ denotes the prefactor of the linear term reminiscent of the zero-flux Kitaev spin liquids. Based on Sec. IV.3, we thus find

[TABLE]

With excitations of flux pairs, the linear scaling factor now combines two pieces of information: the vison gap $\Delta_{v}$ determining the amplitude of $Z^{4}$ and the Majorana fermion gap $\Delta_{f}$ coming into play through $\alpha_{\mathcal{F},0}$.

It may be relevant to mention that the two-spin fluctuations become already nonzero in the perturbed limit. From Eq. (5.2), one gets contributions from the $f$-sector of Majorana fermions. Accompanied by an exponential decay in the spin-spin correlation, we obtain

[TABLE]

where

[TABLE]

When $\Delta_{v}\gg\Delta_{f}$, no excitation of fluxes is allowed. When $f^{z}\to 0$, $\alpha^{\text{TS}}_{\mathcal{F}}\to 0$, we check the result of vanishing two-spin fluctuations in the solvable limit (2.25).

V.1.2 Transition to U(1) gapless spin liquids

If one continues to increase the strength of the uniform magnetic field, from numerical simulations [48], while the Kitaev ferromagnet produces a trivial polarized phase (PL), the Kitaev antiferromagnet might give rise to an intermediate phase supporting U(1) gapless spin liquids (GSL).

In the PL phase, there is no correlation between two subsystems and both the two-spin and valence bond fluctuations vanish

[TABLE]

For the GSL phase in the Kitaev antiferromagnet, one can assume a gapless spinon Fermi surface coupled to a U(1) gauge field. In an effective picture of complex Abrikosov fermions [48], a spin operator is mapped onto the product of fermions $2\vec{S}=f_{\alpha}^{\dagger}\vec{\sigma}_{\alpha\beta}f_{\beta}$, with spin index $\alpha,\beta=\uparrow,\downarrow$ and a U(1) symmetry $f_{\alpha}^{\dagger}\to e^{i\theta}f_{\alpha}^{\dagger}$. From this perspective, the spin and bond correlations follow power-law decays ($r^{-4}$, $r^{-8}$ respectively) in the gapless phase. We predict that an “enhancement” might be observed in the prefactor of the area-law bipartitie fluctuations

[TABLE]

The existence of a similar peak structure in $\alpha_{\mathcal{F}}$ between the gapped Kitaev spin liquids and U(1) gapless spin liquids is possible and can be tested numerically in the future.

V.2 Comparison with the Néel phase

In the end, to make a closer link with quantum materials, 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 [66, 56].

Performing a spin-wave analysis in Appendix F, then we find that the same valence bond correlation shows:

[TABLE]

with $c_{0}=0.131$ and $c_{1}=0.141$.

As a result, the bipartite fluctuations now follow a volume square law:

[TABLE]

arising from the non-vanishing long-range correlation of $c_{0}$. 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 [56], as in the Kitaev spin model. The violation of the lower bound (2.10) originates from the finite-size regularization procedure taken in the modified spin-wave theory.

V.3 Conclusion

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. Valence bond fluctuations appear as a relevant tool to identify phases and phase transitions of Majorana magnetic quantum systems. Application to three-dimensional systems [67, 68] can be studied in the future.

Acknowledgements

We would like to thank two anonymous referees for their valuable suggestions on comparing different types of spin liquids and on developing the non-perturbative regime. This work has also benefitted a lot 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/.

Appendix A SU(2) quantum spin systems

Here, we evaluate the two-spin fluctuations and valence bond fluctuations for the SU(2)-symmetric quantum spin models both in one and two dimensions. We differentiate two cases: the one where the singlets resonate with equal weights and the other where they decay exponentially over distance.

A.1 Finite-size singlet state

We start from one dimension. Fig. 5 shows all the pairing states $\{|\varphi\rangle_{i}\}$ for $N=2,4,6$-site singlets. $N$ can also be interpreted as the maximum resonating range of the VB states on an infinite chain. As $N\to\infty$, the system approaches the critical point.

When $N=2$, we have a unique singlet state $|\Phi_{0}\rangle=|\varphi\rangle$ shown in Fig. 5 (a). Once the boundary between two subsystems $A$ and $B$ is set on the bond, one recovers immediately $\mathcal{F}^{\text{TS}}=1$, $\mathcal{F}^{\text{VB}}=0$, $\mathcal{S}^{\text{VB}}=\ln 2$ and the relation $\mathcal{S}^{\text{VB}}/\mathcal{F}^{\text{TS}}=\ln 2$. Nevertheless, if the boundary is off the bond, $\mathcal{F}^{\text{TS}}=\mathcal{F}^{\text{VB}}=\mathcal{S}^{\text{VB}}=0$. In the following analysis, we always keep the size of the two subsystems to be the same.

For $N=4$ in Fig. 5 (b), a simple normalized singlet state reads $|\Phi_{0}\rangle=(|\varphi_{1}\rangle+|\varphi_{2}\rangle)/\sqrt{3}$ with an overlap between two pairings $\langle\varphi_{1}|\varphi_{2}\rangle=1/2$. Also, we have assumed the system has equal chance to stay in each pairing state. Then

[TABLE]

We find a new relation

[TABLE]

Yet Eq. (A2) does not hold if $N>4$. We can check the case $N=6$ shown in Fig. 5 (c). The system allows $(N/2)!=6$ pairing configurations and the normalization factor $\mathcal{D}$ becomes larger than 6 on account of the overlaps:

[TABLE]

While each $|\varphi_{p}\rangle$ is responsible for $\mathcal{S}^{\text{VB}}$ and $\mathcal{F}^{\text{TS}}$, only $|\varphi_{p=4,5}\rangle$ contribute to $\mathcal{F}^{\text{VB}}$.

[TABLE]

Therefore, ${\mathcal{S}^{\text{VB}}}/({\mathcal{F}^{\text{TS}}+\mathcal{F}^{\text{VB}}})=({5}/{2})\ln 2>\ln 2$.

In two dimensions, for a gapped VB state composed of $N$-site singlets, the generalized result becomes

[TABLE]

$n$ denotes the number of singlets that the boundary crosses. Both relations above take the equality $``="$ if the resonating range $N\leq 4$. Meanwhile, we recover the area law of bipartite fluctuations

[TABLE]

where the coefficients depend on the size of the singlets. For instance, we have $f^{\text{TS}}(2)=1,f^{\text{VB}}(2)=0;f^{\text{TS}}(4)=2/3,f^{\text{VB}}(4)=1/3;f^{\text{TS}}(6)=5/9,f^{\text{VB}}(6)=1/9$.

A.2 General singlet state

Next, we consider a gapped system where the singlet bonds decay exponentially with distance, and $\mathcal{F}^{\text{TS}}$ and $\mathcal{F}^{\text{VB}}$ still show an area law as the entanglement entropy.

We first test the 1D case. A general VB state can be constructed as

[TABLE]

where $\xi$ stands for the correlation length of dimers. The probability of two sites $(i,1)$ and $(j,2)$ being paired reads

[TABLE]

Correspondingly, on a chain with a total length $2L$ ($1\ll\xi\ll L$),

[TABLE]

$\mathcal{F}^{\text{TS/VB}}_{\text{1D}}$ becomes a constant proportional to the square of the correlation length.

In two dimensions, one can apply the general scaling rule of Appendix D.1 on $\mathcal{F}^{\text{TS/VB}}=\sum_{i\in A,j\in B}I(i,j)$. Since the correlators $I(i,j)$ follow an exponential decay (as in one dimension), the bipartite fluctuations between two subsystems should be proportional to the perimeter $x$ of the boundary

[TABLE]

Hence, we find a similar area law in the two types of fluctuations as for the VB entanglement entropy in (2.17).

Appendix B Kitaev spin chain

Here, we give a mathematical proof of Eq. (3.14).

B.1 Quantum critical point $J_{1}=J_{2}$

For the critical Kitaev spin chain at the gapless point $J_{1}=J_{2}$, we first evaluate the bipartite fluctuation within subregion $A$:

[TABLE]

The bond correlation depends only on the difference of two variables. One can thus convert the double sum into a single sum through

[TABLE]

From the expression (3.10) of $I(k)$: $1/[{\pi^{2}}(k^{2}-1/4)]$ for $k\neq 0$; $1-4/\pi^{2}$ for $k=0$, we get

[TABLE]

with $\sigma_{k}=1$ for $k=\text{odd}$ and $\sigma_{k}=0$ for $k=\text{even}$. It is convenient to re-express the finite sum in terms of the difference of two infinite sums

[TABLE]

In the second equality, we have also used the relation: $(\text{odd terms})=(\text{all terms})-(\text{even terms})$. Now we can apply the properties of the digamma function, which shares the series representation related to the Euler’s constant $\gamma\simeq 0.57721$, as well as the asymptotic expansion

[TABLE]

Therefore,

[TABLE]

For a bipartition $l_{A}=l_{B}=L/2$, we then get the critical scaling of $\mathcal{F}_{A}$ and at the same time, $\mathcal{F}_{AB}$ from their relation (2.6)

[TABLE]

B.2 Gapped regime $|J_{1}|>|J_{2}|$

We continue to study the gapped phase of the Kitaev spin chain with negative exchanges such that $|J_{1}|>|J_{2}|$, such that the strong bonds occur on the $x$-links. In Eqs. (3.12) and (3.13), we predict that the bond correlator behaves as $I(i,j)=c_{1}|i-j|^{-2}$ for $|i-j|\leq\xi$ and $I(i,j)=c_{2}e^{-|i-j|/\xi}$ for $|i-j|>\xi$, with a correlation length $\xi\sim\left|J_{2}-J_{1}\right|^{-1}$. For $1<\xi<l_{A}=l_{B}=L/2$, the valence bond fluctuations between two subregions become

[TABLE]

Approximating the single summation by an integral and supposing $\xi\ll L$, one obtains

[TABLE]

Appendix C Asymptotic form of bond correlator in the intermediate gapless phase

Here, we derive the power-law behavior of bond correlation functions in the intermediate gapless spin liquid phase of the Kitaev honeycomb model. We assume a simple case where three Ising couplings share the same strength: $J_{x}=J_{y}=J_{z}=J$.

When $\kappa=0$, the valence bond correlator (4.10) can be re-expressed as the product of two sums

[TABLE]

The main contribution comes from the two Dirac points $\pm\vec{k}^{*}=\pm(k_{x}^{*},k_{y}^{*})=\pm({4\pi}/{3},0)$ which satisfy $|f(\pm\vec{k}^{*})|=0$. It allows us to approximate the summation by an expansion around each Dirac point within a small radius $\xi$: $\vec{k}\in\Omega(\pm\vec{k}^{*},\xi)$.

For the first sum in Eq. (C1), around one Dirac point $\vec{k}^{*}=(4\pi/3,0)$ we get

[TABLE]

Here $\theta^{\prime}$ is the angle between the relative vector around the Dirac cone $\delta\vec{k}=\vec{k}-\vec{k}^{*}$ and the $x$ axis. It is clear to see that $I(i,j)$ is anisotropic. To simplify the exponential, we denote the direction of the two unit cells as $\vec{r}=\vec{r}_{j}-\vec{r}_{i}=(r\cos{\theta^{*}},r\sin{\theta^{*}})$ with $\theta^{*}$ the angle between vectors $\vec{r}$ and $\vec{k}^{*}$. Then,

[TABLE]

with $\theta$ the relative angle between $\delta\vec{k}$ and $\vec{r}_{j}-\vec{r}_{i}$.

Now we can evaluate the summation by taking the continuum limit

[TABLE]

The factor $2$ comes from the contribution of two Dirac points. The other factor $\sqrt{3}/2$ originates from a change of basis from $dk_{1}dk_{2}$ in the Brillouin zone (with unit vectors $\vec{n}_{1}$ and $\vec{n}_{2}$) to $dk_{x}dk_{y}$. We have also used the relation $\theta+\theta^{\prime}=\theta^{*}$.

Via a change of variables $k^{\prime}=kr$, one reaches

[TABLE]

where $t(\Lambda)=\sqrt{3}/(2\pi)\cdot\int_{0}^{\Lambda}J_{1}(k)kdk$ with a cutoff $\Lambda=\xi r$ and inside the integral $J_{1}(k)$ denotes the Bessel function of the first kind.

A similar expansion around the other Dirac point $-\vec{k}^{*}=(-4\pi/3,0)$ would give an additional phase factor $e^{i(-\vec{k}^{*}\cdot\vec{r}-(\pi-\theta^{*}))}$. Thus, the total contribution reads

[TABLE]

For the second sum in Eq. (C1), we only need to change $\vec{r}$ to $-\vec{r}$ and adjust the relative angle from $\theta^{*}$ to $\theta^{*}-\pi$:

[TABLE]

Combining Eqs. (C1), (C6) and (C7), we then recover the $r^{-4}$ scaling of the bond correlator in the gapless phase [47]:

[TABLE]

Furthermore, from our calculations the amplitude $\widetilde{c}_{1}$ retrieves an anisotropic factor $Y(\vec{r})$:

[TABLE]

One can also verify that the forms of $\tilde{c}_{1}$ and of the anisotropic $Y$-function in Eq. (C9) are valid for the whole gapless region.

Appendix D Bipartite fluctuations on honeycomb geometry

We evaluate here the bipartite fluctuations on the honeycomb lattice, involving the lattice summation.

D.1 General scaling rule

Consider a bipartition on the honeycomb lattice shown in Fig. 3 (a). The parallelogram is expanded by two unit vectors $\vec{n}_{1}=(1/2,\sqrt{3}/2)$ and ${\vec{n}_{2}}=(-1/2,\sqrt{3}/2)$ with a total size $\Omega=l_{x}\times l_{y}$ and the subsystems are chosen as $A=B=(l_{x}/2)\times l_{y}={(L/2)}\times L$. For convenience, we adopt new coordinates $\vec{r}=x\vec{n}_{1}+y\vec{n}_{2}:x=1,2,\dots,l_{x},y=1,2,\dots,l_{y}$. The summation in the bipartite fluctuations can then be re-expressed into

[TABLE]

To derive general scaling arguments in a ‘simple’ way, we only consider the case where $I(\vec{r})$ is an isotropic function of the distance $|\vec{r}|=\sqrt{x^{2}+xy+y^{2}}$. By analogy to Eq. (B2), a relation between the double and single sums can be established

[TABLE]

and the same for $\sum_{y,y^{\prime}=1}^{l_{y}}I(x,y^{\prime}-y)$. Then the bipartite fluctuation function can be grouped into four parts:

[TABLE]

with

[TABLE]

The dominant scaling terms in $\mathcal{F}_{\Omega}$ depend on the particular form of $I(r)$. Suppose a general case where $I(r)\propto r^{-\alpha}$ and $l_{x}$, $l_{y}$ are of the same order as $L$:

[TABLE]

where $\mathcal{O}({1}/{L^{0}})\sim\mathcal{O}(\ln L)$. The leading-order scaling in $\mathcal{F}_{\Omega}$ becomes

[TABLE]

When $I(r)\propto e^{-r/\xi}$, $\alpha\to\infty$, $\mathcal{F}_{\Omega}$ still show the volume law: $\mathcal{F}_{\Omega}\propto L^{2}=\mathcal{V}$. Besides, after the subtraction: $\left|\mathcal{F}_{A\cup B}-\mathcal{F}_{A}-\mathcal{F}_{B}\right|/2$, while the higher order terms $L^{4}$ and $L^{3}$ survive in $\mathcal{F}_{AB}$ , the square term $L^{2}$ always vanishes, which leads to an area law in $\mathcal{F}_{AB}\propto L=\mathcal{A}$.

To evaluate $\mathcal{F}_{AB}$ more precisely, we are going to study next-leading order terms in $\mathcal{F}_{\Omega}$ case by case.

D.2 Kitaev model: the gapless phase

For the gapless phase of the Kitaev honeycomb model, first we consider an isotropic form of the valence bond correlator in Eq. (4.11)

[TABLE]

where ${\tilde{c}_{1}}$ is a constant for a given $J_{z}\in[0,0.50]$ in our convention.

As $\alpha=4$, from the general scaling rule in Eq. (D5) we get $I_{i}(i=1,2,3)\propto\mathcal{O}(1)$ and $I_{4}\propto\mathcal{O}(\ln L)$. In particular, due to the convergence of the lattice summations in $I_{i}(i=1,2,3)$ in Eq. (D4), we can introduce a cutoff $l_{x}=l_{y}=\Lambda=10^{3}$ to approximate these pre-factors:

[TABLE]

Here $I(0)$ denotes the on-site contribution and has the value

[TABLE]

$\langle{Q}_{i}\rangle=\langle\sigma^{z}_{i_{1}}\sigma^{z}_{i_{2}}\rangle$ represents the integral over the Brillouin zone:

[TABLE]

with $\cos\theta(k_{1},k_{2})=a/\sqrt{a^{2}+b^{2}}$ and $f(k)=a+bi=J_{x}e^{ik_{1}}+J_{y}e^{ik_{2}}+J_{z}$. Combined with the general expression of $\mathcal{F}_{\Omega}$ in Eq. (D3), we arrive at

[TABLE]

where

[TABLE]

In the special case $J_{x}=J_{y}=J_{z}=1/3$, we have

[TABLE]

As indicated by the inset of Fig. 3 (b), the pre-factor of the $L^{2}$ term in $\mathcal{F}_{A}$ takes the value $\alpha=0.391$ close to $I(0)/2$. We conclude that the major contribution to $\mathcal{F}_{A}$ comes from the on-site interactions.

D.3 Kitaev model: the gapped phases

In the gapped phases, the bond correlation function decays exponentially and there is less anisotropy observed in Fig. 2 (b). We can safely start with an isotropic form

[TABLE]

For $\alpha=\infty$ in Eq. (D5), all of the amplitudes $I_{i}(i=1,2,3,4)\propto\mathcal{O}(1)$. Yet we can still relate them to the powers of the finite correlation length $\xi$. The following assumption is taken

[TABLE]

such that when $n<5$,

[TABLE]

Back to Eq. (D4), we can then replace the summations on the lattice with integrals

[TABLE]

The pre-factors now read

[TABLE]

Correspondingly, the bipartite fluctuations share the quadratic and linear forms respectively

[TABLE]

with

[TABLE]

It indicates that in gapped phases, ${\alpha}_{\mathcal{F}}$ increases at the transition towards the intermediate gapless phase.

Appendix E Entanglement entropy of the Fermion sector

In this part, we extract the linear factor ${\alpha}_{\mathcal{S}}$ in $S_{F}$ applying the same numerical approach as in Refs. [65, 30]. The general idea is to map the total system $A\cup B$ on the honeycomb lattice of Fig. 3 (a) onto a torus geometry depicted in Fig. 6 (a). Periodic boundary conditions (PBC) are imposed along $X$ and $Y$ directions and the boundaries between subsystems now turn into two circles. Fig. 6 (a) also illustrates one of the two zigzag boundaries located on $x$-links. For this shape of the boundary, the generic Hamiltonian of subsystem $A$ reads

[TABLE]

where $\langle\ \ \rangle$ denotes the nearest-neighbor links and $\langle\langle\ \ \rangle\rangle$ the next-nearest-neighbor links. The gauge choice $u_{\langle jk\rangle_{a}}=+1,-u_{\langle ji\rangle_{a}}u_{\langle ik\rangle_{b}}=+1$ is manifested in the direction of the arrows from site $k$ to site $j$ in Fig. 6 (a).

Now along the $X$ direction, we go to the momentum space ($k_{x}=2m\pi/N_{x},m=0,1,\dots,N_{x}-1$) and Fourier transform the Hamiltonian (E1) into

[TABLE]

in the basis of the matter Majorana fermions $\mathcal{C}^{\dagger}(k_{x})=\left(c_{-k_{x},1},c_{-k_{x},2}\right)$. The matrix takes the form

[TABLE]

with matrix elements $r=-J_{x}$, $s=(J_{y}+J_{z})\cos(k_{x}/2)-i(J_{y}-J_{z})\sin(k_{x}/2)$, $\alpha=2\kappa\sin k_{x}$, $\beta=2\kappa\sin(k_{x}/2)$. It should be noted that the intrinsic structure of our $2N_{y}\times 2N_{y}$ Hamiltonian matrix $M(k_{x})$ is distinct from Ref. [30] where only one type of next-nearest-neighbor couplings $J^{\prime}$ is included. In our case, all three types of next-nearest-neighbor couplings between the matter Majorana fermions are taken into account, as they arise naturally from the effects of the uniform magnetic field [21]. Additionally, to be consistent with the bipartition we choose for $\mathcal{F}_{AB}$ in Fig. 3 (a), the boundary position has been switched from $z$-links in Ref. [30] to $x$-links.

The diagonalized Hamiltonian reads

[TABLE]

where $n=1,\dots,2N_{y}$ and $\psi_{nk_{x}}$ represents the standard complex fermionic annihilation operators. For a free fermion system, the energy spectrum $\xi_{n}(k_{x})$ can be used to calculate the entanglement entropy [65]

[TABLE]

Here $\lambda_{n}(k_{x})=(e^{\beta\xi_{n}(k_{x})}+1)^{-1}$ denote the eigenvalues of the single-particle correlation function $\langle\psi^{\dagger}_{nk_{x}}\psi_{n^{\prime}k_{x}}\rangle$. It follows a Fermi-Dirac distribution with an inverse temperature $\beta=1/(k_{B}T)$. Fig. 6 (b) shows the numerical entanglement spectrum $\lambda_{n}(k_{x})$ for the non-Abelian and Abelian phases. We observe two gapless branches in the non-Abelian phase. These two modes are responsible for the peaks in entanglement entropy plotted in Fig. 6 (c). Both features disappear in the Abelian phase.

Since $\mathcal{S}_{F}=\sum_{k_{x}}\mathcal{S}(k_{x})$, through the summation of $\mathcal{S}(k_{x})$ in the momentum space followed by a finite-size scaling with respect to $N_{x}$, we are able to obtain the pre-factor $\alpha_{\mathcal{S}}$ of the linear term in $\mathcal{S}_{F}$. The results are shown in Fig. 4 (d) in the Article.

Appendix F Heisenberg antiferromagnet on honeycomb lattice

First, we give a review of the modified spin-wave theory on a two-dimensional Heisenberg antiferromagnet. We then analyze the asymptotic behavior of the two-spin correlation function on the honeycomb lattice. Later, a closed form of the valence bond correlator is derived and the $r^{-1}$ scaling is verified both analytically and numerically.

F.1 Modified spin-wave theory

For an antiferromagnetic Heisenberg model on the honeycomb lattice, we can reach the Néel state by applying a staggered magnetic field [66, 56].

The Hamiltonian reads

[TABLE]

Two sets of sublattices $\{1,2\}$ are differentiated by $(-1)^{|\vec{r}|}=1$ for $\vec{r}\in 1$ and $(-1)^{|\vec{r}|}=-1$ for $\vec{r}\in 2$. Vectors $\vec{\delta}$ connect each site to its nearest neighbors, with the total number of nearest neighbor sites denoted by $z$. The introduction of the staggered field breaks the $O(3)$ spin-rotational symmetry and helps to repair the divergence in Green functions arising from the zero mode (or the Goldstone mode).

In the modified spin-wave theory [66], one can map spin operators to bosonic operators: $\vec{r}\in 1$, $S_{r}^{+}=\sqrt{2S-a_{r}^{\dagger}a_{r}}\ a_{r}$, $S_{r}^{-}=a_{r}^{\dagger}\sqrt{2S-a_{r}^{\dagger}a_{r}}$, $S_{r}^{z}=S-a_{r}^{\dagger}a_{r}$; $\vec{r}\in 2$, $S_{r}^{+}=b_{r}^{\dagger}\sqrt{2S-b_{r}^{\dagger}b_{r}}$, $S_{r}^{-}=\sqrt{2S-b_{r}^{\dagger}b_{r}}\ b_{r}$, $S_{r}^{z}=-\left(S-b_{r}^{\dagger}b_{r}\right)$ with $[a_{r},a_{r^{\prime}}^{\dagger}]=[b_{r},b_{r^{\prime}}^{\dagger}]=\delta_{r,r^{\prime}}$. Expansion around large $S$ then gives the Hamiltonian of the order $\mathcal{O}(S^{2})$ and $\mathcal{O}(S)$.

Combining Fourier transform and Bogoliubov transformation, it is straightforward to obtain single-particle expectation values

[TABLE]

with

[TABLE]

and others all vanish. Here, $N$ denotes the total number of lattice sites. $\eta$ and $\gamma_{k}$ are functions depending on the geometry of the lattice

[TABLE]

Let us take a closer look at a finite honeycomb lattice. The geometric function reads

[TABLE]

When $h\to 0$, $\eta\to 1$, there exists one zero mode $\vec{k}_{0}=(0,0)$ making $f(\vec{r})$ and $g(\vec{r})$ divergent.

Following Ref. [66], one repairs the divergence by adjusting the strength $h$ of the local staggered magnetic field such that the magnetization becomes zero

[TABLE]

It is noted that only the zero mode is regularized by $h\neq 0$ and the sum over the remaining region can be safely approximated by a finite integral at $h=0$. One arrives at

[TABLE]

Taking $S=1/2$ in Eq. (F6), we get

[TABLE]

In the same manner,

[TABLE]

F.2 Asymptotic behavior of spin-spin correlation

Now, we can evaluate the behavior of the two-spin correlation function. From Wick’s theorem, it differs between sites

[TABLE]

To restore the spin-rotational symmetry at zero magnetic field, we have introduced an extra factor $1/3$ to Eq. (F10).

At large distances, an expansion of $f(\vec{r})$ around $\vec{k}_{0}=(0,0)$ within radius $\xi$ leads to

[TABLE]

Once we take $r\to\infty$,

[TABLE]

with $c_{0}=m_{0}\simeq 0.246$, $c_{1}=3/(8\pi)\simeq 0.119$. Numerically, we check in Fig. 7 (a) the finite-size scaling of the single-particle expectation function $f(\vec{r})$ from Eq. (F3). The coefficients read $c_{0}=0.245$, $c_{1}=0.12$ with a power-law index $\alpha=1.04$. They agree well with the asymptotic form.

Since the approximation $(-\eta\gamma_{k})\simeq-1$ is still valid for $\vec{k}\in\Omega\left(0,\xi\right)$, one finds over a long distance

[TABLE]

The two-spin correlation function then reveal a power-law $r^{-1}$ decay with alternating signs on the same and different sublattices.

F.3 Closed form of the valence bond correlator

Next, we study the response of valence bond correlation functions in the Néel state. We adopt the same definition as the Kitaev honeycomb model

[TABLE]

The bond index $\langle i_{1}i_{2}\rangle$ denotes two sites in the $i$-th unit cell of the sublattice $\{1\}$ and $\{2\}$. In the modified spin-wave theory, $S_{i_{1}}^{z}=1/2-a_{i}^{\dagger}a_{i}$, $S_{i,2}^{z}=-(1/2-b_{i}^{\dagger}b_{i})$. Reassembling different terms, we reach

[TABLE]

To simplify the notation, we have introduced a new set of number operators

[TABLE]

First, taking into account $\langle a_{i}^{\dagger}a_{i}\rangle=\langle b_{i}^{\dagger}b_{i}\rangle=-{1}/{2}+f(\vec{0})={1}/{2}$, the single-particle expectation values become

[TABLE]

.

Then, Wick’s theorem can be applied to calculate the remaining terms involving (2,3,4)-particle expectation values. Before proceeding, from Eq. (F3) we identify three useful functions

[TABLE]

All the other terms like $\langle a_{i}^{\dagger}a_{j}^{\dagger}\rangle$ and $\langle a_{i}^{\dagger}b_{j}\rangle$ vanish in the ground state. Two-particle expectation values then read

[TABLE]

The three-particle expectation value takes the form

[TABLE]

The summation is invariant under the exchange of $(i\leftrightarrow j)$. Combined with the sublattice symmetries, we have

[TABLE]

For the four-particle expectation value, we verify

[TABLE]

Therefore, we arrive at a closed form of the bond-bond correlators

[TABLE]

Considering the asymptotic behaviors of $f(\vec{r})$ and $g(\vec{r})$ functions in Eqs. (F12) and (F13), we establish

[TABLE]

where $c_{0}=4m_{0}^{2}(g_{0}^{2}+m_{0}^{2})\simeq 0.131$ and $c_{1}={m_{0}(3g_{0}^{2}+6m_{0}^{2})}/{\pi}\simeq 0.141$. Fig. 7 (b) confirms numerically the long-range behavior of the valence bond correlator with a power-law index $\alpha=1.07$. The two pre-factors $c_{0}=0.13,c_{1}=0.155$ are also consistent with our analytical predictions.

In the end, we find for the Néel state supported by strong antiferromagnetic Heisenberg exchanges, that the valence bond correlator gives a signature of the $r^{-1}$ scaling accompanied by a non-vanishing constant from finite-size effects (the regularization of the zero mode). This is clearly distinct from the pure $r^{-4}$ scaling in the gapless Kitaev spin liquid phase.