Quenching the Kitaev honeycomb model

I studied the non-equilibrium response of an initial N\'{e}el state under time evolution with the Kitaev honeycomb model. With isotropic interactions ($J_x = J_y = J_z$) the system quickly loses its antiferromagnetic order and crosses over into a steady state valence bond solid, which can be inferred from the long-range dimer correlations. There is no signature of a dynamical phase transition. Upon including anisotropy ($J_x = J_y \neq J_z$), an exponentially long prethermal regime appears with persistent magnetization oscillations whose period derives from an effective toric code.


Introduction
Quantum spin liquids [1] are intriguing forms of matter characterized by the absence of magnetic order and the presence of long-range entanglement. A defining feature is that they cannot be transformed smoothly into a non-entangled magnetic product state, such as the Néel antiferromagnet. One might wonder whether these opposite extremes can be connected under a rapid change of external parameters.
Such a rapid change is known as a quench [2,3], and this set-up has lead to the prediction and observation of dynamical phase transitions. [4][5][6] For example, in the transverse field Ising model, time evolution of an initial magnetic state under a Hamiltonian with a trivial paramagnetic ground state leads to nonanalytic behavior in the return amplitude at certain times after the quench. [4] Also the opposite quench, starting from a paramagnetic spin liquid and time-evolving with a Hamiltonian not supporting spin liquid behavior, has been studied [7]. In this work, I will combine these works to answer the question: what happens when time evolves a magnetic state with a Hamiltonian that has a spin liquid ground state? Will we see a dynamical phase transition or crossover into a spin liquid regime at some finite timescale? Or will signatures of the initial magnetic order remain? The Kitaev honeycomb model [8] provides an ideal playground to answer this question since it is exactly solvable. A slow ramp in this model has been studied before [9,10], but there the dynamics started from an initial spin liquid state. Here, I start from an antiferromagnetic Néel state, the simplest possible non-entangled magnetically ordered state, and time evolve with the Kitaev Hamiltonian with both isotropic (J x = J y = J z ) and anisotropic interactions (J x = J y = J z ). In order to time evolve with the Kitaev model I first express the Néel state as a superposition of different gauge fields configurations. Within each gauge sector, I then compute the exact time evolution of the free Majorana fermions.
As expected, the initial staggered magnetization vanishes (see Fig. 1) after the quench. Unlike quenches in the transverse field Ising model, however, there seem to be no signatures of a dynamical phase transition. At long times the system becomes a steady state valence bond solid. More surprisingly, an exponentially long prethermal regime appears when the interactions are anisotropic, as seen for example in the time evolution of the magnetization (Fig. 1). This prethermal regime is governed by an effective high-temperature toric code.
In Sec. 2 I will first present the Kitaev honeycomb model, its exact solution and an outline of the quench method. The results are discussed in Sec. 3, with a special emphasis on the question of a dynamical phase transition (3.1), the prethermal regime (3.2) and the final steady state (3.3). I conclude with a brief discussion on entanglement and experimental realizations in Sec. 4.

Model, initial state and method
Before presenting the results in detail, let me introduce the set-up of the quench. Consider spin-1 2 degrees of freedom σ i on a honeycomb lattice. The unit cell has two sites, which I will label as the A and B site, shown in Fig. 2. The initial state will be a perfect Néel state polarized along the z-direction, which is an unentangled product state |ψ 0 = i | ↑ iA ⊗| ↓ iB . Starting from this initial state I will compute the time evolution using the Kitaev honeycomb model. In this model the bonds between lattice sites are divided into three types, depending on their direction, as shown in Fig. 2. Each bond-type has an Ising spin interaction along a different spin orientation, Kitaev's key insight was that one can solve this model exactly by representing each spin by four Majorana operators b x , b y , b z and c. This enlarges the Hilbert space, and in the enlarged Hilbert space we can define 'enlarged' spin operators σ x = ib x c, σ y = ib y c, and σ z = ib z c. The projection operator onto the real, physical, subspace is P = 1 2 (1 + b x b y b z c). Therefore, the physical spins are given by σ α = P σ α P , which implies σ In the following, I will use that within the physical subspace, the real spins can also be represented by the operators of the form σ z = −ib x b y and similar expressions hold for σ x and σ y .
In terms of the new Majorana operators, the Hamiltonian reads where j sums over unit cells and u jα = ib α jA b α j+δα,B = ±1 is a static Z 2 gauge field living on the α = x, y, z bond. The product of Z 2 gauge fields along a plaquette is gauge-invariant and is the 'flux' w p = σ x 1 σ y 2 σ z 3 σ x 4 σ y 5 σ z 6 . The remaining c-Majorana's are called 'matter' and are noninteracting.
The spin liquid ground state of the Kitaev honeycomb model is in the zero-flux sector, meaning all gauge fields u jα are the same. In contrast, the Néel state, when expressed in terms of gauge and matter fields, is in a superposition of different flux configurations since ψ 0 |w p |ψ 0 = 0 where w p is the plaquette flux operator. We can show which flux configurations are included in this superposition by repeated use of the fact that the Néel state is an eigenstate of the physical operator σ z i . A good basis to describe the Néel state is by pairing the remaining matter Majorana's along the z-bonds within one unit cell, v j = ic jA c jB = ±1. Any possible state in the enlarged Hilbert space can be written as a superposition of u, v-configurations, |ψ = {u jα ,v j } c {u jα ,v j } | {u jα , v j } , and our task is to find the weight constants c {u jα ,v j } . The fact that the Neel state is physical and therefore must satisfy P j |ψ 0 = |ψ 0 , and that it is an eigenstate of σ z j for every j, leads to two constraints on the possible u, v-configurations. On a lattice consisting of L x × L y unit cells with periodic boundary conditions, we have periodic chains of xy-bonds. The product of all 2L x z-spins along such a xy-chain equals (−1) times the product of all x and y gauge fields. Therefore, this product of gauge fields must equal (−1) Lx+1 . Consequently, the Néel state is an equal-weight superposition of all N c = 2 3LxLy−Ly possible u jα gauge field configurations that satisfy this constraint. The matter content v i is fixed by the constraint σ z jA σ z jB = −1 within each unit cell, which implies u jz = v j . The relative phases between different {u jα , v j }-configurations are fixed by the expectation value of σ z j operators, and are multiples of i.
Note that in principle the gauge freedom allows us to construct the same Néel state with a different set of gauge field configurations. However, the current choice is extremely transparent since it represents the Néel state as an equal superposition of all allowed configurations. This in turn makes the calculation of observables straightforward.
Because the gauge fields are integrals of motion only the matter fields will be changing over time, where {u jα } represents a gauge field configuration that respects the aforementioned constraints, |ψ is the initial matter field configuration determined by v j = u jz and H {u jα } is a free matter Majorana Hamiltonian with hoppings depending on the Z 2 gauge fields. The magnetization on an A lattice site m jA (t) = ψ(t)|σ z jA |ψ(t) can be found using the gaugefield-only representation of spin, σ z jA = −ib x j b y j . Therefore, the magnetization can be written as the return amplitude with two matter Hamiltonians, where the configurations {u jα } and {u jα } differ only by the flip of the two gauge fields u x j and u y j . The sum over exponentially many gauge field configurations can be replaced by a random Monte Carlo sampling over all configurations [11,12] that satisfy the constraints relevant for the initial Néel state. For each such configuration I need to compute these generalized return amplitudes for the matter Hamiltonian, which can be done efficiently using the Balian-Brezin decomposition as outlined in Appendix A. [13,14] Note that an alternative way of deriving my results is by using the 'brick wall'-representation of the Kitaev honeycomb model. [15] Note that in the basis we use, where we pair Majorana matter particles along the zbonds to create complex fermions, it is most natural to compute expectation values of S z and correlations thereof. In contrast, the computation of correlations functions containing S x or S y is more involved. Specifically, such correlation functions can no longer be expressed purely in terms of gauge fields b, and would require keeping track of the time evolution of Majorana matter fermions. Since we start with a z-polarized Néel, we only focus in this manuscript on the dynamics of correlation functions involving only S z operators.
Finally, it is worth mentioning that the sampling over all gauge fields is specific to our choice of initial state. In the extreme case that one has an initial state that lies purely within one flux sector, the corresponding quench dynamics is that of a non-interacting fermion model. This is done in, for example, Refs. [9,10], where they study a ramp within the zero-flux sector. In this case the dynamics follow the general lore of dynamical phase transitions [4][5][6]16]. Much of our results in the following section deviate from this precisely because we intertwined various flux sectors by choosing an initial Néel state.

Results
I will now show results for a quench from an initial Néel antiferromagnet, to the Kitaev honeycomb model. I consider both quenches to the isotropic case (J x = J y = J z ) as well as the anisotropic Kitaev model (J x = J y = J z ). The anisotropy is defined by the ratio J z /J xy where J xy ≡ J x = J y .

Phase transition or crossover?
The dynamics studied here can be viewed as a quench through a quantum critical point separating an antiferromagnetic phase and a paramagnetic spin liquid phase. It has been suggested that a quench from the ferromagnetic to the paramagnetic phase leads to nonanalytic behavior of the return amplitude at given times. [4][5][6]16] Specifically, consider the nonequilibrium free energy density The nonequilibrium free energy f (t) = − 1 N log G(t) for various J xy , fixed J z = 1, N mc = 20000, and system sizes L = 6 (thin lines) and L = 8 (thick lines), where L x = L y = L. No dynamical phase transition is observed in the short times before a steady state plateau emerges. A prethermalization regime appears with increasing anisotropy. Right: The nonequilibrium free energy for the isotropic case J z = J xy as a function of system size L, averaged over N mc = 500, 000 gauge configurations. I extrapolated these results to L = ∞, which still does not show any signatures of nonanalytic behavior.
where N = L x L y = L 2 , and G(t) is the return amplitude or Loschmidt echo In the transverse field Ising model, this quantity is shown to be nonanalytical at several moments after the quench. Such nonanalytic points are called dynamical phase transitions, in analogy to thermal phase transitions where the free energy becomes nonanalytic. In order to study the possible appearance of a dynamical phase transition in our quench model, I computed the nonequilibrium free energy density f (t). The results are shown in Fig. 3. On the left-hand side I show the free energy, for J z = 1 and various J xy , averaged over N mc = 20000 gauge configurations. Independent of J xy , there is an initial growth of free energy. Subsequently, a plateau (discussed in Sec. 3.2) appears for the anisotropic cases. After that, there seem to be severe system-size fluctuations and it is not apparent whether or not a true nonanalyticity appears.
On the right of Fig. 3 I show the system-size dependence of the free energy for the isotropic case, computed using much more gauge configurations (N mc = 500, 000). Around the time t = 10 there seems to be a steady state plateau developing for the free energy, which is strongly system size dependent. The infinite L limit, however, does not seem to suggest any nonanalytic behavior. The evolution of the free energy is likely more accurately described as a crossover.
It is important to note that many dynamical phase transitions have been found in models that are noninteracting, such as the transverse field Ising model, [4][5][6]16] the XY model, [17] or fermionic band insulators. [18] Even though the Kitaev model is exactly solvable, it is not a noninteracting theory. This might be the reason why I do not observe any dynamical phase transition.  Figure 5: The static spin-spin correlations σ z i σ z j for various short-range spins, in the isotropic model, with L = 8 and N mc = 2000. After a short time all spin correlations vanish except the nearest-neighbor correlation along a z-bond. The long-range dimer-dimer correlation function, here measured at the longest possible distance between two sets of zbonds, also obtains a nonzero steady state value. These correlations are indicative of a valence bond solid phase. [19]

Prethermalization
Upon increasing the anisotropy J z /J xy , a plateau emerges in the free energy (Fig. 3) that lasts long in the anisotropy ratio J z /J xy . During this exponentially long timescale, persistent oscillations in the staggered magnetization m(t) = j (−1) j σ z j (t) appear. To emphasize this behavior, we show in Fig. 4 the staggered magnetization for J xy = 0.2J z as a function of system size, including a L = ∞ limit. Even though the anisotropy is only J z /J xy = 5, the time-scale over which the magnetization persists is about 10 5 longer than for the isotropic case. Different measures of a typical time-scale, namely the onset of the free energy plateau or when the magnetization reaches a 0.2 or 0.001 threshold, all display an approximately exponential dependence on the anisotropy t * ∼ e cJz/Jxy , as is shown in the inset of Fig. 1.
Both of these phenomena -the long time-window t * and the persistent magnetization oscillations -can be understood within the framework of prethermalization. [20][21][22][23][24][25] Let us first consider the length of the time-scale t * . In typical quenched systems, dynamical time-scales would depend in a power-law fashion on an anisotropy parameter. For example, in a quench of the transverse field Ising model from the ferromagnetic to the paramagnetic phase, the typical time-scale is set by t * = π/ k * (g 1 ), where k (g) = (g − cos k) 2 + sin 2 k, cos k * = 1+g 0 g 1 g 0 +g 1 and g 0 ,g 1 are the values of the transverse field before and after the quench. [4] This timescale diverges as g 1 , the post-quench transverse field, becomes close to the critical value g c = 1. It is easy to show that this divergence indeed follows a power-law, t * ∼ (g 1 − g c ) −1/2 . So why is the time-scale t * so much longer in the case of quenching the anisotropic Kitaev model? The answer lies in the concept of prethermalization [20][21][22][23][24][25] that occurs in systems close to integrability. In particular, the dynamics here can be understood using the framework of Ref. [23]. For the anisotropic Kitaev model, the coupling along the z-bonds is significantly stronger than along the x, y-bonds. We can therefore treat the x, y-bond coupling as a perturbation. The Kitaev model is written aŝ whereN is the sum of all the z-bond couplings, andŶ contains the couplings along the x, y-bonds. The termN is trivially integrable, it is just a sum of local commuting terms with integer eigenvalues. Following Ref. [23], for J xy < J z , we can perform a unitary transformation such that the Hamiltonian becomeŝ where the new termĤ eff commutes withN . This means thatĤ eff does not affect the relative orientation of two spins along a z-bonds. More importantly, the remaining term is exponentially small, meaning that for an exponentially long time the dynamics preserve the spin correlations along each z-bond. To summarize, an exponentially long timescale t * ∼ e cJz/Jxy appears because for a suitable unitary transformation the Hamiltonian has effectively new conservation laws that constrain dynamics for a long time.
The fact that the spin correlations are locked along a z-bond can be inferred from measuring the static spin correlation function S zz ij (t) = ψ(t)|σ z i σ z j |ψ(t) . As shown in Fig. 5, even in the isotropic case the relative orientation of spins along a z-bond remains nonzero in the infinite time-limit. This can be further corroborated by computing the static spin correlations in the diagonal ensemble (see Appendix B), which indeed yields a steady state with zero magnetization but nonzero spin correlations along the z-bond. Notice that other static spin correlations vanish.  Figure 6: In the anisotropic limit J xy J z , the exponentially long prethermal regime is governed by an effective toric code. This causes persistent oscillations in the magnetization with a period T (J xy ) given by Eq. (10). Here we show the magnetization oscillations in the range J xy /J z = 0.05−0.3 for linear system size L = 6, with N mc = 2000. The time (horizontal axis) is rescaled for each value of J xy to correspond to exactly four oscillation periods. We find that indeed for small J xy , we approach the predicted periodicity.
Having established the z-bond spin-lock for an exponentially long time, we can investigate how this leads to persistent magnetization oscillations during this prethermal regime. The key lays in the effective HamiltonianĤ eff in Eq. (8), which describes the effective interaction between the locked spins along a z-bond. At each z-bond, the configuration must be antiferromagnetic, meaning that only the spin configurations ↑↓ and ↓↑ are allowed. These two states constitute a 'new spin' τ , and using Kitaevs fourth order perturbation theory [8] the effective Hamiltonian becomeŝ where p is every plaquette of the honeycomb lattice, and left/right/etc. refers to the z-bond to the left/right/etc. of this plaquette. Note that this model is equivalent to the toric code. [26] Any attempts to understand the dynamics in terms of topology are futile, as we are very far from the ground state during our quench dynamics and signatures of topology, such as anyonic excitations, are defined only close to the ground state.
Nonetheless, a simple calculation shows what would happen if one starts with an initial Néel state (meaning ↑↓ on every z-bond) and let time evolve with the Hamiltonian of Eq. (9). In such a quench, the staggered magnetization will oscillate according to m(t) = cos 2 (2J eff t) where J eff = To confirm this prediction, we analyzed how the period of persistent oscillations varies with J xy /J z . The results are shown in Fig. 6 for L = 6, N mc = 2000 and varying small J xy . Whenever J xy < 0.2J z , the magnetic oscillations have a period that is well approximated by Eqn. (10). We conclude that indeed, for an exponentially long time, the system undergoes magnetic oscillations as dictated by the toric code.

Steady state valence bond solid
In the isotropic case there is no signature of prethermalization and after a short time of order unity the system equilibrates. While there is zero net staggered magnetization in this steady state, there are remnant nearest-neighbor spin correlations along the z-bond, as shown in Fig. 5. This suggests the steady state is a valence bond solid with the singlets oriented along the z-bonds. To further corroborate this claim, I studied the dimer-dimer correlation function D zz ij (t) = σ z i σ z i+δz σ z j σ z j+δz that has been used before as an indication of valence bond order. [19] Indeed, I find long-range dimer order, even though the state is at relatively high temperatures. Notice that this state does break rotational invariance, since the singlet bonds live on the z-bonds which are inequivalent to the x, y-bonds of the lattice.
Another way to quantify the steady state is through the dynamic two-time spin correlation function S zz j (t, t ) = ψ(t)|σ z jA (t )σ z jB |ψ(t) . [27][28][29] The Fourier transform with respect to t −t can be interpreted as an AC spin conductivity. Specifically, the DC (ω = 0) response measures the antiferromagnetic correlations along a z-bond. For small ω the correlations are reduced over a frequency scale set by the flux-averaged Majorana density of states, see Fig. 7. At later times this correlations gets suppressed in the frequency range between 0 and 6J, which is the flux-averaged bandwidth of the matter Majorana's. Interestingly, for times t > 0.5 a reversal of the dynamic correlations for 4J < ω < 6J appears. This is a signature of the elementary triplet excitation of the valence bond solid.
We thus find a dynamic crossover from a Néel state to a valence bond solid. This transition in equilibrium is known as a deconfined quantum phase transition and falls outside the usual Landau classification of continuous phase transitions. [30] The absence of any finite time singularity is due to the fact that the location of the valence bonds are determined by the orientation of the initial Néel state. There is no dynamical spontaneous symmetry breaking: a Néel state polarized along the x-axis would give rise to a valence bond solid with singlets along the x-bonds, and so forth. It is an interesting open question to study what happens when an initial Néel state is not aligned along one of the principal spin axes.

Conclusion and Discussion
I showed that starting from a Néel state, time evolution with the Kitaev honeycomb model leads to crossover to a steady state valence bond solid. When the interactions are anisotropic (J z /J xy = 1), an exponentially long prethermal regime appears whose dynamics can be effectively described by a toric code.
Note that similar results are expected if one would start with an initial ferromagnetic product state, rather than an antiferromagnetic product state.
To what extent these results remain valid beyond the exactly solvable model, for example by introducing a small Heisenberg term, is an open question. Based on the proof of Ref. [23] I expect that the prethermal regime will persist even in the presence of such perturbations, but quantifying this requires new computational techniques beyond the ones used in this work.
An interesting aspect that was not included in this study is the topological nature of the Kitaev honeycomb model. The ground state of the model has nontrivial entanglement entropy, [31,32] a topological groundstate degeneracy and has anyonic excitations. [8] It is hard, however, to see one of these topological signatures in the quench set-up. After all, the Néel state has a high energy density expectation value in the Kitaev honeycomb model, meaning that the quench dynamics are far away from the ground state. The generated entanglement is therefore volume law, making it almost impossible to detect a possible topological entanglement entropy. The same holds for anyonic excitations, which are well-defined only close to the ground state. Interestingly, topological degeneracy might be observed. On a torus, the Kitaev honeycomb model in the toric code regime has a fourfold degeneracy, meaning that there are orthogonal ground states that are locally indistinguishable. Acting with a string of spin operators allows you to go from one to the other ground state. Now in the final valence bond steady state, acting with a vertical string of σ x and σ y operators will create an orthogonal state that is indistinguishable from the original state on the level of single-spin operator measurements. I leave it for future work to investigate whether indeed this amounts to a topological degeneracy, which might be relevant for quantum computations.
Finally, in recent years some material systems have been proposed to be experimental realizations of the Kitaev honeycomb model [33]. Though straining these materials is unlikely to give rise to the desired anisotropy to observe a prethermal regime, it might be possible to chemically engineer these system to get the anisotropic interactions desired. It will also be interesting to see the dynamic response after a quench with an initial state resembling the spiral magnetic order found in these materials. [34,35] Funding information This research was supported by Perimeter Institute for Theoretical Physics. Research at Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science and Economic Development Canada and by the Province of Ontario through the Ministry of Research, Innovation and Science.

A Matter Hamiltonian time evolution
As described in the main text, the time evolution of the Kitaev honeycomb model is completely due to the Majorana fermions. In each unit cell j, which contains a z-link, we indentify an A and B sublattice site. The c-Majorana's in Kitaev's notation are then paired along the z-link to form complex fermions, The matter Hamiltonian on the full honeycomb lattice reads where j labels a unit cell, and δ connects to the unit cell with center at position δ x = 1 2 (− √ 3x − 3ŷ), and δ y = 1 2 ( √ 3x − 3ŷ)). In each gauge sector, the required initial state is the product state where unit cells with u z j = 1 are occupied with a complex matter fermion. For later purposes it is practical to perform a particle-hole transformation on 'occupied' sites, so that the Hamiltonian becomes With the Hamiltonian Eqn. (14), the initial matter state is nothing but the a-vacuum |0 , defined by a j |0 = 0. This matter Hamiltonian can be brought into a canonical Bogoliubov-De Gennes (BdG) format, where H d is a real-valued symmetric matrix, ∆ a real-valued antisymmetric matrix, and the vector a † a contains all creation and annihilation operators for all unit cells. The 2N × 2N BdG matrix in Eqn. (15) can be diagonalized, H BdG = V ΛV , with real eigenvalues Λ = diag ( 1 , 2 , . . . , N , − 1 , . . . , − N ) and V a real orthogonal matrix of the form V = Q R R Q . This diagonalization allows us to compute the Balian-Brezin decomposition of the time evolution operator, [13,14] e −iHt = e 1 2 a † Xa † e a † Y a e 1 2 aZa det Re −iΛt/2 + Qe iΛt/2 (16) where The simplest quantity to compute is the overlap of the initial state with the time-evolved state, known as the return amplitude G(t) = ψ(t)|e −iHt |ψ 0 . Because different flux sectors are orthogonal to one another, the total return amplitude is a sum of matter Majorana return amplitudes in each gauge sector, Note that due to the particle-hole transformation, the state |ψ {u jα } 0 is equal to the a-vacuum, so a j |ψ {u jα } 0 = 0 for every a j . To simplify notation, from now on I will write |0 for the initial state.
Since the number of gauge field configurations scales exponentially with system size, it is impossible to compute the above sum of Eqn. (17) exactly. Instead, I averaged over N mc random gauge field configurations that satisfy the constraints set by the initial state. It turns out that N mc = 1000 yields sufficient accuracy for the system sizes considered. The staggered magnetization, defined as m(t) = 1 2N j ψ(t)|σ z jA − σ z jB |ψ(t) , will decay over time starting from m(t = 0) = 1. Using the representation σ z jA = −ib x j b y j , valid within the physical subspace, we see that the magnetization can be computed as a sum over return amplitudes involving two Hamiltonians, where H 1 and H 2 only differ through a flip of the u jx and u jy gauge fields neighboring the spin that we want to measure. I proceed by making the Balian-Brezin decomposition for both H 1 and H 2 , R(t) = det R 2 e iΛ 2 t/2 + Q 2 e −iΛ 2 t/2 det R 1 e −iΛ 1 t/2 + Q 1 e iΛ 1 t/2 0|e 1 2 aZ * 2 a e 1 2 a † X 1 a † |0 . (20) state is then |ψ(t) = n c n e −iEnt |n . The diagonal ensemble is a density matrix composed of the time-independent diagonal of the initial state density matrix, ρ D = n |c n | 2 |n n|.
In our case, the eigenstates of our system have the form |{u} ⊗ |{f } where |{f } are Fock states composed of the single-particle wavefunctions diagonalizing the matter BdG Hamiltonian. Since the flux sectors are orthogonal we can construct a diagonal ensemble within each flux sector. The trace carries over to the extended Hilbert space provided we use the physical subspace projector and our initial state is completely embedded in the physical subspace. Any operator that changes the flux sector, such as an isolated σ z j , must have a zero expectation value in the diagonal ensemble. The only (possibly) nonzero expectation values of two-spin operators are of of the form σ z jA σ z jB along a z-bond. Using σ z A σ z B = ib z A c A ib z B c B = −iu z c A c B , and the particle-hole transformation defined above, I find Tr σ z jA σ z jB ρ D = 1 − 2Tr a † j a j ρ D where I used the diagonalization of the matter Hamiltonian defined in the previous section.