Arbitrary-time error suppression for Markovian adiabatic quantum computing using stabilizer subspace codes

TL;DR

This paper demonstrates that error suppression in adiabatic quantum computing can be maintained at arbitrary times under thermal noise, using stabilizer codes with polynomial evolution time and local interactions.

Contribution

It extends initial error suppression results to arbitrary times, showing that efficient quantum computation remains feasible under thermal environments with stabilizer codes.

Findings

Thermal excitations can be suppressed at any time with polynomial evolution.

Error suppression requires at least 4-local stabilizer interactions.

Energy penalty strength grows logarithmically with system size.

Abstract

Adiabatic quantum computing (AQC) can be protected against thermal excitations via an encoding into error detecting codes, supplemented with an energy penalty formed from a sum of commuting Hamiltonian terms. Earlier work showed that it is possible to suppress the initial thermally induced excitation out of the encoded ground state, in the case of local Markovian environments, by using an energy penalty strength that grows only logarithmically in the system size, at a fixed temperature. The question of whether this result applies beyond the initial time was left open. Here we answer this in the affirmative. We show that thermal excitations out of the encoded ground state can be suppressed at arbitrary times under the additional assumption that the total evolution time is polynomial in the system size. Thus, computational problems that can be solved efficiently using AQC in a closed system setting, can still be solved efficiently subject to coupling to a thermal environment. Our construction uses stabilizer subspace codes, which require at least $4$-local interactions to achieve this result.

Equations

$$H(t)=\sum _{l\ge 0}{ϵ}_l(t){\Pi }_l(t),$$

$$p_{\perp }$$

$${\dot{p}}_{\perp }=-{\partial }_t\mathrm{Tr}({\Pi }_0\rho )\equiv -R(t)$$

$$p_{\perp }(t)=-{\int }_0^{t}R(t^{\prime })d{t}^{\prime }+p_{\perp }(0).$$

$$R(t)={\partial }_t\mathrm{Tr}({\Pi }_0\rho )=\mathrm{Tr}({\dot{\Pi }}_0\rho )+\mathrm{Tr}({\Pi }_0\dot{\rho }).$$

$$\mathrm{Tr}({\dot{\Pi }}_0\rho )$$

$$+\mathrm{Tr}({\dot{\Pi }}_0{\Pi }_0^{\perp }\rho {\Pi }_0)+\mathrm{Tr}({\dot{\Pi }}_0{\Pi }_0\rho {\Pi }_0^{\perp }).$$

$${\Pi }_0{\dot{\Pi }}_0+{\dot{\Pi }}_0{\Pi }_0={\dot{\Pi }}_0⟹{\dot{\Pi }}_0{\Pi }_0={\Pi }_0^{\perp }{\dot{\Pi }}_0.$$

$${\Pi }_0{\dot{\Pi }}_0{\Pi }_0=0,{\Pi }_0^{\perp }{\dot{\Pi }}_0{\Pi }_0^{\perp }=0.$$

$$R(t)=\mathrm{Tr}({\Pi }_0\dot{\rho })\mathrm{if}\rho ={\Pi }_0\rho {\Pi }_0+{\Pi }_0^{\perp }\rho {\Pi }_0^{\perp }\forall t.$$

$$R(t)=\mathrm{Tr}({\Pi }_0\dot{\rho })+[\mathrm{Tr}({\dot{\Pi }}_0{\Pi }_0\rho {\Pi }_0^{\perp })+\text{c.c.}].$$

$$H_{\mathrm{tot}}(t)=H_S(t)+H_B+H_{SB},$$

$$H_{SB}=\sum _{\alpha }{A}_{\alpha }\otimes B_{\alpha }$$

$$\dot{\rho }=\mathcal{L}(t)[\rho ]=-i[H(t)+H_{\mathrm{LS}}(t),\rho ]+D(t)[\rho ],$$

$$⟨{\mathcal{B}}_{\alpha {\alpha }^{\prime }}(t)⟩=\mathrm{Tr}({\rho }_B{e}^{-i{H}_{B}t}{B}_{\alpha }{e}^{i{H}_{B}t}{B}_{{\alpha }^{\prime }}),$$

$${\gamma }_{\alpha {\alpha }^{\prime }}(\omega )={\int }_{-\infty }^{\infty }dt{e}^{i\omega t}⟨{\mathcal{B}}_{\alpha {\alpha }^{\prime }}(t)⟩={\gamma }_{{\alpha }^{\prime }\alpha }^{\ast }(\omega ).$$

$$F_{\alpha }=\sum _{{\alpha }^{\prime }}(u^{†}{)}_{\alpha {\alpha }^{\prime }}{A}_{{\alpha }^{\prime }},$$

$$F_{\alpha }(\omega )=\sum _{{ϵ}_{l^{\prime }}-{ϵ}_l=\omega }{\Pi }_l{F}_{\alpha }{\Pi }_{l^{\prime }},$$

$$D[\rho ]=\sum _{\omega \alpha }{\gamma }_{\alpha }(\omega )[F_{\alpha }(\omega )\rho F_{\alpha }^{†}(\omega )-\frac{1}{2}\{F_{\alpha }^{†}(\omega )F_{\alpha }(\omega ),\rho \}],$$

$${\gamma }_{\alpha }(-\omega )=e^{-\beta \omega }{\gamma }_{\alpha }(\omega ),\omega >0.$$

$$|\mathrm{Tr}({\dot{\Pi }}_0{\Pi }_0\rho {\Pi }_0^{\perp })|$$

$$\le \Vert S\dot{H}\Vert \Vert {\Pi }_0\rho {\Pi }_0^{\perp }{\Vert }_1$$

$$\le \frac{\Vert \dot{H}\Vert }{\Delta (H)}\Vert {\Pi }_0\rho {\Pi }_0^{\perp }{\Vert }_1,$$

$$\Vert {\Pi }_0^{\perp }\rho {\Pi }_0{\Vert }_1$$

$$\le \Vert {\Pi }_0^{\perp }\Vert \Vert \rho -{\rho }^{\mathrm{ad}}{\Vert }_1\Vert {\Pi }_0\Vert \le C/t_f,$$

$$|\mathrm{Tr}[{\dot{\Pi }}_0{\Pi }_0\rho {\Pi }_0^{\perp }]|\le \frac{\Vert \dot{H}\Vert C_v}{\Delta (H)t_f^{v+1}},$$

$$\mathrm{Tr}({\Pi }_0\dot{\rho })=-i\mathrm{Tr}({\Pi }_0[H+H_{\mathrm{LS}},\rho ]+\mathrm{Tr}({\Pi }_{0}D[\rho ]).$$

$$\mathrm{Tr}({\Pi }_0[H^{\prime },\rho ])=\mathrm{Tr}([{\Pi }_0,H^{\prime }]\rho )=0.$$

$$\mathrm{Tr}({\Pi }_{0}D[\rho ])$$

$$\displaystyle\qquad-\frac{1}{2}\mathrm{Tr}[\Pi_{0}\{F_{\alpha}^{\dagger}(\omega)F_{\alpha}(\omega),\rho\}]\bigg{)}\ .$$

$$\sum _{\omega }{\gamma }_{\alpha }(\omega )\mathrm{Tr}[{\Pi }_0{F}_{\alpha }(\omega )\rho F_{\alpha }^{†}(\omega ){\Pi }_0]$$

$$=\sum _{\omega }{\gamma }_{\alpha }(\omega )\sum _{{ϵ}_{l^{\prime }}-{ϵ}_l=\omega }\sum _{{ϵ}_{l^{\prime \prime \prime }}-{ϵ}_{l^{\prime \prime }}=\omega }$$

$$\mathrm{Tr}[{\Pi }_0({\Pi }_l{F}_{\alpha }{\Pi }_{l^{\prime }})\rho ({\Pi }_{l^{\prime \prime \prime }}{F}_{\alpha }^{†}{\Pi }_{l^{\prime \prime }}){\Pi }_0]$$

$$=\sum _{\omega }{\gamma }_{\alpha }(\omega )\sum _{{ϵ}_{l^{\prime }}=\omega +{ϵ}_0}\sum _{{ϵ}_{l^{\prime \prime \prime }}=\omega +{ϵ}_0}$$

Full text

I Introduction

In adiabatic quantum computing (AQC), computations are performed using a time-dependent Hamiltonian that evolves smoothly from an initial Hamiltonian with a known and easily preparable ground state, to a final Hamiltonian whose ground state is unknown and encodes the desired result Farhi et al. (2000) (for a review see Ref. Albash and Lidar (2018)). This model appears promising for near-future large scale realization, especially in terms of (non-universal) quantum annealing devices, which already feature a few thousand qubits Harris et al. (2018); King et al. (2018).

Despite enjoying a certain degree of inherent robustness to errors, AQC has four main and well documented failure modes Childs et al. (2001); Sarandy and Lidar (2005); Aberg et al. (2005); S. Ashhab, J. R. Johansson, and F. Nori (2006); Tiersch and Schützhold (2007); Amin et al. (2008, 2009a, 2009b); Deng et al. (2013); Sarovar and Young (2013): (i) diabatic transitions out of the ground state arising from an evolution on a timescale that is faster than that set by the inverse gap, (ii) control errors resulting in the implementation of the wrong final Hamiltonian, (iii) decoherence of the ground state, and (iv) thermal excitations out of the ground state. The first of these is a purely unitary error mode which arises even in the absence of coupling to the environment. It is mitigated by slowing the evolution down in accordance with the adiabatic theorem Kato (1950), i.e., in order to remain in the ground state throughout, the total evolution time is required to be large relative to the timescale set by (a small power of) the inverse of the smallest energy gap from the ground state encountered along the evolution Jansen et al. (2007). The second can be viewed as arising from technical imperfections or from the environment; either way it can be mitigated to some extent by imposing smooth boundary conditions on the interpolation between the initial and final Hamiltonians Lidar et al. (2009); Wiebe and Babcock (2012); Ge et al. (2016); Campos Venuti and Lidar (2018) or by encoding the final Hamiltonian Young et al. (2013a), but the absence of a complete theory of fault tolerance in AQC (despite impressive attempts Mizel (2014)) means that it is not currently known how to scalably and reliably overcome control errors. The third and fourth are entirely environment-induced errors. Decoherence of the ground state (due to decoherence in the computational basis) is a catastrophic failure mode that occurs when the coupling to the environment is too strong for AQC to be meaningfully executed. Quantum error correction methods can be deployed in principle, but at present they are impractical in that they require the use of many-body interactions that scale with the problem size Young et al. (2013b). To avoid decoherence of the ground state, AQC should be performed in systems obeying the weak coupling limit to the environment, where decoherence occurs in the instantaneous energy eigenbasis Albash and Lidar (2015).

In this work we address thermal excitations. This failure mode can be suppressed using a scheme first proposed by Jordan, Farhi, and Shor (JFS) Jordan et al. (2006). In the JFS scheme, an error detecting stabilizer subspace code is chosen, and the system Hamiltonian is encoded using the logical operators of the same code. A penalty Hamiltonian $H_{\mathrm{p}}$ proportional to the sum of the stabilizer generators of the code is added, which suppresses excitations out of the code subspace. This is useful since without encoding thermal excitations are suppressed only by the gap of the system Hamiltonian, but with encoding, thermal excitations are suppressed by the gap of the penalty Hamiltonian (which is a constant for stabilizer codes) times the magnitude $\eta_{\mathrm{p}}$ of the energy penalty.

In their analysis, JFS assumed a particular system of spins weakly coupled to a photon bath and a pure initial state. They then identified the lowest-weight possible subspace stabilizer codes for detecting $1$-local and $2$-local noise compatible with the suppression of thermal excitation errors. Ref. Marvian and Lidar (2017a) generalized the JFS suppression result to arbitrary Markovian master equations and arbitrary subspace (as opposed to subsystem) error detection codes, while allowing for mixed initial states. However, both Refs. Jordan et al. (2006); Marvian and Lidar (2017a) only considered the ultra-short-time performance of this error suppression scheme for Markovian environments. More precisely, they established conditions for the success of the scheme only in terms of the initial thermal excitation rate out of the code subspace.

Here we complete the analysis initiated in Ref. Marvian and Lidar (2017a) and consider the performance of subspace-based error suppression schemes for arbitrary times $t$. We prove that thermal excitation errors can be suppressed for all physically reasonable, local Markovian environments by increasing $\eta_{\mathrm{p}}$ only logarithmically in the number of qubits $n$ at constant bath temperature, provided the total evolution time scales at most polynomially in $n$. Our main technical result is formulated in terms of an upper bound on the excited state population at arbitrary $t$, assuming that the system is initialized in the ground subspace. We show that, provided the conditions mentioned above hold, this bound can be made arbitrarily small by increasing $\eta_{\mathrm{p}}$ in proportion to $\log(n)$. Since we require that the total evolution time $t_{f}\sim\mathrm{poly}(n)$, our result does not guarantee protection against thermal excitation errors for problems with exponentially (or superpolynomially) small gaps, for which, by the adiabatic theorem, we expect $t_{f}$ to have to scale faster than $\mathrm{poly}(n)$.

The structure of this paper is as follows. In Sec. II we provide a general bound on the excitation rate out of the ground subspace at arbitrary time. We observe that the bound involves an off-diagonal component (coherence between the ground and excited subspaces) that did not appear in the earlier initial-time treatment of Refs. Jordan et al. (2006); Marvian and Lidar (2017a). In Sec. III we derive an upper bound on the excited state population after encoding using a subspace-based error detection code and adding an energy penalty term, and show that it can be made arbitrarily small provided the penalty strength $\eta_{\mathrm{p}}\sim\log(n)$ and total evolution time $t_{f}\sim\mathrm{poly}(n)$. This is the content of our main result, Eq. (52). Readers interested primarily in the conclusions can skip many details of the derivation and read the paper starting from this point. We provide a summary and discussion in Sec. IV, and provide a few additional technical details in the Appendix.

II Bounding the excitation rate out of the ground subspace at arbitrary time

Consider the spectral decomposition of a time-dependent Hamiltonian $H(t)$:

[TABLE]

where $\Pi_{l}(t)$ denotes the projection onto the (possibly degenerate) $H(t)$-eigensubspace with eigenvalue $\epsilon_{l}(t)$. The eigenvalues are ordered so that $\epsilon_{0}(t)\leq\epsilon_{1}(t)\leq\dots$ $\forall t$, and we assume that there are no level crossings. The eigenprojectors are orthogonal: $\Pi_{l}(t)\Pi_{l^{\prime}}(t)=\delta_{ll^{\prime}}\Pi_{l}(t)$. From now on we usually drop the explicit time-dependence to simplify the notation. But it important to remember that all our quantities are explicitly time-dependent unless explicitly stated otherwise.

II.1 General expression for the excitation rate out of the ground subspace

Assume that the system is initially prepared in the (possibly degenerate) ground subspace of $H$, with energy $\epsilon_{0}$, i.e., $\rho(0)=\Pi_{0}\rho(0)\Pi_{0}$. Using $\rho=(\Pi_{0}+\Pi_{0}^{\perp})\rho(\Pi_{0}+\Pi_{0}^{\perp})$, the population in the subspace orthogonal to $\Pi_{0}$ is

[TABLE]

so that

[TABLE]

and

[TABLE]

Since we assumed that the initial population is fully in $\Pi_{0}$, i.e., $\mathrm{Tr}[\Pi_{0}\rho(0)\Pi_{0}]=1$, it follows that $p_{\perp}(0)=0$.

In contrast to Ref. Marvian and Lidar (2017a), which focused on the initial excitation rate out of the ground subspace, $R(0)=\partial_{t}\mathrm{Tr}(\Pi_{0}{\rho})|_{t\simeq 0}$, here we are interested in the excitation rate for arbitrary $t$

[TABLE]

Again using $\rho=(\Pi_{0}+\Pi_{0}^{\perp})\rho(\Pi_{0}+\Pi_{0}^{\perp})$, we have

[TABLE]

The terms in line (6a) vanish: differentiating the identity $\Pi_{0}^{2}=\Pi_{0}$ yields

[TABLE]

Multiplying from the left or from the right gives

[TABLE]

The two complex conjugate terms in line (6b) are due to coherence between $\Pi_{0}$ and $\Pi_{0}^{\perp}$. They vanish if we assume that at arbitrary time $t$ the state is in $\Pi_{0}\oplus\Pi_{0}^{\perp}$:

[TABLE]

Since the initial state satisfies $\rho(0)=\Pi_{0}\rho(0)\Pi_{0}$, Eq. (9) holds at $t=0$ without an additional no-coherence assumption. This was the case studied in Ref. Marvian and Lidar (2017a). But since here we are interested in arbitrary $t$, we have:

[TABLE]

One way for the assumption $\rho=\Pi_{0}\rho\Pi_{0}+\Pi_{0}^{\perp}\rho\Pi_{0}^{\perp}\ \forall t$ to hold is if decoherence between $\Pi_{0}$ and $\Pi_{0}^{\perp}$ is fast on the timescale of the evolution, i.e., $T_{2}\ll t_{f}$, where $T_{2}$ is the timescale over which $\Pi_{0}\rho\Pi_{0}^{\perp}$ and $\Pi_{0}^{\perp}\rho\Pi_{0}$ decay, and $t_{f}$ is the final time, i.e., the total evolution time. This is certainly true in the time-independent case (where we expect these coherences to decay at least as fast as $e^{-t/T_{2}}$), but it does not hold in the general time-dependent case, as we discuss below.

II.2 Adiabatic master equation in Davies-Lindblad form

Let us define our open system model. Assuming a total Hamiltonian of the form

[TABLE]

where $H_{S}(t)$ is the time-dependent system Hamiltonian, $H_{B}$ is a general bath Hamiltonian, and

[TABLE]

is a general system-bath interaction Hamiltonian, an adiabatic Markovian master equation in Davies-Lindblad form Davies (1974); Lindblad (1976) can be derived in the weak coupling limit Albash et al. (2012):

[TABLE]

where $H_{\mathrm{LS}}(t)$ is the Lamb shift, which commutes with $H_{S}(t)$, and $D(t)$ denotes the dissipative (non-unitary) part. Henceforth we use units such that $\hbar\equiv 1$. We assume that the system operators $A_{\alpha}$ are $\ell$-local, with $\ell$ a constant that is independent of the number of system particles (e.g., qubits) $n$. The interaction Hamiltonian $H_{SB}$ then has a local structure, and can be expressed as a sum over $\binom{n}{\ell}$ terms, which is polynomial in $n$.

We briefly review the structure of $D(t)$ (see, e.g., Ref. Albash and Lidar (2015) for more details). Let $\rho_{B}$ denote the initial state of the bath. The bath correlation function is

[TABLE]

and its Fourier transform is

[TABLE]

The matrix $\gamma=\{\gamma_{\alpha\alpha^{\prime}}\}$ is positive semi-definite. Therefore it can be diagonalized by a unitary matrix $u$. Define new system operators

[TABLE]

and their transforms

[TABLE]

where the sum is over all pairs of eigenvalues $\epsilon_{l^{\prime}}$ and $\epsilon_{l}$ whose difference is equal to the given Bohr frequency $\omega$. Then, the dissipator $D(t)$ can be written as:

[TABLE]

where $\{\gamma_{\alpha}(\omega)\}$ are the eigenvalues of $\gamma(\omega)$.

If the bath is in thermal equilibrium at inverse temperature $\beta=1/(k_{B}T)$, then the matrix of decay rates satisfies the Kubo-Martin-Schwinger (KMS) condition Haag et al. (1967):

[TABLE]

It relates the excitation rate $\gamma_{\alpha}(-\omega)$ to the relaxation rate $\gamma_{\alpha}(\omega)$, and shows that excitation is exponentially suppressed in $\beta\omega$ relative to relaxation.

II.3 Upper bound on the off-diagonal term

To bound the off-diagonal term $\mathrm{Tr}(\dot{\Pi}_{0}\Pi_{0}\rho\Pi_{0}^{\perp})$ in Eq. (10), we note that $\dot{\Pi}_{0}$ can be replaced by the reduced resolvent $S=\Pi_{0}^{\perp}(H-\epsilon_{0})^{-1}\Pi_{0}^{\perp}$. Namely, it is well known that $\dot{\Pi}_{0}\Pi_{0}=-S\dot{H}\Pi_{0}$ and $\|S\|\leq 1/\Delta(H)$, the inverse minimum spectral gap of $H(t)$, i.e., the gap between $\Pi_{0}$ and the first excited state (see, e.g., Appendices B and F of Ref. Rezakhani et al. (2010)). This gives the following bound:

[TABLE]

where $\|\cdot\|$ is the operator norm (largest singular value) and $\|\cdot\|_{1}$ is the trace-norm, and we used the inequality $|\mathrm{Tr}(AB)|\leq\|A\|\|B^{\dagger}\|_{1}$ (operator norm of $A$ times trace norm of $B^{\dagger}$) R. Bhatia (1997).

How tightly can we bound the $\|\Pi_{0}^{\perp}\rho\Pi_{0}\|_{1}$ factor in general? As mentioned above, for the time-independent case and in the weak coupling limit this quantity decays rapidly due to decoherence between eigenstates. But in the time-dependent case the best we can do in general is the following. Let $\rho^{\mathrm{ad}}$ be the solution of of the Davies-Lindblad master equation in the adiabatic limit. This is known to be the Gibbs state Venuti et al. (2016), which is diagonal in the energy eigenbasis. Hence $\Pi_{0}^{\perp}\rho^{\mathrm{ad}}\Pi_{0}=0$, and using $\|\Pi_{0}\|=\|\Pi_{0}^{\perp}\|=1$:

[TABLE]

where we used the adiabatic theorem for open systems Venuti et al. (2016) for the last inequality. The $t_{f}$-independent constant $C$ is given by Eq. (4) of  Venuti et al. (2016) and depends on a power of $1/\Delta(\mathcal{L})$ where $\mathcal{L}$ is the Lindbladian from Eq. (13) and $\Delta(\mathcal{L})$ is the minimal gap from the zero eigenvalue of $\mathcal{L}(s)$ along the evolution from $s=0$ to $s=1$, where $s=t/t_{f}$.

The bound (21b) can be tightened using boundary cancellation methods, and improved to $C_{v}/t_{f}^{v+1}$, where $v$ is the number of vanishing derivatives of $H(t)$ at $t_{f}$, and the bounded, $t_{f}$-independent constants $C_{v}$ are given in Eq. (11b) of Ref. Campos Venuti and Lidar (2018). It is important to note that $C_{v}$ does depend on the system size $n$, a point we return to below. Thus, altogether we have:

[TABLE]

where $\|\dot{H}\|$ and $\Delta(H)$ are maximized and minimized over the interval $[0,t_{f}]$, respectively.

II.4 Computation of the diagonal term

Now we consider the diagonal term in Eq. (10):

[TABLE]

Let us compute each term in turn.

II.4.1 Computation of $\mathrm{Tr}(\Pi_{0}[H+H_{\mathrm{LS}},\rho])$

Let $H^{\prime}\equiv H+H_{\mathrm{LS}}$. Recall that $[H,H_{\mathrm{LS}}]=0$ and $[\Pi_{0},H^{\prime}]=0$. Now note that $\Pi_{0}[H^{\prime},\rho]=[H^{\prime},\Pi_{0}\rho]+[\Pi_{0},H^{\prime}]\rho$, so that:

[TABLE]

Thus there is no contribution from the unitary part.

II.4.2 Computation of $\mathrm{Tr}(\Pi_{0}D[\rho])$

Eq. (18) gives:

[TABLE]

The term in line (25a) is:

[TABLE]

For the term in line (25b), we first note that:

[TABLE]

An identical result holds for $F_{\alpha}^{\dagger}(\omega)F_{\alpha}(\omega)\Pi_{0}$. Using this, the term in line (25b) becomes $-\sum_{\l}\gamma_{\alpha}(\epsilon_{0}-\epsilon_{l})\mathrm{Tr}[{\Pi_{0}\rho\Pi_{0}F^{\dagger}_{\alpha}}{\Pi_{l}F_{\alpha}}]$. Therefore, after using the KMS condition to write $\gamma_{\alpha}(\epsilon_{0}-\epsilon_{l})=e^{-\beta(\epsilon_{l}-\epsilon_{0})}\gamma_{\alpha}(\epsilon_{l}-\epsilon_{0})$:

[TABLE]

These terms represent opposite processes: the first represents relaxation into $\Pi_{0}$, the second represents excitation out of $\Pi_{0}$. When just the initial excitation rate is accounted for, only the excitation term appears Marvian and Lidar (2017a).

II.4.3 The case without degeneracies: Pauli master equation

In the absence of any degeneracies the projectors are all rank $1$, i.e., $\Pi_{l}=|{\epsilon_{l}}\rangle\langle\epsilon_{l}|$ and we can simplify Eq. (II.4.2) by factoring out the populations

[TABLE]

This yields:

[TABLE]

where we used the KMS condition and defined the Markov transition matrix (whose elements are positive):

[TABLE]

In this case $\mathrm{Tr}(\Pi_{0}\dot{\rho})=\mathrm{Tr}(\Pi_{0}D[\rho])$ is simply the Pauli master equation, expressing repopulation of the ground state $\ket{\epsilon_{0}}$ with transition rate $W_{0l}$, and depopulation with transition rate $W_{l0}=e^{-\beta(\epsilon_{l}-\epsilon_{0})}W_{0l}$ (detailed balance). Both positivity and the detailed balance conditions are proved in Appendix A.

III Excitation rate reduction using an error detecting code

III.1 The encoded Hamiltonian

We now choose a code $C$ that can detect all the errors (system operators) $F_{\alpha}$ Knill and Laflamme (1997):

[TABLE]

where $P_{C}$ projects onto the code space. More explicitly, we choose $C$ to be an $[[n,k,d]]$ stabilizer code, where the number of physical and logical qubits is $n$ and $k$, respectively, and $d$ is the code distance ($d\geq 2$). Let us denote the stabilizer generators by $\{S_{i}\}_{i=1}^{n-k}$. This partitions the Hilbert space into $2^{n-k}$ stabilizer syndrome subspaces, each of which is $2^{k}$-dimensional. Each such subspace is defined by a particular ordered assignment of $\pm 1$ eigenvalues of the stabilizer generators Gottesman (1996).

We construct a penalty Hamiltonian $H_{\mathrm{p}}$ by summing the stabilizer generators Jordan et al. (2006):

[TABLE]

The eigenvalues of $H_{\mathrm{p}}$ are

[TABLE]

with corresponding $r_{m}=\binom{n-k}{m}2^{k}$-dimensional eigensubspaces having the property that exactly $m$ of the stabilizer generators have eigenvalue $-1$. Note that $\sum_{m=0}^{n-k}r_{m}=2^{k}\sum_{m=0}^{n-k}\binom{n-k}{m}=2^{n}$, as required.

We define the codespace $C$ as usual as the linear span of the simultaneous eigenstates of all the stabilizer generators with eigenvalue $+1$ Gottesman (1996). Therefore the codespace is the ground subspace of $H_{\mathrm{p}}$. The codespace ($m=0$) is $r_{0}=2^{k}$-dimensional, supporting $k$ logical qubits. The code therefore has $k$ logical $X$ and $k$ logical $Z$ operators denoted $\bar{X}_{i}$ and $\bar{Z}_{i}$, respectively, and their products form logical $XX$ and $ZZ$ operators. We encode the system Hamiltonian $H_{S}(t)$ using these logical operators. The following encoded Hamiltonian is universal for AQC Biamonte and Love (2008):

[TABLE]

so we may assume this form without loss of generality. The entire time-dependence is in the parameters $\{h_{i}^{x}(t),h_{i}^{z}(t),J^{x}_{ij}(t),J^{z}_{ij}(t)\}$.

The total system Hamiltonian is then

[TABLE]

where the dimensionless quantity $\eta_{\mathrm{p}}>0$ quantifies the strength of the energy penalty. Every error detected by the code anticommutes with at least one stabilizer generator [this is equivalent to Eq. (32)] Gottesman (1996), so every such error “pays” an energy penalty equal to the number of anticommuting generators times $\eta_{\mathrm{p}}$. Adding the encoded system Hamiltonian $\bar{H}_{S}(t)$ to $\eta_{\mathrm{p}}H_{\mathrm{p}}$ splits the $2^{k}$-fold degeneracy of the code space. The ground state of $H(t)$ is the encoded ground state which the error suppression scheme is designed to protect. The protection is against bath-induced errors that excite the system out of the code space ($m=0\mapsto m>0$), but not against errors that induce transitions inside the code space ($m=0\mapsto m=0$). The latter can be either logical errors due to coupling of the bath to system operators with weight $\geq d$, or due to diabatic transitions arising from non-adiabatic evolution. These considerations, as well as additional ones discussed below, are illustrated in Fig. 1.

Since by construction $[H_{\mathrm{p}},\bar{H}_{S}(t)]=0$ $\forall t$, it follows that we may express the two Hamiltonians in terms of the same set of time-dependent eigenprojectors:

[TABLE]

where the $\{\xi_{l}\}_{l}$ and $\{\bar{\omega}_{l}(t)\}_{l}$ are, respectively, the eigenvalues of $H_{\mathrm{p}}$ and $\bar{H}_{S}(t)$. Since $\bar{H}_{S}(t)$ breaks the degeneracy of $H_{\mathrm{p}}$, the index $l$ ranges over a set of values that is at least as large as that of the index $m$ in Eq. (34). I.e., unlike in Eq. (34) where $m\neq m^{\prime}\Longrightarrow\xi_{m}\neq\xi_{m^{\prime}}$, the index $l$ may repeat certain eigenvalues of $H_{\mathrm{p}}$, and $\xi_{l}=\xi_{l^{\prime}}$ is possible when $l\neq l^{\prime}$.111 It may seem surprising that the time-independent $H_{\mathrm{p}}$ can be expressed in terms of a linear combination of time-dependent projectors. To see this explicitly, let $\{\xi_{m},\{\ket{v_{j}^{(m)}}\}_{j=1}^{r_{m}}\}_{m}$ denote the time-independent eigenvalues and corresponding eigenvectors of $H_{\mathrm{p}}$, where $r_{m}$ is the degeneracy of $\xi_{m}$ [recall Eq. (34)]. Since $[H_{\mathrm{p}},\bar{H}_{S}(t)]=0$, $\bar{H}_{S}(t)$ preserves the eigenspaces of $H_{\mathrm{p}}$ but may break the degeneracy within each such subspace, and the eigenvectors of $\bar{H}_{S}(t)$ are time-dependent linear combinations of the eigenvectors of $H_{\mathrm{p}}$. I.e., we can write the eigenvectors of $\bar{H}_{S}(t)$ as $\{\{\ket{w_{k}^{(l)}(t)}=\sum_{j}a_{jk}^{(l)}(t)\ket{v_{j}^{(l)}}\}_{k}\}_{l}$, with corresponding time-dependent eigenvalues $\{\bar{\omega}_{l}(t)\}_{l}$. We can trivially check that the time-dependent vectors $\{\{\ket{w_{k}^{(l)}(t)}\}_{k}\}_{l}$ are also eigenvectors of $H_{\mathrm{p}}$ with eigenvalues $\{\xi_{l}\}_{l}$. Namely, $H_{\mathrm{p}}\ket{w_{k}^{(l)}(t)}=\sum_{j}a_{jk}^{(l)}(t)H_{\mathrm{p}}\ket{v_{j}^{(l)}}=\sum_{j}a_{jk}^{(l)}(t)\xi_{l}\ket{v_{j}^{(l)}}=\xi_{l}\ket{w_{k}^{(l)}(t)}$. Therefore we may express $H_{\mathrm{p}}$ as in Eq. (37), with the projectors explicitly identified as $\Pi_{l}(t)=\sum_{k}|{w_{k}^{(l)}(t)}\rangle\langle w_{k}^{(l)}(t)|$.

We may now choose the $\Pi_{l}$’s of Eq. (1) as the same eigenprojectors and write the eigendecomposition of $H(t)$ as

[TABLE]

As mentioned above, the codespace is the ground space of $H_{\mathrm{p}}$. We may associate a corresponding projection operator

[TABLE]

The fact that we have a sum over $l$ is due to the breaking of the degeneracy of the codespace by $\bar{H}_{S}(t)$; the sum is over all $l$ required to span the $2^{k}$-dimensional codespace. Eq. (38) tells us the ground subspace projector of $H(t)$ is $\Pi_{0}$, i.e., the ground space of $H(t)$ is a subspace of the code space.

III.2 Modified excitation rate after encoding

We assume that the initial state, which is now encoded into $C$ and evolves according to $H(t)$ of Eq. (36), belongs to $\Pi_{0}$, i.e., again $\rho(0)=\Pi_{0}\rho(0)\Pi_{0}$. This means that $p_{\perp}(0)=0$ in Eq. (4), and we can focus exclusively on the excitation rate $R(t)$. Let us collect the results from above [Eqs. (10), (22)-(24), (II.4.2)], as follows.

First, we note that the off-diagonal term $Q(t)\equiv\mathrm{Tr}(\dot{\Pi}_{0}\Pi_{0}\rho\Pi_{0}^{\perp})+\textrm{c.c.}$ satisfies the bound given in Eq. (22):

[TABLE]

Next:

[TABLE]

Thanks to the error detection properties of the code [Eq. (32)] we have

[TABLE]

so that the sum over $l$ reduces to a sum only over terms not in the codespace:

[TABLE]

The term $m_{\alpha}^{0,l}$ represents the excitation out of the codespace associated with error operator $F_{\alpha}$, while $m_{\alpha}^{l,0}$, represents the corresponding relaxation into the codespace.

III.3 The non-codespace population is exponentially suppressed by the energy penalty

Next, let us show that for reasonable models of the bath the excitation rate is exponentially suppressed with increasing energy penalty $\eta_{\mathrm{p}}$.

First, repeating the argument given in Ref. Marvian and Lidar (2017a), let $\Pi_{l}$ denote an eigenprojector of $H$ with energy $\epsilon_{l}=\mathrm{Tr}[\Pi_{l}H]$. Recall that the $\Pi_{l}$’s are simultaneous eigenprojectors of $\bar{H}_{S}$ and $H_{\mathrm{p}}$ as well. We have $\forall l\in C^{\perp}$:

[TABLE]

where $g$ is the ground state gap of $H_{\mathrm{p}}$:

[TABLE]

When $H_{\mathrm{p}}$ is a sum of commuting terms, as is true for the stabilizer construction we consider here, the gap $g$ is guaranteed to be a constant Bravyi and Vyalyi (2005).

Next, note that, using the spectral decomposition $\rho(t)=\sum_{i}\lambda_{i}|{i}\rangle\langle i|$, both the excitation and relaxation terms in Eq. (43a) are positive:

[TABLE]

Using Eqs. (4) and (43a), the non-codespace population is

[TABLE]

We may replace $e^{-\beta(\epsilon_{l}-\epsilon_{0})}$ by the upper bound $e^{-\beta\eta_{\mathrm{p}}g}$, and further increase the RHS by removing the relaxation term, since it is positive. Hence:

[TABLE]

where we defined

[TABLE]

and used the fact that we already know that both ${\tilde{\gamma}_{\max}}$ and $\tilde{m}_{\alpha}^{0,l}$ are positive. To obtain line (48b) we used Eq. (38). Note the appearance of the factor $t$ due to the integration in Eq. (47). This factor will play an important role in our final upper bound considerations below, and was absent from the initial-time-only considerations of Refs. Jordan et al. (2006); Marvian and Lidar (2017a).

The bound on $|p_{\perp}|$ depends on ${\tilde{\gamma}_{\max}}$. To ensure a non-trivial bound this quantity has to be finite, which is a natural assumption. We also assume that $\gamma(\omega)$ is a polynomial (or any subexponential) function of $\omega$ for $\omega>0$; this too is an assumption that is compatible with all commonly used bath models Breuer and Petruccione (2002). Therefore $\tilde{\gamma}_{\max}\sim\mathrm{poly}(\eta_{\mathrm{p}})$.

What remains to be shown is that the sum over all non-code states—which might appear to involve exponentially many terms—does not spoil this conclusion. Now, since we already showed [Eq. (46)] that each $m_{\alpha}^{0,l}\geq 0$, we have:

[TABLE]

where we used $\sum_{l}\Pi_{l}=I$. Using the inequality $|\mathrm{Tr}(AB)|\leq\|A\|\|B^{\dagger}\|_{1}$ again, we have

[TABLE]

since $\|\Pi_{0}\rho\Pi_{0}\|_{1}\leq 1$. Thus:

[TABLE]

The sum over $\alpha$ contains a polynomial number of terms in $n$ due to our earlier assumption of $\ell$-local system operators in the system-bath interaction Hamiltonian $H_{SB}=\sum_{\alpha}{A_{\alpha}\otimes B_{\alpha}}$. Each $\|F_{\alpha}\|$ is itself a sum over $\mathrm{poly}(n)$ $n$-independent terms due to Eq. (16).

Combining all this with Eqs. (40) and (47b) we thus conclude that

[TABLE]

Line (52a) states that the non-codespace population is exponentially suppressed in terms of the energy penalty, and grows at most linearly in time. This contribution to the bound is due to the error suppression strategy, and is similar to the result in Ref. Marvian and Lidar (2017a). The new aspect is the factor of $t_{f}$. The dependence of $t_{f}$ on problem size $n$ is dictated by the adiabatic theorem: $t_{f}$ must scale as a small inverse power ($2$ or $3$) of the minimum gap $\Delta(H)$ encountered along the evolution Jansen et al. (2007); Lidar et al. (2009). The same condition also ensures that the denominator in line (52b) grows with $n$. This line arises purely due to diabatic transitions, which cannot be suppressed using the energy penalty or encoding. Note that one factor of $t_{f}$ present in Eq. (40) cancelled after integration, so that we must enforce $v\geq 1$.

To ensure that the entire contribution of line (52b) decreases with $n$, the number of vanishing derivatives $v$ at $t=t_{f}$ must be sufficiently large to overcome both the scaling of $\|\dot{H}\|$ with $n$ (which is at most quadratic, in the case of all-to-all interactions), and the scaling of $C_{v}$ with $n$. The latter depends on powers of both $\Delta(H)$ and $\Delta(\mathcal{L})$, since it is known that in the Davies-Lindblad adiabatic master equation case and in the presence of a thermal bath satisfying the KMS condition, the constant $C$ appearing in the adiabatic theorem bound without assuming any boundary conditions [Eq. (21b)] satisfies $C=O[\|\mathcal{L}^{\prime}\|/\Delta^{2}(\mathcal{L})]$ (where prime denotes differentiation with respect to $s$) Venuti:2017aa.

Therefore, assuming that all inverse gap dependencies are polynomial in $n$, Eq. (52) implies that as long as $t_{f}\sim\mathrm{poly}(n)$ (with an appropriately high degree) then one only needs to increase the strength of the energy penalty, $\eta_{\mathrm{p}}$, logarithmically in $n$, at any fixed inverse temperature $\beta$.

To see explicitly why, let $\mathrm{poly}(n)\sim n^{p}$, $t_{f}\sim n^{q}$, and $\eta_{\mathrm{p}}\sim r\log(n)$. Then

[TABLE]

so that if $r>(p+q)/(\beta g)$ then this guarantees that line (52a) decreases (polynomially) in the system size $n$. The $\mathrm{poly}(\eta_{\mathrm{p}})$ factor does not change this conclusion since it is a polynomial in $\log(n)$.

If, on the other hand, the inverse gap dependencies are exponential in $n$ then we need $t_{f}\sim\exp(n)$, and then $\eta_{\mathrm{p}}$ must grow at least linearly in $n$ in order for error suppression to be effective, which is unacceptable: the same suppression could be achieved simply by scaling up all the coupling constants linearly without incurring the cost of encoding.

The main conclusion reported in Refs. Jordan et al. (2006); Marvian and Lidar (2017a) therefore remains valid for arbitrary evolution times, namely, that by using error detecting codes built on commuting Hamiltonians (for which $g$ is constant), for physically plausible Markovian models, a logarithmically increasing energy penalty strength suffices for error suppression. The main new caveat is that the proof holds for problems with polynomially small gaps, but not for problems with gaps that decrease superpolynomially.

IV Summary and Discussion

In this paper we studied error suppression for Markovian models. These are interesting, despite the fact that general results of a similar nature have already been established for non-Markovian models Bookatz et al. (2015); Marvian (2016), since they are widely used Breuer and Petruccione (2002), and moreover decay in these models is always exponential R. Alicki and K. Lendi (1987). In this sense error suppression in the non-Markovian case is less challenging than for Markovian AQC, since in the latter case one cannot rely on the use of non-Markovian recurrences, as is commonly done in error suppression techniques such as dynamical decoupling Lidar (2008); Quiroz and Lidar (2012); Ganti et al. (2014) or the Zeno effect Paz-Silva et al. (2012).

Ref. Marvian and Lidar (2017a) generalized the JFS result Jordan et al. (2006), that it suffices to increase the energy penalty logarithmically with system size in order to protect AQC against excitations out of the ground state, to general Markovian dynamics and mixed states. However, these earlier results were only valid for the initial excitation rate, and the natural question of whether they generalize to arbitrary evolution times was left open. Here we settled this question in the affirmative, under the assumption that the problem gap is at most polynomially small in the problem size. While it seems unlikely, we did not rule out the possibility that this method of error suppression could be adapted to work even for problems with an exponentially small gap. Still, the present result establishes that one of the main failure modes of AQC can be overcome under plausible physical assumptions for problems for which AQC is efficient.

We emphasize that our results require the ability to encode both the final and the initial Hamiltonian. Therefore they do not apply to transverse field implementations of quantum annealing, where only the final Hamiltonian can be encoded Pudenz et al. (2014); Vinci et al. (2016). A further caveat is that it is known that for penalty Hamiltonians comprising a sum of commuting Hamiltonians, as is the case here, at least $3$-local interactions are required Marvian and Lidar (2014), and for stabilizer subspace codes at least $4$-local interaction are required for universality Jordan et al. (2006). However, we expect that our entire construction will generalize straightforwardly to the stabilizer subsystem setting Poulin (2005), where the penalty Hamiltonian becomes a sum over the $2$-local generators of the (non-Abelian) gauge group of the code Jiang and Rieffel (2017); Marvian and Lidar (2017b). We expect this to improve upon the locality of the construction presented here as well. The generalization of our results to the subsystem code case is an important problem left for future work.

Acknowledgements.

This work benefited greatly from numerous constructive discussions with Milad Marvian, who carefully read and improved the manuscript. Thanks also to Paolo Zanardi and Lorenzo Campos Venuti for insightful comments. The research is based upon work (partially) supported by the Office of the Director of National Intelligence (ODNI), Intelligence Advanced Research Projects Activity (IARPA), via the U.S. Army Research Office contract W911NF-17-C-0050. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of the ODNI, IARPA, or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright annotation thereon.

Appendix A Properties of the Markov transition matrix

We defined the Markov transition matrix in Eq. (31): $W_{l^{\prime}l}=\sum_{\alpha}\gamma_{\alpha}(\epsilon_{l}-\epsilon_{l^{\prime}})\mathrm{Tr}[{\Pi_{l}F^{\dagger}_{\alpha}\Pi_{l^{\prime}}F_{\alpha}}]$.

Positivity of the matrix elements can be seen by recalling Eq. (16) and that the $A_{\alpha}$ are Hermitian. Then:

[TABLE]

so that:

[TABLE]

To prove that detailed balance holds, let us write:

[TABLE]

where $\gamma_{\alpha^{\prime\prime}\alpha^{\prime}}(\epsilon_{l}-\epsilon_{l^{\prime}})=\sum_{\alpha}u_{\alpha^{\prime\prime}\alpha}\gamma_{\alpha}(\epsilon_{l}-\epsilon_{l^{\prime}})(u^{\dagger})_{\alpha\alpha^{\prime}}$ and $A_{l^{\prime}l,\alpha}\equiv\bra{\epsilon_{l^{\prime}}}A_{\alpha}\ket{\epsilon_{l}}$. Thus

[TABLE]

But using the the KMS condition with $\epsilon_{l^{\prime}}>\epsilon_{l}$, we have from Eq. (56c)

[TABLE]