Efficient computation of topological order

Abstract

We analyze the computational aspects of detecting topological order in a quantum many-body system. We contrast the widely used topological entanglement entropy with a recently introduced operational definition for topological order based on error correction properties, finding exponential scaling with the system size for the former and polynomial scaling for the latter. We exemplify our approach for a variant of the paradigmatic toric code model with mobile particles, finding that the error correction method allows to treat substantially larger system sizes. In particular, the phase diagram of the model can be successfully computed using error correction, while the topological entanglement entropy is too severely limited by finite size effects to obtain conclusive results. While we mainly focus on one-dimensional systems whose ground states can be expressed in terms of matrix product states, our strategy can be readily generalized to higher dimensions and systems out of equilibrium, even allowing for an efficient detection of topological order in current quantum simulation experiments.

Topologically ordered states of matter transcend Landau’s symmetry-breaking paradigm and cannot be described by local order parameters Wen (2019). While such states are both important from a fundamental point of view and for applications in quantum computing, computing the topological properties of a many-body state is a notoriously difficult task. Here, we show that an approach based on the error correction properties of a state allows for an efficient computation of the phase boundaries of topologically ordered states, requiring only polynomial computational time in the size of the system.

While there are many conceptual frameworks to identify topological order Thouless et al. (1982); Kohmoto (1985); den Nijs and Rommelse (1989); Haegeman et al. (2012); Pollmann and Turner (2012); Elben et al. (2020); Chen et al. (2010); Kitaev and Preskill (2006); Levin and Wen (2006); Jiang et al. (2012); Zhang et al. (2012); Zhu et al. (2013); Nussinov and Ortiz (2009); Qiu and Wang (2020), their use in actual calculations or even experiments remains limited. The most pratically relevant quantity is the topological entanglement entropy Kitaev and Preskill (2006); Levin and Wen (2006); Jiang et al. (2012); Zhang et al. (2012), but being a highly nonlinear functional of a reduced density matrix, it can only be evaluated for small cluster sizes Satzinger et al. (2021). This underlines that the computational decision whether a state is topologically ordered or not remains an outstanding problem.

In this Letter, we show that these challenges can be overcome when computing topological properties based on an operational definition of topological order related to the error correction properties of a many-body state Jamadagni and Weimer (2022a). We exemplify this approach for a variant of the toric code, where the consituent spins can move around on a lattice. Using matrix product state (MPS) simulations, we show that the ground state undergoes a series of two distinct phase transitions. We emphasize that this behavior cannot reliably captured using the topological entanglement entropy, which shows a continuous dependence on the strength of the hopping of the particles rather than converging to a discrete value.

Topological Order.— Quantum phase transitions involving local order parameters can be successfully classified within Landau’s symmetry breaking paradigm Sachdev (2011). However, in the last decades, many examples have been discovered that transcend this notion. Most prominently, gapped quantum phases exhibiting topological order Wen (2017) exhibit a wide range of unrivaled properties such as long-range entanglement, robustness against local perturbations and the ability to perform fault-tolerant quantum computing using such many-body states.

One major challenge surrounding topological order is that its inherent nonlocality makes the classification of topologically ordered states of matter very difficult. On the one hand, this refers to conceptual challenges, i.e., identifying which nonlocal quantities should be considered, especially for systems approaching the thermodynamic limit of infinite system sizes. On the other hand lie practical challenges concerning the formalization of a decision problem allowing to compute whether a given many-body state is topologically ordered or not. For the latter, it is of special importance that the decision problem can be solved efficiently, i.e., computational ressources scale at most polynomially in the system size, as otherwise only very small systems can be probed.

Topological entanglement entropy and its limitations.— While there are many different concepts to describe topological order Thouless et al. (1982); Kohmoto (1985); den Nijs and Rommelse (1989); Haegeman et al. (2012); Pollmann and Turner (2012); Elben et al. (2020); Zhang et al. (2012); Zhu et al. (2013), resorting to the topological entanglement entropy (TEE) Kitaev (2006); Levin and Wen (2006) is one of the most widely used approaches. In its essence, the TEE $\gamma$ captures the subleading correction to the entanglement entropy when looking at the reduced density matrix of a system cut along a path of length $L$, i.e.,

$$S=\alpha L-\gamma+O(1/L).$$ (1)

The TEE is connected to the total quantum dimension of a topologically ordered state Kitaev (2006) and hence directly captures topological order. For practical calculations, the TEE can be computed from a cylindrical geometry Jiang et al. (2012), see Fig. 1a. In many cases, it is possible to derive the correct value of the two-dimensional (2D) case in the limit where the $Z$ extension of the cylinder is very small Zou and Haah (2016). One important requirement for this to work is that the model does not include boundary terms that immediately destroy topological order Jamadagni et al. (2018), i.e., the quasi-1D model is a variant of symmetry-protected topological order Chen et al. (2011) corresponding to the intrinsically topologically ordered 2D model.

Figure 1: Sketch of the topology. A cylindric system (a) is mapped onto a quasi-1D ladder geometry in the limit of small $Z$ (b). The cut of length $L$ for computing the the topological entanglement entropy is shown as a dotted line. For the mobile toric code, the system consists of two different sublattices (circles and squares) on which $A_{v}$ and $B_{p}$ operators act to implement the toric code interactions. Mobility is introduced by allowing the particles to hop along the dashed lines.

Unfortunately, actually obtaining the TEE is often prohibitively challenging. In general, computing the entanglement entropy of a many-body system scales exponentially in the system size. This is true even in a quasi-1D situation, where the ground state can be approximately represented by a matrix product state (MPS) of the form Perez-Garcia et al. (2007)

$$\mathopen{|}\psi\mathclose{\rangle}=\sum_{i_{1},....i_{N}}^{d}\text{Tr}[A_{i_{1}}^{[1]}A_{i_{2}}^{[2]}......A_{i_{N}}^{[N]}]\mathopen{|}i_{1},i_{2},....,i_{N}\mathclose{\rangle},$$ (2)

where $\mathopen{|}\psi\mathclose{\rangle}$ is a state of the Hilbert space of $N$ particles of local Hilbert space dimension $d$. The matrices $A_{i}$ have a maximum size of $\chi\times\chi$, with $\chi$ being the maximum bond dimension, which is a measure of complexity of the many-body state. While for many 1D problems, $\chi$ grows at most polynomially with system size, allowing to compute expectation values in polynomial time, this is typically not the case for computing the TEE Hamma et al. (2008). Generically, the cost of computing the entanglement entropy scales exponentially in the number of bonds broken by the cut, resulting in $O(d^{3N/2})$ steps for the cut required for the TEE on a square ladder, see Fig. 1b. In practice, this means that it is only possible to use the TEE for computing topological order for spin $1/2$ systems having $d=2$ and provided that the TEE converges sufficiently fast in the system size.

Operational definition of topological order.— To obtain a computationally viable approach to topological order, we turn to a recently introduced operational definition Jamadagni and Weimer (2022a). In a nutshell, a state is topologically ordered if and only if it can be corrected to a topologically ordered reference state using an error correction circuit of finite depth. Here, the circuit depth refers to the number of steps that are required for the circuit to complete in a fully parallelized implementation. In our approach, the reference state takes a similar role as the order parameter in ordered phases involving spontaneous symmetry breaking. We note that this operational definition is related to efforts to enhance the topological order in quantum states by means of error correction Cong et al. (2024). The error correction algorithm can be formulated as a sequence of (i) syndrome measurements to detect the errors with respect to the reference state and (ii) quantum operations to correct the errors. In sharp contrast to the TEE, the operational definition does not require the computation of a nonlinear functional and hence can be implmented in polynomial time within an MPS framework. For practical purposes, the most convenient implementation of the operational definition corresponds to a Monte-Carlo algorithm Jamadagni and Weimer (2022b): Starting from the MPS state, $O(N)$ projective measurements are being performed on the state, randomly selecting a measurement outcome according to the respective expectation values. After each projection, the MPS has to be reorthogonalized using $O(N)$ singular value decompositions (SVDs) Orus (2014). Each SVD carries a computational cost of $O(d\chi^{3})$, leading to a total computational cost of $O(dMN^{2}\chi^{3})$, with $M$ being the number of Monte-Carlo samples. This polynomial scaling demonstrates a drastic improvement over the exponential complexity of computing the TEE.

The mobile toric code.— Let us now turn to a concrete model to apply our error correction method and compute its phase diagram. Here, we will be interested in an extension to the toric code, where the constituents particles can move around on a square lattice, see Fig. 1b. As an experimental realization, we envision an implementation of the toric code using ultracold atoms in an optical lattice undergoing Rydberg interactions Weimer et al. (2010), where residual movement of the atoms is a natural imperfection stemming from the finiteness of the optical lattice potentials Bloch et al. (2008). The toric code is described by the spin-$1/2$ Hamiltonian

$$H_{\text{toric}}=-\sum_{\text{vertex}}\underset{A_{v}}{\underbrace{\sigma^{i}_{x}\sigma^{j}_{x}\sigma^{k}_{x}\sigma^{l}_{x}}}-\sum_{plaquette}\underset{B_{p}}{\underbrace{\sigma^{m}_{z}\sigma^{n}_{z}\sigma^{o}_{z}\sigma^{p}_{z}}}.$$ (3)

Since all $A_{v}$ and $B_{p}$ operators mutually commute, the ground states satisfy the constraints $A_{v}\mathopen{|}\psi\mathclose{\rangle}=\mathopen{|}\psi\mathclose{\rangle}$ and $B_{p}\mathopen{|}\psi\mathclose{\rangle}$ for all $v$ and $p$ Kitaev (2003). Violations of the constraints correspond to errors that can be detected by measuring the $A_{v}$ and $B_{p}$ operators.

In order to account for the motion of particles, we extend the local Hilbert space by a third possible state $\mathopen{|}0\mathclose{\rangle}$ indicating the absence of a particle. In the following, we assume a hard-core constraint for the particles, i.e., there is at most one atom per lattice site, which is well justified in sufficiently deep optical lattices Bloch et al. (2008). Then, the local basis for each lattice site is given by the set $\{\mathopen{|}\downarrow\mathclose{\rangle},\mathopen{|}0\mathclose{\rangle},\mathopen{|}\uparrow\mathclose{\rangle}\}$, where $\mathopen{|}\uparrow\mathclose{\rangle}$ and $\mathopen{|}\downarrow\mathclose{\rangle}$ refer to an atom in the spin up or down state, respectively. Using this notation, we can define bosonic creation and annhilation operators according to

$$a_{\uparrow}^{\dagger}=\mathopen{|}\uparrow\mathclose{\rangle}\mathopen{\langle}0\mathclose{|},\;\;\;\;\;a_{\downarrow}^{\dagger}=\mathopen{|}\downarrow\mathclose{\rangle}\mathopen{\langle}0\mathclose{|}.$$ (4)

The Pauli spin operators can then be understood in the usual way, e.g., $\sigma_{x}=\mathopen{|}\uparrow\mathclose{\rangle}\mathopen{\langle}\downarrow\mathclose{|}+\text{H.c.}$, with the only difference being that they are now acting on a three-dimensional local Hilbert space. As a consequence, $A_{v}$ and $B_{p}$ operators now also have eigenstates with eigenvalue of zero; this occurs whenever one of the sites taking part in the vertex $v$ or the plaquette $p$ is in the vacuum state $\mathopen{|}0\mathclose{\rangle}$. Then, the total Hamiltonian has the form

$$H=H_{\text{toric}}-t\sum_{\langle ij\rangle\;s}a_{i,s}a^{{\dagger}}_{j,s}+\text{H.c.},$$ (5)

where the sum over $s$ runs over both spin states. Since the square lattice has twice as many lattice sites as the checkerboard lattice of the original toric code, the toric code part has to be doubled to act on both sublattices, see the top and bottom part of Fig. 1b.

As discussed previously, we focus on a quasi-1D realization of the model. Similar to the toric code in a magnetic field Jamadagni and Weimer (2022a), the minimal instance of the mobile toric code is given by a ladder geometry, where the four-body interactions of the original model are replaced by three-body interactions, see Fig. 1b. Before performing the error correction analysis, let us first remark on some general properties of the phase diagram. For zero hopping $t$, the system corresponds exactly to the quasi-1D limit of the toric code and hence exhibits intrinsic topological order when extended back into 2D. Additionally, the doubling of the lattice sites creates an additional $\mathbb{Z}_{2}$ symmetry corresponding to which of the two sublattices are being populated with atoms, with the other half being empty for $t=0$. For very large values of $t$, the toric code part irrelevant and the ground state corresponds to a topologically trivial Luttinger liquid phase. At intermediate $t$, both the topological order and the density-wave order of the broken $\mathbb{Z}_{2}$ symmetry will eventually disappear. However, it is possible that both types of order disappear at the same time or there is an intermediate phase exhibiting density-wave order but no topological order. Hence, we will employ our computationally efficient strategy to determine the entire phase diagram.

Error correction circuits.— For the mobile toric code, the error correction algorithm consists of two parts. First, the positions of the particles have to be corrected to a perfect crystal before the conventional error correction to the toric code can be applied. Immediately applying the error correction for the toric code would result in some $A_{v}$ and $B_{p}$ operators being measured in the zero eigenvalue, which cannot be corrected. Importantly, the depth of the error correction circuit for correcting the positions of the particle is a measure of density-wave ordering Jamadagni and Weimer (2022b). If this part of the error correction circuit has finite depth, the system exibits density-wave ordering; if the circuit depth diverges with the size of the system and the system enters a Luttinger liquid phase. Hence, density-wave ordering is a necessary requirement for the system being able to exhibit topological order.

Figure 2: Error correction circuit for density-wave order. (a) The local particle configuration is measured and the error syndrome for density-wave order is computed along the dashed lines with error variables being located on the dual lattice (squares). Two neighboring particles signal the presence of a particle error (+1), while two empty sites correspond to a hole error (-1). (b) Errors are removed by swapping particle positions until a particle error is fused with a hole error. All measurement and correction operations are carried out without disurbing the spin state.

Let us first describing the error correction algorithm for density-wave order. As the first step, we identify the position of the errors by measuring the positions of the particle without disturbing the spin sector. We perform these measurements on an MPS by applying our Monte-Carlo algorithm. Errors are identified according to the measurement results by transversing the system in the pattern shown in Fig. 2, with errors being located on the dual lattice between the sites of the physical lattice. If there is exactly one particle on two neighboring site, there is no error (0), while two particles correspond to a particle error (+1) and two empty sites correspond to a hole error (-1). Error correction is now being performed by decorating each error with a walker that explores the surroundings of the error, with each movement of the walker by one site corresponds to one timestep of the error correction circuit Jamadagni and Weimer (2022a). Once a walker from a particle error encounters a hole error or vice versa, the errors can be removed by swapping the particle positions in between such that an error-free configuration is obtained. Again, this operation is carried out without affecting the spin state of the particles. Finally, the total number of timesteps required to remove all errors corresponds to the depth of the density-wave circuit.

Figure 3: Depths of error corrections circuits for the mobile toric code. (a) Total circuit depth $d_{tot}$ as a function of the hopping strength $t$ for different system sizes. For large $t$, the circuit depth diverges in the system size indicating the absence of topological order. (b) Finitize size scaling of the peak of $\partial d_{tot}/\partial t$, leading to a critical value of $t_{1}=0.37(1)$ in the thermodynamic limit for the breakdown of topological order. (c) Correcting only the density-wave order leads to a critical value of $t_{2}=0.50(1)$ for the disappearance of density-wave order. Number of trajectories range from $M=5000$ for $N=14$ to $M=2000$ for $N=30$.

Once the positions have been corrected to the perfect crystal, we perform the usual correction of the toric code errors by pairwise fusion of the $A_{v}$ and $B_{p}$ violations Jamadagni and Weimer (2022a). The circuit depth of this part is added to the circuit depth of the density-wave ordering $d_{dw}$ to obtain the total circuit depth $d_{tot}$ for detecting topological order. Figure 3 shows the total circuit depth for systems up to $N=30$ particles. MPS simulations were carried out using the TenPy library Hauschild and Pollmann (2018) using a maximum bond dimension of $\chi=3000$. We note that this bond dimension is required to achieve convergence of the circuit depth in the Luttinger liquid phase. We locate the phase transition points by looking at the peaks of the derivatives $\partial d_{tot}/\partial t$ and $\partial d_{dw}/\partial t$, respectively. The derivatives are obtained by fitting a fourth-order polynomial to the simulation data. Using a finite-size scaling analysis, we find that topological order disappears at a critical hopping strength of $t_{1}=0.37(1)$ and density-wave order disappears at $t_{2}=0.50(1)$. Hence, our numerical simulations provide strong evidence for the existence of an intermediate density-wave ordered phase that is topologically trivial.

Figure 4: Topological entanglement entropy $\gamma$ of the mobile toric code for system sizes up to $N=14$. For $t=0$, the TEE corresponds to the toric code value of $\gamma=\log 2$, but strong finite size effects immediately lead to a reduction in the TEE even though the system is still in its topologically ordered phase.

Comparison to TEE calculations. We now try to investigate the topological properties of the mobile toric code using the topological entanglement entropy. As a first obstacle, we notice that the exponential scaling limits the accessible system size to at most $N=14$, which is much smaller than the error correction computations when using comparable computational resources. For zero hopping, the TEE is the same as for the toric code on a cylinder, i.e., we have $\gamma=\log 2$. For large values of $t$, we again expect a transition to a topologically trivial phase having $\gamma=0$. However, finitize-size effects are expected to appear as a rounding of the step-like behavior of the TEE Morampudi et al. (2014); Jamadagni et al. (2018). Technically, the TEE is computed by calculating the entanglement entropy for large system sizes and extrapolating back to $L=0$ to obtain the subleading term according to Eq. 1. Figure 4 shows the behavior of the TEE as a function of the hopping strength $t$. As can be clearly seen, the results are completely inconclusive once the hopping is switched on, having little resemblance to the expected step function behavior. Even for very small hopping strengths, the TEE starts to deviate from the toric code value of $\gamma=\log 2$, indicating a loss of topological order much earlier than is actually the case. These findings unambigously demonstrate the advantages of the error correction approach to topological order.

Summary and outlook.— In summary, we have presented an efficient numerical strategy to compute topological order in a quantum many-body system. Our approach relies on an operational definition for topological order and is exponentially faster than computing the topological entanglement entropy. We have applied our computational framework to the example of an extension toric code where particles are mobile, finding that we can successfully obtain the phase diagram. Our approach naturally generalizes to higher-dimensional tensor network states such as projected entangled-pair states Orus (2014); Cirac et al. (2021), as well as to topological order in non-equilibrium situations Bardyn et al. (2013); Yarloo et al. (2018); Petiziol et al. (2024). Finally, we would like to point that the quantities required to compute topological order can be measured in state-of-the-art quantum simulator experiments Semeghini et al. (2021); Chen et al. (2023), while experimentally obtaining the topological entanglement entropy requires quantum state tomography and is again exponentially hard.