Thermodynamics of a gauge-frustrated Kitaev spin liquid

TL;DR

This paper explores how geometric frustration in a Kitaev spin liquid suppresses gauge field ordering, leading to residual entropy and a Majorana metal phase, with thermodynamic signatures analyzed via quantum Monte Carlo simulations.

Contribution

It introduces a frustrated gauge sector in a Kitaev model, revealing suppressed gauge ordering and a novel Majorana metal phase, extending understanding of spin liquid thermodynamics.

Findings

Gauge frustration suppresses gauge ordering transition.

Residual entropy persists at low temperatures.

Formation of a Majorana metal phase at very low temperatures.

Abstract

Two- and three-dimensional Kitaev magnets are prototypical frustrated quantum spin systems, in which the original spin degrees of freedom fractionalize into Majorana fermions and a $\mathbb{Z}_2$ gauge field -- a purely local phenomenon that reveals itself as a thermodynamic crossover at a temperature scale set by the strength of the bond-directional interactions. For conventional Kitaev magnets, the low-temperature thermodynamics reveals a second transition at which the $\mathbb{Z}_2$ gauge field orders and the system enters a spin liquid ground state. Here we discuss an explicit example that goes beyond this paradigmatic scenario -- the $\mathbb{Z}_2$ gauge field is found to be subject to geometric frustration, the thermal ordering transition is suppressed, and an extensive residual entropy arises. Deep in the quantum regime, at temperatures of the order of one per mil of the interaction strength, the degeneracy in the gauge sector is lifted by a subtle interplay between the gauge field and the Majorana fermions, resulting in the formation of a Majorana metal. We discuss the thermodynamic signatures of this physics obtained from large-scale, sign-free quantum Monte Carlo simulations.

Equations

$${\mathcal{H}}_{\mathrm{Kitaev}}=\sum _{⟨j,k⟩,\gamma }{J}_{\gamma }{\sigma }_j^{\gamma }{\sigma }_k^{\gamma },$$

$$\mathbf{W}=\left(\begin{array}{c}{W}_x\\ W_y\\ W_z\end{array}\right)\stackrel{(\mathrm{a})}{=}\left(\begin{array}{c}-1\\ -1\\ 1\end{array}\right)\stackrel{(\mathrm{b})}{=}\left(\begin{array}{c}1\\ -1\\ -1\end{array}\right)\stackrel{(\mathrm{c})}{=}\left(\begin{array}{c}-1\\ 1\\ -1\end{array}\right),$$

$$P=\frac{4}{3N}\sum _j⟨{\mathbf{W}}_0\cdot {\mathbf{W}}_j⟩,$$

$$C_{v,\mathrm{MF}}(T)=-\frac{1}{T^2}{⟨\frac{\partial E_f(\{u_{jk}\})}{\partial \beta }⟩}_{\mathrm{MC}},$$

$$C_{v,\mathrm{GF}}(T)=\frac{1}{T^2}\left({⟨E_f^2(\{u_{jk}\})⟩}_{\mathrm{MC}}-{⟨E_f(\{u_{jk}\})⟩}_{\mathrm{MC}}^2\right).$$

$${\text{a}}_1=\left(1,0,0\right),{\text{a}}_2=\left(-\frac{1}{2},\frac{\sqrt{3}}{2},0\right),{\text{a}}_3=\left(0,0,\frac{2}{5}\right)$$

$${\text{b}}_1=2\pi \left(1,\frac{1}{\sqrt{3}},0\right),{\text{b}}_2=2\pi \left(0,\frac{2}{\sqrt{3}},0\right),{\text{b}}_3=2\pi \left(0,0,\frac{5}{2}\right).$$

$${\text{r}}_1$$

$${\text{r}}_5$$

$${\bf P}=\left(\begin{array}[]{ccc}\frac{2}{\sqrt{6}}&-\frac{1}{\sqrt{6}}&-\frac{1}{\sqrt{6}}\\ 0&\frac{1}{\sqrt{2}}&-\frac{1}{\sqrt{2}}\\ \frac{1}{\sqrt{3}}&\frac{1}{\sqrt{3}}&\frac{1}{\sqrt{3}}\end{array}\right)\left(\begin{array}[]{c}W_{x}\\ W_{y}\\ W_{z}\end{array}\right)\,,$$

$${\bf P}_{a}=\left(\begin{array}[]{c}-\sqrt{\frac{2}{3}}\\ -\sqrt{2}\\ -\frac{1}{\sqrt{3}}\end{array}\right),\hskip 5.69046pt{\bf P}_{b}=\left(\begin{array}[]{c}2\sqrt{\frac{2}{3}}\\ 0\\ -\frac{1}{\sqrt{3}}\end{array}\right),\hskip 5.69046pt{\bf P}_{c}=\left(\begin{array}[]{c}-\sqrt{\frac{2}{3}}\\ \sqrt{2}\\ -\frac{1}{\sqrt{3}}\end{array}\right)\,.$$

$$\frac{{\bf P}_{m}\cdot{\bf P}_{n}}{\|{\bf P}_{m}\|\|{\bf P}_{n}\|}=\frac{{\bf W}_{m}\cdot{\bf W}_{n}}{\|{\bf W}_{m}\|\|{\bf W}_{n}\|}=\left\{\begin{array}[]{c}-\frac{1}{3},\hskip 5.69046ptm\neq n\\ 1,\hskip 14.22636ptm=n\end{array}\right.$$

$$P_{abc}=\frac{1}{9}\left(1+1+1-6\cdot \frac{1}{3}\right)=\frac{1}{9}.$$

$$P_{bc}=\frac{1}{4}\left(1+1-2\cdot \frac{1}{3}\right)=\frac{1}{3},$$

$$P_a=1.$$

$$P=\frac{1}{N_s}\sum _j\frac{⟨{\mathbf{W}}_0\cdot {\mathbf{W}}_j⟩}{\Vert {\mathbf{W}}_0\Vert \Vert {\mathbf{W}}_j\Vert }=\frac{4}{3N}\sum _j⟨{\mathbf{W}}_0\cdot {\mathbf{W}}_j⟩$$

Full text

In frustrated magnetism, lattice gauge theories are an ubiquitous tool to capture the physics of quantum spin liquids Savary and Balents (2017). The fundamental distinction in these theories between confined and deconfined regimes corresponds to the formation of trivial, magnetically ordered states versus macroscopically entangled spin liquids, respectively. The principal nature of the underlying gauge theory can further be used to broadly categorize different types of quantum spin liquids such as $\mathbb{Z}_{2}$ spin liquids Read and Sachdev (1991); Senthil and Fisher (2000), $U(1)$ spin liquids Anderson (1973), or chiral spin liquids Kalmeyer and Laughlin (1987) – for which the corresponding gauge theory exhibits either a discrete $\mathbb{Z}_{2}$ or continuous $U(1)$ symmetry or an underlying Chern-Simons action Wen (2002). This classification allows one to immediately draw conclusions about the stability of the corresponding spin liquids, in particular to thermal fluctuations. While the (spontanteous) breaking of time-reversal symmetry in chiral spin liquids mandates their thermal stability and the presence of a finite-temperature phase transition, a more complex picture emerges for $\mathbb{Z}_{2}$ and $U(1)$ spin liquids. Here, spatial dimensionality needs to be taken into account. For the $\mathbb{Z}_{2}$ spin liquid the elementary vison excitations of the underlying gauge structure are point-like objects in two spatial dimensions allowing them to proliferate at finite temperatures and destroy the entangled spin liquid state. In contrast, in three spatial dimensions the $\mathbb{Z}_{2}$ spin liquid is stable to thermal fluctuations, as now the visons form (small) loop-like objects that cannot destroy the spin liquid and break open into extended line-like objects only at a finite-temperature transition. For $U(1)$ spin liquids the elementary instanton excitations of the bare $U(1)$ gauge theory are point-like objects in both two and three spatial dimensions, implying that these spin liquids are generically not stable at finite temperatures 111Possible exceptions, guided by physical intuition, are discussed in the context of quantum spin ice Savary and Balents (2013).

In this manuscript, we consider the explicit example of a three-dimensional $\mathbb{Z}_{2}$ spin liquid, realized in a numerically tractable Kitaev model, that proves to be an exception from these paradigms. At sufficiently low temperature the gauge field is found to be subject to geometric frustration, arising from local constraints that impose a divergence-free condition on the gauge field and ultimately result in an extensive residual entropy. The net result is a suppression of the expected thermal ordering transition of the $\mathbb{Z}_{2}$ gauge field and the emergence of a spin liquid state that, in some sense, is “doubly frustrated”, as it arises from the interplay of exchange frustration on the level of the original spin degrees of freedom and geometric frustration on the level of the emerging fractional degrees of freedom. Extensive quantum Monte Carlo simulations reveal that at ultralow temperatures of the order of one per mil of the interaction strength, i.e., deep in the quantum regime, the degeneracy in the gauge sector is eventually lifted by a subtle interplay with the Majorana fermion degrees of freedom, which emerge in parallel with the gauge field upon fractionalization of the original spin degrees of freedom in Kitaev models. The formation of a collective ground state of these itinerant fermions, a Majorana metal with a distinct nodal line structure feeds back into the gauge sector and leads to the concurrent formation of columnar ordering of the gauge field. Our model system thereby proves to be a principal example of a spin liquid, for which not only the phenomenon of fractionalization, but also of the subsequent non-trivial interplay of the emergent fractional degrees of freedom and the underlying lattice gauge theory can be captured by numerically exact simulations.

Gauge frustration.– The Kitaev model with its characteristic bond-directional spin exchanges of the form

[TABLE]

is well known to be analytically tractable for a class of two-dimensional Kitaev (2006); Hermanns et al. (2018) and three-dimensional O’Brien et al. (2016) lattice geometries. The analytical treatment relies on a parton construction introduced by Kitaev that decomposes the spin degrees of freedom into itinerant Majorana fermions and a static $\mathbb{Z}_{2}$ gauge field. The fact that the gauge field remains static and assumes, at sufficiently low temperatures, an ordered ground state is key for the exact solvability of the model, since it allows one to reduce the problem to one of free Majorana fermions hopping in a fixed background. In fact, a powerful theorem by Lieb Lieb (1994) describes the ground state of the gauge sector in terms of $\mathbb{Z}_{2}$ fluxes through the elementary plaquettes – plaquettes of length $6,10,\ldots$ are flux-free, while plaquettes of length $4,8,\ldots$ carry a $\pi$ flux. Recent classification work of 3D Kitaev models O’Brien et al. (2016) has shown that Lieb’s theorem generically predicts the correct ground state flux assignment, even for lattices that do not fulfill all the mathematical requirements for the theorem to apply. This is also true for the (8,3)c lattice.

The key ingredient for the study at hand is a “frustrated” three-dimensional lattice geometry whose central motif are hexagonal sites at which three plaquettes of length 8 meet, see the illustration in Fig. 1. Following the above intuition based on Lieb’s theorem each of these plaquettes is destined to carry a $\pi$ flux, which however conflicts with the fact that for any closed volume, such as the one spanned by the three neighboring plaquettes, the fluxes must obey a divergence-free condition – if a flux enters the volume through one of the plaquettes, it must leave through another one. This divergence-free condition allows only two of the three plaquettes to carry a $\pi$ flux and leaves one of the plaquettes in a flux-free state. For isotropic coupling strength $J_{x}=J_{y}=J_{z}$ this produces three possible flux arrangements per hexagonal site and an extensive residual entropy for the entire system. It is this formation of an extensive manifold of (almost) degenerate states in the gauge sector that designates the term “gauge-frustrated” for the Kitaev model at hand 222A similar scenario of “gauge frustration” plays out in the 3D Kitaev model on the hypernonagon lattice Kato et al. (2017) where the pattern of $\pm\pi/2$ fluxes arising from elementary 9-bond loops is subject to geometric frustration..

One way to relieve the frustration in the system is to vary the relative coupling strengths of the bond-directional exchange (keeping an overall normalization $J_{x}+J_{y}+J_{z}=1$). To see this, consider that every length-8 plaquette consists of an uneven number of bond-directional coupling types, e.g., $3\times J_{x},3\times J_{y},2\times J_{z}$ for the bottom plaquette illustrated in Figs. 2(a)-(c), while the two upper plaquettes have $3\times J_{z}$ couplings. If the $z$-bond coupling is enhanced, i.e. for $J_{z}>J_{x}=J_{y}$, one finds that the local threefold degeneracy is immediately lifted and only one local gauge configuration, illustrated in Fig. 2(a) is favored. For $J_{z}<J_{x}=J_{y}$ the two flux configurations of Figs. 2(b),(c) remain degenerate, thus only partially lifting the original threefold degeneracy.

To check that this phenomenon of gauge frustration, which in the above line of arguments is primarily motivated by the intuition build on Lieb’s theorem, indeed plays out in the model at hand, we have performed large-scale sign-free Monte Carlo simulations of the model over a wide range of temperatures. To capture the local gauge physics, we define for any given hexagonal site a pseudospin vector

[TABLE]

where the individual components $W_{x,y,z}$ can take values +/- 1 indicating the absence/presence of a $\pi$ flux in the three adjacent plaquettes, thus allowing for eight different possible vectors. For the three states that fulfill the local divergence-free condition, their $\pi$ flux assignments are given on the r.h.s. of the above equation in correspondence with Figs. 2(a)-(c). Using these pseudospin vectors we can define a two-point correlation function

[TABLE]

where [math] and $j$ denote two hexagonal sites of the lattice. $P$ readily reveals the nature of the ground-state manifold and can be directly probed in our Monte Carlo simulations. Its expectation value is $P=1$ for the case of a single ground-state of the gauge field ($J_{z}>J_{x},J_{y}$), and $P<1$ for the extensively degenerate cases, specifically $P=1/3$ for the scenario with a local twofold degeneracy ($J_{z}<J_{x},J_{y}$) and $P=1/9$ for the scenario with a local threefold degeneracy expected for isotropic coupling strengths ($J_{z}=J_{x}=J_{y}$) – see the Supplemental Material for an analytical derivation. Numerical results from Monte Carlo runs are shown in Fig. 2(d) for different strengths of $J_{z}$. The data clearly shows that down to temperatures of the order of $10^{-2}$ (in units of the Kitaev coupling $J$) the pseudospin correlation function goes to zero, indicating a completely disordered state of the gauge fields. Below this temperature scale, the pseudospin correlation function rises and indeed saturates. These simulations thereby unambiguously confirm that the system indeed enters a regime of gauge frustration at low temperatures, with an extensive degeneracy building up in the gauge sector.

Lifting of gauge degeneracy.–

The formation of an “accidential” degeneracy, i.e., a degeneracy that is not protected by any inherent symmetries of the system, is often accompanied by some residual effect that splits this degeneracy, at sufficiently small temperature scales, in favor of a unique (or less degenerate) ground state – an effect that typically goes hand-in-hand with a macroscopic phase transition. Such residual effects can include the energetic or entropic selection of ground states, driven either by otherwise negligible interactions (such as e.g., longer-range interactions) or thermal fluctuations in an order-by-disorder scheme Villain et al. (1980).

For the Kitaev system at hand, we find the particularly intriguing scenario that it is an (energetic) interplay between the emergent fractional degrees of freedom that ultimately lifts the gauge frustration discussed above. From the perspective of the itinerant Majorana fermions, the residual degeneracy in the gauge sector is equivalent to a complex scattering potential, as every individual gauge configuration corresponds to a distinct sign structure of the Majorana hopping amplitudes. In the gauge frustrated regime, the collective state of the Majorana fermions is therefore best described as a thermal metal Chalker and Coddington (1988); Chalker et al. (2001), as the degeneracy in the gauge sector has a similar effect as (thermal) disorder. This observation readily points to a scenario where the formation of a collective Majorana state – a more conventional, disorder-free metallic state – might become favorable at the expense of inducing an ordering in the gauge sector. This is precisely what happens at ultralow temperatures, of the order of $10^{-3}$ of the magnetic coupling strength, in the system at hand – the itinerant Majorana degrees of freedom form a nodal-line semimetal, while simultaneously enforcing a columnar ordering in the gauge sector that lifts the gauge frustration. Schematically, the key signatures of these states are illustrated in Fig. 3, which shows the gapless nodal line in the Majorana band structure for different values of the bond-directional exchange $J_{z}$, and the corresponding columnar ordering patterns of the gauge field.

Thermodynamics.–

To quantitavily probe this physics we have measured a variety of thermodynamic observables in quantum Monte Carlo (QMC) simulations covering four orders of magnitude in temperature. These QMC simulations are performed in the sign-free parton basis Nasu et al. (2014), i.e. we sample configurations of the gauge field, $\{u_{jk}=\pm 1\}$ for every bond $\langle j,k\rangle$ of the lattice, with the change of the Majorana free energy, $F_{f}(\{u_{jk}\})$, being calculated explicitly in every update step (either by exact diagonalization or a Green’s function based kernel polynomial method Weiße (2009); Mishchenko et al. (2017), see the Supplemental Material for further details of this implementation). This procedure also allows us to separately distill the entropic contributions to the specific heat of the Majorana fermions Shimomura et al. (2005)

[TABLE]

and the gauge field

[TABLE]

Results are given in Fig. 4(a) for the isotropic coupling point ($J_{x}=J_{y}=J_{z}$). Some features of the multipeak structure of the specific heat are well known from conventional Kitaev models, such as the crossover feature at temperatures of order 1, where the system releases about half of its entropy, see Fig. 4(b), upon the fractionalization of the local spin degrees of freedom Nasu et al. (2014, 2015), primarily by the Majorana fermions (whose energy scale is set by the hopping/magnetic coupling strength). Below this crossover peak there are two additional features in the specific heat that originate in the gauge sector. At a temperature of about $10^{-2}$ a broad shoulder forms, which does not show any scaling with system size pointing to a local crossover phenomenon. It is at this temperature scale that the system enters the manifold of “gauge frustrated” states, which is evident from (i) the average plaquette flux dropping to a value of $\langle W_{p}\rangle=-1/3=(-2+1)/3$ (expected for local configurations where two out of three plaquettes have a $\pi$ flux, i.e., $W_{p}=-1$, and one plaquette remains flux-free, $W_{p}=+1$), (ii) the fluctuations of the plaquette flux exhibiting a maximum upon entering this constrained ground-state manifold, and (iii) the pseudospin correlator (3) raising and saturating at the expected value of $P=1/9$, as documented in Figs. 4(c)-(e). Below this second crossover peak, at a temperature of the order of $2\times 10^{-3}$, one finds a sharp peak in the specific heat that sharpens with increasing system size – this is a true thermal phase transition, where the system releases entropy by forming a columnar ordering of the gauge field. In this ordered state every column of hexagonal sites exhibits a staggered pattern of the flux-free plaquettes as indicated by the yellow plaquettes in Fig. 3(b). However, since the columns order individually and there are two possible staggered states for each column, the resulting overall order not only still allows for a residual entropy, but also breaks the lattice rotational symmetry. This is a remarkable symmetry-breaking effect as it plays out solely in the gauge sector.

Phase diagram.–

Expanding this analysis of key thermodynamic observables to a range of $J_{z}$ parameters, we have compiled the composite phase diagram of Fig. 5. Plotted here are different indicators for the crossover scale to the constrained gauge manifold. As proxies for this crossover, we have marked (i) the location of the peak in the variance of the fluxes, akin to Fig. 4(d), by the filled green circles and (ii) the line of temperature points at which the pseudospin correlator (3) crosses $P=1/9$ by the filled white circles. The low-temperature phase transition, at which the concurrent formation of a nodal-line Majorana semimetal and columnar order of the gauge field occurs, is indicated by the white squares. Depending on the strength of $J_{z}$, we distinguish two principal scenarios. First, there is a line of transitions (indicated by the filled white squares) where it is the formation of the Majorana metal that lifts the degeneracy in the gauge sector and enforces the columnar gauge order. This is the case for $J_{z}\leq 1/3$. For $J_{z}\gtrsim 0.40$, it is the energetics within the gauge sector that readily selects a single configuration of the constrained gauge field for each hexagonal site, see Fig. 2(a), which results in the columnar ordering depicted in Fig. 3(f). For $1/3<J_{z}\lesssim 0.40$, a more subtle mechanism is at play where the energetics of the gauge field favors the same type of columnar order as for $J_{z}\gtrsim 0.40$, but the minimization of the Majorana energy enforces yet another type of columnar order, depicted in Fig. 3(d), which, in a certain sense, is an intermediate type of order with a staggered flux pattern involving flux-free states on some of the bottom $W_{z}$ plaquettes. The overall phase diagram thereby reveals multiple distinct regimes, in which a subtle interplay between the emergent parton degrees of freedom leads to the formation of different types of collective ground states – including gapless spin liquids with a Majorana nodal line and columnar-ordered $\mathbb{Z}_{2}$ gauge fields.

Conclusions.–

The main results of the study at hand are two advances in the conceptual understanding of quantum spin liquids. First, we have introduced the concept of “gauge frustration”, which we showcased in a three-dimensional Kitaev model where the emergent $\mathbb{Z}_{2}$ gauge degrees of freedom are subject to local constraints resulting in an extensive residual entropy. Second, we showed by large-scale numerical simulations that this residual entropy can be lifted by an interplay of the $\mathbb{Z}_{2}$ lattice gauge theory and the itinerant Majorana fermions, which concurrently emerge with the gauge field upon fractionalization of the original local spin degrees of freedom. As such, the model at hand realizes a scenario intermediate between more conventional Kitaev models where the parton degrees of freedom fully decouple (allowing for an analytical solution where one first identifies the ground state of the gauge field and subsequently solves the Majorana problem), and the scenario of strongly interacting partons as it is the case for, e.g., a $U(1)$ spin liquid, in which the gauge field remains heavily fluctuating to the lowest temperatures and thereby strongly feeds back into the formation of a collective parton state.

Acknowledgements.

Acknowledgments.– T.E., M.H., and S.T. acknowledge partial funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Projektnummer 277146847 – SFB 1238 (projects C02 and C03). M.H. acknowledges partial funding by the Knut and Alice Wallenberg Foundation and the Swedish Research Council. P.A.M, T.A.B, Y.K., and Y.M. acknowledge funding by Grant-in-Aid for Scientific Research under Grant No. 15K13533 and 16H02206. Y.M. and Y.K. were also supported by JST CREST (JPMJCR18T2). The numerical simulations were performed on the JUWELS cluster at the Forschungszentrum Jülich.

I Lattice structure

The lattice vectors of the inversion-symmetric (8,3)c lattice are chosen as

[TABLE]

and the reciprocal lattice vectors become

[TABLE]

It has eight sites per unit cell that are located at

[TABLE]

II Pseudospin correlations

The pseudospin correlations defined in (3) can be used to identify the ground state manifolds of the different (an)isotropic regimes of the system.

We can introduce the unitary transformation Ishizuka and Motome (2013); Fujiki et al. (1984)

[TABLE]

mapping the triplet of 3-flux states that preserve the local divergence-free condition to be rotationally symmetric around the $W_{z}$-axis:

[TABLE]

It can be easily seen that $\|{\bf W}_{m}\|=\|{\bf P}_{m}\|=\sqrt{3}$ for all $m\in\{a,b,c\}$ and

[TABLE]

For large $T$, all eight possible 3-flux states are allowed, so we expect P to average to 0. For the low temperature regime of the isotropic system ($J_{x}=J_{y}=J_{z}$), the pseudospins ${\bf P}_{a}$, ${\bf P}_{b}$, ${\bf P}_{c}$ will be equally distributed in the system, giving the expectation value of the pseudospin correlations

[TABLE]

For the weak-$J_{z}$ limit, only the pseudospins ${\bf P}_{b}$, ${\bf P}_{c}$ are selected, resulting in

[TABLE]

while for the strong-$J_{z}$ limit, all plaquette triplets select the state ${\bf P}_{a}$, leading to

[TABLE]

Since we choose all our systems to have periodic boundary conditions in all spatial directions, each lattice with $N$ sites has $N_{p}=3N/4$ plaquettes and $N_{s}=N/4$ pseudospins. For a finite system, the pseudospin correlator results as

[TABLE]

III Sign-free quantum Monte Carlo simulations

Our large-scale quantum Monte Carlo simulations are performed in the parton basis of the spin system Kitaev (2006), in which they are sign-problem free Nasu et al. (2014). Technically, a single Monte Carlo step consists of a local update of the $\mathbb{Z}_{2}$ gauge degree of freedom on a random lattice bond, followed by an evaluation of the Majorana (free) energy required for the calculation of a Metropolis acceptance probability. While this latter step can be accomplished using exact diagonalization (ED) techniques Nasu et al. (2014), a far more efficient algorithm that avoids the cubic scaling of the ED approach is to employ the Green-Function-Based Kernel Polynomial Method (GF-KPM) Weiße (2009), which was recently introduced for Kitaev systems Mishchenko et al. (2017). Using this apporoach, the free energy change of the flip of a $\mathbb{Z}_{2}$ gauge variable is calculated from a set of four Green functions which can be approximated with a limited number of Chebyshev polynomials (typically with 128 to 256 coefficients in our case). The computational cost of each MC update is thereby reduced to $\mathcal{O}(N)$ scaling – a crucial step which enables us to simulate relatively large system sizes with up to 1728 spins.

Note that for the transformation of the spin model to the Majorana basis, we follow the local approach originally pioneered by Kitaev Kitaev (2006). In this local approach all lattice bonds carry gauge field degrees of freedom and can be addressed by a single-flip MC step, in contrast to an alternative scheme which makes use of a Jordan-Wigner transformation and where the resulting $\mathbb{Z}_{2}$ gauge field degrees of freedom are located only on the $z$-bonds Nasu et al. (2014). A crucial advantage of the local approach is that it allows us to simulate systems with periodic boundary conditions in all directions, thereby minimizing finite-size effects. There remains one subtlety, however: By replacing each spin with four Majorana fermions, the Kitaev ansatz artificially expands the Hilbert space of the system, which would normally require a rigorous projection into the physical subspace. However, it is known that physical and unphysical states for a given lattice and $\mathbb{Z}_{2}$ gauge field configuration are distinguished by their fermionic parity, reducing the effect of states of the expanded Hilbert space to deviations of order $1/N$, which can be neglected in the thermodynamic limit Zschocke and Vojta (2015). In practical terms, we have checked that already for small system sizes of $N=64$ sites, we cannot establish differences between the two approaches up to our numerical efficiency.

Finally, in order to prevent our MC simulation from freezing at low temperatures, we have applied a parallel tempering scheme Hukushima and Nemoto (1996) with one replica exchange proposed between every pair of neighboring temperature points after every MC sweep. This enabled us to perform simulations down to temperatures of $\mathcal{O}(10^{-3})$ in terms of the magnetic coupling strength. To calculate observables on the lower end of this limit and to interpolate between temperature points, we have used the Ferrenberg-Swendsen reweighting method with multiple histograms Ferrenberg and Swendsen (1988, 1989).