An exact chiral amorphous spin liquid

Topological insulator phases of non-interacting particles have been generalized from periodic crystals to amorphous lattices, which raises the question whether topologically ordered quantum many-body phases may similarly exist in amorphous systems? Here we construct a soluble chiral amorphous quantum spin liquid by extending the Kitaev honeycomb model to random lattices with fixed coordination number three. The model retains its exact solubility but the presence of plaquettes with an odd number of sides leads to a spontaneous breaking of time reversal symmetry. We unearth a rich phase diagram displaying Abelian as well as a non-Abelian quantum spin liquid phases with a remarkably simple ground state flux pattern. Furthermore, we show that the system undergoes a finite-temperature phase transition to a conducting thermal metal state and discuss possible experimental realisations.


I. INTRODUCTION
Amorphous materials are condensed matter systems characterised by short-range regularities, and an absence of long-range crystalline order as studied early on for amorphous semiconductors [1,2].The bonds of a wide range of covalent compounds can enforce local constraints around each ion, e.g. a fixed coordination number z, which has enabled the prediction of energy gaps even in lattices without translational symmetry [3,4], the most famous example being amorphous Ge and Si with z = 4 [5,6].Recently, following the discovery of topological insulators (TIs), it has been shown that similar phases can exist in amorphous systems characterized by protected edge states and topological bulk invariants [7][8][9][10][11][12][13].However, research on electronic systems has been mostly focused on non-interacting systems with few exceptions, for example, to account for the observation of superconductivity [14][15][16][17][18] in amorphous materials or very recently to understand the effect of strong electron repulsion in TIs [19].
Magnetism in amorphous systems has been investigated since the 1960s, mostly through the adaptation of theoretical tools developed for disordered systems [20][21][22][23] and with numerical methods [24,25].Research has focused on classical Heisenberg and Ising models, which are able to describe ferromagnetic, disordered antiferromagnetic and widely observed spin glass behaviour [26].However, the role of spin-anisotropic interactions and quantum effects in amorphous magnets has not been addressed.It is an open question whether frustrated magnetic interactions on amorphous lattices can give rise to genuine quantum phases, i.e. to long-range entangled quantum spin liquids (QSL) [27][28][29][30].
The combination of a fixed local coordination number in conjunction with magnetic frustration generated by bond-anisotropic Ising exchanges can lead to stable QSL phases.The seminal Kitaev model on the honeycomb lattice [31] provides an exactly solvable model whose ground state is a QSL characterized by a static Z 2 gauge field and Majorana fermion excitations.Several instances of Kitaev candidate materials have been synthesized in the last decade [32][33][34][35][36] following the suggestion that heavyion Mott insulators formed by edge-sharing octahedra may realize dominant Kitaev interactions [32].In particular, recently it has been shown that the Kitaev material Li 2 IrO 3 can be created with an amorphous structure [37].In fact, with sufficiently fast cooling, any crystalline material can be made amorphous [2,38], opening the possibility for exploring a wide variety of non-crystalline Kitaev materials.
It is by now well known that the Kitaev model on any three-coordinated (z = 3) graph has conserved plaquette operators and local symmetries [39,40] which allow for a mapping onto an effective free Majorana fermion problem in a background of static Z 2 fluxes [41][42][43][44].However, in general this neither means that any z = 3 lattice Kitaev model can be straightforwardly constructed, nor that the QSL properties are obvious.Several obstacles remain.First, the labelling of bonds necessary to create a soluble Hamiltonian can be an NP-complete problem.Second, once the Majorana system has been constructed, determining the ground state out of the exponentially large number of Z 2 flux sectors is generically hard, since Lieb's theorem -which defines the ground state flux configuration for the honeycomb -is not applicable for most lattices.Previous studies have relied on translation and reflection symmetries to reduce the number of sectors that must be checked [43,45,46], which cannot be done in an amorphous system.Third, once the ground state flux sector is found, it needs to be determined whether lattice disorder induces a gapless phase [47][48][49] or whether the fermionic spectrum is gapped, possibly with non-trivial topology [43].
In this article we study the Kitaev model on amorphous lattices and establish it as an example of a topologically ordered amorphous QSL phase.Concentrating on random networks generated via Voronoi tessellation [7,9] with z = 3, we show how to colour the bonds consistently.We find that the presence of plaquettes with an oddnumber of sites lead to a chiral QSL with spontaneously broken time-reversal symmetry (TRS) [43,[50][51][52][53][54][55].We establish via extensive numerics that the ground state Z 2 flux sector follows a remarkably simple counting rule consistent with Lieb's theorem [56].We map out the phase diagram of the model and show that the chiral phase around the symmetric point is gapped and characterized by a quantized local Chern number ν [7,57] as well as protected chiral Majorana edge modes.Finally, we discuss the effect of additional bond disorder and comment on the role of finite temperature fluctuations, showing that the proliferation of flux excitations leads to an Anderson transition to a thermal metal phase [47][48][49].

II. METHODS
We start with a brief review of the Kitaev model on the honeycomb lattice [31] before generalising it to amorphous systems.A spin-1/2 is placed on every vertex and each bond is labelled by an index α ∈ {x, y, z}.The bonds are arranged such that each vertex connects to exactly one bond of each type.The Hamiltonian is given by where σ α j is a Pauli matrix acting on site j, ⟨j, k⟩ α is a pair of nearest-neighbour indices connected by an α-bond with exchange coupling J α .For each plaquette of the lattice, we define a flux operator W p = σ α j σ α k , where the product runs clockwise over the bonds around the plaquette.These commute with one another and the Hamiltonian, so correspond to an extensive number of conserved quantities.This allows us to split the Hilbert space according to the eigenvalues ϕ p = ±1 (±i for odd plaquettes) of {W p }.
The Hamiltonian in eq. ( 1) can be exactly solved by transforming to a Majorana fermion representation [31], see fig. 1.Each spin is represented with four Majorana operators, σ α i = ib α i c i .We define a set of conserved bond operators ûjk = ib α j b α k .As with the W p operators, we may partition the Majorana Hilbert space according to the eigenvalues of these operators, u jk = ±1.For a given choice of these bond variables, eq. ( 1) reduces to a quadratic Majorana Hamiltonian where A jk = 2J α u jk .
The matrix iA determines properties of the fermionic degrees of freedom for a given flux configuration {u jk }.The spectrum is obtained by rotating to a new Majorana basis consisting of pairs of operators c′ j , c′′ j , defined by a matrix cj = R jk c k containing the fermionic eigenstates.The Hamiltonian takes the form H = j ε j ic ′ j c′′ j , and in what follows we refer to fermionic properties of the system as those determined by iA in a fixed flux sector.
The Kitaev Hamiltonian remains exactly solvable on any lattice in which no site connects to more than one bond of the same type [41].Thus, we shall restrict our investigation to lattices in which every vertex has coordination number z ≤ 3.Here we generate such lattices using Voronoi tessellation [58].Once a lattice has been generated, the bonds must be labelled in such a way that no vertex touches multiple edges of the same type, which we refer to as a three-edge colouring.The problem of finding such a colouring is equivalent to the classical problem of four-colouring the faces, which is always solvable on a planar graph [59,60] but can take up to seven colours on the torus [61].In practice, we reduce the colouring to a Boolean satisfiability problem [62] with details described in the supplemental material.One example of a coloured amorphous lattice is shown in fig.1(a).
Once the lattice and colouring has been found, the amorphous Hamiltonian is diagonalised using the same procedure as for the honeycomb model.Note that the Majorana system is only strictly equivalent to the initial spin system after a parity projection [63,64], details of which for the amorphous case are described in the supplemental material.Nevertheless, one can still use eq.( 2) to evaluate the expectation values of operators that conserve ûjk in the thermodynamic limit [65,66].The ground state energy of a given flux sector is the sum of the negative eigenvalues of iA/4 in eq. ( 2), and excitation energies are given by the positive eigenvalues.

III. RESULTS
We first investigate which flux patterns minimize the ground state energy on the amorphous lattice.When represented in the Majorana Hilbert space, flux operators W p = σ α j σ α k correspond to ordered products of link variables ûjk , and their eigenvalues describe the Z 2 flux through each plaquette, where the product is taken over the u jk values going clockwise around the border ∂p of each plaquette.We refer to a particular choice of a set of {ϕ p } as a flux sector.The spin Hamiltonian is real, thus it has TRS.However, the flux ϕ p through any plaquette with an odd number of sides has imaginary eigenvalues ±i.Thus, states with a fixed flux sector spontaneously break TRS, which in the context of crystalline Kitaev models was first described by Yao and Kivelson [67].All flux sectors come in degenerate pairs, where time reversal is equivalent to inverting the flux through every odd plaquette [43,46].For a system with n p plaquettes in periodic boundaries, there are 2 np−1 possible flux sectors, and in general it is a nontrivial task to determine which pair of flux sectors has the lowest energy.On the honeycomb lattice, the ground state was shown by Lieb to be flux free, ϕ p = +1 [56], however no such proof exists for amorphous lattices, since all lattice symmetries are broken.
To numerically determine the ground state flux sector, we first test a large number of finite size lattices (∼ 25, 000 lattices with 16 plaquettes), directly enumerating all possible flux configurations to find the lowest energy.In practice, care must be taken to account for finite size effects, as well as to ensure that the results hold as system size is increased -detailed in the supplemental material.Remarkably, we find that the energy is always minimised by setting the flux through each plaquette p to where n sides is the number of edges that form the plaquette and the global choice of the sign of i gives the two TRSdegenerate ground state flux sectors.The conjecture is consistent with results found on other regular lattices for which Lieb's theorem is not applicable [42].Having identified the ground state, any other flux sector can be characterized by the configuration of vortices, i.e. by the plaquettes whose flux is flipped with respect to ϕ g.s.
The ground state phase diagram can then be determined by varying the strength of each bond type, J α while remaining in the ground state flux sector, and we numerically calculate the ternary phase diagram shown in fig. 1 (c).The diagram contains two distinct phases: close to the corners of the triangle, e.g.|J z | ≫ |J x |, |J y |, the (A) phase is equivalent to the toric code on an amorphous lattice [68].The phase has a fermionic gap and supports Abelian excitations.Around the isotropic point J x = J y = J z , the (B) phase is also gapped in contrast to the honeycomb case as a consequence of TRS breaking from the finite density of odd plaquettes.All lattices studied in this work were generated from a voronoi lattice with completely random seed points, and so had on average equal proportions of odd and even plaquettes.We will confirm below that the (B) phase is indeed a chiral spin liquid.
As the values of J α are varied, the fermionic gap closes at the boundary between the two phases.In the honeycomb model, the phase boundaries are located on the straight lines for any permutation of α, β, γ ∈ {x, y, z}.We find that on the amorphous lattice these boundaries exhibit an inward curvature similar to honeycomb Kitaev models with flux [69] or bond [66] disorder.
Note, the presence of the gapped B-phase is non-trivial and related to our choice of homogeneous couplings for each colour of the bonds.In the supplemental material we study the robustness of the B phase with respect to bond disorder, e.g. a bond-length dependence of the interaction strength.In general one might expect disorder to lead to a gap closing, however we find that the gap is reduced but remains robust up to sizeable bond disorder.
A fundamental tool for understanding the distinction between the two phases is the Chern number.The original definition relies on momentum space, and so cannot be used here, where the system lacks any translational symmetry.However, recently methods have been developed for evaluating a real-space analogue of the Chern number [70,71].Here we shall use a slight modification of Kitaev's definition [7,31,57].For a choice of flux sector, we calculate the projector P onto the negative energy eigenstates of the matrix iA defined in eq. ( 2).The local Chern number around a point R in the bulk is given by where θ Rx is a step function in the x-direction, with the step located at x = R x , θ Ry is defined analogously.The trace is taken over a region around R in the bulk of the material, where care must be taken not to include any points close to the edges.Provided that the point R is sufficiently far from the edges, this quantity will be very close to quantised to the Chern number.Using this local Chern marker, we determine that the (A) phase has Chern number ν = 0, whereas the two TRSdegenerate ground state flux sectors in the (B) phase have Chern number ν = ±1 respectively.In closed boundaries, this leads to the appearance of gap-crossing protected edge modes, in accordance with the bulk-boundary correspondence [72], an example is shown in fig. 2. The edge modes are exponentially localised to the boundary of the system, and can be qualitatively distinguished from bulk states by their large inverse participation ratio, where ψ(r) denotes an eigenmode of the free Majorana Hamiltonian derived in eq. ( 2).Finally, we note that the closing of the fermionic gap on the boundary between the two phases is necessary in order to transition between states with different Chern numbers.
Anderson Transition to a Thermal Metal -Having understood the spontaneous formation of a chiral amorphous QSL ground state, we are now in a position to discuss the finite temperature behavior of the model.In general, an Ising-like thermal phase transition into the chiral QSL phase is expected akin to the one observed for the Yao-Kivelson model [73] but a full Monte-Carlo sampling, which is further complicated by the inherent disorder in the amorphous lattice, is beyond the scope of this letter.Nevertheless, the main effect of increasing temperature is the proliferation of fluxes which allow us to gain a qualitative understanding of the finite temperature behavior [69].
On the honeycomb Kitaev model with explicit TRS breaking, Majorana zero modes bind to fluxes forming Ising non-Abelian anyons [74].Their pairwise interaction decays exponentially with separation [48,49,75].As temperature is increased, the proliferation of vortices in the system produces a finite density of anyons and their hybridization leads to an Anderson transition to a macroscopically degenerate state known as a thermal metal phase [47,48,76].This exotic phase has two key signatures.Firstly, the metallic phase is defined by a closing of the fermion gap -that is, it is driven by vortex configurations with a gapless fermionic spectrum.Secondly, we expect the density of states in a thermal metal to diverge logarithmically with energy and display characteristic low energy oscillations predicted by random matrix theory [47,77].In the supplemental material we present numerical evidence showing that all of the above features carry over to the amorphous QSL with spontaneous TRS breaking, giving strong evidence for the transition to the thermal metal phase.

IV. DISCUSSION
We have studied an extension of the Kitaev honeycomb model to amorphous lattices with coordination number z = 3.We found that it is able to support two quantum spin liquid phases that can be distinguished using a realspace generalisation of the Chern number.The presence of odd-sided plaquettes results in a spontaneous breaking of TRS, and the emergence of a chiral spin liquid phase.Furthermore we found evidence that the amorphous system undergoes an Anderson transition to a thermal metal phase, driven by the proliferation of vortices with increasing temperature.Our exactly soluble chiral QSL provides a first example of a topologically quantum many-body phase in amorphous magnets, which raises a number of questions for future research.
First, a numerically challenging task would be a study of the full finite temperature phase diagram via Monte-Carlo sampling and possible violations of the Harris criterion for the Ising transition stemming from the inherent lattice disorder [78][79][80].Second, it would be worthwhile to search for experimental realisations of amorphous Kitaev materials, which can possibly be created from crystalline ones using standards method of repeated liquifying and fast cooling cycles [3,21,23].The putative QSL behavior of the intercalated Kitaev compound H 3 LiIr 2 O 6 [81,82] could possibly be related to amorphous lattice disorder.Moreover, metal organic frameworks are promising platforms forming amorphous lattices [83] with recent proposals for realizing strong Kitaev interactions [84] as well as reports of QSL behavior [85].We expect that an experimental signature of a chiral amorphous QSL is a half-quantized thermal Hall effect similar to magnetic field induced behavior of honeycomb Kitaev materials [86][87][88][89].Alternatively, it could be characterized by local probes such as spin-polarized STM [90][91][92] and the thermal metal phase displays characteristic longitudinal heat transport signatures [74].Third, it would be interesting to study the stability of the chiral amorphous Kitaev QSL with respect to perturbations [93][94][95][96][97] and, importantly, to investigate whether QSL may exist for spin-isotropic Heisenberg models on amorphous lattices.
Overall, there has been surprisingly little research on amorphous quantum many body phases albeit material candidates aplenty.We expect our exact chiral amorphous spin liquid to find many generalisation to realistic amorphous quantum magnets and beyond.

Appendix A: Lattice Generation
We generate tri-coordinate lattices by taking the Voronoi partition of the unit square with respect to uniformly sampled seed points [58].This partitions the space into polyhedral volumes enclosing the region closest to each seed point.In two dimensions, the vertices and edges of these polygons form a tri-coordinate lattice, exactly what is necessary for the Kitaev model.To produce lattices with periodic boundary conditions we first tile the seed points into a repeating lattice.The Voronoi partition of the tiled seed points can then be converted into a lattice embedded onto the torus by connecting corresponding edges that cross the unit square boundaries, see fig. 3. Finally, to somewhat regularise bond lengths in our lattice, a single step of Lloyd's algorithm is performed, where every vertex is shifted to the center of mass of the three plaquettes that it touches.This is done to improve the readability of the lattice, reducing the number of extremely short bonds appearing, and has no effect on the physics.
Once a lattice has been generated, the bonds must be labelled in such a way that no vertex touches multiple edges of the same type, which we refer to as a three-edge colouring.The problem of finding such a colouring is equivalent to the classical problem of four-colouring the faces, which is always solvable on a planar graph [59,60].On the torus, a face colouring can require up to seven colours [61], so not all graphs can be assumed to be 3-edge colourable.However, such exceptions seem rare for graphs generated via voronisation -every graph generated in this study admitted multiple distinct 3-edge colourings.In practice, the problem of finding a colouring for a given graph can be reduced to a Boolean satisfiability problem [62], which we then solve using the open-source solver MiniSAT [99].
Care must be taken in the definition of open boundary conditions, simply removing bonds from the lattice leaves behind unpaired b α operators that need to be paired in some way to arrive at fermionic modes.In order to fix a pairing we always start from a lattice defined on the torus and generate a lattice with open boundary conditions by defining the bond coupling J α ij = 0 for sites joined by bonds (i, j) that we want to remove.This creates fermionic zero modes u ij associated with these cut bonds which we set to 1 when calculating the projector.All our code is available online [98].

Appendix B: The Projector
Closely following the derivation of [63] we can extend the projector from Majorana wavefunctions to physical spin states to the amorphous case.In the standard way, we define normal mode operators from there we form fermionic operators and their associated number operators n i = f † i f i .The many body ground state within a vortex sector is then defined by the set of occupation numbers n m = 0, 1. Lastly we need to define the fermion parity π = i N (1 − 2n i ).
The projector can be written as where D i are the local projectors.S symmetrises over gauge equivalent states while P 0 is responsible for annihilating unphysical states, see [63] for details.
To extend this to the amorphous case we calculate the product of the local projectors The Voronoi partition (lines) splits a region up into polyhedra closer to one of the seed points (points) than any other.In two dimensions this yields a tri-coordinate lattice.Dotted lines go off to infinity.(b) To create a lattice on the torus, we tile the seed points into a 3x3 grid and compute a Voronoi partition.By identifying pairs of edges (dotted lines) that cross the unit square (in grey) as the same we turn this lattice into on defined on the torus.(c) The final tri-coordinate lattice in periodic boundary conditions, coloured such that all three colours meet at every vertex.
for a tri-coordinate lattice with N faces, 2N vertices and 3N edges.The operators can be ordered by bond type without utilising any property of the lattice.
The product over c i operators reduces to a determinant of the Q matrix and the fermion parity.The only problem is to compute the factors p x , p y , p z = ±1 that arise from reordering the b operators such that pairs of vertices linked by the corresponding bonds are adjacent.
This is simple the parity of the permutation from one ordering to the other and can be computed quickly with a cycle decomposition.The final form is almost identical to the honeycomb case with the addition of the lattice structure factors p x , p y , p z where det(Q u ) and u ij depend on the lattice and the particular vortex sector.

Appendix C: Numerical Evidence for the Ground State Flux Sector
In this section we detail the numerical evidence collected to support the claim that, for an arbitrary lattice, a gapped ground state flux sector is found by setting the flux through each plaquette to ϕ g.s.= −(±i) n sides .This was done by generating a large number (∼ 25,000) of lattices and exhaustively checking every possible flux sector to find the configuration with the lowest energy.We checked both the isotropic point (J α = 1), as well as in the toric code (J x = J y = 0.25, J z = 1).
The argument has one complication: for a graph with n p plaquettes, there are 2 np−1 distinct flux sectors to search over, with an added factor of 4 when the global fluxes Φ x and Φ y wrapping around the cylinder directions are taken into account.Note that the −1 appears in this counting because fluxes can only be flipped in pairs.To be able to search over the entire flux space, one is necessarily restricted to looking at small system sizeswe were able to check all flux sectors for systems with n p ≤ 16 in a reasonable amount of time.However, at such small system size we find that finite size effects are substantial.In order to overcome these effects we tile the system and use Bloch's theorem (a trick that we shall refer to as twist-averaging for reasons that shall become clear) to efficiently find the energy of a much larger (but periodic) lattice.Thus we are able to suppress finite size effects, at the expense of losing long-range disorder in the lattice.
Our argument has three parts: First we shall detail the techniques used to exhaustively search the flux space for a given lattice.Next, we discuss finite-size effects and explain the way that our methods are modified by the twist-averaging procedure.Finally, we demonstrate that as the size of the disordered system is increased, the effect of twist-averaging becomes negligible -suggesting that our conclusions still apply in the case of large disordered lattices.
Testing All Flux Sectors -For a given lattice and flux sector, defined by {u jk }, the fermionic ground state The energy of every flux sector explored for a system of 16 plaquettes, the order is arbitrary.The two ground state flux sectors can be identified as the points with lowest energy.(b) The fermion gap for each of the flux sectors explored.Note that the largest fermion gap coincides with the ground state flux sector.This occurred in ∼ 85% of cases tested.(c) Average energy of the systems tested over a range of system sizes from np = 9 to np = 1600.The region between the upper and lower quartiles is shown in red, and the full range of energies obtained is shown in orange.(d) Average fermion gap as a function of system size.Again, the region between the upper and lower quartiles is shown in red, and the full range is in orange.As can be seen, no gapless systems were found for np > 20.
energy is calculated by taking the sum of the negative eigenvalues of the matrix The set of bond variables u jk , which we are free to choose, determine the Z 2 gauge field.However only the fluxes, defined for each plaquette according to have any effect on the energies.Thus, there is enormous degeneracy in the u jk degrees of freedom.Flipping the bonds along any closed loop on the dual lattice has no effect on the fluxes, since each plaquette has had an even number of its constituent bonds flipped -as is shown in the following diagram: where the flipped bonds are shown in red.In order to explore every possible flux sector using the u jk variables, we restrict ourselves to change only a subset of the bonds in the system.In particular, we construct a spanning tree on the dual lattice, which passes through every plaquette in the system, but contains no loops.
The tree contains n p − 1 edges, shown in red, whose configuration space has a 1 : 1 mapping onto the 2 np−1 distinct flux sectors.Each flux sector can be created in precisely one way by flipping edges only on the tree (provided all other bond variables not on the tree remain fixed).Thus, all possible flux sectors can be accessed by iterating over all configurations of edges on this spanning tree.
Finite Size Effects -In our numerical investigation, the objective was to test as many example lattices as possible.We aim for the largest lattice size that could be efficiently solved, requiring a balance between lattice size and cases tested.Each added plaquette doubles the number of flux sectors that must be checked.25,000 lattices containing 16 plaquettes were used.However, in his numerical investigation of the honeycomb model, Kitaev demonstrated that finite size effects persist up to much larger lattice sizes than we were able to access [31].
In order to circumvent this problem, we treat the 16plaquette amorphous lattice as a unit cell in an arbitrarily large periodic system.The bonds that originally connected across the periodic boundaries now connect adjacent unit cells.This infinite periodic Hamiltonian can then be solved using Bloch's theorem, since the larger system is diagonalised by a plane wave ansatz.For a given crystal momentum q ∈ [0, 2π) 2 , we are left with a Bloch Hamiltonian, is identical to the original Hamiltonian aside from an extra phase on edges that cross the periodic boundaries in the x and y directions, where q jk = q x for a bond that crosses the x-periodic boundary in the positive direction, with the analogous definition for y-crossing bonds.We also have q jk = −q kj .Finally q jk = 0 if the edge does not cross any boundaries at all -in essence we are imposing twisted boundary conditions on our system.The total energy of the tiled system can be calculated by summing the energy of M (q) for every value of q.In practice we constructed a lattice of 50 × 50 values of q spanning the Brillouin zone.The procedure is called twist averaging because the energyper-unit cell is equivalent to the average energy over the full range of twisted boundary conditions.Evidence for the Ground State Ansatz For each lattice with 16 plaquettes, 2 15 = 32,768 flux sectors are generated.In each case we find the energy (averaged over all twist values) and the size of the fermion gap, which is defined as the lowest energy excitation for any value of q.We then check if the lowest energy flux sector aligns with our ansatz (given by ϕ g.s.p = −(±i) n sides ) and whether this flux sector is gapped.
In the isotropic case (J α = 1), all 25,000 examples conformed to our guess for the ground state flux sector.A tiny minority (∼ 10) of the systems were found to be gapless.As we shall see shortly, the proportion of gapless systems vanishes as we increase the size of the amorphous lattice.An example of the energies and gaps for one of the systems tested is shown in fig. 4 (a).For the anisotropic phase (we used J x , J y = 0.25, J z = 1) the overwhelming majority of cases adhered to our ansatz, however a small minority (∼ 0.5%) did not.In these cases, however, the energy difference between our ansatz and the ground state was at most of order 10 −6 .Further investigation would need to be undertaken to determine whether these anomalous systems are a finite size effect due to the small amorphous system sizes used or a genuine feature of the toric code phase on such lattices.
A Gapped Ground State -Now that we have collected sufficient evidence to support our guess for the ground state flux sector, we turn our attention to checking that this sector is gapped.We no longer need to exhaustively search over flux space for the ground state, so it is possible to go to much larger system size.We generate 40 sets of systems with plaquette numbers ranging from 9 to 1600.For each system size, 1000 distinct lattices are generated and the energy and gap size are calculated without phase twisting, since the effect is negligible for such large system sizes.As can be seen in fig. 4 (c-d), for very small system size a small minority of gapless systems appear, however beyond around 20 plaquettes all systems had a stable fermion gap in the ground state.Finally, the energy and gap difference between the phase twisted and non-phase twisted results vanished exponentially with the system size, supporting the claim that the results can be straightforwardly extrapolated to large systems.At small ρ, the states populating the gap possess τ ≈ 0, indicating that they are localised states pinned to individual fluxes -the system remains insulating.At larger ρ, the in-gap states merge with the bulk band and become extensive, closing the gap -the system transitions to a metallic phase.
Finally, the averaged density of states in the ρ = 0.5 case is shown in fig.6 (b) for both the Honeycomb model and our amorphous lattice.Note that only the honeycomb model is calculated in the presence of an effective magnetic field explicitly breaking TRS.In both cases we see the logarithmic scaling alongside the characteristic RMT oscillations at low energy, giving strong evidence that the amorphous model displays a finite temperature transition to a thermal metal phase.

FIG. 1 .
FIG. 1. Construction details for the amorphous lattice model.(a) Amorphous lattice generated via Voronoi tessellation of a uniformly distributed random point set on the unit square.Periodic boundary conditions are imposed by tiling the unit square before Voronoi tessellation.(b) Magnified portion of the amorphous lattice.Arrows from site j to site k indicate direction where the bond variable u jk = 1.An arbitrary flux sector is shown, where shaded plaquettes have Z2 flux flipped with respect to the ground state.Colours correspond to a valid assignment of the bond colourings, α jk .The inset demonstrates the Majorana construction on a tri-coordinate motif, which allows for the exact solution of the model.(c) Ternary phase diagram of the amorphous Kitaev model with varying exchange coupling.The isotropic regime |Jx| ≈ |Jy| ≈ |Jz| (B), exhibits a topologically non-trivial chiral QSL ground state with Chern number ν = ±1.The fermion gap of the ground state flux sector closes at the phase boundary (solid black lines), and a transition occurs to a ν = 0 phase (A) for anisotropic couplings.The phase boundary was obtained by averaging over 20 amorphous lattice realisations with ∼ 400 sites.Dotted black lines indicate the corresponding phase boundaries in the honeycomb model.
FIG. 2. Ground-state flux sector wavefunctions and spectrum.(a) In-gap fermionic wavefunction drawn from the ground state flux sector in open boundary conditions, showing a topological edge mode.The number density for this state along a line of lattice sites spanning the system (black line) is shown in the bottom subfigure on a logarithmic scale, demonstrating the characteristic exponential decay of topological edge modes with distance from the edge.(b) Ground-state flux sector fermionic density of states in open boundary conditions, coloured by inverse participation ratio.The increased inverse participation ratio of the in-gap states signifies their localisation to the edges of the system.

FIG. 4 .
FIG. 4. (a)The energy of every flux sector explored for a system of 16 plaquettes, the order is arbitrary.The two ground state flux sectors can be identified as the points with lowest energy.(b) The fermion gap for each of the flux sectors explored.Note that the largest fermion gap coincides with the ground state flux sector.This occurred in ∼ 85% of cases tested.(c) Average energy of the systems tested over a range of system sizes from np = 9 to np = 1600.The region between the upper and lower quartiles is shown in red, and the full range of energies obtained is shown in orange.(d) Average fermion gap as a function of system size.Again, the region between the upper and lower quartiles is shown in red, and the full range is in orange.As can be seen, no gapless systems were found for np > 20.

FIG. 6 .
FIG.6.(a) Density of states (top) and inverse participation ratio scaling exponent (bottom) of the fermionic spectrum as a function of flux defect density, ρ, for isotropic couplings.Each pixel is averaged over 10 independent lattice realisations, in flux sectors sampled from an ensemble with a proportion ρ of fluxes flipped with respect to the ground state sector.White pixels correspond to bins containing no fermionic states.At low defect density, the fermionic spectra are gapped.As the defect density increases, in-gap states appear with a small τ , indicating that they are strongly localized around defects.At large defect density, τ increases for the in-gap sites, indicating that they are delocalised and the system becomes gapless.Using defect density as a proxy for temperature, this demonstrates the thermal phase transition from a chiral spin liquid to a thermal metal phase.(b) A histogram of fermionic density of states sampled from the thermodynamic ensemble of flux sectors for T → ∞, i.e. all gauge configurations equally likely.The oscillations at low E are characteristic of a thermal metal phase[47], demonstrated for the Kitaev honeycomb lattice model subject to a magnetic field (top) and the amorphous Kitaev model (bottom).L corresponds to the linear extent of the system -L ∼ √ N with N sites -for both lattice types.