Subdiffusive dynamics and critical quantum correlations in a disorder-free localized Kitaev honeycomb model out of equilibrium

Disorder-free localization has recently emerged as a mechanism for ergodicity breaking in homogeneous lattice gauge theories. In this work we show that this mechanism can lead to unconventional states of quantum matter as the absence of thermalization lifts constraints imposed by equilibrium statistical physics. We study a Kitaev honeycomb model in a skew magnetic field subject to a quantum quench from a fully polarized initial product state and observe nonergodic dynamics as a consequence of disorder-free localization. We find that the system exhibits a subballistic power-law entanglement growth and quantum correlation spreading, which is otherwise typically associated with thermalizing systems. In the asymptotic steady state the Kitaev model develops volume-law entanglement and power-law decaying dimer quantum correlations even at a finite energy density. Our work sheds light onto the potential for disorder-free localized lattice gauge theories to realize quantum states in two dimensions with properties beyond what is possible in an equilibrium context.


I. INTRODUCTION
It is the general expectation that realistic isolated quantum many-body systems driven out of equilibrium eventually thermalize such that the relaxed long-time steady states become locally indistinguishable from thermal ensembles [1][2][3][4][5][6]. Two types of celebrated exceptions beyond this paradigm are quantum integrable models [7][8][9][10][11] and the Anderson or many-body localization (MBL) mechanism imposed by strong disorder [12,13]. In two dimensions, the exploration of ergodicity breaking dynamics in interacting systems remains a challenge especially in view of the argued instability of MBL in two dimensions [14]. Recent years have witnessed a new type of mechanism for nonergodic dynamics unique to lattice gauge theories where static local gauge charge or flux serves as a source for an effective internal disorder [15][16][17]. Importantly, this so-called disorder-free localization scenario does not rely on breaking translational invariance and can even occur in interacting two dimensional models [18,19], opening up a promising route targeting the challenge of realizing quantum states in two dimensional nonergodic systems with properties beyond any equilibrium counterpart.
In this work we show that the Kitaev honeycomb model driven to highly-excited states by a nonequilibrium quantum quench enters a peculiar disorder-free localized phase exhibiting subdiffusive dynamics towards a critical state exhibiting an algebraically decaying dimer quantum correlation function. Specifically, we investigate the nonequilibrium dynamics in the Kitaev honeycomb model in a weak skew magnetic field starting from a spin polarized initial state. The problem can be mapped to a weakly interacting Majorana fermion model coupled to a static Z 2 gauge field [20], which for the considered dynamics becomes effectively disordered. * timeexplorer1991@gmail.com Although a number of previous works have considered the intertwined physics between fermion and flux in the Kitaev model [21][22][23][24][25][26][27][28][29], the central open question has remained as to whether this model can break ergodicity and can potentially host non-thermal quantum order. In the noninteracting limit, we find that the gauge flux disorder localizes most of the Majorana fermions but fails to freeze the metallic and critical modes, leading to the observed subdiffusive dynamics although the system is overall nonergodic [30][31][32][33][34]. We identify the subdiffusive dynamics in both an algebraic spread of quantum correlations and the power-law growth of entanglement. At late times, the system relaxes to a steady state with dimer quantum correlation functions decaying algebraically in space, which is characteristic of quasi-long range order not accessible in thermal equilibrium. We argue that this quasi-long range order implies a divergent multipartite entanglement as quantified by the quantum Fisher information. We find evidence that our main findings are robust against the leading order perturbative Majorana fermion interactions induced by the skew magnetic field according to our numerical calculations for up to 128 spins on long time scales. Our results can be extended to any Z 2 lattice gauge theory coupled to chiral Majorana fermions as long as the gauge flux can be considered static and disordered on the considered time scales.

II. MODEL
The Kitaev model consists of spin- 1 2 degrees of freedom on the honeycomb lattice with spin-orbital locking Ising interactions D µ j = −σ µ j σ µ j+eµ , where j labels spin site and e µ denotes nearest neighbour vector of different orientations µ = x, y, z (Fig. 1a). In the presence of a weak [111] skew magnetic field the Hamiltonian iŝ For h = 0 the productŴ ≡ D x D y D z D x D y D z surrounding a hexagon plaquette commutes withĤ K , implying an extensive number of local integrals of motion. In the targeted limit h ≪ J µ we take into account the magnetic field perturbatively to the leading order that preserves these local symmetries [20]: whereh ∝ h 3 J 2 . The perturbative interaction acts on three spins that live on any wedges ∧, or the end of any Y junction. Either by introducing the gauge redundancy [20] and fixing the gauge, or by a Jordan-Wigner transformation [35,36], one can mapĤ onto an interacting Majorana fermion minimally coupled with Z 2 gauge field on the links: H = ⟨j→l⟩ J µ iu j,l β j α l +h ⟪j→l⟫ iu j,k u k,l (α j α l + β j β l ) +h Y u i,j u i,k u i,l (β i α j α k α l − α i β j β k β l ) , where α(β) denotes Majorana fermion on A(B) sublattice marked with open(close) circle in Fig. 1a. The static gauge field on link is pinned to u j,j−e x(y) = 1, u j,j−ez = ±1. W on plaquettes are transformed to be locally conserved Z 2 gauge fluxes. The last four fermion term is a chiral and gauged Majorana Hubbard interaction [37], where j, k, l are arranged in a counter-clockwise order around i.

III. QUANTUM QUENCH PROTOCOL
We prepare a simple initial state as a Néel state such that σ z Ψ 0 ⟩ = ± Ψ 0 ⟩ on A B sublattice respectively, which is to be evolved byĤ later on. In the fermion representation the initial state becomes a gauged fermion vacuum coupled to a disordered gauge field background: where N is the number of unit cells(z-links), and the Fock state satisfies iu j,j−ez α j β j−ez ψ {u} ⟩ = ψ {u} ⟩. In the sector ∏ j σ z j = 1 the anti-periodic boundary condition in spin Hamiltonian is mapped to periodic boundary in Majorana Hamiltonian. In this main-text we'll mainly focus on the isotropic coupling J x = J y = J z ≡ J,h = 0.25J.
For observables that preserve the gauge fieldÔ = ∑ {u} O {u} [15,16], where the average over gauge-field configurations can be performed via Monte-Carlo sampling. The typical {u} configuration is random, making the dynamical problem equivalent to Majorana fermions subject to Z 2 gauge (π) flux disorder, although our model is overall translational invariant [15,16].
Overall we target the description of the nonequilibrium dynamics through a sequence of two steps. First, we will study in detail the exact solvable point i.e. the noninteracting limit of Eq. (3), where we find that the system becomes nonergodic due to disorder-free localization, and afterwards explore the influence of interactions.

IV. EXACTLY SOLVABLE POINT
When the Majorana interactions are neglected, the model becomes exactly solvable. For each gauge configuration the dynamics is governed by a free Majorana fermion Gaussian Hamiltonian that can be computed efficiently. By randomly sampling gauge fields u = ±1 on the z-links, we compute the real-time evolution of various physical observables that are natural in both the spin and the fermion language.
First, we consider the spin dimer expectation values D µ j = iu j,j−eµ α j β j−eµ (µ = z, x), which relax exponentially fast to the same constant loosing the memory of initial anisotropy (see Fig. 1b). The observables along x and y directions can be related by mirror symmetry. As we will show, the dimer quantum correlation functions j+r Ψ(t)⟩ c , exhibit a much slower and intricate dynamics, which quantifies the correlation of gauged local fermion parity. Remarkably, we find that C zz(xx) (r, t) exhibits an algebraic light-cone in spacetime, see Fig. 1(c), with the the wave-front following a power-law Jt ∝ r z behavior. The dynamical exponent is obtained as z = 2.5(2) for C zz while z = 2.7(3) for C xx . Notice that the dynamical exponent z here is associated with information transport instead of particle or energy transport, and z > 1 signals subdiffusion [38]. In Fig. 1d we further corroborate this by achieving a data collapse upon rescaling the time axis Jt r z . While such subballistic behavior in systems with conventional disorder is typically observed on the ergodic side close to the MBL transition lying between diffusive and glassy limit [30-34], here we observe such dynamics for a disorder-free localized model, as we will argue in more detail below.
At long times the system settles to a steady state, which, as we find, is of nonergodic critical nature with correlations decaying algebraically in space, as seen in Fig. 2. We observe that the decay of C zz(xx) (r, t) is consistent with a power-law in space, whose exponent increases for larger system sizes and appears to converge near 0.63 (2). The power-law decaying correlation function in all directions x, y, z are reminiscent of the Kosterlitz-Thouless phase with quasi-long range order [39], without spontaneously breaking the spin-orbital three-fold rotation symmetry. However, even when this symmetry is explicitly broken in the anisotropic regime, we still observe critical correlation [40]. It is the effective disorder that partially inhibits the finite energy density fluctuations and stabilizes the quasi-long range spin dimer order [41].
These critical quantum correlations further have an immediate impact onto the entanglement content of the reached steady state, as the dimer quantum correlation function can be directly linked to a quantum Fisher infor- mation density via f ∼ N 1−∆ diverges in the thermodynamic limit. As a consequence the steady state exhibits strong multipartite entanglement.
For a more detailed quantification of the entanglement properties we further consider the projective bipartite entanglement entropy that serves as an entanglement diagnostics for quantum disentangled liquids [45][46][47]. Namely, we measure the von Neumann entanglement entropy for half of the Majorana fermions, when the gauge field is projected onto the diagonal ensemble: with the reduced density matrix  Fig. 3, at early time the entanglement grows with an area law. At a second stage, the entanglement entropy exhibits a further growth according to a subballistic power-law S ∝ t z . From a fit to the data we obtain the entanglement dynamical exponent z = 2.4(1), which within the accuracy of our simulations aligns with the exponent appearing for the subballistic spreading in C zz(xx) (r, t). In a system of finite size, we find that the entanglement entropy saturates to a volume-law state S v ∝ L x L y typical for the highly excited free fermion states [53,54], as shown in inset of Fig. 3. These numerical findings again highlight the unconventional nonequilibrium dynamics that we observe in the disorder-free localized Kitaev model.
One may ask what if we deform the initial state. By tuning the gauge flux density in the initial state by applying an operator ∏ q ( 1 2 +( 1 2 −p)Ŵ q ), the exponent ∆ as well as z change continuously, as visible in Fig. 4, which corroborates the robustness of the critical dynamical phase reminiscent of Kosterlitz-Thouless phase. However, notice that our critical dynamical phase at late-time steady state should be contrasted with the one exhibiting critical initial slip in short-time relaxation [55].

V. LOCALIZATION ANALYSIS
The peculiar coexistence of subdiffusive dynamics at transient time and the quasi-long-range order at late time implies a subtle localization scenario in behind. Indeed we find a mixture of localized and critical modes from the standard numerical diagnostics [40] including level spacing statistics [56, 57], localization length in two dimensions [58-60], and Chern number [20,[61][62][63].
The localization length is calculated by the retarded Green's function using iterative Dyson's equation for a semi-infinite quasi-one-dimensional geometry, followed by a one-parameter-scaling-collapse for varying narrow width. As shown in Fig. 5a, the localization length and the level spacing ratio for a finite size system are consistent in showing three energy windows with delocalization tendency. The delocalization at zero energy was known to be responsible for a low-energy Majorana thermal metal state in the class D dirty superconductors [64][65][66] or Majorana lattice model [67][68][69], which entails logarithmic divergent density of states and weak multi-fractal nature as we numerically verify [40,70]. Intuitively, the low-energy delocalized Majorana mode arises from percolating through an extensive number of resonating Majorana zero modes trapped in Z 2 gauge fluxes in a weak pairing topological superconductor [68,71].
To gain more insight into the topology of the fermion, we calculate the Chern number of the fermion eigenstates [72][73][74] of varying energy by using the real-space formula based on concept of non-commutative Brillouin zone: Here x, y are the real-space coordinate operators which label the first quantized orbitals and generate the translation of crystal momenta, ⟩ is the single-particle eigenstate of the first quantized Hamiltonian matrix with energy . P (E) is the spectral projector where the singleparticle mode with energy smaller than −E is occupied, mimicking the Fermi level in the complex fermion system with number conservation. In a fermionic system with only fermion parity conservation, half of the singleparticle eigenstates are redundant, so we consider only E ≥ 0. The change of C(E) reveals the Berry flux carried by the fermion mode at the corresponding energy. From Fig. 5b, the delocalized mode near E ≃ 2.5(1)J is clearly associated with a topological quantum critical point separating two distinct Chern plateaus, that is robust against perturbation and weak disorder [75,76]. As for the energy window 1.5J ≲ E ≲ 2.0J, it is unclear whether it would maintain a finite mobility edge or shrink to a singular point or become fully localized in the thermodynamic limit [77,78].
A final comment is that our result is consistent with the argument that non-Abelian topological phases cannot be fully localized [79]. While disorder tends towards localization, i.e., a divergent dynamical exponent z → ∞, this tendency competes with the metallic and topology induced critical modes favoring ballistic propagation with z = 1, leading to the observed subdiffusive dynamics with z > 1.

VI. BEYOND EXACT SOLVABLE LIMIT
Now we aim to address the robustness of our observations upon the influence of interactions present in Eq. (3).
Here, we will focus on the leading order resonant contributions responsible for an eventual destabilisation, by utilizing the approach introduced in Ref. [80], where it has been shown that these resonant contributions can capture the essential non-perturbative effects of interactions such as the logarithmic entanglement growth in MBL phases not only on a qualitative but also on a quantitative level. The Hamiltonian is then expressed in the canonical fermion basiŝ where the canonical Majorana fermions γ ′ , γ ′′ are related to the original local Majorana fermions by an orthogonal transformation obtained from diagonalizing the non-interacting fermion part. The leading order resonant interaction preserves the parity of the canonical fermion mode ⟨iγ ′ n γ ′′ n ⟩ but induces a dephasing effect, analogous to the l-bit theory in MBL systems [12,13] in which it leads to dramatic non-perturabative effect [80]. In our case of Kitaev model with weak magnetic field, we also find a special structure for this interaction strength V m,n , which endows a hierarchy of dephasing timescales (see Fig. 6a). It is the key observation in Ref. [80] that the dynamics of any fermion correlation function can be effectively written as a sum over a number O(N 4 ) of Gaussian evolution trajectories, schematically abbreviated as which can be further factorized into the product of a Loschmidt amplitude quantity and an effective correlation function [40, [81][82][83][84][85]. As shown in Fig. 6d, we calculate the system with 16×4 unit cells (128 spins) up to the timescale Jt ≤ 10 4 , where the interacting scenario turns out to collapse with the non-interacting case within numerical accuracy. The critical quantum correlations are therefore stable up to this timescale. We estimate the validity of the perturbative approach by statistics of resonances in first-order perturbation theory of the omitted terms. We find that they are off-resonant with probability ≳ 99.5% for the considered parameter regime [40], above typical thresholds [86], so that they can become relevant only via higher-order processes manifesting on longer timescales.

VII. CONCLUDING DISCUSSION
While thermalization may occur eventually for the full Hamiltonian in Eq. (1) on long time scales, we emphasize that our findings imply a long intermediate time window with nonergodic behavior leading to exotic quantum dynamics and correlations. On the other hand, it is tempting to ask whether the exotic dynamics and non-equilibrium quantum order are related to the lowenergy non-Abelian Ising topological order [20,87]. However, in the strongly anisotropic coupling regime which in zero temperature exhibit a distinct Abelian Z 2 topological order [88], we find similar high energy critical mode and subdiffusive dynamics as well as critical correlation, which goes beyond the low energy universality class [40]. Our findings of subdiffusive dynamics and critical quantum correlations may emerge universally in general Z 2 lattice gauge theories coupled to chiral Majorana matter fermions, provided two essential ingredients: (i) nontrivial fermion topology; (ii) static disordered gauge flux [89][90][91][92]. Above all, our observation of quasi-long range order with associated divergent multipartite entanglement in a non-equilibrium high-energy steady state marks a concrete first step towards yet unexplored unconventional phase structures in ergodicity breaking two dimensional quantum models. Subdiffusion might be present also in a more general context of Majorana spin liquids as long as the effective Z 2 gauge flux (vison) dynamics is much slower than that of the fermions, which will be a challenging but valuable scope for future research. Furthermore, motivated by the recent developments showing that the gauge charge disorder in a 1D unconstrained gauge theories can stabilize a time crystal order [93], it would also be interesting to generalize this idea to two dimensions in search of nontrivial spatiotemporal order from a driven Kitaev model [94,95]. Finally, beyond the conceptual interest, the nonequilibrium quench dynamics can in principle be realized in various quantum architectures including ultracold atoms [96,97], superconducting qubits [98,99] or topological nanowires [100,101], as well as via ultrafast pump-probe techniques in the Kitaev candidate materials at sufficiently low temperatures with suppressed phonon influence [29, [102][103][104][105][106]. The quantum quench protocol with a prequench spin product state excites all the allowed gauge configurations like in a thermal ensemble of infinite temperature, where the typical gauge configuration is maximally random. It is also interesting to see the dynamics interpolating between the clean limit and this dirty limit, analogous to tuning the temperature from zero to infinite in a thermal ensemble. In general there could be three related but inequivalent ways of interpolating between the zero gauge flux configuration and the typical maximally random gauge flux configuration. The first way is to endow any specific gauge field configuration with a thermal weight depending on the corresponding Majorana fermion free energy, which is commonly used in Monte Carlo calculations for finite temperature thermal ensemble [23]. The second way is simply to tune the density of π link: p(u = −1) ∈ [0, 0.5] [67]. But notice that the gauge field on a link is not gauge invariant quantity. The third interpolation way is to tune the average gauge-invariant flux [24] by applying the spin ring exchange interaction operator ∏ q (1 2 + (1 2 − p)Ŵ ) to the initial state, which is equivalent to tuning the density of π flux on a hexagon plaquette from 0 to 0.5. One could even further extend this range to p ∈ [0, 1] to interpolate between the absolute 0 flux and fully packed π flux gauge configurations [68], with the dirty limit lying in between. This manner of controlling the density of π flux can be realized by deforming the quench protocol and hence the initial state Ψ 0 (p)⟩, where the fermion maintains in a gauged vacuum in each gauge configuration ⟨Ψ 0 (p) iuαβ Ψ 0 (p)⟩ = 1. As shown in Fig. 7a, at intermediate time, the projected Majorana bipartite entanglement entropy grows algebraically with non-monotonically varying exponent depending on the density of random π fluxes. In Fig. 7bc, the late time steady state exhibits algebraic dimer correlation functions in space, with a non-monotonically continuously varying exponent. Overall, the cleaner system with p close to 0 or 1 have faster information propagation, but less prominent critical correlation in steady state. Notice that the maximally random case p = 0.5 does not yield the slowest entanglement propagation nor the largest correlation exponent, which might be due to the asymmetry between the clean 0 flux sector and π flux sector. From this discrete sequence of disorder density, we do not find divergent tendency for the dynamical exponent like in the MBL phase transition, which is consistent with the fact that non-Abelian topological phases of matter cannot be fully localized [79]. In Fig. 7d, we show the spatiotemporal profiles of the dimer quantum correlation functions, where a continuously varying light-cone is witnessed. The clean limit without disorder exhibits oscillation inside the lightcone due to the finite size non-interacting nature, which is damped by the onset of disorder.

Appendix B: Other Majorana fermion correlations
Throughout the main text we mainly consider the physical observable of spin dimers, which is equivalent to the next nearest neighbouring gauge invariant Majorana fermion bilinear term. This is the simplest and most natural choice in both spin and fermion representation. However, the spin dimer correlation might be much more difficult to measure in experiments than the spin correlation. The single bare spin operator is composed of single matter Majorana fermion and an auxiliary Majorana operator that flips the gauge connection and subsequently the fluxes. Nevertheless, it was shown in Ref. [22] that the flux-conservation-breaking perturbation could generally dress the spin operator and lead to a contribution of pure matter Majorana fermion bilinear term that preserves the flux. For example, to leading order σ z ≈ σ z +f σ x σ z σ y +⋯, where the second term is just the three-spin interaction in an neighbouring wedge, transformed to the next nearest neighbour hopping term of the Majorana fermions on the same sublattice. In low energy, this term contributes to the chiral mass of the Majorana fermion. This pure matter contribution qualitatively changes the low temperature low frequency spin spin correlation function [22]. Here we also perform a modest calculation for the spreading of (i) correlation between two Majorana fermions on the same sublattice, separated at a distance; (ii) connected correlation between two next-nearest-neighbour-Majorana-fermion-bilinears on the same sublattice. The latter one should contribute to the spin spin correlation that is accessible by experimental probe. As shown in Fig. 8, the spread of this correlation appears to be consistent with the subdiffusive fermion entanglement growth we show in the main-text. In the late time steady state, the correlations also appears to show a power-law signature. Therefore, the subdiffusive dynamics and critical correlation we obtain may leave signatures in the pump-probe experiments for transient spin dynamics. (f) It is tempting to ask whether the critical dynamical phase we obtain out of a nonequilibrium quantum quench is associated with the low energy topological phases. To answer this question we may tune the Kitaev interaction to be strongly anisotropic, by chooisng J z = 2J, J x = J y = 0.5J,h = 0.25J. In this case at zero temperature it belongs to the Abelian Z 2 topological phase, in the same universality class with the celebrated toric code model. In fact, the low energy physics of the anisotropic Kitaev honeycomb model can be mapped to the toric code model [20]. However, at high energy density the connection between the Kitaev honeycomb model and the toric code model is not a priori known.
The perturbed toric code model was argued to be many-body localized in the presence of strong disorder [41], since it can be dual to the Ising model when charge is absent [88]. The anisotropic Kitaev honeycomb model in the presence of magnetic field, on the other hand, is found to exhibit subdiffusive spreading of quantum correlation and  entanglement, see Fig. 9 and Fig. 10. Notice that the spreading of dimer correlation is anisotropic, and that along z direction has the same exponent as the spreading of entanglement. Behind that we also find the topological critical fermion modes at high energy for the random flux configurations or π-flux configuration, even though the zero energy ground state is topological trivial, see Fig. 11. It means that unlike the 0 flux sector, the random flux sectors cannot be adiabatically connected to the toric code limit without closing the mobility gap at certain energy. The late-time steady state also exhibits critical dimer correlation along z-direction, see Fig. 12. These evidences show that the critical dynamics phase goes beyond the low energy effective theory. The topological critical Majorana modes at high energy are responsible for the anomalous subdiffusive dynamics and critical correlation. For completeness, we also compute the dynamics for the time reversal symmetric Kitaev model in the absence of magnetic field h = 0, as shown in Fig. 13. In such case, the problem is equivalent to a free Majorana fermion with only nearest neighbour hopping on the honeycomb lattice, that is exposed to a random π flux penetrating the plaquette. Notice that without magnetic field there is not only time reversal symmetry but also the sublattice (chiral) symmetry (sign change to one sublattice changes the sign of the Hamiltonian), which connects the positive energy to negative energy but acts as a unitary symmetry distinct from the anti-unitary particle hole symmetry. The zero energy single particle mode has exact chiral symmetry, and therefore often behaves qualitatively distinct from the other energy  states in the presence of disorder. A closely related problem is the gapless Dirac fermion in the presence of random U (1) magnetic flux disorder, see Ref. [89][90][91]. Here we are dealing with a Majorana version of the Dirac fermion in the presence of π flux disorder, with additional particle hole symmetry. As shown in Fig. 14, our numerical result suggests a diverging localization length towards zero energy. Besides, there seems to be multiple singular points.
For comparison, we summarize the dynamical exponents fit from different physical observables and for different parameter regime as in Table. I.  Below we review two alternative fermionization approaches for the Kitaev model, and draw the connections between them. Generally speaking, Kitaev's original gauge theory approach is physically transparent and looks more symmetric in putting gauge field on every link, while the Jordan Wigner transformation approach [35] is practically convenient without gauge redundancy. The latter can be achieved from the former by gauge fixing, with proper care taken for the boundary terms.

Jordan-Wigner transformation
As a warmup let us first consider only one zig-zag row composed of x-links and y-links: ∑ j∈A σ x j σ x j+1 + ∑ j∈B σ y j σ y j+1 . It can be readily solved by performing a zigzag Jordan-Wigner transformation: such that the Ising exchange interactions on x(y)-links are simplified to Majorana fermion hopping terms Now consider the honeycomb lattice, where spins are ordered row by row to be Jordan-Wigner transformed, and we label the sites by row and column indices (i, j). The above intra-row interactions remain unchanged, while the inter-row couplings are simply where a static link variable u z i,j = ±1 emerges on every z-link. Therefore the three anisotropic Ising interactions become free Majorana hopping terms coupled with classical Z 2 variables. Since Jordan-Wigner transform is a nonlocal unitary, the boundary term in each row yields a non-local term: where u y i commutes with the Hamiltonian therefore serves as a conserved quantum number, separating the fermion system into disconnected sectors with distinct boundary condition. In other words, the spin model with a definite boundary condition contain two fermion sectors with distinct boundary condition for fermions, which is reminiscent of the 1D quantum Ising critical point.
Here we briefly comment that the boundary condition in fermion basis is relevant to the topological ground state degeneracy, which is generated by moving two Z 2 local gauge fluxes around the torus before annihilation, equivalent to pumping a global Z 2 flux for fermions [87]. In the Abelian phase (toric code), one obtains 4-fold degeneracy by the freedom of pumping global Z 2 flux along either direction. In the non-Abelian Ising phase, there are only 3-fold degeneracy, because due to the Majorana zero mode trapped by the flux, pumping global flux along both directions would leave a global fermionic excitation behind and take the state away from ground state manifold.
Under such Jordan Wigner transformation, the prequench initial spin polarized product state that satisfies is transformed to a product state of local Majorana dimers satisfying ib j∈A α j = 1, iβ j b j∈B = 1. By grouping (b, b) and (α, β), it is a maximally entangled EPR state between the local Majorana bilinears and the gauge field on z-link: Therefore the prequench state is a gauged fermion vacuum state with disordered gauge configurations. The postquench Hamiltonian can be viewed as a Bogoliubov-de-Gennes superconducting Hamiltonian with gauged pairing and hopping terms, which would create fermion pairs from the gauged vacuum and allow them to propagate. On the other hand, such initial state projects onto the sector that locks opposite boundary condition between spin and fermion representations i.e. −σ y σ y = −iαβ on boundary while σ y σ y = −iαβ in the bulk. We'll work with the anti-periodic boundary condition for the spin model, and hence periodic boundary condition for the Majorana fermions.

Redundant gauge theory
Kitaev's original solution [20] is to extend the Hilbert space by rewriting each spin operator by four Majorana fermions with a physical constraint: Then the spin Hamiltonian rewritten in terms of Majorana fermions are always trivially conserved by e iπG ≡ b x j b y j b z j c j . In other words, the subsequent Hamiltonian in terms of Majorana fermions by definition has a Z 2 gauge symmetry, generated by e iπG . Notice that G is the Z 2 analogue of the Gauss law operator n j −∑ k E j,k as in the conventional U(1) gauge theory. The constraint b x j b y j b z j c j = 1 actually picks up the gauge neutral sector among all other super-selection sectors to be the physical space, as a superposition of all gauge equivalent configurations.
In the exact solvable regime, the Z 2 gauge field on every link u j−e x(y)(z) is conserved. To avoid extensive redundancy, one can always use the gauge transformation to fix the gauge such that u x(y) j = 1 on every link in the bulk. In other words, the gauge field is oriented only towards z-direction, analogous to the Landau gauge. However, one should pay special attention when the lattice is not on a planar but a compact manifold such as torus. In a torus, there could be global fluxes which can neither be detected locally nor be locally gauged away. Therefore, there must remain u x(y) j ≠ 1 on certain boundary links that account for the global fluxes. One can also see this by taking the product of gauge neutral constraints along one zig-zag row: ∏ j∈row u x j u y j = − ∏ j∈row σ z j . Therefore one may fix u x(y) j = 1 in the bulk but must allow one boundary x(y)-link on every row fluctuating, which is equal to − ∏ j∈row σ z j . In this way, we reach exactly the same fermionic Hamiltonian with a nonlocal boundary term, as derived from Jordan-Wigner transformation. And in this way we shall see that the Z 2 variable emergent in Jordan Wigner transformation has the physical meaning of a gauge field.
Last but not least, one should notice that the onset of magnetic field does not explicitly spoils the Z 2 gauge symmetry of fermion since that is just a gauge redundancy in the extended Hilbert space. Instead, the magnetic field breaks another set of local Z 2 symmetries which stand for the gauge flux conservation.
which is skew-symmetric matrix. The matrix elements immediately give the spin dimer expectations expressed in terms of gauged Majorana fermion bilinears: The disconnected dimer correlation functions expressed in terms of the four-Majorana-fermion correlation functions can be decomposed using Wick's theorem In each gauge configuration, the time evolved density operator is Gaussian. Therefore the reduced density operator of the half-partitioned fermion system is also a Gaussian density operator that can be disentangled into reduced canonical fermions denoted as d n , n = 1, ⋯, N 2: This Gaussian density operator is faithfully encoded in the covariant matrix which is just the subsystem block of the total covariant matrix Γ(t). Therefore one can extract the eigenvalue of the disentangled density operator (occupation of the reduced canonical fermion mode) and the entanglement entropy: The upper-bound is reached if and only if each canonical fermion is maximally entangled with a flat spectrum for the subblock of covariant matrix i.e. ρ n = 1 2.
that satisfy . For convenience here we use i, j, k, l to label the spin site instead of unit-cell. Consider two inversion related Y junctions tied to the same z-link, with site labeled in counter-clockwise ordering, where cycl.perm denotes the cyclic permutation counterparts. One could qualitatively verify the above result by the fact that the chiral three spin interaction σ x A(B) σ y A(B) σ z A(B) respects three-fold rotation symmetry C 3 but behaves odd under either time reversal or mirror reflection symmetry And M z maps β 1 α 1 α 2 α 3 to α 1 β 1 β 2 β 3 , which enforces the sign difference. Inside the weight, taking the imaginary part of the wave-function overlap between the same sublattice enforces the mirror reflection odd condition. These diagonal part contributes to the leading order resonant interaction. Generally, the omitted off-diagonal terms can be Schrieffer-Wolff rotated to yield higher order correction to the resonant interaction, which usually only plays a role in exponential longer time. Here we just take the leading order term with the gauge coupling to get the symmetric interaction coupling matrix: i,m ψ j,m Im ψ * k,n ψ l,n + cycl.perm + (m ↔ n), in which j, k, l are nearest neighbours arranged in counter-clockwise order surrounding site i, and i is summed over both A and B sublattices. In this way, we arrive at our final effective Hamiltonian to leading order: For completeness we here show the conserved canonical fermion parity in each random gauge configuration ⟨iγ ′ n γ ′′ n ⟩, and the corresponding mean-field shift of energy due to the resonant interaction, see Fig. 15.

Statistics of omitted off-diagonal interaction
Here we do a statistical analysis for the off-diagonal interaction terms that were omitted. For simplicity let's fix the gauge such that nonzero gauge connection occurs only at z-link, and take a unit-cell and label the 6 sites involved as in inset of Fig. 16. The Majorana Hubbard interaction is the product of certain Majorana fermion bilinears on bondŝ V =h ∑ unit−cell u 0 (iα 0 β 1 iβ 2 β 3 + iα 0 β 1 iα 4 α 5 ), which can be straightforwardly verified by σ x 2 σ z 1 σ y 3 = (σ z 0 σ z 1 ) (σ x 2 σ z 0 σ y 3 ) = (−u 0 iα 0 β 1 ) (iβ 3 β 2 ), likewise for σ x 4 σ z 0 σ y 5 . Using α j = √ 2 ∑ n (−iψ j,n c n + h.c.) , β j = √ 2 ∑ n (ψ j,n c j,n + h.c.), the involved fermion bilinears can be decomposed into off-diagonal and diagonal parts: A rough overview: the product of two diagonal fermion bilinears γ ′ p γ ′′ p γ ′ q γ ′′ q contributes to the l-bit type resonant interaction terms that have been shown previously. The product of one diagonal bilinear γ ′ p γ ′′ p and one off-diagonal bilinear γ m γ n contributes to the off-diagonal interaction term that flips the fermion parity of two canonical fermions, which can be further decomposed into an assistant pairing term ∼ 1 − 2c † p c p c † m c † n and an assistant hopping term ∼ 1 − 2c † p c p c † m c n , mediated by the presence or absence of other canonical fermion modes. Since the density term depending on index p is to be integrated out, the assistant pairing/hopping interaction is roughly proportional to the overlap of canonical fermion wave-function of m and n on the same given unit-cell. In contrast, the product of two off-diagonal bilinears γ m γ n γ p γ q (m ≠ n ≠ p ≠ q) that changes the fermion parity of four canonical fermion modes requires the overlap of four localized fermion wave-function on the same unit-cell, which are less dominant. In the following we do statistical analysis for the assistant hopping and pairing terms.
In detail, for the Kitaev Hamiltonian at solvable point, we sample 200 random gauge configurations {u j } where j labels the sites. For each gauge configuration we randomly sample 100 eigenstates, which are determined by the canonical fermion pairity array {ν n }, with ν n ≡ iγ ′ n γ ′′ n = ±1being randomly sampled. We consider the assistant hopping(pairing) interactions connecting these random states to the target states with canonical fermion modes m, n FIG. 16. Cumulative probability distribution function(CDF) for the off-diagonal interaction matrix element divided by the corresponding energy difference, in the eigenstate basis of exact solvable Hamiltonian in random gauge sectors. Varying system sizes are shown, as well as the assistant hopping and pairing types. The assistant hopping term is relatively easier to reach resonance than pairing term, with probability ≲ 0.5%. In comparison, in the zero flux sector without effective gauge flux disorder, the assistant hopping interaction shows strongly resonant contribution with ln V ∆E > 10 6 with probability ≳ 4% among the same randomly sampled eigenstates. Therefore it is the disorder that suppresses the strong resonance. Parameters: Jx = Jy = Jz = J,h = 0.25J, sampled over 200 disorder realization of random gauge configurations, with 100 randomly sampled eigenstate in each configuration.

Correlations with dephasing interaction
In the presence of dephasing interaction, we need to transform the physical observables into the canonical fermion basis: The two-dimer disconnected correlation function is equivalent to the four-Majorana correlation functions as follows: 24A[f zz m,n,p,q ]⟨γ m (t)γ n (t)γ p (t)γ q (t)⟩.

(G13)
Notice that due to the orthogonality of Q matrix ∑ m Q j,m Q i,m = δ i,j , the terms with overlapping indices drops out such that only off-diagonal Majorana correlation functions with m ≠ n ≠ p ≠ q contribute. A[f zz m,n,p,q ] means antisymmetrizing the form factor f zz with respect to indices m, n, p, q, using the anticommutation relation of Majorana fermions. Likewise for ⟨D x ⟩ and C xx .
Unlike the free fermion case, the four-point Majorana correlation function cannot be immediately factorized by Wick's theorem into simple product of two-point Majorana correlation functions. Each canonical Majorana fermion doublet γ n ≡ (γ ′ n , γ ′′ n ) T is evolved effectively by a particle-hole superposition of Gaussian operators: Then the two-point Majorana fermion correlation functions can be expressed as The key is to evaluate the expectation of e 1 4 γ(µAm+νAn)γ γ m γ T n over the initial state ψ {u} ⟩. This can be generally done because the initial density matrix and the exponential operator are both Gaussian. For simplification in a given gauge configuration, we have the A sublattice absorb the neighboring gauge field on z-link α j → u j α j such that initial state is rotated to a clean Fock vacuum, while the gauge field dependence is absorbed by Q → (u, 1)Q matrix. Turning back to the original fermion basis γ T Aγ = ζ T QAQ T ζ where µA m + νA n ≡ A and ζ j labels the Majorana doublet in unit-cell j, we have where iα j β j vac⟩ = vac⟩, and every eigenvector matrix Q has been regauged accordingly. The quantity in the middle is a vacuum expectation of the product of Gaussian operator and Majorana fermions, which can be evaluated using the generic formula that is to be derived later. Likewise, the four-Majorana-fermion disconnected correlation function can be reduced to the expectation of the product of Gaussian operator and four Majorana fermions: (G17) By denoting µA m + νA n + κA p + λA q ≡ A, we can factorize the term similarly and evaluate the effective correlation function part by Wick's theorem: Likewise, the Loschmidt amplitude Z and the effective two-point correlation functionΓ can be evaluated using the general result we are going to derive in the next section. In numerical computation, although each four-Majorana-fermion correlation function with fixed m, n, p, q reduces to effective Gaussian evolution and can be efficiently calculated, there are ∼ O(N 4 ) number of independent Gaussian evolving trajectories in total, which is a huge computation complexity to keep track of. Fortunately, it is easy for parallelization. Our computation amounts to N = 16 × 4 unit-cells, which accounts for 128 spins, being averaged over 200 randomly generated disorder samples. For completeness, here we show the slices of correlation growth in fixed distances, and compare the interacting case with non-interacting case, for both isotropic coupling and anisotropic coupling respectively, see Fig. 17.

Extended Majorana correlation function
In this section we are going to derive an independent formula for a generic complex extended Majorana correlation function as below: (G20) where the vac⟩ is defined by ⟨vac iγγ T vac⟩ = i(1 + τ y ), γ stands for a column vector of Majorana fermion operators. The complex correlation function is factorized into a product of a Loschmidt amplitude and an effective correlation function. In general, the input argument A can be any generic complex antisymmetric matrix that has the canonical form, and K is the matrix that would appear as the Gaussian exponent coupled to the pairing term in the Balian-Brezin decomposition: This formula can be applied to evaluate multiple expressions in the former section, where we turned back to the original fermion basis ζ, and regauge the initial state to be a fermion vacuum of ζ. The Pfaffian can be numerically evaluated efficiently using the algorithms and codes provided by Ref. [83].
Detailed derivation for this formula are given as follows. To evaluate the Loschmidt amplitude without sign ambiguity, we first derive a generic formula for the expectation of generic Majorana fermion Gaussian operator over the fermion vacuum, by virtue of the Baker-Campbell-Hausdorff formula for tracing out sequence of Majorana Gaussian operators as relevant in the generalized Levito's formula Tr(e n cosh a n 2 . (G22) Notice that when getting rid of the square root there could be a sign factor, which is fixed to be positive in this case by analytic continuity such that Z(A) → 1 when A → 0, as long as Z(A) remains analytic along the path of A variation and does not trespass a zone with condensed Lee-Yang-Fisher zeroes in thermodynamic limit. In finite size system, Z(A) is a finite order polynomial complex function where isolated Fisher zeroes can always be avoided by infinitesimal variation. Notice that the final expression is rather complex compared to a simple √ det at the beginning lines of derivation, but being rewritten in terms of Pfaffian has the advantage of being free from sign ambiguity [81]. Otherwise, summation over different trajectories with a sign ambiguity would lead to a sign problem like in quantum Monte Carlo calculations typically for fermion systems.
Next we derive the generic correlation matrix defined with respect to an effective Gaussian operator vac⟩⟨vac e where we denote T ≡ e A BdG and the subscript labels particle-hole blocks. Notice that the above can also be directly derived by using the Balian-Brezin decomposition for ⟨vac e 1 4 γAγ c † c † vac⟩ = ⟨vac e cKc c † c † vac⟩, where the anti-symmetric matrix Turning back to the Majorana fermion basis, we have that completes the derivation. Finally we also comment that the expression above can also be consistently obtained by using the general formula of tracing out the product of two generic Gaussian operators, and the formula of the product of two Gaussian correlation functions [82] up to some further simplification. Namely, for any generic square skew-symmetric matrices A and B, Tr(e One way is to generalize the Brillouin zone using the concept of non-commutative geometry, which manifest in replacing summation by matrix trace in calculation, and is guaranteed to be convergent and topological if in the presence of a mobility gap [20,61,63]. The formula we need to calculate is where P is the spectral projector to be concretely defined below, the real space continuum coordinate operator x(y) is to be expanded by the lattice coordinate in a properly defined series [62] such that the finite size error is exponentially small. Our numerical result in the main text is based on this method. The other way is to twist the boundary phase instead [73], which takes the spirit of Laughlin's gedanken experiment of flux pumping, and can be understood as generating a super-lattice with effective quasi-momentum. In numerical calculation for a finite size lattice of discretized twisting boundary phase, it has great advantage to regularize the problem into the form of a lattice gauge theory, whose global flux over the twisting phase space is guaranteed to be quantized [74]: where φ is the twisting boundary phase, ψ φ ⟩ is the corresponding single particle or manybody wave-function. The other advantage of this approach is that for a free fermion problem it results in eigen-state resolved Chern numbers, from which one can determine the density of extended states [75,77,78], as shown in Fig. 19. In the following we briefly review the technical details for these two methods. To avoid confusion, all the notations work only within the following sub-section independently and does not apply to the rest of the supplemental material.

Calculation using real space formula
Even without translation symmetry, one can still perform a Fourier transform and define momentum, although it is no longer good quantum number, because the single particle Hamiltonian in momentum space carries off-diagonal elements connecting different momenta. Naively, one can simply generalize the summation over momentum to matrix trace over of the spectral projector. Concretely, first let's denote the single particle spectral projector to the ground state in real space basis: P = 1−sgn(H) 2 , where sgn(H) is obtained by flattening the spectrum of the single particle Hamiltonian H but the sign. The spectral projector in momentum space and its derivative arẽ One may generalize the Chern number formula (I4) In the end, it is simply a commutator of the real space coordinate operator being projected into the Fermi sea. One could immediately verify that when P = 1 ⇒ C ∝ [x, y] = 0, the contribution over the total Hilbert space is guaranteed to vanish. The essential questions are whether it converges to integer and whether it is of topological nature. First, it was elaborated early by Bellisard using the concept of non-commutative geometry [61], and later Kitaev also gave an intuitive argument for the topological nature of the 2-current of the spectral projector, and its quantization as the flow of a quasi-diagonal unitary matrix reminiscent of the Laughlin's flux pumping gedanken experiment F e i2πP xP = Tr e i2πP xP P yP e −i2πP xP − P yP N = ∫ 2π 0 Tr ∂ φ e iφ[P xP,] P yP dφ N = 2πi[P xP, P yP ] N [20]. Note that the quasi-diagonal condition of the unitary matrix is readily satisfied when there is a spectral or mobility gap so that correlation exponentially decays. Note also that Bianco et.al. also derived a similar formula from the linear response theory for continuum real space system [63], and they even took out the local contribution in each unit-cell i.e. the diagonal entry of the commutator 2πi[P xP, P yP ] and defined a so-called local Chern marker, which in a clean system is equivalent to the total Chern number. However, one should notice that the trade of k-derivative with commutator of coordinate ∂ k → [ir, ] is exact only at thermodynamic limit N → ∞ such that ∆k → 0 can be infinitesimally small. Prodan et.al. proposed an efficient numerical algorithm to exponentially suppress the finite size error in this line as follows [62]. In a finite size system with L units along x-direction, one could approximate the real-space coordinate by an expansion x ≃ ∑ L 2 n=1 a n e in 2π L x − c.c. . As what we need is not to cover the whole lattice coordinate, but to approach the continuum limit of small x, a direct discrete Fourier transform for the lattice coordinate is apparently not the optimal choice, which yields polynomial O(1 L 2 ) error. One optimal solution is chosen by solving the linear matrix equation ∑ L 2 n=1 a n n 2j−1 = L 2π δ j,1 , (j = 1, ⋯, L 2) such that the deviation between the continuous x and the expanded series is exactly 0 up to L 2-th order, leading to exponentially small error with the system size: O 1 L L . As coordinate operator generates the translation of momentum, we have

Calculation by twisting boundary phase
It is akin to the momentum space formula for the Chern number [72,73,77]. By discretizing the twisted boundary phase space into M x × M y sites, it is equivalent to treat the L x × L y finite-size system as a giant unit-cell and generate a periodic super-lattice of size L x × L y × M x × M y , thus effectively generating a mini-Brillouin-zone(BZ) for the quasimomentum in unit of 2π M x(y) . In the mini-BZ, the L x × L y large unit cell folds into L x × L y mini-bands. Since the twisting phase φ , conjugate to the real space coordinate, plays the role of a quasi-momentum ranging between [0, 2π L] [73]: where ψ(φ) can either be the single particle eigenstate or manybody eigenstate for the Hamiltonian with twising boundary phase φ. Physically, the twisting boundary phase also concurs with Laughlin's flux pumping gedanken experiment. In a finite-size lattice system, one could regularize the above formula like regularizing the continuum gauge theory into a lattice gauge theory [74]. Namely, regularize the local flux integrated over a plaquette as the Wilson loop of the U (1) link variable, which is physically the parallel transportation of the wave-function: where arg is defined as taking the phase angle from the principal branch (−π, π]. Notice that the parallel transporter circulating the BZ should be equivalent to identity due to the periodicity (close manifold without boundary), which means the total flux must be quantized in units of 2π, therefore guaranteeing the quantization of Chern number defined on this lattice. Whether this global flux coincides with the continuum limit is a question. While there is no limitation to the flux on each plaquette in the continuum limit as a noncompact gauge theory, there is 2π ambiguity of flux on each plaquette in the lattice gauge theory with compact gauge group. This is the source of the discrepancy between the lattice and the continuum limit. Generally, when the phase grid of the lattice is fine enough such that the corresponding flux on each plaquette in continuum limit is controlled within the principal branch window, the finite size result should unambiguously concurs with the continuum limit. Assuming a smooth Berry curvature configuration, by uniformly distributing the Berry flux onto each plaquette approximately, one can estimate a lower-bound for the required size of the phase space: M x M y > 2 C . As we know, for a low-energy Dirac fermion with light mass, the Berry curvature is sharply peaked around the Dirac point, but the distribution is flattened when the mass is large. The quantization of Chern number in continuum limit relies on the adiabatic evolution of the ground state wavefunction circulating around the twisting phase space. In other words, it also requires a mobility gap at the Fermi-level, or a (exponentially) fast decay of the correlation function, which excludes the situation with Fermi-level crossing a band in a clean system. In cases with band overlapping or level crossing, the Berry curvature could have singularity. However, in a generally disordered finite size system, this generally holds because either the states near the Fermilevel is localized without contribution to Berry curvature, or the delocalized states near Fermi-level are generally experiencing level repulsion due to scattering. Notice that the Bloch state in clean system without momentum scattering and manybody interaction is a very special case in this sense.
In practice, one can either use the manybody ground state (such as the slater determinant for free fermion system) to construct the link variable, or resolve the Berry phase for each single particle eigen-state, which gives biproduct of determining the localization nature of a wave-function [75,77]. Further, one can perform scaling analysis for the density of extended states versus the total states to deduce the existence of extended states in thermodynamic limit that escapes localization [78].

Appendix J: Zero energy thermal metallic state
In this section we show numerical verification for the delocalized nature of the single particle mode at zero energy for the isotropic coupling model for larger system size, which is qualitatively consistent with the infinite temperature limit of a moderate size numerical calculations in Ref. [69]. To resolve the density of states at extremely small energy, we use the standard kernel polynomial method [70], to expand the density of states, defined as a collection of δ functions, in the orthogonal basis of Chebyshev polynomials, Tr (T n (H E max )) T n (E E max ) .
where T n is the first kind of Chebyshev polynomials, and the coefficient at each order can be evaluated iteratively in terms of sparse matrix. 10 random complex states are initiated to evaluate the trace of matrix stochastically. In this way we can go to very large system size with finer energy resolution. When truncated at finite orders, a Jackson kernel is attached to damp the Gibbs oscillations, resulting in a smooth regular function. Results are shown in Fig. 20ab. A logarithmic divergent behaviour close to zero energy is visible, consistent with the theoretical prediction [65].
To verify the weak multi-fractal nature of the single fermion wave-function at zero energy, we first look at the fractal dimension of the generalized inverse participation ratio [76]: where p(r) ≡ ψ(r) 2 is the zero energy single particle wave-function probability distribution, L is the coarse grained linear system size. In practical numerical calculation, to reduce the interference of the lattice length scale, we typically do a coarse graining for p(r), where r is the coarse grained coordinate. To calculate the singularity spectrum f (α) that is related to D q by the Legendre transform D q (1 − q) = f (α) − αq, f ′ (α) = q, we first define a set of normalized weight µ q (r) = p(r) q ∑ r p(r) q . Then the scaling dimension of q-th moment weighted average probability, and the associated volume scaling dimension are respectively given by α(q) = ∑ r µ q (r) ln p(r) − ln L , f (α(q)) = ∑ r µ q (r) ln µ q (r) − ln L . (J3) It can be checked that f (α(q = 0)) = 2. Similar to a 2+ dimensional metal near the Anderson localization, our results show weak multi-fractal bahaviour D q = 2 − γq, when γq ≪ 2, and α(q = 0) = d − γ, as shown in Fig. 20cd.