1 Introduction↩︎

Gibbs state preparation is a key computational technique used in physics, statistics, and many other scientific fields. Given a local Hamiltonian \(H\) and an inverse temperature \(\beta\), the Gibbs state \(\rho_{\beta H}\sim\exp(-\beta H)\) describes the thermal equilibrium properties of quantum systems at finite temperature, making them essential for studying the phase diagram, stability of topological quantum memory [5], [6] as well as the thermalization process [7], [8]. In addition to physics, Gibbs state preparation also has found various applications in optimization [9], [10] and Bayesian Inference [11], [12]. Various Gibbs state preparation algorithms (or Gibbs samplers) have been proposed, including approaches inspired by the Davies generator [13]–[16], the Metropolis-Hasting type method [17], [18], and ones based on Grover amplification [19] and the Quantum Singular Value Transform [20].

The key requirement for a good Gibbs sampler is fast mixing, that is, the algorithm prepares the Gibbs state in polynomial time. Gibbs samplers for classical Hamiltonians have been studied for decades and fast mixing algorithms have been successfully designed for various scenarios. In particular, Glauber dynamics yield fast mixing Gibbs samplers for 1D systems at any constant temperature [21]–[23] and for 2D systems at high temperature [24]. On the other hand, for 2D systems like the Ising model, Glauber dynamics-based samplers are known to suffer critical slow downs and are slow mixing at low temperature [25]–[27]. Instead, the Swendsen-Wang algorithm [28], [29] and approaches based on Barvinok’s method [30] were proved to achieve fast mixing for low temperature 2D systems.

Recent efforts on developing fast mixing Gibbs samplers for quantum Hamiltonians have largely focused on commuting local Hamiltonians (CLHs). CLHs are an important subclass of quantum Hamiltonian which exhibits non-trivial quantum phenomenon. Different from classical Hamiltonians whose eigenstates are computational basis, the eigenstates of a CLH instance can be highly entangled and cannot be prepared by any constant depth quantum circuit, as is true for the famous example of Kitaev’s Toric code [31]. Besides, it was shown that Gibbs sampling of CLHs at constant temperature remains classically hard [32], [33]. Nonetheless, several aforementioned fast mixing results for classical Hamiltonians have been successfully generalized to the CLH case. In particular, multiple results utilize the Davies generator [34], [35], which represents a quantum Markov chain (Lindbladian) for thermalization process in the weak coupling limit. This has yielded fast mixing Gibbs samplers for 1D CLH at any constant temperature  [2], [3] and 2D CLH at high temperature [1], [3]. The fast mixing proofs are obtained by generalizing classical techniques for analyzing the mixing time of transition matrices [21]–[24], [24] to analyzing the mixing time of the Davies generator. These generalizations are highly non-trivial and very technical since analyzing the spectrum of a quantum operator (in this case the Davies generator) is generally hard. In the low temperature regime, fast mixing Gibbs samplers for CLHs on two or higher dimensions are only known for the standard 2D Toric code [36], [37].

In this work, we introduce a new approach which does not use the Davies generator. Instead, we design new Gibbs samplers for various CLHs by giving a reduction from quantum Gibbs state preparation to classical Gibbs sampling. Combined with the existing fast mixing results for classical Hamiltonians, our Gibbs sampler is able to replicate the state-of-the-art performances mentioned above. Furthermore, our algorithm can prepare low temperature Gibbs states as long as there exists a fast mixing Gibbs sampler for the corresponding classical Hamiltonians. Our reductions are summarized in the following two theorems. More details and comparisons between previous results and the performance of our Gibbs sampler are contained in 1.1.

Roughly speaking, we say that there is a Gibbs sampling reduction from a quantum Hamiltonian \(H\) to a classical Hamiltonian \(H^{(c)}\) if, given the existence of an algorithm that performs Gibbs sampling for \(( H^{(c)},\beta)\) in time \(T\), there exists a quantum algorithm preparing the quantum Gibbs state for \((H,\beta)\) in time \(T\) plus a small overhead polynomial in the number of qubits.³ First we notice that the Structure Lemma, which is the key technique used in studying the complexity of CLHs [38]–[41], directly gives the desired reduction for 2-local CLHs.

For more physically motivated 4-local CLHs such as the Toric code, the Structure Lemma can no longer transform 4-local CLHs to classical Hamiltonians. Instead, by leveraging a symmetry in the eigenspaces, we demonstrate that an oblivious randomized correction approach yields the desired reduction for qubit CLHs in 2D, via generalizing a ground state preparation algorithm [40] to Gibbs state preparation. The locality of the resulting classical Hamiltonian depends on whether there are classical qubits in the CLHs. Roughly speaking a qubit is classical if by choosing proper basis of this qubit, all terms look like \(\vert 0\rangle\!\langle 0\vert\otimes ... + \vert 1\rangle\!\langle 1\vert\otimes ...\) on this qubit.

Note that in the above theorem the quantum Hamiltonian is on qubits while the classical Hamiltonian is on qudits. An example of a qubit 2D CLH without classical qubits is the defected Toric code, a generalization of the Toric code with arbitrary complex coefficients. We will give a technical overview based on the defected Toric code in 1.3.

Our reduction also gives a quantum analogy of Stockmeyer’s result [4] for the complexity of quantum approximate counting. In particular, a fundamental result of Stockmeyer [4] states that classical approximate counting (approximating the partition function of a classical Hamiltonian) is contained in the complexity class \(\textsf{BPP}^{\textsf{NP}}\). It is natural to conjecture that the quantum approximate counting is upper bounded by a complexity class like \(\textsf{BQP}^{\textsf{QMA}}\), but few results are known. By the connection between quantum approximate counting and the Gibbs state preparation [42]⁴, our reduction shows that for various CLHs, the corresponding quantum approximate counting problem is contained in \(\textsf{BQP}^{\textsf{CS}}\), where \(\textsf{CS}\) is an oracle which can perform arbitrary classical Gibbs sampling.

1.1 Comparison to previous work.↩︎

Recall that most of the previous work on Gibbs samplers for CLHs are based on simulating the Davies generator, which is a Lindbladian closely related to the thermalization process. In this section, we give a detailed comparison between previous results and our result, demonstrating that instead of using the Davies generator, our reduction gives a new Gibbs sampler for CLHs directly utilizing fast mixing Gibbs samplers for classical Hamiltonians. In particular, our reduction is able to replicate state-of-the-art results and also derive new results. Our results are summarized in 1. In this section we discuss related results, and a more thorough discussion on Gibbs state preparation can be found in [1].

1.1.0.1 Review of Markov chains, Lindbladians and mixing time

We briefly review some key concepts. Consider an \(n\)-qubit local Hamiltonian \(H\) and an inverse temperature \(\beta\). We assume \(\beta \in \mathcal{O}(1)\) unless further specified. We will first assume \(H\) is classical and introduce key concepts for classical Gibbs sampling. Then we will generalize to quantum Hamiltonians.

Suppose \(H\) is a classical Hamiltonian which diagonalizes in the computational basis, and the task is performing classical Gibbs sampling for \((H,\beta)\). The target is the Gibbs distribution \(\pi\) which samples computational basis states \(\left\vert x \right\rangle\) with probability proportional to \(\exp\left(-\beta \langle x|H|x\rangle\right)\). The goal of classical Gibbs sampling is to design a classical process which drives any distribution \(\nu\) to \(\pi\). The commonly used method is the classical Metropolis algorithm [43], which is a discrete-time Markov chain described by a transition matrix \(P\), such that \(\pi\) is the unique fixed point of \(P\), i.e. \(P\pi=\pi\). The mixing time \(t(\epsilon)\) is the time needed to get \(\epsilon\)-close to the invariant distribution \(\pi\) with respect to \(1\)-norm (the total variation distance), that is \[\begin{align} t(\epsilon):=\min\{ t: \|P^t\nu-\pi\|_1\leq \epsilon,\forall \text{ distribution } \nu \}. \end{align}\] In addition to this discrete-time Markov chain, one can also design a continuous-time Markov chain, described by a generator matrix \(G\) such that \(\pi\) is the unique invariant distribution of \(G\), i.e. \(G\pi=0\) or equivalently \(e^{Gt}\pi=\pi\), \(\forall t\). Similarly to above, the mixing time is defined to be \[\begin{align} t(\epsilon):=\min\{ t: \|e^{Gt}\nu-\pi\|_1\leq \epsilon,\forall \text{ distribution } \nu \}. \end{align}\]

• The Markov chain is poly-time mixing, or fast mixing, if the the spectral gap of \(P\) or \(G\) is \(\Omega(1/poly(n))\), which implies \(t(\epsilon)=poly(n)\times \log \frac{1}{\epsilon}\).

• The Markov chain is rapid mixing if it reaches the invariant distribution in a time which scales logarithmically with the system size, that is \(t(\epsilon)= poly(\log n) \times \log\frac{1}{\epsilon}\). Rapid mixing is typically proved by bounding the log-Sobolev constant [21] for the continuous-time chain.

When \(H\) is a quantum Hamiltonian, we wish to prepare a quantum Gibbs state, defined as \[\rho_{\beta H}:= \rho( H,\beta ):= \frac{1}{\text{tr}[\exp(-\beta H)]}\exp(-\beta H).\] The goal is to design a quantum process which drives any quantum state \(\sigma\) to \(\rho_{\beta H}\). One commonly used method is to design a Lindbladian \(\mathcal{L}\) such that \(\rho_{\beta H}\) is the unique fixed point of \(\mathcal{L}\), i.e. \(\mathcal{L}(\rho_{\beta H})=0\) or, equivalently, \(e^{\mathcal{L}t}(\rho_{\beta H})=\rho_{\beta H}, \forall t\). The Lindbladian is the quantum analogy of a continuous-time Markov chain. The mixing time \(t(\epsilon)\) is defined to be \[\begin{align} t(\epsilon):=\min\{ t: \|e^{\mathcal{L}t}(\sigma)-\rho_{\beta H}\|_1\leq \epsilon,\forall \sigma \}. \end{align}\] The notion of poly-time mixing and rapid mixing is defined similarly to the classical setting. One can prepare the quantum Gibbs state on a quantum computer by Lindbladian simulation techniques [13], [44]. For simplicity we will assume that \(\epsilon=1/poly(n)\) for the remainder of the section.

1.1.0.2 Replication of poly-log time mixing algorithm for 1D commuting Hamiltonians.

For 1D classical Hamiltonians, it is well-known that there is no computational phase transition [21]–[23]. As a result, for any constant temperature, Glauber dynamics is rapid mixing for all translation-invariant, 1D classical Hamiltonians with finite-range interactions.

A large body of previous work has focused on generalizing the rapid mixing results from classical Hamiltonians to quantum CLHs. In particular, for 1D CLHs, [3] proved that the Davies generator has a constant spectral gap and thus is fast (poly-time) mixing. Then Bardet et.al. [2] strengthened the result to obtain rapid mixing. Specifically, they proved that for any constant temperature, the Davis generator \(\mathcal{L}\) for any finite-range, translation-invariant 1D CLHs is rapid mixing. This is proved by generalizing the classical technique of bounding the log-Sobolev constant  [21]–[23], [45] to CLHs. Note that this generalization is highly non-trivial since \(\mathcal{L}\) is a quantum operator and analyzing its spectral gap and log-Sobolev constant are complicated. Combining this bound with quantum simulations of the Lindbladian, Bardet et.al. gives a quantum Gibbs sampler with runtime \[poly(\log n) \times f_1,\] where \(f_1\) is the overhead brought by simulating the Lindbladian evolution.

In contrast to [2], which obtained a rapidly mixing Gibbs sampler by developing sophisticated techniques to bound the log-Sobolev constant of the Davies generator, our reduction gives a Gibbs sampler of similar performance by directly using classical results:

Here \(f_2\) is the overhead incurred by simulating the classical Markov chain. \(\mathcal{O}(n)\) is the time needed to prepare a constant depth quantum circuit arising from the Structure Lemma which implements the quantum-to-classical reduction.

1.1.0.3 High and low temperature Gibbs state for 2-local CLHs on 2D.

Recall that \(\beta\) is the inverse temperature, thus low temperature corresponds to large \(\beta\). We begin with a literature review for classical Hamiltonians. Unlike 1D classical Hamiltonians where the Glauber dynamics is rapid mixing for any constant temperature, 2-local 2D classical Hamiltonians exhibit a constant-temperature computational phase transition. For example, the ferromagnetic 2D Ising model has a constant critical inverse temperature \(\beta_c\), such that

• For \(\beta< \beta_c\) the Glauber Dynamics is poly-time mixing [24], [24].

• For \(\beta\geq \beta_c\), the Glauber dynamics meets a critical slow down where the spectral gap of the Glauber Dynamics is smaller than \(\exp(-\alpha\sqrt{n})\) for \(\alpha>0\) [25]–[27].

Similar results also hold for the Potts model [46], [47]. To understand this phase transition intuitively, note that the Glauber dynamics is a Markov chain with local update rules. Intuitively an algorithm using local updates is good at solving a “local” problem. In the high temperature region, most spins will interact effectively weakly thus the Gibbs state has little entanglement [48] and a local optimization suffices. However, in the low temperature region, there are strong correlations in large regions and thus the Gibbs state is highly non-local.⁵ Thus, to prepare low temperature Gibbs state , one needs to carefully design Markov chains with non-local update rules, such as the cluster updates in the Swendsen-Wang algorithm [28]; or uses other methods such as the Barvinok’s method [30].

In the quantum case, to the best knowledge of the authors, all previous work on Gibbs state preparation for 2D CLHs has focused on the high temperature region. In particular,

• [3] showed that there is a constant \(\beta_1\) such that for \(\beta\leq \beta_1\) the Davies generator is poly-time mixing.

• [1] showed that there is a constant \(\beta_2\) such that for \(\beta\leq \beta_2\)⁶, the Schmidt generator defined in [1] is rapid mixing.

Both the Davies generator and the Schmidt generator for 2D CLHs with respect to local jump operators are local Lindbladians. An adaption of the classical proofs [25]–[27] will show that they are slow mixing for 2D systems at low temperature [53]. That is, there exists a constant \(\beta_3\) such that for \(\beta\geq \beta_3\), the spectral gap of any \(O(\log n)\)-local Lindbladian (not necessarily the Davies generator) which fixes the Gibbs state of 2D Ising model at inverse temperature \(\beta\) has an exponentially-small spectral gap.

Our Gibbs sampler improves on the existing results in two main ways. First, in the high temperature region our reduction again gives a way to directly utilize classical results [50] and obtain a Gibbs sampler of similar performance as the best prior work [1] (i.e., rapid mixing), without involving heavy proofs for analyzing the log-Sobolev constant of the Schmidt generator like [1].

Our second contribution is in the low temperature regime. Unlike [1], [3] which only work for the high temperature region, our reduction allows us to prepare low-temperature Gibbs states by utilizing classical techniques such as the Swendsen-Wang algorithm. To the best knowledge of the authors, prior work has not addressed low-temperature Gibbs samplers for 2-local CLHs. As an example, with our reduction we can obtain the following result,

A key feature of our reduction is that it is agnostic to the underlying classical Gibbs sampler. At the critical temperature, when applied to qudit 2-local 2D CLHs with large constant qudit dimension \(d\), the Swendsen-Wang algorithm will also mix slowly and have a spectral gap that is exponentially small in the side length of the lattice [54]. However, our reduction allows us to substitute in other samplers, such as the Gibbs sampler for the low temperature Potts model based on Barvinok’s method, which remains poly-time for large \(d\) [30].

1.1.0.4 4-local, 2D commuting Hamiltonian.

The best prior work is due to [3], who proved that for high enough temperature, the Davies generator is poly-time mixing for 4-local 2D CLHs. Unlike their result, the mixing time of our algorithm is dependent on the classical Hamiltonian produced by the reduction and thus our results are not directly comparable.

For the standard Toric code, a concurrent work [37] showed that for any inverse temperature \(\beta<+\infty\) (not necessarily constant), Lindbladian dynamics with nonlocal jump operators prepares the Gibbs state efficiently (for very low temperature, the mixing time is approximately \(\mathcal{O}(n^3)\)). Our Gibbs sampler is based on different techniques and gives a \(\mathcal{O}(n^ 2poly(\log n))\)-time Gibbs state preparation algorithm for the general defected Toric code at any non-zero temperature, where the defected Toric code is the Toric code with arbitrary coefficients. Our Gibbs sampler is based on generalizing the standard ground state preparation algorithm for the Toric code (which measures all stabilizers) to the task of Gibbs state preparation via an oblivious randomized correction technique. We will give a technical overview based on the example of defected Toric code in Section 1.2 and 1.3. We remark that since our Gibbs sampler is not based on Lindbladian, our results do not offer additional insights into Lindbladian dynamics, unlike in [37]. Another related work [55] uses classical Monte Carlo methods to simulate the Gibbs states of t-doped stabilizer Hamiltonians, although without discussing convergence guarantees.

In addition to the defected Toric code, Theorem [thm:intro4local] also works for more general families of qubit CLHs and can prepare the corresponding Gibbs state as long as there exists an efficient algorithm for the corresponding classical Gibbs sampling task.

1.2 Technical overview↩︎

Recall that in [thm:intro2local], we construct a Gibbs sampling reduction from 2-local qudit CLHs to 2-local qudit classical Hamiltonians, The proof is primarily based on the Structure Lemma [38]–[41]. The Structure Lemma has been the principle tool in studying the complexity of CLHs and, intuitively, says that one can transform a 2-local qudit CLH \(H^{(2)}\) to a 2-local qudit classical Hamiltonian \(H^{(2c)}\) via a constant depth quantum circuit \(\mathcal{C}_H\). In other words, there is a one-to-one correspondence between the computational basis of the classical Hamiltonian \(H^{(2c)}\) and the eigenstates of the quantum Hamiltonian \(H^{(2)}\). By this observation, there is a simple procedure to sample from the Gibbs state of \(H^{(2)}\), i.e. sample an eigenstate of \(H^{(2)}\) according to the Gibbs distribution. First, we sample a computational basis state \(\left\vert \psi \right\rangle\) from the Gibbs distribution of the classical Hamiltonian \(H^{(2c)}\) (via a classical Gibbs sampler). Then, applying a constant depth quantum circuit to \(\left\vert \psi \right\rangle\) yields an eigenstate of \(H^{(2)}\), distributed according to the Gibbs state of \(H^{(2)}\).

For Hamiltonians of higher locality, the exact correspondence present in the 2-local case does not hold. Nonetheless, we show in [thm:intro4local] that we can extend our techniques beyond 2-local Hamiltonians. We demonstrate a reduction from Gibbs sampling of 4-local qubit CLHs in 2D to classical Gibbs sampling.

The case “without classical qubits” is the simpler setting. Still, even in this case, we can no longer straightforwardly apply the Structure Lemma as is possible for 2-local Hamiltonians. This is not due to a deficiency in our techniques; rather, 4-local Hamiltonians can exhibit topological order [31] and there cannot be a constant depth quantum circuit \(\mathcal{C}_H\) as in the previous theorem. However, we observe that the eigenspace of qubit CLH is symmetric in some sense and via an oblivious randomized correction technique, we can adapt an algorithm for preparing ground state (as given in [40]) to preparing a Gibbs state. In particular, [40] proves that any 2D qubit CLH without classical qubits is equivalent to a defected Toric code permitting boundaries. That is, the “interior” terms look like Pauli X or Pauli Z terms and terms on the “boundary” have more freedom. The presence of boundaries makes it non-trivial to utilize this equivalence to design a Gibbs sampler for general 2D qubit CLH. We will use the defected Toric code as an example to explain our Gibbs sampler in Section 1.3.

When the initial Hamiltonian has classical qubits, the situation becomes more complex, as the connection from [40] between 2D qubit CLH and the Toric code only applies when there are no classical qubits. This does not pose a problem in [40] as they simply want to verify ground energy, and an \(\textsf{NP}\) prover can provide a recursive restriction of classical qubits consistent with some ground state. This effectively removes all classical qubits and the resulting Hamiltonian can be translated into a defected Toric code permitting boundaries.

In our case, we would like to recover the distribution over eigenstates, and thus we cannot perform the same recursive restriction. Additionally, we need our reduction to be efficient and should not depend on the power of an prover. We develop a propagation lemma to characterize the limits of the recursive restriction. Combined with an assumption that all fully quantum terms⁷ are uniformly correctable (see 1), we will argue that the statement of [40] can be modified to obtain a Gibbs sampling reduction from any 2D qubit CLH to a constant-locality classical Hamiltonian.

1.2.0.1 Extensions to more general 2D Hamiltonians

For simplicity, the theorems above are proved in the setting when the underlying Hamiltonian is placed on a \(2D\) lattice. However, in [40] the authors consider a more general setting of Hamiltonians on polygonal complexes. A straightforward generalization of our proofs works in this setting as well.

1.3 Case study: the punctured defected Toric code↩︎

To illustrate our Gibbs sampler for 2D qubit CLH, that is [thm:intro4local], we consider the restricted setting of the defected Toric code. In the remainder of this section, we assume that the inverse temperature is a constant \(\beta <+\infty\). We will formally define the defected Toric code and first for the punctured Toric code (defined later) we will give a \(\mathcal{O}(n^2)\)-time algorithm to prepare its Gibbs states via an oblivious randomized correction idea. The algorithm for Toric code on torus use similar ideas and is of runtime \(\mathcal{O}(n^ 2poly(\log n))\), which is put in Appendix 8. Those algorithms are specific to the defected Toric code. Then we describe a slightly different algorithm which is not as fast as the first algorithm, but by using the tools from [40] it can be extended to prepare Gibbs state for general qubit 2D CLHs.

The defected Toric code \(H_{DT}\) is embedded on a 2D, \(L \times L\) square lattice, with qubits placed on the vertices. Terms are grouped into “black” terms \(\mathcal{B}\) and “white” terms \(\mathcal{W}\), as in [fig:intro95fig95defected]. Formally, we define \[\begin{align} &H_{DT} =\sum_{p\in \mathcal{B}} c_pX^p + \sum_{p\in \mathcal{W}} c_pZ^p, \label{eq:HDT} \end{align}\tag{1}\] where \(X\) and \(Z\) are the standard Pauli \(X\) and \(Z\) operators. For a given term \(p\) acting on qubits \(q_1, \dots, q_4\), \(X^p\) denotes \(X_{q_1} \otimes \dots \otimes X_{q_4}\) (and same for \(Z^p\)). The coefficients \(c_p\) can be any real number (whereas \(c_p = -1\) in the standard Toric code).

1.3.1 A \(\mathcal{O}(n^2)\)-time algorithm for the punctured defected Toric code.↩︎

First, we explain the \(\mathcal{O}(n^2)\)-time algorithm to prepare the Gibbs states of the punctured defected Toric code. Specifically, the defected Toric code is punctured with a white term \(p_w\) and a black term \(p_b\) missing from the Hamiltonian. Thus, for any term \(p'\), there is a correction operator \(L_{p'}\), which anti-commutes with \(p'\) and commutes with all other terms. If \(p'\) is a white term, \(L_{p'}\) is realized as a string of Pauli \(X\) operators, starting in the support of \(p'\) and ending adjacent to the missing white term \(p_w\). Similarly if \(p'\) is a black term, \(L_{p'}\) is a string of \(Z\)’s from \(p'\) to \(p_b\). To prepare the Gibbs states for the punctured Hamiltonian \(H'_{DT}\), we initialize our state as the maximally mixed state, then sequentially measure and correct each plaquette term \(p\) in \(H'_{DT}\). That is, if we measure the plaquette operator \(p\) and get measurement outcome \(\lambda\in \{+c_p,-c_p\}\), then we perform the following oblivious randomized correction:

• With probability \(prob:=\frac{\exp(-\beta \lambda)}{\exp(\beta \lambda)+\exp(-\beta \lambda)}\) we do nothing.

• With probability \(1-prob\) we apply the correction operation \(L_{p}\).

The algorithm correctly prepares the Gibbs states of \(H'_{DT}\) because \(L_p\) bijectively maps the eigenspace associated with measurement outcome \(+c_p\) to that of \(-c_p\) and vice versa. A more detailed description of this procedure and proof of its correctness can be found in 4.2. Furthermore, the above process only performs the correction once for each plaquette term, and thus should not be interpreted as an iterative, randomized Accept/Reject process as done in the general Metropolis algorithm.

In the case where the Toric code has no punctured terms—such as the standard Toric code defined on a torus in the error-correcting code literature—a slightly different algorithm is required. In this case, correction operators are strings of Paulis between two measured terms \(p_1, p_2\). The above outline does not work directly for the following reason. Suppose \(p_1, p_2\) are measured and we obtain eigenvalues \((\lambda_1, \lambda_2) \in \{\pm c_p\}^{\otimes 2}\). Although, as above, we can bijectively map between \((c_p, c_p) \rightarrow (-c_p, -c_p)\) and for \((c_p, -c_p) \rightarrow (-c_p, c_p)\), the same is not possible for \((c_p, c_p) \rightarrow (c_p, -c_p)\) and vice versa. Instead, by adapting similar ideas, we show that there exists a Gibbs sampling reduction from the standard Toric code to 1D Ising model. Further details are provided in 8.

1.3.2 A generalizable algorithm↩︎

For general 2D qubit CLH (without classical qubits), not every plaquette term has a correction operator; in fact, the presence of a correction operator qualitatively characterizes the correctable interior terms and the non-correctable boundary terms. The proof that the interior terms are correctable uses techniques from [40], and the exterior terms are handled by a reduction to classical Gibbs sampling. The details are as below. For simplicity, we assume in this section that \(H_{DT}\) is embedded on a plane rather than torus, and the initial boundary of the lattice naturally plays the roles of the puncture terms in \(H_{DT}\).

1.3.2.1 Reduction to Classical Hamiltonian

As in [40] the first step is to remove enough terms such that the resulting Hamiltonian can be viewed as \(2\)-local. In the case of \(H_{DT}\), we can simply remove alternating rows of white terms, as in [fig:intro95fig95295local]. Finally, grouping all qubits on a white term as a single \(2^4\)-dimensional qudit, we see that the white terms become \(1\)-local and the black terms are all \(2\)-local. In [fig:intro95fig95295local], we group qubits \(u, v, w, \tau\) to form the qudit \(p_1\). Similarly we form the qudits \(p_2,p_3,p_4\). Then, the black term to the right of qudit \(p_1\) becomes 2-local, acting on \(p_1\) and \(p_3\). Call this \(2\)-local Hamiltonian \(H_{DT}^{(2)}\). The Structure Lemma of [38] gives a way to transform \(H_{DT}^{(2)}\) to a 2-local classical Hamiltonian. By working out the details (see 7.0.0.2), it turns out that in this 2-local classical Hamiltonian we obtain three distinct “types” of terms:

• \(h_\text{vert}\), 2-local terms corresponding to black terms acting on vertically arranged white terms (e.g. between \(p_1\) and \(p_2\)),

• \(h_\text{horiz}\), 2-local terms corresponding to black terms acting on horizontally arranged white terms (e.g. between \(p_1\) and \(p_3\)), and

• \(h_w\), 1-local terms corresponding to the white terms.

The final classical Hamiltonian is then \[H^{(2c)}_{DT} = \sum_{\text{vertical} p_i, p_j} (h_\text{vert})_{p_i,p_j} + \sum_{\text{horizontal} p_i,p_j} (h_\text{horiz})_{p_i,p_j} + \sum_p (h_w)_{p}.\]

1.3.2.2 Preparation of Quantum Gibbs State.

So far, we’ve removed terms from \(H_{DT}\) to obtain \(H^{(2)}_{DT}\), then argued that we can view this as a classical Hamiltonian \(H^{(2c)}_{DT}\). Assume we are able to perform classical Gibbs sampling at a given temperature on \(H^{(2c)}_{DT}\). To obtain a sampler for our original Hamiltonian, we need to reverse each step of the reduction. First, since the transformation from \(H_{DT}^{(2)}\) to \(H_{DT}^{(2c)}\) is via a low-depth quantum circuit, we can easily obtain a Gibbs state of \(H_{DT}^{(2)}\) from the Gibbs state on \(H_{DT}^{(2c)}\). The primary challenge is to correct for the terms we have removed to make \(H_{DT}\) \(2\)-local.

Suppose we want to correct for a remove white term \(p \in \mathcal{W}\). We first measure the current state \(\psi\) with respect to the term \(p\). If we only need to obtain some state with the correct eigenvalues, whenever we obtain an incorrect outcome, we could simply perform the correction operator \(L_p\) depicted in [fig:intro95fig95correct]. To obtain \(L_p\), we find a path from a corner of \(p\) to the boundary of the lattice, and apply a Pauli \(X\) on each qubit along the path. However, this does not immediately work when we are trying to sample from the Gibbs distribution.

Denote \(\Pi^p_{+c_p} \left\vert \psi \right\rangle\) be state if we get measurement outcome \(+c_p\) when measuring \(p\). Similarly for \(\Pi^p_{-c_p} \left\vert \psi \right\rangle\). There are two challenges in preparing the Gibbs state. First, we need to maintain the proper distribution over \(\Pi^p_{+c_p} \left\vert \psi \right\rangle\) and \(\Pi^p_{-c_p} \left\vert \psi \right\rangle\). Second, applying the correction \(L_p\) after measuring \(\Pi^p_{+c_p}\) may not yield \(\Pi^p_{-c_p} \left\vert \psi \right\rangle\), i.e. \[\label{eq:eigenstate95not95equal} L_{p} \Pi^p_{+c_p} \left\vert \phi(\boldsymbol{y}) \right\rangle \not\propto \Pi^p_{-c_p} \left\vert \phi(\boldsymbol{y}) \right\rangle .\tag{2}\] Nonetheless, we show that this can be done via an oblivious randomized correction technique. That is, based on the measurement outcome, we apply the correction operation \(L_p\) with some probability \(\mu\). At a high level, the correctness of the idea comes from the symmetry of the eigenspace; despite 2 , we do have that \[L_p \Pi^p_{+c_p} \Pi_{\lambda(\psi)} \Pi^p_{+c_p} L_p = \Pi^p_{-c_p} \Pi_{\lambda(\psi)} \Pi^p_{-c_p}\] where \(\Pi_{\lambda(\psi)}\) is an eigenspace of the non-removed operators corresponding to \(\left\vert \psi \right\rangle\). We will leverage this fact by applying a correction uniformly across this eigenspace, and the correction probability \(\mu\) will depend only \(\Pi_{\lambda(\psi)}\) rather than \(\left\vert \psi \right\rangle\) itself.

1.4 Conclusion and future work↩︎

In this manuscript, we give a reduction from Gibbs state preparation for various families of CLHs to the task of Gibbs sampling for classical Hamiltonians. In particular, based on the Structure Lemma we show that there is a Gibbs sampling reduction from 2-local qudit CLHs to 2-local qudit classical Hamiltonians. Based on the symmetry in qubit CLH and the idea of oblivious randomized correction we give a Gibbs sampling reduction from various 2D qubit CLH to qudit classical Hamiltonians. This approach yields a Gibbs sampler based on techniques very different than those traditionally used, such as analyzing the Davies generator. We also demonstrate that combined with existing fast mixing results for classical Hamiltonians, our Gibbs sampler matches the performance of state-of-the-art results in [1]–[3]. It will be interesting to apply our framework to more examples.

A natural direction to explore is whether our reduction can be generalized to other CLHs, especially those for which the complexity of proving ground energy is in \(\textsf{NP}\), such as the factorized qudit CLH on 2D lattice [39], the factorized CLH on any geometry [38] and the qutrit CLH on 2D [39].

Our work also gives an interesting characterization for the complexity of quantum approximate counting for specific CLHs. It is well-known that classical approximate counting is in \(\textsf{BPP}^{\textsf{NP}}\) [4]. Due to the connection between quantum approximate counting and the Gibbs state preparation [42], our work shows that quantum approximate counting for various CLHs is contained in the complexity class \(\textsf{BQP}^{\textsf{CS}}\), where \(\textsf{CS}\) is an oracle which can perform arbitrary classical Gibbs sampling. It would be interesting to explore whether there exist other families of quantum Hamiltonians where one can also upper bound the complexity of quantum approximate counting by \(\textsf{BQP}^{\mathcal{O}}\) for some oracle \(\mathcal{O}\) which is weaker than quantum approximate counting.

2 Preliminary↩︎

2.1 Notation↩︎

For two operators \(h\) and \(h'\), we use \([h,h']\) to denote the commutator \(hh'-h'h\). We say that \(h\) and \(h'\) commutes if \([h,h']=0\). Given two \(n\)-qubit quantum states \(\rho\) and \(\sigma\), we use \(\|\rho-\sigma\|_1:= \frac{1}{2}\text{tr}(|\rho-\sigma|)\) to denote their trace distance. Given two probability distributions \(\mathcal{D}_1\) and \(\mathcal{D}_2\) over \(\{0,1\}^n\), we use \(\|\mathcal{D}_1-\mathcal{D}_2\|_1\) to denote the total variation distance, that is \(\|\mathcal{D}_1-\mathcal{D}_2\|_1 = \frac{1}{2}\sum_{x\in\{0,1\}^n} |\mathcal{D}_1(x)-\mathcal{D}_2(x)|\), where \(\mathcal{D}_i(x)\) is the probability of sampling \(x\) in distribution \(\mathcal{D}_i\). Given a graph \(G=(V,E)\), for any vertex \(v\in V\), we use \(N(v)\) to denote the set of vertices which are adjacent to \(v\) (excluding \(v\)). For \(v,w\in V\), we use \(\{v,w\}\) and \(\langle v, w \rangle\) for unordered set and ordered set respectively. For a positive integer \(m\in \mathbb{N}\), we use \([m]\) to denote \(\{1,2...,m\}.\)

2.2 Formal problem definitions↩︎

\(k\)-local Hamiltonians. We say an \(n\)-qudit Hermitian operator \(H\) is a \(k\)-local Hamiltonian, if \(H=\sum_{i=1}^{m} h_i\) for \(m=poly(n)\), and each \(h_i\) only acts non-trivially on at most \(k\) qudits. We allow different qudits to have different dimensions.

For the special case of \(k=2\), one can define \(H\) via a graph \(G=(V,E)\). That is, on each vertex there is a qudit, and on each edge \(\{v,w\}\) there is a Hermitian term \(h^{vw}\), such that \[H_G=\sum_{\{v,w\}\in E} h^{vw}.\]

2.2.0.1 2D Hamiltonians.

Consider a 2D lattice \(G=(V,E)\) as in 3, where qudits are placed on vertices. As above, we can define a \(2\)-local Hamiltonian on \(G\) via \(H_G = \sum_{\{v,w\} \in E} h^{v,w}\). We can also define a \(4\)-local Hamiltonian over the 2D lattice by associating a Hermitian term to each plaquette \(P\). We use \(v\in P\) to denote that vertex \(v\) is in the plaquette \(P\). With some abuse of notations, we also use \(v\) and \(P\) to denote the corresponding qubit and the Hermitian term. The 2D, 4-local Hamiltonian on the lattice is given by \[H = \sum_P P\,.\]

For many of the proofs in this work, it will useful to partition the plaquette terms \(\{P\}_P\) into a set of “black” terms \(\mathcal{B}\) and “white” terms \(\mathcal{W}\) by viewing the 2D lattice as a chess board, as in 3.

There is another natural notion of a 2D local Hamiltonian where the qudits are placed on the edges and Hermitian terms are corresponding to plaquettes and stars; this is the version which is primarily considered in [40]. However, the two settings (where qudits are on vertices or are on edges) are equivalent when the underlying graph is the 2D square lattice (See Appendix C in [39]).

2.2.0.2 Commuting and classical Hamiltonians.

We say a \(k\)-local Hamiltonian \(H = \sum_{i=1}^m h_i\) is a commuting local Hamiltonian (CLH) if \([h_i,h_j]=0,\forall i,j\). Whenever we have a Hamiltonian \(H=\sum_P P\) defined over the plaquettes of a 2D lattice, we have \([p,p']=0\) \(\forall p,p'\). We say a \(k\)-local Hamiltonian \(H = \sum_{i=1}^m h_i\) is classical, if each \(h_i\) is diagonalized in the computational basis.

2.2.0.3 Defected Toric code and 2D Ising model.

Here we give two examples of qubit CLHs on 2D. For a vertex \(v\) in plaquette \(p\), we use \(Z^p_v,X^p_v\) to denote the Pauli Z and Pauli X operator on the qubit \(v\). When \(v\) uniquely identifies a vertex we abbreviate \(Z^p_v\) as \(Z_v\) and similarly for other Pauli operators. For a plaquette term \(p\), we define \(Z^p:= \otimes_{v\in p} Z_v\) and \(X^p\) similarly.

As shown in [fig:toric95code], the defected Toric code is defined as \[H=\sum_{p\in \mathcal{W}} c_p Z^p + \sum_{p\in \mathcal{B}} c_p X^p,\] where \(c_p\in \mathbb{R}\) can be arbitrary. The standard Toric code is a special case when all \(c_p=-1\).

Denote the 2D lattice as \(G=(V,E)\). As in [fig:ising95model], the ferromagnetic 2D Ising model is a 2-local Hamiltonian \[H=\sum_{\{u,v\}\in E} Z_u\otimes Z_v.\]

The 2D Ising model is a classical Hamiltonian whose eigenstates are all computational basis. The defected Toric code is not a classical Hamiltonian. The ground state of the standard Toric code is highly entangled and cannot be prepared by any constant depth quantum circuit.

2.2.0.4 Gibbs state preparation.

Given a \(k\)-local Hamiltonian \(H=\sum_i h_i\), an inverse temperature \(\beta<+\infty\), the Gibbs state with respect to \(\left(H,\beta\right)\) is defined as \[\begin{align} \rho( H,\beta)=\frac{1}{\text{tr}(\exp(-\beta H))}\exp(-\beta H). \end{align}\] Given \(\epsilon> 0\), we say an algorithm \(\mathcal{A}\) prepares \(\rho( H,\beta)\) within precision \(\epsilon\), if \(\mathcal{A}\) outputs a state \(\rho\) such that \[\begin{align} \|\rho-\rho( H,\beta)\|_1 \leq \epsilon. \label{eq:rho} \end{align}\tag{3}\] When \(H\) is classical, the classical Gibbs distribution with respect to \(\left(H,\beta\right)\) is denoted as \(\mathcal{D}_{\beta H}\), which samples a classical string \(x\in\{0,1\}^n\) with probability \(\exp(-\beta \langle x|H|x\rangle)/ \text{tr}(\exp(-\beta H))\). We say an algorithm \(\mathcal{A}\) performs classical Gibbs sampling \(\mathcal{D}_{\beta H}\) with precision \(\epsilon\), if \(\mathcal{A}\) outputs a distribution \(\mathcal{D}\) such that \[\begin{align} \|\mathcal{D}-\mathcal{D}_{\beta H}\|_1\leq \epsilon. \label{eq:D} \end{align}\tag{4}\] Note that Eq. (4 ) is equivalent to Eq. (3 ) when \(\rho\) and \(\rho(H, \beta)\) are diagonal matrices.

3 Reduction for 2-local qudit CLHs↩︎

In this section, we prepare the Gibbs state for qudit 2-local CLHs. In particular, in Section 3.1 we will prove the Gibbs sampling reduction for general 2-local qudit CLHs in Theorem 1. Then in Section 3.2 we will give several examples, whose proofs are put into Appendix 9.

3.1 General case↩︎

Recall that a 2-local CLHs \(H_G^{(2)}\) is defined on a graph \(G=(V,E)\), where \[H_G^{(2)}=\sum_{\{v,w\}\in E} h^{vw} .\] Here \(\{h^{vw}\}_{\{v,w\}\in E}\) are Hermitian terms and commute with each other. The superscript \((2)\) is to emphasize that the Hamiltonian is 2-local. In this section, we assume on each vertex there is a qudit rather than a qubit, and we allow \(G\) to be an arbitrary graph rather than just a 2D lattice.

The Gibbs sampling reduction for 2-local CLHs comes from the Structure Lemma, which was originally developed by [38] to study the computational complexity of commuting Hamiltonians. A constructive proof of the Structure Lemma can be found in Section 7.3 of [56]. Intuitively, the Structure Lemma says that one can decouple all commuting 2-local terms. This will allow us to identify eigenstates of a 2-local CLH \(H_G^{(2)}\) with the computational basis of a classical Hamiltonian \(H_G^{(2c)}\) defined later. In addition, each such eigenstate can be prepared by a constant depth quantum circuit. Thus, to prepare the Gibbs state of \(H_G^{(2)}\), it suffices to first do classical Gibbs sampling for \(H_G^{(2c)}\), yielding a distribution over computational bassi states, then prepare the corresponding eigenstate of \(H_G^{(2)}\) indexed by a sampled basis state via a constant depth quantum circuit.

We first give the formal statement of the Structure Lemma.

The Structure Lemma can be understood via 4. We can understand ?? as the following: by choosing a proper local basis for the Hilbert space \(\mathcal{H}_{j_v}^v\), it is equivalent to the Hilbert space of \(|N(v)|\) distinct new qudits.

If for every qudit \(v\), one applies 4 and chooses an index \(j_v\in [J_v]\) and corresponding subspace \(\mathcal{H}^v_{j_v}\), this will decouple all terms in \(H_G\). Each term \(h^{vw}\) restricted to the subspaces \(\mathcal{H}^v_{j_v} \otimes \mathcal{H}^v_{j_w}\) will act on distinct qudits and \[h^{vw} \in \mathcal{L}\left(\mathcal{H}_{j_v}^{\langle v,w\rangle}\otimes \mathcal{H}_{j_w}^{\langle w,v\rangle}\right).\]

To define the corresponding classical Hamiltonian \(H_G^{(2c)}\), we use the indices \(\{j_v\}_{v \in V}\) to index the eigenstates of \(H_G^{(2)}\). Denote \[h^{vw}|_{j_vj_w}:= \text{ restriction of h^{vw} onto \mathcal{H}_{j_v}^{\langle v, w \rangle} \otimes \mathcal{H}_{j_w}^{\langle w,v \rangle}.}\] Note that eigenstates of \(h^{vw}|_{j_vj_w}\) might not be computational basis states (and in particular could be entangled). Nonetheless, we have shown that under the restriction corresponding to \(\{j_v\}_{v \in V}\), all terms \(h^{vw}|_{j_v, j_w}\) act on distinct qudits and we can use the computational basis to index the eigenstates. As shown in Figure 4, let \(D^{\langle v, w \rangle}_{j_v} := \dim (\mathcal{H}^{\langle v, w \rangle}_{j_v})\) and write the basis of each subspace \(\mathcal{H}_{j_v}^{\langle v, w \rangle}\) as \(\left\vert b_{j_v}^{\langle v, w \rangle} \right\rangle\) where \(b_{j_v}^{\langle v, w \rangle}\) ranges over \([D_{j_v}^{\langle v, w \rangle}]\). Then \(\bigotimes_{w\in N(v)} \left\vert b_{j_v}^{\langle v, w \rangle} \right\rangle\) is a computational basis state of \(\mathcal{H}_{j_v}^{v}\). For each edge \(\{v,w\}\in E\), the term \(h^{vw}|_{j_vj_w}\) is Hermitian and thus can be diagonalized; the computational basis states \(\left\vert {\boldsymbol{b}}_{j_vj_w}^{vw} \right\rangle\) are used to index the eigenstates. Thus, a basis for the full eigenspace of \(h^{v,w}|_{j_v,j_w}\) is given by \[\left\vert {\boldsymbol{b}}_{j_vj_w}^{vw} \right\rangle : = \left\vert b_{j_v}^{\langle v, w \rangle},b_{j_w}^{\langle w,v \rangle} \right\rangle,\quad {\boldsymbol{b}}^{v,w}_{j_v j_w} \in [D^{\langle v,w \rangle}_{j_v} \times D^{\langle w,v \rangle}_{j,w}]\,.\]

Given an index \({\boldsymbol{b}}^{v,w}_{j_v j_w}\), the corresponding eigenstate is denoted \(\psi({\boldsymbol{b}}_{j_v j_w}^{vw})\) and the eigenvalue \(\lambda({\boldsymbol{b}}_{j_v j_w}^{vw})\). The classical Hamiltonian is defined by substituting the eigenstate with its index, that is \[\begin{align} \label{eq:classical95ham} H_G^{(2c)} := \sum_{\{v,w\}\in E} \quad \sum_{j_v,j_w} \quad \sum_{{\boldsymbol{b}}_{j_vj_w}^{vw}} \lambda({\boldsymbol{b}}_{j_vj_w}^{vw}) \left\vert {\boldsymbol{b}}_{j_vj_w}^{vw} \right\rangle \left\langle {\boldsymbol{b}}_{j_vj_w}^{vw} \right\vert. \end{align}\tag{5}\] Following the usual convention, each term in the summand is implicitly padded with identities as necessary.

In this way, the eigenstates of \(H^{(2c)}_G\) are given by specifying an index \(j_v\) for each vertex \(v \in V\), then a basis state \(b^{\langle v, w \rangle}_{j_v j_w}\) for each of the decoupled Hilbert spaces \(\mathcal{H}^{\langle v, w\rangle}_{j_v} \otimes \mathcal{H}^{\langle w,v\rangle}_{j_w}\) \[\begin{align} &{\boldsymbol{j}}:=\{j_v\}_v,\label{eq:14}\\ & {\boldsymbol{b}}_{\boldsymbol{j}}:=\{b_{{\boldsymbol{j}}_v, {\boldsymbol{j}}_w}^{\langle v, w \rangle}\}_{\langle v, w \rangle}. \end{align}\tag{6}\] Thus, by construction, we have the following.

We can also easily map eigenstates of \(H^{(2c)}_G\) to eigenstates of \(H^{(2)}_G\) via \[\begin{align} &\lambda({\boldsymbol{b}}_{\boldsymbol{j}}): = \sum_{\{v,w\}\in E} \lambda({\boldsymbol{b}}_{j_vj_w}^{vw})\\ & \left\vert \psi({\boldsymbol{b}}_{\boldsymbol{j}}) \right\rangle : = \bigotimes_{\{v,w\}\in E } \left\vert \psi({\boldsymbol{b}}_{j_vj_w}^{vw}) \right\rangle\,,\label{eq:17} \end{align}\tag{7}\] and any classical Gibbs sampling procedure for \(H^{(2c)}_G\) yields a Gibbs sampler for \(H^{(2)}_G\).

3.2 Examples↩︎

In this section, we write down the Gibbs sampling reduction for some specific Hamiltonians as illustrative examples. All the proofs are deferred to 9.

We first consider an \(n\)-qudit Hamiltonian on a 1D chain \[\begin{align} H=\sum_{i} h_i,\label{eq:hi} \end{align}\tag{8}\] we say that \(H\) is \(r\)-range if \(h_i\) only acts non-trivially on qudits \(i,i+1,..,i+r-1\). For simplicity we assume \(n\) is an integer multiple of \(r\) and the qudit dimension \(d\) is a power of \(2\). We say that \(H\) is finite-range if \(r\) is a constant, and \(H\) is translation-invariant if all the terms \(h_i\) are the same.

By coarse-graining \(H\), we can always assume \(H\) is 2-local: group each consecutive set of \(r\) qudits as a new qudit so that each \(h_i\) acts non-trivially on at most two (grouped) qudits. For each pair of new qudits \(\{j, j+1\}\), we associate the new term \(H_{j,j+1}\), which is a sum of all terms from \(H\) acting on the corresponding qudits. Thus \(H\) can be viewed as a 2-local qudit Hamiltonian on 1D written as \(H=\sum_j H_{j,j+1}\). Note the terms \(H_{j,j+1}\) can also be made translation-invariant.

We also give another example of a Hamiltonian on a 2D lattice. Recall our first definition of a 2-local Hamiltonian on 2D (see 2.2) \[H_{2D}=\sum_{\{v,w\}\in E} h^{vw}.\] As usual, \(H_{2D}\) is translation-invariant if all terms \(h^{vw}\) are the same. We say that a 2D lattice has periodic boundary condition if it can be embedded onto torus; we will assume a periodic boundary for simplicity.

4 Reduction for 2D 4-local qubit CLH without classical qubits↩︎

In this section we describe how to prepare the Gibbs state of 2D qubit CLHs. In 4.1 we first review the canonical form of the 2D qubit CLH as developed in [40], who establishes a connection between 2D qubit CLHs and the defected Toric code. Based on this connection and an observation on the symmetry of the ground space, we use an oblivious randomized correction technique to generalize the Gibbs state preparation algorithm for the defected Toric code presented in 1.2 to prepare the Gibbs state for the more general family of 2D qubit CLHs without classical qubits.⁸

4.1 A canonical form↩︎

This section is primarily a review of the results in [40]. In that work, the authors prove that 2D qubit CLH⁹ without classical qubits (which we will define shortly) is in some sense equivalent to the defected Toric code. [40] used this connection to show that one can prepare the ground state of 2D qubit CLHs similar to the way ground states of the defected Toric code are prepared; this is via the measure and correct approach mentioned in 1.3.1.

We summarize necessary definitions and theorems which will be used in later sections. Recall that a 2D qubit CLH is defined as \(H=\sum_{p\in P} p\), where \(P\) is the set of plaquettes of the lattice.

In other words a qubit \(q\) is classical if under some basis for the qubit, all terms look like \(\vert 0\rangle\!\langle 0\vert_q \otimes ... + \left\vert 1 \right\rangle 1_q \otimes ...\).

Even when all terms are in the interior, Lemma 8 does not imply all the terms should be a tensor product of Pauli \(Z\). To be written in the form of Eq. (?? ), adjacent plaquettes may need different choices of basis. An example is the Toric code, where plaquettes are alternate \(X^{\otimes 4}\) and \(Z^{\otimes 4}\).

Recall the chessboard partition of the 2D lattice into black plaquettes \(\mathcal{B}\) and white plaquettes \(\mathcal{W}\).

In [40] the authors use Lemma 9 to reduce the task of preparing the ground state of any (4-local) 2D qubit CLH to the task of preparing the ground state of a 2-local qudit CLH. This reduction is characterized by the following Corollary.

For instance, in 2 from the technical overview, \(\mathcal{P}\) is the set of all black terms \(\mathcal{B}\) and non-removed white terms \(\mathcal{O}\). Then, \(H_{\mathcal{P}}\) is exactly the 2-local Hamiltonian \(H_{DT}^{(2)}\).

4.2 Reduction to a classical Hamiltonian↩︎

In this section, we describe how to reduce the task of Gibbs state preparation for qubit CLH without classical qubits to the task of classical Gibbs sampling. An explicit and canonical example is the Gibbs state preparation for the defected Toric code, as described in 7 and 8, with respect to the punctured defected Toric code (e.g. embedded on a planar lattice) and on a torus, respectively. For general qubit CLHs without classical qubits, our result is summarized in the following theorem.

We begin with some notation. Recall that the set of all plaquette terms \(P\) is partition into black and white terms, i.e. \(P = \mathcal{B}\cup\mathcal{W}\). In 1 we defined \(\mathcal{R}\) as the set of terms which have access to the boundary. Let \(Q\subseteq \mathcal{B}\cup\mathcal{W}\) be an arbitrary subset of the plaquette terms. We write \(\mathcal{R}\setminus Q\) as the set difference of \(\mathcal{R}\) and \(Q\). For \(p\in \mathcal{R}\) define the correction operator \(L_p\) as in 4.

Define \({\boldsymbol{\lambda}}_Q:=\{\lambda_p\}_{p\in Q}\) to be a set of real values, where each \(\lambda_p\) coresponds to an eigenvalue of \(p\in Q\). Let \(\lambda(Q) :=\sum_{p\in Q} \lambda_p\). Recall that all plaquette terms are commuting. Thus, the terms are simulataneously diagonalizable and the common eigenspace of each \(p \in Q\) with eigenvalue \(\lambda_p\) is well defined; we denote this as \(\mathcal{H}^Q_{{\boldsymbol{\lambda}}_Q}\). Formally, \[\begin{align} \mathcal{H}^Q_{{\boldsymbol{\lambda}}_Q}:= \{\left\vert \phi \right\rangle \, | \, p\left\vert \phi \right\rangle = \lambda_p \left\vert \phi \right\rangle , \forall p\in Q\}. \end{align}\] Let \(\Pi^Q_{{\boldsymbol{\lambda}}_Q}\) be the projection onto \(\mathcal{H}^Q_{{\boldsymbol{\lambda}}_Q}\). Again by commutation, \(\Pi^Q_{{\boldsymbol{\lambda}}_Q}\) is equal to the product of the individual projectors: \[\begin{align} \Pi^Q_{{\boldsymbol{\lambda}}_Q} = \prod_{p\in Q} \Pi^p_{\lambda_p}\,.\label{eq:proj} \end{align}\tag{9}\]

Note that for a general 2D qubit CLH, a plaquette term \(p\) can be any arbitrary 4-qubit operator. For example one can set \(H=p_0 + \mathcal{I}\) where \(p_0\) is an arbitrary operator on one plaquette. Nonetheless, one can show that all terms in \(\mathcal{R}\) are in a sense quite regular.

Additionally, the eigenspaces corresponding to each eigenvalue of \(p\) are symmetric. 11 is the key observation which leads to the oblivious randomized correction idea.

A consequence of this lemma is that the the the eigenspace corresponding to \({\boldsymbol{\lambda}}_Q\) is balanced across any \(p\)’s \(+c_p\) and \(-c_p\) eigenspaces.

Underyling the 2 is the following algorithm. We will prove 2 by proving correctness via 13. Recall that by 1 the Hamiltonian \(H_G^{(2)}\) is 2-local after we group some qubits into qudits. Using the notation from 3, \(H_G^{(2)}\) and \(H_G^{(2c)}\) denote the \(2\)-local CLH and the corresponding classical Hamiltonian, and \(\{{\boldsymbol{b}}_{\boldsymbol{j}}\}_{{\boldsymbol{b}}_{\boldsymbol{j}}}\) denote the computational basis of the grouped qudits. Since \(G\) is planar (1 item (1)) the number of edges \(m\) is \(\mathcal{O}(n)\), with \(n\) being the number of vertices. As usual, we assume access to a classical Gibbs sampler which obtains a computational basis state \(\left\vert \psi({\boldsymbol{b}}_{\boldsymbol{j}}) \right\rangle\) with probability \(p({\boldsymbol{b}}_{\boldsymbol{j}})\) in time \(T+ \mathcal{O}(m)=T+ \mathcal{O}(n)\). That is, we prepare a state \[\begin{align} &\rho(\mathcal{P}):= \sum_{{\boldsymbol{b}}_{\boldsymbol{j}}} p({\boldsymbol{b}}_{\boldsymbol{j}}) \left\vert \psi({\boldsymbol{b}}_{\boldsymbol{j}}) \right\rangle\left\langle \psi({\boldsymbol{b}}_{\boldsymbol{j}}) \right\vert\\ \text{such that }&\left\| \rho(\mathcal{P}) - \rho(H_G^{(2)},\beta)\right\|_1 \leq \epsilon\,. \label{eq:gibbs} \end{align}\tag{11}\] Here we did not write down the explicit formula for \(p({\boldsymbol{b}}_{\boldsymbol{j}})\) since we will not use it. In the second step of the algorithm, we sequentially measure the current state with respect to each removed term \(p \in \mathcal{R}\) and perform an oblivious randomized correction. The details are in 5.

We prove correctness by induction. Define \(H_Q:=\sum_{p\in Q} p\). We claim the following.

5 Reduction for 2D (4-local) qubit CLH with classical qubit↩︎

5.1 A review of the techniques of Aharonov et. al [40]↩︎

Recall that [40] studied the structure of qubit 2D commuting local Hamiltonians to argue that preparing the ground state can be done in \(\textsf{NP}\). The two main technical results of our work ([thm:4local2,thm:gibbs_with_classical]) require opening up this result so that we are not only able to prepare a single ground state, but sample from the Gibbs state of the Hamiltonian. This requires a good understanding of their correction operators, as well as the role classical qubits play in their proof. As such, we dedicate this section to reviewing the techniques used in their paper.

5.1.1 \(\mathrm C^\star\)-algebras and the Structure Lemma↩︎

The primary technical tool used in [40] (and in nearly all other works on commuting local Hamiltonians [38], [39], [41]) is the Structure Lemma for \(\mathrm C^\star\) algebras.

The connection between local Hamiltonians and algebras is made through the concept of an “induced algebra”.

The following lemma tells us that the induced algebra is independent of the choice of basis for \(q_2\).

The induced algebra gives us a tool to analyze whether two terms commute.

Finally, we recall the Structure Lemma, which was first applied to understanding commuting local Hamiltonians by [38]. Since then, it has been the principal tool used by subsequent works on the CLH [39], [41], [57], [58]. At a high level, the Structure Lemma says that algebras can be block-diagonalized, and that within each of these blocks, the algebra takes on a tensor product structure which identifies its commutant. For a proof of the lemma, see [56].

We remark that 16 is equivalent to 4 in 3.

A corollary of 16, and the reason why it is so useful for characterizing the properties of commuting local Hamiltonians, is that in order to commute with an algebra, another algebra must live entirely within the \(\mathcal{H}^q_{(i,2)}\) subspaces. Formally, we have the following.

5.1.2 Restrictions on 2D commuting local Hamiltonians↩︎

With the Structure Lemma in hand, we can see how [40] apply these tools to show that 2D qubit CLH is in NP.

We start with some definitions. First, we review the notion of a classical qubit, defined originally in 3. For a qubit \(q\) and a commuting Hamiltonian \(H\), let \(\mathcal{N}(q)\) be the set of terms acting non-trivially on \(q\). We say that \(q\) is classical if there is a non-trivial decomposition of \(\mathcal{H}^q\) (i.e. \(\mathcal{H}^q = \bigoplus_{i \in \ell} \mathcal{H}^q_i\), with \(\ell > 1\)) and each term \(h \in \mathcal{N}(q)\) keeps this decomposition invariant. We define \(\mathcal{C}_0(H)\) as the set of classical qubits of \(H\). If the Hamiltonian \(H\) is clear from context, we abbreviate \(\mathcal{C}_0(H)\) as \(\mathcal{C}_0\).

A key point is that an initial restriction \(H|_s\), with \(s \in \{0,1\}^{\mathcal{C}_0}\) can lead to the creation of new classical qubits in \(H|_s\). For instance, in 6, setting \(s_c = \left\vert 1 \right\rangle\) removes \(h\) from the Hamiltonian and causes \(q\) to become classical.

When we refer to a terminating restriction, that restriction is implictly assumed to be valid.

In other words, \(\mathcal{C}\) is the set of qubits which might become classical when we sequentially choose an assignment for the current classical qubits. For example in Figure 6 \(q\) is a “possible classical qubit” and both \(c,q\in \mathcal{C}\).

Now that we have characterized all the possible classical qubits, we define the notion of a “fully quantum” qubit (and term).

With this notation, we can understand the first step of the algorithm of [40] as applying a valid, terminating restriction \(\boldsymbol{s} = (s_0,\dots, s_i)\) to the Hamiltonian.

5.2 Classical qudits↩︎

In this section we show a reduction from a CLH instance \(H\) to a classical Hamiltonian with constant locality, albeit with some (fully quantum) terms removed, but without needing to first remove classical qubits. In general, the choice of classical qubits can affect the correction operators for a removed term. To deal with this issue, we make the following assumption.

Formally, our result is stated as follows.

We prove this theorem in two steps. First, in 5.2.2 we show how to modify the proof of [40] so that classical qubits do not need to be removed when preparing a single ground state. Then in 5.2.3 we show how to extend this idea to Gibbs state sampling and get a proof for 3.

5.2.1 Characterization of classical qubits↩︎

Proving 3 will require characterizing the set of possible classical qubits \(\mathcal{C}\) (12). We begin with a lemma about algebra on qubits in the interior (recall 2).

This lemma provides strong restrictions on the possible set of classical qubits after “propagating” the initial set of classical qubits.

5.2.2 Uniform preparation of a ground state↩︎

In this section, we describe how to avoid restricting classical qubits in the proof of [40].

Let us first recall the high level proof of [40], which can be boiled down to:

1. The prover provides a terminating restriction \(\boldsymbol{s}\) such that \(H|_{\boldsymbol{s}}\) contains a ground state of \(H\). This removes all classical qubits.

2. Remove a set of terms \(\mathcal{R}\) of \(H|_{\boldsymbol{s}}\) so that the Hamiltonian becomes two local.

3. Prepare a ground state of the two-local Hamiltonian via [38] (with the help of the prover).

4. Correct the removed operators in \(\mathcal{R}\) to get a ground state for \(H\).

Our first observation is that [item:removal] can be performed without first removing classical qubits, as this step depends only on the geometric structure of the Hamiltonian. However, the issue comes in [item:correct]; depending on the choice of classical qubits, the set of “correctable” terms may be different. Therefore, if we perform [item:removal] without considering the classical restriction and the resulting set of correctable terms, we may not be able to correct all terms in \(\mathcal{R}\).

More precisely, [40] characterize the set of correctable terms via paths to the boundary.

These “paths to the boundary” yield correction operators, and the choice to remove \(h\) or \(h'\) depends on which term is correctable. For instance in 8, we see that our choice for the classical qubit \(c\) changes which of the two boxed terms may be removed.

To avoid this issue, we take 1, which states that if some \(h\) is correctable under some valid, terminating restriction \(H|_{\boldsymbol{s}}\), then it is correctable under any other valid, terminating restriction \(H|_{\boldsymbol{t}}\). This implies that the classical restriction can be deferred until after [item:prepare95gs] is performed. We will now formalize this idea.

5.2.2.1 Choice of Removed Terms.

In the proof of [40], the authors triangulate the 2D complex, then construct a “co-triangulation”, dividing the surface in tiles \(T \in \mathcal{T}\) such that each qubit within a single tile \(T\) becomes a new qudit \(q_T\) in the transformed Hamiltonian (see 9). Under this transformation, terms in the original Hamiltonian are either 1) internal to a single tile \(T\) (and are now 1-local), 2) cross between two adjacent tiles \(T, T'\) (and are now 2-local), or 3) on the corner of three tiles (and are now 3-local). Since the goal is to produce a \(2\)-local Hamiltonian, the natural strategy is to simply set the set of removed terms \(\mathcal{R}\) to be precisely these corner terms. However, we need to be a bit careful, since any given term is not necessarily correctable; all we know from 19 is that either the corner term \(h\) or its neighbor \(h'\) is correctable.¹¹ This is easy to handle. As long as the original triangulation has sufficient girth, the co-triangulation can be “shifted” so that its corners lie entirely in correctable terms.

In our case, we no longer have the ambiguity of which terms is correctable, but we instead need to deal with classical terms carefully. Consider a triangulation and corresponding tiling \(\mathcal{T}\). Let \(h\) be an arbitrary term containing a corner of the tiling \(\mathcal{T}\). We consider the following set of cases:

• (Case 1) \(h\) is a fully quantum term (14).

  • (Case 1a) \(h\) is an interior term.

  • (Case 1b) \(h\) is a boundary term.

• (Case 2) \(h\) is supported on an originally classical qubit \(q \in \mathcal{C}_0\).

• (Case 3) \(h\) is supported on a qubit \(q \in \mathcal{C} \setminus \mathcal{C}_0\) (i.e. a qubit which is classical under some restriction but which is not classical in the original Hamiltonian).

In the first case, we apply the same logic as [40]: for Case 1a, \(h\) is an interior term, it will be correctable (in our case by 1) and we added it to the removed set \(\mathcal{R}\); in Case 1b, \(h\) is a boundary term which implies a neighbor acts trivially on one of its qubits. This neighbor results in a “hole” in the original surface. Place the corner of the co-triangulation in the hole.

In Case 2, we construct the co-triangulation such that one tile contains only the classical qubits of \(h\); this will ensure that any assignment to the classical qubit renders this term \(2\)-local as well. For Case 3, we appeal to 18; since \(q\) is not originally classical but becomes classical, it must be on a boundary term or a term with a classical qubit. In the boundary case, we apply the same argument as in Case 1b. In the classical case, apply the argument from Case 2. Finally, we denote the Hamiltonian with terms in \(\mathcal{R}\) removed as \(\widetilde{H}\).

5.2.2.2 Construction of the Classical Hamiltonian.

Now that we have fixed the tiling \(\mathcal{T}\) and chosen the set of removed terms \(\mathcal{R}\), it remains to describe how to translate this into a classical Hamiltonian. To do so, we introduce some notation. For each tile \(T \in \mathcal{T}\), we denote \(Q(T)\) to refer to the qubits contained entirely within \(T\) and \(H(T)\) as the set of original Hamiltonian terms in \(\widetilde{H}\) overlapping \(T\) (i.e. acting non-trivially on some \(q \in Q(T)\)). In addition to the qubits internal to \(T\), the qubits “nearby” \(T\) are also important. Define \[\overline{Q}(T) = \bigcup_{h \in H(T)} Q(T),\] Additionally, let \(\mathcal{N}(T)\) be the tiles neighboring \(T\). Then, we define the Hilbert space of the tile qudit \(q_T\) as \(\mathcal{H}_T = \otimes_{q \in T} \mathcal{H}_q\). To define the terms in the grouped Hamiltonian, we define the following sets:

• \(S_T\) is the set of Hamiltonian terms which act non-trivially on \(Q(T) \setminus \mathcal{C}_0\) (i.e. on a qudit in \(T\) which is not originally classical).

• \(S_{T,T'} := S_T \cap S_{T'}\) are the terms acting on both \(T\) and \(T'\).

• \(h_{T,T'} = \sum_{h \in S_{T,T'}} h\) will be the 2-local terms in the grouped Hamiltonian.

• \(h_T = \sum_{\substack{h \in S_T\\\forall T' \neq T, h \not\in S_{T'}}} h\) are the 1-local terms.

Then, the grouped Hamiltonian is denoted as \(\widetilde{H}_{\mathcal{T}} = \sum_{T} h_T + \sum_{T' \neq T} h_{T,T'}\). Again, note that \(h_{T,T'}\) is technically not \(2\)-local; it only becomes \(2\)-local once a classical restriction to \(\mathcal{C}_0\) is fixed. But given such a restriction, we have the following simple consequence of the Structure Lemma [38].

We can say something slightly stronger than this; for a tile \(T\), the above decomposition only depends on terms \(h_{T,T'}\) which intersect \(T\) (and the internal terms as well). Therefore, if we consider a restriction \(s'\) where a bit outside of \(\overline{Q}(T)\) is flipped, this induces the same decomposition, i.e. \(\mathcal{D}^{s'}_T = \mathcal{D}^s_T\). As a result, when specifying a restriction in the context of a triangle, we implicitly imagine \(s\) is defined over \(\{0,1\}^{|Q(T) \cap \mathcal{C}_0|}\).

Recall our goal is to argue that \(h_{T,T'}|_s\) is classical. Recalling 3, we see that the classical terms constructed for a term \(h_{v,w}\) in 5 are only a function of the decomposition on vertices \(v\) and \(w\) (in our case \(q_T\) and \(q_{T'}\)). In particular, for a fixed assignment \(s\) and term \(h_{T,T'}|_s\), we obtain the classical term, \[\label{eq:dep95classical95ham} h_{T,T'}|_s := \sum_{\substack{\boldsymbol{j} = (j_1,j_2)\\ \in [\ell_T] \times [\ell_{T'}]}} \sum_{\boldsymbol{b}_{\boldsymbol{j}}^{T,T'}} \lambda(\boldsymbol{b}^{T,T'}_{\boldsymbol{j}}) \vert\psi(\boldsymbol{b}^{T,T'}_{\boldsymbol{j}})\rangle\!\langle\psi(\boldsymbol{b}^{T,T'}_{\boldsymbol{j}})\vert,\tag{21}\] where \(\boldsymbol{b}^{T,T'}_{\boldsymbol{j}}\) and \(\lambda(\boldsymbol{b}^{T,T'}_{\boldsymbol{j}})\) are defined as in , except that we use \(T,T'\) to refer to qudits, rather than \(v,w\). Then, we have that the following \(2\)-local Hamiltonian is equivalent to \(\widetilde{H}|_s\), \[\widetilde{H}^{(c,s)} := \sum_{T,T'} h_{T,T'}|_s + \sum_T h_T|_s.\] Now, we have that by construction, any ground state \(\left\vert \psi \right\rangle\) to \(\widetilde{H}^{(c,s)}\) is (after some 1-local unitary transformations) equal to a ground state \(\left\vert \phi \right\rangle\) of \(H|_s\); in turn, we can obtain a ground state of the original Hamiltonian \(\widetilde{H}\) as \(\left\vert s \right\rangle \otimes \left\vert \phi \right\rangle\). In particular, we can write down a Hamiltonian equivalent to \(\widetilde{H}\) as \[\sum_{s \in \{0,1\}^{|\mathcal{C}_0|}} \vert s\rangle\!\langle s\vert \otimes \widetilde{H}^{(c,s)} = \sum_{T,T'}\sum_{s \in \{0,1\}^{|\mathcal{C}_0|}} \vert s\rangle\!\langle s\vert \otimes h_{T,T'}|_s + \sum_T\sum_{s \in \{0,1\}^{|\mathcal{C}_0|}} \vert s\rangle\!\langle s\vert \otimes h_T|_s.\] Unfortunately, since \(|\mathcal{C}_0|\) can be as large as \(\textsf{poly}(n)\), this yields an exponentially-large Hamiltonian. But we use our earlier observation, which is that the decompositions on \(T\) and \(T\) only depend on \(\overline{Q}(T) \cap \mathcal{C}_0\). Since the local decompositions are the same, so too are classical terms obtained in 21 . This means that for any \(T,T'\), \[\sum_{s \in \{0,1\}^{|\mathcal{C}_0|}} \vert s\rangle\!\langle s\vert \otimes h^{(c,s)}_{T,T'} = \sum_{s' \in \{0,1\}^{|\mathcal{C}_0 \cap (S_T \cup S_{T'})|}} \vert s'\rangle\!\langle s'\vert \otimes h^{(c,s')}_{T,T'},\] and the resulting Hamiltonian is \[\label{eq:final95classical95ham} \widetilde{H}^{(c)} = \sum_{T,T'} \sum_{s' \in \{0,1\}^{|\mathcal{C}_T \cup \mathcal{C}_{T'}|}} \vert s'\rangle\!\langle s'\vert \otimes h^{(c,s')}_{T,T'} + \sum_T \sum_{s' \in \{0,1\}^{|\mathcal{C}_T |}} \vert s'\rangle\!\langle s'\vert \otimes h^{(c,s')}_{T}.\tag{22}\] Since each tile \(T\) is of constant size, each set \(\mathcal{C}_T \cup \mathcal{C}_{T'}\) has constant size as well. This implies we only get a constant blow-up in the number of terms of \(H\), yielding a polynomial-size instance. The locality of the instance depends on the max size of the set \(\mathcal{C}_T \cup \mathcal{C}_{T'}\). Assuming this is bounded by \(k\), we obtain an \(2 + k\)-local classical Hamiltonian.

Correcting Removed Terms \(\mathcal{R}\). As in the original proof, a ground state \(\left\vert \psi \right\rangle\) for \(\widetilde{H}^{(c)}\) can be constructed in NP. To obtain a ground state for the original Hamiltonian \(H\), we need to correct the terms removed terms \(\mathcal{R}\). The first step is to remove the classical qubits. We do this by successively measuring the classical qubits to obtain a terminating restriction \(\boldsymbol{s} = (s_0, \dots, s_\ell)\).

Therefore, the above procedure yields a Hamiltonian \(H' = \widetilde{H}|_{\boldsymbol{s}}\) and a ground state \(\left\vert \psi \right\rangle\) of \(H'\). Now, recall that the removed term \(h \in \mathcal{R}\) is an interior, fully quantum term of the Hamiltonian \(H\). 1 implies that there is a correction unitary operator \(L_h\) anti-commuting with \(h\) and commuting with every other term in \(H'\). Thus, we may perform the measurement \(\mathcal{M} = \{\tfrac 1 2(\mathcal{I}+ h), \tfrac 1 2 (\mathcal{I}-h)\}\) and we can apply \(L_h\) to flip the measurement outcome if required.

5.2.3 Gibbs state sampling↩︎

In this section, we extend the previous proof so that rather than obtaining a single ground state of \(H\), we obtain a sample from the Gibbs distribution. Initially, we apply exactly the same steps as in the previous section, until we obtain a classical Hamiltonian \[\widetilde{H}^{(c)} = \sum_{T,T'} \sum_{s' \in \{0,1\}^{|\mathcal{C}_T \cup \mathcal{C}_{T'}|}} \vert s'\rangle\!\langle s'\vert \otimes h^{(c,s')}_{T,T'} + \sum_T \sum_{s' \in \{0,1\}^{|\mathcal{C}_T |}} \vert s'\rangle\!\langle s'\vert \otimes h^{(c,s')}_{T}\] equivalent to the original Hamiltonian with terms \(\mathcal{R}\) removed. Now, rather than obtaining a ground state of \(\widetilde{H}^{(c)}\), we assume that we can sample from the Gibbs distribution of \(\widetilde{H}^{(c)}\), and then argue that we can recover the Gibbs distribution of the original Hamiltonian, by correcting for the terms \(\mathcal{R}\). We will show how to correct these terms inductively. Let \(Q\) be the set of terms that were either not removed or already corrected, and let \(H^Q = \sum_{h\in Q} h\). Assume we have the Gibbs state \(\rho(Q)\) of \(H^Q\) \[\rho(Q) := \frac{1}{Z} \sum_{s, \lambda^Q_s} e^{-\beta \lambda^Q_s} \vert s\rangle\!\langle s\vert \otimes \Pi^Q_{\lambda_s},\] where we use the fact that the Hamiltonian \(H^Q\) is diagonal w.r.t. \(s \in \{0,1\}^{|\mathcal{C}_0|}\). To correct terms \(h \in \mathcal{R}\), we use the same observation from 21 that given some classical restriction \(s\), the resulting Hamiltonian \(H|_s\) is block-diagonal w.r.t. the qubits in \(\mathcal{C}_1(s)\). Thus, we may assume we see the state \[\rho^{(s,t)}(Q) = \sum_{\lambda^Q_{s,t}} e^{-\beta \lambda^Q_{s,t}} \vert s,t\rangle\!\langle s,t\vert\Pi^Q_{\lambda_{s,t}}\] with probability proportional to \(\sum_{\lambda^Q_{s,t}} e^{-\beta \lambda^Q_{s,t}}\). Here, \(\Pi^Q_{\lambda_{s,t}}\) is the projector onto the \(\lambda^Q_{s,t}\)-eigenspace of \(H^Q|_{s,t}\). The idea is to apply the proof of 13 within the subspace \(\vert s,t\rangle\!\langle s,t\vert \otimes \mathcal{I}\). Specifically, we make the following observations.

• 10 holds since each \(h \in \mathcal{R}\) is a fully quantum term (for any classical restriction, \(h\) is in the interior).

• By 1 there is a correction operator \(L_h\) anti-commuting with \(h\) and commuting with all other terms of \(H^Q|_{s,t}\), for any choice of \(s,t\). This assumption means that \(L_h\) itself is only defined over non-classical qubits and thus for a particular choice of \(s,t\), we can think of this operator as being \(\vert s,t\rangle\!\langle s,t\vert \otimes L_h\). Then, 11 follows by tensoring both sides of [eq:sym1,eq:sym2] with \(\vert s,t\rangle\!\langle s,t\vert\).

• Similarly, 12 holds, again by tensoring both matrices with \(\vert s,t\rangle\!\langle s,t\vert\).

These observations plus the proof of 13 implies that we can prepare the state, \[\hat{\rho}(Q \cup \{h\}) = \frac{1}{Z(Q \cup \{h\})}\sum_{\lambda^Q_{s,t}} \sum_{\lambda_h \in \{\pm c_p\}} \exp(-\beta \lambda^Q_{s,t}) \cdot \exp(-\beta \lambda_h) \vert s,t\rangle\!\langle s,t\vert\Pi^{Q \cup \{h\}}_{\lambda^Q_{s,t} + \lambda_h}\] It’s easy to see that within the subspace \(\vert s,t\rangle\!\langle s,t\vert\) this is the correct Gibbs state. Moreover, across different \((s,t), (s',t')\) pairs, we retain the correct distribution. This is because correcting additional fully quantum terms \(h\) only applies an identical multiplicative factor across each of the \(\vert s,t\rangle\!\langle s,t\vert\)-subspaces.

6 Acknowledgement↩︎

We thank Dominik Hangleiter, Sandy Irani, Anurag Anshu, and Yongtao Zhan for the helpful discussion.

YH is supported by the National Science Foundation Graduate Research Fellowship under Grant No. 2140743; YH is also supported by NSF grant CCF-2430375. Any opinion, findings, and conclusions or recommendations expressed in this material are those of the authors(s) and do not necessarily reflect the views of the National Science Foundation. Jiaqing Jiang is supported by MURI Grant FA9550-18-1-0161 and the IQIM, an NSF Physics Frontiers Center (NSF Grant PHY-1125565).

7 Gibbs sampling reduction for punctured defected Toric code↩︎

As described in the proof overview section (Section 1.3), our Gibbs sampler gives an \(\mathcal{O}(n^2)\)-time algorithm for preparing the Gibbs states of the punctured defected Toric code \(H_{DT}\). In this section we instead describe a different and slower Gibbs sampler for \(H_{DT}\), which achieves the following Claim 16 and captures the key ideas for our Gibbs sampler for general qubit 2D CLHs.

We recall the defected Toric code, \(H_{DT}\). As shown in 11, imagine partitioning the plaquettes in the 2D lattice as Black \(\mathcal{B}\) and White \(\mathcal{W}\). The defected Toric code is putting \(Z\) and \(X\) terms in white and black plaquetes respectively: \[\begin{align} &H_{DT} =\sum_{p\in \mathcal{B}} c_pX^p + \sum_{p\in \mathcal{W}} c_pZ^p, \end{align}\] where \(c_p\) can be an arbitrary real number. In the standard Toric code \(c_p=-1\) for any \(p\). One can check that \(H_{DT}\) is a qubit 4-local CLH on 2D.

In this section we describe the algorithm and the reduction in details, while as the correctness proof is omitted, as is the same as in Section 4. The algorithm will require removing a set \(\mathcal{R} \subseteq \mathcal{W}\) of terms from \(H_{DT}\). In particular, if we remove alternating rows of white terms (as in [fig:fig2b]) then treat the \(4\) qubits in the support of each remaining white term as a single grouped qudit, then the resulting Hamiltonian \(H^{(2)}_{DT}\) is \(2\)-local. We’ll denote the remaining white terms as \(\mathcal{O} = \mathcal{W} \setminus \mathcal{R}\).

7.0.0.1 Notation.

For simplicity, here we assume the 2D lattice is a square \(L\times L\) lattice embedded on a plain which has boundaries. The number of qubits is \(n=L\times L\). For better illustration, here we denote the computational basis as \(\left\vert x \right\rangle\in\{\pm 1\}^n\) rather than the conventional notation \(\left\vert x \right\rangle\in \{0,1\}^n\), where the Pauli \(X\) and \(Z\) operators act as \[\begin{align} &Z \left\vert 1 \right\rangle=\left\vert 1 \right\rangle, Z \left\vert -1 \right\rangle=-\left\vert -1 \right\rangle\\ & X \left\vert 1 \right\rangle = \left\vert -1 \right\rangle, X \left\vert -1 \right\rangle = \left\vert 1 \right\rangle. \end{align}\]

This notion is only used in this section. We say \(\left\vert x \right\rangle\) is of even Hamming weight if there are even 1s in \(x\).

7.0.0.2 The classical Hamiltonian \(H_{DT}^{(2c)}\).

To map \(H^{(2)}_{DT}\) to the classical Hamiltonian \(H_{DT}^{(2c)}\), we change the basis of the qudit. More specifically, for any term \(p\in \mathcal{O}\), consider the 4 qubits on \(p\), as shown in [fig:fig2b], name them as \(u,v,w,\tau\) . A natural basis for the 4 qubit Hilbert space is the computational basis, which are common eigenvalues of \(\{Z_u\}_{u\in p}\). An alternative labeling is choosing another set of 4 independent stabilizers \[\begin{align} S^p_1:=Z^p, S^p_2:=X_u\otimes X_v, S^p_3:=X_v\otimes X_w, S^p_4 := X_w\otimes X_\tau, \end{align}\] The 4 new stabilizers will specify a new basis for the 4-qubit Hilbert space indexed by \(s\in\{\pm 1\}^4\), which can be viewed as the computational basis for 4 virtual qubits. The 4 new stabilizers act as Pauli Z on the 4 virtual qubits. That is if we use \((Z_{*i})^p\) to denote Pauli Z on virtual qubit \(i\) in plaquette \(p\), then \[(Z_{*i})^p = S^p_i.\] Each \(s\in\{\pm 1\}^4\) indicates the common eigenvector of \(S^p_i\) w.r.t eigenvalues \(s_i\). Denote the corresponding vector as \(\left\vert s^p \right\rangle^*\) where \(*\) means \(\left\vert s^p \right\rangle^*\) is the computational basis for the virtual qubits rather than the original qubits. As an example one can verify that \(\left\vert \psi^p \right\rangle\) in Eq. (24 ) corresponds to \(\left\vert 1111^p \right\rangle^*\). One can check that in this virtual qubit basis, we have that the Hermitian terms¹² can be re-phrase as \[\begin{align} &\forall p\in \mathcal{O}, p = c_p\cdot (Z_{*1})^p, \end{align}\] Besides, note that \(X_{u}\otimes X_{\tau}= S_2^p\times S_3^p\times S_4^p\), and \(X^{\otimes 4} = X^{\otimes 2}\otimes X^{\otimes 2}\), we have \[\begin{align} &\text{For horizontal p connecting p_1,p_3, }p = c_p \cdot (Z_{*3})^{p_1} \otimes \left(Z_{*2}\otimes Z_{*3} \otimes Z_{*4} \right)^{p_3}.\\ &\text{For vertical p connecting p_1,p_2, }p = c_p \cdot (Z_{*4})^{p_1} \otimes (Z_{*2})^{p_2}. \end{align}\] Thus the resulting Hamiltonian is exactly \(H_{DT}^{(2c)}\) as shown in [fig:fig2c] and Eq. (23 ), where we omitting the subscript * for simplicity, \[\begin{align} H_{DT}^{(2c)} : =\sum_{\text{horizontal (p,p_i,p_j)}} c_p \cdot (Z_3)^{p_i}\otimes (Z_2Z_3Z_4)^{p_j} + \sum_{\text{vertical (p,p_i,p_j)}} c_p \cdot (Z_4)^{p_i}\otimes (Z_2)^{p_j} +\sum_{p\in O} c_p Z_1^p \label{eq:classical} \end{align}\tag{23}\] In summary, by choosing another basis each qudit can be viewed as four virtual qubits, and the 2-local Hamiltonian \(H_{DT}^{(2)}\) becomes a classical Hamiltonian \(H_{DT}^{(2c)}\) w.r.t those virtual qubits. The computational basis of the virtual qubits \(\left\vert \boldsymbol{y} \right\rangle^*\) corresponds to an eigenvector of \(H_{DT}^{(2)}\), denoted as \(\left\vert \phi(\boldsymbol{y}) \right\rangle\). Denote the eigenvalue w.r.t \(\left\vert \boldsymbol{y} \right\rangle^*\) and \(p\in \mathcal{B}\cup\mathcal{O}\) as \(\lambda_p(\boldsymbol{y})\), and define \[\lambda_{\mathcal{B}\cup\mathcal{O}}(\boldsymbol{y}) :=\sum_{p\in \mathcal{B}\cup\mathcal{O}} \lambda_p(\boldsymbol{y}).\]

7.1 Preparation of a ground state↩︎

Before describing our Gibbs state sampler, we first review a folklore quantum algorithm for preparing the ground state of the standard Toric code.

For ease of presentation, for any subset \(Q\subseteq \mathcal{W}\) or \(Q\subseteq\mathcal{B}\), we use “terms in \(Q\)" to denote the operators \((-X^p)\) if \(p\in \mathcal{B}\), and \((-Z^p)\) if \(p\in \mathcal{W}\). Note that the ground state of the standard Toric code is the common \((-1)\)-eigenvector of all terms in \(\mathcal{W}\) and \(\mathcal{B}\). The first step in preparing the ground state of \(H_{DT}\) is preparing the ground state of \(H^{(2)}_{DT}\), where \[H_{DT}^{(2)}:=\sum_{p\in \mathcal{B}\cup\mathcal{O}} p\] is two local.

The ground state of \(H_{DT}^{(2)}\) is easy to describe. For each term \(p\in \mathcal{O}\), denote \(\left\vert \psi^p \right\rangle\) as the uniform superposition of basis states with even Hamming weights, that is \[\begin{align} \left\vert \psi^p \right\rangle =\frac{1}{2}\sum_{x\in \{\pm 1\}^p, |x| \text{ even}} \left\vert x^p \right\rangle,\label{eq:psip} \end{align}\tag{24}\] where \(\{\pm 1\}^p\) is the computational basis of the four qubits of \(p\). One can check that \[\begin{align} \left\vert \psi^{\mathcal{B}\cup \mathcal{O}} \right\rangle:= \otimes_{p\in \mathcal{O}} \left\vert \psi^p \right\rangle, \end{align}\] is the common \((-1)\)-eigenstate of all terms in \(\mathcal{B}\cup\mathcal{O}\).

Next, we need to correct for the removed terms \(p \in \mathcal{R}\). Starting with \(\left\vert \psi^{\mathcal{B}\cup\mathcal{O}} \right\rangle\), we sequentially measure the current state w.r.t the measurement \(-Z^p\), for each \(p \in \mathcal{R}\).

• If we get \(- 1\), we know the current state is the common \((-1)\)-eigenstate for all terms in \(\mathcal{B}\cup\mathcal{O}\cup\{p\}\) and we move to the next \(p\in\mathcal{R}\).

• Otherwise we perform a deterministic correction: as shown in [fig:intro95fig95correct], we can connect term \(p\) to the boundary by a path \(\gamma_p\). The deterministic correction is done by applying a sequence of \(X\) operator, denoted as \[\begin{align} L_{p}=\otimes_{v\in \gamma} X_v. \label{eq:correction} \end{align}\tag{25}\]

where \(X_v\) is the Pauli X on qubit \(v\). Note that \(L_p\) anti-commutes with \(p\), and commutes with all other terms. Thus the final state will again be the common \((-1)\)-eigenstate for all terms in \(\mathcal{B}\cup\mathcal{O}\cup\{p\}\).

In summary, to prepare the ground state, one first prepares the ground state of a 2-local Hamiltonian \(\sum_{p\in \mathcal{B}} -X^p + \sum_{p\in \mathcal{O}} -Z^p.\) Then performs a deterministic correction to modify the ground state of the 2-local Hamiltonian to the ground state of the Toric code. The ground state preparation for the defected Toric code is similar.

7.2 Preparation of the Gibbs state↩︎

We now describe how to generalize the ground state preparation algorithm to prepare Gibbs state for the defected Toric code. At a high level, the idea is to:

• Step 1: First prepare the Gibbs state of the Hamiltonian \[H_{DT}^{(2)}:=\sum_{p\in \mathcal{B}\cup\mathcal{O}} p.\] by mapping it to the classical 2-local Hamiltonian \(H_{DT}^{(2c)}\) defined in Eq. (23 ).

• Step 2: Perform a randomized correction, to modify the current Gibbs state for \(H_{DT}^{(2)}\) to the final Gibbs state for \(H_{DT}\).

7.2.0.1 Step 1: Gibbs state for \(H_{DT}^{(2)}\).

As in [fig:fig2b], \(H^{(2)}_{DT}\) is 2-local in the sense that: if we group every 4 qubits in a plaquette \(p\in \mathcal{O}\) as one qudit, then every \(p\in \mathcal{O}\) only acts on one qudit, every \(p\in \mathcal{B}\) acts on two qudits.

Moreover, since we only change local basis for every qudit, \(\left\vert \phi(\boldsymbol{y}) \right\rangle\) is in fact a tensor product of single-qudit state. Thus we can prepare the quantum Gibbs state of \(H_{DT}^{(2)}\) by a simple algorithm: do classical Gibbs sampling for \(H_{DT}^{(2c)}\), get a string \(\boldsymbol{y}\), and prepare the tensor product state \(\left\vert \phi(\boldsymbol{y}) \right\rangle\).

7.2.0.2 Step 2: Correcting removed terms.

Step 2 is more tricky. By Step 1, we assume that we can sample the string \(\boldsymbol{y}\) with probability \(\exp(-\beta \lambda_{\mathcal{B}\cup\mathcal{O}}(\boldsymbol{y}))/Z_{\mathcal{B}\cup\mathcal{O}}\), where \(Z_{\mathcal{B}\cup\mathcal{O}}\) is the normalization factor.

Then we take a term in the removed set \(p\in \mathcal{R}\), and try to prepare the Gibbs state w.r.t. \(\left(H_{\mathcal{B}\cup\mathcal{O}\cup\{p\}},\beta\right)\). Consider measuring \(\left\vert \phi(\boldsymbol{y}) \right\rangle\) w.r.t \(p\). Denote the measurement outcome as \(\lambda \in\{\pm c_p\}\). Denote the projectors to be \[\begin{align} \Pi^{p}_{+c_p} : = \frac{1}{2}\left(I+\frac{p}{c_p}\right), \quad \Pi^{p}_{-c_p} : = \frac{1}{2} \left(I-\frac{p}{c_p}\right). \end{align}\]

We have \[\begin{align} Pr\left( \text{outcome is } \lambda \right) = \exp(-\beta \lambda_{\mathcal{B}\cup\mathcal{O}}(\boldsymbol{y}))/Z \cdot \langle \phi(\boldsymbol{y})| \Pi^{p}_{\lambda} | \phi(\boldsymbol{y})\rangle. \label{eq:prac} \end{align}\tag{26}\] Recall that at this moment, the ideal distribution we want is \[\begin{align} Pr_{ideal} \left( \text{outcome is }\lambda\right) \text{ is proportional to } \exp(-\beta \lambda_{\mathcal{B}\cup\mathcal{O}}(\boldsymbol{y})) \cdot \exp(-\beta \lambda).\label{eq:ideal} \end{align}\tag{27}\] Compared to the preparation of ground state, which only needs to correct the eigenvalue, in the task of Gibbs state preparation, to get Eq. (27 ) from Eq. (26 ), it seems that one needs to correct

• The probability incurred by measurement, that is \(\langle \phi(\boldsymbol{y})| \Pi^{p}_{\lambda} | \phi(\boldsymbol{y})\rangle\).

• The probability incurred by the new energy \(\exp(-\beta \lambda)\).

This first correction is tricky because it depends on \(\phi(\boldsymbol{y})\), thus the probability to be corrected is different for every \(\boldsymbol{y}\). We circumvent this problem by an intuitive oblivious randomized correction. That is if we get outcome \(\lambda\) then

• With probability \(prob:=\frac{\exp(-\beta \lambda)}{\exp(\beta \lambda)+\exp(-\beta \lambda)}\) we do nothing.

• With probability \(1-prob\) we apply the correction operation \(L_{p}\) defined in Eq. (25 ).

One may doubt whether the above algorithm successfully prepares the Gibbs state or not, since the correction only brings us back to the right eigenspace, but does not really help us get back to the right state, as noted in 18. That is \[\begin{align} & L_{p} \Pi^p_{+c_p} \left\vert \phi(\boldsymbol{y}) \right\rangle \not\propto \Pi^p_{-c_p} \left\vert \phi(\boldsymbol{y}) \right\rangle, \end{align}\] where here \(\not\propto\) means the two vectors are not proportional to each other. In fact, \(L_p \Pi_{+c_p} \left\vert \phi(\boldsymbol{y}) \right\rangle\) might not even be (proportional to) one of the eigenstates \(\{ \Pi^p_{+c_p}\left\vert \phi(\boldsymbol{y}) \right\rangle\}_{\boldsymbol{y}}\cup \{ \Pi^p_{-c_p}\left\vert \phi(\boldsymbol{y}) \right\rangle\}_{\boldsymbol{y}}\). The key fact that makes this oblivious randomized correction work, is the symmetry of the eigenspace and the fact that \(L_p\) is a unitary. More specifically, let \(\Pi^{\mathcal{B}\cup\mathcal{O}}_{\boldsymbol{y}}\) be the common eigenspace of \(p'\in \mathcal{B}\cup\mathcal{O}\) w.r.t eigenvalue \(\lambda_{p'}(\boldsymbol{y})\). Then \[\begin{align} L_p \Pi^{p}_{+c_p} \Pi^{\mathcal{B}\cup\mathcal{O}}_{\boldsymbol{y}} \Pi^{p}_{+c_p} L_p = \Pi^{p}_{-c_p} \Pi^{\mathcal{B}\cup\mathcal{O}}_{\boldsymbol{y}}\Pi^{p}_{-c_p}. \end{align}\] Finally, note that all \(p\in \mathcal{R}\) can be corrected independently. Thus it suffices to perform the same measure and randomized correction procedure sequentially.

For general qubit 2D CLH \(H\) without classical qubits, denote the removable terms \(\mathcal{R}\) as the set of terms whose measurement value can be corrected without changing the value of other terms. We prepare the Gibbs state of \(H\) in a similar way as for the defected Toric code. That is we first remove terms in \(\mathcal{R}\) and transform the remaining 2-local Hamiltonian \(H^{(2)}\) to a classical Hamiltonian \(H^{(2c)}\). Then we first do classical Gibbs sampling for \(H^{(2c)}\), then prepare the Gibbs state of \(H^{(2c)}\), then perform measurement and randomized correction to prepare the Gibbs state of \(H\).

8 Gibbs Sampler for defected Toric Code embedded on a closed 2D surface↩︎

As described in 1.3.1, the oblivious correction algorithm works only when there are punctured terms that act as a boundary in the defected Toric code. If the defected Toric code is defined on a torus—as in the standard formulation in the error-correcting code literature—then correction operators correspond to strings of Pauli operators connecting two measured terms of the Hamiltonian, and the oblivious correction approach does not applies directly. Nonetheless, we show that there is a simple algorithm which prepares the Gibbs state of such defected Toric code on torus. One can similarly prepare the Gibbs state for the defected Toric code embedded on other closed 2D surface, such as a sphere, or for the models defined in [40].

As before, we partition the terms into “white” terms \(\mathcal{W}\) and “black” terms \(\mathcal{B}\), which correspond to \(Z\) and \(X\)-type terms respectively.

The algorithm proceeds in two phases. In the first phase, we prepare the Gibbs state for the \(Z\)-type terms, i.e., \[\label{eq:z95gibbs}\rho(H_{\mathcal{W}}, \beta) = \sum_{\lambda \in \Lambda_{\mathcal{W}}} \mu_{\mathcal{W}}(\lambda) \Pi_{\lambda}.\tag{28}\] Here \(H_{\mathcal{W}} = \sum_{p \in \mathcal{W}} p\), and \(\Lambda_{\mathcal{W}}\) denote the set of length-\(|\mathcal{W}|\) vectors with entries \(\pm 1\) whose product equals 1, that is, \(\lambda=(...,\lambda_p,...)_{p\in W}\) where \(\lambda_p\in\{+1,-1\}\) and \(\prod_{p\in \mathcal{W}} \lambda_p =1\). The vector \(\lambda\) is used to label the eigenspaces of \(H_{\mathcal{W}}\), denoted as \(\Pi_\lambda\), which satisfies \(\forall p\in \mathcal{W}, p \Pi_\lambda =\lambda_p \Pi_\lambda.\) The \(\mu_{\mathcal{W}}(\lambda)\) is the corresponding Boltzmann weights, \[\mu_{\mathcal{W}}(\lambda)= \prod_{p\in \mathcal{W}} exp(-\beta \lambda_p c_p)/Z_{\mathcal{W}} \text{ where Z_{\mathcal{W}}= \sum_{\lambda \in \Lambda_{\mathcal{W}}} \prod_{p\in \mathcal{W}} exp(-\beta \lambda_p c_p)}.\] In the second phase, begin with \(\rho(H_{\mathcal{W}}, \beta)\) we further consider the \(X\)-type terms, prepares the full states for \(H=H_{\mathcal{W}}+H_{\mathcal{B}}\), that is \[\rho(H, \beta) = \sum_{\lambda \in \Lambda_{\mathcal{W}}} \mu_{\mathcal{W}}(\lambda) \Pi_{\lambda} \sum_{\lambda'\in \Lambda_{\mathcal{B}}} \mu_{\mathcal{W}}(\lambda') \Pi_{\lambda'}\,, \label{eq:Hbeta}\tag{29}\] Note that 29 holds since the correction operators for \(X\) terms commute with the correction operators for \(Z\) terms, thus the correction operators commute with each \(\Pi_\lambda\).

More specifically, for both phases, we use 1D Ising chains to guide the correction operations. For the first phase, we construct the clasiscal Ising chain \[\begin{align} H_\text{Ising} = \sum_{i=1}^k c_i h_i,\label{eq:Toric95Ising} \end{align}\tag{30}\] with \(h_i = Z_i Z_{i+1}\), \(k = |\mathcal{W}|\), \(k+1\) identified with \(1\), and the coefficients \(c_i\) will be specified later. Start with the uniform distribution over bit strings. In each step of Glauber dynamics for a state \(\left\vert x \right\rangle\), a site \(i \in [k]\) is picked randomly, and the change in energy \(\Delta E = \text{Tr}\left[H_{\mathrm{Ising}} (\vert x\rangle\!\langle x\vert - X_i \vert x\rangle\!\langle x\vert X_i)\right]\) is computed. The bit flip \(X_i\) is accepted with probability \(\min\{1,exp(-\beta \Delta E)\}\). To map this Gibbs sampling for Ising model onto Gibbs states preparation for \(H_{\mathcal{W}}\), we order the terms \(p \in \mathcal{W}\) in any arbitary order \(p_1, \dots, p_i, \cdots, p_{k}\), with \(k = |\mathcal{W}|\). For convenience, we also write the corresponding coefficient \(c_p\) in Eq. (1 ) as \(c_i\). Note that each \(p_i\in \mathcal{W}\) in the Toric code will correspond to \(h_i\) in the aforementioned Ising model, i.e. Eq. (30 ). Imagine the strings \(x \in \{0,1\}^k\) as labeling the eigenvalues of \((p_1, \dots, p_k)\) via \[\label{eq:eig95mapping} \lambda_i = (-1)^{x_i \oplus x_{i+1}}\,.\tag{31}\] We argue that this mapping is valid. \(H_{\mathcal{W}}\)’s eigenspaces are labeled by \(\lambda \in \{\pm 1\}^k\) with the parity constraint \(\prod_{i \in [k]} \lambda_i = 1\). This constraint is obeyed by 31 : \[\begin{align} \prod_{i=1}^k \lambda_i = (-1)^{x_1 \oplus x_2} \cdot (-1)^{x_2 \oplus x_3} \dots (-1)^{x_{k-1} \oplus x_k} \cdot (-1)^{x_k \oplus x_1} = (-1)^{x_1 \oplus x_k} \cdot (-1)^{x_k \oplus x_1} = 1\,.\label{eq:TC95prod1} \end{align}\tag{32}\] Moreover, the update rule of Glauber dynamics maps cleanly to \(H_{\mathcal{W}}\). The classical bit flip \(X_i\) is mapped to the string operator \(X^s\), where \(s\) is a path from \(p_i\) to \(p_{i+1}\) in the 2D Toric code lattice. This operator is defined so that it anti-commutes with \(p_{i-1}\) and \(p_{i}\) and commutes with all other terms \(p \in \mathcal{W}\), and commutes with all terms in \(\mathcal{B}\). This retains the equivalency in 31 since applying \(X_i\) to the \(i\)-th bit toggles the parity, and the operation takes \(\lambda_i \rightarrow -\lambda_i\) and \(\lambda_{i-1} \rightarrow -\lambda_{i-1}\). This is exactly the effects of applying the correction operators \(X^s\) on the Toric code. One can see from 31 that the Gibbs distribution of \(H_\text{Ising}\) corresponds to 28 since \[\begin{align} \rho(H_\text{Ising}, \beta) &= \frac{1}{Z} \sum_{x \in \{0,1\}^k} \exp(-\beta \text{Tr}[H_\text{Ising} \vert x\rangle\!\langle x\vert]) \vert x\rangle\!\langle x\vert\\ &= \frac{1}{Z} \sum_{x \in \{0,1\}^k} \exp\left(-\beta \sum_i c_i (-1)^{x_i \oplus x_{i+1}}\right) \vert x\rangle\!\langle x\vert\\ &= \frac{1}{Z} \sum_{\lambda \in \{\pm 1\}^k} \exp\left(-\beta c_i \lambda_i \right) \vert x\rangle\!\langle x\vert\\ \end{align}\] where each \(x\) label the eigenspace projector \(\Pi_\lambda\) in \(\rho(H_{\mathcal{W}},\beta)\).

Finally, to prepare the full Gibbs state for \(H = H_{\mathcal{W}} + H_{\mathcal{B}}\), we begin with the state \(\rho(H_{\mathcal{W}}, \beta)\), and repeat the above mapping—this time with respect to the \(X\)-type terms instead of the \(Z\)-type terms. Specifically, starting from \(\rho(H_{\mathcal{W}}, \beta)\), we run a classical Glauber dynamics corresponding to the terms in \(\mathcal{B}\), and use this Ising dynamics to guide the correction procedure for the \(X\)-type terms on the state \(\rho(H_{\mathcal{W}}, \beta)\).

Note that Gibbs sampling for 1D Ising model is known to be rapid mixing and the corresponding Gibbs samplers have runtime \(O(n \,poly(\log n))\)  [21]–[23], [45]. Thus the runtime of our algorithm is \(O(npoly(log n) \times n)=\mathcal{O}(n^ 2poly(\log n))\), where \(n\) is the cost for each correction operation.

Finally, note that the algorithm described above for preparing the Gibbs state of the defected Toric code on a torus can be directly generalized to defected Toric code defined on any closed 2D surfaces. This is because the algorithm does not rely on any special properties of the torus, but only on the constraint that all syndromes of the same type multiply to one, and on the fact that in the defected toric code the value of a syndrome can be flipped by string operators.

9 Proofs of reductions for specific Hamiltonians↩︎

The proof of Lemma 7 involves the notations of induced algebra, which should be read after reading section 5.1.1.

Author Contribution↩︎

J.J. conceived the study. J.J and Y.H. performed the theoretical analysis. All authors discussed the results and contributed to writing the manuscript.

References↩︎

---

1. yeongwoohwang@g.harvard.edu↩︎

2. jiaqingjiang95@gmail.com↩︎

3. In our case, we will obtain scaling like \(T+\mathcal{O}(n)\) or \(T+\mathcal{O}(n^2)\) or \(T\times O(n)\).↩︎

4. Lemma 12 in [42], where the k-QMV can be estimated by measuring the Gibbs state .↩︎

5. More precisely there is an equivalence between the mixing time and the spatial decay of correlation in the Gibbs measure [49]–[52]. In some classical literature, decay of correlation is referred to as mixing condition [21] or spatial mixing [49].↩︎

6. We did not check that whether \(\beta_1,\beta_2\) and the later mentioned high temperatures are equal.↩︎

7. Intuitively, this is the set of terms which remains quantum under any recursive restriction of the classical qubits, see 14.↩︎

8. The formal definition of classical qubits is given in Definition 3.↩︎

9. The definition of a “2D qubit CLH” in [40] is slightly different; qubits are put on the edges of 2D lattice while we put qubits on the vertices. However, the two settings are equivalent, as explained in Appendix C of [39].↩︎

10. That is, the defected Toric code \(H=\sum_{p\in \mathcal{W}} c_p Z^p + \sum_{p\in \mathcal{B}} c_p X^p\) is embedded on a closed 2D surface and the coefficient \(c_p\neq 0\), for all \(p\).↩︎

11. Recall that in [40] classical qubits have already been removed so \(h\) and \(h'\) are automatically quantum.↩︎

12. With some abuse of notations, we use \(p\) to denote both the plaquette and the Hermitian term on the plaquette.↩︎

13. Note that there are constructive proofs for the structure Lemma, as in Section 7.3 of [56].↩︎