title: Revealing divergent length scales using quantum Fisher information in the Kitaev honeycomb model
authors: James Lambert Erik S. Sørensen
text: 
  
  
  
  
  lambej3@mcmaster.ca
  
  sorensen@mcmaster.ca
  Department of Physics & Astronomy, McMaster University
  1280 Main St. W., Hamilton ON L8S 4M1, Canada
  
    We compute the quantum Fisher information (QFI) associated with two different local
    operators in the ground state of the Kitaev honeycomb model, 
    and find divergent behaviour in the
    second derivatives of these quantities with respect to the driving parameter at
    the quantum phase transition between the gapped and gapless phases for both
    fully anti-ferromagnetic and fully ferromagnetic exchange couplings, thus
    demonstrating that the second derivative a locally defined, experimentally
    accessible, QFI can detect
    topological quantum phase transitions.
    The QFI associated with a local magnetization operator behaves differently
    from that associated with a local bond operator
    depending on whether the critical point is approached from the gapped or
    gapless side. We show how the behaviour of the second derivative of the QFI at
    the critical point can be understood in terms of the diverging length scales
    associated to the two and four point correlators of
    the Majorana degrees of freedom. We present critical exponents associated with the 
    divergences of these length scales. 
  
  
  Revealing divergent length scales using 
  quantum Fisher information in the Kitaev honeycomb model
      Erik S. Sørensen
      December 30, 2023
  =================================================================================================
  
  
  
  
  
  § INTRODUCTION
  
  
  
   §.§ Overview
  
  
  
  
  The quantum Fisher information (QFI), ℱ, arises naturally in quantum metrology
   <cit.>. Given a general
  state ρ(θ) where θ is some parameter, the QFI 
  bounds the precision with which θ may be extracted in any
  M measurements through the
  Cramér-Rao bound, Var_ρ̂(θ_i) ≥ 1
  /√(Mℱ). 
  In other words the QFI quantifies the extent to which a parameterized state ρ(θ) may
  be distinguished from a neighbouring state ρ(θ+θ). By
  quantifying the distinguishability of neighbouring states, 
  the QFI furnishes a natural notion of distance on the Hilbert space, with more
  easily distinguishable states separated by a greater distance. Formally, the QFI
  quantifies the local change in the Bures distance under the aforementioned
  parameterization. <cit.>
  This geometrical interpretation of the QFI expands its scope of
  application to probing the physics of condensed matter phases and phase
  transitions. <cit.> The QFI also exhibits interesting behaviour during a quantum
  quench in spin chain systems <cit.> In fact, a special case of the QFI is already ubiquitous in theoretical studies
  of condensed matter systems. The fidelity
  susceptibility <cit.> (FS) is directly
  proportional to the QFI <cit.> for an appropriate
  parameterization. In particular 
  One often considers parametrizations that have been
  generated unitarily (though this is not the only choice) 
  by a Hermitian operator Ô according to,
  ρ(θ) = e^iθÔρ e^-iθÔ, and we will
  restrict ourselves to this case here. The operator Ô is usually
  expressed as a sum over sub-operators Ô^α,
  
      Ô = ∑_r Ô^α_r_r
  
  A local operator, is one for which all Ô_r^α_r depend on a
  contiguous sublattice that is small relative to the total lattice. One may also
  consider non-local parameterizations, such as the string operators
  considered in Ref. <cit.>. Non-local
  parameterizations reveal remarkable behaviour in topological phases as
  demonstrated in Ref. <cit.>, where a characteristic, 
  super-extensive scaling
  of the non-local QFI is demonstrated in the topologically non-trivial phases of the
  Kitaev wire.  
  
  In a many-body state containing N degrees of freedom, the QFI density
  F=ℱ/N quantifies the degree of multipartite entanglement when the
  state ρ is projected into the eigenbasis of the operator Ô that
  generates the parametrization. For F>m where m|N, we say the state is
  (m+1) partite entangled. <cit.>
  Specifically for the case of pure states, ψ, and unitary parametrizations, the QFI is
  proportional to the variance of the generator, <cit.>
  
      ℱ = 4Var_ψ(Ô)
  
  While one could perform an interesting study looking only at the variances, we
  prefer to work within the context of the QFI because it continues to be well
  defined at finite temperature. This link allows the critical properties of the
  ground state to be inferred from the thermal scaling of the QFI. <cit.> 
  While we do not consider finite temperature behaviour in this study, the
  connection offers a path forward for future work.
  Recently it was shown that the QFI can be detected experimentally in inelastic
  scattering measurements. <cit.> Thus working within the context of
  the a locally defined QFI also allows for connection with experiment. This
  contributes to a growing body of research on experimental approaches to extract
  multipartite
  entanglement. <cit.>
  We emphasize that we do not detect genuine multipartite entanglement in this
  study, but mention the connection for completeness. 
  
  Given that the QFI is defined at finite temperature and that the zero temperature
  QFI is proportional to the variance, one might ask, what is the generalization of the notion of
  variance to the finite temperature case? By imagining that the variance
  of an observable contains a
  quantum contribution and a thermal contribution,  
  quantum variances (QV) may be defined which are proportional to an upper and
  lower
  bound of the QFI. <cit.> The QFI at zero temperature can be
  viewed as the zero temperature limit of the quantum contribution to the
  variance (the thermal contribution being zero at zero temperature).
  
  
  
  The QFI has now been studied in a wide range of
  models. <cit.> Of
  particular interest for our purposes is the work done on the Kitaev wire
  in Ref. <cit.>, where the first derivative of the QFI associated with a
  local generator was shown to exhibit a divergence at the topological phase
  transition of that model, and where the topologically non-trivial phase
  exhibits super-extensive scaling of the QFI associated with a non-local
  generator.  
  
  Quantum spin liquids (QSL) are characterized by a lack of any form of long range
  magnetic order down to zero temperature. <cit.> Such phases are
  thought to exhibit instead subtle forms of quantum ordering, along with
  topologically non-trivial anyonic excitations. <cit.> In so far
  as these phases are characterized by a lack of order, their detection in
  experiment presents a substantial challenge. In this work we examine the
  behaviour of the QFI in the Kitaev honeycomb model
  (KHM) <cit.>, which presents two spin liquid phases (one
  gapped and one gapless), induced by exchange coupling anisotropy. 
  
  The KHM has been studied from an information theoretic perspective before, with
  studies examining the Jensen-Shannon
  divergence <cit.>, and the mutual
  information. <cit.> Of particular interest is the Fidelity
  susceptibility, which was studied in Ref. <cit.> and the study of the
  Bures distance in. <cit.> In the case of an n parameter
  estimation scenario (or an n dimensional unitary parametrization), the Bures
  distance is locally equivalent to the QFI Matrix which is a Riemannian metric on
  the Hilbert space. <cit.> The Fidelity susceptibility is recovered by examining the
  particular parametrization of the Hilbert space corresponding to the driving
  operator of the phase transition. The physics of the KHM have also been studied
  using SU(2) parton approaches <cit.>. Details of the
  dynamical response of the model in the presence of magnetic fields may be found
  in Ref. <cit.>
  
  For the remainder of this section we introduce the KHM and explain its key
  features. In Sec. (<ref>) we discuss the 
  relationship between the scaling of the second
  derivative of the QFI density (hereafter called the QFI susceptibility), ∂_u^2F, 
  where u drives the phase transition,
  and the correlation functions of the generator.
  In Sec. (<ref> we analyze the behaviour of these quantities
  for the magnetization operator, ∑_j S_j^α, (where j represents
  both a unit cell position and sublattice index), in
  Sec. (<ref>), and the bond, ∑_r S_r,A^α
  S_r,B^α, where r indicates a unit cell, in
  Sec. (<ref>). Finally, we conclude our discussion in
  Sec. (<ref>), where we discuss the relevance of this work to
  studies of  the geometric phase. 
  
  
  
  
  
   §.§ Kitaev Honeycomb Model
  
    
    The Kitaev honeycomb model (KHM) is given by,
    
      H = ∑_⟨ j,k⟩ K^γ_j,k
          S^γ_j,k_jS^γ_j,k_k
  
    where the sum is over nearest neighbour bonds and γ∈{x,y,z} denotes a
    bond-dependent Ising exchange. 
    If the exchange couplings are sufficiently
    isotropic (|K^γ| ≤|K^α| + |K^β|, for all
    choices of 
    α,β,γ∈{x,y,z}), the spectrum is gapless. In the regime
    where one exchange coupling is dominant (the opposite inequality), the model
    is gapped. This phase transition between two topologically different spin
    liquid phases presents no local order parameter. It is instead associated with
    a subtle kind of symmetry breaking to do with the structure of the gauge
    fields themselves. <cit.> 
    On the gapped side of the transition, the model is mapped
    onto the lattice gauge Ising model <cit.>, 
    with alternative rows of hexagon plaquettes
    becoming associated with one of the two excitations in that model
    (conventionally called e and m  excitations). In both phases, the spin-spin 
    correlation functions are identically zero beyond nearest neighbour. 
    The model also possesses an extensive number of conserved charges defined by
    the plaquette operators.
  
    Remarkably, the KHM is analytically solvable. <cit.> 
    By mapping each spin operator
    into the space of four Majorana fermions, {c,b^x,b^y,b^z} via
    
      S^γ_j = 1/4ic_jb^γ_j
  
    an extensive number of conserved charges can be constructed, given by u_j,k =
    ib_j^γ_j,k b_k^γ_j,k. These operators take eigenvalues ±1.   
    Using the above mapping, the KHM becomes,
    
      H = i/4∑_jk K^γ_j,ku_jk c_jc_k
  
    Since the u_j,k commute with the Hamiltonian, we may fix a particular
    configuration of eigenvalues on each bond, and the problem is reduced to free
    Majoranas hopping in the gauge fields. 
    The lowest energy configuration will be the flux free configuration, as
    follows from
    Lieb's theorem. <cit.> We therefore choose to work in the configuration where all
    u_i,j have eigenvalue +1 (hereafter referred to as the standard gauge).
    Once we fix a gauge configuration, the model is a simple hopping Hamiltonian,
    which may be diagonalized by Fourier transforming and then performing a
    Bogoliubov rotation, where the mixing angle is defined implicitly via,
    
      tan(2θ_q) = ϵ_q/Δ_q,
  
    where, 
    
    
      ϵ_q    = K_αcos(q_x) + K_βcos(q_y) + K_γ
      Δ_q    = K_αsin(q_y) + K_βsin(q_y)
  
     
    Where α,β,γ∈{x,y,z} depending on the choice of which bond
    acts as the unit cell. Here q_x=𝐚·𝐪 and
    q_y=𝐚·𝐪 where 𝐚_1 and 𝐚_2 are
    any choice of translation vectors on the principle lattice and
    q=n/L_x𝐛_1 + m/L_y𝐛_2, with L_x and
    L_y the side length of the lattice, is a general vector in the reciprocal
    space. 
  
    We note that the true, physical ground state, must be the
    symmetrized product over all physically equivalent choices of the gauge
    fields (i.e. all choice of the gauge fields resulting in zero flux).
    Following arguments described in Ref. <cit.>, 
    the operators we consider are not dependent on projection into the physical
    subspace at large system sizes.
    Details of the solution to the Kitaev model are provided
    in <ref>. 
  
  
  
  
  
  § QFI SUSCEPTIBILITY AND DIVERGING LENGTH
  SCALES
  
  
  In order to interpret the divergences at the critical point, 
  consider a generator, Ô=∑_rÔ_r  where the generators are 
  given by a sum over local products of spins,
  𝒮̂_r = ∏_j∈Λ_rS_r+ℓ_j^α_j, where Λ_r is
  some local, contiguous sublattice. The associated QFI density in a pure state is given by, 
  
      F{Ô}   =
        1/N∑_r_1,r_2𝒮_r_1𝒮_r_2
        -
        𝒮_r_1𝒮_r_2
  
  Through Kitaev's mapping, we may decompose our spin blocks into a component
  operating on the flux sector, ℬ_r=∏_j∈Λ_r
  b_r+ℓ_j^α_j, 
  and a component operating on the matter sector,
  𝒞_r=∏_j∈Λ_rc_r+ℓ_j,
  
      F{Ô} =
        1/N∑_r_1,r_2ℬ_r_1ℬ_r_2𝒞_r_1𝒞_r_2
        - 
        ℬ_r_1ℬ_r_2𝒞_r_1𝒞_r_2
  
  Now there are three possible values for the flux sector expectation values. If
  ℬ_r is diagonal in the gauge sector, then the contribution from the
  gauge sector factorizes and gives an overall prefactor of ± 1. If
  ℬ_r is strictly off-diagonal, but
  ℬ_r_1ℬ_r_2 has diagonal entries, then the situation is
  the same. Finally, it may be the case that ℬ_r_1ℬ_r_2
  has non-zero diagonal elements only for certain separations. Regardless of which
  scenario is realized, the contribution to the QFI from the flux sector will be
  independent of u, since the gauge fields commute at all points in the phase
  diagram. The QFI is then given by a sum over the correlation functions  in the
  matter sector Majorana fermions with some prefactor (which might be ±1 or
  0 as a function of the separation), determined by the situation
  above. We adopt the standard ansatz for the matter sector correlations 
  
      𝒞_r_1𝒞_r_2 = Ξ(r,u)r^-ae^-r/ξ(u)
  
  where a is determined by the phase and does not depend explicitly on u, and
  ξ(u) is a length scale associated with correlations between the Majorana
  operators which depends on the position in the phase diagram.
  
  Taking the assumption that, near the critical point,
  
      ξ(u)∼|u-u_c|^-ν
  
  one can show that the second derivative of the QFI density must diverge at the
  critical point like, 
  
      ∂_u^2F{Ô}∼|u-u_c|^ν-2
  
  Therefor Δ_Ô = ν-2. The QFI and, by extension, the QFI
  susceptibility, are, in principle, experimentally accessible probes. In
  particular for the case of local generator such as the total magnetization
  operator. The QFI can there be used to extract experimentally the scaling of the
  correlation length associated with the matter sector of the KHM. This analysis
  is similar to the analysis performed in the supplementary materials of
   <cit.>, where the authors examined the effects of coarse-graining
  transformation on the QFI to arrive a scaling hypothesis for the near field and
  finite temperature regimes. 
  
  In practice,
  experimentally relevant models will not be amenable to the above treatment, as
  the KHM acquires additional terms in real materials that break the
  integrability of Kitaev's  original solution (e.g. Heisenberg terms and
  symmetric off-diagonal terms). <cit.> In these cases the analysis may
  instead be applied to the correlation length of the spin degrees of freedom
  directly, and divergences in the QFI susceptibility may still be
  linked to the critical exponent for the divergence of a correlation length. 
  
  
  
  
  § QFI SUSCEPTIBILITY FOR MAGNETIZATION AND BOND
  OPERATORS
  
  
  Motivated by the results of  <cit.>, we compute the QFI
  associated with two local operators and examine the second derivatives of those
  operators with respect to the driving parameter of the phases transition. In
  both cases the second derivative of the QFI is found to diverge. We term this
  the QFI susceptibility.
  
  Throughout this section we consider a path through the space of exchange
  couplings parameterized by u, 
  
    
      K^z    = 1 + u, 
      
          K^x    = 1, 
      
          K^y    = 1
  
  
  For this parametrization, u_c=1 represents the critical point between the
  gapless phase (u<1)  and the gapped phase (u>1). The ground state of the
  KHM is a function of u and is hereafter denoted as ψ_0(u).
  We consider the case of fully ferromagnetic and fully anti-ferromagnetic
  exchange couplings for both parametrizations (in which case K^γ→
  -K^γ). Unless otherwise noted, calculations are carried out for L_x=L_y=10^4 with periodic
  boundary conditions in a rhombic geometry. 
  
  
  
  
    
    
  
  
   §.§ Magnetization Operator
  
  
  First we examine the QFI in the Kitaev honeycomb as generated by the
  magnetization operator,
  
      Ô_M^α = ∑_r(S_r,A^α + S_r,B^α).
  
  here α∈{x,y,z}, r denotes a unit cell in the two site
  basis, and A,B denotes the sublattice. The corresponding QFI is given by, 
  
      F_M,α(u) ≡ F{Ô_M^α; ψ_0(u)}   = 4Var_ψ_0(O_M^α) 
         =1 - 1/N∑_qcos(2θ_q)
  
  with ψ_0(u) defined at the start of this section.
  
  Fig. <ref> shows this quantity plotted along 
  the path defined by <ref> for the fully antiferromagnetic (a)
  and ferromagnetic (b) cases
  respectively. 
  In the AFM case the ground state possesses F_M,α > 1 for each spin
  component. This is indicative of at least bipartite entanglement. In a pure
  state any non-zero QFI is indicative of the presence of quantum correlations.
  Nonetheless, the QFI in the fully ferromagnetic case is insufficient to witness
  even bipartite entanglement, indicating that quantum correlations are reduced
  for the FM coupling. 
  
  
  The absolute values of the derivatives of F_M,α will be the same in the 
  AFM and FM cases.
  This can be seen by considering the fact that the functional dependence of
  F_M,α on the driving parameter u enters through the nearest neighbour
  correlation functions, which are the same in both cases up to a negative sign.
  
  The first and second derivatives of the QFI are given in
  Fig. <ref>. We observe that the QFI susceptibility associated
  with the magnetization operator exhibits a power law divergence when
  approaching the critical point from the gapped side. When approaching the
  critical point from the gapless side the transition appears first order. The
  behaviour of the transition from the gapped side can be understood in light of the analysis in
  Sec. <ref>. 
  Using the scaling hypothesis,
  
      ∂_u^2 F_M,α (u) ∼|u-u_c|^Δ_M,α
  
  we extract the following critical exponents for the second derivative of the QFI
  for each spin component of the magnetization operator,
  
    
      Δ_M,x = Δ_M,y   ≈ -0.52(1) 
      Δ_M,z   ≈ -0.50(1)
  
  
  and can be seen in Fig. <ref> over a region from 10^-1.7 >
  u >10^-3.0. At this point finite, finite size effects enter, and the scaling
  ansatz is no longer valid. This leads to plateaus in the QFI susceptibility
  which occur closer to the critical point for larger system sizes as seen from
  the data in Fig. (<ref>) and
  Fig. (<ref>).
  
  
  
  
   §.§ Bond Correlation Operator
  
  
  We now turn to the QFI as parameterized by the bond correlation operator,
  
      Ô_B^α = ∑_r S_r,A^α S_r,B^α,
  
  with the corresponding QFI density given by, 
  
      F_B,α = 1/2N∑_q sin^2(2θ_q)
  
  In this case the AFM and FM cases are identical. 
  We repeat the same analysis as for the magnetization operator as shown in
  Fig. <ref> and  <ref>. The QFI associated with the
  bond operator along the x and y components converges to a constant value
  immediately following the phase transition, while the QFI associated with the
  z component bond operator falls towards zero. This behaviour can be understood
  by the fact that the Hamiltonian in the gapped phase is dominated by the Ising
  exchange on the z bonds. Consequently, the commutator between the Hamiltonian
  and the bond operator approaches zero in the limit of u→∞. 
  
  Assuming the same scaling ansatz as for the second derivative F_M,α we
  find a crossover between two scaling regimes. The first regime is given by the critical
  exponents, 
  
    
      Δ_B,x^(I)=Δ_B,y^(I)   ≈ -0.52(1) 
      Δ_B,z^(I)   ≈ -0.58(1)
  
  
  which appears valid on the interval 10^-1.9 > u >
  10^-2.8 and a second regime
  characterized by the exponents,
  
    
      Δ_B,x^(II)=Δ_B,y^(II)   ≈ -0.62(1) 
      Δ_B,z^(II)   ≈ -0.65(1)
  
  
  which appears to be valid on the interval 10^-3.0> u>
  10^-3.3. While the magnetization operator exhibits a divergence when approaching the
  critical  point from the gapped phase, the bond operator exhibits a divergence
  approaching the critical point from the gapless phase. 
  
  The oscillatory behaviour on the gapless side of the transitions for both
  QFI's is related to divergences in the QFI susceptibility due to points where
  the denominator of the integrand goes to zero. In the gapped phase these points
  are necessarily absent. 
  
  
  
   §.§ Diverging length scales
  
  
  Using the results of the previous section we can determine the scaling of the
  divergence in the correlation length for c_r_1,Ac_r_2,B 
  (using the divergence in the
  magnetization operator), and for c_r_1,Ac_r_1,Bc_r_2,Ac_r_2,B.
  
  In light of Eq. (<ref>), we can now understand that the QFI
  susceptibility associated with the magnetization operator diverges from the
  gapped side due specifically to the divergence in that correlation function of
  the matter sector Majorana's. On the gapless side of the transition, the
  correlation function for the matter sector Majorana's is critical, and
  consequently the second derivative of Eq. (<ref>) is given
  specifically by Ξ(r,u) and contains no divergence.
  
  Using
  Eq. (<ref>), we may extract the scaling exponents for the
  correlation length of the matter sector correlation functions in the x and
  y, and z channels for ∂_u^2F_M,α,
  
    
      ν_M,x=ν_M,y   ≈ 1.48(1) 
      ν_M,z   ≈ 1.50(1)
  
  
  and for the two scaling regimes of ∂_u^2F_B,α
  for u<u_c. The first given by,
  
    
      ν_B,x^(I) = ν_B,y^(I)   ≈ 1.48(1) 
      ν_B,z^(I)   ≈ 1.42(1),
  
  
  and the second by,
  
    
      ν_B,x^(II) = ν_B,y^(II)   ≈ 1.38(1) 
      ν_B,z^(II)   ≈ 1.35(1),
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § CONCLUSION
  
  
  We have examined the QFI for the bond and magnetization operators in both the
  gapped and gapless phases of the KHM and at the transition between these two
  phases for fully ferromagnetic and fully anti-ferromagnetic couplings. The
  second derivative of the QFI with respect to the magnetization operator is shown
  to diverge when approaching the phase transition from the gapped side like at a
  second order transition, while the QFI susceptibility approaching the critical
  point from the gapless side appears first order. Conversely we find that the QFI 
  susceptibility associated with the bond operator
  diverges like a second order transition when approaching the critical point from
  the gapless side, and like a first order transition when approaching the
  critical point from the gapped side. 
  
  In both cases, the divergences in the QFI susceptibility can be associated to
  diverging length scales in the two point correlators of the local generators of
  the QFI. For the particular case of the KHM, these divergences can be linked to
  diverging length scales in the matter sector Majorana's, even when the physical
  spin-spin correlation functions are truncated (as in the case of the two point
  correlation function). The implication is that presence of the topological phase
  transition between the gapped and gapless phases may be detected experimentally
  at low temperatures.  
  
  The has been related work examining the Geometric phase associated with a twist operator
  acting on both sites. <cit.> We note that the critical
  exponents presented in Eq. <ref> for scaling
  regime (II) of the bond operators are within the margin of error of those in
  Ref. <cit.>. Geometrically, the QFI that we compute with
  respect to the bond operator is the diagonal component of the quantum geometric
  tensor <cit.>. The imaginary component of this tensor
  corresponds to the Berry curvature, while the real component corresponds to the
  notion of distance induced by the distinguishability of states. In
  Ref. <cit.>, the connection between these two components
  of metric is discussed. The implication is that the geometry detected by the
  Berry phase is intimately related to the geometry of distinguishability, opening
  the prospect of experimentally measuring the Berry phase in condensed matter
  systems. Extracting the full quantum metric tensor has recently been achieved in
  cold atom systems. <cit.>
  
  The QFI associated with the magnetization operator in the fully anti-ferromagnetic
  phase is shown to be greater than for the fully ferromagnetic phase, as one
  would expect from the tendency of the anti-ferromagnetic coupling to produce
  spin singlets on the bonds. In the gapped phase defined by large K^γ
  coupling, the QFI associated with the
  bond operator converges to a constant value for the QFI generated by the
  transverse spin components (specifically the x and y components in our
  analysis).
  
  It is shown in  <cit.> that the QFI is proportional to the
  FS if the operator parameterizing the QFI is the same as the
  operator that generates the change in parameter for the ground state. 
  This implies that the QFI for the bond operator from the gapless side of the
  transition is proportional the FS calculated
  in Ref. <cit.>, however we do not find this to be the case.
  This may be due to the fact that the ground state used in
  Ref. <cit.> differs from that used by Kitaev in
  Ref. <cit.>, which is the one we employ here. 
  Understanding the details  of the connection between the QFI and the fidelity
  susceptibility warrants further investigation. We also note that the definition of the QFI may not be unique when a
  Hamiltonian posseses a degenerate ground state manifold. In particular one can
  imagine a situation where the generator of the QFI lifts the degeneracy of this
  manifold, affecting the results. While in our calculation this ambiguity is not present
  in the gapped phase, it may affect the results in the gapless phase. 
  
  Future research is warranted to examine the behaviour of the QFI at finite
  temperatures around the critical point, where the ground state scaling will be
  modified by finite temperature effects. The connection between the finite
  temperature scaling and the length scale of the Majorana fermions in this case may
  offer insight into the details of candidate Kitaev spin liquid phases in
  materials where the pure Kitaev Hamiltonian is modifed by material relevant
  terms. <cit.>
  
  
  
  
  
  § SCALING BEHAVIOUR OF QFI SUSCEPTIBILITY
  
  
   Let's work specifically on the case of pure states and unitary QFI. 
  The generator of the QFI is most generally given by,
  
      Ô = ∑_rÔ_r
  
  where Ô_r is an operator associated with the site located at r. We
  assume that r is contiguous and local, that is it encompasses a finite number of degrees of
  freedom all lying within a distance ℓ from the site r. We consider a
  state ψ  that depends on some parameter u that drives a phase transition
  at a value u_c = 1. 
  
      f{Ô,ψ(u)} =1/N∑_r_1,r_2⟨Ô_r_1Ô_r_2⟩_ψ -
        ⟨Ô_r_1⟩_ψ⟨Ô_r_2⟩_ψ
        = 1/N∑_r_1,r_2C_r_1,r_2(u)
  
  Let's assume that the model is translation invariant and define r|r_1-r_2|. In general, we may assume that the connected correlation functions can be
  fit to the following form,
  
      C_r(u) = Ξ(r,u)r^-ae^-r/ξ(u)
  
  where a depends on the phase (i.e. is assumed independent of the driving
  parameter), and ξ is the correlation length, taken to be a function of the
  parameter u (we hereafter drop the explicit dependence). The function
  Ξ(r,u) is assumed to be a smooth function of r and of the parameter u
  within a particular phase (though not necessarily smooth at the phase boundary). 
  The divergence in the
  second derivative of the QFI must emerge from a divergence in the two point
  correlation functions. We therefore consider the second derivative of
  Eq. (<ref>)
  
      ∂_u^2C_r(u)    = 
        ∂_u
        (
          Ξ(r,u)r^1-aξ^-2∂_uξ e^-r/ξ 
          + ∂_uΞ(r,u)r^-ae^-r/ξ) 
         = 
        -2r^1-aΞ(r,u)ξ^-3(∂_uξ)^2 e^-r/ξ
        + r^1-aΞ(r,u)ξ^-2∂_u^2ξ e^-r/ξ 
        + r^2-aΞ(r,u)ξ^-4(∂_uξ)^2
        e^-r/ξ
         +∂_u^2Ξ(r,u)r^-ae^-r/ξ
        + ∂_uΞ(r,u)r^1-aξ^-2∂_uξ e^-r/ξ
  
  Naively, the correlation length is expected to diverge at the critical point.
  let ũ=|u-u_c| be the distance from the critical point. Then the
  correlation length goes as,
  
      ξ∼ũ^-ν
  
  This ansatz may be used to infer the scaling relations for the derivatives of
  the correlation length,
  
  
      ∂_uξ   ∼ -νũ^-(ν+1)
      ∂_u^2ξ   ∼ν(ν+1) ũ^-(ν+2)
  
  
  Substituting this into Eq. (<ref>), gives,
  
      ∂_u^2C_r(u)    = 
        - 2Ξ(r,u)r^1-aũ^3νν^2ũ^-2(ν+1)e^-r/ξ
        + Ξ(r,u)r^1-aũ^2νv(v+1)ũ^-(ν+1) e^-r/ξ 
        + Ξ(r,u)r^2-aũ^4νν^2ũ^-2(ν+1)e^-r/ξ
         +∂_u^2Ξ(r,u)r^-ae^-r/ξ
        + ∂_uΞ(r,u)r^1-aũ^2ν(-ν)ũ^-(ν+1)
        e^-r/ξ
         = e^-r/ξ(
          -2Ξ(r,u)r^1-aũ^ν - 2ν^2 
          + Ξ(r,u)r^1-aũ^ν - 2ν(ν+1)
          + Ξ(r,u)r^2-aũ^2ν-2.
         +.∂_u^2Ξ(r,u)r^-a 
          -ν∂_uΞ(r,u)r^1-aũ^ν-2)
  
  We can now pull out the divergence associated with the proximity to the critical
  point. 
  
      ∂_u^2 C_r(u)    = ũ^ν-2e^-r/ξ(
          -2Ξ(r,u)r^1-aν^2 
          + Ξ(r,u)r^1-aν(ν+1)
          + Ξ(r,u)r^2-aũ^ν.
         +.ũ^2-ν∂_u^2Ξ(r,u)r^-a 
          -ν∂_uΞ(r,u)r^1-a)
  
  The scaling behaviour of the QFI susceptibility is thus given by,
  
      ∂_u^2f{Ô,ψ}=ũ^ν-2ζ(r,u)
  
  We define Δ_Ô = ν- 2 as the scaling of the QFI with proximity
  to the critical point. The y-intercept on the log-log plot will be given by the
  non-universal function ζ(r,u).
  
  
  
  § CALCULATING VARIANCES
  
  
  
  
   §.§ Solution of the Kitaev Model
  
  
  We adopt the approach of <cit.>, where the Majorana degrees of
  freedom are recombined into Dirac fermions, with three bond fermions,
  
    
      b_r,A^γ   = 1/2(β_r^γ +
          (β_r^γ)^†) 
      
          b_r,B^γ   = 1/2i(β_r^γ -
          (β_r^γ)^†)
  
  
  and one matter fermion
  
  
      c_r,A   = 1/2(f_r + f_r^†) 
      
          c_r,B   = 1/2(f_r - f_r^†)
  
  
  The bond fermions are not present in the Hamiltonian, since we simply replace
  the bond operators with the eigenvalues of the standard gauge configuration
  (u_j,k=1). The resulting Hamiltonian is quadratic in the matter fermions and
  translation invariant. It can be diagonalized first by mapping each matter
  fermion to momentum space, f_r = 1/√(N)∑_q e^iq· r f_q,
  and then applying the Bogoliubov rotation, f_q = cos(θ_q)a_q +
  isin(θ_q)a_-q^†, where θ_q is defined by,
  
      tan(2θ_q) = K^xcos(q_x) + K^ycos(q_y) + K^z/K^xsin(q_x) +
        K^ysin(q_y).
  
  
  
  
  
  
  
   §.§ Magnetization Operator
  
  Begin with, 
  
      = ∑_rα_r
  
  The variance is given generally by, 
  
      Ô = O^2 - O^2,
  
  which, for the magnetization operator gives,
  
      = ∑_r_1,r_2α_r_1α_r_2
  
  Using translation invariance and converting the Majorana representation, this
  expression can be given as,
  
      Ô^α_Mag   = N ∑_rα_0α_r
  
  In the Kitaev model, the two  point correlator is zero for all values of r
  except nearest neighbours. Thus the sum above can be reduced to, 
  
      Ô^α_Mag = N ( 1/4 +
        α_0,Aα_0,B)
  
  Thus we only need to calculate the nearest neighbour correlation function, 
  
      α_0,Aα_0,B   = 1/4σ^α_0,Aσ^α_0,B
         = 1/4(ib_0,A^α c_0,A)(ib_0,B^α c_0,B)
         = 1/4⟨ℱ|b_0,A^α b_0,B^α|ℱ⟩⟨ℳ|c_0,Ac_0,B|ℳ⟩
         = 1/4⟨ℱ|(-i)(2n̂^β^α_0 - 1)|ℱ⟩⟨ℳ|(-i) (2n̂^f_0 - 1)|ℱ⟩
         = -1/4(2n̂^f_0 - 1) 
         = -1/4(21/N∑_q_1,q_2f^†_q_1f_q_2 -
        1)
         = -1/4(2/N∑_qsin^2(θ_q) - 1)
  
  The QFI density is four times the variance divided by the system size. Thus,
  
      f{} = 1 +
        σ_0,A^ασ_0,B^α
  
  
  
  
   §.§ Bond Operator
  
  
      Var()    = ()^2 - ^2 
         = ∑_r_1,r_2α_r_1 Aα_r_1 Bα_r_2 Aα_r_2 A
          - (∑_r α_r_1 Aα_r_1 B)^2 
         = 1/16(∑_r_1,r_2b_r_1 A^α b_r_1 B^α b_r_2 A^α b_r_2 B^αc_r_1 Ac_r_1 Bc_r_2 Ac_r_2 B - (∑_r
          -b_r A^α b_r B^αc_r Ac_r B)^2
          )
  
  We can compute the flux sector expectation values easily, 
  
      b_r_1 A^α b_r_1 B^α b_r_2 A^α b_r_2 B^α   = (-i)^2(β_r_1 + β_r_1^†)(β_r_1-β_r_1^†)
                 (β_r_2 + β_r_2^†)(β_r_2-β_r_2^†)
         = -(2n_r_1 - 1)(2n_r_2 - 1)
         = -1
  
  Where the last line
  follows from the fact that the ground state in the standard flux configuration
  is defined by u_r = 2n_r - 1 = 1. Similarly we find,
  
      b_rA^α b_rB^α   = (-i)(2n_r -1)
         = (-i)
  
  The variance is therefore,
  
      Var()    = 1/16(∑_r_1,r_2
        (-1)c_r_1Ac_r_1Bc_r_2Ac_r_2B
        -( ∑_r(i)c_rAc_rB)^2  
        )
        
         = 1/16(
        -∑_r_1,r_2c_r_1 Ac_r_1 Bc_r_2 Ac_r_2 B
        + 
        ∑_r_1,r_2c_r_1 A c_r_1 Bc_r_2 A c_r_2 B) 
         = 1/16(
         ∑_r_1,r_2 
         -c_r_1 Ac_r_1 Bc_r_2 Ac_r_2 B 
         +c_r_1 Ac_r_2 Ac_r_1 Bc_r_2 B
         -c_r_1 Ac_r_2 Bc_r_1 Bc_r_2 A. 
         +.c_r_1 Ac_r_1 Bc_r_2 Ac_r_2 B) 
         = 1/16(
          ∑_r_1,r_2c_r_1 Ac_r_2 Ac_r_1 Bc_r_2 B
         -c_r_1 Ac_r_2 Bc_r_1 Bc_r_2 A)
  
  We now need only evaluate the two point correlators above. For the first term we
  have,
  
      c_r_1 A c_r_2 A   = 
        1/N∑_q_1,q_2 e^iq_1r_1 e^iq_2r_2(f_q_1 + f_q_1^†)(f_q_2 + f_q_2^†)
         = 1/N∑_q_1,q_2e^iq_1r_1 e^iq_2r_2f_q_1f_q_2+f_q_1f_q_2^† + f_q_1^†
                  f_q_2 + f_q_1^† f_q_2^†
         =1/N∑_q_1,q_2e^iq_1r_1 e^iq_2r_2f_q_1f_q_2 + f_q_1^† f_q_2^† +
                  δ_q_1,q_2
  
  We can see that,
  
  
      f_q_1f_q_2   = icos(θ_q_1)sin(θ_q_2)a_q_1
        a_-q_2^†
         = icos(θ_q_1)sin(θ_q_2) δ_q_1,-q_2
      f_q_1^† f_q_2^†   =
        (-i)cos(θ_q_1)sin(θ_q_2) a_-q_1a_q_2^†
         = (-i) cos(θ_q_1)sin(θ_q_2)δ_-q_1,q_2
  
  
  Consequently, 
  
      c_r_1 A c_r_2 A   = 
        1/N∑_q_1,q_2 e^iq_1r_1e^iq_2r_2δ_q_1,q_2
         = 1/N∑_q e^iq(r_1-r_2)
         = δ_r_1,r_2
  
  Similarly we may show that, 
  
      c_r_1 B c_r_2 B = δ_r_1,r_2
  
  For the second term we begin with, 
  
      c_r_1Ac_r_2B   =
        1/N∑_q_1,q_2e^iq_1r_1e^iq_2r_2
        (-i)(f_q_1+f_q_1^†)(f_q_2-f_q_2^†
         = 1/N∑_q_1,q_2e^iq_1r_1e^iq_2r_2(-i)
        f_q_1f_q_2-f_q_1f_q_2^† + f_q_1^† f_q_2
        -f_q_1^† f_q_2^†
         = (-i)/N∑_q_1,q_2 e^iq_1r_1e^iq_2r_2(f_q_1f_q_2 - f_q_1^† f_q_2^† 
                      +f_q_1^† f_q_2 - f_q_1f_q_2^†)
                      
         = (-i)/N∑_q_1,q_2 e^iq_1r_1e^iq_2r_2(f_q_1f_q_2-f_q_1^† f_q_2^† 
             +2f_q_1^† f_q_2 - 1) 
         = (-i)/N∑_q_1,q_2e^iq_1r_1e^iq_2r_2(icos(θ_q_1)sin(θ_q_2)(δ_q_1,-q_2 +
             δ_-q_1,q_2) +
             2sin(θ_q_1)sin(θ_q_2)δ_q_1,q_2 - 1)
         =(-i)/N∑_q 2ie^-iq(r_1-r_2)cos(θ_q)sin(θ_q) 
             +(-i)/N∑_q (2sin^2(θ_q) - 1)
         = 1/N∑_q e^-iq(r_1-r_2)2cos(θ_q)sin(θ_q)
             +i/N∑_qe^iq(r_1+r_2)cos(2θ_q) 
         = 1/N∑_q e^-iq(r_1-r_2)sin(2θ_q) 
                + i/N∑_qe^iq(r_1+r_2)cos(2θ_q)
  
  The final two point correlator is given by,
  
      c_r_1Bc_r_2A   = 1/N∑_q_1,q_2e^iq_1r_1e^iq_2r_2 
        (-i)(f_q_1-f_q_1^†)(f_q_2+f_q_2^†)
         =(-i)/N∑_q_1,q_2e^iq_1r_1e^iq_2r_2(
          f_q_1f_q_2 - f_q_1^† f_q_2^†
          -2f_q_1^† f_q_2 - 1)
         =(-i)/N∑_q_1,q_2e^iq_1r_1e^iq_2r_2(icos(θ_q_1)sin(θ_q_2)(δ_q_1,-q_2 +
        δ_-q_1,q_2) - (2sin(θ_q_1)sin(θ_q_2)δ_q_1,q_2 -
        1)
        ) 
         = (-i)/N∑_q 2ie^-iq(r_1-r_2)cos(θ_q)sin(θ_q)
        -(-i)/N∑_q e^iq(r_1+r_2)(2sin^2(θ_q) - 1) 
         = 1/N∑_qe^-iq(r_1-r_2)sin(2θ_q) -
        i/N∑_qe^iq(r_1+r_2)cos(2θ_q)
  
  The final term is therefore a product of differences,
  
      c_r_1Bc_r_2Ac_r_1Ac_r_2B   =
        1/N^2∑_q_1,q_2(e^-i(q_1+q_2)(r_1-r_2)sin(2θ_q_1)sin(2θ_q_2)
        +e^i(q_1+q_2)(r_1+r_2)cos(2θ_q_1)cos(2θ_q_2)) 
         = 1/N^2∑_q_1,q_2(e^ir_1(q_1+q_2)
        e^-ir_2(q_1+q_2)sin(2θ_q_1)sin(2θ_q_2)
        + e^ir_1(q_1+q_2)e^ir_2(q_1+q_2)cos(2θ_q_1)cos(2θ_q_2)
       )
  
  Under the summation, we may extract the delta functions,
  
      ∑_r_1,r_2c_r_1Bc_r_2Ac_r_1Ac_r_2B   = 1/N∑_r_2∑_q_1,q_2δ_q_1,-q_2e^-ir_2(q_1+q_2)sin(2θ_q_1)sin(2θ_q_2) + δ_q_1,-q_2
        e^ir_2(q_1+q_2)cos(2θ_q_1)cos(2θ_q_2)  
         = ∑_q cos^2(2θ_q) - sin^2(2θ_q)
  
  We are now ready to return to our original expression for the variance  which
  reads,
  
      Var()    = 1/16(N - ∑_q (cos^2(2θ_q) -
        sin^2(2θ_q)))  
         = 1/16∑_q(1-cos^2(2θ_q) +
        sin^2(2θ_q)) 
         = 1/8∑_qsin^2(2θ_q)
  
  Thus the associated QFI density is,
  
      ℱ() = 4Var()/N = 
        1/21/N∑_qsin^2(θ_q)
  
  passing to the continuum limit we find, 
  
      ℱ() = 1/2∫_BZsin^2(2θ_q) ^2q
  
  
  
  
  § DERIVATIVES OF THE QFI SUSCEPTIBILITY MOMENTUM DENSITY
  
  
  In order to gain more insight into the divergence in F_M,α, we define
  the QFI momentum density for the magnetization operator, f_M,α
  via,
  
      F_M,α = ∑_q(1/N - cos(2θ_q)) = ∑_q
        f_M,α(q;u)
  
  We may explicitly evaluate the first and second derivative of this quantity for
  each spin component, giving, 
  
    
      ∂_uf_M,x   = 
          Δ_x^2cos(q_y) - ε_xΔ_xsin(q_y)
          /(ε_x^2 + Δ_x^2)^3/2
      ∂_uf_M,y   = 
          Δ_y^2cos(q_x) - ε_yΔ_ysin(q_x)
          /(ε_y^2 + Δ_y^2)^3/2
      ∂_uf_M,z   =
          Δ_z^2/(ε_z^2 + Δ_z^2)^3/2
  
  
  for the first derivatives, and
  
    
      ∂_u^2 f_M,x   =
        
          -3ε_xΔ_x^2cos^2(q_y)
          +Δ_x(3(ε_x^2-Δ_x^2) +
          ε_x+Δ_x)sin(q_y)cos(q_y)
          +ε_x(2Δ_x^2-ε_x^2)sin^2(q_y)
        /(ε_x^2+Δ_x^2)^5/2
      ∂_u^2 f_M,y   =
        
          -3ε_yΔ_y^2cos^2(q_x)
          +Δ_y(3(ε_y^2-Δ_y^2) +
          ε_y+Δ_y)sin(q_x)cos(q_x)
          +ε_y(2Δ_y^2-ε_y^2)sin^2(q_x)
        /(ε_y^2+Δ_y^2)^5/2
      ∂_u^2 f_M,z   =
        
          -3ε_zΔ_z^2
        /(ε_z^2+Δ_z^2)^5/2
  
  
  for the second derivatives.
  We again compute explicitly the first and second derivatives of the QFI with
  respect to the driving parameter by rewriting the QFI in terms of an integral
  over a QFI density,
  
      F_B,α = ∫_BZ f_B,α(q;u)
  
  finding,
  
  
      ∂_u f_B,x   = 
        Δ_xε_x(ε_xsin(q_y) -
          Δ_xcos(q_y))
        /
          (ε_x^2+Δ_x^2)^2
        
      ∂_u f_B,y   = 
        Δ_yε_y(ε_ysin(q_x) -
          Δ_ycos(q_x))
        /
          (ε_y^2+Δ_y^2)^2
        
      ∂_u f_B,z   = 
        Δ_z^2ε_z
        /
          (ε_z^2+Δ_z^2)^2
  
  
  for the first derivatives and
  
  
      ∂_u^2 f_B,x   =
        (ε_xsin(q_y)-Δ_xcos(q_y))
          (
            (
              Δ_x^2-3ε_x^2
            )Δ_xcos(q_y)
            +
            (
              ε_x^2-3Δ_x^2
            )ε_xsin(q_y)
          )
        /(ε_x^2+Δ_x^2)^3
        
      ∂_u^2 f_B,y   =
        (ε_ysin(q_x)-Δ_ycos(q_x))
          (
            (
              Δ_y^2-3ε_y^2
            )Δ_ycos(q_x)
            +
            (
              ε_y^2-3Δ_y^2
            )ε_ysin(q_x)
          )
        /(ε_y^2+Δ_y^2)^3
        
      ∂_u^2 f_B,z   =
        Δ_z^2(Δ_z^2-3ε_z^2)
        /(ε_z^2+Δ_z^2)^3
  
  
  
  
  
  
  
  
  