Partial flux ordering and thermal Majorana metals in (higher-order) spin liquids

In frustrated quantum magnetism, chiral spin liquids are a particularly intriguing subset of quantum spin liquids in which the fractionalized parton degrees of freedom form a Chern insulator. Here we study an exactly solvable spin-3/2 model which harbors not only chiral spin liquids but also spin liquids with higher-order parton band topology -- a trivial band insulator, a Chern insulator with gapless chiral edge modes, and a second-order topological insulator with gapless corner modes. With a focus on the thermodynamic precursors and thermal phase transitions associated with these distinct states, we employ numerically exact quantum Monte Carlo simulations to reveal a number of unconventional phenomena. This includes a heightened thermal stability of the ground state phases, the emergence of a partial flux ordering of the associated $\mathbb{Z}_2$ lattice gauge field, and the formation of a thermal Majorana metal regime extending over a broad temperature range.


I. INTRODUCTION
The emergence of topological phases from local constraints induced by competing interactions in frustrated quantum magnets has fascinated researchers for decades [1,2]. In a groundbreaking conceptual work, Kalmeyer and Laughlin [3] in the late 80s put forward the formation of a bosonic analogue of the fractional quantum Hall state in what has since been termed a chiral spin liquid. First envisioned as resonating valence bond (RVB) ground states of geometrically frustrated Heisenberg antiferromagnets and relevant to high-temperature superconductivity [4,5], such chiral spin liquids have remained elusive for many years. In the past decade, however, tremendous progress has been achieved in firmly establishing chiral spin liquids as ground states of microscopic models [6][7][8][9][10][11] and their eventual experimental observation in quantized thermal Hall measurements [12,13]. On the theoretical side, conceptual insight has been gained from the unambiguous numerical detection of chiral spin liquids via the calculation of modular matrices [14,15] from the ground state entanglement structure [16] of a number of kagome models [10,11]. Analytically, Kitaev showed that a time-reversal symmetry breaking magnetic field gives rise to a non-Abelian chiral spin liquid in his eponymous model [6]. The advent of Kitaev materials [17] has since produced a direct experimental observation of this state in measurements of a half-quantized thermal Hall effect in the spin-orbit entangled Mott insulator α-RuCl 3 [12].
The physical mechanism underlying the formation of a chiral spin liquid can be elegantly formulated using Wen's parton construction [18]: At low temperatures the original spin degrees of freedom fractionalize into a parton coupled to an emergent lattice gauge field [2]. For the Kitaev model, this parton decomposition is typically done [6] in terms of Majorana fermions coupled to a Z 2 gauge field [19] -thereby casting the original interacting spin model to a free (Majorana) fermion problem with a static Z 2 gauge order (at zero temperature). One can then resort to the tools of topological band theory to classify the band structure of the emergent Majorana fermions, and thereby classify the fundamental topological properties of the spin liquid state. For the case of the Kitaev honeycomb model, this perspective results in an understanding of the field-induced topological spin liquid as the formation of a Chern insulator (in the Majorana band structure) with all bands carrying a non-trivial Chern number |ν| = 1. Such a state -a Chern insulator of the emergent Majorana fermions -is gapped in the bulk, but exhibits chiral gapless modes at its boundary that are topologically protected, as depicted in Fig. 1(b). That such a Majorana Chern insulator indeed forms in the Kitaev material α-RuCl 3 is substantiated by the observation of a half-quantized thermal Hall conductance, which is direct evidence for Majorana fermions (and not conventional electrons) carrying the thermal edge currents. Subsequent measurements [13] of the anomalous character of the quantized thermal Hall effect, which arises even without any perpendicular magnetic field component, have brought verification that the observed state is indeed arising from the formation of topological Chern bands (and not the formation of Landau levels).
In parallel the field of topological band theory has been developed [20]. There one distinguishes strong topological insulators (TIs), whose topological features are protected by time reversal and/or charge conjugation symmetries, from crystalline TIs whose topology is endowed by certain lattice symmetries. A particular example of such crystalline TIs are second-order topological insulators, an instance of higherorder topology [21]. An n th order TI in d spatial dimensions exhibits (d − n)-dimensional topologically protected gapless modes that are localized at the intersection of n boundary planes, while the boundaries of codimension less than n remain fully gapped. A second-order TI in two spatial dimensions thus exhibits topologically protected corner modes, i.e. zero-dimensional gapless modes at the intersection of two boundaries giving rise to a corner. Such higher-order topology can also play out in the context of a quantum spin liquid -with the itinerant fractionalized parton degrees of freedom forming a non-trivial topological band structure. The example of a second-order spin liquid with topologically protected Majorana corner modes has been discussed in the context of an exactly solvable spin-3/2 generalization of the Kitaev model [22]. Within this same framework, the chiral spin liquid can be thought of as a first-order spin liquid, since a Chern insulator can be considered an example of a first-order topological insulator [23].
In this manuscript, we study the thermodynamic precursors and symmetry-breaking thermal phase transitions leading to the formation of a family of spin liquid ground states which arXiv:2007.04332v1 [cond-mat.str-el] 8 Jul 2020 FIG. 1. Variable-order spin liquids can be conceptualized via the formation of (topological) band structures of emergent partons. Shown here is the reflection of these band structures in the ground-state parton wavefunction |ψij| 2 for different settings with (a) trivial topology, (b) conventional (first-order) topology with gapless (chiral) edge modes, and (c) second-order topology with localized corner modes. The actual calculations have been performed in the Majorana representation of model (1) on a 16 × 16 lattice. exhibit the full range of parton band topology, second-order, conventional (first-order), and trivial topology, in a generalized Kitaev model. We employ sign-problem free quantum Monte Carlo simulations in the parton basis [24]. These numerically exact calculations allow us to track the fractionalization of the original spin degrees of freedom, the formation of gauge order, and the spontaneous breaking of time-reversal symmetry upon entering the different flavors of spin liquid ground states. Our main results include (i) the observation that the thermal stability of the chiral spin liquid in our model is enhanced by almost an order of magnitude in comparison with other chiral spin liquid models, with the highest transition temperatures reaching about 1/10 of the bare coupling strength; (ii) the emergence of partial flux order in an intermediate temperature range, accompanied by a characteristic 3-peak signature in the specific heat; and (iii) the formation of a gapless phase at finite temperatures that is best described as a thermal Majorana metal.
Our discussion of these results in the remainder of the manuscript is structured as follows. In Section II we briefly introduce a generalized spin-3/2 Kitaev model, its Γ-matrix representation, and the underlying five-coordinated Shastry-Sutherland lattice. In discussing its analytical solution at zero temperature, we also introduce the parton basis relevant to our sign-free QMC simulations to explore the thermodynamics at finite temperatures. The formation of a conventional (firstorder) chiral spin liquid is discussed in Section III. Our main results on thermal stability, partial flux ordering, and thermal metal formation are all discussed in detail here. In Section IV we then turn to the formation of a second-order spin liquid, whose zero-temperature properties we previously discussed in Ref. 22. Our focus here is on its thermodynamic properties. We conclude with an outlook in Section V.

II. THE SHASTRY-SUTHERLAND KITAEV MODEL
We start our discussion with a brief review of the generalization of the Kitaev model to the Shastry-Sutherland lattice [22,25], its fundamental (lattice) symmetries, the formation of spin liquid ground states of various levels of topology, and its numerical representation in sign-free quantum Monte Carlo simulations. A. The model: spin-3/2 and Gamma matrices The Kitaev honeycomb model is the paradigmatic example of an exactly solvable quantum spin liquid model. The model consists of spin-1/2 degrees of freedom on the sites of a honeycomb lattice interacting via bond-dependent Ising interactions. By representing the spin operators in terms of Majorana fermions, the model can be reduced to a nearestneighbor hopping model of non-interacting fermions coupled to a static Z 2 gauge field. The Kitaev model and its exact solution can be straightforwardly generalized to other lattices with an odd coordination number, z = 2n − 1, wherein the local "spins" are decomposed into 2n Majorana fermions.
Here, we study such a generalization of Kitaev's honeycomb model to the pentacoordinated Shastry-Sutherland lattice, previously introduced in Refs. [22,25] (see Fig. 2). The lattice is most well-known for the orthogonal dimer model, which was solved by Shastry and Sutherland [26] and serves as an effective low-temperature model for the transition metal oxide SrCu 2 (BO 3 ) 2 [27]. The generalized Kitaev model is described by the Hamiltonian where γ = 1, . . . 5 labels the bond direction and we have a set of five 4 × 4 anticommuting matrices Γ γ j for each site. Physically, the Γ-matrices can be interpreted as acting on either j = 3 2 spins or two coupled spin-1 2 degrees of freedom, such as spin and orbital degrees of freedom, on each lattice site [22].
To solve this model exactly, we represent the Γ-matrices on each site in terms of six Majorana operators c j and {b γ j } by setting [28] Γ γ j = ib γ j c j . The "bond Majoranas" b γ j are recombined into bond operatorsû jk = ib γ j b γ k .. Since allû jk commute with the Hamiltonian, they can be replaced by their eigenvalues u jk = ±1. We are thus left with a hopping model of Majorana fermions {c i } coupled to a static Z 2 -gauge field u jk , described by the Hamiltonian The spectrum of H must be invariant under Z 2 gauge transformations, and can thus depend only on the Z 2 fluxes W p degrees of freedom. The horizontal and vertical (blue) bonds carry a coupling J (which we set to unity), and the diagonal (red) bonds a coupling Jd. For a suitable choice of these parameters, the resulting higher-order spin liquid system is shown to possess a topologically non-trivial ground state. associated with plaquettes p, defined as where, as a convention, the product here is taken with a clockwise orientation. The Shastry-Sutherland lattice has two types of elementary plaquettes, viz, square plaquettes with flux W = ±1, and triangular plaquettes with flux W = ±i.

B. Sign-problem free sampling of flux configurations
To determine the ground state, one primary task is to find the flux configuration that minimizes the total energy. This is generally a non-trivial matter and one that can almost never be resolved in an analytically exact fashion -with the honeycomb Kitaev model being the most notable exception, for which a theorem by Lieb [29] can be invoked. For the general case one can, however, rely on a numerically exact treatment by performing quantum Monte Carlo (QMC) sampling of the gauge field, which in the Majorana decomposition introduced above is possible without encountering a sign problem [22,24] (and yet keeping track of all the Majorana physics).
In such a sign-problem-free QMC approach, the simulation is performed in the Majorana basis, where the Z 2 gauge degrees of freedom u jk can be sampled as classical Ising variables, with the Majorana quantum physics entering in the weights of the Markov chain sampling. For a fixed gauge configuration, the noninteracting Majorana Hamiltonian (Eq. (2)) can be numerically diagonalized, a computation which scales as N 3 , where N is the number of lattice sites [24,30]. All thermodynamic observables can then be obtained without needing to introduce an additional imaginary-time dimension or other non-trivial mapping schemes for the quantum problem. The QMC method is thus guaranteed to be sign-problemfree on any lattice geometry. A comprehensive discussion of the technicalities associated with this QMC approach to the Shastry-Sutherland Kitaev model can be found in Ref. [22].
Performing such a quantum Monte Carlo (QMC) simulation readily demonstrates that the ground state flux order corresponds to W = −1 for all square plaquettes. However, though W = W 2 , the value of W itself is not fixed. The two possibilities W = ±i are related by time-reversal symmetry (TRS) and degenerate in energy. At zero temperature, the ground state must thus spontaneously break TRS by selecting either W = i or W = −i for all triangular plaquettes. This spontaneous breaking of time-reversal symmetry is generally true for any Kitaev-type model on a lattice containing plaquettes with an odd number of bonds, as first noted by Kitaev [6], and later elucidated within a concrete spin-1/2 model by Yao and Kivelson [7]. This is precisely what is seen in our thermodynamic QMC data for finite temperatures, as we will discuss in detail below.

C. Symmetries and the ground-state phase diagram
On a conceptual level, the possible ground-state phases of the (Majorana) Hamiltonian can be inferred from the tenfold classification of topological insulators and superconductors [31,32]. The Majorana Hamiltonian, by design, obeys a particle-hole symmetry (PHS) which squares to +1, while time-reversal symmetry is broken spontaneously, as discussed earlier. Thus, this Majorana Hamiltonian belongs to symmetry class D. In two spatial dimensions, this symmetry class D allows for a Z invariant, viz, a Chern number, signaling the possibility of the formation of topological order. Such a topological state is precisely the chiral spin liquid (with a gapless chiral edge mode) discussed in the introduction.
If we further restrict to the case where all the "square lattice" bonds are of equal strength (which we set to J = 1), the only remaining parameter is J d , the strength of the diagonal "dimer" bonds of the lattice. As a function of varying J d the Majorana band structure exhibits two gapped phases: a topological phase with Chern number The two gapped phases are separated by a gap closing at (π, π), as illustrated in the lower panel of the schematic phase diagram in Fig. 3. The C = ±1 phase is a chiral spin liquid with a chiral Majorana edge mode and bulk Ising anyon topological order. On the other hand, the C = 0 phase, while still spontaneously breaking TRS, is fully gapped and possesses the same Abelian topological order as the toric code. We refer to this phase as a "trivial" spin liquid.

III. TRIVIAL AND CHIRAL SPIN LIQUIDS
Coming to the actual results of our finite-temperature analysis of the Shastry-Sutherland Kitaev model, we first concentrate on the thermodynamic behavior above the phase transition from the topological (first-order) spin liquid (with Chern number ν = ±1) to the trivial spin liquid (with Chern number ν = 0).

A. Thermodynamics
A central quantity distinguishing the various finitetemperature phases of our model is the specific heat C v (T ) and its characteristic multi-peak structure, as illustrated in Fig. 4. For the closely related Kitaev spin liquids, it is well established [24] that one finds two well-separated peaks in the specific heat -a smooth high-temperature peak, indicating a thermal crossover, corresponding to the (local) fractionalization of spins and a second low-temperature peak associated with the freezing of the Z 2 gauge field. We observe a similar two-peak structure of the specific heat in our model for a broad range of parameters J d 1 as well. Plotted in Fig. 4 is a characteristic C v (T ) trace in comparison with data for the Kitaev honeycomb model for similar system sizes. Both systems show a shallow high-temperature crossover around the value of the elementary coupling strength, in our case T ∼ 2, whose shape and height are essentially independent of the system size -indicating a purely local crossover. In fact, this is where the spin fractionalization happens and the system releases precisely half of its entropy [22,24]. A more pointed distinction is found in the low-temperature peak -this is a sharp peak for the model at hand (that sharpens with increasing system size), while it is a more shallow feature in the Ki- taev honeycomb model. This is an immediate reflection of the fact that in the model at hand there is a true phase transition occurring at this lower temperature -the spontaneous breaking of time-reversal symmetry upon entering the chiral spin liquid regime, while in contrast the Kitaev honeycomb model exhibits a finite-temperature crossover at this lower temperature scale at which, for a given system size, the Z 2 gauge field freezes into its ground-state configuration. While the latter is a true phase transition in three spatial dimensions [33], it remains a thermodynamic crossover in two spatial dimensions [34,35].
The key distinction between these model systems is that the Shastry-Sutherland lattice is non-bipartite and exhibits elementary (triangular) plaquettes with an odd number of bonds. For the emergent Majorana fermions this results in an ambiguous situation in which they can pick up a phase e ±iπ/2 upon hopping around such a triangular plaquette, endowing the latter with a flux W = ±i as discussed in the previous section. By spontaneously breaking time-reversal symmetry, one of the two possible signs is chosen, with the system simultaneously undergoing an Ising-type phase transition. Such a time-reversal symmetry breaking thermal phase transition was first observed in the Yao-Kivelson model (on a decorated honeycomb lattice) [36] at a temperature scale T c ∼ 10 −2 J. This is also the typical temperature scale for thermal phase transition in 3d Kitaev models. In comparison, the critical temperature scale of the model at hand is elevated, see Fig. 4, with its maximum close to one-tenth of the Kitaev coupling (at J d ∼ 1.2). One might speculate that this enhanced transition temperature is a reflection of the higher coordination number of the Shastry-Sutherland lattice (z = 5) in comparison to conventional Kitaev models on tricoordinated lattice geometries. This idea, however, does not hold up when further generalizing our model to a 7-coordinated lattice (by placing additional diagonal bonds on the lattice), which has a transition temperature of the same order of magnitude as the Shastry-Sutherland case.

B. Partial flux ordering
Upon closer inspection, the formation of flux order at low temperatures turns out to be slightly more intricate. As noted above, the spontaneous breaking of time-reversal symmetry is intimately connected with the flux ordering of the triangular plaquettes. With two such triangular plaquettes constituting a single square plaquette, this also implies ordering for the latter. The corresponding flux satisfies W = W 1 · W 2 = +1, independent of the actual assignment of the triangular plaque-FIG. 6. Partial flux ordering. For Jd 1, the square plaquettes assume an ordered π-flux configuration at T . In this regime, the triangle plaquettes remain disordered, assuming pairwise fluxes W = ±i (darker / lighter yellow). At Tc ∼ 0.1Jd, also the triangle plaquettes order into one of the two homogeneous flux configurations.
ttes [37]. Thus, an ordering of the 3-plaquettes implies ordering of at least half the 4-plaquettes. The converse is, however, not true; we can have ordering of the 4-plaquettes while the 3-plaquettes remain disordered, resulting in a partial flux ordering. This is exactly what we observe for a limited parameter range 0 < J d 1, where the 4-plaquette fluxes order at a higher temperature T than the critical temperature T c for the 3-plaquettes fluxes. This is illustrated in Fig. 5 where we plot the phase diagram of our model by color-coding the flux of the 3-plaquettes (top panel) and 4-plaquettes (bottom panel) as function of J d . A schematic rendering of the indermediate partial flux ordering is provided in Fig. 6.
The evolution of this partial flux ordering with varying coupling strengths can be seen in the sequence of data sets for vertical cuts through the phase diagram provided in Fig. 7. The top row shows the gauge contribution to the specific heat (i.e. omitting the Majorana contribution resulting in the higher temperature crossover). In the partial flux ordering regime 0 < J d 1, one finds that the low-temperature specific heat peak actually splits into two parts -with the lower temperature peak indicating the true thermal phase transition associated with time-reversal symmetry breaking and 3-plaquette flux ordering, see also the medium row of panels. While this transition quickly moves to zero temperature as J d → 0, there remains a signature in the specific heat around O(1/10), which upon closer inspection is the thermal crossover associated with the ordering of the remaining half of the 4-plaquettes, see the lower row of panels.
Taking a step back, we conclude that the emergence of a partial flux ordering is, for a limited range of parameters, a precursor phenomenon to the formation of a low-temperature topological chiral spin liquid.

C. Thermal Majorana metal
The thermodynamic signatures discussed so far -the finitetemperature, time-reversal symmetry breaking phase transition as well as the emergence of a partial flux ordering and its associated thermal crossover -are both closely connected to the underlying lattice gauge theory. From a more conceptual perspective, these aspects are an interesting variation to the gauge physics that has been intensely studied in the context of two-and three-dimensional Kitaev models [24,33]. We now turn to additional thermodynamic aspects of our model at hand, which are genuinely rooted in the physics of the emergent Majorana fermions. The most notable feature here is the formation of a thermal metal regime above the transition to the topological chiral spin liquid as illustrated in the schematic phase diagram of Fig. 3, and as evidenced in the numerical observations of Fig. 8. This thermal metal regime, principally located in the intermediate temperature regime T c T T (i.e. between spin fractionalization and the time-reversal symmetry breaking transition), does not extend over the entire parameter space of our model, but marks a relatively sharp transition for a critical value of J d that is considerably shifted in comparison with the zero-temperature transition between the chiral and trivial spin liquids. We rationalize these numerical observations using analytical arguments based on the analytical (self-consistent) Born and T-matrix approximations along with a numerical computation using transfer matrices. We thus explain how the entire gapless thermal metal phase emanates from the quantum critical point between the chiral and trivial spin liquid at zero temperature.

Numerical observations
In our numerics, the existence of two distinct regimes above the low-temperature chiral spin liquid phases is most evident in calculations of the average Chern number |ν| of the Majorana band structures encountered when sampling (and averaging over) the different gauge configurations for a given temperature. The behavior of this average Chern number as a function of J d is plotted in Fig. 8(a), which immediately reveals three distinct regimes: For a horizontal cut in the lowest temperature regime, the average Chern number jumps from |ν| = 0 in the trivial phase to |ν| = 1 in the topological phase, as expected. Unanticipated is probably the behavior for a horizontal cut in the intermediate-temperature regime, where our simulations also reveal a relatively sharp (vertical) boundary at J d ≈ 1.9, at which the average Chern number mimics the low-temperature behavior, jumping from |ν| = 0 for large couplings J d 1.9, to |ν| ≈ 0.6 for J d 1.9 (see also the scans of |ν| shown in Figs. 19 and 20 of the Appendix). This clearly separates two different thermodynamic regimes, which seemingly persist not only up to the fractionalization crossover scale but beyond into the paramagnetic regime extending all the way up to infinite temperatures.
One important aspect in interpreting this numerical result, in particular the non-quantization of the average Chern num- ber for J d 1.9, is to ask whether the Majorna band structure actually remains gapped -a prerequisite for the proper calculation of a Chern number -also at finite temperatures. In fact, this is precisely what distinguishes the two regimes: while the system remains gapped for J d 1.9, it becomes gapless for smaller couplings. That the phase for J d 1.9 is in fact gapless can be shown in a straightforward manner by computing the gap in the Majorana spectrum averaged over flux configurations sampled with a uniform distribution, or equivalently, at infinite temperature. As shown in the top panel of the schematic phase diagram of Fig. 3, the gap in the Majorana spectrum indeed vanishes for small J d and slowly opens only for J d 1.9, as compared to the scenario of a gap closing and re-opening for J d = 2 √ 2 ≈ 2.8 in the zero-temperature case. The Chern number results are thus only well-defined in the gapped regimes of the phase diagram, J d 1.9 and T < T c . For detailed scans of the average Chern number we refer to Appendix A.
Having established the principal gapless character of the finite-temperature spectrum for J d 1.9, a more meaningful quantity to calculate is the effective number of states accessible to the system at a given temperature where ρ(E) is the Majorana density of states and n F (E, T ) is the Fermi function, whose derivative has a peak at E = 0 of width ∼ T . Plotting N 0 (T ) as a function of J d for various temperatures in Fig. 9, we indeed see that the number of available state becomes finite for small J d , indicating a gapless state, and vanishes only near the transition point J d ≈ 1.9, indicating a truly gapped state for larger J d . We note in passing that the partial flux ordered phase also appears to be gapped, as N 0 (T ) → 0 in the corresponding parameter regime, i.e., for J d 0.3 and T 0.04. This gap opens at the crossover to the thermal metal phase, as indicated by the increase in N 0 (T ). This is also indicated in the Chern number plot of Fig. 8(a), where a narrow blue stripe above the partially flux-ordered phase indicates the thermal crossover points associated with the ordering of square plaquettes. Thus, for J d < 1, the Kitaev Shastry-Sutherland system starts in the gapped chiral spin liquid ground state at T = 0, undergoes a thermal phase transition into a gapped partial flux-ordered phase, and the gap finally collapses at the crossover to the intermediate temperature phase.

Symmetry-class considerations
To interpret these numerical observations, it is useful to remind oneself of the symmetry class classification of free fermion systems [38,39] and its application to topological states of matter [31,32]. When considering the Majorana fermion representation of our model, the system naturally exhibits particle-hole symmetry (due to the self-conjugation of Majorana fermions) and broken time-reversal symmetry, which is implied by the non-bipartite geometry of the lattice as discussed earlier. This puts our model into symmetry class D. In the 10-fold way classification of topological states of matter [31,32], this symmetry class allows, in two spatial dimensions, for the formation of both a trivial and a topological insulator (with vanishing and integer Chern number, respectively). In addition, it is well established that in the presence of disorder this symmetry class allows for the formation of a gapless "thermal metal" phase [40][41][42][43][44][45].
In the Kitaev Shastry-Sutherland model, all three of these phases are found to be realized. In the low-temperature regime, where time-reversal symmetry is spontaneously broken and the gauge field is ordered, we realize the two "clean" Going to the intermediate-temperature regime, i.e. above the ordering transition at T c , a more subtle situation arises: Here the Z 2 gauge field does not exhibit any order, effectively providing a static disorder potential for the Majorana system. And while time-reversal symmetry (TRS) is not broken on a global level, which strictly speaking does not put the system into symmetry class D anymore, TRS is broken for any given static disorder configuration: Each particular flux configuration must have a fixed value of flux (±i) across the 3-plaquettes, and hence breaks time-reversal symmetry. The ensemble of Majorana Hamiltonians thus always corresponds to the symmetry class D at any temperature, even if the true state of the Kitaev model is itself time-reversal invariant. This allows the system to principally form the thermal metal phase of symmetry class D.
This thermal metal phase is indeed realized in the intermediate temperature regime for J d 1.9, i.e. in the parameter regime that our numerical experiments indicated as gapless. The very nature of the thermal metal phase can be probed explicitly via a calculation of the low-energy density of states ρ(E) of the Majorana fermions, which is supposed to show a characteristic "ringing" [40]. For E → 0, the latter takes the universal form with a single fitting parameter α. This is indeed what we find upon numerical inspection, as shown in Fig. 10 for a parameter set deep in the thermal metal regime (J d = 1.5 and T = 1.86). Similar behavior has first been reported for the honeycomb Kitaev model in a magnetic field in a broad temperature regime above the thermal crossover into its chiral spin liquid ground state [46].

Analytical perspective
The essential physics of the formation of the thermal metal phase can be understood analytically by analyzing our Majorana hopping model with uncorrelated Z 2 flux disorder [47]. Explicitly, this implies that arbitrary disorder configurations can be obtained by starting from the ground state flux configuration and flipping each plaquette flux with a probability p, independent of all other fluxes. This disorder density p then roughly corresponds to the temperature in the full Kitaev model. In particular, p = 0 is the ground state (T = 0), while p = 1/2 is the state with completely random fluxes (T = ∞).
We analytically study the effect of this flux disorder by computing the self-energy as a power series in the disorder density -which corresponds to a series in the number of scattering events -using the self-consistent Born approximation (SCBA) as well as the self-consistent T-matrix approximation (SCTA). We find that the net effect of the disorder is a renormalization of the bare coupling strength J (which we have set to 1 so far), so that the new phase boundary occurs at where k 0 = (π, π) is the gap closing point at J d = 2 √ 2 in the clean limit (see Appendix B for a detailed derivation). In Fig. 11, we plot this phase boundary as a function of the disorder density. We find that for p 0.25, the disorder-induced renormalization of J leads to a continuous shift of the phase boundary between the trivial and topological gapped phases towards smaller J d , i.e., to a net suppression of the chiral spinliquid phase [48].
The most intriguing result of this computation is that the entire thermal metal phase can be thought of as emerging from the renormalization of the quantum critical point at J d = 2 √ 2. This can be understood by a simple scaling argument: The leading term in the perturbation series for Σ p (0, k 0 ) corresponding to scattering events from n impurities scales as 4pJ 2 n . Thus, for small p, we neglect the n > 1 terms to write J eff as a linear function of p. Within this "noncrossing" approximation, Eq. (6) has a single solution for a given p, leading to a single gapless point for each disorder density which we can approximate using the self-consistent Born (SCBA) and T-matrix approximations (SCTA). On the other hand, for p 1/4, we can no longer ignore the higherorder diagrams, so that Eq. (6) has, in fact, an infinite number of solutions, accounting for an entire gapless region in the phase space, which is precisely the thermal metal phase.
Since the perturbative approach is no longer useful in the high disorder density regime, we complement it with a numerical computation of the transmission coefficient. In particular, we compute the disorder-averaged generalized transfer matrix [49] for the Majorana model on a cylinder geom- etry with L x L y . The eigenvalues of this transfer matrix are related to the inverse localization lengths ξ i , which can be used to compute the transmission coefficient as g = 2Ly i=1 sech 2 (L x /ξ i ). This can be thought of as the effective number of gapless channels, and is proportional to conductance for complex fermions [50].
This approach effectively captures the contributions from all n-impurity scattering events, and thus yields a reliable estimate of the phase boundary also for large disorder densities. In Fig. 11 we plot g as a function of J d and p, overlaid with the phase boundaries obtained from SCBA/SCTA. We see a qualitative agreement between the numerical and analytical computations in the low density regime p 1/4, verifying the effective renormalization of J. For p 1/4, we clearly see the spreading of the phase boundary into a gapless region extending up to J d ≈ 1.8, which is also in good agreement with the results obtained using QMC for the full Kitaev model (Fig. 8).

Physical interpretation
In closing our discussion of the thermal metal, we need to address the question whether this thermal metal is a legitimate physical phase of our spin model or a mere artefact of our calculations relying on a Majorana decomposition. The relevance of this question acutely presents itself in the observation that the putative thermal metal regime seems to possibly extend beyond the intermediate-temperature range all the way up to infinite temperatures, as possibly suggested by the calculation of the average Chern number plotted in Fig. 8.
The answer to this question is two-fold. It sensitively depends on whether the Majorana fermions correspond to actual physical degrees of freedom of the system (after spin fractionalization) or whether they remain mere mathematical objects that can always be invoked in a spin decomposition. This fine distinction is well known from the solution of the honeycomb Kitaev model [6], whose intricacy lies precisely in the fact that the spin decomposition "becomes real" in the lowtemperature phase and describes the emergent fractionalized degrees of freedom. It is also exactly this distinction which distinguishes the intermediate-temperature regime from the high-temperature paramagnet in our model. The intermediatetemperature regime is defined as the temperature regime sandwiched between the fractionalization crossover at T and the low-temperature gauge ordering transition at T c , in which the emergent fractionalized degrees of freedom of itinerant Majorana fermions moving in the background of a static, but still disordered Z 2 gauge field, form. In this regime, our line of arguments outlined above fully applies and we conclude that in this temperature regime the thermal metal is a true physical phase.
The high-temperature paramagnet, on the other hand, is different. Here the Majorana decomposition of the spin operators is a mathematical possibility, but there is no deeper physical meaning associated with it (or any other spin decomposition using, e.g., complex fermions or other parton degrees of freedom). As such our numerical observation of a finite average Chern number and other class D physics in the Majorana representation is a direct reflection of the choice of decomposition, but not physically relevant.
We return to the question of how the thermal metal regime in the intermediate-temperature regime can be probed in the original spin system in our discussion section at the end of the manuscript.

IV. SECOND-ORDER SPIN LIQUID
The physics of the Kitaev Shastry-Sutherland model becomes even richer when considering a staggering of the plaquette couplings, which has been shown [22] to induce another variant of spin liquid ground state. This "second-order" spin liquid (SOSL) exhibits a Majorana band structure in the ground state that, akin to a second-order topological insulator [21], is gapped in the bulk but exhibits gapless corner modes (i.e., d − 2 dimensional zero modes, as opposed to the usual d − 1 dimensional zero energy modes). These corner modes are a manifestation of the formation of a symmetry-enriched topological order that is protected by two mirror symmetries of the lattice (indicated in Fig. 12).
In the following, we will turn to the thermodynamics accompanying the formation of such a second-order spin liquid. In doing so, we will concentrate on a representative choice of coupling parameters deep in the SOSL regime. Specifically, we introduce the plaquette staggering δJ = 0.7 (resulting in weakly and strongly coupled plaquettes with coupling strength J −δJ and J +δJ, respectively) and vary the relative coupling of the diagonal versus plaquette bonds J d /J. At zero FIG. 12. Kitaev Shastry-Sutherland mode with staggered plaquette couplings. The dashed/solid lines correspond to staggered couplings J ± δJ, respectively, which generates a hierarchy among the square plaquettes: A quarter of the squares possess either only weak (green) or strong (blue) couplings on the edge bonds, while the remaining half of plaquettes has two weak and two strong bonds (cyan). The dotted gray lines denote the two mirror axes. temperature, the phase diagram again consists of two phases [22]: a chiral spin liquid for J d /J > 2 √ 2 and a second-order spin liquid phase for J d /J < 2 √ 2.

A. Thermodynamics
The principle thermodynamic signatures above the formation of this SOSL are similar to what we have seen for the more conventional spin liquids discussed in the previous section: The specific heat C v (T ) exhibits a characteristic multipeak structure, with a high-temperature local crossover at T ∼ 2J d in this case, and a sharp low-temperature peak, at T c , marking a true thermal phase transition associated with spontaneous breaking of time-reversal symmetry. For J d /J 6.7 the lower peak splits into two peaks, the lower of which marks the symmetry-breaking transition and the upper one, at T , marking a crossover into another partial flux ordered state. This is summarized in the finite-temperature cuts of Fig. 13, which clearly indicate the multi-peak structure of the specific heat, and the thermal phase diagrams of Fig. 14, which shows the extent and boundaries of the different regimes.
The position of the low-T transition in temperature space monotonously decreases if J d /J is decreased. In the chiral spin liquid phase, it shows the particularly high value T c ∼ 0.1J d for large J d /J, a phenomenon which is discussed above. For J d /J → J c = 2 √ 2, the transition temperature reaches the order of magnitude T c ∼ 10 −2 J d , and, in the SOSL phase, it is further lowered to T c → 10 −3 J d . We note that for J d /J < 1.8, the transition temperature T c moves below the temperature range of our QMC simulations. In this limit, where the coupling J − δJ on half of the lattice bonds approaches 0, it is expected that the transition temperature rapidly decreases to lower temperature scales [22]. . The specific heat Cv shows a three-peak structure (a). While the high-temperature crossover is associated with spin fractionalization and the low-T phase transition with spontaneous breaking of time-reversal symmetry, the intermediate crossover indicates a partial flux-ordering of square plaquettes (b), which is a consequence of staggered bond couplings J ± δJ. This choice of coupling generates a hierarchy between the square plaquettes of the lattice, which results in different ordering temperature scales.

B. Partial flux ordering
The emergence of a partial flux ordering in an intermediate temperature range is again a precursor phenomenon for the formation of spin liquid ground states that, in the presence of staggered plaquette couplings, comes in even more variations. This is readily illustrated by measurements of the three-and four-plaquette fluxes, overlaid as color-coding in the thermal phase diagrams of Fig. 14, which reveal a partially flux ordered regime in the parameter range J d /J 6.7 and above the thermal phase transition T > T > T c . However, in this case, the pattern of partial ordering is very different to the one previously discussed in Section III for the original model. While there, we observed a region in which the square plaquettes were fully ordered, and the triangular plaquettes remained disordered, here the partial flux order is characterized by partial order of the square plaquettes, with the triangular plaquettes again remaining disordered. This can be clearly seen in Here, the partial flux ordering is determined by the staggering of bond couplings (dashed / solid lines).At T < Tc, all square plaquettes have a π-flux (green). At Tc, the plaquettes with only weak edge couplings J − δJ are the first to become disordered, with πand 0-fluxes (green / red). Within this partially flux-ordered phase at Tc < T < T , also the squares with two weak and two strong bond couplings become disordered. Finally, at T , also the square plaquettes with strong couplings disorder.
the thermal phase diagram with Fig. 14(a) showing the threeplaquette flux |W | = 0 and Fig. 14(b) showing the squareplaquette flux −0.75 W −0.25 within the partial flux ordered regime.
What is the reason for this new behavior? When looking at the model with staggered couplings, J ± δJ, we see that the lattice is now composed of three different kinds of square plaquettes: (i) one quarter are "strong plaquettes", which contain a diagonal J d bond and four "strong" bonds with coupling J + δJ, (ii) one quarter are "weak" plaquettes, which contain a diagonal J d bond and four "weak" bonds with coupling J − δJ, (iii) while the remaining half are "mixed" plaquettes, which do not contain any diagonal bond and are made up of two "strong" bonds and "two" weak bonds. This hierarchy of couplings, and thus vison gaps, for the three different kinds of square plaquettes is precisely what underlies the emergence of the partial flux ordering.
If we consider a lattice of N p total square plaquettes, the behavior seen within the partially flux-ordered regime can be explained as follows. Starting from the lowest temperatures, for 0 < T < T c , all square plaquettes are in an ordered πflux phase, with W = −1 for all N p square plaquettes. At T c the triangular plaquettes become disordered, triggering the recovery of time-reversal symmetry, and, also at T c , the N p /4 "weak" square plaquettes similarly become disordered. This loss of N p /4 plaquettes explains the drop of W from −1 to −3/4 at T c , shown in the lower panel of Fig. 13. Further increasing the temperature within the regime T c < T < T , the N p /2 "mixed" plaquettes gradually disorder for higher temperatures, resulting in a smooth change of W from −3/4 just above T c to −1/4 just below T , see again the lower panel of Fig. 13. Finally, at T , the remaining N p /4 "strong" plaquettes disorder, resulting in a fully disordered flux state with W = 0 for all plaquettes for T > T .

C. Thermal Majorana metal
Turning to signatures of Majorana physics in the thermodynamic behavior, we find evidence for the formation of a thermal metal regime also for the staggered model, similar to what we discussed extensively for the original Kitaev Shastry-Sutherland model in Sec. III C above. The numerical evidence for the emergence of such a phase again is the observation of the vanishing of the Majorana gap, which is reflected in the finite (but not quantized) average Chern number illustrated in Fig. 16. The broad temperature regime where this average Chern number does not vanish (beyond the topological ground state phases) sharply sets in above the (partially) flux ordered phases, i.e. in the regime where the Z 2 fluxes are effectively disordered and thereby create a disorder potential for the Majorana fermions (which form a nearly gapless band structure). Like in the original model, the average Chern number remains finite also in the high-temperature paramagnet and we refer to our previous discussion on the physical interpretation of this observation at the end of section III C.

V. SUMMARY
One of the most intriguing phenomena associated with the formation of quantum spin liquids is the emergence of novel, FIG. 16. Thermal metal. The color coding shows the average Chern number |ν| across the thermal phase diagram of the staggered Kitaev Shastry-Sutherland model. For the low-temperature groundstate phases the Chern number is |ν| = 1, in the topological spin liquid phase and, as expected, |ν| = 0 in the second-order spin liquid (SOSL). In the intermediate temperature regime, we find that the average Chern number vanishes in the partially flux-ordered phase in the temperature range Tc < T < T for the respective range of couplings. Above T (and Tc for large Jd, respectively) the average Chern number is finite, but not quantized around |ν| ≈ 0.8. The black data points in indicate the respective transition temperatures of the system as a function of the coupling ratio Jd/J, extracted from specific heat traces akin to the ones shown in Fig. 13. fractionalized quantum mechanical degrees of freedom -a quasiparticle, generally referred to as a parton, coupled to the gauge field of a deconfined lattice gauge theory. For the spin-3/2 Kitaev Shastry-Sutherland model studied in this manuscript, these emergent degrees of freedom are itinerant Majorana fermions and a Z 2 lattice gauge field, with all spin liquids coming in the form of varying levels of topology -trivial, first-and second-order, akin to the classification of higher-order topological insulators.
With a focus on the thermodynamic signatures of this fractionalization, we have observed in our numerically exact (sign-free) quantum Monte Carlo simulations characteristic fingerprints of both the underlying gauge physics and the emergent Majoranas. Below the thermodynamic crossover where these fractionalized degrees of freedom (locally) come to live, the gauge physics manifests itself through a sequence of (partial) flux ordering transitions, with the transition into the ground-state manifold being accompanied by (global) time-reversal symmetry breaking. The latter is also mandated (already on a local level) by the Majorana physics for nonbipartite lattice geometries.
Another striking manifestation of the interplay of Majorana physics and the Z 2 gauge structure is the formation of a thermal metal regime in an intermediate temperature range where the gauge field is intrinsically disordered. Notably, this gapless phase emerges only for a limited parameter regime of the Kitaev Shastry-Sutherland model, with a sharp transition to a gapped regime. This principle distinction might also make the thermal metal phase observable in experimental studies, e.g. in heat transport measurements which have been demonstrated, both theoretically [51][52][53] and experimentally [12,13,54,55], to be sensitive probes of the Majorana physics. The gapless versus gapped character of the intermediate-temperature phases should also reflect itself in four-spin correlation functions that probe the algebraic versus exponential decay of the bond-energy bond-energy correlations.