Polynomial-Time Preparation of Low-Temperature Gibbs States for 2D Toric Code

The ability (or the lack of it) to efficiently prepare Gibbs states has far-reaching implications in quantum information theory, condensed matter physics, quantum chemistry, statistical mechanics, and optimization. Given a quantum Hamiltonian $H\in\mathbb{C}^{2^{N}\times 2^{N}}$, we would like to prepare the associated thermal state $\sigma_{\beta}\propto e^{-\beta H}$, where $\beta$ is the inverse temperature. A number of quantum algorithms have been designed to efficiently prepare high-temperature Gibbs states with a small $\beta=\mathcal{O}(\mathrm{poly}(N^{-1})$) [44, 18, 52, 26, 1]. However, as the temperature lowers (i.e., $\beta$ becomes large), the complexity of these algorithms can scale exponentially in the number of qubits $N$, rendering them impractical for low-temperature regimes where the Gibbs state has a significant overlap with the ground state of $H$.

Recent advancements have rekindled interest in designing quantum Gibbs samplers based on Lindblad dynamics [41, 47, 14, 16, 17, 55, 21]. These algorithms rely on a specific form of open quantum system dynamics to drive the system toward its thermal equilibrium, an idea pioneered by Davies in the 1970s [19, 20]. The efficiency of such algorithms largely depends on the mixing time of the underlying dynamics, which can vary significantly across different systems and different forms of Lindbladians.

The computational complexity of preparing quantum Gibbs states, computing partition functions, and the potential for establishing quantum advantage in these tasks is a topic of ongoing debate in the literature. On one hand, at high enough temperatures, there exists polynomial-time classical algorithms to sample from Gibbs states and to estimate partition functions [7, 40, 56, 10]. On the other hand, in the low-temperature regime, preparing classical Gibbs states is already \NP-hard in the worst case [4, 48], and we do not expect efficient quantum algorithms in these cases. The development of these new Gibbs samplers has also contributed to advancements in our complexity-theoretic understanding [45, 8, 46]. [45] proved that simulating the Lindbladian proposed in [17] to time $T=\poly(N)$ at $\beta=\Omega(\log(N))$ for a $k$-local Hamiltonian is \BQP-complete. [8, 46] constructed a family of $k$-local Hamiltonians such that quantum Gibbs sampling at constant temperatures (lower than the classically simulatable threshold) can be efficiently achieved with the block-encoding framework [16], and the task is classically intractable assuming no collapse of the polynomial hierarchy.

None of these constructions imply efficient preparation of Gibbs states at low temperatures. Indeed, when the temperature is sufficiently low, the Gibbs state can exhibit a high overlap with the ground state, and cooling to these temperatures is expected to be \QMA-hard in the worst case. However, it is important to recognize that \QMA-hardness does not preclude the possibility of developing efficient Gibbs samplers for specific Hamiltonians. Consequently, understanding and controlling the mixing time for specific systems (or specific classes of systems) is a fundamental open question in this field and can provide valuable insights into the practical performance of these Gibbs samplers [2, 50, 49, 22, 11, 23, 51, 36, 31, 5, 30, 6, 45].

In this work, we propose a novel Gibbs sampler with nonlocal jump operators, and analyze its convergence rate for preparing low-temperature Gibbs states of the 2D toric code [33], a paradigmatic model in quantum information theory, quantum error correction, and condensed matter physics. In the context of thermalization, when the toric code is exposed to thermal noise modeled by a specific form of Lindbladians called the Davies generator, the seminal work of Alicki et al. [2] showed that the inverse spectral gap (which leads to an upper bound of the mixing time) grows exponentially with inverse temperature $\beta$. Using a different argument based on energy barriers, Temme et al. [50, 49] confirmed that the thermalization time (i.e., mixing time of the Lindblad dynamics) of the 2D toric code with local noise should indeed scale exponentially with $\beta$. However, it is unknown whether the $\exp(\beta)$ factor in the mixing time is unavoidable for all Lindblad dynamics on this problem. This leads to the central question of this work:

Can we design a quantum algorithm, based on Lindblad dynamics, that efficiently prepares the low-temperature Gibbs state of the 2D toric code with a polynomial runtime in both the inverse temperature $\beta$ and the number of qubits $N$?

Let $H$ be a quantum many-body Hamiltonian on $\mathcal{H}\cong\mathbb{C}^{2^{N}}$ and $\sigma_{\beta}=e^{-\beta H}/\mathcal{Z}_{\beta}$ be the associated thermal state, where $\beta$ is the inverse temperature and $\mathcal{Z}_{\beta}=\operatorname{Tr}(e^{-\beta H})$ is the partition function. In this section, we shall recall the canonical form of the Davies generator with $\sigma_{\beta}$ being the invariant state and some basic facts for its spectral gap analysis.

A superoperator $\Phi:\mathcal{B}(\mathcal{H})\to\mathcal{B}(\mathcal{H})$ is a quantum channel if it is completely positive and trace preserving (CPTP). Lindblad dynamics is a $C_{0}$-semigroup of quantum channels with the generator defined by $\mathcal{L}^{\dagger}(\rho):=\lim_{t\to 0^{+}}t^{-1}(\mathcal{P}^{\dagger}_{t}(\rho)-\rho)$ for $\rho\in\mathcal{D}(\mathcal{H})$. Here and in what follows, we adopt the convention that the adjoint operators (i.e., those with ${\dagger}$) are the maps in Schrödinger picture acting on quantum states. Both generators $\mathcal{L}$ and $\mathcal{L}^{\dagger}$ are usually referred to as Lindbladian. Davies generator is a special class of Lindbladians derived from the weak coupling limit of open quantum dynamics with a large thermal bath.

We first introduce the Bohr frequencies of $H$ by

$$B_{H}:=\{\omega=\lambda_{i}-\lambda_{j}\,:\ \forall\leavevmode\nobreak\ \lambda_{i},\lambda_{j}\in{\rm Spec}(H)\}\,,$$

where ${\rm Spec}(H)$ is the spectral set of $H$. Let $\{S_{a}\}_{a\in\mathcal{A}}$ be a set of coupling operators with $\mathcal{A}$ being a finite index set that satisfies

$$\{S_{a}\}_{a\in\mathcal{A}}=\{S_{a}^{\dagger}\}_{a\in\mathcal{A}}\,.$$

(2.1)

The jump operators $\{S_{a}(\omega)\}_{a,\omega}$ for the Davies generator are given by the Fourier components of the Heisenberg evolution of $S_{a}$:

$$e^{iHt}S_{a}e^{-iHt}=\sum_{\omega\in B_{H}}e^{i\omega t}\sum_{\lambda_{i}-\lambda_{j}=\omega}P_{i}S_{a}P_{j}:=\sum_{\omega\in B_{H}}e^{i\omega t}S_{a}(\omega)\,.$$

(2.2)

where $P_{i/j}$ is the projection into the eigenspace $\lambda_{i/j}$. By (2.1), we have $S_{a}(\omega)^{\dagger}=S_{a}(-\omega)$ for any $\omega\in B_{H}$.

We then introduce the Davies generator in the Heisenberg picture:

$$\mathcal{L}_{\beta}(X):=\sum_{a\in\mathcal{A}}\mathcal{L}_{S_{a}}(X)\,,\quad X\in\mathcal{B}(\mathcal{H})\,,$$

(2.3)

with

$$\mathcal{L}_{S_{a}}(X):=\sum_{\omega\in B_{H}}\gamma_{a}(\omega)\left(S_{a}(\omega)^{\dagger}\,X\,S_{a}(\omega)-\frac{1}{2}\left\{S_{a}(\omega)^{\dagger}S_{a}(\omega),X\right\}\right)\,,$$

(2.4)

where the transition rate function $\gamma_{a}(\omega)>0$ is given by the Fourier transform of the bath auto-correlation function satisfying the KMS condition [32]:

$$\gamma_{a}(-\omega)=e^{\beta\omega}\gamma_{a}(\omega)\,.$$

In this work, we always choose the transition rate function $\gamma_{a}(\omega)$ as the Glauber form:

$$\gamma_{a}(\omega)=\frac{2}{e^{\beta\omega}+1}\,,$$

(2.5)

For later use, we define $g_{a}(\omega):=e^{\beta\omega/2}\gamma_{a}(\omega)$ and find $g_{a}(\omega)=g_{a}(-\omega)$. Then, letting

$$L_{a}(\omega)=\sqrt{g_{a}(\omega)}S_{a}(\omega)\,,$$

(2.6)

we reformulate the Davies generator (2.3)–(2.4) as follows:

	$\displaystyle\mathcal{L}_{\beta}(X)$	$\displaystyle=\sum_{a\in\mathcal{A}}\sum_{\omega\in B_{H}}e^{-\beta\omega/2}\Big{(}L_{a}(\omega)^{\dagger}XL_{a}(\omega)-\frac{1}{2}\left\{L_{a}(\omega)^{\dagger}L_{a}(\omega),X\right\}\Big{)}$		(2.7)
		$\displaystyle=\frac{1}{2}\sum_{a\in\mathcal{A}}\sum_{\omega\in B_{H}}e^{-\beta\omega/2}L_{a}(\omega)^{\dagger}[X,L_{a}(\omega)]+e^{\beta\omega/2}[L_{a}(\omega),X]L_{a}(\omega)^{\dagger}$
		$\displaystyle=\frac{1}{2}\sum_{a\in\mathcal{A}}\sum_{\omega\in B_{H}}\gamma_{a}(\omega)S_{a}(\omega)^{\dagger}[X,S_{a}(\omega)]+\gamma_{a}(-\omega)[S_{a}(\omega),X]S_{a}(\omega)^{\dagger}\,.$

where the second step follows from

$\displaystyle L_{a}(\omega)^{\dagger}=e^{\beta\omega/4}\sqrt{\gamma_{a}(\omega)}S_{a}(-\omega)=e^{\beta\omega/4}\sqrt{e^{-\beta\omega}\gamma_{a}(-\omega)}S_{a}(-\omega)=L_{a}(-\omega)\,.$

We now define the GNS inner product associated with the Gibbs state $\sigma_{\beta}$:

$$\langle Y,X\rangle_{\sigma_{\beta}}=\operatorname{Tr}(Y^{\dagger}X\sigma_{\beta})\,.$$

(2.8)

It is known [32, 21] that the Davies generator satisfies the GNS detailed balance:

$$\langle Y,\mathcal{L}_{\beta}(X)\rangle_{\sigma_{\beta}}=\langle\mathcal{L}_{\beta}(Y),X\rangle_{\sigma_{\beta}}\,,$$

and thus the associated Lindblad dynamics $e^{t\mathcal{L}_{\beta}^{\dagger}}$ admits $\sigma_{\beta}$ as an invariant state, i.e.,

$$\mathcal{L}_{\beta}^{\dagger}(\sigma_{\beta})=0\,.$$

It follows that $\mathcal{L}_{\beta}$ is similar to a self-adjoint operator for the HS inner product, called the master Hamiltonian, and has only real spectrum. To be precise, we introduce the transform $\varphi_{X}:=X{\sigma_{\beta}}^{1/2}$, which gives $\langle Y,X\rangle_{\sigma_{\beta}}=\langle\varphi_{Y},\varphi_{X}\rangle$. Then, the master Hamiltonian $\widetilde{\mathcal{L}}_{\beta}$ is given by the similar transform of $\mathcal{L}_{\beta}$ via $\varphi_{X}$:

$$\widetilde{\mathcal{L}}_{\beta}:=\varphi\circ\mathcal{L}_{\beta}\circ\varphi^{-1}\,,$$

(2.9)

satisfying $\langle Y,\mathcal{L}_{\beta}(X)\rangle_{\sigma_{\beta}}=\langle\varphi_{Y},\widetilde{\mathcal{L}}_{\beta}\varphi_{X}\rangle$. It is easy to see that $-\widetilde{\mathcal{L}}_{\beta}$ is positive semi-definite and $\sqrt{\sigma_{\beta}}$ is the zero-energy ground state of $-\widetilde{\mathcal{L}}_{\beta}$. Thus, the spectral gap of $\mathcal{L}_{\beta}$ is the same as the ground state spectral gap of the Hamiltonian $-\widetilde{\mathcal{L}}_{\beta}$.

We say that $\mathcal{L}_{\beta}$ is primitive if $\sigma_{\beta}$ is the unique invariant state; equivalently, the kernel $\ker(\mathcal{L}_{\beta})$ is of one dimension, spanned by ${\bf 1}$. In this case, we have

$\displaystyle\lim_{t\to\infty}e^{t\mathcal{L}_{\beta}^{\dagger}}(\rho)=\sigma_{\beta}\,,\quad\forall\rho\in\mathcal{D}(\mathcal{H})\,,$

and the spectral gap ${\rm Gap}(\mathcal{L}_{\beta})$ of the primitive Davies generator can be characterized by the variational form:

$${\rm Gap}(\mathcal{L}_{\beta})=\inf_{X\neq 0\,,\operatorname{Tr}(\sigma_{\beta}X)=0}\frac{\langle X,\mathcal{L}_{\beta}(X)\rangle_{\sigma_{\beta}}}{\langle X,X\rangle_{\sigma_{\beta}}}\,.$$

Moreover, thanks to the detailed balance, the operator norm of $\mathcal{L}_{\beta}$ can be computed by

$$\left\lVert\mathcal{L}_{\beta}\right\rVert_{\sigma_{\beta}\to\sigma_{\beta}}=\sup_{X\neq 0}\frac{\langle X,\mathcal{L}_{\beta}(X)\rangle_{\sigma_{\beta}}}{\langle X,X\rangle_{\sigma_{\beta}}}\,.$$

By [54, 57], a sufficient and necessary condition for the primitivity of $\mathcal{L}_{\beta}$ is the $\mathbb{C}$-algebra generated by all the jump operators $\{S_{a}(\omega)\}_{a,w}$ is the whole algebra $\mathcal{B}(\mathcal{H})$. By this condition, one can check that for the choice of $\{S_{a}\}_{a\in\mathcal{A}}=\{\sigma_{i}^{x},\sigma_{i}^{y},\sigma_{i}^{z}\}_{i=1}^{N}$, the associated Davies generator is always primitive. Without loss of generality, when discussing the Gibbs samplers in Appendix A and Section 4, we always let $\{\sigma_{i}^{x},\sigma_{i}^{y},\sigma_{i}^{z}\}_{i=1}^{N}$ be a subset of $\{S_{a}\}_{a\in\mathcal{A}}$ to guarantee the primitivity.

The following lemmas are collected from [2, Lemmas 1 and 2] with proof omitted, which shall be used repeatedly in our subsequent spectral gap analysis.

In this section, we provide a technical overview of the proof of 1. In the analysis, there are three main steps:

1. (a)

   Decompose $\mathcal{H}$ into logic and syndrome subspaces and the observable algebra $\mathcal{B}(\mathcal{H})$ correspondingly.

2. (b)

   Demonstrate efficient transition between logic subspaces.

3. (c)

   Demonstrate fast mixing inside the syndrome subspace.

The decomposition in the first step leverages the special structure of the stabilizer Hamiltonian, following the approach outlined in previous work by [2]. Specifically, for the 2D toric code, we decompose $\mathcal{H}=\mathcal{H}_{\rm logic}\otimes\mathcal{H}_{\rm syndrome}$ according to the eigendecomposition of $H^{\rm toric}$. Since $H^{\rm toric}$ has a four-dimensional ground state space, it encodes two logical qubits, making $\dim(\mathcal{H}_{\rm logic})=4$. The syndrome subspace $\mathcal{H}_{\rm syndrome}$ is then spanned by the electric and magnetic excited states, which are characterized by the bond configurations of the local observables in $H^{\rm toric}$. That is, it can be further decomposed into the electric and magnetic excited subspaces $\mathcal{H}_{\rm syndrome}=\mathcal{H}^{\rm m}_{\rm b}\otimes\mathcal{H}^{\rm e}_{\rm b}$. According to the decomposition of $\mathcal{H}$, we can naturally decompose the observable algebra $\mathcal{B}(\mathcal{H})=\mathcal{Q}_{1}\otimes\mathcal{Q}_{2}\otimes\mathcal{A}^{\rm full}_{\rm m}\otimes\mathcal{A}^{\rm full}_{\rm e}$, where $\mathcal{Q}_{1}\otimes\mathcal{Q}_{2}$ is generated by the logical operators, i.e., the global operators $\mathsf{X}_{1},\mathsf{Z}_{1},\mathsf{X}_{2},\mathsf{Z}_{2}$ appearing in (1.1). The syndrome algebras $\mathcal{A}^{\rm full}_{\rm m}$ and $\mathcal{A}^{\rm full}_{\rm e}$ are spanned by linear transformations acting on the syndrome subspaces $\mathcal{H}^{\rm m}_{\rm b}$ and $\mathcal{H}^{\rm e}_{\rm b}$, respectively. Here, the $\otimes$ symbol represents multiplication between commuting matrices, and every element in $\mathcal{Q}_{1},\mathcal{Q}_{2},\mathcal{A}^{\rm full}_{\rm m}$, and $\mathcal{A}^{\rm full}_{\rm e}$ is understood as a matrix defined over the entire Hilbert space $\mathcal{H}$. We put the detailed discussion of the above decomposition in Section 4.1.

After decomposing $\mathcal{B}(\mathcal{H})$, we can further decompose the Davies generator (1.1):

$$\mathcal{L}_{\beta}=\underbrace{\mathcal{L}_{\rm local}+\mathcal{L}_{\rm global}}_{:=\mathcal{L}^{\rm gapped}}+\mathcal{L}^{\rm rest}\,,$$

with $\mathcal{L}^{\rm gapped}$ defined in (4.13) and $\mathcal{L}^{\rm rest}=\mathcal{L}_{\text{\rm local full}}-\mathcal{L}_{\rm local}$. Here we can show that $\mathcal{L}^{\rm gapped}$ is ergodic, i.e., $\mathrm{Ker}(\mathcal{L}^{\rm gapped})={\rm Span}\{{\bf 1}\}$. By Lemma 3 (item 2), we can lower bound $\mathrm{Gap}(\mathcal{L}_{\beta})$ by $\mathrm{Gap}(\mathcal{L}^{\rm gapped})$. More importantly, the generator $\mathcal{L}^{\rm gapped}$ is block diagonal with respect to the following decomposition:

$$\mathcal{B}(\mathcal{H})=\bigoplus_{B_{1}\in\{{\bf 1},\mathsf{X}_{1},\mathsf{Y}_{1},\mathsf{Z}_{1}\},B_{2}\in\{{\bf 1},\mathsf{X}_{2},\mathsf{Y}_{2},\mathsf{Z}_{2}\}}\mathcal{B}_{B_{1},B_{2}}\,,$$

where $\mathcal{B}_{B_{1},B_{2}}=B_{1}\otimes B_{2}\otimes\mathcal{A}^{\rm full}_{\rm m}\otimes\mathcal{A}^{\rm full}_{\rm e}$ and the operators $\mathsf{X}_{i},\mathsf{Y}_{i},\mathsf{Z}_{i}$ are given in (4.4). Because $\mathrm{Ker}(\mathcal{L}^{\rm gapped})={\rm Span}\{{\bf 1}\}$ and $\mathcal{L}^{\rm gapped}$ is block diagonal, we have

$$\mathrm{Gap}\left(-\mathcal{L}_{\beta}\right)\geq\mathrm{Gap}\left(-\mathcal{L}^{\rm gapped}\right)\\ \geq\min\left\{\mathrm{Gap}\left(-\mathcal{L}^{\rm gapped}\middle|_{\mathcal{B}_{{\bf 1},{\bf 1}}}\right),\lambda_{\min}\left(-\mathcal{L}^{\rm gapped}|_{\mathcal{B}_{B_{1},B_{2}}}\right)\middle|B_{1}\neq{\bf 1}\ \text{or}\ B_{2}\neq{\bf 1}\right\}\,.$$

(3.1)

Thus, for a lower bound estimate of $\mathrm{Gap}(-\mathcal{L}^{\rm gapped})$, it suffices to consider $\mathrm{Gap}(-\mathcal{L}^{\rm gapped}|_{\mathcal{B}_{{\bf 1},{\bf 1}}})$ and $\lambda_{\min}(-\mathcal{L}^{\rm gapped})|_{\mathcal{B}_{B_{1},B_{2}}}$ with $B_{1}\neq{\bf 1}\ \text{or}\ B_{2}\neq{\bf 1}$. The latter term characterizes the transition rate between different logical subspaces. If $\mathcal{L}_{\beta}$ contains only local coupling operators, the system requires a long evolution time to overcome the energy barrier and transit between different logical subspaces, which leads to a slow mixing time scaling as $\exp(\Theta(\beta))$ according to [2]. In our work, an important observation is that the global logical operators $\mathsf{X}_{1},\mathsf{Z}_{1},\mathsf{X}_{2},\mathsf{Z}_{2}$ can directly flip the logical qubits, enabling transitions between different logical subspaces without the need to overcome the energy barrier. The efficient transition in logic density space implies a fast decaying of $\exp(t\mathcal{L}^{\rm gapped})$ in the logic subspace $\{\mathcal{B}_{B_{1},B_{2}}\}_{B_{1}\neq{\bf 1}\ \text{or}\ B_{2}\neq{\bf 1}}$ and a lower bound of $\{-\mathcal{L}^{\rm gapped}|_{\mathcal{B}_{B_{1},B_{2}}}\}_{B_{1}\neq{\bf 1}\ \text{or}\ B_{2}\neq{\bf 1}}$. Specifically, in 8 in Section 4.2, we show

$$-\mathcal{L}^{\rm gapped}|_{\mathcal{B}_{B_{1},B_{2}}}\succeq\mathrm{Gap}\left(-\mathcal{L}_{\rm local}\middle|_{\mathcal{B}_{{\bf 1},{\bf 1}}}\right)\,.$$

(3.2)

According to the above analysis, the remaining thing to lower bound the spectral gap of $\mathcal{L}_{\beta}$ is to study the spectral gap of $\mathcal{L}_{\rm local}$ on $\mathcal{B}_{{\bf 1},{\bf 1}}$. Roughly speaking, this requires us to prove 2. Similar task is done in [2]. However, we emphasize that the lower bound and proof technique in [2] is not suitable for our purpose. In [2, Proposition 2], while the spectral gap is independent of the system size, it decays as $\exp(-\Theta(\beta))$, indicating slow mixing in the syndrome subspace at low temperatures. A main contribution of this work is to show that this lower bound is not tight when $\beta\gg 1$, and the system actually mixes fast in the syndrome subspace at low temperatures. To this end, we develop a new iteration argument and a decomposition trick to prove that $\mathrm{Gap}\left(\mathcal{L}_{\rm local}\middle|_{\mathcal{B}_{{\bf 1},{\bf 1}}}\right)\geq\min\left\{\exp(-\Theta(\beta)),\mathrm{poly}(1/N)\right\}$, which provides a much sharper lower bound of the spectral gap in the lower temperature regime. This result is summarized in 9 in Section 4.2.1, which also provides a proof of 2. Plugging this into (3.1) and (3.2), we can conclude 1.

Let $N=2L^{2}$ spins be on the edges of the toroidal lattice, modeled by the Hilbert space $\mathcal{H}\cong\mathbb{C}^{2^{N}}$. The Hamiltonian is given by

$$H^{\rm toric}=-\sum_{s}\textbf{X}_{s}-\sum_{p}\textbf{Z}_{p}\,,$$

(4.1)

where indices $s$ and $p$ denote a star and plaquette that consist of four sites around a node of the lattice and the center of a cell, respectively, as drawn in Fig. 1, and the associated observables $\textbf{X}_{s}$ and $\textbf{Z}_{p}$ are given by

$$\textbf{X}_{s}=\prod_{i\in s}\sigma^{x}_{i}\,,\quad\textbf{Z}_{p}=\prod_{i\in p}\sigma^{z}_{i}\,,$$

which commute with each other: $[\textbf{X}_{s},\textbf{Z}_{p}]=[\textbf{X}_{s},\textbf{X}_{s^{\prime}}]=[\textbf{Z}_{p},\textbf{Z}_{p^{\prime}}]\equiv 0$ for any stars $s,s^{\prime}$ and plaquettes $p,p^{\prime}$. In addition, due to the periodic boundary condition, it holds that

$$\prod_{s}\textbf{X}_{s}=1\,,\quad\prod_{p}\textbf{Z}_{p}=1\,,$$

(4.2)

This section is devoted to the spectral gap analysis of the Davies generator (Gibbs sampler) in (1.1) and the proof of 1, building on the discussion in Section 3. We will first discuss the ground states and the decomposition of the observable algebra of $H^{\rm toric}$ in Section 4.1. Then, in Section 4.2, we first address the step (b) outlined in Section 3, reducing the problem to step (c): the study of the spectral gap of the local Davies generator on the syndrome space. This will be explored in detail in Section 4.2.1.