Valence bond fluctuations in the Kitaev spin model

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.


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 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 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 [29,30]. 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 arXiv:1901.03973v2 [cond-mat.str-el] 8 Dec 2019 many-body Hamiltonians characterized by different symmetries [49][50][51][52]. In Sec. II A, 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 B, 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 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 C we verify that the valence bond fluctuations represent non-vanishing lower bound for the entropy.

A. Generalities
We decompose a generic quantum system into two parts A ∪ B. For the subsystem A, the entanglement is measured by the von Neumann entropy [31] where ρ A = Tr B ρ represents the reduced density matrix of sub-system A. Once given two density matrices ρ and ρ , the distance between two states can be probed by the relative entropy with a norm bound [53] S(ρ, ρ ) ≥ Here the norm stands for ||ρ|| = Tr ρ † ρ and we have assumed that Trρ = Trρ = 1. Making an analogy with vectors, one may also write S(ρ, ρ ) = S(ρ ρ ). 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

5)
Here, Q is a chosen operator for targeted systems: charge, particle number, one spin or two spins on a valence bond and Q i Q j c = Q i Q j − Q i Q j denotes the reduced correlation function. It is easy to notice while F A measures the fluctuations in subsystem A, F AB covers the correlations between A and B. There is an equality between the two quantities: An important finding so far established is to relate F AB with S A by mutual information [54] From the definition of relative entropy (2.2), the mutual information has an alternative expression (2.9) The numerator recovers F AB . Correspondingly, for F AB = 0, we arrive at 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.

B. 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 B 2, for the SU(2)-symmetric RVB state, that the bond correlator is also a good option for operator Q.
Without any calculation, one can already see the similarities between charge fluctuations and the entropy. If we takeQ A =N A in F A (2.4), whereN A represents the number of particles in sub-region A, as a result of total charge conservationN A − N A = −(N B − N B ), the fluctuations also inherit the symmetric and sub-additive characteristics In fact, one can relate the two quantities more rigorously by the cumulant expansion of entropy [49] (2.12) Here, the coefficients are all positive and related to the unsigned Stirling number of the first kind: α n (K) = 2 K k=n−1 S 1 (k, n−1)/(k!k). The cumulants C n are given by the generating function χ(λ) = e iλN A according to (2.13) By definition, one verifies 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 (2.12) 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, S A can also be truncated by a K = 1 upper-bound. Thus one gets 14) The constant proves to encode rich information, for instance [49] cst · 3 π 2 = K, Luttinger liquids; c/g, U(1) CFTs, (2.15) 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 = πvκ consists of the velocity v and compressibility κ = ∂n/∂µ.

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 twodimensional. 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 (2.14), subsists if we replace F A with F AB , by analogy to mutual information (2.10). Furthermore, we would like to extend the result of Ref. [49] and show how the twospin 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: 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 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]: 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 B. Suppose our system comprises N sites on even and odd sublattices. Dimer coverings between different sublattices sharing the form (2.20) Formally, acting on the singlet state |Φ 0 , the fluctuations and the VB entropy can be obtained from the decomposition (2.19) Eq. (2.21) indicates in general there is no simple correspondence between F TS/VB and S VB . Yet, we may try to simplify (2.21) as For F TS , the relation (2.22) is exact owing to the sublattice symmetry of two-spin correlation functions. For F VB , it is a redefinition in the sense one counts the bond fluctuations inside each pairing pattern with probability |λ p | 2 and at the same time, ignores the contributions from the overlaps of different pairing patterns. This redefinition is crucial for F 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 (λ p = λ), we have the following inequality reminiscent of the lower bound for mutual information (2.10) (2.23) Here n denotes the number of singlets that the boundary crosses. Both relations above take the equality " = " if the maximum resonating range N ≤ 4. As N → ∞, the system approaches the gapless critical point.
On the other hand, as soon as the singlet bonds decay exponentially with distance (λ p ∼ e −r/ξ ), (2.24) A similar area law scaling is revealed in two types of fluctuations alongside the VB entanglement entropy (2.17).

C. Generalization to Kitaev Z2 spin liquids
Valence bond states in SU(2) spin systems and Kitaev 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 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) 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 (2.26) 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 twospin 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 (2.27) 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 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 The sum acts on odd sites only such that 1 ≤ m ≤ M is an integer with M being the total number of unit cells.

A. Valence bond correlator
Applying the Jordan-Wigner transformation, one can then map quantum spins-1/2 onto spinless fermionic operators with occupancy 1 and 0 at each site: We further consider the following representation with two Majorana fermions per site: The key point is that in this basis all the d Majorana fermions decouple from the chain and encode the doubledegeneracy 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] 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) thus forming a dual lattice with site index m. In momentum space, choosing the basis Ψ † = (ψ † k , ψ −k ), the Hamiltonian becomes The matrix elements read From now on, we set the lattice spacing l of the original chain to unity 1. Matrix M has two eigenvalues reflecting a gap in the spectrum if J 1 = J 2 . We check that the gap closes at k F = π/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 A, we introduce the bondbond correlation functions I(i, j) = Q i Q j c , with here: 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, α) represents the j-th unit cell of the sublattice α = {1, 2}. In the dual lattice, Q j relates to the density of bond fermions ψ † j ψ j . At the gapless point, we get from Wick's theorem for i = j. As a comparison, take the usual two-spin 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 bondfermion model and from the Ising symmetry of the spin chain, we predict that the correlation length ξ of the bond operator is proportional to the inverse of the gap (3.12) The valence bond correlations share the behavior We then perform numerical calculations based on DMRG which verify these predictions with the associated critical exponent ν = 0.94 ∼ 1 in the inset of Fig. 1 (c).

Logarithmic growth
Below, we focus on F AB , which defines the fluctuations between two subsystems (2.5) associated to bondbond correlations (3.9). Here the valence bond fluctuations 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 total /2 = L. A direct lattice summation in Appendix B leads to with γ 0.57721 the Euler constant. On the contrary, 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 F AB a logarithmic scaling with respect to the length of subregion A: F AB ∝ 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 α/π 2 is recovered with α = 0.95. Here, the fact that α 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 | |J 2 |, we check that F AB goes to zero reflecting the crystallization of the dimers. Slowly closing the gap ∆, near the phase transition point, |c 1 | |c 2 |. we check that the logarithmic behavior ∝ ln ξ dominates in F AB .

Relation to entanglement entropy
It is interesting to go beyond the lower bound and reveal the relation between F AB and the original entanglement entropy. Deep in the gapped phase driven by |J 1 | |J 2 |, eigenstates are formed on strong x-links. 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 The same response is observed in valence bond fluctuations 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 [59,62]: where g counts the boundary entropy and s 1 stands for a non-universal constant. In inset of Fig. 1  At the gapless point, both F AB and 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 ψ 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 (3.17)

IV. MODEL ON THE HONEYCOMB LATTICE
As discussed briefly in Sec. II C and as will be described below, on the two-dimensional honeycomb lattice, due to the protection of 0-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].  Here we present the case Jx = Jy = Jz = 1/3. The direction of the relative vector is chosen on d(i, j) = (i − j)( n1 + n2).

A. Majorana representation
We start from the Hamiltonian where ij 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 | |J a |, the cubic term in perturbation theory breaks time-reversal symmetry, and the effective Hamiltonian is simplified to [21] Here κ ∼ h x h y h z /(J a J b ) and ik describes nextnearest neighbors i and k connected by site j. Below, we study the effect of κ in the perturbative regime of Eq.  in terms of four Majorana fermions per site [21]: with the constraint for a physical state Here, c j are matter Majorana fermion operators determining the fermionic spectrum and c a j belong to gauge Majorana fermions which in the ground state, can be fixed artificially. Noticing that [ic a i c a j , H] = 0 and (ic a i c a j ) 2 = 1, we choose the gauge u ij a = ic a i c a j := +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 onedimensional 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): To obtain the eigenvectors, again it is more convenient to use the basis of complex bond fermions: 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 Here, f ( k) and g( k) are functions of the form with two unit vectors n 1 = (1/2, √ 3/2) and n 2 = (−1/2, √ 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] where (a, b, c) involve all permutations of (x, y, z).

B. 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: 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 σ z j,1 = ic j,1 c z j,1 influences the gauge structure [24]. Defining the bond gauge fermion χ = (c z j,1 + ic z j,2 )/2, we find that the number operator N χ = χ † χ = (u j1j2 z + 1)/2 takes the value 1 and 0 depending on the gauge choice of u j1j2 z = ±1. In either gauge, σ z j,1 = ic j,1 (χ † + χ) changes the occupation number of the bond gauge fermion χ. It is equivalent to flip the linking number u j1j2 z to −u j1j2 z and excite one π-flux pair in two neighbouring plaquettes. Therefore, the two-spin correlation is totally suppressed by the static 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 nonvanishing valence bond correlation. From Wick's theorem, one obtains a compact form with N the total number of lattice unit cells.

Zero field
In the absence of magnetic field κ = 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 ± 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] 11) and establish that the c 1 coefficient depends on the cutoff function t(Λ) and on the anisotropic function Y ( r). More precisely, Here, θ * is the angle between the vectors r = r i − r j and k * . The space variable r refers to | r i − r j |. We can start from the simplest case by making the directions of r and k * perpendicular to each other: r ⊥ = r j − r i = (j − i)( n 1 + n 2 ). The spatial oscillations disappear in the bond correlator with Y ( r ⊥ ) = 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 where inside the integral J 1 (k) denotes the Bessel function of the first kind. Here, the cutoff Λ = ξr can be further approximated by setting the radius of the momentum integration ξ = 1 and taking r r max to be the total system size L = 100. Analytically, we obtain log |I(i, j)| = c 1 −α log |i−j| with c 1 = log( c 1 /9)| Λ=100 = −4.63, α = 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, α = 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 ( r) = cst, c 1 = log(| cos 2 ( k * · r) − cos 2 (θ * )|) − 4.63. Accordingly, as verified by Fig. 2 (b) and (c), the sampling points of log |I(i, j)| along the nonperpenticular direction oscillate rapidly. We also emphasize here the forms ofc 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 thê z direction, an identical power-law behavior (including the same angular dependence) emerges in the dynamical correlation function [25]: where h z is the strength of the magnetic field and the parameter h 0 can be estimated as h 0 ∼ J. We find at long distances and at a finite time, Both observables are proportional to the density-density correlation function of the bond fermions ψ r in (4.5).

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 < κ 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 κ = 0, the valence bond correlator (4.10) can be expressed in an alternative form (4.16) 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 When κ > 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 κ parameter. In this circumstance, I(i, j) is still governed by an expansion |δ k| ∈ Ω(0, 1) around two original Dirac points ± k * . The denominator in Eq. (4.16) turns out to be When the parameter λ = J/(6κ) > 1, κ < κ c = 0.055, the F ( r) 2 term arising from the Ising interactions is dominant and I(i, j) keeps a positive sign. Otherwise λ 1, κ κ c , then the −G( 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 κ > κ c .

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 ( 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 Ω = l x × l y with l x and l y of the same order as L.
For the Kitaev honeycomb model, from Sec. IV B we see the valence bond correlator reveals a power law decay (α = 4) in the gapless phase and an exponential decay (α → ∞) in the gapped phases. Therefore, in all phases F A shows the volume law: F A ∝ L 2 = V. As usual, we can extract F AB from the equality (2.6): With a subsystem size A = B = (L/2) × L, it is noticeable that the volume term vanishes after the subtraction, leading to an area law in F AB ∝ L = A, where 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: (4.21) In the gapless phase, we obtain α F = 3.84c 1 wherec 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 α F ∝ ξ 3 with ξ the correlation length. This approach then implies that with a rapidly growing correlation length, α 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 Yisotropic form approximation. In Fig. 3 (b), we recover the linear scaling of F AB in the gapless phase (J x = J y = J z ). The inset shows the scaling of 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, F AB contains more information about the entanglement properties.

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 (α F ∝ ξ 3 ), we check that the anisotropy effects in the Y -function are responsible for the decrease of α F when the system goes deeper into the gapless phase.
It is noted that after the double summation in 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 ( r) in the bond correlator (4.12) along the direction perpendicular to the boundary. Illustrated by Fig. 4 (b), at short distances, Y ( r) reaches the maximum value when J z evolves to the phase transition line. Consequently, the amplitude of the linear scaling α F in 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 α F with different magnetic strengths. The signature of the peak structure in α F across the phase transition line is robust against small fields (κ ≤ 0.10). By increasing κ, the gap is enlarged. The anisotropic effects of the Y ( r) function originally dominant in the gapless region become reduced significantly, thus making the cusp of α F more smooth. Stronger magnetic field effects are discussed qualitatively in Sec. V A.

Relation to Fermion entropy
Now, we address the behavior of the entanglement entropy in the Kitaev model. As pointed out in Ref. [30], Since the measurement of valence bond fluctuations preserves the Z 2 gauge field structure, F AB probes the entanglement properties of the fermion sector. Fortunately, S F is responsible for all the essential differences between the Abelian and non-Abelian phases. We then extract the linear factor α S from S F following the methods of Refs. [30,65]. 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 κ = 0.01, the linear scaling factor α S shares the same response as α 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 α S in Fermi entropy disappears slightly more quickly than the one in bond fluctuations. In addition, it is relevant to observe that α S remains zero for J z > 0.5. Another interesting observation is that as a function of J z , the magnitudes of α F and α 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: 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 [27,37], 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.

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] Here, ∆ v denotes the gap from the creation of a pair of fluxes (or visons) and ∆ f ∝ 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 ∆ v ∆ 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 |ϕ 0 : |ϕ = U |ϕ 0 . The transformed spin operator takes the form [27] U † σ a i U = iZc i c a i + · · · + f a ijk ic j c k + · · · , (5. One immediately notices for the valence bond operator Q i = σ z i,1 σ z i,2 , that the leading order contribution turns into Concerning the gauge Majorana fermions c a i , we have used the gauge convention form u i1i2 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) = Q i Q j c within the distance shorter than the correlation length (r < ξ). 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 The linear scaling factor is rescaled according to where α F ,0 denotes the prefactor of the linear term reminiscent of the zero-flux Kitaev spin liquids. Based on Sec. IV C, we thus find With excitations of flux pairs, the linear scaling factor now combines two pieces of information: the vison gap ∆ v determining the amplitude of Z 4 and the Majorana fermion gap ∆ f coming into play through α 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 When ∆ v ∆ f , no excitation of fluxes is allowed. When f z → 0, α TS F → 0, we check the result of vanishing twospin fluctuations in the solvable limit (2.25).

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 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 S = f † α σ αβ f β , with spin index α, β =↑, ↓ and a U(1) symmetry f † α → e iθ f † α . From this perspective, the spin and bond correlations follow powerlaw 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 (5.10) The existence of a similar peak structure in α F between the gapped Kitaev spin liquids and U(1) gapless spin liquids is possible and can be tested numerically in the future.

B. 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 [56,66].
Performing a spin-wave analysis in Appendix F, then we find that the same valence bond correlation shows: with c 0 = 0.131 and c 1 = 0.141. As a result, the bipartite fluctuations now follow a volume square law: 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.

C. 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.

Finite-size singlet state
We start from one dimension. Fig. 5 shows all the pairing states {|ϕ 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 → ∞, the system approaches the critical point.
When N = 2, we have a unique singlet state |Φ 0 = |ϕ shown in Fig. 5 (a). Once the boundary between two subsystems A and B is set on the bond, one recovers immediately F TS = 1, F VB = 0, S VB = ln 2 and the relation S VB /F TS = ln 2. Nevertheless, if the boundary is off the bond, F TS = F VB = S 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 |Φ 0 = (|ϕ 1 + |ϕ 2 )/ √ 3 with an overlap between two pairings ϕ 1 |ϕ 2 = 1/2. Also, we have assumed the system has equal chance to stay in each pair- ing state. Then We find a new relation S VB F TS + F VB = ln 2. (A2) 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 D becomes larger than 6 on account of the overlaps: While each |ϕ p is responsible for S VB and F TS , only |ϕ p=4,5 contribute to F VB .
In two dimensions, for a gapped VB state composed of N -site singlets, the generalized result becomes S VB ∝ ln 2 · n, S VB F TS + F VB ≥ ln 2.

General singlet state
Next, we consider a gapped system where the singlet bonds decay exponentially with distance, and F TS and F VB still show an area law as the entanglement entropy.
We first test the 1D case. A general VB state can be constructed as where ξ stands for the correlation length of dimers. The probability of two sites (i, 1) and (j, 2) being paired reads Correspondingly, on a chain with a total length 2L (1 ξ L), 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 F TS/VB = i∈A,j∈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 Hence, we find a similar area law in the two types of fluctuations as for the VB entanglement entropy in (2.17).

Quantum critical point J1 = J2
For the critical Kitaev spin chain at the gapless point J 1 = J 2 , we first evaluate the bipartite fluctuation within subregion A: The bond correlation depends only on the difference of two variables. One can thus convert the double sum into a single sum through i,j∈A From the expression (3.10) of I(k): 1/[π 2 (k 2 − 1/4)] for k = 0; 1 − 4/π 2 for k = 0, we get with σ k = 1 for k = odd and σ k = 0 for k = even. It is convenient to re-express the finite sum in terms of the difference of two infinite sums In the second equality, we have also used the relation: (odd terms) = (all terms) − (even terms). Now we can apply the properties of the digamma function, which shares the series representation related to the Euler's constant γ 0.57721, as well as the asymptotic expansion For a bipartition l A = l B = L/2, we then get the critical scaling of F A and at the same time, F AB from their relation (2.6) 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| ≤ ξ and I(i, j) = c 2 e −|i−j|/ξ for |i − j| > ξ, with a correlation length ξ ∼ |J 2 − J 1 | −1 . For 1 < ξ < l A = l B = L/2, the valence bond fluctuations between two subregions become Approximating the single summation by an integral and supposing ξ L, one obtains When κ = 0, the valence bond correlator (4.10) can be re-expressed as the product of two sums The main contribution comes from the two Dirac points ± k * = ±(k * x , k * y ) = ±(4π/3, 0) which satisfy |f (± k * )| = 0. It allows us to approximate the summation by an expansion around each Dirac point within a small radius ξ: k ∈ Ω(± k * , ξ).
For the first sum in Eq. (C1), around one Dirac point Here θ is the angle between the relative vector around the Dirac cone δ k = k− 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 r = r j − r i = (r cos θ * , r sin θ * ) with θ * the angle between vectors r and k * . Then, e i k·( rj − ri) = e i k * · r e iδ k· r = e i k * · r e iδkr cos θ , with θ the relative angle between δ k and r j − r i . Now we can evaluate the summation by taking the continuum limit The factor 2 comes from the contribution of two Dirac points. The other factor √ 3/2 originates from a change of basis from dk 1 dk 2 in the Brillouin zone (with unit vectors n 1 and n 2 ) to dk x dk y . We have also used the relation θ + θ = θ * .
Via a change of variables k = kr, one reaches where t(Λ) = √ 3/(2π) · Λ 0 J 1 (k)kdk with a cutoff Λ = ξr and inside the integral J 1 (k) denotes the Bessel function of the first kind.
A similar expansion around the other Dirac point − k * = (−4π/3, 0) would give an additional phase factor e i(− k * · r−(π−θ * )) . Thus, the total contribution reads For the second sum in Eq. (C1), we only need to change r to − r and adjust the relative angle from θ * to θ * − π: Combining Eqs. (C1), (C6) and (C7), we then recover the r −4 scaling of the bond correlator in the gapless phase [47]: Furthermore, from our calculations the amplitude c 1 retrieves an anisotropic factor Y ( r): One can also verify that the forms ofc 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.

General scaling rule
Consider a bipartition on the honeycomb lattice shown in Fig. 3 (a). The parallelogram is expanded by two unit vectors n 1 = (1/2, √ 3/2) and n 2 = (−1/2, √ 3/2) with a total size Ω = l x × l y and the subsystems are chosen as A = B = (l x /2) × l y = (L/2) × L. For convenience, we adopt new coordinates r = x n 1 + y n 2 : x = 1, 2, . . . , l x , y = 1, 2, . . . , l y . The summation in the bipartite fluctuations can then be re-expressed into To derive general scaling arguments in a 'simple' way, we only consider the case where I( r) is an isotropic function of the distance | r| = x 2 + xy + y 2 . By analogy to Eq. (B2), a relation between the double and single sums can be established The dominant scaling terms in F Ω depend on the particular form of I(r). Suppose a general case where I(r) ∝ r −α and l x , l y are of the same order as L: where O(1/L 0 ) ∼ O(ln L). The leading-order scaling in F Ω becomes When I(r) ∝ e −r/ξ , α → ∞, F Ω still show the volume law: F Ω ∝ L 2 = V. Besides, after the subtraction: |F A∪B − F A − F B | /2, while the higher order terms L 4 and L 3 survive in F AB , the square term L 2 always vanishes, which leads to an area law in F AB ∝ L = A.
To evaluate F AB more precisely, we are going to study next-leading order terms in F Ω case by case.

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) wherec 1 is a constant for a given J z ∈ [0, 0.50] in our convention. As α = 4, from the general scaling rule in Eq. (D5) we get I i (i = 1, 2, 3) ∝ O(1) and I 4 ∝ 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 = Λ = 10 3 to approximate these pre-factors: [I(x, y) + I(−x, y)] = I(0) + 7.32 c 1 , Here I(0) denotes the on-site contribution and has the value represents the integral over the Brillouin zone: with cos θ(k 1 , k 2 ) = a/ √ a 2 + b 2 and f (k) = a + bi = J x e ik1 + J y e ik2 + J z . Combined with the general expression of F Ω in Eq. (D3), we arrive at where α = I 1 /2 = I(0)/2 + 3.66c 1 , β = 3I 2 /2 = −11.5c 1 , α F = |I 2 |/2 = 3.84 c 1 . (D12) In the special case J x = J y = J z = 1/3, we have As indicated by the inset of Fig. 3 (b), the pre-factor of the L 2 term in F A takes the value α = 0.391 close to I(0)/2. We conclude that the major contribution to F A comes from the on-site interactions.

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 For α = ∞ in Eq. (D5), all of the amplitudes I i (i = 1, 2, 3, 4) ∝ O(1). Yet we can still relate them to the powers of the finite correlation length ξ. The following assumption is taken such that when n < 5, Back to Eq. (D4), we can then replace the summations on the lattice with integrals The pre-factors now read I 1 = I(0) + 7.26c 2 ξ 2 + 4c 2 ξ, Correspondingly, the bipartite fluctuations share the quadratic and linear forms respectively It indicates that in gapped phases, α F increases at the transition towards the intermediate gapless phase. free fermion system, the energy spectrum ξ n (k x ) can be used to calculate the entanglement entropy [65] S(k x ) = − 1 2 Here λ n (k x ) = (e βξn(kx) + 1) −1 denote the eigenvalues of the single-particle correlation function ψ † nkx ψ n kx . It follows a Fermi-Dirac distribution with an inverse temperature β = 1/(k B T ). Fig. 6 (b) shows the numerical entanglement spectrum λ 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 S F = kx S(k x ), through the summation of S(k x ) in the momentum space followed by a finite-size scaling with respect to N x , we are able to obtain the prefactor α S of the linear term in S F . The results are shown in Fig. 4 (d) in the Article. geometry of the lattice Let us take a closer look at a finite honeycomb lattice. The geometric function reads γ k = 1 3 cos (k y / √ 3) + 2 cos (k x /2) cos( √ 3k y /6) . (F5) When h → 0, η → 1, there exists one zero mode k 0 = (0, 0) making f ( r) and g( 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 1/4, r = r ; f 2 ( r − r )/3, r = r and r, r ∈ same sublattice; −g 2 ( r − r )/3, otherwise. (F10) To restore the spin-rotational symmetry at zero magnetic field, we have introduced an extra factor 1/3 to Eq. (F10). Once we take r → ∞, with c 0 = m 0 0.246, c 1 = 3/(8π) 0.119. Numerically, we check in Fig. 7 (a) the finite-size scaling of the single-particle expectation function f ( r) from Eq. (F3). The coefficients read c 0 = 0.245, c 1 = 0.12 with a powerlaw index α = 1.04. They agree well with the asymptotic form.
Since the approximation (−ηγ k ) −1 is still valid for k ∈ Ω (0, ξ), one finds over a long distance g( r) = −f ( r). (F13) The two-spin correlation function then reveal a powerlaw r −1 decay with alternating signs on the same and different sublattices.

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 The bond index i 1 i 2 denotes two sites in the i-th unit cell of the sublattice {1} and {2}. In the modified spinwave theory, S z i1 = 1/2 − a † i a i , S z i,2 = −(1/2 − b † i b i ). Reassembling different terms, we reach 16I(i, j) (F15) To simplify the notation, we have introduced a new set of number operatorŝ First, taking into account a † i a i = b † i b i = −1/2 + f ( 0) = 1/2, the single-particle expectation values become n 1 = n 2 = m 1 = m 2 = 1. (F17) .
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 finitesize effects (the regularization of the zero mode). This is clearly distinct from the pure r −4 scaling in the gapless Kitaev spin liquid phase.