Variational preparation of finite-temperature states on a quantum computer

Abstract

The preparation of thermal equilibrium states is important for the simulation of condensed matter and cosmology systems using a quantum computer. We present a method to prepare such mixed states with unitary operators and demonstrate this technique experimentally using a gate-based quantum processor. Our method targets the generation of thermofield double states using a hybrid quantum-classical variational approach motivated by quantum-approximate optimization algorithms, without prior calculation of optimal variational parameters by numerical simulation. The fidelity of generated states to the thermal-equilibrium state smoothly varies from 99 to 75% between infinite and near-zero simulated temperature, in quantitative agreement with numerical simulations of the noisy quantum processor with error parameters drawn from experiment.

INTRODUCTION

The potential for quantum computers to simulate other quantum mechanical systems is well known¹, and the ability to represent the dynamical evolution of quantum many-body systems has been demonstrated². However, the accuracy of these simulations depends on efficient initial state preparation within the quantum computer. Much progress has been made on the efficient preparation of non-trivial quantum states, including spin-squeezed states³ and entangled cat states⁴. Studying phenomena like high-temperature superconductivity⁵ requires preparation of thermal equilibrium states, or Gibbs states. Producing mixed states with unitary quantum operations and measurements is not straightforward, and has only recently begun to be explored^6,7. In this work, we demonstrate the use of a variational quantum-classical algorithm to realize Gibbs states using (ideally unitary) gate control on a transmon quantum processor.

Our approach is mediated by the generation of thermofield double (TFD) states, which are pure states sharing entanglement between two identical quantum systems with the characteristic that when one of the systems is considered independently (by tracing over the other), the result is a mixed state representing equilibrium at a specific temperature. TFD states are of interest not only in condensed matter physics but also for the study of black holes^8,9 and traversable wormholes^10,11. We use a variational protocol¹² motivated by quantum-approximate optimization algorithms (QAOA) that relies on alternation of unitary intra- and inter-system operations to control the effective temperature, eliminating the need for a large external heat bath. Other methods have been studied for generation of Gibbs states, such as quantum metropolis sampling¹³ and imaginary time evolution using variational quantum simulation^14,15. However, the advantage of QAOA compared to these proposals is that the form of the ansatz is relatively straightforward and low-depth, whereas the metropolis sampling involves phase estimation which leads to a high-depth circuit, and the imaginary time evolution proposal does not have a clear proposal for the form of the ansatz. Recently, verification of TFD state preparation was demonstrated on a trapped-ion quantum computer⁶. Our work experimentally demonstrates the generation of finite-temperature states in a superconducting quantum computer by variational preparation of TFD states in a hybrid quantum-classical manner.

RESULTS

Theory

Consider a quantum system described by Hamiltonian H with eigenstates $\left|j\right\rangle$ and corresponding eigenenergies E_j:

$$H\left|j\right\rangle ={E}_{j}\left|j\right\rangle .$$

(1)

The Gibbs state ρ_Gibbs of the system is

$${\rho }_{{{{\rm{Gibbs}}}}}(\beta )=\frac{1}{Z}\mathop{\sum}\limits_{j}\exp \left(-\beta {E}_{j}\right)\left|j\right\rangle \left\langle j\right|,$$

(2)

where β = 1/k_BT is the inverse temperature, k_B is the Boltzmann constant, and

$$Z=\mathop{\sum}\limits_{j}\exp \left(-\beta {E}_{j}\right)$$

(3)

is the partition function. Except in the limit β → ∞, the Gibbs state is a mixed state and thus impossible to generate strictly through unitary evolution. To circumvent this, we define the TFD state¹² on two identical systems A and B as

$$\left|{{{\rm{TFD}}}}(\beta )\right\rangle =\frac{1}{\sqrt{Z}}\mathop{\sum}\limits_{j}\exp \left(\frac{-\beta {E}_{j}}{2}\right){\left|j\right\rangle }_{{{{\rm{B}}}}}\ {\left|j\right\rangle }_{{{{\rm{A}}}}}.$$

(4)

Tracing out either system yields the desired Gibbs state in the other.

To prepare the TFD states, we follow the variational protocol proposed by¹² and consider two systems each of size n. In the first step of the procedure, the TFD state at β = 0 is generated by creating Bell pairs $\left|{{{\Phi }}}_{i}^{+}\right\rangle =\left({\left|0\right\rangle }_{{{{\rm{B}}}}i}{\left|0\right\rangle }_{{{{\rm{A}}}}i}+{\left|1\right\rangle }_{{{{\rm{B}}}}i}{\left|1\right\rangle }_{{{{\rm{A}}}}i}\right)/\sqrt{2}$ between corresponding qubits i in the two systems. Tracing out either system yields a maximally mixed state on the other, and vice versa. The next steps to create the TFD state at finite temperature depend on the relevant Hamiltonian. Here, we choose the transverse field Ising model in a one-dimensional chain of n spins¹⁶, with n = 2 [Fig. 1(a)]. We map spin up (down) to the computational state $\left|0\right\rangle$$(\left|1\right\rangle )$ of the corresponding transmon. The Hamiltonian describing system A is

$${H}_{{{{\rm{A}}}}}={ZZ}_{{{{\rm{A}}}}}+g{X}_{{{{\rm{A}}}}},$$

(5)

where Z Z_A = Z_A2Z_A1, X_A = X_A2 + X_A1, and g is proportional to the transverse magnetic field. The Hamiltonian for system B is the same. We focus on g = 1, where a phase transition is expected in the transverse field Ising model at large n¹⁷. We use a QAOA-motivated variational ansatz^12,18, where intra-system evolution is interleaved with a Hamiltonian enforcing interaction between the systems:

$${H}_{{{{\rm{BA}}}}}={XX}_{{{{\rm{BA}}}}}+{ZZ}_{{{{\rm{BA}}}}},$$

(6)

where X X_BA = X_B2X_A2 + X_B1X_A1, and analogously for ZZ_BA. For single-step state generation, the unitary operation describing the TFD protocol is

$$U\left({{{\boldsymbol{\alpha }}}},{{{\boldsymbol{\gamma }}}}\right)={U}_{{{{\rm{inter}}}}}\left({{{\boldsymbol{\alpha }}}}\right){U}_{{{{\rm{intra}}}}}\left({{{\boldsymbol{\gamma }}}}\right),$$

(7)

where

$$\begin{array}{l}{U}_{{{{\rm{intra}}}}}({{{\boldsymbol{\gamma }}}})=\exp \left[-{{{\rm{i}}}}{\gamma }_{2}\left({ZZ}_{{{{\rm{B}}}}} + {ZZ}_{{{{\rm{A}}}}}\right)/2\right]\\ \qquad\qquad\quad\times \exp \left[-{{{\rm{i}}}}{\gamma }_{1}\left({X}_{{{{\rm{B}}}}}+{X}_{{{{\rm{A}}}}}\right)/2\right],\end{array}$$

(8)

$${U}_{{{{\rm{inter}}}}}({{{\boldsymbol{\alpha }}}})=\exp \left(-{{{\rm{i}}}}{\alpha }_{2}{ZZ}_{{{{\rm{BA}}}}}/2\right)\exp \left(-{{{\rm{i}}}}{\alpha }_{1}{XX}_{{{{\rm{BA}}}}}/2\right).$$

(9)

The variational parameters ${{{\boldsymbol{\gamma }}}}=\left({\gamma }_{1},{\gamma }_{2}\right)$, ${{{\boldsymbol{\alpha }}}}=\left({\alpha }_{1},{\alpha }_{2}\right)$ are optimized by the hybrid classical-quantum algorithm to generate states closest to the ideal TFD states. A single step of intra- and inter-system interaction ideally produces the state $\left|\psi ({{{\boldsymbol{\alpha }}}},{{{\boldsymbol{\gamma }}}})\right\rangle =U\left({{{\boldsymbol{\alpha }}}},{{{\boldsymbol{\gamma }}}}\right)\left(\left|{{{\Phi }}}_{2}^{+}\right\rangle \otimes \left|{{{\Phi }}}_{1}^{+}\right\rangle \right)$¹⁹.

Fig. 1: Principle and generation of Thermofield Double state.

a Two identical systems A and B are variationally prepared in an ideally pure, entangled joint state such that tracing out one system yields the Gibbs state on the other. b Corresponding qubits in the two systems are first pairwise entangled to produce the β = 0 TFD state. Next, intra- and inter-system Hamiltonians are applied with optimized variational angles $\left({{{\boldsymbol{\alpha }}}},{{{\boldsymbol{\gamma }}}}\right)$ to approximate the TFD state corresponding to the desired temperature.

The variational algorithm extracts the cost function after each state preparation. We engineer a cost function ${{{\mathcal{C}}}}$ to be minimized when the generated state is closest to an ideal TFD state¹⁹. Following recent work on the concentration of control parameters for QAOA^20,21, we expect engineering the cost function based on the target state for small-sized systems to lead to a general expression for the cost function of an arbitrary-sized system. The engineered cost function is given by:

$$\begin{array}{l}C(\beta )=\langle {X}_{{{{\rm{A}}}}}\rangle +\langle {X}_{{{{\rm{B}}}}}\rangle +1.57\left(\langle {ZZ}_{{{{\rm{A}}}}}\rangle +\langle {ZZ}_{{{{\rm{B}}}}}\rangle \right)\\ -{\beta }^{-1.57}(\langle {XX}_{{{{\rm{BA}}}}}\rangle +\langle {ZZ}_{{{{\rm{BA}}}}}\rangle ).\end{array}$$

(10)

We compare the performance of this engineered cost function ${{{{\mathcal{C}}}}}_{1.57}$ to that of the non-optimized cost function ${{{{\mathcal{C}}}}}_{1.00}$, using the reduction of infidelity to the Gibbs state as the ultimate metric of success [see Supplementary Note 1]. The engineered cost function achieves an average improvement of 54% across the β range covered ([10⁻², 10²] in units of 1/g), as well as a maximum improvement of up to 98% for intermediate temperatures (β ~ 1). Our choice of the class of cost functions to optimize lets us trade off a slight decrease in low-temperature performance with a significant increase in performance at intermediate temperatures. See ref. ¹⁹ for further details on the theory.

The quantum portion of the algorithm prepares the state according to a given set of angles $\left({{{\boldsymbol{\alpha }}}},{{{\boldsymbol{\gamma }}}}\right)$, performs the measurements, and returns these values to the classical portion. The classical portion then evaluates the cost function according to the returned measurements, performs classical optimization, generates and returns the next set of variational angles to evaluate the quantum portion.

Experiment

We implement the algorithm using four of seven transmons in a monolithic quantum processor [Fig. 2(a)]. The four transmons (labeled A₁, A₂, B₁, and B₂) have square connectivity provided by coupling bus resonators, and are thus ideally suited for implementing the circuit in Fig. 1(b). Each transmon has a microwave-drive line for single-qubit gating, a flux-bias line for two-qubit controlled-Z (CZ) gates, and a dispersively coupled resonator with dedicated Purcell filter^22,23. The four transmons can be simultaneously and independently read out by frequency multiplexing, using the common feedline connecting to all Purcell filters. All transmons are biased to their flux-symmetry point (i.e., sweetspot²⁴) using static flux bias to counter residual offsets. Device details and a summary of measured transmon parameters are provided in Supplementary Note 3. Details on the experimental setup can be found in Supplementary Note 4.

Fig. 2: Device and optimized quantum circuit.

a Optical image of the transmon processor was used in this experiment, with false color highlighting the four transmons employed and the dedicated bus resonators providing their nearest-neighbor coupling. Scale bar corresponds to 1 mm. b Optimized circuit equivalent to that in Fig. 1(b) and conforming to the native gate set in our architecture. All variational parameters are mapped onto rotation axes and angles of single-qubit gates. Tomographic pre-rotations R₁–R₄ are added to reconstruct the terms in the cost function ${{{\mathcal{C}}}}$ and to perform two-qubit state tomography of each system following optimization.

In order to realize the theoretical circuit in Fig. 1(b), we first map it to the optimized depth-13 equivalent circuit shown in Fig. 2(b), which conforms to the native gate set in our control architecture. This gate set consists of arbitrary single-qubit rotations about any equatorial axis of the Bloch sphere, and CZ gates between nearest-neighbor transmons. Conveniently, all variational angles are mapped to either the axis or angle of single-qubit rotations. Further details on the compilation steps are reported in Methods and Supplementary Note 2. Bases pre-rotations are added at the end of the circuit to first extract all the terms in the cost function ${{{\mathcal{C}}}}$ and finally to perform two-qubit state tomography of each system.

Prior to implementing any variational optimizer, it is helpful to build a basic understanding of the cost-function landscape. To this end, we investigate the cost function ${{{\mathcal{C}}}}$ at β = 0 using two-dimensional cuts Fig. 3: we sweep γ while keeping α = 0 to study the effect of U_intra and vice versa to study the effect of U_inter. Note that owing to the β^−1.57 divergence, the cost function reduces to −〈H_BA〉 in the β = 0 limit. Consider first the landscape for an ideal quantum processor, which is possible to compute for our system size. The γ landscape at α = 0 is π-periodic in both directions due to the invariance of $\left|{{{\rm{TFD}}}}(\beta =0)\right\rangle$ under bit-flip (X) and phase-flip (Z) operations on all qubits. The cost function is minimized to −4 at even multiples of π/2 on γ₁ and γ₂: $\left|{{{\rm{TFD}}}}(\beta =0)\right\rangle$ is a simultaneous eigenstate of X X_BA and Z Z_BA with eigenvalue +2 due to the symmetry of the constituting Bell states $\left|{{{\Phi }}}_{i}^{+}\right\rangle$. In turn, the cost function is maximized to +4 at odd multiples of π/2, at which the $\left|{{{\Phi }}}_{i}^{+}\right\rangle$ are transformed to singlets $\left|{{{\Psi }}}_{i}^{-}\right\rangle =\left({\left|0\right\rangle }_{{{{\rm{B}}}}i}{\left|1\right\rangle }_{{{{\rm{A}}}}i}-{\left|1\right\rangle }_{{{{\rm{B}}}}i}{\left|0\right\rangle }_{{{{\rm{A}}}}i}\right)/\sqrt{2}$. The α landscape at γ = 0 is constant, reflecting that $\left|{{{\rm{TFD}}}}(\beta =0)\right\rangle$ is a simultaneous eigenstate of XX_BA and ZZ_BA and thus also of any exponentiation of these operators. The corresponding experimental landscapes show qualitatively similar behavior. The γ landscape clearly shows the π periodicity with respect to both angles, albeit with reduced contrast. The α landscape is not strictly constant, showing weak structure particularly with respect to α₂. These experimental deviations reflect underlying errors in our noisy intermediate-scale quantum (NISQ) processor, which include transmon decoherence, residual ZZ coupling at the bias point, and leakage during CZ gates. We discuss these error sources in detail further below.

Fig. 3: Landscape of the cost function for infinite temperature.

Panels (a) and (c): Landscape cuts obtained in noiseless simulation and in experiment, respectively, when varying control angles γ while keeping α = 0 (the ideal solution value). Panels (b) and (d): Corresponding cuts obtained when varying control angles α while keeping γ = 0 (the ideal solution value). See text for details. These landscape cuts are sampled at 100 points and interpolated using the Python package adaptive⁴¹.

The task of the variational algorithm is to balance the mixture of the states at each β, in order to generate the corresponding Gibbs state. Although thermal states are well understood, it is challenging to accurately generate them in NISQ devices for studies of finite temperature systems. When working with small systems, it is possible and tempting to predetermine the variational parameters at each β by a prior classical simulation and optimization for an ideal or noisy quantum processor. We refer to this common practice^6,25 as cheating, since this approach does not scale to larger problem sizes and skips the main quality of variational algorithms: to arrive at the parameters variationally. Here, we avoid cheating altogether by starting at β = 0, with initial guess the obvious optimal variational parameters for an ideal processor (γ = α = 0), and using the experimentally optimized $\left({{{\boldsymbol{\alpha }}}},{{{\boldsymbol{\gamma }}}}\right)$ at the last β as an initial guess when stepping β in the range $\left[0,5\right]$ (in units of 1/g). This approach only relies on the assumption that solutions (and their corresponding optimal variational angles) vary smoothly with β. At each β, we use the Gradient-Based Random-Tree optimizer of the scikit-optimize²⁶ Python package to minimize ${{{\mathcal{C}}}}$, using 4096 averages per tomographic pre-rotation necessary for the calculation of ${{{\mathcal{C}}}}$. After 200 iterations, the optimization is stopped. The best point is remeasured two times, each with 16384 averages per tomographic pre-rotation needed to perform two-qubit quantum state tomography of each system. A new optimization is then started for the next β, using the previous solution as the initial guess.

To begin comparing the optimized states ${\rho }_{{{{\rm{Exp}}}}}$ produced in an experiment to the target Gibbs states ρ_Gibbs, we first visualize their density matrices (in the computational basis) for a sampling of the β range covered (Fig. 4). Starting from the maximally-mixed state I I/4 at β = 0, the Gibbs state monotonically develops coherences (off-diagonal terms) between all states as β increases. Coherences between states of equal (opposite) parity have 0 (π) phase throughout. Populations (diagonal terms) monotonically decrease (increase) for even (odd) parity states. By β = 5, the Gibbs state is very close to the pure state $\left|{{\Upsilon }}\right\rangle \left\langle {{\Upsilon }}\right|$, where $\left|{{\Upsilon }}\right\rangle \approx \sqrt{0.36}\left(\left|01\right\rangle +\left|10\right\rangle \right)-\sqrt{0.14}\left(\left|00\right\rangle +\left|11\right\rangle \right)$. The noted trends are reproduced in ${\rho }_{{{{\rm{Exp}}}}}$. However, the matching is evidently not perfect, and to address this we proceed to a quantitative analysis.

Fig. 4: Qualitative comparison of optimized states to the Gibbs state.

Visualization of the density matrices (in the computational basis) for the targeted Gibbs states ρ_Gibbs (left) and the optimized experimental states ${\rho }_{{{{\rm{Exp}}}}}$ (right) at (a, b) β = 0, a–d β = 1 and (e–f) β = 5. As β increases, the Gibbs state monotonically develops coherence between all states, with phase 0 (π) for states with the same (opposite) parity. Populations in even (odd) parity states decrease (increase). The optimized experimental states show qualitatively similar trends.

We employ two metrics to quantify experimental performance: the fidelity F of ${\rho }_{{{{\rm{Exp}}}}}$ to ρ_Gibbs and the purity P of ${\rho }_{{{{\rm{Exp}}}}}$, given by

$$F=\,{{\mbox{Tr}}}\,\ \left(\sqrt{{\rho }_{{{{\rm{Gibbs}}}}}^{1/2}{\rho }_{{{{\rm{Exp}}}}}^{\phantom{1/2}}{\rho }_{{{{\rm{Gibbs}}}}}^{1/2}}\right),$$

(11)

$$P=\,{{\mbox{Tr}}}\,\ \left({\rho }_{{{{\rm{Exp}}}}}^{2}\right).$$

(12)

At β = 0, F = 99% and P = 0.262, revealing a very close match to the ideal maximally-mixed state. However, F smoothly worsens with increasing β, decreasing to 92% at β = 1 and 75% by β = 5. Simultaneously, P does not closely track the increase of purity of the Gibbs state. By β = 5, the Gibbs state is nearly pure, but P peaks at 0.601.

In an effort to quantitatively explain these discrepancies, we perform a full density-matrix simulation of a four-qutrit system using quantumsim²⁷. Our simulation incrementally adds calibrated errors for our NISQ processor, starting from an ideal processor (model 0): transmon relaxation and dephasing times at the bias point (model 1), increased dephasing from flux noise during CZ gates (model 2), crosstalk from residual ZZ coupling at the bias point (model 3), and transmon leakage to the second-excited state during CZ gates (model 4). The experimental input parameters for each increment and details of modeling are described in Methods and Supplementary Notes 5–9. The added curves in Fig. 5 clearly show that model 4 quantitatively matches the observed dependence of F and P over the full β range, and identifies leakage from CZ gates as the dominant error.

Fig. 5: Performance of the variational algorithm.

a Fidelity to the Gibbs state as a function of inverse temperature β for experimental states obtained by optimization and cheating. b Purity of experimental states as a function of β, and comparison to the purity of the Gibbs state. Added curves in both panels are obtained by numerical simulation of a noisy quantum processor with incremental error models based on calibrated error sources for our device: qubit relaxation and dephasing times, increased dephasing from flux noise during CZ gates, residual ZZ crosstalk at the bias point, and leakage during CZ gates. Leakage is identified as the dominant error source.

DISCUSSION

The power of variational algorithms relies on their adaptability: the optimizer is meant to find its way through the variational parameter space, adapting to mitigate coherent errors as allowed by the chosen parametrization. For completeness, we compare in Fig. 5 the performance achieved with our variational strategy to that achieved by cheating, i.e., using the pre-calculated optimal $\left({{{\boldsymbol{\alpha }}}},{{{\boldsymbol{\gamma }}}}\right)$ for an ideal processor. Our variational approach, whose sole input is the obvious initial guess at β = 0, achieves comparable performance at all β. This aspect is crucial when considering the scaling with problem size, as classical pre-simulations will require prohibitive resources beyond ~50 qubits, but variational optimizers would not. Given the dominant role of leakage as the error source, which cannot be compensated by the chosen parametrization, it is unsurprising in hindsight that both approaches yield nearly identical performance.

In summary, we have presented the generation of finite-temperature Gibbs states in a quantum computer by variational targeting of TFD states in a hybrid quantum-classical manner. The algorithm successfully prepares mixed states for the transverse field Ising model with Gibbs-state fidelity ranging from 99% to 75% as β increases from 0 to 5/g. The loss of fidelity with decreasing simulated temperature is quantitatively matched by a numerical simulation with incremental error models based on experimental input parameters, which identifies leakage in CZ gates as dominant. This work demonstrates the suitability of variational algorithms on NISQ processors for the study of finite-temperature problems of interest, ranging from condensed-matter physics to cosmology. Our results also highlight the critical importance of continuing to reduce leakage in two-qubit operations when employing weakly-anharmonic multi-level systems such as the transmon.

During the preparation of this manuscript, we became aware of related experimental work²⁵ on a trapped-ion system, applying a non-variationally prepared TFD state to the calculation of a critical point.

METHODS

Quantum circuit

We map the theoretical circuit in Fig. 1(b) to an equivalent circuit conforming to the native gate set in our control architecture and exploiting virtual Z-gate compilation²⁸ to minimize circuit depth. Single-qubit rotations R_XY(ϕ, θ), by arbitrary angle θ around any equatorial axis $\cos (\phi ){{{\bf{x}}}}+\sin (\phi ){{{\bf{y}}}}$ on the Bloch sphere, are realized using 20ns DRAG pulses^29,30. Two-qubit CZ gates are realized by baseband flux pulsing^31,32 using the Net Zero scheme^33,34, completing in 80ns. In the optimized circuit [Fig. 2(b)], CZ gates only appear in pairs. These pairs are simultaneously executed and tuned as one block. Single-qubit rotations R₁-R₄ are used to change the measurement bases, as required to measure ${{{\mathcal{C}}}}$ during optimization and to perform two-qubit tomography³⁵ in each system to extract F and P. A summary of single- and two-qubit gate performance [see Supplementary Note 5] and a step-by-step derivation of the optimized circuit is provided [see Supplementary Note 2].

Modeling and simulations

Noiseless simulations were performed prior to experiments, for verification of algorithm convergence. During simulations, experimental conditions were maintained exactly for the algorithm and the control software, while the outcome from readout hardware (Zurich Instruments UHFQC) was replaced with a simulated readout [see Supplementary Note 10 for details].

The models used to simulate the performance of the algorithm are incremental: model k contains all the noise mechanisms in model k − 1 plus one more, which we use for labeling in Fig. 5. Model 0 corresponds to an ideal quantum processor without any error. Model 1 adds the measured relaxation and dephasing times measured for the four transmons at their bias point. These times are tabulated in Supplementary Table 1. Model 2 adds the increased dephasing that flux-pulsed transmons experience during CZ gates. For this, we extrapolate the echo coherence time ${T}_{2}^{{{{\rm{echo}}}}}$ to the CZ flux-pulse amplitude using a 1/f noise model^36,37 with amplitude $\sqrt{A}=1\mu {{{\Phi }}}_{0}$. This noise model is implemented following³⁸. Model 3 adds the idling crosstalk due to residual ZZ coupling between transmons. This model expands on the implementation of idling evolution used for coherence times: the circuit gates are simulated to be instantaneous, and the idling evolution of the system is trotterized. In this case, the residual ZZ coupling operator uses the measured residual ZZ coupling strengths at the bias point [see Supplementary Note 6 for details]. Finally, model 4 adds leakage to the CZ gates, based on randomized benchmarking with modifications to quantify leakage^33,39, and implemented in simulation using the procedure described in ref. ³⁸.

Leakage in transmons

Leakage to transmon second-excited states is found essential to quantitatively match the performance of the algorithm by simulation. To reach this conclusion it was necessary to thoroughly understand how leakage affects the two-qubit tomographic reconstruction procedure employed. The readout calibration only considers computational states of the two transmons involved. Moreover, basis pre-rotations only act on the qubit subspace, leaving the population in leaked states unchanged. Using an overcomplete set of basis pre-rotations for state tomography, comprising both positive $\left(X,Y,Z\right)$ and negative $\left(-X,-Y,-Z\right)$ bases for each transmon, leads to the misdiagnosis of a leaked state as a maximally mixed state qubit state for that transmon. See Supplementary Note 8 for further details.