$\hat{a}_{1}\hat{3}\hat{b}^{\dagger}$ are comparable to $g_{2}\lesssim
2\pi\times 5\textrm{MHz}$. Incorporating the bosonic enhancement factor due to
the average excitation number $\alpha^{2}$, we require the detuning $\Delta$
to be much larger than $g_{2}\alpha^{2}$, e.g., $\Delta=10g_{2}\alpha^{2}$.
Then the achievable coupling strength of the quartic interaction is given by
$g_{2}^{2}/(10g_{2}\alpha^{2})=g_{2}/(10\alpha^{2})\lesssim 2\pi\times
60\textrm{kHz}$ assuming $\alpha^{2}=8$.

##  Appendix G Measurement

In this appendix we discuss measurement schemes for high fidelity readout in
both the $X$ and $Z$ basis. Compared to gates where the optimal gate errors
can be cleanly described in terms of dimensionless constants like
$\kappa_{1}/\kappa_{2}$, for readout absolute timescales matter more. We
enumerate the parameter choices for the different schemes at the end of the
corresponding sections.

###  G.1 $X$-basis measurement

Here we discuss in more detail the $X$-basis readout scheme used to generate
the infidelities used in most of the error correction simulations. Note that
throughout this appendix, when we refer to measurement infidelities or error
probabilities, we are referring to misassignment probabilities

$\displaystyle\epsilon_{s}=1-\text{P}(s|s)$ (262)

where $\text{P}(s|s)$ is the probability of reading out the state s given that
the cavity was in state $s$.

$X$-basis measurement refers to determining the parity of a phononic mode or
equivalently readout in the basis of even and odd cat states, i.e.,
$|\pm\rangle\propto|\alpha\rangle\pm|-\alpha\rangle$. To realize such a
measurement with minimal impact on the length of an error correction cycle we
utilize an additional phononic mode which we refer to as the readout mode.
This mode is interrogated by a transmon in parallel with the next error
correction cycle. As is pictured in Fig. 2, every unit cell contains this
additional readout mode connected to a transmon. Pictured in Fig. 39 is the
circuit we use for measuring an $X$ stabilizer which we now walk through in
more detail. To perform an $X$ stabilizer measurement first the ancilla qubit
$\hat{a}_{1}$ is entangled with the data qubits. Subsequently we “deflate” the
ancilla qubit on a timescale comaparable $1/\kappa_{2}$, mapping the even
parity state to $|\hat{n}=0\rangle$ and the odd parity state to
$|\hat{n}=1\rangle$ [27]. Deflation is achieved under evolution with the two-
phonon dissipator,

$\displaystyle\frac{d\hat{\rho}(t)}{dt}$
$\displaystyle=\kappa_{2}\mathcal{D}[\hat{a}_{1}^{2}-\alpha(t)^{2}]\hat{\rho}(t),$
(263)

by taking $\alpha(t)$ from $\alpha_{0}$ to $\alpha_{1}<\alpha_{0}$. In our
case we rapidly take $\alpha(t)$ from its initial value to $\alpha_{1}=0$
where we evolve for a time on the order of $1/\kappa_{2}$. The deflation is
not required to be adiabatic since we do not need to maintain phase coherence
between the even and odd parity states. The utility of the deflation is that
it makes the state of the cavity less susceptible to single phonon loss events
which change its parity.

After the deflation we perform a SWAP between the ancilla phononic mode
($\hat{a}_{1}$) and the readout mode ($\hat{a}_{2}$). The SWAP is performed
using the Hamiltonian

$\displaystyle\hat{H}_{rot}$
$\displaystyle=g_{r}(\hat{a}_{1}^{\dagger}\hat{a}_{2}+\hat{a}_{1}^{\dagger}\hat{a}_{2})$
(264)

Evolution under this Hamiltonian for a duration ${\pi}/{2g_{r}}$ realizes a
SWAP gate (there is a rotation of the swapped state by 90 degrees). Physical
implementation of this Hamiltonian is discussed in Eq. 250.

An advantage of this readout scheme is that after the exchange has occurred
the next cycle of quantum error correction can continue in parallel with the
measurement of the readout mode. With the idling time now set only by the
deflation + SWAP steps, we can spend more time measuring the readout mode
without compromising on idling error. This simple layout choice could be
generally useful in other architectures. In the specific case of this
proposal, we perform repeated QND parity measurements which we majority vote
to get our final measurement outcome [65, 66, 64]. In general more advanced
methods than majority voting will give higher fidelity [65].

Measurement of the readout mode parity is done using a dispersive coupling
with a transmon qubit [64] $H=-\chi\hat{\sigma}_{z}\hat{a}^{\dagger}\hat{a}$
where $\hat{a}$ corresponds to a bosonic mode and $\hat{\sigma}_{z}$
corresponds to a transmon qubit. Evolution under this Hamiltonian for a time
$t={\pi}/{\chi}$ realizes the unitary $\hat{U}=I\otimes|g\rangle\langle
g|+e^{i\hat{a}_{2}^{\dagger}\hat{a}_{2}\pi}|e\rangle\langle e|$ in which is a
controlled parity gate. As pictured in Fig. 39, when the controlled parity
gate is placed between an initialization in $|+\rangle$ and measurement of the
transmon in the $X$-basis this realizes a QND measurement of the readout mode
parity. Importantly our use of repeated measurement suppresses the effect of
transmon error mechanisms (i.e. transmon readout error, transmon T1…) on the
final readout fidelity. The readout is still QND in the presence of such
errors as they commute with the dispersive coupling.

We have performed simulations of this measurement scheme to determine the
rough infidelities for different single phonon loss rates
$\kappa_{1}=\kappa_{2}*(\kappa_{1}/\kappa_{2})$. We start by performing master
equation evolution under the deflation and swap to determine $P(\text{even})$
($P(\text{odd})$) which is the probability that the final state after the
deflation and swap steps is even (odd). We include the effects of single
phonon loss during the deflation and swap. Then we sample from these
probabilities to determine the state of the readout mode after the first
measurement. After the measurement the readout mode is in the
$|\hat{n}=0\rangle$ and $|\hat{n}=1\rangle$ manifold which is a good
approximation since after the deflation step the population of the readout
mode is very close to 0 or 1 depending on the initial state. Starting from the
state the readout mode is projected onto after the first measurement we
perform master equation evolution to include the effects of the single phonon
loss, gain, and dephasing on the readout mode during the inter measurement
period and during measurement
($T_{\mathrm{entangle}}+T_{\mathrm{meas}}+T_{\mathrm{reset}}$). We then repeat
this projection and evolution for the remaining number of measurements that
are used, giving us one sequence of projections. To include the effect of
transmon errors such as loss, dephasing, and incorrect measurement we add
additional randomness associated with a fixed transmon error probability
($\epsilon_{q}$) giving us the final measurement sequence. We have performed
Monte Carlo sampling of these measurement sequences to determine the
infidelities of this measurement process. A plot of the infidelities are
pictured in Fig. 40. The assumed numbers in the simulations are listed at the
end of this section.

In other circumstances, such as stabilization of four component cat codes,
decay during entanglement with the transmon is problematic because it induces
dephasing of the cavity. In our case since we are only concerned about
measuring the parity this dephasing is not important. As a result we are
justified in lumping the effect of this transmon decay into our fixed
parameter representing the transmon infidelity mechanisms. We also note that
recent advances in transmon measurement would allow more aggressive transmon
measurement fidelities than what we assume [66]. This would allow us to use
fewer measurements to achieve the same fidelities we currently expect to
achieve.

We can also get an approximate form for the measurement infidelities to expect
with this repeated measurement procedure. Defining $N$ (odd) to be the total
number of measurements and $k\equiv(N+1)/2$, to leading order the error
probability for the majority voting of the repeated measurements for initial
even and odd cat states in the case of no gain are

$\displaystyle\epsilon_{\text{even}}$ $\displaystyle=\epsilon_{\text{(deflate
+ SWAP)}}+\binom{N}{k}\epsilon_{q}^{k}(1-\epsilon_{q})^{N-k}$
$\displaystyle\epsilon_{\text{odd}}$ $\displaystyle=\epsilon_{\text{(deflate +
SWAP)}}+\binom{N}{k}\epsilon_{q}^{k}(1-\epsilon_{q})^{N-k}+\kappa_{1}\text{T}_{p}$
(265)

where $\text{T}_{p}$ is the amount of time after the SWAP and before the kth
measurement and $\epsilon_{\text{(deflate + SWAP)}}$ is the error from the
deflation and SWAP for the given initial state. In the above expressions the
first term is the contribution to the error from the deflation and SWAP steps.
The second term is due to transmon error where k measurements are incorrect.
The last term in the case of an odd initial state is the probability of a
$T_{1}$ event before the kth measurement which will lead with high probability
to all the remaining measurements giving 0. Note that this is the reason that
for an odd initial state and larger $\kappa_{1}/\kappa_{2}$ values majority
voting 5 measurements underperforms majority voting 3 measurements. Using a
more advanced procedure than majority voting would mitigate this problem.

Assumptions–The properties of the measurement were chosen to make the error
probabilities dependent only on $\kappa_{1}/\kappa_{2}$ and the measurement
times scale as $1/\kappa_{2}$ to follow the gates which have the same
dependence. This is convenient for the error correction simulations because it
means that the logical failure rates are independent of the absolute scale of
$\kappa_{2}$.

The properties that feed into our measurement error probabilities and times
are a transmon related error probability ($\epsilon_{q}$) of $1\%$, a
deflation time of $3/{\kappa_{2}}$, transmon related times of
$\text{T}_{\text{entangle}}=\text{T}_{\text{readout}}=\text{T}_{\text{reset}}=200\mathrm{ns}*(1/((100\
\mathrm{ns})*\kappa_{2})$, $\alpha^{2}=8$, and $g/{2\pi}=1\mathrm{MHz}*((100\
\mathrm{ns})*\kappa_{2})$[111, 112]. For simplicity, in the error correction
simulations we have used up to 3 measurements in the repetition code and up to
5 measurements in the surface code. This choice is well justified in the
important regime for REGIME 3. For example, for the surface code the duration
of 4 CNOT gates is $\sim 28\mu s$, larger than the relevant measurement
duration of
$5*T_{\mathrm{entangle}}+5*T_{\mathrm{measure}}+5*T_{\mathrm{reset}}\sim
17.05\mu s$. In REGIME 1 and REGIME 2 this simplification breaks down and
fewer measurements or more time for the measurements would be needed. However
in these regimes fault tolerant quantum computation is infeasible even with
this optimistic measurement error model. In the overhead calculations the
idling duration corresponding to the $X$-basis measurement is
$3/\kappa_{2}+\pi/2g$.

Figure 39: Circuit used for an $X$-basis measurement in the context of an
$X$-type stabilizer measurement. The first step consists of entangling the
ancilla qubit with the data qubits. Afterwords, the ancilla qubit is deflated
followed by a SWAP with a readout mode. Lastly, the readout mode is repeatedly
measured using a transmon qubit. The duration’s for the parts of the
measurement procedure are labeled at the bottom of the figure below each
circuit element. While these repeated parity measurements are occurring, the
CNOT gates of the next error correction cycle can begin. Also included is a
diagram of the physical layout of the stabilizer to give context to the
measurement circuit. Figure 40: Simulated (markers) and analytical (curves)
error probabilities for $X$ measurement for the cases of both even parity and
odd parity initial states. We have taken the parity measurement to be QND and
projective. The dependence on $\kappa_{1}/{\kappa_{2}}$ is stronger for the
case of an odd initial state since the cavity is mapped to $|1\rangle$ after
the deflation. The plotted curve is the leading order analytic model for the
case of $\kappa_{\phi}=0$ Eq. 265. Numerical imprecision predominantly due to
the deflation simulation can have about a 10 percent effect on the simulated
infidelities for the smaller $\kappa_{1}/\kappa_{2}$.

###  G.2 $Z$-basis measurement

For $Z$ measurement we use a beamsplitter interaction between the buffer mode
and a phononic storage mode

$\displaystyle\hat{H}_{r}=g_{r}(\hat{a}^{\dagger}\hat{b}+\hat{b}^{\dagger}\hat{a})$
(266)

where $\hat{a}$ is an annihilation operator on a storage mode and $\hat{b}$ is
an annihilation operator on the buffer mode. By homodyning the output of the
buffer mode we determine the state of the storage mode. We perform this
readout scheme with the two-phonon dissipation off. A similar scheme has been
realized for Kerr-Cat qubits in [27].

To realize this interaction we drive the two terms
$\hat{a}^{\dagger}\hat{b}^{2}$ and $\hat{b}$ off-resonantly. In particular, we
use a pump $\epsilon_{p}(t)=\epsilon_{p}\cos(\omega_{p}t)$ with frequency
$\omega_{p}=2\omega_{b}-\omega_{a}+\Delta$ to off-resonantly drive the term
$\hat{a}\hat{b}^{\dagger 2}$ and directly drive the $\hat{b}$ mode at
frequency $\omega_{d}=\omega_{b}+\Delta$ to produce the term

$\displaystyle\hat{H}_{d}$
$\displaystyle=\epsilon_{d}(\hat{b}^{\dagger}e^{-i\omega_{d}t}+\textrm{h.c.}).$
(267)

Then the complete Hamiltonian in the rotating frame of all of the modes is

$\displaystyle\hat{H}_{\textrm{rot}}$
$\displaystyle=\frac{1}{2}E_{J}\epsilon_{p}\varphi_{a}\varphi_{b}^{2}\hat{a}^{\dagger}\hat{b}^{2}e^{i\Delta
t}+\textrm{h.c.}$ $\displaystyle\quad+\epsilon_{d}\hat{b}^{\dagger}e^{-i\Delta
t}+\textrm{h.c.}+\hat{H}^{\prime}$ (268)

where $\hat{H}^{\prime}$ contains rapidly rotating terms. Now let
$\chi_{1}\equiv E_{J}\epsilon_{p}\varphi_{a}\varphi_{b}^{2}/2$. Then,
neglecting $\hat{H}^{\prime}$ and constants, average Hamiltonian theory yields
[56, 55]

$\displaystyle\hat{H}_{\textrm{eff}}$
$\displaystyle=\frac{1}{\Delta}[\chi_{1}\hat{a}^{\dagger}\hat{b}^{2}+\epsilon_{d}\hat{b},\chi_{1}\hat{a}\hat{b}^{\dagger
2}+\epsilon_{d}\hat{b}^{\dagger}]$
$\displaystyle=\frac{1}{\Delta}\Big{[}\chi_{1}^{2}[2(1+2\hat{b}^{\dagger}\hat{b})\hat{a}^{\dagger}\hat{a}-\hat{b}^{\dagger
2}\hat{b}^{2}]+2\chi_{1}\epsilon_{d}(\hat{a}^{\dagger}\hat{b}+\hat{a}\hat{b}^{\dagger})\Big{]}$
$\displaystyle\xrightarrow{\hat{b}^{\dagger}\hat{b}\ll
1}g(\hat{a}^{\dagger}\hat{b}+\hat{a}\hat{b}^{\dagger})+g_{b}\hat{a}^{\dagger}\hat{a}.$
(269)

The coupling constant is given by
$g=E_{J}\epsilon_{p}\beta\varphi_{a}\varphi_{b}^{2}$ where
$\beta=\epsilon_{d}/\Delta$ and there is an energy shift. The strength of the
coupling is on the order of $g_{2}$ since it depends twice on
$\varphi_{b}>\varphi_{a}$ and $\beta<1$ to ensure
$\hat{b}^{\dagger}\hat{b}<1$.

(a) Error vs. Time (b) Error Fit

Figure 41: a.) Measurement error probability for $Z$-basis readout as a
function of time. The colored lines correspond to the analytic formula for the
separation error Eq. 277 for $\alpha^{2}=1,2,...,10$. The black lines
correspond to infidelities from simulations of the corresponding stochastic
master equation Eq. 278 (QuTiP) in the interaction picture for a few thousand
trajectories with the initial condition $|\alpha\rangle$ for the cases
$\alpha^{2}=1,\alpha^{2}=3,\alpha^{2}=5$, and $\alpha^{2}=7$. The simulated
curved which include $\kappa_{1}/2\pi=1\ \mathrm{KHz}$ and analytic curves
agree well indicating the small effect of the additional single phonon loss.
In order to get concrete time numbers the simulations use $1/\kappa_{2}=100$
ns but as discussed this can be scaled to any $\kappa_{2}$. b.) Plot of the
minimum infidelities vs. $\alpha^{2}$ and the fit line.

Next we proceed by finding the measurement SNR for this readout scheme. The
coupled Langevin equations governing the evolution of the storage and buffer
modes in the interaction picture are

$\displaystyle\dot{\hat{a}}$
$\displaystyle=-i[\hat{a},\hat{H}_{r}]=-ig_{r}\hat{b},$
$\displaystyle\dot{\hat{b}}$
$\displaystyle=-i[\hat{b},\hat{H}_{r}]-\frac{\kappa_{b}}{2}\hat{b}-\sqrt{\kappa_{b}}\hat{b}_{in}=-ig_{r}\hat{a}-\frac{\kappa_{b}}{2}\hat{b}-\sqrt{\kappa_{b}}\hat{b}_{in}$
(270)

Here $\kappa_{b}$ is the single photon loss rate of the buffer mode and we
have neglected the single phonon loss rate of the storage mode under the
assumption that it is far slower than the readout timescale. These equations
can be straightforwardly integrated [113] to give

$\displaystyle\hat{a}(t)$
$\displaystyle=\frac{\hat{a}(0)}{\beta}e^{-\frac{\kappa_{b}t}{4}}(\beta\cosh{\frac{\beta
t}{4}}+\kappa_{b}\sinh{\frac{\beta t}{4}}),$ $\displaystyle\hat{b}(t)$
$\displaystyle=-i\frac{4g\hat{a}(0)}{\beta}e^{-\kappa_{b}t/4}\sinh{\frac{\beta
t}{4}}$ (271)

where $\beta=\sqrt{\kappa_{b}^{2}-(4g_{r})^{2}}$. Here we have not included
the mean zero terms with $\hat{b}_{\text{in}}$ since they are not relevant for
computing the signal. The measurement operator with a uniform readout window
is defined to be [114, 115]

$\displaystyle\hat{M}(\tau)=\sqrt{\kappa_{b}}\int_{0}^{\tau}dt[\hat{b}_{out}^{\dagger}(t)e^{i\phi_{h}}+\hat{b}_{out}(t)e^{-i\phi_{h}}]$
(272)

Using the input-output boundary condition that
$\hat{b}_{out}=\hat{b}_{in}+\sqrt{\kappa_{b}}\hat{b}$. We can determine the
average of the measurement operator to be

$\displaystyle\langle\hat{M}(\tau)\rangle$ $\displaystyle=$
$\displaystyle\frac{2\kappa_{b}\langle\hat{a}(0)\rangle\sin{\phi_{h}}}{g_{r}}*$
$\displaystyle\left[1-e^{-\kappa_{b}\tau/4}\left[\cosh{\frac{\beta\tau}{4}}+\frac{\kappa_{b}}{\beta}\sinh{\frac{\beta\tau}{4}}\right]\right].$
(273)

We have also taken the input to be the vacuum with the property
$\langle\hat{b}_{in}(t^{\prime})\hat{b}_{in}^{\dagger}(t)\rangle=\delta(t-t^{\prime})$.
There is a weak drive on $\hat{b}$ to realize the Hamiltonian which we neglect
and could be replaced with a flux pump. From the average of the measurement
signal $\langle\hat{M}(\tau)\rangle$ we can determine the measurement SNR
using

$\displaystyle\text{SNR}^{2}=\frac{|\langle\hat{M}\rangle_{\alpha}-\langle\hat{M}\rangle_{-\alpha}|^{2}}{\langle\hat{M}_{N(\alpha)}^{2}\rangle+\langle\hat{M}_{N(-\alpha)}^{2}\rangle}.$
(274)

We take the noise terms to be given by

$\displaystyle\langle\hat{M}_{N(\pm\alpha)}^{2}\rangle=\langle(\hat{M}(\tau)-\langle\hat{M}(\tau)\rangle_{\pm\alpha})\rangle=\kappa_{b}\tau$
(275)

which we have checked in numerics.

Solving for the SNR and optimizing the phase we get

$\displaystyle\text{SNR}_{\alpha}(\tau)$ $\displaystyle=$
$\displaystyle\alpha\sqrt{8\kappa_{b}}\frac{\left[1-e^{-\kappa_{b}\tau/4}\left[\cosh{\frac{\beta\tau}{4}}+\frac{\kappa_{b}}{\beta}\sinh{\frac{\beta\tau}{4}}\right]\right]}{g\sqrt{\tau}}.$
(276)

As is expected since this readout scheme is not QND and does not preserve the
state of the cavity at long times the readout SNR goes as $1/\sqrt{\tau}$ as
we are only integrating noise. The readout separation error, which is the
dominant source of error for this readout scheme, will be given by

$\displaystyle\epsilon_{\text{sep},\alpha}(\tau)=\frac{1}{2}\text{Erfc}(\frac{\text{SNR}_{\alpha}(\tau)}{2}).$
(277)

From these equations we can then determine the fidelity as a function of time
for different alpha for our measurement scheme as is shown in the colored line
in Fig. 41 a.). In Fig. 41 b.) we fit the optimal readout error to an
exponential decay. We used a more conservative relation
$\epsilon=e^{-1.5-0.9|\alpha|^{2}}$ for the error correction simulations.

Figure 42: Distance-5 repetition code, where one round of stabilizer
measurements is performed (green measurements) followed by a direct
measurement of the data qubits (red measurements). The data qubits are the
yellow rectangles, and the ancilla qubits (prepared in $|+\rangle$) are the
grey rectangles. On the left, a measurement error on the third data qubit
occurs during the direct measurement of the data, which is equivalent to
having a $Z$ data qubit error immediately before the measurement (shown on the
right). Such settings illustrate the importance of the STOP algorithm, where
one might have to correct errors prior to applying a non-Clifford gate, and a
round of perfect error correction (which in practice is achieved by directly
measuring the data) cannot be performed. In such settings, a single
measurement error during the last round of stabilizer measurements (green
measurements in the figure) can lead to a logical failure if the syndrome
measurement is repeated a fixed number of times (say $d$) rather than using
the STOP algorithm.

We have simulated trajectories for this measurement procedure using a
stochastic master equation for confirmation. The stochastic master equation
simulations include $\kappa_{1}/2\pi=1$ KHz. The integrated and classified and
measurement results from the stochastic master equation compared to the
analytic expression are pictured in Fig. 41 a.). The stochastic master
equation simulations were performed with the evolution

$\displaystyle\frac{d\hat{\rho}(t)}{dt}=$
$\displaystyle-i[g_{r}(\hat{a}^{\dagger}\hat{b}+\hat{b}^{\dagger}\hat{a},\hat{\rho}(t)]$
$\displaystyle+\kappa_{b}(1-\eta)D[\hat{b}]\hat{\rho}(t)+\kappa_{b}\eta
D[\hat{b}]\hat{\rho}(t)$
$\displaystyle+\kappa_{1}D[\hat{a}]\hat{\rho}(t)+n_{\mathrm{th}}\kappa_{1}D[\hat{a}^{\dagger}]\hat{\rho}(t)+\kappa_{\phi}D[\hat{a}^{\dagger}\hat{a}]\hat{\rho}(t)$
(278)

where the second term on the second line is the detection part of the master
equation that is unraveled. We see exponential suppression of the error
probability with increasing $\alpha^{2}$. In the future we expect to be able
to improve the performance by optimizing the window function for the readout
and using the confidence of the measurement result to feed back and improve
the matching. These advances in addition to the robustness of the EC to larger
measurement errors than those currently assumed would allow us to make looser
assumptions.

Assumptions– The chosen parameters for this work are
$\kappa_{b}/2\pi=20*((100\ \mathrm{ns})*\kappa_{2})$ MHz and $g/2\pi=4*((100\
\mathrm{ns})*\kappa_{2})$ MHz. These rates depend on $\kappa_{2}$ so the
readout time will scale as $1/\kappa_{2}$. We also include the effect of a
non-unity quantum efficiency $\eta=.5$ ($\eta$ is the proportion of the signal
detected). For $\kappa_{2}=2\pi*280e3$ used in the paper the optimal duration
of the measurement is roughly $850\ ns$. Note that $\kappa_{b}$ can be made
larger with the main effect of lengthening the readout time.

##  Appendix H Alternative architecture with ATS-based basis readout

In this appendix, we describe an alternative version of our architecture where
the $X$-basis readout is performed directly using the ATS. As discussed below,
this modification favorably impacts practicality and performance across
several different levels of our architecture. For example, this modification
obviates the need for a transmon in each unit cell, and allows us to reduce
the number of modes coupled to each reservoir from five to four, see Figure
43. As a result, crosstalk is reduced, easing limitations on logical lifetimes
posed by correlated errors. Though such benefits are highly desirable, the
ATS-based X-basis readout that underlies this alternative architecture may be
difficult to perform in practice. Because the practical feasibility of the
readout scheme is more speculative, we do not use the ATS-based readout scheme
in any of the analyses in the main text. Instead, we describe the ATS-based
readout scheme and enumerate its favorable impacts here in this appendix. The
ATS-based readout scheme is presented and analyzed in Section H.1, while
Section H.2 catalogs the beneficial impacts of this modification across
different levels of the architecture.

Figure 43: Alternative architecture. The $X$-basis readout is performed
directly using the ATS [contained within the reservoir (RESVR)], as opposed to
using an ancillary readout mode and transmon qubit, cf. Figure 2. By
eliminating the readout mode, the number of modes per unit cell is reduced
from five to four.

###  H.1 ATS-based $X$-basis measurement scheme

Here we present a modified $X$-basis readout scheme that can be implemented
using only the ATS and the buffer mode. The readout scheme is implemented by
engineering a coupling Hamiltonian of the form

$\displaystyle\hat{H}_{r}=ig_{r}\hat{a}^{\dagger}\hat{a}(\hat{b}^{\dagger}-\hat{b}).$
(279)

Here $\hat{a}$ is the annihilation operator for a storage mode and $\hat{b}$
is the annihilation operator for the buffer mode. Note that this Hamiltonian
is equivalent to the longitudinal readout discussed with transmons [115].
Practically speaking, engineering this Hamiltonian may be challenging because
it rotates at the frequency of the buffer mode. That is, we cannot engineer
this Hamiltonian simply by pumping the system at the resonance frequency of
the buffer mode, because the same pump would also bring corresponding
Hamiltonians for other storage modes on resonance simultaneously (e.g.,
$\hat{a}_{2}^{\dagger}\hat{a}_{2}(\hat{b}+\hat{b}^{\dagger})$).

To circumvent this frequency-selectivity problem, we sketch how a Hamiltonian
of the form (279) could be engineered via off-resonant pumping. Consider the
following Hamiltonian (in the rotating frame),

$\displaystyle\hat{H}$ $\displaystyle=(g_{a}\hat{a}\hat{b}^{\dagger
2}+ig_{b}\hat{a}\hat{b}^{\dagger})e^{-i\Delta t}+h.c.$ (280)

This Hamiltonian can be realized in a number of ways. For example, the
$\hat{a}{{\hat{b}}^{\dagger 2}}$ term can be realized by pumping the ATS at
frequency $2\omega_{b}-\omega_{a}+\Delta$, while the
$\hat{a}\hat{b}^{\dagger}$ term can be realized by synthesizing two odd terms
(as described in the vicinity of Equation 250). For simplicity, we assume
access to Hamiltonians of the form (280) and leave the precise details of how
they are to be engineered to future work. For $g_{a,b}\ll\Delta$, evolution
generated by the Hamiltonian (280) is well described by the corresponding
time-averaged effective Hamiltonian [56],

$\displaystyle\hat{H}_{\text{eff}}$
$\displaystyle=\frac{g_{a}^{2}}{\Delta}[\hat{a}^{\dagger}\hat{b}^{2},\hat{a}{{\hat{b}}^{\dagger
2}}]+\frac{g_{b}^{2}}{\Delta}[\hat{a}^{\dagger}\hat{b},\hat{a}\hat{b}^{\dagger}]+$
$\displaystyle\frac{ig_{a}g_{b}}{\Delta}([\hat{a}^{\dagger}\hat{b}^{2},\hat{a}\hat{b}^{\dagger}]-[\hat{a}^{\dagger}\hat{b},\hat{a}{{\hat{b}}^{\dagger
2}}])$
$\displaystyle=\frac{g_{a}^{2}}{\Delta}(2\hat{a}^{\dagger}\hat{a}(1+2\hat{b}^{\dagger}\hat{b})-{\hat{b}^{\dagger
2}}\hat{b}^{2})+\frac{g_{b}^{2}}{\Delta}(\hat{a}^{\dagger}\hat{a}-\hat{b}^{\dagger}\hat{b})+$
$\displaystyle\frac{ig_{a}g_{b}}{\Delta}(2\hat{a}^{\dagger}\hat{a}-\hat{b}^{\dagger}\hat{b})(\hat{b}-\hat{b}^{\dagger})$
$\displaystyle\xrightarrow{\left<\hat{b}^{\dagger}\hat{b}\right>\ll
1}\frac{2ig_{a}g_{b}}{\Delta}\hat{a}^{\dagger}\hat{a}(\hat{b}-\hat{b}^{\dagger})+\frac{2g_{a}^{2}+g_{b}^{2}}{\Delta}\hat{a}^{\dagger}\hat{a}$
(281)

where in the last line we have neglected many of the terms due to the small
occupation of the buffer mode.

To implement an $X$-basis readout using this Hamiltonian, we first deflate the
storage mode. As described in Appendix G, this deflation is achieved by
abruptly setting the two-phonon dissipator for the storage mode equal to
$\mathcal{D}[\hat{a}^{2}]$, as opposed to
$\mathcal{D}[\hat{a}^{2}-\alpha^{2}2]$, and waiting for a timescale comparable
to $1/\kappa_{2}$. When the storage mode begins in an even (odd) parity state,
this deflation procedure maps it to the $|0\rangle$ ($|1\rangle$) phonon Fock
state. Subsequent to this deflation, we perform a homodyne measurement of the
buffer mode while the system evolves under the Hamiltonian (H.1). In effect,
this Hamiltonian drives the buffer mode conditioned on whether the parity of
the storage mode was initially even or odd.

We compute the fidelity of this homodyne readout, where we aim to distinguish
$|\hat{n}_{a}=0\rangle$ from $|\hat{n}_{a}=1\rangle$. We closely follow the
derivation of SNR in [115]. The Langevin equation for the evolution of the
buffer mode in the interaction picture is

$\displaystyle\dot{\hat{b}}$
$\displaystyle=-i[\hat{b},\hat{H}_{\text{eff}}]-\frac{\kappa_{b}}{2}\hat{b}-\sqrt{\kappa_{b}}\hat{b}_{in}$
$\displaystyle=g_{r}\hat{a}^{\dagger}\hat{a}-\frac{\kappa_{b}}{2}\hat{b}-\sqrt{\kappa_{b}}\hat{b}_{in}$
(282)

where $\hat{b}_{in}$ is the input field and $g_{r}\equiv 2g_{a}g_{b}/\Delta$.
In the following calculations we will neglect the single phonon loss of the
ancilla mode which in this simple case will add an average readout error
probability of roughly $\kappa_{a}t/{4}$. We integrate this equation to get
the expected value of the buffer mode

$\displaystyle\langle\hat{b}(t)\rangle=\frac{2g_{r}}{\kappa_{b}}\langle\hat{a}^{\dagger}\hat{a}\rangle(1-e^{-\frac{\kappa_{b}t}{2}}).$
(283)

The measurement operator for integration up to time $\tau$ and with homodyne
angle $\phi_{h}$ is defined as [114, 115]

$\displaystyle\hat{M}(\tau)=\sqrt{\kappa_{b}}\int_{0}^{\tau}dt[\hat{b}_{out}^{\dagger}(t)e^{i\phi_{h}}+\hat{b}_{out}(t)e^{-i\phi_{h}}].$
(284)

Evaluating the average of this integral with the optimal phase gives

$\displaystyle\langle\hat{M}(t)\rangle=\frac{4g_{r}\langle\hat{a}^{\dagger}\hat{a}\rangle}{\kappa_{b}}(-2+2e^{-\frac{\kappa_{b}t}{2}}+\kappa_{b}t).$
(285)

Here we have used the standard input-output condition that
$\hat{b}_{out}=\hat{b}_{in}+\sqrt{\kappa_{b}}\hat{b}$ and the conditions on
$\hat{b}_{in}$ that $\langle\hat{b}_{in}\rangle=0$ and
$\langle\hat{b}_{in}(t)\hat{b}_{in}^{\dagger}(t^{\prime})\rangle=\delta(t-t^{\prime})$.
There is a drive on $\hat{b}$ to realize the Hamiltonian which we neglect and
could be replaced by an appropriate pump. Next we compute the SNR which is
defined as

$\displaystyle\text{SNR}^{2}=\frac{|\langle\hat{M}\rangle_{1}-\langle\hat{M}\rangle_{0}|^{2}}{\langle\hat{M}_{N(1)}^{2}\rangle+\langle\hat{M}_{N(0)}^{2}\rangle}.$
(286)

where $\hat{M}_{N(x)}=\hat{M}-\langle\hat{M}\rangle_{x}$ so
$\langle\hat{M}_{N(0)}^{2}\rangle=\langle\hat{M}_{N(1)}^{2}\rangle=\kappa_{b}t$.
Thus the SNR is

$\displaystyle\text{SNR}(\tau)=\frac{4g_{r}}{\kappa_{b}\sqrt{2\kappa_{b}t}}(-2+2e^{-\frac{\kappa_{b}t}{2}}+\kappa_{b}t).$
(287)

The separation error for this readout, which will be in addition to the effect
of the single phonon loss mentioned earlier, will be given by [67]

$\displaystyle\epsilon_{\text{sep}}(\tau)=\frac{1}{2}\text{Erfc}(\frac{\text{SNR}(\tau)}{2}).$
(288)

There will also be an additional contribution to the readout error associated
with the deflation procedure (see Appendix G)

###  H.2 Benefits of ATS-based X-basis readout

The analysis of the previous section demonstrates that, in principle, high-
fidelity $X$-basis readout can be performed directly using the ATS. We now
show how this readout scheme can be exploited to improve the practicality and
performance of our architecture. The following analysis demonstrates that
improved $X$-basis readout represents one potential path to improved future
designs. At the same time, it also illustrates a more general point:
improvements at low levels of the architecture can propagate into significant
savings at higher levels. These findings thus underscore the importance and
utility of comprehensive architectural analyses.

At the hardware level, the modified $X$-basis readout scheme enables us to
simplify the unit cell of our architecture. Specifically, in comparison to the
unit cell described in the main text (Figure 2), the unit cell of Figure 43
has no transmon and has four phononic modes instead of five. Removing the
transmon reduces the number of control lines and is helpful for device layout
because the transmon requires significant additional space in comparison to
the phononic modes

Reducing the number of modes per unit cell from five to four significantly
reduces crosstalk within each unit cell. We can quantify this reduction
following the approach described in Appendix B. To briefly summarize: we
optimize the frequencies of the storage modes within each unit cell in order
to minimize the effects of crosstalk. More specifically, the frequencies are
chosen so that the effects of coherent _crosstalk_ errors are minimized
subject to the constraint that all _incoherent_ crosstalk errors must be
exponentially suppressed by the filter. The impacts of the residual coherent
errors are quantified via two parameters, $p_{\text{double}}$ and
$p_{\text{triple}}$, that respectively describe the probabilities that either
two data qubits, or two qubits and an ancilla qubit, suffer correlated $Z$
errors.

In Table 10, we compare the results of this mode-frequency optimization for
four- and five-mode unit cells. There are two important takeaways from these
optimization results. The first takeaway is simply that the magnitude of the
residual coherent errors is substantially reduced. Indeed, the total
correlated error probability, $p_{\text{double}}+p_{\text{triple}}$, is over
an order of magnitude smaller for the four-mode case in than the five-mode
case. As we discuss below, this reduction can significantly impact logical
error rates.

# modes | $4J$ | $\omega_{\alpha},\omega_{\beta},\omega_{\gamma},\omega_{\delta},\omega_{\rho}$ | $p_{\text{double}}$ | $p_{\text{triple}}$ | $p_{\text{double}}+p_{\text{triple}}$
---|---|---|---|---|---
4 | 180 | 0, 1000, 798, 101, - | $1.22*10^{-9}\left[\frac{|\alpha|^{2}g_{2}}{2\pi\text{MHz}}\right]^{4}$ | $3.87*10^{-10}\left[\frac{|\alpha|^{2}g_{2}}{2\pi\text{MHz}}\right]^{4}$ | $1.60*10^{-9}\left[\frac{|\alpha|^{2}g_{2}}{2\pi\text{MHz}}\right]^{4}$
5 | 100 | 0, 1000, 242, 879, 61 | $1.83*10^{-8}\left[\frac{|\alpha|^{2}g_{2}}{2\pi\text{MHz}}\right]^{4}$ | $5.20*10^{-10}\left[\frac{|\alpha|^{2}g_{2}}{2\pi\text{MHz}}\right]^{4}$ | $1.88*10^{-8}\left[\frac{|\alpha|^{2}g_{2}}{2\pi\text{MHz}}\right]^{4}$
Table 10: Frequency optimization results. The parameters $4J$ and $\omega$ are
given in units of $2\pi\times$ MHz. The correlated error probabilities in the
last three columns are expressed in terms of $\alpha$ and $g_{2}$. For
realistic choices of $|\alpha|=\sqrt{8}$ and $g_{2}/2\pi=2$ MHz, the
probabilities evaluate to $C=1.05*10^{-4}$ and $C=1.23*10^{-3}$ for the four-
and five-mode configurations respectively. See Appendix B for further details.

The second important takeaway from Table 10 is that the four-mode unit cell
can accommodate a larger filter bandwidth, which ultimately eases requirements
on the lifetimes of the phononic resonators to achieve a desired value of
$\kappa_{1}/\kappa_{2}$. As described in Appendix B, the constraint that all
incoherent crosstalk errors be suppressed by the filter places an upper bound
on the filter bandwidth, and this upper bound decreases as the number of modes
in the unit cell increases. Reducing the number of modes from five to four
thus enables us to increase the filter bandwidth. In particular, we find that
the bandwidth can be increased by nearly a factor of 2 (from $4J=2\pi*100$ MHz
to $4J=2\pi*180$ MHz). As described in Section II, the filter bandwidth
constrains the achievable value of $\kappa_{2}$, meaning that an increased
bandwidth enables an increased $\kappa_{2}$. Equivalently, an increased filter
bandwidth eases the coherence requirement on the phononic resonators to
achieve a desired value of $\kappa_{1}/\kappa_{2}$. From Figure 3, we see that
the factor of $\sim 2$ increase in bandwidth correspondingly lowers the
$T_{1}$ requirement on the phononic resonators by roughly a factor of 2.

Having described the effects of the ATS-based $X$-basis readout scheme at the
hardware level, we now consider effects at the logical level. We begin by
quantifying the impact of the increased measurement error rate of the ATS-
based readout scheme (see Table 3) on the logical error rates. In Fig. 44, we
plot the logical $Z$ error rate of the thin rotated surface code, both for the
main-text architecture of (five-mode unit cell, transmon-based readout) and
for the alternative architecture of this appendix (four-mode unit cell, ATS-
based readout). For the former architecture, we used the same data as shown in
Fig. 8, while for the latter architecture (labeled high measurement), we fixed
the measurement error rate to $2*10^{-3}$. As can be seen, even though the
measurement error rate can be more than an order of magnitude larger for the
ATS-based readout, the logical $Z$ error rates increase by a small amount in
the low $\kappa_{1}/\kappa_{2}$ regime. The reason the logical failure rate is
not greatly affected by the large increase in measurement failure rates is
that CNOT failures are the dominant source of noise. As such, we do not expect
the overhead results of Section VIII to increase when using the alternative
architecture of Figure 43 because the same code distances can be used for
implementing the algorithms of interest.

Figure 44:  Logical $Z$ error rates of the thin rotated surface code with five
modes coupled to an ATS (which includes a transmon qubit and an additional
readout mode in each unit cell and corresponds to the data in Fig. 8) and four
modes coupled to an ATS (which excludes the transmon qubit and performs direct
$X$-basis measurements). We labeled measurements with the transmon qubit as
low measurement and the direct $X$-basis measurement as high measurement. For
the direct $X$-basis measurement, the measurement error rate is fixed at
$2*10^{-3}$ for all values of $\kappa_{1}/\kappa_{2}$. Measurement error rates
with the transmon qubit were obtained from Table 3 with five parity
measurements.

Finally, we consider how the reduction in crosstalk (Table 10) impacts logical
error rates. Logical $Z$ error rates of the thin rotated surface code for the
main-text architecture and the alternative architecture are plotted in Figure
45 (a) and (b), respectively, where correlated errors are simulated with
probabilities given in Table 10. The results reveal that the effects of
crosstalk on logical errors are less severe in the alternative architecture.
In particular, the main-text architecture can accommodate values of
$g_{2}/(2\pi)$ up to $\sim 2$ MHz before the impacts of correlated errors
become significant, while the alternative architecture can accommodate values
twice as large $g_{2}\sim 4$ MHz before the impacts of correlated errors
become significant. This increase in the maximum allowable $g_{2}$ can enable
larger stabilization rates $\kappa_{2}$, or equivalently can reduce the
requirement on the phonon-mode lifetime required to reach a given
$\kappa_{1}/\kappa_{2}$.

(a)  (b)

Figure 45: Logical $Z$ failure rates for a $d_{x}=3$ and $d_{z}=11$ thin
surface code in the presence of the residual crosstalk errors (given in Table
10). The $X$-basis measurement error rates are obtained from Table 3 with five
parity measurements. We compute the logical $Z$ error rates for different
values of $g_{2}$ shown in the legend, and compare such results to the case
where the crosstalk errors are not present. The results in panel (a) [resp.
(b)] are obtained with crosstalk error rates taken from the five-mode (four-
mode) row of Table 10 with $|\alpha|^{2}=8$.

##  Appendix I STOP algorithm

When performing physical non-Clifford operations in between rounds of error
correction (EC), in order to maintain the full effective code distance, it is
crucial to use a fault-tolerant error correction protocol which satisfies the
following definition (taken from [68, 116]):

######  Definition 1.

Fault-tolerant error correction

For $t=\lfloor(d-1)/2\rfloor$, an error correction protocol using a
distance-$d$ stabilizer code $C$ is $t$-fault-tolerant if the following two
conditions are satisfied:

  1. 1.

For an input codeword with error of weight $s_{1}$, if $s_{2}$ faults occur
during the protocol with $s_{1}+s_{2}\leq t$, ideally decoding the output
state gives the same codeword as ideally decoding the input state.

  2. 2.

For $s$ faults during the protocol with $s\leq t$, no matter how many errors
are present in the input state, the output state differs from a codeword by an
error of at most weight $s$.

Apart from being useful for proving thresholds of fault-tolerant error
correction schemes based on code concatenation [117], such a definition of
fault-tolerant error correction is also relevant when performing physical non-
Clifford operations on encoded qubits before directly measuring the data
qubits. In particular, if one implements a minimum weight perfect matching
(MWPM) decoder (see Ref. [118]) with $\mathcal{O}(d)$ rounds of stabilizer
measurements (where $d$ is the code distance of the error correcting code
protecting the data), a measurement error in the last round can lead to a
logical failure and Definition 1 would not be satisfied. In many fault-
tolerant implementations, such a problem can be avoided by implementing non-
Clifford operations via gate injection and stabilizer operations, followed by
direct measurement of the data qubits (hence physical non-Clifford gates are
never directly applied to the data qubits). An example of direct measurements
of the data qubits after performing one round of stabilizer measurements for a
$d=5$ repetition code is given in Fig. 42. By measuring the data qubits,
measurement errors can be treated as data qubit errors arising prior to
performing the measurement [72]. As such, measuring the data directly acts as
a round of perfect error correction.

As was shown in Section VI, post-selection can be avoided when preparing the
logical $|0\rangle_{L}$ and $|1\rangle_{L}$ states (used to obtain the state
$|\psi\rangle_{1}$) if we have a decoder that is robust to measurement errors
in the last syndrome measurement round prior to applying the physical Toffoli
gates. For the BUTOF protocol, we cannot directly measure the data prior to
applying the physical Toffoli gates. Using ideas from Ref. [116], in this
section we propose an algorithm which tells us when to terminate the sequence
of error syndrome measurements, which we call the STOP algorithm, and which
satisfies Definition 1 when using the syndrome measurement from the last round
to correct errors. Further, in Appendix J, we show how the STOP algorithm can
be used with gate injection to perform all stabilizer operations of the
repetition code.

Figure 46: Example of a single controlled-$Z$ failure resulting in the error
$Z_{n}\otimes Z_{n+1}$ (where $Z_{n+1}$ acts on the ancilla qubit) when
measuring the operator $Z^{\otimes n}$. Here $n$ is the number of data qubits.
This single fault can cause three consecutive syndrome measurements to yield
three distinct outcomes. Here $E_{\text{in}}$ is an input error with syndrome
$s(E_{\text{in}})=s_{1}$.

The goal of the STOP algorithm is to track consecutive syndrome outcomes
$s_{1},s_{2},\cdots,s_{r}$ and to compute the minimum number of faults which
could have caused this sequence of syndromes. In particular, let
$n_{\text{diff}}$ be a counter which tracks the minimum number of faults
causing changes in syndrome outcomes, and consider the consecutive syndromes
$s_{k-1},s_{k}$ and $s_{k+1}$. Given that a single fault can lead to two
syndrome changes as in the example below, suppose we obtain different
syndromes in rounds $k$ and $k+1$ (so that $s_{k}\neq s_{k+1}$). In order to
decide whether to increment $n_{\text{diff}}$ by one, we must first check
whether $n_{\text{diff}}$ was incremented after measuring the $k$’th error
syndrome. If $n_{\text{diff}}$ didn’t increase after the $k$’th round, then we
increment $n_{\text{diff}}$ by one. Otherwise, $n_{\text{diff}}$ remains
unchanged.

As an example, suppose a single fault occurs during the second round of
stabilizer measurements of an EC protocol adding a weight-one error to the
data qubits while also flipping the measurement outcome of one of the
stabilizers (in this case $Z^{\otimes n}$ as shown in Fig. 46). Further,
suppose the input error to the second round of the EC protocol $E_{\text{in}}$
has the error syndrome $s(E_{\text{in}})=s_{1}$, and that the error
$Z_{n}E_{\text{in}}$ has error syndrome $s_{3}\neq s_{1}$ (here $Z_{n}$ is the
$Z$ error added to the data qubit arising from the two-qubit gate failure).
Since the $Z$ error flipped the measurement outcome of $Z^{\otimes n}$, the
syndrome $s_{2}$ measured during the second round can differ from both $s_{1}$
and $s_{3}$.

Result: Final syndrome $s_{r}$ for $r$ repetitions of the syndrome
measurement.

initialize: $t=(d-1)/2$; $n_{\text{diff}}=0$; $\text{countSyn}=1$;
$\text{SynRep}=1$; $n_{\text{diff}}\text{Increase}=0$; $\text{test}=0$;

while  _$\text{test}=0$_ do

if  _$n_{\text{diff}}=t$_ then

$\text{test}=1$

end if

Measure the error syndrome $s_{j}$. Store the error syndrome $s_{j-1}$ from
the previous round in synPreviousRound and the current syndrome $s_{j}$ in
synCurrentRound.;

if  _$\text{countSyn} >1$_ then

if  _$\text{synPreviousRound}=\text{synCurrentRound}$_ then

$\text{SynRep}=\text{SynRep}+1$;

$n_{\text{diff}}\text{Increase}=0$;

else

$\text{SynRep}=0$;

if  _$n_{\text{diff}}\text{Increase}=0$_ then

$n_{\text{diff}}=n_{\text{diff}}+1$;

$n_{\text{diff}}\text{Increase}=1$;

else

$n_{\text{diff}}\text{Increase}=0$;

end if

end if

end if

if  _$\text{SynRep}=t-n_{\text{diff}}+1$_ then

$\text{test}=1$;

end if

$\text{countSyn}=\text{countSyn}+1$;

end while

Algorithm 1 STOP algorithm

With the above example in mind, the STOP algorithm is given by Algorithm 1. To
see why a decoding algorithm based on Algorithm 1 satisfies Definition 1,
consider the case where the total number of input errors and faults during the
EC is $t=(d-1)/2$ for a distance $d$ error correcting code. If at any time
during the EC the same syndrome $s_{j}$ is measured $t-n_{\text{diff}}+1$
times in a row, then it must have been the correct syndrome (with very high
probability). The reason is that given the value of $n_{\text{diff}}$, which
counts the minimum number of faults compatible with the syndrome history since
the beginning of the current cycle of error correction, there would need to be
more than $t$ faults to cause all $t-n_{\text{diff}}+1$ consecutive syndromes
to be incorrect due to failures resulting in flipped measurement outcomes. As
such one could use the syndrome $s_{j}$ to correct errors and terminate the
protocol. Doing so, there could only be $\leq t$ residual leftover errors that
went undetected in the last measurement round.

Similarly, if after measuring the $r-1$’th syndrome $n_{\text{diff}}$ gets
incremented to $n_{\text{diff}}=t$, then we know that at least $t$ faults must
have occurred during the EC. As such, by repeating the syndrome measurement
one more time (resulting in the syndrome $s_{r}$) and using that syndrome to
decode, there would need to have been more than $t$ faults for $s_{r}$ to be
the wrong syndrome (due to faults flipping some of the stabilizer measurement
outcomes in the last round). Hence using $s_{r}$ to decode would result in
residual errors with weight $v\leq t$ (where, as stated at the beginning of
the previous paragraph, the total number of input errors and faults during the
EC is $t$).

Given the above, we conclude that when using Algorithm 1, the sequence of
syndrome measurements will terminate if one of the following conditions is
satisfied:

  1. 1.

The syndrome $s_{r}$ is repeated $t-n_{\text{diff}}+1$ times in a row.

  2. 2.

The counter $n_{\text{diff}}$ gets incremented to $n_{\text{diff}}=t$. Measure
the syndrome one more time resulting in the syndrome $s_{r}$. Use $s_{r}$ to
decode.

Decoding will succeed if the total number of input errors and faults during
the EC cycle is $\leq t$.

We now provide a few remarks. Firstly, given a particular error correcting
code and decoder along with the STOP algorithm for repeating the syndrome
measurement, one can satisfy Definition 1 by using the last measured syndrome
$s_{r}$ to decode while ignoring the entire syndrome history. Hence in such
settings, one can use a simple code-capacity-type decoder to decode with
$s_{r}$ (i.e. a decoder which ignores measurement and space-time correlated
errors). As an example, one can decode with the surface code using a MWPM or
Union Find decoder (see Ref. [119]) on a two-dimensional graph instead of a
three-dimensional graph tracking the entire syndrome history. Doing so could
significantly reduce the overall decoding time. In general however, the
approach where the STOP algorithm is used to ignore the entire syndrome
history apart the last syndrome $s_{r}$ does not have a threshold 444There are
cases where families of error correcting codes have thresholds when using the
STOP algorithm to decode with the syndrome from the last round. One such
example includes concatenated codes using the methods of [117]. To see this,
consider a distance $d$ repetition code and a stochastic noise model where
fault locations fail with probability $p$. After $d$ rounds of repeating the
syndrome measurement, there will be approximately $pd^{2}$ failures. For large
distances $d$, with high probability the error syndrome will change in every
round. Hence the probability of a measurement error in the final round will
not depend on the past syndrome history and the decoder will fail to correct
the errors with high probability.

On the other hand, tracking the entire syndrome history and using Algorithm 1
to decide when to terminate the rounds of repeated syndrome measurement
generally leads to lower failure rates and has a threshold. Indeed, when
computing the memory failure rates of the repetition code using Algorithm 1
for deciding when to terminate the syndrome measurements), we found that
performing MWPM on the entire syndrome history leads to lower logical failure
rates compared to performing MWPM on a one-dimensional graph using only the
final syndrome $s_{r}$ (note that the logical $Z$ failure rates for the
repetition code in Fig. 6 were computed by applying MWPM to the full syndrome
history of the measured syndromes using the STOP algorithm). As such, the EC
protocols used in this work when considering repetition codes implement MWPM
on the entire syndrome history in conjunction with Algorithm 1 to decide when
to stop measuring the error syndrome.

We conclude this section by providing a lower and upper bound on the maximum
number of syndrome measurement repetitions that can be performed using the
STOP algorithm. For the case where there are no failures, it is
straightforward to see that the syndrome measurement will be repeated $t+1$
times. To find the upper bound, we consider the worst case scenario, where
(starting with $n_{\text{diff}}=0$) there are no failures in the first $t$
rounds of syndrome measurement, so that the same syndrome is repeated $t$
times. However, in round $t+1$, a measurement error occurs and
$n_{\text{diff}}$ gets incremented to $n_{\text{diff}}=1$. Now again, suppose
there are no failures in the next $t-1$ rounds (so the same syndrome is
repeated $t-1$ times) and a measurement error occurs in the $t$’th round, so
that $n_{\text{diff}}$ is incremented to $n_{\text{diff}}=2$. Suppose the same
pattern repeats itself until all $t$ faults are exhausted resulting in
$n_{\text{diff}}=t$. By the protocol of the STOP algorithm, we must repeat the
syndrome measurement one more time. For such a fault pattern, the total number
of syndrome measurements $s_{\text{tot}}$ is then given by

$\displaystyle
s_{\text{tot}}=\sum_{k=0}^{t-1}(t-k)+t+1=\frac{1}{2}(t^{2}+3t+2)=\binom{t+2}{2}.$
(289)

For low code distances and low noise rate regimes, the average number of
repetitions will approach $t+1$. However for large code distances, with high
probability, the syndrome measurement outcome will change every round and thus
$n_{\text{diff}}$ changes every other round. Thus after $2t$ rounds,
$n_{\text{diff}}=t$ and the syndrome must be repeated one more time resulting
in a total number of $2t+1=d$ rounds. It should then be expected that for
large code distances, the performance of MWPM decoders based on a fixed $d$
rounds will perform similarly to a decoder which uses the STOP algorithm to
terminate while implementing MWPM over the full syndrome history.

##  Appendix J Stabilizer operations with the repetition code

In this section we describe how to do all stabilizer operations with the
repetition code. However, the methods we provide apply to any family of
Calderbank-Shor-Steane (CSS) codes.

###  J.1 Computational basis states

We begin by describing how to prepare the logical computational basis states
of the repetition code. Doing so, we provide two schemes for preparing
$|0\rangle_{L}$.

Scheme 1: Using the fact that for an $n$-qubit repetition code
$|+\rangle_{L}=|+\rangle^{\otimes n}$, preparing $|+\rangle^{\otimes n}$
followed by a logical $Z_{L}=Z^{\otimes n}$ measurement (see Fig. 47) projects
the state to $|0\rangle_{L}$ given a $+1$ outcome and $|1\rangle_{L}$ given a
$-1$ outcome. Since a measurement error on the ancilla results in a logical
$X_{L}=X_{1}$ error applied to the data, fault-tolerance can be achieved by
repeating the measurement of $Z_{L}$ using the STOP algorithm (where the
syndrome corresponds to the ancilla measurement outcome) and applying the
appropriate $X_{L}$ correction given the final measurement outcome. For
instance, if $|0\rangle_{L}$ is the desired state and the final measurement
outcome at the termination of the STOP algorithm is $-1$, $X_{1}$ would be
applied to the data. Lastly, note that only $X$ errors can propagate from the
ancilla to the data but these are exponentially suppressed by the cat-qubits.

Figure 47: Circuit for preparing logical computational basis states of the
repetition code. A measurement error results in a logical $X$ error applied to
the data. Fault-tolerance is achieved by repeating the measurement using the
STOP algorithm.

Scheme 2: Here we present a more conventional approach for preparing the
computational basis states which only involves stabilizer measurements (see
for instance Refs. [121, 122, 123]). Starting with the state
$|\psi_{1}\rangle=|0\rangle^{\otimes n}$ which is a $+1$ eigenstate of
$Z_{L}$, measure all stabilizers of the repetition code (each having a random
$\pm 1$ outcome) resulting in the state

$\displaystyle|\psi_{2}\rangle=\prod_{i=1}^{n-1}\Big{(}\frac{I\pm
X_{i}X_{i+1}}{2}\Big{)}|0\rangle^{\otimes n}.$ (290)

If the measurement outcome of $X_{k}X_{k+1}$ is $-1$, the correction
$\prod_{j=1}^{k}Z_{j}$ can be applied to the data to flip the sign back to
$+1$. However given the possibility of measurement errors, the measurement of
all stabilizers $\langle X_{1}X_{2},X_{2}X_{3},\cdots,X_{n-1}X_{n}\rangle$
must be repeated. If physical non-Clifford gates are applied prior to
measuring the data, then the STOP algorithm can be used to determine when to
stop measuring the syndrome. Subsequently, MWPM is applied to the full
syndrome history to correct errors and apply the appropriate $Z$ corrections
to fix the code-space given the initial stabilizer measurements. After
performing numerical simulations, we found that scheme 2 achieves lower
logical failure rates compared to scheme 1. Further, since physical Toffoli
gates are applied to the data qubits in order to prepare a
$|\text{TOF}\rangle$ magic state (see Section VI) and given the constraints
imposed by our ATS architecture (which make performing global $Z$ measurements
very challenging using a single ancilla qubit), we always use scheme 2 along
with the STOP algorithm when preparing logical computational basis states.

Lastly, we remark that although the logical component of an uncorrectable
error $E^{(z)}Z_{L}$ (where $E^{(z)}$ is correctable) can always be absorbed
by $|0\rangle_{L}$ resulting in an output state
$|\psi_{\text{out}}\rangle=E^{(z)}|0\rangle_{L}$, it is still important to
have a fault-tolerant preparation scheme for $|0\rangle_{L}$ (and thus to
repeat the measurement of all stabilizers enough times). For instance, if a
single fault results in a weight-two correctable $Z$ error (assuming $n\geq
5$), a second failure adding one or more data qubit errors during a subsequent
part of the computation can combine with the weight-two error resulting in an
uncorrectable data qubit error. Hence, such a preparation protocol would not
be fault-tolerant up to the full code distance.

###  J.2 Implementation of logical Clifford gates

Figure 48: Circuit for implementing a logical $S$ gate. The circuit requires
the preparation of $|i\rangle_{L}$, and the CNOT gate is transversal. A
logical $|Z\rangle_{L}$ operator is applied when the measurement outcome of
the ancilla is $-1$.

Since the CNOT gate is transversal for the repetition code, we focus on
implementing a generating set of single-qubit Clifford operations. Recall that
the Clifford group is generated by $\mathcal{P}_{n}^{(2)}=\langle
H_{i},S_{i},\text{CNOT}_{ij}\rangle$, where

$\displaystyle
H=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}\begin{array}[]{ccc}1&1\\\
1&-1\end{array}\\\ \end{array}\right),\
S=\left(\begin{array}[]{c}\begin{array}[]{ccc}1&0\\\ 0&i\\\ \end{array}\\\
\end{array}\right),$ (297)

are the Hadamard and phase gate operators. In what follows we show how to
implement $S$ and $Q=SHS$ which also forms a generating set for single-qubit
Clifford operations. A key to the implementation of such gates will be the
injection of the state $|i\rangle=\frac{1}{\sqrt{2}}(|0\rangle+i|1\rangle)$,
which is a $+1$ eigenstate of the Pauli $Y$ operator. The logical state
$|i\rangle_{L}$ can be prepared using scheme 1 of Section J.1 by replacing
$Z_{L}$ with $Y_{L}=Y_{1}Z_{2}\cdots Z_{n}$.

Figure 49: Circuit for implementing a logical $Q=SHS$ gate. The circuit
requires the preparation of $|i\rangle_{L}$, and the CNOT gate is transversal.
A logical $|Y\rangle_{L}$ operator is applied when the measurement outcome of
the ancilla is $-1$.

In Fig. 48 we provide a circuit for implementing $S_{L}$ which requires
$|i\rangle_{L}$ as an input state, a transversal CNOT gate, and a logical
$Z$-basis measurement. If a $-1$ measurement outcome is obtained, we apply a
$Z_{L}$ correction to the data. Note however that a measurement error can
result in a logical $Z_{L}$ being applied incorrectly to the data. As such, to
guarantee fault-tolerance, one can repeat the circuit of Fig. 48 and use the
STOP algorithm to decide when to terminate. The final measurement outcome is
then used to determine if a $Z_{L}$ correction is necessary. The
implementation of $S_{L}$ can thus be summarized as follows:

Figure 50: Efficient circuit for implementing a $CZ$ gate given the higher
cost of logical $H$ gates compared to logical $S$ gates.

$S_{L}$ gate implementation:

  1. 1.

Implement the circuit in Fig. 48 and let the measurement outcome be $s_{1}$.

  2. 2.

Repeat the circuit in Fig. 48 and use the STOP algorithm to decide when to
terminate.

  3. 3.

If the final measurement outcome $s_{r}=+1$, do nothing, otherwise apply
$Z_{L}=Z_{1}Z_{2}\cdots Z_{n}$ to the data.

The circuit for implementing the logical $Q=SHS$ gate is given in Fig. 49. The
circuit consists of an injected $|i\rangle_{L}$ state, a transversal CNOT gate
and a logical $X$-basis measurement is applied to the input data qubits. If
the measurement outcome is $-1$, $Y_{L}$ is applied to the data. As with the
$S_{L}$ gate, we repeat the application of the circuit in Fig. 49 to protect
against measurement errors. The full implementation of $Q_{L}$ is given as
follows:

$Q_{L}$ gate implementation:

  1. 1.

Implement the circuit in Fig. 49 and let the measurement outcome be $s_{1}$.

  2. 2.

Repeat the circuit in Fig. 49 and use the STOP algorithm to decide when to
terminate.

  3. 3.

If the final measurement outcome $s_{r}=+1$, do nothing, otherwise apply
$Y_{L}=Y_{1}Z_{2}\cdots Z_{n}$ to the data.

Note that the logical Hadamard gate can be obtained from the $S_{L}$ and
$Q_{L}$ protocols using the identity
$H=S^{\dagger}SHSS^{\dagger}=S^{\dagger}QS^{\dagger}$. Hence ignoring
repetitions of the circuits in Figs. 48 and 49, the implementation of $H_{L}$
requires three logical CNOT gates, two $|-i\rangle_{L}$ and one
$|i\rangle_{L}$ state, two logical $Z$ basis measurements, and one logical $X$
basis measurement. Instead of using two logical Hadamard gates and one CNOT
gate to obtain a $CZ$ gate, we provide a more efficient circuit in Fig. 50.

Lastly, we point out that since the circuits in Figs. 48 and 49 contain only
stabilizer operations and injected $|i\rangle_{L}$ states, using the STOP
algorithm to repeat the measurements is not strictly necessary. For instance,
one could repeat the measurement a fixed number of times and majority vote
instead of using the STOP algorithm. However in low noise rate regimes, the
STOP algorithm can potentially be much more efficient since the average number
of repetitions for the measurements can approach $t+1$ where $t=(d-1)/2$.

##  Appendix K Growing encoded data qubits to larger code distances with the
repetition code

In this section, we provide a simple protocol for growing a state
$|\overline{\psi}\rangle_{d_{1}}=\alpha|0\rangle_{d_{1}}+\beta|1\rangle_{d_{1}}$
encoded in a distance $d_{1}$ repetition code to a state
$|\overline{\psi}\rangle_{d_{2}}=\alpha|0\rangle_{d_{2}}+\beta|1\rangle_{d_{2}}$
encoded in a distance $d_{2}>d_{1}$ repetition code. We emphasize that the
protocol presented in this section is applicable to arbitrary states and will
be used for growing $|\text{TOF}\rangle$ states prepared using the fault-
tolerant methods of Section VI to larger code distances.

Let $\mathcal{S}_{d_{1}}=\langle
X_{1}X_{2},X_{2}X_{3},\cdots,X_{d_{1}-1}X_{d_{1}}\rangle$ be the stabilizer
group for the distance $d_{1}$ repetition code with cardinality
$|\mathcal{S}_{d_{1}}|=d_{1}-1$. Similarly, we define
$\mathcal{S}_{d^{\prime}_{1}}=\langle
X_{d_{1}+1}X_{d_{1}+2},\cdots,X_{d_{2}-1}X_{d_{2}}\rangle$ with
$|\mathcal{S}_{d^{\prime}_{1}}|=d_{2}-d_{1}-1$. Furthermore, the stabilizer
group for the distance $d_{2}$ repetition code is given by
$\mathcal{S}_{d_{2}}=\langle
X_{1}X_{2},X_{2}X_{3},\cdots,X_{d_{2}-1}X_{d_{2}}\rangle$.

In the remainder of this section, we define $g_{i}^{(d_{1})}$ to be the $i$’th
stabilizer in $\mathcal{S}_{d_{1}}$ and $g_{i}^{(d^{\prime}_{1})}$ to be the
$i$’th stabilizer in $\mathcal{S}_{d^{\prime}_{1}}$, so that
$g_{i}^{(d_{1})}=X_{i}X_{i+1}$ and
$g_{i}^{(d^{\prime}_{1})}=X_{d_{1}+i}X_{d_{1}+i+1}$.

Protocol for growing $|\overline{\psi}\rangle_{d_{1}}$ to
$|\overline{\psi}\rangle_{d_{2}}$:

  1. 1.

Prepare the state $|\psi_{1}\rangle=|0\rangle^{\otimes(d_{2}-d_{1})}$.

  2. 2.

Measure all stabilizers in $\mathcal{S}_{d^{\prime}_{1}}$ resulting in the
state
$|\psi_{2}\rangle_{d^{\prime}_{1}}=\prod_{i=1}^{d_{2}-d_{1}-1}\Big{(}\frac{I\pm
g_{i}^{(d^{\prime}_{1})}}{2}\Big{)}|0\rangle^{\otimes(d_{2}-d_{1})}$.

  3. 3.

Repeat the measurement of stabilizers in $\mathcal{S}_{d^{\prime}_{1}}$ and
apply MWPM to the syndrome history to correct errors and project to the code-
space. If $g_{i}^{(d^{\prime}_{1})}$ is measured as $-1$ in the first round,
apply the correction $\prod_{k=d_{1}+1}^{d_{1}+i}Z_{k}$ to the data.

  4. 4.

Prepare the state
$|\psi_{3}\rangle=|\psi\rangle_{d_{1}}\otimes|\psi_{2}\rangle_{d^{\prime}_{1}}$
and measure $X_{d_{1}}X_{d_{1}+1}$.

  5. 5.

Repeat the measurement of all stabilizers if $\mathcal{S}_{d_{2}}$ and use
MWPM over the syndrome history to correct errors. If in the first round the
stabilizer $X_{d_{1}}X_{d_{1}+1}$ is measured as $-1$, apply the correction
$\prod_{i=1}^{d_{1}}Z_{i}$.

As remark, the corrections stated in step 3 and 5 can be postponed to a later
time after the growing protocol is completed. The reason is that one can use
the entire syndrome history from each step, in addition to the syndromes
measured after the states have merged to apply the appropriate corrections.

Figure 51: Diagram illustrating our protocol for growing the state
$|\psi\rangle_{d_{1}}$ to $|\psi\rangle_{d_{2}}$ with the ATS layout by
starting with the two blocks stabilized by $\mathcal{S}_{d_{1}}$ and
$\mathcal{S}_{d^{\prime}_{1}}$ . The yellow vertices are the data qubits, and
the gray vertices correspond to the ancilla qubits used to measure the
stabilizers of the repetition code. The measurement of $X_{d_{1}}X_{d_{1}+1}$
(with random $\pm 1$ outcome) is highlighted by the purple semi-circle. After
performing the appropriate corrections, the final block is stabilized by
$\mathcal{S}_{d_{2}}$.

The growing scheme involves two blocks, the first being the state
$|\overline{\psi}\rangle_{d_{1}}$ which we want to grow to
$|\overline{\psi}\rangle_{d_{2}}$. The second block involves the set of qubits
which are prepared in the state $|\psi_{2}\rangle_{d^{\prime}_{1}}$ and
stabilized by $\mathcal{S}_{d^{\prime}_{1}}$ (steps 1-3). The key is to
measure the boundary operator $X_{d_{1}}X_{d_{1}+1}$ between the two blocks,
which effectively merges both blocks into the encoded state
$|\overline{\psi}\rangle_{d_{2}}$ and constitutes a simple implementation of
lattice surgery [76, 124, 78, 125]. To see this, consider the state prior to
step 4:

$\displaystyle|\psi_{3}\rangle$
$\displaystyle=|\overline{\psi}\rangle_{d_{1}}\otimes|\psi_{2}\rangle_{d^{\prime}_{1}}$
$\displaystyle=\alpha|0\rangle_{d_{1}}\otimes|\psi_{2}\rangle_{d^{\prime}_{1}}+\beta|1\rangle_{d_{1}}\otimes|\psi_{2}\rangle_{d^{\prime}_{1}}$
$\displaystyle=\alpha\prod_{i=d_{1}+1}^{d_{2}-1}\Big{(}\frac{I+g_{i}^{(d^{\prime}_{1})}}{2}\Big{)}|0\rangle_{d_{1}}\otimes|0\rangle^{\otimes(d_{2}-d_{1})}$
$\displaystyle+\beta
X_{1}\prod_{i=d_{1}+1}^{d_{2}-1}\Big{(}\frac{I+g_{i}^{(d^{\prime}_{1})}}{2}\Big{)}|0\rangle_{d_{1}}\otimes|0\rangle^{\otimes(d_{2}-d_{1})},$
(298)

where we used $|1\rangle_{d_{1}}=X_{1}|0\rangle_{d_{1}}$. When measuring
$X_{d_{1}}X_{d_{1}+1}$ and performing the correction
$\prod_{i=1}^{d_{1}}Z_{i}$ if the measurement outcome is $-1$,
$|\psi\rangle_{3}$ is projected to

$\displaystyle|\psi\rangle_{f}$
$\displaystyle=\alpha\prod_{i=d_{1}+1}^{d_{2}-1}\Big{(}\frac{I+g_{i}^{(d^{\prime}_{1})}}{2}\Big{)}\Big{(}\frac{I+X_{d_{1}}X_{d_{1}+1}}{2}\Big{)}\prod_{j=1}^{d_{1}-1}\Big{(}\frac{I+g_{j}^{(d_{1})}}{2}\Big{)}|0\rangle^{\otimes
d_{2}}$ $\displaystyle+\beta
X_{1}\prod_{i=d_{1}+1}^{d_{2}-1}\Big{(}\frac{I+g_{i}^{(d^{\prime}_{1})}}{2}\Big{)}\Big{(}\frac{I+X_{d_{1}}X_{d_{1}+1}}{2}\Big{)}\prod_{j=1}^{d_{1}-1}\Big{(}\frac{I+g_{j}^{(d_{1})}}{2}\Big{)}|0\rangle^{\otimes
d_{2}}$
$\displaystyle=\alpha\prod_{i=1}^{d_{2}-1}\Big{(}\frac{I+X_{i}X_{i+1}}{2}\Big{)}|0\rangle^{\otimes
d_{2}}$ $\displaystyle+\beta
X_{1}\prod_{i=1}^{d_{2}-1}\Big{(}\frac{I+X_{i}X_{i+1}}{2}\Big{)}|0\rangle^{\otimes
d_{2}}$ $\displaystyle=\alpha|0\rangle_{d_{2}}+\beta X_{1}|0\rangle_{d_{2}}$
$\displaystyle=|\psi\rangle_{d_{2}}.$ (299)

The rounds of repeated stabilizer measurements in steps 3 and 5 are required
due to the random outcomes and measurement errors which can occur when
performing the appropriate projections. A pictorial representation for the
growing scheme is shown in Fig. 51.

##  Appendix L Toffoli simulation twirling approximation

To simulate the fault-tolerant preparation of the $|\text{TOF}\rangle$ state
taking into account all fault-locations, we implement Monte-Carlo methods
using a Gottesman-Knill type simulation [126] to avoid running into
scalability issues. However, since the circuit in Fig. 14b contains physical
Toffoli gates, some type of approximation is necessary to perform a Gottesman-
Knill type simulation with such circuits.

In order to determine the most appropriate type of approximation, writing a
Toffoli gate as CCX, we first observe that

$\displaystyle(\text{CCX})(I\otimes I\otimes Z)|\psi\rangle=(CZ_{A,B}\otimes
Z)(\text{CCX})|\psi\rangle,$ (300)

for some arbitrary state $|\psi\rangle$. In other words, propagating a $Z$
error through the target qubit of the Toffoli gate results in a $CZ$ error on
the two control qubits. Recall that we label the three logical qubits by
$\\{A,B,C\\}$.

In what follows, we will consider the transformation of the
$|\text{TOF}\rangle$ state with input data qubit $Z$ errors on the third block
when using a single $|+\rangle$ ancilla to measure $g_{A}$. Note that all
conclusions remain unchanged if instead we used the GHZ state of Fig. 14b.

Let $A_{k}$ be a subset of $k$ qubits and consider $k\geq 1$ data qubit errors
on the third block expressed as $E^{(C)}=I\otimes I\otimes\prod_{j\in
A_{k}}Z_{j}\equiv\prod_{j=1}^{k}Z^{(C)}_{j}$. We have that

$\displaystyle|\psi\rangle_{\text{in}}=|+\rangle\prod_{j=1}^{k}Z_{j}^{(C)}|\text{TOF}\rangle.$
(301)

After applying $g_{A}$ and propagating the $Z$ errors through the Toffoli
gates, $|\psi\rangle_{\text{in}}$ becomes

$\displaystyle\prod_{j=1}^{k}Z_{j}^{(C)}\Big{(}|0\rangle|\text{TOF}\rangle+|1\rangle\prod_{j=1}^{k}Z_{j}^{(A)}|\text{TOF}\rangle\Big{)}$
$\displaystyle=\prod_{j=1}^{k}Z^{(C)}_{j}\Big{[}|+\rangle\Big{(}\frac{I+\prod_{j=1}^{k}Z^{(A)}_{j}}{\sqrt{2}}\Big{)}|\text{TOF}\rangle$
$\displaystyle+|-\rangle\Big{(}\frac{I-\prod_{j=1}^{k}Z^{(A)}_{j}}{\sqrt{2}}\Big{)}\Big{]},$
(302)

where $\prod_{j=1}^{k}Z_{j}^{(A)}$ are products of $Z$ errors on the first
data block which have identical support with the $Z$ errors on the third
block. After measuring the ancilla in the $X$ basis, a $\pm 1$ measurement
outcome results in the state $|\psi\rangle_{f}$ given by

$\displaystyle|\psi\rangle_{f}=\prod_{j=1}^{k}Z_{j}^{(C)}\Big{(}\frac{I\pm\prod_{j=1}^{k}Z^{(A)}_{j}}{\sqrt{2}}\Big{)}|\text{TOF}\rangle.$
(303)

From Eq. 303, we see that when performing one round of error detection of the
first block $A$, the error
$\Big{(}\frac{I\pm\prod_{j=1}^{k}Z^{(A)}_{j}}{\sqrt{2}}\Big{)}$ will project
either to $I$ or $\prod_{j=1}^{k}Z^{(A)}_{j}$ with $50\%$ probability each
unless $\prod_{j=1}^{k}Z^{(A)}_{j}=Z^{(A)}_{L}$ in which case the state
remains unchanged.

Given the above, when performing our Gottesman-Knill type simulations when
measuring $g_{A}$, if the input $Z$ errors to the third block are
$\prod_{j=1}^{k}Z^{(C)}_{j}$, we flip the GHZ ancilla measurement outcome with
$50\%$ probability and do the following: If $k<d$, we add the $Z$ errors
$\prod_{j=1}^{k}Z^{(A)}_{j}$ to the first block with $50\%$ probability. If
$k=d$, we add $Z_{L}$ to the first block with $100\%$ probability.

Note that such a simulation method is exact when $k<d$ and only introduces a
discrepancy when $k=d$. Since such events are rare, our approximation method
differs from an exact simulation of the bottom-up $|\text{TOF}\rangle$ state
preparation scheme only by a small amount.

##  Appendix M Fitting procedure for memory and lattice surgery

Here we extend the discussion of lattice surgery presented in Section V as
well as describe and justify the fitting procedures used in our error
correction simulations. These fits enable us to reliably extrapolate to larger
code sizes than simulated, which is required for our analysis of resource
costs for large scale quantum computations (see Section VIII).

In addition, to presenting results for memory errors we also consider lattice
surgery errors. Lattice surgery is the primary technique we consider for
performing Clifford gates and magic state injection. It is a procedure for
measuring multi-qubit logical Pauli operators such as $X_{L}^{\otimes m}$ with
$m\geq 2$. It can be regarded as a code deformation where the $m$ logical
qubits are temporally merged into a code of $m-1$ logical qubits, and then
split into their constituent $m$ logical qubits. For the simple $m=2$ case, we
illustrate the space-time diagram for this process in Fig. 11 of the main
text. Here we present a more detailed in Fig. 52.

Figure 52: The three stages of lattice surgery corresponding to cross sections
(time slices) of the lattice surgery spacetime diagram in Fig. 11. Step 1
Prepare: data qubits between the surface code blocks are prepared in the
$|0\rangle$ state. Step 2 Merge: start measuring the $Z$ and $X$ stabilizers
indicated. The product of the $X$ stabilizers (highlighted with white
vertices) yields the outcome $X_{L1}X_{L2}$. However, a measurement error on a
white vertex will flip the outcome and so these stabilizer measurements must
be repeated $d_{m}$ times, with $d_{m}$ chosen sufficiently large to suppress
time-like errors to the desired probability. Step 3 Split: The qubits between
the initial surface code blocks are measured in the $Z$ basis. Note that it is
not possible to use $X$ basis measurements to disentangle as this would
measure $X_{L1}$ and $X_{L2}$. If the parity of the single-qubit $Z$
measurements is $``-1"$ then we must apply a Pauli correction $X_{L1}$ (or
equivalently $X_{L2}$) as a correction. Both the measurement of $X_{L1}X_{L2}$
and the estimated Pauli correction must be done fault-tolerantly after having
decoded the syndrome. In the case of $X_{L1}X_{L2}$, we choose 1 particular
timeslice $t_{p}$ and make an initial guess by multiplying all the white
vertices at time $t_{p}$. If the decoder assigns a measurement error to any
white vertex at time $t_{p}$, then we must account by flipping the
$X_{L1}X_{L2}$. If the accumulated physical $Z$ errors before time $t_{p}$
anticommute with $X_{L1}X_{L2}$ then we flip the outcome. For a similar
discussion of lattice surgery see Ref. [125]. Compared to Figs. 2 and 5 we use
a similar graphical representation but for simplicity omit the location of the
transmon, readout qubit and ATS.

An incredibly powerful and beautiful feature of lattice surgery is that
decoding via matching naturally extends over this 3D spacetime structure
without being interrupted by lattice surgery. However, some care is needed to
correctly account for boundaries and assess different failure modes. For a
planar surface code, it is well known that one must allow defects to match
with the appropriate boundaries in the space direction. When performing
lattice surgery, it is also important to match to appropriate boundaries in
the time directions.

To understand boundary effects, consider the more detailed explanation of
lattice surgery in Fig. 52. The procedure starts and ends with $Z$ basis state
preparations and measurements. A bit-flipped single qubit measurement or
preparation will yield a pair of $Z$ syndrome defects. That is, the initial
and final rounds of $Z$ stabilizer measurements are semi-ideal as they are
reconstructed from single-qubit information so that any defects occur in
pairs. These short $X$ strings are then easily matched. In contrast, an $X$
syndrome measurement error (at the start/end of lattice surgery) can lead to
an isolated defect and is potentially harmful as it flips the outcome of the
lattice surgery operation. However, for such an isolated defect near a time
boundary, the best explanation is clearly an isolated measurement error.
Therefore, we match these defects to red boundaries in the time direction.

As a warm-up to discussing the probability of time-like errors, we first recap
the error scaling properties of memory and logical $Z$-errors. Consider a
$d_{x}$ by $d_{z}$ surface code patch storing a logical qubit for $t$ surface
code cycles. We expect the total logical error probability to scale as
$(1-\exp(-\lambda t))/2$ for some constant rate $\lambda$, which for small
lambda is approximately $\sim\lambda t/2$. Furthermore, as $d_{x}$ increases
the number of paths across the code increases linearly, so we expect that
$\lambda\propto d_{x}$ and the total $Z$-logical probability to scale as

$P_{Z}=d_{x}tF(d_{z},p_{1},\ldots,p_{k}),$ (304)

for some function $F$ of $d_{z}$ and relevant hardware parameters
$(p_{1},\ldots,p_{k})$. Note that $d_{x}t$ corresponds to the area of the
vertical red boundaries in Fig. 11. For fixed parameters
$(d_{z},p_{1},\ldots,p_{m})$ the value of $F(d_{z},p_{1},\ldots,p_{m})$ can be
estimated by Monte Carlo simulation and evaluating $P_{Z}/(d_{x}t)$. For
simulation purposes, standard practice is to assume: at time zero, the system
is in a $``+1"$ eigenstates of all stabilizers; at time $t$, the round of
stabilizer measurements is ideal. This assumption introduces a finite size
effect error into $P_{Z}/(d_{x}t)$. This is suppressed by taking $t$ large,
and community folklore suggests that $t=\mathrm{max}[d_{z},d_{x}]$ will
suffice though one could push higher. The exact form of function $F$ can be
quite involved, though we know it will be exponentially suppressed by the
relevant distance $d_{z}$. Taking our sole experimental parameter to be
$\kappa_{1}/\kappa_{2}$ we find good fits of the form

$P_{Z}=d_{x}ta_{z}(b_{z}\kappa_{1}/\kappa_{2})^{c_{z}d_{z}},$ (305)

where $a_{z},b_{z},c_{z}$ are fitted parameters. For small $d_{x}$, there will
be a finite size effect so the scaling is not linear in $d_{x}$. However, we
can still use such a fit when $d_{x}$ is held constant provided we do not
attempt to extrapolate to larger $d_{x}$. Note that Eq. 305 is not necessarily
a leading order fit of the classical form $O(p^{(d-1)/2})$. Since the
probability of logical failures has a entropic/combinatorial component, it is
dominated by errors with a weight much larger than $(d-1)/2$. As such, we do
not attempt a leading order fit but rather it is appropriate to fit the
scaling exponent $c_{z}$.

We present the result of this fitting procedure in Fig. 53 and observe that it
works well over the interval $10^{-5}\leq\kappa_{1}/\kappa_{2}\leq 10^{-4}$.
At higher values of $\kappa_{1}/\kappa_{2}$, higher order contributions to a
logical $Z$ failure become important. Similarly, at lower values of
$\kappa_{1}/\kappa_{2}$, lower order contributions become important. Even if a
more sophisticated fitting function of $\kappa_{1}/\kappa_{2}$ is assumed, we
expect a finite range of applicability since there are other relevant
experimental parameters in the noise model.

Similar reasoning can be applied to timelike errors. The relevant boundary has
an area $\ell d_{x}$ where $\ell$ is the distance between the codeblocks. As
with $Z$-logical errors, the exponential decay of timelike errors follows from
a percolation theory analysis [127, 128] of a strings connecting the timelike
boundaries. As always in percolation problems, the probability of a
percolation event is exponentially suppressed in the distance between the
boundaries (whenever below some threshold). The relevant boundaries are
separated by a distance $d_{m}$, which we call the measurement distance, and
physically corresponds to the number of repeated rounds of syndrome
measurements during the merge step. Therefore, we fit to the ansatz

$P_{M}=\ell d_{x}a_{z}(b_{m}\kappa_{1}/\kappa_{2})^{c_{m}d_{m}},$ (306)

where $a_{m},b_{m},c_{m}$ are fitted parameters. To obtain an estimate of
$P_{M}$ we simulate the middle group of qubits in Fig. 52. We wish to isolate
the timelike errors and so freeze out $Z$-logical errors by assuming that the
left-most and right-most qubits are ideal and error-free. This is analogous to
the assumption of ideal measurements in a memory simulation. Furthermore,
since the $d_{z}$ distance is temporally extended during lattice surgery, such
errors will be rare in comparison. Again, this idealization introduces a
finite size effect that vanishes as $\ell$ grows relative to $d_{m}$.

Figure 53: Fitted results for simulation of $d_{x}=3$ surface code for logical
$Z$ and timelike errors. We fit according to the anstaz of Eqs. 305 and 306.
(Top) The logical $Z$ simulations for which we set $t=d_{z}$ and plot the
error probability divided by $t$. (Bottom) the probability of a timelike error
during lattice surgery for which we set $\ell=d_{m}-1$. All data points shown
are used in fitting. This is a truncated data set eliminating points above
$10^{-4}$ on the error rate axis and eliminating points outside the relevant
range of $\kappa_{1}/\kappa_{2}$.

We present the result of this fitting procedure for thin surface codes in Fig.
53 and observe that it works well over the interval
$10^{-5}\leq\kappa_{1}/\kappa_{2}\leq 5*10^{-4}$. We did not collect data for
$\kappa_{1}/\kappa_{2}\geq 5*10^{-4}$ as we had already identified that the
surface code overhead will be prohibitively large in this regime.

To the best of our knowledge, there have not been previous simulations that
investigate time-like errors in codes with boundaries and/or using circuit-
level noise. For instance, timelike errors were accounted for by Raussendorf
and Harrington [129] but using a toy, phenomenological noise model and
periodic boundary conditions in both space and time.

Widespread practice is to set $d_{m}=d_{x}=d_{z}$ but there is no a prior
reason to believe this is optimal. Indeed, just as physical bias in $X$ and
$Z$ noise leads to an asymmetry in our choice of $d_{x}$ and $d_{z}$, a
realistic noise model will influence the optimal choice of $d_{m}$. In later
calculations we find that $d_{m}=d_{z}-2$ is the most common optimal choice
for the main algorithm. Furthermore, in the design of magic state distillation
factories, the time-like errors are not critically important (see Table 14)
and so inside the factory $d_{m}$ can be set much smaller (by about a factor
1/2) than one would otherwise expect.

##  Appendix N Edge weights and decoding graphs for the repetition and surface
codes

In this section we provide the decoding graphs used to implement MWPM with the
repetition and surface codes considered in this paper. We also provide details
for computing the edge weights of all edges in a given graph.

###  N.1 Repetition code decoding graphs

(a)  (b)

Figure 54: (a) Circuit for measuring the stabilizers of the $d=5$ repetition
code. The dark rectangular boxes correspond to idling qubit locations. (b)
MWPM decoding graph for the $d=5$ repetition code where the syndrome
measurement is repeated five times.

The circuit for measuring the stabilizer of the $d=5$ repetition code is shown
in Fig. 54a and can straightforwardly be generalized to arbitrary code
distances. The corresponding graph for decoding the $d=5$ repetition code
using MWPM with five rounds of syndrome measurements is shown in Fig. 54b. The
purple vertices correspond to the measurement outcome of each ancilla qubit
(prepared in $|+\rangle$ and measured in the $X$-basis), and the horizontal
edges correspond to the physical data qubits. A vertex is highlighted if the
measurement outcome of the corresponding ancilla is non-trivial. We also add
space-like boundary vertices and edges (shown in red). For a given syndrome
measurement round (corresponding to a one-dimensional slice of the graph in
Fig. 54b), a boundary vertex is highlighted if on odd number of bulk vertices
in the corresponding one-dimensional slice are highlighted. To deal with
measurement errors, dashed grey vertical edges are added and connect vertices
of two one-dimensional graphs. Lastly, cross-diagonal edges (shown in green)
are added to deal with space-time correlated errors arising from failures at
CNOT gate locations (see below for explicit examples). More details for
implementing graph based decoding using MWPM can be found in Refs. [130, 131,
83].

We now describe how to compute the edge weights for each edge type of the
graph in Fig. 54b. For a given edge $e$, we must first compute the probability
$P_{e}$ of all failure events resulting in $e$ being highlighted. The weight
$w_{e}$ for the edge $e$ is then given by $w_{e}=-\log(P_{e})$ (see for
instance Refs. [130, 100, 132]). Such a prescription ensures that edges
arising from more likely failure events are chosen with higher probability
when finding the lowest weight path between two highlighted vertices. In what
follows, we will refer to $P_{e}$ as the edge weight probability for the edge
$e$.

The first and last data qubits in Fig. 54a have an extra idling location
compared to all other data qubits, and their edge weight probabilities are
labeled $P_{h_{1}}$ and $P_{h_{2}}$, whereas the other data qubits have edge
weight probabilities labeled by $P_{h}$. The dashed grey vertical edges
connecting are labeled $P_{v}$ and the green space-time correlated egdes are
labeled $P_{d}$. Next we define $P_{s}$ to be the probability of a $|+\rangle$
state preparation error, $P_{m}$ the probability of a measurement error,
$P_{i}$ the probability of an idle error and $P_{Z_{1}}$, $P_{Z_{2}}$ and
$P_{Z_{1}Z_{2}}$ the probability of a $Z\otimes I$, $I\otimes Z$ and $Z\otimes
Z$ CNOT failure (where the first qubit is the control qubit of the CNOT).

We now show how to compute $P_{h}$ and $P_{d}$ to leading order (the other
edge weight probabilities can be obtained using analogous methods). In the
case of a single failure, a bulk horizontal edge (say corresponding to an
error on the data qubit $q_{j}$) can be highlighted if either a $Z$ error
occurs at the idle location during the preparation of $|+\rangle$, a $Z\otimes
Z$ failure on the CNOT gate at the second time step with $q_{j}$ as a target
qubit, or an $I\otimes Z$ failure on the CNOT gate on the third time step
occuring in the previous syndrome measurement round. Hence we have

$\displaystyle
P^{(t_{1})}_{h}=P_{i}(1-P_{Z_{1}Z_{2}})+P_{Z_{1}Z_{2}}(1-P_{i}),$ (307)

and

$\displaystyle P^{(t>t_{1})}_{h}$
$\displaystyle=2P_{i}(1-P_{Z_{1}Z_{2}})(1-P_{Z_{2}})(1-P_{i})$
$\displaystyle+P_{Z_{1}Z_{2}}(1-P_{i})^{2}(1-P_{Z_{2}})$
$\displaystyle+P_{Z_{2}}(1-P_{i})^{2}(1-P_{Z_{1}Z_{2}}),$ (308)

where $P^{(t_{1})}_{h}=P_{h}$ in the first syndrome measurement round, and
$P^{(t>t_{1})}_{h}=P_{h}$ in all subsequent syndrome measurement rounds.

Now suppose a $Z\otimes Z$ error occurs on a CNOT gate in the third time step
of the syndrome measurement round $t$ resulting in a $Z$ data qubit error on
qubit $q_{j}$ while also flipping the measurement outcome of the ancilla
$a_{k}$. Note that if a $Z$ error had occurred on qubit $q_{j}$ prior to
applying the two CNOT gates, both ancillas $a_{k}$ and $a_{k+1}$ would be
measured non-trivially. Hence in round $t+1$ (and assuming no other failures),
the measurement outcome of $a_{k}$ will not change whereas the measurement
outcome of $a_{k+1}$ will change. To ensure that such an event is treated to
leading order, we introduce a green cross-diagonal edge as seen in Fig. 54b.
Also note that a $I\otimes Z$ error on a CNOT in the second time step also
results in such an edge. Hence we have that

$\displaystyle P_{d}=P_{Z_{1}Z_{2}}(1-P_{Z_{2}})+P_{Z_{2}}(1-P_{Z_{1}Z_{2}}).$
(309)

A similar analysis results in the following expressions for the remaining edge
weight probabilities

$\displaystyle P_{v}$
$\displaystyle=P_{m}(1-P_{s})(1-P_{Z_{1}})^{2}+P_{s}(1-P_{m})(1-P_{Z_{1}})^{2}$
$\displaystyle+2P_{Z_{1}}(1-P_{Z_{1}})(1-P_{s})(1-P_{m}),$ (310)
$\displaystyle P^{(t_{1})}_{h_{1}}=P^{(t_{1})}_{h},$ (311) $\displaystyle
P^{(t>t_{1})}_{h_{1}}$
$\displaystyle=3P_{i}(1-P_{i})^{2}(1-P_{Z_{1}Z_{2}})(1-P_{Z_{2}})$
$\displaystyle+P_{Z_{1}Z_{2}}(1-P_{i})^{3}(1-P_{Z_{2}})$
$\displaystyle+P_{Z_{2}}(1-P_{i})^{3}(1-P_{Z_{1}Z_{2}}),$ (312) $\displaystyle
P^{(t_{1})}_{h_{2}}=2P_{i}(1-P_{i})(1-P_{Z_{1}Z_{2}})+P_{Z_{1}Z_{2}}(1-P_{i})^{2},$
(313)

and

$\displaystyle P^{(t>t_{1})}_{h_{2}}=P^{(t>t_{1})}_{h_{1}}.$ (314)

###  N.2 Surface code decoding graphs

(a)  (b)  (c)

Figure 55: (a) Surface code lattice with $d_{x}=5$ and $d_{z}=7$. (b) Graph
used for decoding $X$ stabilizer measurement outcomes with both bulk and
boundary edge weight probability labels. (c) Graph used for decoding $Z$
stabilizer measurement outcomes with both bulk and boundary edge weight
probability labels.

The two-dimensional graphs for decoding the $X$ and $Z$ stabilizer measurement
outcomes of a $d_{x}=5$ and $d_{z}=7$ surface code, along with their
corresponding edge weight probability labels, are shown in Fig. 55. We will
show below the edges that need to be added when considering measurement errors
and space-time correlated errors arising from CNOT gate failures. However, we
first provide edge weight probabilities for the edges of the two-dimensional
graphs.

Let $G^{(2D)}_{(d_{x})}$ and $G^{(2D)}_{(d_{z})}$ be the two-dimensional
graphs corresponding to the $X$ and $Z$ stabilizer measurement outcomes. For
the graph $G^{(2D)}_{(d_{x})}$, we label the bulk edge weight probabilities by
$P^{(2D)}_{BLTRX}$ and $P^{(2D)}_{TLBRX}$. All other labels in Fig. 55b are
used for boundary edges. Similarly, for the graph $G^{(2D)}_{(d_{z})}$, we
label the bulk edge weight probabilities by $P^{(2D)}_{BLTRZ}$ and
$P^{(2D)}_{TLBRZ}$ with all other labels in Fig. 55c representing boundary
edge weight probabilities. In order to simplify the expressions for the edge
weight probabilities, we define the following function

$\displaystyle\Gamma(P_{1},P_{2},\cdots,P_{j};n_{1},n_{2},\cdots,n_{j})\equiv$
$\displaystyle\sum_{k=1}^{j}n_{k}P_{k}(1-P_{k})^{n_{k}-1}\prod_{i=1,i\neq
k}^{j}(1-P_{i})^{n_{i}}.$ (315)

In what follows, we define $P^{(P_{i}P_{j})}_{\text{CNOT}}$ to be the
probability that a CNOT gate failure results in a two-qubit Pauli error of the
form $P_{i}\otimes P_{j}$. We also define $P^{(P_{i})}_{\text{Id}}$ to be the
probability that a single-qubit idling location results in a $P_{i}$ Pauli
error on that qubit. To further simplify the edge weight probability
polynomials, we define the following probabilities:

$\displaystyle
P^{(1)}_{ZZCX}=P^{(ZZ)}_{\text{CNOT}}+P^{(ZY)}_{\text{CNOT}}+P^{(YZ)}_{\text{CNOT}}+P^{(YY)}_{\text{CNOT}},$
(316) $\displaystyle
P^{(1)}_{IZCX}=P^{(IZ)}_{\text{CNOT}}+P^{(XZ)}_{\text{CNOT}}+P^{(IY)}_{\text{CNOT}}+P^{(XY)}_{\text{CNOT}},$
(317) $\displaystyle P^{(1)}_{ZICX}$
$\displaystyle=P^{(ZI)}_{\text{CNOT}}+P^{(ZX)}_{\text{CNOT}}+P^{(ZY)}_{\text{CNOT}}+P^{(ZZ)}_{\text{CNOT}}+P^{(YI)}_{\text{CNOT}}$
$\displaystyle+P^{(YX)}_{\text{CNOT}}+P^{(YZ)}_{\text{CNOT}}+P^{(YY)}_{\text{CNOT}},$
(318) $\displaystyle P^{(2)}_{IZCX}$
$\displaystyle=P^{(IZ)}_{\text{CNOT}}+P^{(XZ)}_{\text{CNOT}}+P^{(IY)}_{\text{CNOT}}+P^{(XY)}_{\text{CNOT}}+P^{(ZI)}_{\text{CNOT}}$
$\displaystyle+P^{(ZX)}_{\text{CNOT}}+P^{(YI)}_{\text{CNOT}}+P^{(YX)}_{\text{CNOT}},$
(319) $\displaystyle P^{(3)}_{IZCX}$
$\displaystyle=P^{(IZ)}_{\text{CNOT}}+P^{(IY)}_{\text{CNOT}}+P^{(ZZ)}_{\text{CNOT}}+P^{(ZY)}_{\text{CNOT}}+P^{(XZ)}_{\text{CNOT}}$
$\displaystyle+P^{(XY)}_{\text{CNOT}}+P^{(YZ)}_{\text{CNOT}}+P^{(YY)}_{\text{CNOT}},$
(320)

and

$\displaystyle
P^{(2)}_{ZICX}=P^{(ZI)}_{\text{CNOT}}+P^{(YI)}_{\text{CNOT}}+P^{(ZX)}_{\text{CNOT}}+P^{(YX)}_{\text{CNOT}},$
(321)

$\displaystyle P^{(1)}_{d}=P^{(Z)}_{\text{Id}}+P^{(Y)}_{\text{Id}}.$ (322)

Using Eqs. 315, 316, 317, 318, 319, 320, 321 and 322 and the same methods as
in Section N.1, the leading order edge weight probabilities for the graph
$G^{(2D)}_{(d_{x})}$ are given by:

$\displaystyle
P^{(2D)}_{BLTRX}=\Gamma(P^{(1)}_{ZZCX},P^{(1)}_{IZCX},P^{(1)}_{d};1,1,1),$
(323) $\displaystyle P^{(2D)}_{TLBRX}=$
$\displaystyle\Gamma(P^{(1)}_{ZZCX},P^{(1)}_{IZCX},P^{(1)}_{ZICX},P^{(2)}_{IZCX},P^{(1)}_{d};2,2,1,1,1),$
(324) $\displaystyle
P_{C1X}=\Gamma(P^{(3)}_{IZCX},P^{(2)}_{ZICX},P^{(1)}_{d};1,1,1),$ (325)
$\displaystyle P_{TB2X}=P^{(2D)}_{BLTRX},$ (326) $\displaystyle
P_{TB1X}=\Gamma(P^{(1)}_{ZZCX},P^{(1)}_{IZCX},P^{(2)}_{IZCX},P^{(1)}_{d};2,1,1,1),$
(327) $\displaystyle
P_{C2X}=\Gamma(P^{(3)}_{IZCX},P^{(2)}_{IZCX},P^{(1)}_{ZZCX},P^{(1)}_{d};1,1,1,1),$
(328) $\displaystyle
P_{MRX1}=\Gamma(P^{(3)}_{IZCX},P^{(2)}_{ZICX},P^{(1)}_{d};1,2,1),$ (329)
$\displaystyle P_{MRX2}=$
$\displaystyle\Gamma(P^{(3)}_{IZCX},P^{(1)}_{IZCX},P^{(1)}_{ZICX},P^{(2)}_{IZCX},P^{(1)}_{d};1,2,1,1,1),$
(330) $\displaystyle
P_{C3X}=\Gamma(P^{(3)}_{IZCX},P^{(1)}_{IZCX},P^{(1)}_{ZICX},P^{(1)}_{d};1,1,1,1),$
(331) $\displaystyle P_{BB2X}=P^{(2D)}_{BLTRX},$ (332) $\displaystyle
P_{BB1X}=\Gamma(P^{(1)}_{IZCX},P^{(1)}_{ZZCX},P^{(1)}_{ZICX},P^{(1)}_{d};2,1,1,1),$
(333) $\displaystyle P_{C4X}=P_{C3X},$ (334) $\displaystyle P_{MLX1}=$
$\displaystyle\Gamma(P^{(3)}_{IZCX},P^{(1)}_{IZCX},P^{(2)}_{ZICX},P^{(1)}_{ZICX},P^{(1)}_{d};1,1,1,1,1),$
(335)

and

$\displaystyle P_{MLX2}=P_{MLX1}.$ (336)

For the graph $G^{(2D)}_{(d_{z})}$, we first define the following
probabilities:

$\displaystyle
P^{(1)}_{XXCX}=P^{(XX)}_{\text{CNOT}}+P^{(XY)}_{\text{CNOT}}+P^{(YX)}_{\text{CNOT}}+P^{(YY)}_{\text{CNOT}},$
(337) $\displaystyle
P^{(1)}_{XICX}=P^{(XI)}_{\text{CNOT}}+P^{(YI)}_{\text{CNOT}}+P^{(XZ)}_{\text{CNOT}}+P^{(YZ)}_{\text{CNOT}},$
(338) $\displaystyle
P^{(1)}_{IXCX}=P^{(IX)}_{\text{CNOT}}+P^{(ZX)}_{\text{CNOT}}+P^{(IY)}_{\text{CNOT}}+P^{(ZY)}_{\text{CNOT}},$
(339) $\displaystyle
P^{(2)}_{IXCX}=P^{(IX)}_{\text{CNOT}}+P^{(IY)}_{\text{CNOT}}+P^{(ZX)}_{\text{CNOT}}+P^{(ZY)}_{\text{CNOT}}+P^{(XX)}_{\text{CNOT}}$
$\displaystyle+P^{(XY)}_{\text{CNOT}}+P^{(YX)}_{\text{CNOT}}+P^{(YY)}_{\text{CNOT}},$
(340) $\displaystyle
P^{(3)}_{IXCX}=P^{(IX)}_{\text{CNOT}}+P^{(IY)}_{\text{CNOT}}+P^{(ZX)}_{\text{CNOT}}+P^{(ZY)}_{\text{CNOT}}+P^{(XI)}_{\text{CNOT}}$
$\displaystyle+P^{(XZ)}_{\text{CNOT}}+P^{(YI)}_{\text{CNOT}}+P^{(YZ)}_{\text{CNOT}},$
(341) $\displaystyle
P^{(2)}_{XICX}=P^{(XI)}_{\text{CNOT}}+P^{(YI)}_{\text{CNOT}}+P^{(XX)}_{\text{CNOT}}+P^{(YX)}_{\text{CNOT}}+P^{(XZ)}_{\text{CNOT}}$
$\displaystyle+P^{(YZ)}_{\text{CNOT}}+P^{(XY)}_{\text{CNOT}}+P^{(YY)}_{\text{CNOT}},$
(342)

and

$\displaystyle P^{(2)}_{d}=P^{(X)}_{\text{Id}}+P^{(Y)}_{\text{Id}}.$ (343)

Using Eqs. 337, 338, 339, 340, 341, 342 and 343, the leading order edge weight
probabilities for the graph $G^{(2D)}_{(d_{z})}$ are given by:

$\displaystyle
P^{(2D)}_{BLTRZ}=\Gamma(P^{(1)}_{XXCX},P^{(1)}_{XICX},P^{(2)}_{d};1,1,1),$
(344) $\displaystyle P^{(2D)}_{TLBRZ}=$
$\displaystyle\Gamma(P^{(1)}_{XXCX},P^{(1)}_{XICX},P^{(2)}_{IXCX},P^{(3)}_{IXCX},P^{(2)}_{d};2,2,1,1,1),$
(345) $\displaystyle
P_{C1Z}=\Gamma(P^{(2)}_{XICX},P^{(1)}_{IXCX},P^{(2)}_{d};1,1,1),$ (346)
$\displaystyle P_{TB1Z}=$
$\displaystyle\Gamma(P^{(2)}_{XICX},P^{(1)}_{XICX},P^{(2)}_{IXCX},P^{(3)}_{IXCX},P^{(1)}_{XXCX},P^{(2)}_{d}$
$\displaystyle;1,1,1,1,1,1),$ (347) $\displaystyle
P_{TB2Z}=\Gamma(P^{(2)}_{XICX},P^{(1)}_{IXCX},P^{(2)}_{d};1,2,1),$ (348)
$\displaystyle P_{C2Z}=P_{C1Z},$ (349) $\displaystyle
P_{MRZ1}=\Gamma(P^{(1)}_{XXCX},P^{(1)}_{XICX},P^{(3)}_{IXCX},P^{(2)}_{d};1,2,1,1),$
(350) $\displaystyle P_{MRZ2}=P^{(2D)}_{BLTRZ},$ (351) $\displaystyle
P_{C3Z}=P_{C1Z},$ (352) $\displaystyle P_{BB1Z}=$
$\displaystyle\Gamma(P^{(2)}_{XICX},P^{(1)}_{IXCX},P^{(2)}_{IXCX},P^{(1)}_{XICX},P^{(2)}_{d};1,1,1,1,1),$
(353) $\displaystyle P_{BB2Z}=P_{BB1Z},$ (354) $\displaystyle
P_{C4Z}=\Gamma(P^{(2)}_{XICX},P^{(1)}_{XICX},P^{(2)}_{IXCX},P^{(2)}_{d};1,1,1,1),$
(355) $\displaystyle
P_{MLZ1}=\Gamma(P^{(1)}_{XXCX},P^{(1)}_{XICX},P^{(2)}_{IXCX},P^{(2)}_{d};1,2,1,1),$
(356)

and

$\displaystyle P_{MLZ2}=P^{(2D)}_{BLTRZ}.$ (357)

(a)  (b)

Figure 56: (a) Graph used for decoding $X$-type stabilizer measurements which
include vertical edges (dashed gray edges) for dealing with measurement errors
and space-time correlated edges for correcting errors arising from CNOT gate
failures causing two different syndrome measurement outcomes in consecutive
rounds. (b) Same as in (a) but for $Z$-type stabilizer measurements.

We now consider the three-dimensional version of the graphs in Figs. 55b and
55c (which we label $G^{(3D)}_{(d_{x})}$ and $G^{(3D)}_{(d_{z})}$) to deal
with measurement errors in addition to space-time correlated errors arising
from CNOT gate failures. As an example, consider an $I\otimes Z$ error arising
from a CNOT gate failure in the second time-step of an $X$-type (red)
plaquette during the $k$’th syndrome measurement round. Such a failure adds a
$Z$ data-qubit error which propagates through the CNOT in the fifth time-step
of the top right red $X$-type plaquette. Let $v_{j}$ and $v_{k}$ be the
vertices corresponding to the measurement outcomes of the two ancilla qubits
which would detect the $Z$ error. Assuming there were no other failures, only
one of the two vertices (say $v_{j}$) changes from rounds $k-1$ to round $k$.
In the next syndrome measurement round, both $X$-type plaquettes will detect
the $Z$ data qubit error the ancilla qubits in both $X$-type plaquettes will
be highlighted. Hence only the vertex $v_{k}$ changes from round $k$ to $k+1$.
In order to ensure that the highlighted ancillas arising from failures as in
the example considered here can be reached by a single edge when implementing
MWPM, the dark green edges in the graph of Fig. 56a (labeled
$P^{(\text{bulk})}_{d_{1},X}$) are added to the graph of Fig. 55b. The other
types of space-time correlated edges are distinguished by their color and
associated label (all edges of the same color have identical edge-weight
probabilities). Similarly, we add the dashed grey vertical edges in Figs. 56a
and 56b connecting identical vertices from two consecutive syndrome
measurement rounds to deal with measurement errors. The edge weight
probabilities of such edges are labeled $P^{X}_{V}$ and $P^{Z}_{V}$. Note that
there are also solid dark vertical edges at some of the boundaries of the
graphs where weight-two $X$-type and $Z$-type stabilizers occur in Fig. 55a.
These vertical edges have different edges weight probabilities which are
labeled $P^{X,\text{bound}}_{V}$ and $P^{Z,\text{bound}}_{V}$.

In order to avoid making the visualization of the three-dimensional graphs too
cumbersome, in Figs. 56a and 56b we only included vertices corresponding to
the first two syndrome measurement rounds. Further, the two-dimensional edges
from the second round were omitted in order to maintain focus on the vertical
and space-time correlated edges connecting vertices from two consecutive
syndrome measurement rounds.

Let

$\displaystyle
P_{VCX}=P^{(ZI)}_{\text{CNOT}}+P^{(ZX)}_{\text{CNOT}}+P^{(YI)}_{\text{CNOT}}+P^{(YX)}_{\text{CNOT}},$
(358)

and

$\displaystyle
P_{VCZ}=P^{(IX)}_{\text{CNOT}}+P^{(IY)}_{\text{CNOT}}+P^{(ZX)}_{\text{CNOT}}+P^{(ZY)}_{\text{CNOT}}.$
(359)

Further, let $P_{s}$ be the probability of preparing $|-\rangle$ instead of
$|+\rangle$ and $P_{m}$ be the probability that a $X$-basis measurement
outcome is flipped. The edge weight probabilities corresponding to the dashed
grey edges in Figs. 56a and 56b (i.e. the vertical edges of
$G^{(3D)}_{(d_{x})}$ and $G^{(3D)}_{(d_{z})}$) are given by

$\displaystyle P^{X}_{V}=\Gamma(P_{VCX},P_{s},P_{m};4,1,1),$ (360)
$\displaystyle P^{X,\text{bound}}_{V}=\Gamma(P_{VCX},P_{s},P_{m};2,1,1),$
(361) $\displaystyle P^{Z}_{V}=\Gamma(P_{VCZ},P_{s},P_{m};4,1,1),$ (362)

and

$\displaystyle P^{Z,\text{bound}}_{V}=\Gamma(P_{VCZ},P_{s},P_{m};2,1,1).$
(363)

(a)  (b)

Figure 57: (a) Fictitious identity gates illustrating the possible correlated
errors arising before the $X$-basis measurement of the $X$-type ancilla
qubits. Grey squares correspond to the first qubit, blue triangles to the
second qubit and green circles to the third qubit. (b) $X$-type decoding graph
with added edges to correct correlated errors. The edge-weight probabilities
of the orange cross-edges are labeled $P_{\text{cross}}$. We also add red
edges with edge-weight probabilities labelled $P_{d,\text{corr}}$ at the
bottom row of the graph.

Next we consider the edge-weight probabilities for the space-time correlated
edges of $G^{(3D)}_{(d_{x})}$. The dark green edges labeled by
$P^{(\text{bulk})}_{d_{1},X}$ have different values at the boundaries (dashed
dark green edges in the first and last column of Fig. 56a) and are labeled by
$P^{(\text{bound})}_{d_{1},X}$. We have that

$\displaystyle
P^{(\text{bulk})}_{d_{1},X}=\Gamma(P^{(1)}_{IZCX},P^{(1)}_{ZZCX},P^{(2)}_{ZICX};1,1,2),$
(364)

and

$\displaystyle
P^{(\text{bound})}_{d_{1},X}=\Gamma(P^{(1)}_{IZCX},P^{(1)}_{ZZCX},P^{(2)}_{ZICX};1,1,1).$
(365)

The edge weight probability $P_{d_{2},X}$ (represented by the light green
edges in Fig. 56a) is given by

$\displaystyle P_{d_{2},X}=\Gamma(P^{(1)}_{IZCX},P^{(1)}_{ZZCX};1,1).$ (366)

Lastly, the edge weight probability $P_{d_{3},X}$ (represented by the yellow
edges in Fig. 56a) is given by

$\displaystyle P_{d_{3},X}=P_{d_{2},X}.$ (367)

Similarly, for the graph $G^{(3D)}_{(d_{z})}$, the edge weight probability
$P_{d_{1},Z}$ (represented by the light green edges) is given by

$\displaystyle P_{d_{1},Z}=\Gamma(P^{(1)}_{XICX},P^{(1)}_{XXCX};1,1).$ (368)

The bulk and boundary edge weight probabilities $P^{(\text{bulk})}_{d_{2},Z}$
(dark green edges) and $P^{(\text{bound})}_{d_{2},Z}$ (dashed dark green
edges) are given by

$\displaystyle
P^{(\text{bulk})}_{d_{2},Z}=\Gamma(P^{(1)}_{XICX},P^{(1)}_{XXCX},P^{(1)}_{IXCX};1,1,2),$
(369)

and

$\displaystyle
P^{(\text{bound})}_{d_{2},Z}=\Gamma(P^{(1)}_{XICX},P^{(1)}_{XXCX},P^{(1)}_{IXCX};1,1,1).$
(370)

Lastly, the edge weight probability $P_{d_{3},Z}$ (represented by the yellow
edges) is given by

$\displaystyle P_{d_{3},Z}=P_{d_{1},Z}.$ (371)

###  N.3 Adding edges for dealing with correlated errors

In this section we provide a modified version of the graph
$G^{(3D)}_{(d_{x})}$ (described in Section N.2) which includes extra edges to
deal with two-qubit and three-qubit correlated errors arising from the micro
oscillations described in Section B.5.

For the purposes of the edge weight analysis, in Fig. 57a, we illustrate
fictitious two-qubit and three-qubit gates which act as the identity and which
are applied immediately prior to the $X$-basis measurements of the red
plaquettes. The two-qubit correlated errors can be viewed as an $Z\otimes
I\otimes Z$-type error at a three-qubit gate location, where the $Z$ errors
act on the qubits adjacent to the grey squares and green circles of such
gates. Such errors occur with probability $P_{\text{cd}}$. Similarly, the
three-qubit correlated errors can be viewed as an $Z\otimes Z\otimes Z$-type
error at a three-qubit gate location. Such errors occur with probability
$P_{\text{ct}}$. Additionally, there can be correlated errors occurring
between the ancilla and data qubits at the top and bottom boundaries of the
lattice in Fig. 57a. Hence, we add fictitious two-qubit gate locations at such
boundaries as shown in the figure.

In order to incorporate the different types of correlated errors mentioned
above into our MWPM decoding protocol, extra edges are added to the graph
$G^{(3D)}_{(d_{x})}$ as shown in Fig. 57b. The first type of extra edges are
two-dimensional cross edges shown in orange that deal with two and three-qubit
correlated errors arising at the three-qubit fictitious gate locations of Fig.
57a. The edge-weight probabilities of such edges are labeled
$P^{(\text{bulk})}_{\text{cross}}$. Due to boundary effects, we also add
dashed orange edges with edge-weight probabilities labeled
$P^{(\text{bound})}_{\text{cross}}$. Additionally, extra space-time correlated
edges (shown in red) are added at the bottom row of the graph in Fig. 57b with
edge weight probabilities labeled by $P_{d,\text{corr}}$. Note that the two-
qubit correlated errors arising at the top boundary of Fig. 57a result in
space-time correlated edges which are already included in
$G^{(3D)}_{(d_{x})}$.

Figure 58: Example of a decoding graph for correcting timelike errors using a
$d=5$ repetition code with $d_{m}=4$. The top and bottom boundary edges (with
zero weight) and vertices are shown in blue and are connected by a blue edge
with zero weight. As explained in Appendix M, we have removed the left and
two-dimensional black edges (which correspond to the left and rightmost
qubits) to isolate timelike errors.

(a)  (b)

Figure 59: (a) Implementation of the timelike decoding protocol in the
presence of a single measurement error when measuring the stabilizer
$X_{2}X_{3}$ during the first round. The minimum weight path matches to the
bottom boundary going through the vertex $v^{(1)}_{2}$ whose outcome is
correctly flipped (illustrated by the yellow star). (b) Same as in (a) but
with an additional measurement error occurring in the second round when
measuring $X_{2}X_{3}$. In this case, the minimum weight path matches to the
top boundary and fails to flip the measurement outcome of $v^{(1)}_{2}$ (which
is incorrect given the measurement error in the first round) resulting in a
logical failure.

In addition to the extra edges added to $G^{(3D)}_{(d_{x})}$, the edge-weight
probabilities of a subset of the edges already included in
$G^{(3D)}_{(d_{x})}$ need to be renormalized. The edge-weight probabilities of
the added edges in addition to the renormalized edges are given by:

$\displaystyle
P^{(\text{bulk})}_{\text{cross}}=\Gamma(P_{\text{ct}}P_{\text{cd}};2,2),$
(372) $\displaystyle
P^{(\text{bound})}_{\text{cross}}=\Gamma(P_{\text{ct}}P_{\text{cd}};1,1),$
(373) $\displaystyle P_{d,\text{corr}}=P_{\text{ct}},$ (374) $\displaystyle
P_{TB2X}=\Gamma(P^{(1)}_{ZZCX},P^{(1)}_{IZCX},P^{(1)}_{d},P_{\text{cd}};1,1,1,1),$
(375) $\displaystyle
P_{TB1X}=\Gamma(P^{(1)}_{ZZCX},P^{(1)}_{IZCX},P^{(2)}_{IZCX},P^{(1)}_{d},P_{\text{cd}};2,1,1,1,1),$
(376) $\displaystyle P_{C2X}=$
$\displaystyle\Gamma(P^{(3)}_{IZCX},P^{(2)}_{IZCX},P^{(1)}_{ZZCX},P^{(1)}_{d},P_{\text{cd}},P_{\text{ct}};1,1,1,1,1,1),$
(377) $\displaystyle
P_{BB2X}=\Gamma(P^{(1)}_{ZZCX},P^{(1)}_{IZCX},P^{(1)}_{d},P_{\text{cd}};1,1,1,1),$
(378) $\displaystyle
P_{BB1X}=\Gamma(P^{(1)}_{IZCX},P^{(1)}_{ZZCX},P^{(1)}_{ZICX},P^{(1)}_{d},P_{\text{cd}};2,1,1,1,1),$
(379) $\displaystyle P_{C4X}=$
$\displaystyle\Gamma(P^{(3)}_{IZCX},P^{(1)}_{IZCX},P^{(1)}_{ZICX},P^{(1)}_{d},P_{\text{cd}},P_{\text{ct}};1,1,1,1,1,1),$
(380) $\displaystyle
P^{X}_{V}=\Gamma(P_{VCX},P_{s},P_{m},P_{\text{ct}};4,1,1,1),$ (381)

For the space-time correlated edges, at the top row of the graph in Fig. 57b,
we have

$\displaystyle
P^{(\text{bound,top})}_{d_{1},X}=\Gamma(P^{(1)}_{IZCX},P^{(1)}_{ZZCX},P^{(2)}_{ZICX},P_{\text{ct}};1,1,1,1),$
(382)

whereas at the bottom boundary $P^{(\text{bound,bottom})}_{d_{1},X}$ is given
by Eq. 365. Similarly, at the top row of Fig. 57b, we have

$\displaystyle
P^{(\text{top})}_{d_{2},X}=\Gamma(P^{(1)}_{IZCX},P^{(1)}_{ZZCX},P_{\text{ct}};1,1,1),$
(383)

whereas anywhere else in the graph $P_{d_{2},X}$ is given by Eq. 366.

###  N.4 Decoding time-like errors

In this section, we show how the decoding graphs in addition to the MWPM
decoding protocols need to be modified for correcting timelike errors
discussed in Appendix M. Since visualizing three-dimensional graphs can be
challenging, we focus on correcting timelike errors in the context of the
repetition code, even though timelike errors occur in surface code patches
when implementing our lattice surgery schemes. However the main techniques
discussed in the context of the repetition can straightforwardly be applied to
the rotated surface code.

An example of a decoding graph for timelike errors occurring in a $d=5$
repetition code with $d_{m}=4$ is given in Fig. 58. Note that unlike Fig. 54b,
the boundary edges and vertices (shown in blue) are at the top and bottom of
the graph since we follow the matching protocol explained in Fig. 11. In
particular, we are considering a setting analogous to Fig. 52, where data
qubits between two repetition code patches are initially prepared in the
$|0\rangle$ state, and the product of the $X$-stabilizers yields the outcome
$X_{L1}X_{L2}$. Note that although the measurement of each $X$-stabilizer in
the first round is random, the parity of the product of all measurement
outcomes gives the outcome of $X_{L1}X_{L2}$. Due to possibility of
measurement errors, measurements of the $X$-stabilizers are repeated $d_{m}$
times. MWPM is then performed over the entire syndrome history in order to
determine if measurement errors occurred during the measurement of
$X$-stabilizers in the first round. We thus summarize the decoding protocol as
follows:

  1. 1.

Repeat the measurement of all $X$-stabilizers $d_{m}$ times.

  2. 2.

Implement MWPM using a timelike decoding graph (such as the one in Fig. 58).
If there is an odd number of highlighted vertices (purple vertices in Fig.
58), highlight a boundary vertex (the particular choice is irrelevant).

  3. 3.

Let $v^{(1)}_{j}$ correspond to the $j$’th $X$-stabilizer measurement outcome
in the first round (represent by the $j$’th purple vertex, starting from the
left, in the first layer of the graph in Fig. 58). If there are highlighted
timelike edges (i.e. vertical edges) incident to $v^{(1)}_{j}$, flip the
measurement outcome of $v^{(1)}_{j}$.

  4. 4.

Let $\tilde{v}^{(1)}_{j}$ correspond to the values of $v^{(1)}_{j}$ after
implementing MWPM and performing the appropriate measurement flips described
in the previous step. The outcome $v_{f}$ of $X_{L1}X_{L2}$ is given by
$v_{f}=\sum^{d-1}_{j=1}\tilde{v}^{(1)}_{j}\text{mod}(2)$.

In Fig. 59, we provide two examples for the implementation of the timelike
decoder. In Fig. 59a, we consider the case where a single measurement error
occurs in the first round when measuring the stabilizer $X_{2}X_{3}$. Since
the syndrome changes between the first and second round, the second vertex
(starting from the left to right) of the second two-dimensional layer is
highlighted. A boundary vertex is also highlighted to ensure the total number
of highlighted vertices is even. After implementing MWPM, the minimum weight
path connecting the two highlighted vertices correctly passes through
$v^{(1)}_{2}$ in the timelike direction. The decoder then flips the
measurement outcome of $X_{2}X_{3}$ in the first round resulting in the
correct parity for the he outcome of $X_{L1}X_{L2}$. In Fig. 59b, we consider
a similar setting but with two consecutive measurement errors of the
stabilizer $X_{2}X_{3}$ occurring in the first and second round. In this case,
the syndrome only changes between the second and third round resulting in the
red highlighted vertex shown in Fig. 59b. After implementing MWPM, the minimum
weight path connect to the top boundary and so the measurement outcome of
$X_{2}X_{3}$ in the first round is incorrectly left unchanged resulting in a
logical failure.

We conclude this section with an important remark. Suppose a measurement error
occurs in the first round when measuring the $X$-type stabilizer $S^{(x)}_{j}$
of a given code. In order to prevent highlighted timelike edges from being
incident to the vertex $v^{(1)}_{j}$, one requires additional measurement
errors such that minimum weight paths are matched in the top timelike portion
of the decoding graph (as in Fig. 59b). By increasing $d_{m}$ to $d_{m}+2$,
one requires an additional measurement error to guarantee that the minimum
weight path is not incident to $v^{(1)}_{j}$, thus explaining the scaling in
Eq. 306.

##  Appendix O Toffoli state distillation (TDTOF)

###  O.1 Prior state of the art

Here we give a high-level comparison of how our TDTOF protocol compares to the
prior art in terms of magic state conversion rates.

Early protocols for fault-tolerant quantum computation focused on TOF state
preparation in concatenated codes [133] or they protected against 1 type of
error [134]. However, none of these protocols are suitable for protecting
against generic noise in topological (e.g. repetition or surface) codes.

A more modern approach to magic state distillation uses a supply of low
fidelity $T$ magic states. There are many protocols for distillation of noisy
$T$-states to purer $T$-states [40, 37, 135, 38, 136, 137]. One can also use
$T$-states as input to protocols that output other types of magic states,
including TOF states [138, 139, 140, 141, 142, 143]. For instance, there were
parallel discoveries of protocols [138, 139] that distill 1 TOF state from 8
noisy $T$ states, which we will write as $8T\rightarrow 1\text{TOF}$. This was
later generalized using synthillation [141, 142] to a family of protocols
$(6k+2)T\rightarrow k\text{TOF}$ for any integer $k$. However, in some
settings, the supply of noisy TOF states can be prepared with better
fidelities than the noisy $T$ states. For instance, in this paper we have
shown that in system with highly-biased noise we can use a repetition encoding
and the BUTOF protocol to realize TOF state at better fidelities than physical
TOF gates, with only a mild additional resource cost.

It has been previously noted [144] that triorthogonal codes enabling
$(6k+8)T\rightarrow(2k)T$ state distillation can also be lifted to perform
$(6k+8)\text{TOF}\rightarrow(2k)\text{TOF}$. The conversion rate of these
protocols is $2k/(6k+8)$, which is poor when $k$ is small (starting at $1/7$
for $k=1$) but improving when $k$ is larger (approaching $1/3$ for
$k\rightarrow\infty$). However, the ratio of inputs to outputs is not the sole
metric of importance; also crucial is the space-time complexity of the
Clifford circuit implementing the distillation protocol. Previous analysis has
found that the space-time complexity of Clifford distillation circuits tends
to be more favorable for simpler protocols using smaller block sizes [136, 79,
91] and that this effect can outweigh the improvement of conversion rate in
the asymptotic regime. In other words, the desiderata for distillation
protocols converting $n\rightarrow k$ magic states, are that: the protocol has
a good rate, so $k/n$ is large; the protocol is compact so $n$ is as small as
possible. These desiderata are in tension since rates tends to improve
asymptotically as block size $n$ is increased. A protocol satisfying these
desiderata, will likely have a small space-time footprint when compiled down
to physical qubits and gates. In this work, we present a
$8\text{TOF}\rightarrow 2\text{TOF}$ protocol that protects against any single
location fault (of $X$, $Y$ or $Z$ type), so it has a relatively high
conversion rate of $1/4$ without needing to scale to large blocks. In
contrast, to achieve the same conversion rate using the ideas of Ref. [144]
would require a much larger $32\text{TOF}\rightarrow 8\text{TOF}$ protocol.

###  O.2 Transversality proofs

Here we prove that the trio of $[[8,2,2]]$ codes introduced in Section VII
have the required CCZ tranversality properties. Recall that CCZ is a 3-qubit
gate that adds a “$-1$” phase to the state $|111\rangle$ and “$+1$” to all
other computational basis states. The corresponding magic state
$|\text{CCZ}\rangle$ differs from $|\text{TOF}\rangle$ by a single Hadamard
gate. For reasons of mathematical elegance, it is simpler to work mostly in
terms of $|\text{CCZ}\rangle$ state distillation, but our final description of
the distillation protocol will be presented in terms of $|\text{TOF}\rangle$
states.

For our trio of codes, each block encodes $k=2$ logical qubits into $n=8$
physical qubits and can detect an error on any single qubit. In the main text,
we define the code using Fig. 16 and for completeness we give an algebraic
definition here. All three code blocks one $X$-stabilizer $X^{\otimes 8}$ but
different logical $X$ operators

$\displaystyle X_{L1A}$ $\displaystyle=(X\otimes X\otimes
X\otimes{\mathchoice{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}\otimes{\mathchoice{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}\otimes X\otimes{\mathchoice{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}\otimes{\mathchoice{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm 1\mskip-5.0mul}})_{A},$ (384)
$\displaystyle X_{L2A}$ $\displaystyle=(X\otimes X\otimes{\mathchoice{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}\otimes{\mathchoice{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}\otimes X\otimes
X\otimes{\mathchoice{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}\otimes{\mathchoice{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm
1\mskip-5.0mul}})_{A},$ $\displaystyle X_{L1B}$ $\displaystyle=(X\otimes
X\otimes{\mathchoice{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}\otimes{\mathchoice{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}\otimes X\otimes X\otimes{\mathchoice{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}\otimes{\mathchoice{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm 1\mskip-5.0mul}})_{B},$ $\displaystyle
X_{L2B}$ $\displaystyle=(X\otimes{\mathchoice{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}\otimes
X\otimes{\mathchoice{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}\otimes X\otimes{\mathchoice{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}\otimes X\otimes{\mathchoice{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm 1\mskip-5.0mul}})_{B},$ $\displaystyle
X_{L1C}$ $\displaystyle=(X\otimes{\mathchoice{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}\otimes
X\otimes{\mathchoice{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}\otimes X\otimes{\mathchoice{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}\otimes X\otimes{\mathchoice{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm 1\mskip-5.0mul}})_{C},$ $\displaystyle
X_{L2C}$ $\displaystyle=(X\otimes X\otimes X\otimes{\mathchoice{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}\otimes{\mathchoice{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}\otimes
X\otimes{\mathchoice{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}\otimes{\mathchoice{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm
1\mskip-5.0mul}})_{C},$

where the index $\\{A,B,C\\}$ labels the three different codeblocks and the
numerical index labels the 2 logical qubits in this code block. We write
$(\ldots)_{D=A,B,C}$ to emphasize that the operator acts non-trivially on
codeblock $D$ and trivially on other codeblocks. While the code blocks share
the same $X$-stabilizer, they will have different $Z$-stabilizers as a
consequence of having different logical $X$ operators.

We say a set of $[[n,k,d]]$ codes is CCZ transversal whenever
$\text{CCZ}^{\otimes n}$ performs a logical $\text{CCZ}^{\otimes k}$ gate.
Note that if we take three copies of a CSS code that has a transversal $T$
gate (so that $T^{\otimes n}=T_{L}$ or similar, then it must also be CCZ
transversal). This is simply because CSS codes have transversal CNOT gates and
we can synthesize CCZ gates from CNOT and $T$ gates. Essentially, this is the
observation exploited to construct $(6k+8)\text{TOF}\rightarrow(2k)\text{TOF}$
protocols [144]. However, it is possible for a trio of codes to be CCZ
transversal, but not be $T$ transversal. To the best of our knowledge this was
first shown for the 3D surface codes by showing an equivalence (via unfolding)
to 3D colour codes [145]. Later, Vasmer and Brown gave a more direct proof
that the 3D surface codes are CCZ transversal [146].

Here, we use similar proof techniques to Vasmer and Brown [146], though
generalized (to $k>1$) and with a new code construction that code not appear
to be a surface code. From just the $X$-stabilizer and logical operator
information, we will prove that our trio of codes are CCZ transversal. The key
element of the proof is a lemma relating tranversality to the support of
logical $X$ operators and $X$ stabilizers. The lemma requires that for
$j=1,2$, the operators $X_{LjA}$, $X_{LjB}$ and $X_{LjC}$ share support on an
odd number of qubit indices. Furthermore, we need that for any other choice of
three $X$ operators (either logical $X$ or $X$ stabilizer) with one selected
from each code block, they must share support on an even number of qubit
indices. It is easy to verify the operators provided above have this property.

We define codes here using slightly different notation from the main text.
Given an $n$-qubit bit string $\mathbf{s}=(s_{1},s_{2},\ldots,s_{n})$, we use
$X[\mathbf{s}]:=\otimes_{j}X^{s_{j}}$. For example, if

$\mathbf{u}=(1,1,1,0,0,1,0,0),$ (385)

then

$X[\mathbf{u}]=X\otimes X\otimes X\otimes{\mathchoice{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm
1\mskip-5.0mul}}\otimes{\mathchoice{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}\otimes
X\otimes{\mathchoice{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm
1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}\otimes{\mathchoice{\rm
1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm 1\mskip-5.0mul}}.$
(386)

With this notation we can define an $[[n,k,d]]$ CSS code using a binary
$G$-matrix representation as follows.

Let $G$ be a binary matrix that is row-wise linearly independent and
partitioned as follows

$G=\left(\begin{array}[]{c}G_{1}\\\ \hline\cr G_{0}\end{array}\right),$ (387)

where $G$ has $n$ columns and $G_{1}$ has $k$ rows. Letting $m$ denote the
number of rows in $G_{0}$, then for a non-trivial ($d\geq 2$) code we know
$m\geq 1$. Here, we review the relevant facts for $G$-matrices, but for
additional details and proofs refer the reader to Refs. [143, 142, 141, 90].
This allows us to define a CSS code with all-zero logical state

$|(0,\ldots,0)\rangle_{L}=2^{-m/2}\sum_{\mathbf{u}\in\mathbb{F}_{2}^{m}}|\mathbf{u}G_{0}\rangle.$
(388)

Note we use bold-font for row vectors. The notation $\mathbf{u}G_{0}$
represents left multiplication of matrix $G_{0}$ by the row vector
$\mathbf{u}$, performed modulo 2, which will produce a length $n$ row-vector
describing a physical, computational basis state. The set of all
$\mathbf{u}G_{0}$ corresponds to the row-span of $G_{0}$ and form a group
under addition modulo two.

Furthermore, logical computation basis states can be represented by a $k$-bit
string $\mathbf{x}=(x_{1},\ldots,x_{k})$ as follows

$|\mathbf{x}\rangle_{L}=\frac{1}{\sqrt{|\mathcal{G}_{0}|}}\sum_{\mathbf{u}\in\mathbb{F}_{2}^{m}}|\mathbf{u}G_{0}+\mathbf{x}G_{1}\rangle,$
(389)

where $\mathbf{x}G_{1}$ is again obtained by matrix multiplication (modulo 2)
and is a constant shift identifying a coset of the group generated by addition
(modulo 2) of rows of $G_{0}$. We can compress this notation slightly by
noting

$\mathbf{u}G_{0}+\mathbf{x}G_{1}=(\mathbf{x},\mathbf{u})G,$ (390)

where $(\mathbf{x},\mathbf{u})$ is the row-vector resulting from joining
$\mathbf{u}$ and $\mathbf{x}$. Again, note that Eq. 390 should be read as
modulo two and this will be the convention for such expressions throughout the
remainder of this appendix.

The $j^{\mathrm{th}}$ logical $X$ operator, denoted $X_{Lj}$, ought to flip
the $|\mathbf{0}\rangle_{L}$ state to $|(\hat{\mathbf{o}}_{j})\rangle_{L}$
state, where $\hat{\mathbf{o}}_{j}$ is a unit vector with a single “1” entry
at the $j^{\mathrm{th}}$ location. It is straightforward to verify that
$X_{Lj}=X[\hat{\mathbf{o}}_{j}G_{1}]$ performs the required flip and that
$\hat{\mathbf{o}}_{j}G_{1}$ is equal to the $j^{\mathrm{th}}$ row of $G_{1}$.
Therefore, the logical operators of the code are given by the row vectors of
$G_{1}$. Furthermore, for every $\mathbf{g}$ in the row-span of $G_{0}$, the
operator $X[\mathbf{g}]$ is an $X$-stabilizer of the codespace, and this
enumerates all the $X$-stabilizers.

As a final notational preliminary, we will make use of a triple dot product
between triples of vectors. If $\mathbf{a}$, $\mathbf{b}$ and $\mathbf{c}$ are
binary vectors of equal length, we define

$|\mathbf{a}\wedge\mathbf{b}\wedge\mathbf{c}|=\sum_{j}a_{j}b_{j}c_{j}\pmod{2},$
(391)

which we again evaluate modulo 2. It is useful to note that this counts the
parity of the number of locations where operators $X[\mathbf{a}]$,
$X[\mathbf{b}]$ and $X[\mathbf{c}]$ all act non-trivially.

This $G$-matrix representation was also used for triorthogonal codes [37] and
quasi-triorthogonal codes [142, 141] except we are interested in different
transversality properties and so we will require different constraints on the
weight of rows in $G_{0}$ and $G_{1}$. The additional constraints determine
the transversality properties that we summarise with the following result,
which is a slight generalization (beyond $k=1$) of the proof techniques used
by Vasmer and Browne [146]

######  Lemma 1.

Let $\\{G^{A},G^{B},G^{C}\\}$ be a trio of $G$-matrices that represent a trio
of $[[n,k,d]]$ codes. Additionally, assume the following triple dot product
conditions (recall Eq. 391)

$|\hat{\mathbf{o}}_{p}G^{A}\wedge\hat{\mathbf{o}}_{q}G^{B}\wedge\hat{\mathbf{o}}_{r}G^{C}|=\begin{cases}1\mbox{
if }p=q=r\leq k\\\ 0\mbox{ otherwise }\\\ \end{cases}$ (392)

where $\hat{\mathbf{o}}_{p}$ is a binary unit vector with $1$ in location $p$
and 0 everywhere else. Then it follows that a physical $\text{CCZ}^{\otimes
n}$ realizes a transveral, logical $\text{CCZ}^{\otimes k}$.

Let us remark on what Eq. 392 means in terms of operators. Observe that when
$p\leq k$, the operator $X[\hat{\mathbf{o}}_{p}G^{D}]$ is the
$p^{\mathrm{th}}$ logical operator for codeblock $D\in\\{A,B,C\\}$. Therefore,
the condition of Eq. 392 tells us that the $X_{Lp}$ logical operators must
share an odd number of qubit indices where they all act non-trivially. All
other combinations of logical operators and stabilizers have an even number of
such locations.

###### Proof.

To determine the phase acquired from acting on the codespace with
$\text{CCZ}^{\otimes n}$, we first ask how this operator acts on an arbitrary
computational basis state. Recall that $\text{CCZ}^{\otimes
n}=\prod_{j=1}^{n}\text{CCZ}_{j}$ where $\text{CCZ}_{j}$ acts on qubit $j$ in
each block. Given a triple of $n$-qubit binary vectors $\mathbf{a}$,
$\mathbf{b}$ and $\mathbf{c}$ we have

$\text{CCZ}_{j}|\mathbf{a}\rangle|\mathbf{b}\rangle|\mathbf{c}\rangle=(-1)^{a_{j}b_{j}c_{j}}|\mathbf{a}\rangle|\mathbf{b}\rangle|\mathbf{c}\rangle,$
(393)

and so

$\text{CCZ}^{\otimes
n}|\mathbf{a}\rangle|\mathbf{b}\rangle|\mathbf{c}\rangle=(-1)^{|\mathbf{a}\wedge\mathbf{b}\wedge\mathbf{c}|}|\mathbf{a}\rangle|\mathbf{b}\rangle|\mathbf{c}\rangle,$
(394)

where $|\mathbf{a}\wedge\mathbf{b}\wedge\mathbf{c}|=\sum_{j}a_{j}b_{j}c_{j}$
as we introduced earlier. Next, we ask how this acts on the codespace.

Consider a trio of computational basis states
$|\mathbf{x}\rangle_{L}|\mathbf{y}\rangle_{L}|\mathbf{z}\rangle_{L}$ encoded
in blocks $A$, $B$ and $C$ respectively. Using Eqs. 389 and 390, we see that

$\displaystyle|\mathbf{x}\rangle_{L}|\mathbf{y}\rangle_{L}|\mathbf{z}\rangle_{L}$
(395)
$\displaystyle=2^{-3m/2}\sum_{\mathbf{u},\mathbf{v},\mathbf{w}\in\mathbb{F}_{2}^{m}}|(\mathbf{x},\mathbf{u})G^{A}\rangle|(\mathbf{y},\mathbf{v})G^{B}\rangle|(\mathbf{z},\mathbf{w})G^{C}\rangle.$

To determine the phase acquired from acting on
$|\mathbf{x}\rangle_{L}|\mathbf{y}\rangle_{L}|\mathbf{z}\rangle_{L}$ with
$\text{CCZ}^{\otimes n}$, we consider its action on each term in the
superposition using Eq. 394. Each term acquires a phase

$\displaystyle\text{CCZ}^{\otimes
n}|(\mathbf{x},\mathbf{u})G^{A}\rangle|(\mathbf{y},\mathbf{v})G^{B}\rangle|(\mathbf{z},\mathbf{w})G^{C}\rangle$
(396)
$\displaystyle=(-1)^{\lambda}|(\mathbf{x},\mathbf{u})G^{A}\rangle|(\mathbf{y},\mathbf{v})G^{B}\rangle|(\mathbf{z},\mathbf{w})G^{C}\rangle,$

where the phase exponent is

$\displaystyle\lambda$
$\displaystyle=|(\mathbf{x},\mathbf{u})G^{A}\wedge(\mathbf{y},\mathbf{v})G^{B}\wedge(\mathbf{z},\mathbf{w})G^{C}|.$
(397)

Using linearity of the triple dot-product and expanding the vectors in terms
of unit-vectors, e.g
$(\mathbf{x},\mathbf{u})=\sum_{p}\hat{\mathbf{o}}_{p}(\mathbf{x},\mathbf{u})_{p}$,
we have

$\displaystyle\lambda$
$\displaystyle=\sum_{p,q,r}(\mathbf{x},\mathbf{u})_{p}(\mathbf{y},\mathbf{v})_{q}(\mathbf{z},\mathbf{w})_{r}|\hat{\mathbf{o}}_{p}G^{A}\wedge\hat{\mathbf{o}}_{q}G^{B}\wedge\hat{\mathbf{o}}_{r}G^{C}|.$
(398)

Next, using the assumption of Eq. 392, we see almost all these terms vanish
except a few when $p=q=r\leq k$

$\displaystyle\lambda$ $\displaystyle=\sum_{p\leq
k}(\mathbf{u},\mathbf{x})_{p}(\mathbf{v},\mathbf{y})_{p}(\mathbf{w},\mathbf{z})_{p}.$
(399)

Notice that if $p\leq k$, $(\mathbf{x},\mathbf{u})_{p}=(\mathbf{x})_{p}$ since
$\mathbf{x}$ is length $k$. Therefore,

$\displaystyle\lambda$ $\displaystyle=\sum_{p\leq
k}(\mathbf{x})_{p}(\mathbf{y})_{p}(\mathbf{z})_{p}$ (400)