Methods
The practical implemntation of what is described in this section is available as a Python package called Koala (Kitaev On Amorphous LAttices)1 most of the figures shown were generated with Koala.
Voronisation
In order to study the properties of the amorphous Kitaev model we need a way to sample from the space of possible trivalent graphs.
A very simple way to do this is to use a Voronoi partition of the torus2–4. We start by sampling seed points uniformly (or otherwise) on the torus. We then compute the partition of the torus into regions closest (with a Euclidean metric) to each seed point. The straight lines (if the torus is flattened out) at the borders of these regions become the edges of the new lattice and the points where they intersect beceme the vertices.
The graph generated by a Voronoi partition of a two dimensional surface is always planar meaning that no edges cross eachother when the graph is embedded into the plane. It is also trivalent in the sense that every vertex is connected to exactly three edges cite.
Ideally we might instead sample uniformly from the space of possible trivalent graphs, and indeed there has been some work on how to do this using a Markov Chain Monte Carlo approach5, however it does not gurantee that the resulting graph is planar which we will need to ensure that the edges can be 3-coloured.
In practice, we then use a standard algorithm6 from scipy7 which actually computes the Voronoi partition of the plane. In order to compute the Voronoi partition of the torus, I take the seed points and replicate them into a repeating grid, either 3x3 (or for very small numbers of seed points 5x5). I then identify edges in the output to construct a lattice on the torus.
Graph Representation
We represent the graph structure with an ordered list of edges \((i,j)\) so we can represent both directed and undirected graphs which is useful for defining the sign of bond operators \(u_{ij} = - u_{ji}\).
Coloring the Bonds
The Kitaev model requires that each edge in the lattice be assigned a label \(x\), \(y\) or \(z\) such that each vertex has exactly one edge of each type connected to it. Let \(\Delta\) be the maximum degree of a graph which in our case is 3. If \(\Delta > 3\) it is obviously not possible to 3 color the edges but the general theory of when this is and isn’t possible for graphs with \(\Delta \leq 3\) is more subtle.
In the graph theory literature, graphs where all vertices have degree 3 are commonly called cubic graphs, there is no term for graphs with maximum degree 3. Planar graphs are those that can be embedded onto the plane without any edges crossing. Bridgeless graphs do not contain any edges that, when removed, would partition the graph into disconnected components.
It’s important to be clear that this problem is different from that considered by the famous 4 color theorem8 . The 4 color thorem is concerned with assiging colours to the vertices of a graph such that no vertices that share an edge are the same colour. Here we are concerned with an edge colouring.
The four color theorem applies to planar graphs, those that can be embedded onto the plane without any edges crossing. Here we are actually concerned with Toroidal graphs which can be embedded onto the torus without any edges crossing. In fact toroidal graphs require up to 7 colors9 . The complete graph \(K_7\) is a good example of a toroidal graph that requires 7 colours.
\(\Delta + 1\) colours are enough to edge-colour any graph and there is an \(\mathcal{O}(mn)\) algorithm to do it for a graph with \(m\) edges and \(n\) vertices10. Restricting ourselves to graphs with \(\Delta = 3\) like ours, those can be 4-edge-coloured in linear time11 .
It’s trickier if we want to 3-edge-colour them however. Cubic, planar bridgeless graphs can be 3-edge-coloured if and only if they can be 4-face-coloured12 . For which there is an \(\mathcal{O}(n^2)\) algorithm robertson1996efficiently . However it is not clear whether this extends to cubic, toroidal bridgeless graphs.
4-face-colourablity implies 3-edge-colourability
The proof of that 4-face-colourablity implies 3-edge-colourability can be sketched out quite easily: 1. Assume the faces of G can be 4-coloured with labels (0,1,2,3) 2. Label each edge of G according to \(i + j \mathrm{mod} 3\) where i and j are the labels of the face adjacent to that edge. For each edge label there are two face label pairs that do not share any face labels. i,e the edge label \(0\) can come about either from faces \(0 + 3\) or \(1 + 2\).
\[\begin{aligned} 0 + 3 \;\mathrm{or}\; 1 + 2 &= 0 \;\mathrm{mod}\; 3\\ 0 + 1 \;\mathrm{or}\; 2 + 3 &= 1 \;\mathrm{mod}\; 3\\ 0 + 2 \;\mathrm{or}\;1 + 3 &= 2 \;\mathrm{mod}\; 3\\ \end{aligned} \]
- In a cubic planar G, a vertex v in G is always part of 3 faces and the colors of those faces determines the colors of the edges that connect to v. The three faces must take three distinct colors from (0,1,2,3).
- From there’s easy to convince yourself that those three distinct face colours can never produce repeated edge colours according to the \(i+j \;\mathrm{mod}\; 3\) rule.
This implies that all cubic planar graphs are 3-edge-colourable. It does not apply to toroidcal graphs, however I have not yet generated a voronoi lattices on the torus that is not 3-edge-colourable. This suggests that perhaps voronoi lattices have additional structure that makes them 3-edge-colourable. Intuitively, the kinds of toroidal graphs that cannot be 3-edge-coloured look as if they could never be generated by a voronoi partition with more than a few seed points.
Finding Lattice colourings in practice (unfinished)
Some things are harder in theory than in practice. 3-edge-colouring cubic toroidal graphs appears to be one of those things.
The approach I take is relatively standard in the computer science community for solving NP problems computationally. I don’t believe this problem to be in NP but I tried it anyway.
The trick is to map the problem on into a Boolean Satisfiability
‘SAT’ problem13, use an
off the shelf solver, MiniSAT14, and finally to map the problem
back to the original domain. While SAT solvers are very general, they
are also highly optimised and they do seem to yield good results for
this problem.
SAT solvers encode problems as constraints on some number of boolean variables \(x_i \in {0,1}\). The constraints must Conjunctive Normal Form (CNF). CNF means the constraints are encoded as a set of clauses of the form \[x_1 \;\textrm{or}\; \bar{x}_3 \;\textrm{or}\; x_5\] that containt logical ORs of some subset of the variables where any of the variables may also be logical NOT’d which I represent by over bars here.
A solution of the problem is one that makes all the clauses simultaneously true.
I encode the edge colouring problem as a set of statements about a set of boolean variables \(x_i \in {0,1}\). For \(B\) bonds we take the \(3B\) variables \(x_{i\alpha}\) where \(x_{i\alpha} = 1\) indicates that edge \(i\) has colour \(\alpha\).
For edge colouring graphs we need two kinds of constraints: 1. Each edge is exactly one colour. 2. No neighbouring edges are the same color.
The first constraint is a kind of artifact of doing this mapping over to boolean variables, the solver doesn’t know anything about the structure of the problem unless it is encoded into the variables.
The second constraint encodes the structure of the graph itself and can be constructed easily from the adjacency matrix.
I’ll fill in the encoding later but the gist is that we can give this to a solver and get back: whether the problem is solveable, a solution or all the possible solutions. Finding a solution is relatively fast, while finding all the solutions is slower since there appear to be exponentially many of them. Fig \(\ref{fig:multiple_colourings}\) shows some examples.
Mapping between flux sectors and bond sectors
Constructing the Majorana representation of the model requires the particular bond configuration \(u_{jk} = \pm 1\). However the large number of gauge symmetries of the bond sector make it unwieldly to work with. We therefore need a way to quickly map between bond sectors and flux sectors.
Going from the bond sector to flux sector is easy since we can compute it directly by taking the product of \(i u_{jk}\) around each plaquette \[ \phi_i = \prod_{(j,k) \; \in \; \partial \phi_i} i u_{jk}\]
Going from flux sector to bond sector requires more thought however. The algorithm I use is this:
Fix the gauge by choosing some arbitrary \(u_{jk}\) configuration. In practice I use \(u_{jk} = +1\). This chooses an arbitrary one of the 4 topological sectors.
Compute the current flux configuration and how it differs from the target one. Let’s call an plaquette that differs from the target a defect.
Find any adjacent pairs of defects and flip the \(u_jk\) between them. This leaves a set of isolated defects.
Pair the defects up using a greedy algorithm.
Compute paths along the dual lattice between each pair of plaquettes. Flipping the corresponding set of \(u_{jk}\) transports one flux to the other and anhilates them.
Amorphous materials are glassy condensed matter systems characterised by short-range constraints in the absence of long-range crystalline order as first studied in amorphous semiconductors15,16. In general, the bonds of a whole range of covalent compounds 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 symmetry17,18, the most famous example being amorphous Ge and Si with \(z=4\)19,20. 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 invariants2,3,21–25. However, research on electronic systems has been mostly focused on non-interacting systems with a few notable exceptions for understanding the occurrence of superconductivity26–28 in amorphous materials and recently the effect of strong repulsion in amorphous TIs29.
Magnetic phases in amorphous systems have been investigated since the 1960s, mostly through the adaptation of theoretical tools developed for disordered systems and numerical methods~. Research focused on classical Heisenberg and Ising models which have been shown to account for observed behavior of ferromagnetism, disordered antiferromagnetism and widely observed spin glass behaviour~. However, the role of spin-anisotropic interactions and quantum effects has not been addressed. Similarly, it is an open question whether magnetic frustration in amorphous quantum magnets can give rise to long-range entangled quantum spin liquid (QSL) phases.
%Broad constraints to the possible phases hosted by Heisenberg amorphous magnets were provided by the phenomenological theory developed by Andreev and Marchenko30–32. The phases in this theory are described by a set of macroscopic magnetic vectors that transform according to the irreducible representations of the group of spatial symmetries of the system30. Amorphous magnets are treated, on average, as homogeneous and isotropic, being thus symmetric under three-dimensional rotations and spatial inversion31. Only three types of phases are consistent to this group of symmetries, corresponding to ferromagnets, disordered antiferromagnets, or spin glasses31,32.
Two intentional simplifications of Andreev’s and Marchenko’s theory were the neglect of spin-orbit coupling induced anisotropies and the effects arising from the local structure of amorphous lattices. It is then expected that their theory is invalid for amorphous compounds generated from crystalline magnets with strong spin-orbit coupling with tight geometrical arrangements. Several instances of these magnets were synthesized in the last decade, among which we highlight the Kitaev materials33–37. It was suggested (and later observed) that heavy-ion Mott insulators formed by edge-sharing octahedra could be good platforms for the celebrated Kitaev model on the honeycomb lattice33, an exactly solvable model whose ground state is a quantum spin liquid (QSL)38–41 characterized by a static \(\mathbb Z_2\) gauge field and Majorana fermion excitations42. The model displays bond-dependent Ising-like exchanges that give rise to local symmetries, which are essential to its mapping onto a free fermion problem43,44. Such a mapping is rigorously extendable to any three-coordinated graph in two or three dimensions satisfying a simple geometrical condition45–48. Thus, it reasonable to suppose that the Kitaev model is also analytically treatable on certain amorphous lattices, therefore becoming a realistic starting point to study the overlooked possibility of QSLs in amorphous magnets.
In this letter, we study Kitaev spin liquids (KSLs) stabilized by the \(S=1/2\) Kitaev model on coordination number \(z=3\) random networks generated via Voronoi tessellation . On these lattices, the KSLs generically break time-reversal symmetry (TRS), as expected for any Majorana QSL in graphs containing odd-sided plaquettes . An extensive numerical study showed that the \(\mathbb Z_2\) gauge fluxes on the ground state can be described by a conjecture consistent with Lieb’s theorem . In contrast to the honeycomb case, the amorphous KSLs are gapless only along certain critical lines. These manifolds separate two gapped KSLs that are topologically differentiated by a local Chern number \(\nu\) in analogy with the KSLs on the decorated honeycomb lattice . The \(\nu=0\) phase is the amorphous analogue of the abelian toric-code QSL , whereas the \(\nu=\pm1\) KSLs is a non-Abelian chiral spin liquid (CSL). We study two specific features of the latter liquid: topologically protected edge states and a thermal-induced Anderson transition to a thermal metal phase .% The Kitaev spin liquids are classified by their Majorana fermion dispersion and topological properties. On the honeycomb lattice, tuning the exchange couplings \(J^\alpha\) can change the ground state from a gapped QSL with Abelian anyonic excitations (e.g., when \(J^z\gg J^x,J^y\)) or gapless (e.g., when \(J^z=J^x=J^y\)). In the latter case, breaking time reversal symmetry (TRS) opens a gap that signals the onset of a chiral spin liquid (CSL) phase supporting non-Abelian excitations and protected edge modes. On the honeycomb lattice, CSLs are only obtained by perturbing a Hamiltonian with, for example, magnetic fields or Dzyaloshinskii-Moriya exchanges . CSLs on the pure Kitaev model can be obtained on \(z=3\) lattices containing odd-sided plaquettes, for which any Majorana QSL displays spontaneous TRS breaking , as confirmed on decorated honeycomb and non-Archimedean lattices.
{} We start with a brief review of the Kitaev model on the honeycomb lattice . Here, a spin-1/2 is placed on every vertex and each bond is labelled by an index \(\alpha \in \{ x, y, z\}\). The bonds are arranged such that each vertex connects to exactly one bond of each type. The Hamiltonian is given by \[\begin{equation} \label{eqn:kitham} \mathcal{H} = - \sum_{\langle j,k\rangle_\alpha} J^{\alpha}\sigma_j^{\alpha}\sigma_k^{\alpha}, \end{equation}\] where \(\sigma^\alpha_j\) is a Pauli matrix acting on site \(j\), (j,k_) is a pair of nearest-neighbour indices connected by an \(\alpha\)-bond with exchange coupling \(J^\alpha\). For each plaquette of the lattice, we can define the conserved operator $ W_p = _j^{}_k^{}$, where the product runs clockwise over the bonds around the plaquette. This provides an extensive number of conserved plaquettes that allow us to split the Hilbert space in terms of the eigenvalues of \(W_p\). The KSL is uncovered by transforming eqn.~\(\ref{eqn:kitham}\) to a four-Majorana representation of the spin operators, \(\sigma_i^\alpha = i b_i^\alpha c_i\) , where the Hamiltonian takes the form \[\begin{equation}\label{eqn:majorana_hamiltonian} \mathcal{H} = \frac{i}{4}\sum_{j,k}A_{jk}^{(\alpha)}c_jc_k. \end{equation}\] Here, \(A_{jk}^{(\alpha)}=2J^{\alpha}u_{jk}\) with \(\hat u_{jk} = ib_j^{\alpha}b_k^{\alpha}\) being conserved \(\mathbb Z_2\) bond operators. Once the \(\hat u_{jk}\) eigenvalues are fixed, the Kitaev model becomes equivalent to a fermionic problem that can be diagonalized with standard methods . The Kitaev Hamiltonian remains exactly solvable on any graph in which no site connects to more than one bond of the same type . Thus, we are restricted to lattices in which every vertex has coordination number \(z \leq 3\). Here, such graphs are generated with Voronoi tessellation . A set of points are sampled uniformly from the unit square and cells are generated as the region of space closer to a given point than any other. The lattice is given by the boundaries between cells with edges at the interface of two cells and vertices at the point where three edges meet. Periodic boundary conditions are imposed by tiling the initial set of points and then connecting corresponding edges that cross the unit square boundaries - see for technical details. One example of such an amorphous lattice is shown in~(a). Once a random network has been generated, the bonds types must be assigned in a way that is consistent with our condition, which we refer to as a . The problem of finding such a colouring was shown to be equivalent to the classical problem of four-colouring the faces, which is always solvable in planar graphs~. On the torus, a face colouring can require up to seven colours , and so not all graphs can be assumed to be 3-edge colourable. However, such exceptions are rare – every graph generated in this study admitted multiple distinct 3-edge colourings. The problem of finding a colouring for a given graph can be reduced to a Boolean satisfiability problem , which we then solve using the open-source solver ~.Once the three-edge colouring has been found, the Kitaev Hamiltonian is mapped onto eqn.~\(\ref{eqn:majorana_hamiltonian}\), which corresponds to the spin fractionalization in terms of a static \(\mathbb Z_2\) gauge fields and \(c\) matter as indicated in ~(b) . Strictly speaking, the Majorana system is equivalent to the original spin system after applying a projector operator , whose form is presented in . Despite this caveat, one can still use eqn.~\(\ref{eqn:majorana_hamiltonian}\) to evaluate the expectation values of operators conserving \(\hat u_{jk}\) in the thermodynamic limit . This type of operator is exemplified by the Hamiltonian itself, for which the ground state energy of a fixed sector is the sum of the negative eigenvalues of \(iA/4\) in eqn.~\(\ref{eqn:majorana_hamiltonian}\), and whose excitations are extracted from the positive eigenvalues of the same matrix.
{} Let us now consider the conserved operators $ W_p = _j^{}_k^{}$ on amorphous lattices. When represented in the Majorana Hilbert space, these operators correspond to ordered products of \(\hat u_{jk}\), and their fixed eigenvalues are written as \[\begin{equation} \label{eqn:flux_definition} \phi_p = \prod_{(j,k) \in \partial p} (-iu_{jk}), \end{equation}\] where the pairs \(j,k\) are crossed around the border \(\partial p\) of the plaquette on the orientation. In periodic boundaries there is an additional pair of global \(\mathbb{Z}_2\) fluxes \(\Phi_x\) and \(\Phi_y\), which are calculated along an arbitrary closed path that wraps the torus in the \(x\) and \(y\) directions respectively. The energy difference between distinct flux sectors decays exponentially with system size, so that the ground state of any flux sector in the thermodynamic limit displays a fourfold topological degeneracy . We now need to determine the ground state flux sectors. First, let us recall that Majorana QSLs emerging on graphs containing odd-sided plaquette undergo a spontaneous TRS breaking . Therefore, there will be always a twofold ground state degeneracy due to time-reversal, in which one ground state is related to the other by inversion of imaginary \(\phi_p\) fluxes . An insight pointing to the ground state sectors come from the model on the honeycomb lattice, for which a theorem proved by Lieb sets that the ground state sector to be \(\phi_p=+1\), \(\forall p\) . Although Lieb’s theorem is not extendable to amorphous lattices, it is suggested the ground state energy for a sufficiently large system is minimised by setting \[\begin{align} \label{eqn:gnd_flux} \phi_p^{\textup{g.s.}} = -(\pm i)^{n_{\textup{sides}}}, \end{align}\] where \(n_{\textup{sides}}\) is the number of edges that form \(p\) and the global choice of the sign of \(i\) gives each of the two TRS-degenerate ground state flux sectors. Such a conjecture is consistent with Lieb’s theorem on regular lattices and is supported by numerical evidence as detailed in . Once we have identified the ground state, any other sector can be characterized by the configuration of vortices, i.e. by the plaquettes whose flux is flipped with respect to \(\left\{\phi_p^{\textup{g.s.}}\right\}\). {} We numerically found that the amorphous KSLs are generally gapped, except along the critical lines displayed in \(\ref{fig:example_lattice}\)(c). The QSLs separated by these lines are distinguished by a real-space analogue of the Chern number . A similar topological number was discussed by Kitaev on the honeycomb lattice that we shall use here with a slight modification . For a choice of flux sector, we calculate the projector \(P\) onto the negative energy eigenstates of the matrix \(iA\) defined in eqn.~\(\ref{eqn:majorana_hamiltonian}\). The local Chern number around a point \(\bf R\) in the bulk is given by \[\begin{align} \nu (\bf R) = 4\pi \Im \Tr_{\textup{Bulk}} \left ( P\theta_{R_x} P \theta_{R_y} P \right ), \end{align}\] where \(\theta_{R_x}\) is a step function in the \(x\)-direction, with the step located at \(x = R_x\), \(\theta_{R_y}\) is defined analogously. The trace is taken over a region around \(\bf 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 \(\bf R\) is sufficiently far from the edges, this quantity will be very close to quantised to the Chern number. The local Chern marker distinguishes between an Abelian phase (A) with \(\nu = 0\), and a non-Abelian (B) phase characterized by \(\nu = \pm 1\). The (A) phase is equivalent to the toric code on an amorphous system . { Since the (A) phase displays the “topological” degeneracy described above, I think that “topologically trivial” is not a good term to describe it. Another thing that I think it should be considered here. The abelian phase is expected to have 2x4 degeneracy, where the factor of 2 comes from time-reversal. On the other hand, the non-Abelian phase should display 2x3 degeneracy, as discussed by . Did you get any evidence of this?} By contrast, the (B) phase is a , the magnetic analogue of the fractional quantum Hall state. Topologically protected edge modes are predicted to occur in these states on periodic boundary conditions following the bulk-boundary correspondence . The probability density of one such edge mode is given in (a), where it is shown to be exponentially localised to the boundary of the system. The localization of these modes can be quantified by their inverse participation ratio (IPR), \[\begin{equation} \textup{IPR} = \int d^2r|\psi(\mathbf{r})|^4 \propto L^{-\tau}, \end{equation}\] where \(L\sim\sqrt{N}\) is the characteristic linear dimension of the amorphous lattices and \(\tau\) dimensional scaling exponent of IPR. Finally, the CSL density of states in open boundary conditions indicates the low-energy modes within the gap of Majorana bands in (b). { Could you plot the dimensional scaling exponent \(\tau\) in (a)?}The phase diagram of the amorphous model in \(\ref{fig:example_lattice}\)(c) displays a reduced parameter space for the non-Abelian phase when compared to the honeycomb model. Interestingly, similar inward deformations of the critical lines were found on the Kitaev honeycomb model subject to disorder by proliferating flux vortices or exchange disorder .
{} An Ising non-Abelian anyon is formed by Majorana zero-modes bound to a topological defect . Interactions between anyons are modeled by pairwise projectors whose strength absolute value decays exponentially with the separation between the particles, and whose sign oscillates in analogy to RKKY exchanges . Disorder can induce a finite density of anyons whose hybridization lead to a macroscopically degenerate state known as . One instance of this phase can be settled on the Kitaev CSL. In this case, the topological defects correspond to the \(W_p \neq +1\) fluxes, which naturally emerge from thermal fluctuations at nonzero temperature . We demonstrated that the amorphous CSL undergoes the same form of Anderson transition by studying its properties as a function of disorder. Unfortunately, we could not perform a complete study of its properties as a function of the temperature as it was not feasible to evaluate an ever-present boundary condition dependent factor for random networks. Instead, we evaluated the fermionic density of states (DOS) and the IPR as a function of the vortex density \(\rho\) as a proxy for temperature. This approximation is exact in the limits \(T = 0\) (corresponding to \(\rho = 0\)) and \(T \to \infty\) (corresponding to \(\rho = 0.5\)). At intermediate temperatures the method neglects to include the influence of defect-defect correlations. However, such an approximation is enough to show the onset of low-energy excitations for \(\rho \sim 10^{-2}-10^{-1}\), as displayed on the top graphic of \(\ref{fig:DOS_Oscillations}\)(a). We characterized these gapless excitations using the dimensional scaling exponential \(\tau\) of the IPR on the bottom graphic of the same figure. At small \(\rho\), the states populating the gap possess \(\tau\approx0\), indicating that they are localised states pinned to the defects, and the system remains insulating. At large \(\rho\), the in-gap states merge with the bulk band and become extensive, closing the gap, and the system transitions to a metallic phase. { Maybe being a bit more quantitative about \(\tau\) can enrich the discussion by allowing us to discuss a bit about the multifractality of these low-energy states} The thermal metal DOS displays a logarithmic divergence at zero energy and characteristic oscillations at small energies. . These features were indeed observed by the averaged density of states in the \(\rho = 0.5\) case shown in (b) for amorphous lattice. We emphasize that the CSL studied here emerges without an applied magnetic field as opposed to the CSL on the honeycomb lattice studied in Ref. { I have the impression that (b) on the top is very similar to Fig. 3 of . Maybe a more instructive figure would be the DOS of the amorphous toric code at the infinite temperature limit. In this case, the lack of non-Abelian anyons would be reflected by a gap on the DOS, which would contrast nicely to the thermal metal phase}% This high temperature phase of the amorphous model is known as a thermal metal. The signature of the thermal metal phase is characteristic oscillations in the low energy density of states, as seen in~(b).
{} 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 real-space generalisation of the Chern number. The presence of odd-sided plaquettes on these lattices let to a spontaneous breaking of time reversal symmetry, leading to 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. The next step is to search for an experimental realisation in amorphous Kitaev materials, which can be created from crystalline ones using several methods . Following the evidence for an induced chiral spin liquid phase in crystalline Kitaev materials , it would be interesting to investigate if a similar state is produced on its amorphous counterpart. Besides the usual half-quantized signature on thermal Hall effect , such a CSL could be also characterized using local probes such as spin-polarized scanning-tunneling microscopy . The same probes would also be useful to manipulate non-Abelian anyons , thereby implementing elementary operations for topological quantum computation. Finally, the thermal metal phase can be diagnosed using bulk heat transport measurements .This work can be generalized in several ways. Introduction of symmetry allowed perturbations on the model . Generalizations to higher-spin models in random networks with different coordination numbers
% In the present work, we have avoided the need for a rigorous Monte Carlo study of the thermal phase transition. As a consequence, the thermodynamic nature of the transition between the chiral QSL and thermal metal states has not been elucidated. { insert some guff about the Imri-Ma argument}.
{ Probably one way to make this theory experimentally relevant is to do experiments on amorphous phases of Kitaev materials. These phases can be obtained by liquifying the material and cooling it fast. Apparently, most of crystalline magnets can be transformed into amorphous ones through this process. } %Metal-organic frameworks (MOFs) are a promising candidate for realising Kitaev physics in an amorphous system. Yamada et al. propose a realisation of the Kitaev honeycomb model in a crystalline Ru-oxalate MOF~, and Misumi et al.~have demonstrated potential signatures of a resonating valence bond quantum spin liquid state in MOFs with Kagome geometry~. Amorphous MOFs can be generated by introducing disorder into crystalline MOFs through mechanical processes~, suggesting a natural route to realising amorphous Kitaev physics. Assuming it is possible to realise Kitaev physics in a crystalline MOF, it is unclear what superexchange couplings would be retained when disorder is introduced to the lattice. Because it is unlikely one would cleanly reproduce the exact model described in the present work, future work should examine how robust the CSL ground state of the amorphous Kitaev model is to additional disorder in the Hamiltonian, for example random recoloring of the bonds, additional bond forming and breaking, and disorder in coupling strengths.
% Produces the bibliography via BibTeX. %
A random pointset is used to partition space into polyhedral volumes enclosing the region closest to each point in the set. In two dimensions, the vertices and edges of these polygons form a tri-coordinate lattice.
%% The Kitaev honeycomb lattice model (HLM) is composed of spin-()
particles interacting anisotropically along the edges of a lattice: %
\begin{equation} % =
-{(i,j)}J^{{ij}}i^{{ij}}j^{{ij}} +
_{(i,j,k)}_i^{x}j^{y}k^{z} % \end{equation} % where the two spin
term runs over pairs of nearest neighbours and the three spin term runs
over consecutive triplets around a plaquette. The Pauli matrices \(\sigma^\alpha\) in each term are chosen
according to the type, or
coloring', of the bond $i\to j$, $\alpha_{ij}\in\{x,y,z\}$. The bond coloring is chosen such that exactly one bond of each type is connected to each vertex. % This Hamiltonian is exactly solvable by introducing a Majorana representation \(\widetilde{\sigma}_i^{\alpha} = i b^{\alpha}_i c_i\) which the partitions the Hilbert space into a classical $\mathrm{Z}_2$ gauge degree of freedom, \(u_{jk} = ib_j^{\alpha_{jk}}b_k^{\alpha_{jk}}\), on the bonds, and Majorana fermions, $c_j$, living on the vertices. It also doubles the size of the Fock space, necessitating calculating a projector \(P\) from the Majorana Fock space \(\mathcal{\widetilde{M}}\) onto the physical subspace \(\mathcal{M}\)~\cite{pedrocchiPhysicalSolutionsKitaev2011}. We refer to a choice of gauge configuration, $\{u_{jk}\}$, as theflux
sector’. The problem then reduces to solving a free-fermion Hamiltonian
within each flux sector (u) % \begin{equation} % ^u =
{j,k}A{jk}c_jc_k % \end{equation} % where \(A_{jk}=2J^\alpha_{jk}u_{jk}\) for \((j, k)\) nearest-neighbours, \(A_{jk}=2\kappa\sum_l u_{jl}u_{kl}\) for
\((j,k)\) second-nearest-neighbours,
and \(A_{jk}=0\) otherwise. Finally the
Majorana modes can be found with a transformation (Q) % \begin{equation}
% (b^{’}_1, b^{’’}_1, … ;b^{’}_N, b^{’‘}N) = (c_1, c_2, …
;c{2N}) Q % \end{equation} % from which we create the fermionic
operators (a_i = (b^{’}_i + ib^{’’}_i)), bringing (H) to the form %
\begin{equation} % ^u = _m _m (n_m - ) % \end{equation} % with ground
state energy (E_0 = -_m m). The projector has the effect of removing
many body states with either even or odd parity (= i (1 - 2n_i)), an
effect which typically leads to a correction of order (). The gauge
symmetries of \(\{u_{jk}\}\) can be
removed by defining plaquette operators (P_i = {(i,j) P_i}
u{ij}) that wind the plaquettes (faces) of the lattice.
% The ground state flux sector of the HLM in the isotropic phase (\(J^x = J^y = J^z\)) at zero field (\(\kappa=0\)) possesses a gapless fermionic spectrum. A non-zero field (\(\kappa\neq0\)) opens a gap, and the resulting fermionic insulator is known to host non-Abelian anyonic excitations and possess a non-zero Chern number~. This non-abelian phase has been shown to undergo a finite-temperature phase transition to a so-called `thermal metal’ phase, which exhibits multifractility~.
For a lattice with (B) bonds, (V) vertices, (P) plaquettes and Euler characteristic () (0 for the torus) the Euler equation states that (B = P + V + ). This corresponds to the (2^{B}) gauge configurations being composed of (2^{P - 1}) physically distinct vortex states each of which is composed of (2^{V - 1}) gauge equivalent states that correspond to flipping three (u_{ij}) around a vertex, along with (2 - ) non-contractible loop operators. The term (2 - ) is perhaps more easily understood by relating () to the genus of the surface (g), i.e the number of holes with (= 2 - 2g) showing that there are two non-contractible loops for each hole in the surface.
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.
{ Add brief mention of fermions and many body ground state} Closely following the derivation of~ we can extend to the amorphous case relatively simply. The main quantity needed is the product of the local projectors (D_i) [_i^{2N} D_i = _i^{2N} b^x_i b^y_i b^z_i c_i ] for a lattice with (2N) vertices and (3N) edges. The operators can be ordered by bond type without utilising any property of the lattice. [_i^{2N} D_i = _i^{2N} b^x_i _i^{2N} b^y_i _i^{2N} b^z_i _i^{2N} c_i] 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 = ) that arise from reordering the b operators such that pairs of vertices linked by the corresponding bonds are adjacent. [i^{2N} b^i = p{(i,j)}b^_i b^_j] This is simple the parity of the permutation from one ordering to the other and can be computed easily 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) [P^0 = 1 + p_x;p_y;p_z (Q^u) ; ; {{i,j}} -iu{ij}]
((Q^u)) is the determinant of the matrix that takes ((c_1, c_2… c_{2N}) Q = (b_1, b_2… b_{2N})). This along with (u_{ij}) depend on the lattice and the particular vortex sector.
( = ^{N} (1 - 2_i)) is the parity of the particular many body state determined by fermionic occupation numbers (n_i). The Hamiltonian is (H = _i (n_i - 1/2)) in this basis and this tells use that the ground state is either an empty system with all (n_i = 0) or a state with a single fermion in the lowest level.
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 \(\phi_{\textup{g.s.}} = -(\pm i)^{n_{\textup{sides}}}\). This was done by generating a large number (\(\sim\) 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^\alpha = 1\)), as well as in the toric code phase (\(J^x = J^y = 0.25, J^z = 1\)). The argument has one complication: for a graph with \(n_p\) plaquettes, there are \(2^{n_p - 1}\) distinct flux sectors to search over, with an added factor of 4 when the global fluxes \(\Phi_x\) and \(\Phi_y\) 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 sizes – we were able to check all flux sectors for systems with \(n_p \leq 16\) in a reasonable amount of time. However, at such small system size we find that finite size effects are substantial enough to destroy our results. In order to overcome these effects we tile the system and use Bloch’s theorem (a trick that we shall refer to as 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.
{} For a given lattice and flux sector, defined by \(\{ u_{jk}\}\), the fermionic ground state energy is calculated by taking the sum of the negative eigenvalues of the matrix \[\begin{align} M_{jk} = \frac{i}{2} J^{\alpha} u_{jk}. \end{align}\] The set of bond variables \(u_{jk}\), which we are free to choose, determine the \(\mathbb Z_2\) gauge field. However only the fluxes, defined for each plaquette according to eqn.~\(\ref{eqn:flux_definition}\), 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^{n_p - 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.
{} 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 . In order to circumvent this problem, we treat the 16-plaquette 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 $\bf q\(. We then check if the lowest energy flux sector aligns with our ansatz (eqn.~\ref{eqn:gnd_flux}) and whether this flux sector is gapped. \par 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.~\(\ref{fig:energy_gaps_example}\). 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.{} 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, 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.
% Thus, we shall begin with a discussion of how finite size affects the eigenvalues of the Majorana Hamiltonian, followed by our solution to this problem. Evidence for the ground state solution was collected by searching over all possible flux sectors for the lowest energy states. This is repeated for various values of \(J\) over a large number of randomly generated lattices.% For a given lattice and flux sector, defined by \(\{ u_{jk}\}\), the fermionic ground state energy is found by taking the sum of the negative eigenvalues of the matrix % \begin{align} % M_{jk} = J^{} u_{jk}. % \end{align} % A gauge transformation involves flipping the value of \(u_{jk}\) for the three bonds connected to the point at \(j\). Under a gauge transformation, the matrix \(M\) transforms according to \(M \rightarrow D_j M D_j\), where the matrix \(D_j\) is a diagonal matrix with \(-1\) on the \(j\)’th entry, and \(+1\) on all others. This represents a unitary transformation, so the spectrum of \(M\) is invariant under gauge transformations. As demonstrated in , the spectrum is determined entirely by the flux through all circuits in the system, which we define analogously to \(\ref{eqn:flux_definition}\). In this case we include not only plaquettes, but circuits that encircle several plaquettes. In periodic boundaries we must also consider
% In the language of graph theory, this matrix may be interpreted as representing a weighted, directed digraph, with weights determined by the individual entries of \(M\). The Harary-Sachs theorem states that the characteristic polynomial of such a matrix may be written in terms of the weights of the cycles of the graph, defined as the product of the elements of \(M\) around some closed path \(\mathcal C\) on the lattice, % \begin{align} % w_{} = {} M{jk}. % \end{align} % These weights are similar to the fluxes defined in the bulk text, with two important differences. Firstly, the cyclic weights include the factor of \(J^\alpha\) in the product. Secondly, unlike tthe fluxes, which are defined for individual plaquettes, the weights are calculated for every closed path on the lattice. The takeaway is that the characteristic polynomial, and thus all eigenvalues, is determined only by the values of these weights. Any change to the set of \(u_{jk}\) that does not affect the weight of any cycles will have no effect on the energies of the system. For example a gauge transformation, where \(u_{jk}\) is flipped on the three edges connected to a chosen site, cannot affect the energies, as every cycle passing through the chosen site must contain two of the flipped edges.\end{document}