Fault-tolerant preparation of approximate GKP states

TL;DR

This paper develops a fault-tolerant protocol for preparing approximate GKP states, crucial for quantum error correction, and analyzes its implementation challenges and robustness under realistic noise conditions.

Contribution

It introduces a rigorous fault-tolerance definition for GKP state preparation and proposes a fault-tolerant phase estimation protocol suitable for circuit QED.

Findings

The protocol uses one flag qubit and selective acceptance to prevent large shift errors.

Analytic expressions describe effects of non-linear dispersive shift and Kerr non-linearity.

Numerical results show accepted states retain error correction capabilities despite noise.

Abstract

Gottesman-Kitaev-Preskill (GKP) states appear to be amongst the leading candidates for correcting errors when encoding qubits into oscillators. However the preparation of GKP states remains a significant theoretical and experimental challenge. Until now, no clear definitions for fault-tolerantly preparing GKP states have been provided. Without careful consideration, a small number of faults can lead to large uncorrectable shift errors. After proposing a metric to compare approximate GKP states, we provide rigorous definitions of fault-tolerance and introduce a fault-tolerant phase estimation protocol for preparing such states. The fault-tolerant protocol uses one flag qubit and accepts only a subset of states in order to prevent measurement readout errors from causing large shift errors. We then show how the protocol can be implemented using circuit QED. In doing so, we derive analytic expressions which describe the leading order effects of the non-linear dispersive shift and Kerr non-linearity. Using these expressions, we show that to mitigate the non-linear dispersive shift and Kerr terms would require the protocol to be implemented on time scales four orders of magnitude longer than the time scales relevant to the protocol for physically motivated parameters. Despite these restrictions, we numerically show that a subset of the accepted states of the fault-tolerant phase estimation protocol maintain good error correcting capabilities even in the presence of noise.

Figures

Tables

Table 1. Table 1: The first row displays the accepted measurement strings (set A 4 subscript 𝐴 4 A_{4} ) for the four round non-adaptive phase estimation protocol. If the protocol does not output an element of A 4 subscript 𝐴 4 A_{4} , the output is rejected and the protocol begins anew. The second row displays the largest shift difference that can occur when applying the phase estimated shift in the presence of a single measurement readout error. As can be seen, these are all less than π / 6 ≈ 0.295 𝜋 6 0.295 \sqrt{\pi}/6\approx 0.295 . The third row gives the probabilities of obtaining each state at the output of the phase estimation protocol.

Four round phase estimation protocol acceptance set $A_4$	1111	1010	0101	1000	0010
Largest shift difference	0.235	0.225	0.225	0.225	0.225
Probability of obtaining output state	0.1258	0.1287	0.1279	0.0505	0.0499

Table 2. Table 2: Parameter values for the two simulations of the non-adaptive noisy phase estimation protocol. The values for κ 𝜅 \kappa are chosen based on state of the art 3D cavities. The parameters in the second column results in a T 1 subscript 𝑇 1 T_{1} time given by T 1 = 2 ​ π κ = 10 ​ ms subscript 𝑇 1 2 𝜋 𝜅 10 ms T_{1}=\frac{2\pi}{\kappa}=10\text{ms} . For a resonator frequency f r = 7 ​ GHz subscript 𝑓 𝑟 7 GHz f_{r}=7\text{GHz} , the quality factor is given by Q = f r κ / ( 2 ​ π ) = 7 × 10 7 𝑄 subscript 𝑓 𝑟 𝜅 2 𝜋 7 superscript 10 7 Q=\frac{f_{r}}{\kappa/(2\pi)}=7\times 10^{7} . For all simulations, the squeezing level of the input squeezed vacuum state is chosen to be 0.2.

Parameter values	Simulation 1	Simulation 2
$p$	$5\times {10}^{-3}$	${10}^{-3}$
$\frac{\kappa }{2\pi }$	$1.59\times {10}^{-5}$MHz	$1.59\times {10}^{-6}$MHz
$\frac{{\gamma }_1}{2\pi }$	$1.06\times {10}^{-3}$MHz	$1.06\times {10}^{-4}$MHz
$\frac{{\gamma }_{\varphi }}{2\pi }$	$7.96\times {10}^{-4}$MHz	$7.96\times {10}^{-5}$MHz
$g$	$8.92$MHz	$5.09$MHz
$\Delta$	$0.32$GHz	$0.32$GHz
$K$ (Kerr)	${10}^{-4}$MHz	${10}^{-5}$MHz

Table 3. Table 3: The second row corresponds to the probabilities of obtaining the output states of A 4 subscript 𝐴 4 A_{4} for the parameters chosen from the second column of Table   2 conditioned on all flag measurements being 0. The third row is identical but for the parameters chosen from the third column of Table   2 . The probabilities are smaller for noisier circuits since the flag qubit has a higher chance of being measured as 1 causing the protocol to be aborted.

Four round phase estimation protocol acceptance set $A_4$	1111	1010	0101	1000	0010
Probability of obtaining output state noisy simulation 1	0.1297	0.0727	0.0737	0.0474	0.0429
Probability of obtaining output state noisy simulation 2	0.1468	0.0982	0.0978	0.0514	0.0448

Equations

$$|\overline{0}⟩=\sum _{k=-\infty }^{\infty }{S}_p^k|q=0⟩$$

$$|\overline{1}⟩=\sum _{k=-\infty }^{\infty }{S}_p^k|q=\sqrt{\pi }⟩,$$

$$|\text{sq}⟩=\int \frac{dq}{(\pi {\Delta }^2{)}^{1/4}}{e}^{-q^2/(2{\Delta }^2)}|q⟩.$$

$$|\overline{0}{⟩}_{\text{approx}}=C\sum _{k=-\infty }^{\infty }{e}^{-2\pi {\tilde{\Delta }}^2{k}^2}D(k\sqrt{2\pi })|\text{sq}⟩,$$

$$\sum _{k=-\infty }^{\infty }{\int }_{2\sqrt{\pi }k-\delta }^{2\sqrt{\pi }k+\delta }|⟨q|\tilde{0}⟩{|}^{2}dq>1-ϵ.$$

$$\sum _{k=-\infty }^{\infty }{\int }_{\sqrt{\pi }k-\delta }^{\sqrt{\pi }k+\delta }|⟨p|\tilde{0}⟩{|}^{2}dp>1-ϵ.$$

$$\sum _{k=-\infty }^{\infty }{\int }_{2\sqrt{\pi }k-\delta }^{2\sqrt{\pi }k+\delta }|⟨p|\widetilde{+}⟩{|}^{2}dp>1-ϵ.$$

$$\sum _{k=-\infty }^{\infty }{\int }_{\sqrt{\pi }k-\delta }^{\sqrt{\pi }k+\delta }|⟨q|\widetilde{+}⟩{|}^{2}dq>1-ϵ.$$

$$c_{\text{max}}^q=\underset{c}{\mathrm{argmax}}\sum _{k=-\infty }^{\infty }{\int }_{2\sqrt{\pi }k+c-\delta }^{2\sqrt{\pi }k+c+\delta }|⟨q|\tilde{0}⟩{|}^{2}dq.$$

$$c_{\text{max}}^p=\underset{c}{\mathrm{argmax}}\sum _{k=-\infty }^{\infty }{\int }_{\sqrt{\pi }k+c-\delta }^{\sqrt{\pi }k+c+\delta }|⟨p|\tilde{0}⟩{|}^{2}dp.$$

$$\displaystyle\gamma_{m}=\operatorname*{arg\,max}_{\gamma}{\sum_{x_{m}=0,1}\Big{|}\int d\theta e^{i\theta}P_{\gamma}(x[m]|\theta)\Big{|}}.$$

$$\displaystyle P_{\gamma=\gamma_{m}}(x[m]|\theta)=\prod_{i=1}^{m}\cos^{2}{\Big{(}\frac{\theta+\gamma_{i}}{2}+x_{i}\frac{\pi}{2}\Big{)}}.$$

$$\tilde{\theta }=\text{arg}{\int }_{-\pi }^{\pi }d\theta e^{i\theta }P(x[M]|\theta ).$$

$$H_s={\tilde{\omega }}_r{a}^{†}a+{\tilde{\omega }}_{a}Z+\chi a^{†}aZ-\varphi (a^{†}a{)}^{2}Z,$$

$$H_s={\tilde{\omega }}_r{a}^{†}a+{\tilde{\omega }}_{a}Z+\chi a^{†}aZ-\varphi (a^{†}a{)}^{2}Z-\frac{K}{2}(a^{†}a{)}^2.$$

$$H_d(t)=\mathcal{E}(t)a^{†}{e}^{-i{\omega }_{d}t}+{\mathcal{E}}^{\ast }(t)a{e}^{i{\omega }_{d}t},$$

$$H(t)=H_s+H_d(t).$$

$$V_R(0,T)$$

$$=R_{+}(T)D(A_{+})e^{i{B}_{+}}|0⟩⟨0|+R_{-}(T)D(A_{-})e^{i{B}_{-}}|1⟩⟨1|,$$

$$R_{\pm }(T)=e^{-iT{\omega }_{\pm }{a}^{†}a}(1\pm iT{\varphi }_{\pm }(a^{†}a{)}^2),$$

$$A_{\pm }=-i{\int }_0^T\mathcal{E}(t)e^{i{\omega }_{\pm }t}dt\mp {\varphi }_{\pm }(2a^{†}a-1){\int }_0^T\mathcal{E}(t)t{e}^{i{\omega }_{\pm }t}dt+\mathcal{O}({\varphi }_{\pm }^2),$$

$$B_{\pm }$$

$$\displaystyle\pm\frac{\phi_{\pm}(2a^{\dagger}a+1)}{2}\int_{0}^{T}dt\int_{t}^{T}dt^{\prime}\Big{(}\mathcal{E}^{*}(t^{\prime})\mathcal{E}(t)e^{i\omega_{\pm}(t^{\prime}-t)}+\mathcal{E}(t^{\prime})\mathcal{E}^{*}(t)e^{i\omega_{\pm}(t^{\prime}-t)}\Big{)}(t^{\prime}-t)+\mathcal{O}(\phi_{\pm}^{2}).$$

$$\displaystyle\mathcal{E}(t)\approx-2.09562\Big{(}\frac{\sqrt{2\pi}}{T}\Big{)}e^{-\frac{(t-\mu)^{2}}{2\sigma^{2}}},$$

$$A_{\pm }\approx \mp \frac{\sqrt{2\pi }}{2}\mp {\varphi }_{\pm }(2a^{†}a-1){\int }_0^T\mathcal{E}(t)t{e}^{i{\omega }_{\pm }t}dt.$$

$$\dot{\rho }=-i[H(t),\rho ]+\kappa \mathcal{D}[a]\rho +{\gamma }_1\mathcal{D}[{\sigma }_{-}]\rho +{\gamma }_{\varphi }\mathcal{D}[{\sigma }_z]\rho ,$$

$$\displaystyle|\psi^{(1)}\rangle=\frac{1}{\sqrt{2}}\Big{\{}|00\rangle|\varphi\rangle+\sqrt{1-p}|11\rangle D(\sqrt{2\pi})|\varphi\rangle+\sqrt{p}|01\rangle|\varphi\rangle\Big{\}}.$$

$$\displaystyle|\psi^{(2)}\rangle=\frac{1}{\sqrt{2}}\Big{\{}|00\rangle|\varphi\rangle+\sqrt{1-p}|10\rangle D(\sqrt{2\pi})|\varphi\rangle+\sqrt{p}|01\rangle|\varphi\rangle\Big{\}}.$$

$$\displaystyle|\psi^{(3)}\rangle=\frac{1}{\sqrt{2-p}}\Big{\{}|0\rangle|\varphi\rangle+\sqrt{1-p}|1\rangle D(\sqrt{2\pi})|\varphi\rangle\Big{\}}.$$

$$\displaystyle|\psi^{(4)}\rangle=\frac{1}{\sqrt{2(2-p)}}\Big{\{}|0\rangle(1+e^{i\gamma}\sqrt{1-p}D(\sqrt{2\pi}))|\varphi\rangle+|1\rangle(1-e^{i\gamma}\sqrt{1-p}D(\sqrt{2\pi}))|\varphi\rangle\Big{\}}.$$

$$\displaystyle|\psi^{\text{out,0}}\rangle=\frac{1}{\sqrt{2-p}}\Big{\{}(1+e^{i\gamma}\sqrt{1-p}D(\sqrt{2\pi}))|\varphi\rangle\Big{\}}.$$

$$\displaystyle|\psi^{\text{out,1}}\rangle=\frac{1}{\sqrt{2-p}}\Big{\{}(1-e^{i\gamma}\sqrt{1-p}D(\sqrt{2\pi}))|\varphi\rangle\Big{\}}.$$

Full text

I Introduction

Fault-tolerant quantum computing will be essential for implementing large scale quantum algorithms that offer provable speed-ups over the best known classical algorithms. Currently there are many proposals for encoding qubits into error correcting codes in order to perform universal fault-tolerant quantum computation. Depending on the underlying physical architecture, some encoding schemes are more suitable than others.

One method proposed by Gottesman, Kitaev and Preskill is to encode a qubit into an oscillator such that small shift errors in both position and momentum can be corrected. Although some bosonic codes have been designed to correct realistic errors arising from noise models encountered in the experiment (e.g. photon loss), recently it has been shown that GKP codes have better error correction capabilities than such codes under the assumption of perfect encoding and decoding AlbertBosonicCodePerformance18 ; GaussianChannel18 ; BarbaraDanielGKP . In addition, it has been shown how GKP codes can be concatenated with the toric code in order to achieve larger threshold values compared to toric codes with bare physical qubits ToricGKPBarbara ; GKPFaultToleranceJapan ; NohChambsGKPSurface . Lastly, given a supply of GKP-encoded Pauli eigenstates, universal fault-tolerant quantum computation can be achieved using only Gaussian operations MagicStateGKP19 .

Given the above, it is clear that the fault-tolerant preparation of encoded GKP states is an important problem that needs to be addressed. Various proposals for preparing GKP states have been outlined FirstGKPPrep02 ; GlancyKnillCondition ; BarbaraDanielGKP ; GKPStatePrepHome16 ; HomeNatureGKP ; TaiwanGKPPrep18 ; GKPStatePrepHomePhotonCatalysis ; ProbFaultGKP19 ; ShrutiStabilizedCat19 . However to our knowledge, no clear definitions for fault-tolerantly preparing GKP states using qubit-cavity couplings have been proposed. As such, without careful consideration, it is possible that a small number of faults lead to large uncorrectable shift errors. Inspired by AGP06 , in this work we propose new fault-tolerant definitions for preparing GKP states which tolerate small shift errors and a small number of faults occurring on ancilla qubits during the protocol. We then show how the phase estimation protocol proposed in this work satisfies our fault-tolerance criteria. In particular, the protocol is robust to a single fault occurring on the ancilla qubits in addition to shift errors of magnitude at most $0.06$. In order to be fault-tolerant, the protocol uses one flag qubit (and thus requires a total of two ancilla qubits) to prevent damping errors and imperfect implementations of the required gates from causing large shift errors. In addition, we provide an algorithm which only accepts a subset of all output states of phase estimation in order to prevent a single measurement readout error from causing large uncorrectable shift errors. We then proceed to show how our protocol can be implemented using circuit QED. We first analytically derive expressions describing the effects of the non-linear dispersive shift and Kerr non-linearity on the evolution of the cavity. We then numerically show that certain states output from the phase estimation protocol are robust to noise processes found in current 2D and 3D cavities since these can still correct small shift errors with high probability.

Our paper is structured as follows. In Section II, we provide new metrics which we use throughout the remainder of the paper to characterize the error correction capabilities of approximate GKP states. The fault-tolerance definitions used throughout this paper are given in Section III. In Section IV.1, we briefly review the phase estimation protocol described in BarbaraDanielGKP . In Section IV.2, we obtain a new phase estimation protocol and prove that it is fault-tolerant under our definitions. In Section IV.3, we compare the error correction capabilities of various states obtained from the phase estimation protocol in the noise free case. Section V is devoted to the implementation and error analysis of our protocol in circuit QED. In Section V.1, we provide analytic expressions for the time evolution of the qubit-cavity coupling when implementing a controlled-displacement gate. The expressions are derived in Appendix B. In Section V.2, we numerically solve a master equation which includes all considered noise processes, such qubit damping and dephasing, photon loss in addition to measurement, ancilla state-preparation and gate errors which arise from a depolarizing noise channel. In Section VI we summarize our results and discuss possible future directions.

II Goodness of approximate GKP states

As was explained in GKP2001 ; BarbaraDanielGKP , preparing perfect GKP states would require an infinite amount of squeezing. In a realistic setting, one can only prepare approximate GKP states with finite squeezing. Perfect GKP states, which are $+1$ eigenstates of the mutually commuting operators $S_{p}=e^{-2i\sqrt{\pi}p}$ and $S_{q}=e^{2i\sqrt{\pi}q}$, can correct shift errors of size at most $\frac{\sqrt{\pi}}{2}$. Note that for the displacement operator $D(\alpha)=e^{\alpha a^{\dagger}-\alpha^{*}a}$, we can write $S_{p}=D(\sqrt{2\pi})$ and $S_{q}=D(i\sqrt{2\pi})$. The goal of this paper will be to fault-tolerantly prepare approximate GKP states which can correct small shift errors correctable by perfect GKP states with high probability (in Section III we will specify what we mean by correctable shift errors). Therefore, it is important to have a metric which allows us to compute the "goodness" of an approximate GKP state.

Recall that for perfect GKP states, the logical $|\overline{0}\rangle$ and $|\overline{1}\rangle$ states are given by

[TABLE]

up to normalization. In practice, approximate GKP states analogous to those in Eq. 1 can be prepared by first preparing finitely squeezed states in $q$ space and approximately projecting these states onto the $S_{p}=1$ eigenspace. For instance, the $|q=0\rangle$ state can be written as a squeezed vacuum state $|\text{sq}\rangle$ with squeezing parameter $\Delta$ which is given by

[TABLE]

One can then apply a sum of displacements with a Gaussian filter to obtain

[TABLE]

where $C$ is a normalization coefficient. As was shown in BarbaraDanielGKP , it is natural to have $\tilde{\Delta}=\Delta$. In Section IV, we will present a fault-tolerant version of phase estimation for approximately projecting the state in Eq. 2 onto the $+1$ eigenspace of $S_{p}$ (see BarbaraDanielGKP which provides the first description for using phase estimation to prepare approximate GKP states).

Naturally because of the finite width of the peaks of approximate GKP states, it will not be possible to correct a shift error in $p$ or $q$ of magnitude at most $\frac{\sqrt{\pi}}{2}$ with certainty. For example, suppose we have an approximate $|\overline{0}\rangle$ GKP state with a peak at $q=0$ subject to a shift error $e^{-ivp}$ with $|v|\leq\frac{\sqrt{\pi}}{2}$. The finite width of the Gaussian peaks will have a non-zero overlap in the region $\frac{\sqrt{\pi}}{2}<q<\frac{3\sqrt{\pi}}{2}$ and $\frac{-3\sqrt{\pi}}{2}<q<\frac{-\sqrt{\pi}}{2}$. Thus with non-zero probability the state can be decoded to $|\overline{1}\rangle$ instead of $|\overline{0}\rangle$ (see Fig. 1 for an illustration).

In general, if an approximate GKP state is afflicted by a correctable shift error, we would like the probability of decoding to the incorrect logical state to be as small as possible. A smaller overlap of the approximate GKP state in regions in $q$ and $p$ space that lead to decoding the state to the wrong logical state will lead to a higher probability of correcting a correctable shift error by a perfect GKP state. These remarks motivate the following definition

Definition 1**.**

Let $|\tilde{0}\rangle$ be an approximate logical $|\overline{0}\rangle$ GKP state. We say that $|\tilde{0}\rangle$ is $(\frac{\sqrt{\pi}}{2}-\delta,\epsilon)_{q}$-GKP correctable if and only if for a given tuple $(\delta,\epsilon)$ with $\delta\leq\frac{\sqrt{\pi}}{2}$, and $0\leq\epsilon\leq 1$, we have

[TABLE]

We say that $|\tilde{0}\rangle$ is $(\frac{\sqrt{\pi}}{2}-\delta,\epsilon)_{p}$-GKP correctable if and only if for a given tuple $(\delta,\epsilon)$ with $\delta\leq\frac{\sqrt{\pi}}{2}$, and $0\leq\epsilon\leq 1$, we have

[TABLE]

Note that the bounds in Eqs. 4 and 5 are different since GKP states have peaks defined on a rectangular lattice. Similarly, for an approximate logical $|\overline{+}\rangle$ state, we have the following definition

Definition 2**.**

Let $|\tilde{+}\rangle$ be an approximate logical $|\overline{+}\rangle$ GKP state. We say that $|\tilde{+}\rangle$ is $(\frac{\sqrt{\pi}}{2}-\delta,\epsilon)_{p}$-GKP correctable if and only if for a given tuple $(\delta,\epsilon)$ with $\delta\leq\frac{\sqrt{\pi}}{2}$, and $0\leq\epsilon\leq 1$, we have

[TABLE]

We say that $|\tilde{+}\rangle$ is $(\frac{\sqrt{\pi}}{2}-\delta,\epsilon)_{q}$-GKP correctable if and only if for a given tuple $(\delta,\epsilon)$ with $\delta\leq\frac{\sqrt{\pi}}{2}$, and $0\leq\epsilon\leq 1$, we have

[TABLE]

As an example, suppose we have two approximate GKP states $|\tilde{0}\rangle_{1}$ and $|\tilde{0}\rangle_{2}$ which are $(\frac{\sqrt{\pi}}{2}-\delta,\epsilon_{1})_{q}$ and $(\frac{\sqrt{\pi}}{2}-\delta,\epsilon_{2})_{q}$ correctable for a shift of size $e^{i(\frac{\sqrt{\pi}}{2}-\delta)p}$ (which is correctable by a perfect $|\overline{0}\rangle$ GKP state). If $\epsilon_{1}<\epsilon_{2}$, then $|\tilde{0}\rangle_{1}$ will correct a shift of size $\frac{\sqrt{\pi}}{2}-\delta$ with greater probability than $|\tilde{0}\rangle_{2}$. This is due to the fact that $|\tilde{0}\rangle_{1}$ has a smaller overlap in regions which result in decoding the logical $|\overline{0}\rangle$ state to the logical $|\overline{1}\rangle$. In this sense we say that $|\tilde{0}\rangle_{1}$ is better than $|\tilde{0}\rangle_{2}$ at correcting shift errors in $q$ space.

III Fault-tolerant definitions

Given a protocol for preparing an approximate GKP state, errors that occur during the protocol can potentially accumulate resulting in a large shift error in addition to deforming the output state. This can then result in an output state which significantly differs from a good approximate GKP state. By good, we mean a state that for a desired value of $\delta$, the state is $(\frac{\sqrt{\pi}}{2}-\delta,\epsilon_{1})_{q}$ and $(\frac{\sqrt{\pi}}{2}-\delta,\epsilon_{2})_{q}$ correctable with $\epsilon_{1},\epsilon_{2}\ll 1$.

The desired values for $\delta$ depend on the particular fault-tolerant error correction protocol used to correct shift errors on encoded data qubits (see for instance CDT09 ; AliferisCross07 ; CJL16 ; CJL16b ; ChamberlandNeuralNet ). For example, one can use a version of the error correction scheme (which reduces to Steane error correction for qubit CSS stabilizer codes) as proposed by Glancy and Knill in GlancyKnillCondition . In this scheme, logical $|\overline{0}\rangle$ and $|\overline{+}\rangle$ ancillas are prepared and interact with the encoded data qubit via a CNOT gate to correct the shift errors afflicting the data qubit. It is shown that the largest shift error that can be corrected is $\frac{\sqrt{\pi}}{6}$. The threshold of $\frac{\sqrt{\pi}}{6}$ arises from how shift errors that are initially afflicting the encoded data qubit and ancilla states propagate through the CNOT gates and combine prior to the measurement. In practice, if additional shift errors occur during the the error correction scheme (say after applying the CNOT gates), the largest correctable shift error could potentially be smaller than $\frac{\sqrt{\pi}}{6}$ 111Using the optimizations considered in CDT09 , using Knill error correction could potentially increase the threshold of $\frac{\sqrt{\pi}}{6}$ to a larger value.. In what follows, we will assume that after preparing the desired ancilla states using some state preparation protocol, the ancillas will be used in the error correction scheme of GlancyKnillCondition (assumed to be fault-free) to correct shift errors on encoded data qubits. We also point out that due to the propagation of shift errors in the error correction scheme of GlancyKnillCondition , it is important to prepare approximate $|\tilde{0}\rangle$ and $|\tilde{+}\rangle$ states that have small shift errors in both $p$ and $q$.

If a small number of errors that occur during the preparation of an approximate GKP state result in large shift errors (or linear combinations of large shift errors) on the output state, then clearly the protocol used would not be practical. Thus it is important that a given state preparation protocol be fault-tolerant. In what follows we will define what we mean by fault-tolerant. We start with the following two definitions:

Definition 3**.**

A shift error is said to be correctable if the magnitude of the shift is less than or equal to $\frac{\sqrt{\pi}}{6}$. Otherwise, we will say that the shift error is uncorrectable.

Definition 4**.**

Suppose we have a protocol for preparing an approximate GKP state. We will say that the output state is an ideal approximate GKP state if no faults occur during the protocol.

Thus by Definition 4, any approximate GKP state obtained from a state preparation protocol will be called an ideal approximate GKP state if the protocol is implemented fault free, even if the output state is $(\frac{\sqrt{\pi}}{2}-\delta,\epsilon)_{p,q}$ correctable for some desired $\delta$ with large $\epsilon$ (so that the probability of correcting a shift $\frac{\sqrt{\pi}}{2}-\delta$ is small).

Note that the notion of correctable in Definition 3 assumes that only the data and ancilla qubits used in the error correction scheme of GlancyKnillCondition have shift errors. If other operations such as the CNOT gates introduce additional errors, the correctable threshold would be smaller than $\sqrt{\pi}/6$.

With the above definitions we are now ready to define what it means for a state preparation protocol of an approximate GKP state to be fault-tolerant. We point out that in Section IV we consider a fault-tolerant state preparation protocol based on phase estimation. Hence the definitions given below are specific to the case where an approximate GKP state is obtained by coupling a qubit to an oscillator.

Definition 5**.**

$(m,\tilde{\delta})$-Fault-tolerant state preparation of an approximate GKP state*:*

Suppose we have a protocol for preparing an approximate GKP state which is obtained by coupling qubits to a harmonic oscillator. Suppose also that at most $m$ faults occur during the protocol on the qubit Hilbert space and in addition, a correctable shift error in either $p$ or $q$ of size at most $\tilde{\delta}$ occurs on the oscillator Hilbert space. We will say that the protocol is an $(m,\tilde{\delta})$-Fault-tolerant state preparation of an approximate GKP state protocol if the output state differs from an ideal approximate GKP state by a correctable shift error.

A few clarifications are necessary. Firstly, a fault on the qubit Hilbert space corresponds to a location where an error can occur (see for instance AGP06 ). By location we are referring to a particular time step where either a gate is implemented, a qubit is prepared, a qubit is measured or a qubit is idling. On the oscillator Hilbert space, if an error occurs, we can always expand that error into shift errors (see for instance Eq. (7.12) in AlbertBosonicCodePerformance18 and also BarbaraDanielGKP ). By performing the error correction scheme of GlancyKnillCondition , measuring the $q$ and $p$ quadratures to perform error correction will always project the state onto a state with a single shift error in $q$ and $p$. Lastly, we point out that $\tilde{\delta}$ in Definition 5 relates to the largest allowed size of the shift error which occurs during the protocol. This should not be confused with $\delta$ in Definitions 1 and 2 which relates to the size of a shift error that can be corrected by an approximate GKP state222The probability $1-\epsilon$ of correcting a shift of size $\frac{\sqrt{\pi}}{2}-\delta$ depends on the location and width of the peaks..

Intuitively, a fault-tolerant protocol should have the property that if both the number of qubit errors and the size of shift errors are small, the resulting shift error on the output state should be correctable. Depending on the desired value for $\tilde{\delta}$, a particular fault-tolerant error correction protocol might be less tolerant to the size of the shift errors that occurred during the preparation of the approximate GKP state (in practice a different error correction scheme than Steane error correction could be used).

Suppose a state preparation protocol satisfies Definition 5. If during the preparation of the state, $m$ faults occur on the qubit Hilbert space in addition to a shift error of size at most $\tilde{\delta}$ on the oscillator Hilbert space, the remaining shift error on the output state will be corrected in a perfect version of the error correction protocol in GlancyKnillCondition . As we will see in Section IV, due to measurement errors, the phase estimation protocol presented in BarbaraDanielGKP needs to be modified in order to be a $(1,\tilde{\delta})$ fault-tolerant protocol. There are additional fault locations (apart from measurement errors) which can result in large shift errors that need to be treated with care.

Lastly we describe how we will evaluate the error correction properties of an approximate GKP state obtained from a noisy state preparation protocol. It is important to compare the goodness of a GKP state prepared from a noisy circuit to that of an ideal circuit. If the output state of a noisy state preparation protocol is correctable, the shift error will be removed when performing error correction. Therefore we will compare the $(\frac{\sqrt{\pi}}{2}-\delta,\epsilon)_{p,q}$ correctable properties of output states after performing an optimal shift back correction (as long as the shift error is correctable) to the output state. If the shift error is not correctable, the protocol will be deemed too noisy.

The optimal shift back correction is found as follows. In $q$ space, the optimal shift back $c^{q}_{\text{max}}$ is computed as

[TABLE]

Similarly, in $p$ space we have

[TABLE]

For protocols preparing $|\tilde{+}\rangle$, the metric can be defined analogously, but with the integral bounds for $p$ and $q$ switched.

IV Fault-tolerant phase estimation protocol

In Section IV.1 we will give a brief review of the phase estimation protocol presented in BarbaraDanielGKP . In Section IV.2, we will show how the protocol can be modified so that it becomes a $(1,\tilde{\delta})$-Fault-tolerant state preparation of an approximate GKP state and we will provide the value of $\tilde{\delta}$.

IV.1 Brief review of the phase estimation protocol for preparing approximate GKP states

The phase estimation circuit for preparing an approximate $|\tilde{0}\rangle$ GKP state is given in Fig. 2. The $H$ gate is the Hadamard gate and $\Lambda(e^{i\gamma})=\text{diag}(1,e^{i\gamma})$. After applying several rounds of the circuit in Fig. 2, the input squeezed vacuum state (given in Eq. 2) is projected onto an approximate eigenstate of $S_{p}$ with some random eigenvalue333The phase $\theta$ that is obtained depends on the measurement outcome of the ancilla qubit in each round. $e^{i\theta}$. Additionally, an estimated value for the phase $\theta$ is obtained. After computing the phase, the state can be shifted back to an approximate $+1$ eigenstate of $S_{p}$.

There are a variety of ways to choose the phases $\gamma$ at each round in the gate $\Lambda(e^{i\gamma})$. In a non-adaptive phase estimation protocol, half of the values for $\gamma$ can be chosen to be 0 and the other half can be chosen to be $\pi/2$. In an adaptive phase estimation protocol, if the phase to be estimated is $\theta$ and assuming no prior knowledge of $\theta$, during the first round the phase can be chosen as $\gamma_{1}=0$. For later rounds, it was shown in BarbaraDanielGKP that the next phases can be chosen as

[TABLE]

In Eq. 10, the probability $P_{\gamma}(x[m]|\theta)$ is the probability of obtaining measurement outcomes $x_{1},\cdots,x_{m}$ when the produced output eigenstate $|\psi_{\theta}\rangle$ has eigenvalue $e^{i\theta}$. Since the measurement results for each round are independent, the estimated probability is given by

[TABLE]

By analytically computing the evolution of the state, the exact probability of obtaining the outcomes $x[M]$ is derived in Appendix D.

Given the final measurement record $x[M]$, the estimated phase $\tilde{\theta}$ is chosen as

[TABLE]

The shift back correction based on the estimated phase $\tilde{\theta}$ is given by $e^{i\frac{\tilde{\theta}}{2\sqrt{\pi}}q}$. Note that for either the adaptive or non-adaptive protocol, the estimated phase and probabilities are computed using Eqs. 11 and 12.

IV.2 $(1,\tilde{\delta})$ fault-tolerant state preparation using phase estimation

From the circuit of Fig. 2 and from Eq. 12, a measurement readout error can result in the wrong estimated phase which in turn results in an incorrect shift-back correction (application of $e^{ivq}$ where $v$ is computed using Eqs. 10 and 11). Now suppose that the output state is afflicted by a shift error of the form $e^{-iwq}$. In this case, the final output state will have a total shift error of the form $e^{-i(v-w)q}$. If $|v-w|\mod\sqrt{\pi}>\sqrt{\pi}/6$, then the shift error will be uncorrectable and thus the phase estimation protocol will not be fault-tolerant 444Note that $|v-w|$ should be taken modulo $\sqrt{\pi}$ since a perfect $|\overline{0}\rangle$ GKP state in $p$ space has peaks at any integer multiple of $\sqrt{\pi}$.. In what follows we will identify an output state of the phase estimation protocol as $x[m]$ if it arises from the measurement outcomes $x_{1},\cdots,x_{m}$ in the fault-free case.

Consider the case where the phase estimation protocol is implemented in $m$ rounds and let us assume that a single measurement readout error occurs and that all other operations are fault free. In this case the output state $x_{\text{out}}[m]$ will have one bit which differs from the bit string $x_{\text{corr}}[m]$ that would have been obtained in the fault free case (so that the Hamming distance $d_{H}(x_{\text{out}}[m],x_{\text{corr}}[m])=1$). The shift-back operation of Eq. 10 applied to the output state will be the shift correction of an output state that is one Hamming distance away from the actual state that was produced. Thus to ensure that the protocol is $(1,\tilde{\delta})$ fault-tolerant for some $\tilde{\delta}$, we should only accept output states $x_{\text{out}}[m]$ which have the property that applying the shift back correction corresponding to any other state $x^{\prime}[m]$ with $d_{H}(x^{\prime}[m],x[m])=1$ results in a correctable left-over shift error. These remarks motivate the following protocol to prepare an approximate $|\tilde{0}\rangle$ state which is fault-tolerant to a single measurement readout error.

Calculation of acceptance set $A_{m}$ and $\tilde{\delta}$.

Consider all output states obtained from an $m$ round fault-free phase estimation protocol of Section IV.1. Let $A_{m}=\emptyset$ (which we call the acceptance set) and $\Gamma_{m}=\emptyset$. For each output state $x[m]$, do the following

Compute the shift correction $v_{m}=\frac{1}{2\sqrt{\pi}}\text{arg}\int_{-\pi}^{\pi}d\theta e^{i\theta}P(x[m]|\theta)$.

Let $j[m]$ be a bit string of size $m$. For all $j\in\{1,\cdots,m\}$ such that $d_{H}(x[m],j[m])=1$, compute the shift correction $v^{(j)}_{m}=\frac{1}{2\sqrt{\pi}}\text{arg}\int_{-\pi}^{\pi}d\theta e^{i\theta}P(j[m]|\theta)$.

If $\nexists\thinspace\thinspace j$ such that $|v_{m}-v^{(j)}_{m}|\mod\sqrt{\pi}>\frac{\sqrt{\pi}}{6}$, append $x[m]$ to $A_{m}$ and $\forall\thinspace j\in\{1,\cdots,m\}$, append $|v_{m}-v^{(j)}_{m}|$ to $\Gamma_{m}$.

If $A_{m}\neq\emptyset$, $\tilde{\delta}=\max_{m}\Gamma_{m}$.

$(1,\delta)$** fault-tolerant $m$ round phase estimation protocol for measurement readout errors.**

Consider the acceptance set $A_{m}\neq\emptyset$ and $\tilde{\delta}$ computed from the procedure described above. During an implementation of the phase estimation protocol subject to measurement readout errors, if the obtained output state $x[m]\notin A_{m}$, abort the protocol and start anew.

Note that during the protocol, there can be more than one measurement readout error. However if more than one readout error occurs, it is possible that the output state is afflicted by a shift error of magnitude greater than $\sqrt{\pi}/6$. Our protocol only guarantees protection against a single readout error.

As an example, consider the four round phase estimation protocol. Applying the procedure described above for the adaptive phase estimation protocol, we find that $A_{4}$ is empty when choosing the initial phase to be $\gamma_{0}=0$ or $\gamma_{0}=\frac{\pi}{2}$. This indicates that the adaptive phase estimation protocol of four rounds described in Section IV.1 is not fault-tolerant to single measurement readout errors555However it is possible that the adaptive phase estimation protocol described in Section IV.1 could be modified to be fault-tolerant to measurement readout errors. We leave such analysis to future work.. If the phases are computed using the non-adaptive phase estimation protocol and applying the protocol described above, we find that $\tilde{\delta}=\frac{\sqrt{\pi}}{6}-0.235\approx 0.0604$. The set $A_{4}$ of accepted measurement strings is given in Table 1 and has $|A_{4}|=5$. The total probability of obtaining a state in $A_{4}$ is roughly $48.3\%$. For the phase estimation protocol with more than four rounds, the acceptance set and the acceptance probability is very sensitive to the initial phase $\gamma_{0}$ and the domain of the optimized $\gamma$. For example, if we choose the initial phase $\gamma_{0}$ to be $\frac{\pi}{2}$ and the range of $\gamma$ to be $[0,2\pi]$, we get 18 accepted states (for the adaptive protocol). In this case the total probability of obtaining a state in $A_{8}$ is $6.25\%$. If we choose the initial phase $\gamma_{0}$ to be 0, the acceptance set has four states and the probability of obtaining a state in $A_{8}$ is $1.3\%$. For the non-adaptive phase estimation protocol, $|A_{8}|=66$ but the total probability of obtaining a state in $A_{8}$ is roughly $79.2\%$. These results indicate that although the adaptive phase estimation protocol outperforms the non-adaptive protocol in the fault-free case (BarbaraDanielGKP ; DispSensorBarbara ), the adaptive phase estimation protocol cannot be used in the presence of measurement readout errors for four rounds, and has significantly fewer states in $A_{8}$ for eight rounds. Therefore in a noisy implementation of phase estimation, the non-adaptive protocol is preferable. A list of the states belonging to the acceptance set $A_{8}$ for both the adaptive and non-adaptive protocols can be found at https://github.com/godott/GKP_phase_estimation.git.

Suppose now that for $m$ rounds the set $A_{m}$ is not empty. We then know there exists a $\tilde{\delta}$ such that the protocol is a $(1,\tilde{\delta})$ fault-tolerant $m$ round phase estimation protocol for measurement readout errors. However there are other fault locations where a single fault resulting in a qubit error could potentially lead to an uncorrectable shift error on output states. Note that an $X$ error prior to applying the controlled-$D(\sqrt{2\pi})$ gate will do nothing since the qubit is in the $|+\rangle$ state. A $Z$ error will just change the sign of one of the peaks of the output state in $p$ space. However a damping event occurring on the ancilla qubit before or during the application of the controlled-$D(\sqrt{2\pi})$ gate can result in incorrectly applying the $D(\sqrt{2\pi})$ displacement on the cavity. If several damping events occur during the protocol, the output state in $p$ space could potentially be badly deformed. In BarbaraDanielGKP it was proposed to replace the ancilla qubit by a $k$-qubit cat state so that if a single qubit undergoes damping, a single shift error of size $D(\frac{\sqrt{2\pi}}{k})$ would occur which can be made small for large $k$. However this approach would require the fault-tolerant preparation of a large $k$-qubit cat state which would significantly increase the required resources in addition to adding substantially more locations where errors can occur. Similar issues were observed in Rosenblum266 for preparing cat codes. However large cavity displacements were mitigated by interacting the cavity with a three level system $|f\rangle$, $|e\rangle$ and $|g\rangle$ instead of a qubit. A transition from $|f\rangle$ to $|e\rangle$ would not cause the incorrect gate from being applied. Thus two damping events would be required to cause the incorrect gate from being applied.

Here we propose an alternative solution for mitigating damping errors during the application of the controlled-$D(\sqrt{2\pi})$ gate which uses a single extra flag qubit CR17v1 ; CR17v2 ; CB17 ; TCD18Flag ; ReichardtFlag18 ; CC18MagicFlag . Consider the circuit in Fig. 3. If a bit flip error occurs before or during the application of the controlled-$D(\sqrt{2\pi})$ gate, the flag qubit will be measured as 1 instead of 0. In such a case the strategy will be simply to abort the protocol and start anew. The analysis for a general damping event using the flag qubit is given in Appendix A. It is shown that if ever a damping event on the $|+\rangle$ ancilla qubit causes the wrong $D(\sqrt{2\pi})$ gate to be applied, the flag qubit will be measured as 1 instead of 0 in which case we abort the protocol and start anew. Thus with the flag qubit, the phase estimation protocol is robust to a damping error during the application of the controlled displacement gate. Note however that if a fault occurs in addition to a damping event, the shift error resulting from the damping event can potentially go undetected. For instance, if in addition to damping, the flag qubit is subject to a measurement error, the measurement outcome can be 0 instead 1 resulting in acceptance when the protocol should have been aborted.

Suppose that only a single damping event occurs during the four round phase estimation protocol, and that all other operations are fault-free. Using the analytic expressions derived in Appendix A and for a damping rate of $p=0.5$, we performed simulations to numerically characterize the effects of qubit damping on the prepared GKP state. We found that the shift error resulting from a single damping event was negligible compared to the largest shift error resulting from a single measurement readout error.

Lastly, an $X$ error prior to the $\Lambda(e^{i\gamma})$ gate will change the sign of $\gamma$. After propagating through the Hadamard gate, the $X$ error will become a $Z$ error which will not affect the measurement outcome of the ancilla qubit. For the four round non-adaptive phase estimation protocol, we performed a simulation which showed that the shift error resulting from a single $X$ error prior to the $\Lambda(e^{i\gamma})$ gate is much smaller than the largest shift error arising from a measurement readout error on the ancilla qubit.

Let us assume that the noise model afflicting the oscillator can cause a shift error of size at most $\tilde{\delta}$ in $p$ space (i.e. a shift of the form $e^{-i\tilde{\delta}q}$). From the above, the largest shift error that can arise from a single fault afflicting the qubit space during the four round fault-tolerant phase estimation protocol is $0.235<\sqrt{\pi}/6$. Following Definition 5, we conclude that the protocol described in this section is a $(1,\tilde{\delta})$-fault-tolerant state preparation of an approximate logical $|\overline{0}\rangle$ GKP state with $\tilde{\delta}=\sqrt{\pi}/6-0.235\approx 0.06$.

In Section V we will discuss a physical implementation of the phase estimation protocol using circuit QED. A much more detailed analysis of the noise afflicting the cavity and the resulting shift errors will be provided.

IV.3 Goodness of approximate $|\tilde{0}\rangle$ states obtained from the noise free phase estimation protocol

In Section V.2, the fault tolerant implementation of phase estimation described in this section is analyzed for a noise model which introduces gate noise, measurement readout errors and ancilla state preparation errors in addition to photon loss, amplitude damping and dephasing. Here we consider a noiseless implementation of the four round non-adaptive phase estimation protocol described in this section which will be used to benchmark the noisy implementation.

Since the shift correction obtained from non-adaptive phase estimation is in $p$ space, for a range of values for $\delta$, we computed $\epsilon$ using the integral in Definition 1 after performing a shift back correction using Eq. 9. This allowed us to compute the probability of correcting shift errors in $p$ of size $\frac{\sqrt{\pi}}{2}-\delta$. Plots for the states 0101 and 1000 which belong to $A_{4}$ (see Table 1) are given in Fig. 4. The $(\frac{\sqrt{\pi}}{2}-\delta,\epsilon)_{p}$ correctable properties of the other states in $A_{4}$ are similar to the ones shown in Fig. 4. It can be seen that for small values of size $\frac{\sqrt{\pi}}{2}-\delta$, both states can correct the shift error with probability close to 1.

V Circuit QED implementation

V.1 State evolution in the dispersive regime

In this section we will describe a direct implementation of the controlled-$D(\sqrt{2\pi})$ gate. From BlaisFirstPaper ; BlaisSecondPaper , the Hamiltonian describing the coupling between a qubit and a cavity in the dispersive regime is given by

[TABLE]

where $\phi=\frac{g^{4}}{\Delta^{3}}$, $\chi=\frac{g^{2}}{\Delta}-\phi$, $\tilde{\omega}_{r}=\omega_{r}+\phi$ and $\tilde{\omega}_{a}=\frac{\omega_{a}+\chi}{2}$. Here $g$ is the coupling strength between the qubit and the cavity and $\Delta=|\omega_{r}-\omega_{a}|$. The term $\phi(a^{\dagger}a)^{2}Z$ corresponds to the non-linear dispersive shift. The Hamiltonian in Eq. 13 can be derived by performing an exact diagonalization of the Jaynes-Cummings Hamiltonian and keeping only leading order terms in $\phi$ BlaisSecondPaper . Note that a more systematic treatment of the qubit as an anharmonic oscillator leads to an additional term in Eq. 13 given by $-\frac{K}{2}(a^{\dagger}a)^{2}$ which is referred to as the Kerr non-linearity OptimalControl2 . Hence, in our analysis, we choose the system Hamiltonian to be given by

[TABLE]

The direct implementation of the controlled-$D(\sqrt{2\pi})$ gate can be achieved using the drive Hamiltonian

[TABLE]

where $\mathcal{E}(t)$ describes the pulse shape of the drive and $\omega_{d}$ is the drive frequency. Thus the total Hamiltonian describing the evolution of the qubit-cavity system during the implementation of the controlled-$D(\sqrt{2\pi})$ gate is given by

[TABLE]

Going into the rotating frame of the qubit and the cavity and defining $\phi_{\pm}=\phi\pm K$ and $\omega_{\pm}=\pm\chi$, in Appendix B we show that

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

In deriving Eqs. 18, 19 and 20, we kept only leading order terms in $\phi_{\pm}$ since the Hamiltonian in Eq. 14 neglects higher order terms in $\phi$ and $K$. We keep only leading order terms since with current experimental parameter values, $\phi_{\pm}$ is roughly three orders of magnitude smaller than $\chi$.

Firstly, notice that the terms $R_{\pm}(T)$ introduce a relative rotation between the $|0\rangle$ and $|1\rangle$ state of the qubit. Neglecting the terms proportional to $\phi_{\pm}$, it was shown in BarbaraDanielGKP that by choosing the total interaction time to be $T=\pi/\chi$, the relative rotation between $|0\rangle$ and $|1\rangle$ can be set to one. However due to the presence of the non-linear dispersive shift and Kerr terms in Eq. 18, it is clear that the relative rotation cannot be completely removed. Further, $R_{\pm}(T)$ does not depend on the pulse shape $\mathcal{E}(t)$ of the drive term. Therefore even with an optimized pulse shape (see below), the effects from the non-linear dispersive shift and the Kerr will cause an undesired relative rotation between the qubit states. Regardless, we will still choose the interaction time $T$ to be $\pi/\chi$ to ensure that we eliminate the relative rotation from the leading order terms in Eq. 18. In Section B.4, we provide further details using the analytic expressions to show that to mitigate the effects due to the non-linear dispersive shift and Kerr terms would require the protocol to take place on time scales of order $\frac{1}{\phi_{\pm}}$ which (for the parameter values in Table 2) is four orders of magnitude larger than $\frac{\pi}{\chi}$. For such long time scales, noise terms such as photon loss, dephasing and damping would render the protocol impractical.

The terms in Eq. 17 that perform the desired controlled displacement gate are given by $A_{\pm}$ in Eq. 19. The goal is to choose a pulse shape that implements the desired gate while at the same time minimizes the contributions from the non-linear dispersive shift and the Kerr term (terms proportional to $\phi_{\pm}$ in Eq. 19). In order to be experimentally relevant, it is important to choose a pulse shape that is accessible to near term experiments using 2D and 3D cavities. We chose a Gaussian pulse with the following parameters

[TABLE]

where $\mu=\frac{\pi}{2\chi}$ and $\sigma=\frac{\pi}{10\chi}$. With this Gaussian pulse we obtain

[TABLE]

The amplitude of the Gaussian, given by -2.09562\Big{(}\frac{\sqrt{2\pi}}{T}\Big{)}, is chosen numerically to ensure that the appropriate gate is being performed.

Since both the $A_{+}$ and $A_{-}$ terms are symmetric with opposite signs, the gate D\Big{(}-\sqrt{\frac{\pi}{2}}\Big{)} in Fig. 3 is not required. Furthermore, the first term in Eq. 22 is independent of $\chi$. However the contribution from the second term in Eq. 22 will depend on the coupling strength and relative frequencies between the qubit and the cavity. Thus reducing the value of $\chi$ will minimize the undesired effects arising from the non-linear dispersive shift and the Kerr term while still producing the desired gate.

The $B_{\pm}$ terms of Eq. 20 result in a phase difference between the $|0\rangle$ and $|1\rangle$ qubit states. However, computing the first integrals in Eq. 20 (i.e. terms before $\phi_{\pm}$), we find that for our chosen pulse, the phase difference remains fixed during every round of the phase estimation protocol. Therefore after applying the controlled displacement gate, we apply an additional phase gate which removes the phases introduced by the left most integrals of $B_{\pm}$. Note however that the $\phi_{\pm}$ terms in $B_{\pm}$ will not be cancelled and will thus introduce additional shift errors.

V.2 Full noise analysis and master equation results.

In this section we perform a numerical analysis for the noisy implementation of the phase estimation protocol described in Section IV.2. The simulation is performed in several steps that we now describe. First, the controlled displacement gate is modeled using the following master equation

[TABLE]

where $H(t)$ is given by Eq. 16 after going into the rotating frame and $\mathcal{D}[L]\rho=(2L\rho L^{\dagger}-L^{\dagger}L\rho-\rho L^{\dagger}L)/2$. The density matrix corresponds to the joint state of the cavity, ancilla qubit and flag qubit. The parameters $\kappa$, $\gamma_{1}$ and $\gamma_{\phi}$ correspond to the photon loss rate, the qubit decay rate and qubit dephasing. The pulse shape of the drive is given by Eq. 21. The total evolution time of the controlled displacement is given by $\frac{\pi}{\chi}$ with $\chi=\frac{g^{2}}{\Delta}-\phi$.

Both before and after the controlled displacement gate, the state of the cavity is subject to photon loss and its evolution is computed by solving a master equation. Based on current gate, state preparation and measurement times, we chose the evolution time from the preparation of the ancilla to the first CNOT gate (after which the controlled-displacement gate is performed) to be $0.14\mu$s. Similarly, the evolution time after the controlled-displacement gate to the time the ancilla is measured is also chosen to be $0.14\mu$s. Thus during one round, the cavity freely evolves with photon loss for a total of $0.28\mu$s. In addition, before and after the controlled displacement gate, we allow all qubit locations to fail with the following depolarizing noise model

With probability $p$, each two-qubit gate is followed by a two-qubit Pauli error drawn uniformly and independently from $\{I,X,Y,Z\}^{\otimes 2}\setminus\{I\otimes I\}$. 2. 2.

With probability $\frac{2p}{3}$, the preparation of the $|0\rangle$ state is replaced by $|1\rangle=X|0\rangle$. 3. 3.

With probability $\frac{2p}{3}$, any single qubit measurement has its outcome flipped. 4. 4.

With probability $p/10$, each Hadamard gate is followed by a Pauli error drawn uniformly and independently from $\{X,Y,Z\}$.

We chose $p/10$ for Hadamard gate failures since for current superconducting architectures, single qubit gate fidelities are about an order of magnitude higher than two-qubit gate fidelities. In practice, the gate errors applied to the qubit Hilbert space will depend strongly on the circuit QED architecture and should also be modeled using a master equation as in Eq. 23. However, we chose a depolarizing model in order to reduce the computation time and simplify the analysis. We also mention that in our analysis we assumed that $g$ is tunable. Thus both before and after the controlled-displacement gate, the cavity and qubit system is decoupled and can be treated separately.

The master equation was numerically solved using Qutip JOHANSSON20131234 , our code can be accessed at https://github.com/godott/GKP_phase_estimation.git. Due to the long computation time during the controlled displacement gate, performing a full Monte Carlo simulation to take into account gate errors with the depolarizing model was unfeasible (i.e. gate locations both before and after the controlled displacement gate). Instead, we analytically computed the error probabilities for each Pauli operator (using the depolarizing noise model described above) immediately after the first CNOT gate, before both measurement locations and before the phase gate666The error probabilities were computed by considering all possible single-fault locations resulting in a given error at the considered location. For instance, a $Z\otimes I$ error after the first CNOT gate can arise from a $Z\otimes I$ error from the faulty CNOT gate, but also from a $Z$ error after the application of the Hadamard gate prior to the CNOT gate.. At each location, all possible Pauli operators based on their associated probabilities were added. For a given Pauli error, the probabilities (which are expressed as functions of $p$) were then used (in addition to the state evolution obtained from the master equation) to compute the final probability of obtaining a given output state. Instances with two or more faults occurring on the qubit Hilbert space before and after the controlled-displacement gate were neglected. However our analysis is still more complete than previous implementations which considered only measurement errors and errors during the controlled-displacement gate. More details of the Pauli simulation can be found in Appendix D.

We performed two different simulations where for each simulation, the chosen parameter values are given in Table 2. The parameters chosen for the first simulation (middle column in Table 2) are based on current experimental values for 2D and 3D cavities GKPparams1 ; Rosenblum266 ; GKPparams2 ; LongCavityLifetimes ; 2DCavityLifetimes ; 2DCavityLifetimesv2 . The parameters chosen for the second simulation are based on values that might be obtained with improved future technologies. Plots showing the probability of correcting shift errors $\frac{\sqrt{\pi}}{2}-\delta$ after performing a shift correction of the output states of the noisy phase estimation protocol using Eq. 9 are given in Fig. 5. In what follows we will refer to these plots as $\epsilon-\delta$ plots. Note that noise during the phase estimation protocol introduces more shift errors in $p$ space. Therefore in Fig. 5, only $\epsilon-\delta$ plots in $p$ space are shown.

For the four round phase estimation protocol, we find that the state 1010 has a similar $\epsilon-\delta$ plot to the one for the state 0101 whereas the state 0010 has a similar $\epsilon-\delta$ plot to the one for the state 1000. The probability of obtaining each state for the two noisy simulations (parameters in Table 2) are given in Table 3. It can be seen that the probability of the state 0101 to correct shift errors is significantly more affected by the noise than the state 1000 or 1111. To understand why this is the case, it is useful to look at the wave function densities for these three states in the $q$ basis when no noise is present (see Fig. 6). Comparing the wave function densities, it can be observed that the state 1000 has three peaks with similar amplitudes, whereas the states 0101 and 1111 have two smaller peaks compared to the peak at the center. Performing a numerical analysis, it is found that when only damping is present (with all other noise terms including Kerr and non-linear dispersive shift set to zero), damping has a negligible effect on the height of the peaks for all considered states. However performing a numerical simulation with only the Kerr and non-linear dispersive shift terms present, these terms significantly reduce the height of the peaks for the state 0101 and 1010. Therefore, these states are left with one large peak in $q$ with all other peaks close to zero which results in a wave function density with very low resolution in $p$ space (which thus affects the ability for these states to correct shift errors in $p$). Interestingly, the state 1111 is more robust to the Kerr and non-linear dispersive shift contributions since there is a much smaller reduction in the amplitudes of the smaller peaks compared to those for the states 0101 and 1010 (see Fig. 7). Furthermore, the states 1000 and 0010 have three large peaks in $q$ and thus even with a reduced amplitude, these states have a higher wave function density resolution in $p$ space.

VI Conclusion

In this work we presented a fault-tolerant state preparation protocol for preparing GKP states using phase estimation. In Section II we provided metrics for comparing how good approximate GKP states are at correcting shift errors in both $q$ and $p$ space. In Section III, we provided a definition for the $(m,\tilde{\delta})$-Fault-tolerant state preparation of an approximate GKP state where $m$ is the maximum number of allowable faults which can occur on the qubit Hilbert space and $\tilde{\delta}$ is the maximum allowed shift error that can affect the oscillator. Using a non-adaptive phase estimation protocol with one ancilla qubit and one flag qubit, in Section IV we showed how the protocol can be made into a $(1,0.06)$-fault-tolerant state preparation of an approximate GKP state. The flag qubit is used to detect damping errors. In addition, it was shown how the adaptive phase estimation protocol of Section IV.1 can not be made fault-tolerant in the presence of measurement readout errors. For the four round non-adaptive phase estimation protocol, 5 of the 16 output states can be used in the presence of measurement readout errors. The total probability of obtaining the accepted states is approximately 0.48. In Section V, we considered how the protocol can be implemented in circuit QED. We first provided (to leading order) analytic expressions of the non-linear dispersive shift and Kerr terms during the evolution of the qubit and cavity in the fault-tolerant phase estimation protocol. We used these expressions to find a Gaussian pulse shape that allows one to implement the desired gates. However to mitigate effects due to the non-linear dispersive shift and Kerr terms would require the protocol to be implemented on time scales four orders of magnitude larger than those that were considered in this work. Due to the noise processes afflicting the system, such long time scales would render the protocol impractical. Performing two different simulations for both current and futuristic parameter values found in 2D and 3D cavities, we numerically solved a master equation to study the affects of qubit damping, dephasing and photon loss on the cavity in addition to gate and measurement errors during the protocol. It is shown numerically that 3 of the 5 accepted states (for the four round non-adaptive phase estimation protocol) are much more robust to noise arising from the non-linear dispersive shift and Kerr terms and maintain good error correction capabilities even in the presence of noise.

The pulse shape used in this work to implement the controlled-displacement gate of the protocol was obtained from the analytical expressions describing the time evolution of the qubit-cavity system. An important direction for future work is to find a pulse shape using methods such as optimal control OptimalControl1 ; OptimalControl2 in order to obtain a pulse which can further mitigate the effects from the non-linear dispersive shift and Kerr terms. Since the non-linear dispersive shift and Kerr terms are the dominant source of noise that reduce the error correcting capabilities of approximate GKP states, using optimal control could potentially allow the protocol to be implemented using a larger number of rounds in order to obtain better approximate GKP states.

The fault-tolerant state preparation of approximate GKP states presented in this work is tailored to protocols that use phase estimation. An interesting direction for future work would be to find fault-tolerant implementations for preparing approximate GKP states that apply to broader schemes such as those found in GKPStatePrepHome16 ; GKPStatePrepHomePhotonCatalysis . In addition, fault-tolerant state preparation protocols for hexagonal GKP codes could be analyzed since these offer better error correction capabilities than GKP codes on a square lattice. It would also be interesting to extend the ideas presented in this work beyond state preparation. For instance, fault-tolerant protocols for the implementation of logical gates would be of great interest.

VII Acknowledgements

We thank Barbara Terhal, Daniel Weigand and Daniel Gottesman for useful feedback and discussions. We give special thanks to Zlatko Minev for the many useful discussions regarding circuit QED. Y.S. performed the majority of the numerical simulations. C.C. performed the majority of the analytical calculations, developed the fault-tolerant schemes and wrote the majority of the manuscript. Y.S. and A.C. acknowledges scientific discussions with Fred Chong that took place during the completion of this project. Y.S. is funded in part by EPiQC, an NSF Expedition in Computing, under grant CCF-1730449 and NSF award STAQ, PFCQC-1818914. Y.S. is also funded in part by the NSF QISE-NET fellowship under grant number 1747426. This work was completed in part with resources provided by the University of Chicago Research Computing Center.

Appendix A Shift error arising from amplitude damping before the controlled-displacement gate

Consider amplitude damping on the $|+\rangle$ ancilla state before the controlled displacement gate as shown in Fig. 8. Let $|\varphi\rangle$ be the input state part of the cavity Hilbert space. Defining $p$ as the damping rate, after the controlled-displacement gate, the state of the system can be written as

[TABLE]

Here we have omitted writing the state of the environment interacting with the first ancilla qubit. However, this does not affect the result of the calculations that follow.

After the second CNOT gate, the state becomes

[TABLE]

If the flag qubit is measured as 1, the protocol is aborted. So assume that the flag is measured as 0. In this case the state becomes (tracing out the flag qubit)

[TABLE]

Applying the remaining gates in Fig. 8, the final state prior to the measurement of the ancilla qubit is

[TABLE]

If the measurement of the ancilla is 0, the output state will be

[TABLE]

If the measurement of the ancilla is 1, the output state will be

[TABLE]

Both outcomes occur with probability $1/2$.

Appendix B Analytic derivation of the unitary operator describing the evolution of the qubit-cavity system during the implementation of the controlled-$D(\sqrt{2\pi})$ gate.

B.1 Implementation of the controlled-$D(\sqrt{2\pi})$ gate in the lab frame.

In this section, we will derive the unitary operator describing the evolution of the qubit-cavity state during the application of the control-displacement gate. In the dispersive regime, the qubit-cavity interaction can be described as:

[TABLE]

where

[TABLE]

and

[TABLE]

In Eq. 31, $\phi=\frac{g^{4}}{\Delta^{3}}$, $\tilde{\omega}_{r}=\omega_{r}+\phi$, $\chi=\frac{g^{2}}{\Delta}-\phi$ and $\tilde{\omega}_{a}=\frac{\omega_{a}+\chi}{2}$. The term $\phi(a^{\dagger}a)^{2}Z$ corresponds to the non-linear dispersive shift. Note that we have not included the Kerr term $-\frac{K}{2}(a^{\dagger}a)^{2}$. At the end of this section we will modify our results to take into account its effect. We also note that in writing the drive Hamiltonian in Eq. 32, we neglected a term of the form $\lambda X$ (where $\lambda=(g/\Delta$) and all higher powers in $\lambda$. For parameter values considered in this paper, we found the effect of this term to be negligible. Additionally, we chose a pulse shape which is real valued.

In Eq. 32, we represent the drive pulse $\mathcal{E}(t)$ as

[TABLE]

where $\omega_{d}=\tilde{\omega}_{r}-\chi$ is the drive frequency.

Since the Hamiltonian does not commute at different times, the unitary operator describing the time evolution under the Hamiltonian of Eq. 30 is given by

[TABLE]

The right-hand side of Eq. 34 can be computed using the Suzuki-Trotter decomposition, so that

[TABLE]

where $\tilde{t}\equiv T/n$, $\omega_{\pm}\equiv\tilde{\omega_{r}}\pm\chi$.

Now, the two terms on the right-hand side of Eq. 35 can be decomposed as

[TABLE]

where

[TABLE]

and

[TABLE]

Hence from the above, we see that we need to compute terms of the form

[TABLE]

In order to compute the products appearing in Eq. 39, we will use the identity $Ae^{B}A^{-1}=e^{ABA^{-1}}$. Defining,

[TABLE]

and

[TABLE]

with

[TABLE]

we have that

[TABLE]

The term $e^{iH^{\prime}\tilde{t}}(a+a^{\dagger})e^{-iH^{\prime}\tilde{t}}$ in Eq. 43 can be computed using Heisenberg’s equation of motion. We obtain

[TABLE]

with

[TABLE]

where we defined

[TABLE]

We thus have the following set of coupled differential equations:

[TABLE]

and

[TABLE]

The solutions to Eqs. 47 and 48 are given by

[TABLE]

and

[TABLE]

To see this, we take a derivative of Eq. 49 to obtain

[TABLE]

Comparing Eqs. 47 and 51, we must have that

[TABLE]

Using Eqs. 49 and 50 and defining $H_{te}\equiv-kt(\phi^{\prime}+2\phi(a^{\dagger}a-1))$, we have that

[TABLE]

as desired. A similar calculation can be done for $a^{\dagger}(\tilde{t})$. Hence we have that

[TABLE]

for

[TABLE]

Writing $V(0,T)$ (Eq. 35) as

[TABLE]

using Eq. 54 and the fact that $R_{n}^{n}=e^{-iT(\omega_{+}-\phi a^{\dagger}a)a^{\dagger}a}\equiv R_{+}(T)$, we can write

[TABLE]

where we defined

[TABLE]

Notice that we can write $V_{+}(0,T)$ as products of displacements $D(A_{k})$. However, now $A_{k}$ is an operator instead of a complex number. In order to compute the products in Eq. 57, we will need to obtain a relation for terms of the form $D(A_{k})D(A_{j})$ (where $A_{k}$ is given in Eq. 58). Since $A_{k}$ is proportional to $T/n$ and remembering that we will take the limit where $n\to\infty$, using Baker-Campbell-Hausdorff, we have the exact relation

[TABLE]

The commutator on the right hand side of Eq. 59 has four terms

[TABLE]

As we will show, two of these terms vanish when taking the limit $n\to\infty$.

First, we compute

[TABLE]

with

[TABLE]

Now, terms such as $e^{-\frac{iT}{n}\Phi_{j}}a^{\dagger}e^{\frac{iT}{n}\Phi_{j}}$ can be computed by defining

[TABLE]

and invoking Heisenberg’s equation of motion. Using $[H_{\text{temp}},a^{\dagger}]=-j\phi^{\prime}\tilde{\delta}a^{\dagger}$ and solving the differential equation, we obtain

[TABLE]

Using the result of Eq. 64 into Eq. 62, we have

[TABLE]

Using $\phi^{\prime}\tilde{\delta}=2\phi$ and expanding Eq. 65 to leading order in $\phi$, we obtain

[TABLE]

Inserting Eq. 66 into Eq. 61, we obtain

[TABLE]

A similar calculation shows that

[TABLE]

Next we compute the cross terms. We first have

[TABLE]

with

[TABLE]

Expanding the term proportional to $a^{\dagger}a$ to leading order in $\phi$, Eq. 70 becomes

[TABLE]

Inserting Eq. 71 into Eq. 69, we obtain

[TABLE]

The last commutator to compute is

[TABLE]

with

[TABLE]

so that

[TABLE]

Inserting Eqs. 67, 68, 72 and 75 into Eq. 59, we have

[TABLE]

Since the product $D(A_{n})D(A_{n-1})\cdots D(A_{2})D(A_{1})$ will produce sums in the exponents proportional to $n^{2}$, all terms in Eq. 76 proportional to $(\frac{T}{n})^{3}$ will vanish in the limit where $n\to\infty$. Hence Eq. 76 simplifies to

[TABLE]

where $\overline{\Phi}_{(k-j)}$ is defined as

[TABLE]

We thus have that the product of the modified displacement operators is given by

[TABLE]

We now wish to compute

[TABLE]

First notice that for any operators $A$ and $B$, to leading order $e^{A}e^{B}=e^{B}e^{A}e^{[A,B]}$. Now for terms of the form \text{exp}\Big{[}i\mathcal{E}(t_{k})\mathcal{E}(t_{j})\Big{(}\frac{T}{n}\Big{)}^{2}\sin\Big{(}\frac{T}{n}\overline{\Phi}_{(k-j)}\Big{)}\Big{]}D(\sum_{m}A_{m}), the commutator of the two exponents will be proportional to $(\frac{T}{n})^{3}$ which will vanish in the limit where $n\to\infty$. Hence, we can commute all terms \text{exp}\Big{[}i\mathcal{E}(t_{k})\mathcal{E}(t_{j})\Big{(}\frac{T}{n}\Big{)}^{2}\sin\Big{(}\frac{T}{n}\overline{\Phi}_{(k-j)}\Big{)}\Big{]} to the right hand side of $P_{n}$. Hence we have

[TABLE]

Now, taking the limit where $n\to\infty$ of $P_{n}$, we obtain

[TABLE]

We conclude that

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

Hence to conclude, the unitary evolution of the Hamiltonian described in Eq. 30 is given by

[TABLE]

Recall that in deriving the expression for $V_{+}(0,T)$, in several steps of the calculation we expanded to leading order in $\phi$. Hence the expressions obtained in Eqs. 84, 85 and 86 are only valid to leading order in $\phi$. Hence to leading order in $\phi$, we have

[TABLE]

[TABLE]

and

[TABLE]

B.2 Including the Kerr non-linearity

When including the Kerr non-linearity, the Hamiltonian during the control-displacement gate is now given by

[TABLE]

In this case, we have that

[TABLE]

Hence, we see that the analysis of Section B.1 is exactly the same, with $\phi\to\phi\pm\frac{K}{2}$ in $R_{\pm}(T)$, $A_{\pm}$ and $B_{\pm}$.

B.3 Controlled-displacement gate in the rotating frame

In this section, we consider the same Hamiltonian as in Eq. 30 but with a drive term of the form

[TABLE]

We will go into the rotating frame of both the qubit and the cavity. We define

[TABLE]

With $U(t)=e^{iH_{s}^{\prime}t}$ and applying the transformation $U(t)H(t)U^{-1}(t)-iU(t)\frac{\partial}{\partial t}U^{-1}(t)\equiv H_{R}(t)$ to the Hamiltonian in Eq. 30 with $H_{d}(t)$ given in Eq. 93, we obtain

[TABLE]

The frequency dependence in the drive term of Eq. 95 can be eliminated by choosing $\omega_{d}=\tilde{\omega}_{r}$. Again, we wish to compute the unitary evolution

[TABLE]

Using the Suzuki-Trotter decomposition

[TABLE]

where now we have that $\omega_{\pm}=\pm\chi$.

Comparing Eq. 97 to Eq. 35, we see that we can follow the same steps as in Section B.1 by using the new values for $\omega_{\pm}$ and replacing $\mathcal{E}$ by $\mathcal{E}^{*}$ in the conjugate expressions. Doing so, we obtain

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

B.4 Effects of the non-linear dispersive shift and Kerr term on the unitary evolution of the qubit-cavity system

In this section we will show that regardless of the chosen pulse shape (even with numerical tools such as optimal control), to leading order in $\phi_{\pm}$, the changes to the unitary evolution of the qubit-cavity system due to the non-linear dispersive shift and Kerr terms cannot be completely removed for time scales $T\ll\frac{1}{\phi_{\pm}}$. Using Eqs. 84, 89 and 90, we can write (for instance choosing the terms affecting the $|0\rangle\langle 0|$

[TABLE]

where we define

[TABLE]

and

[TABLE]

Let $A=i\omega_{+}T\phi_{+}(a^{\dagger}a)^{2}$ and $B=I_{1}a^{\dagger}-I_{1}^{*}a-I_{2}\phi_{+}(2a^{\dagger}a-1)a^{\dagger}+I_{2}^{*}\phi_{+}a(2a^{\dagger}a-1)$. Using the Baker-Campbell-Hausdorff lemma and keeping only leading order terms in $\phi_{+}$, we have

[TABLE]

Computing the commutators, we find

[TABLE]

Notice that the last term in Eq. 106 is proportional to $(a^{\dagger}a)^{2}$ and independent of the pulse shape (performing the same calculation as above including the $B_{+}$ term will not change the conclusion). The terms dependent on the pulse shape are expressed as lower powers of $a$ and $a^{\dagger}$. One cannot eliminate all terms proportional to $\phi_{+}$ unless one chooses a time scale comparable to $\frac{1}{\phi_{+}}$. For time scales on the order of $\frac{1}{\phi_{+}}$, higher order terms in $\phi_{+}$ will be relevant and therefore it might possible to choose pulse (with numerical techniques such as optimal control) which can eliminate effects from the non-linear dispersive shift and Kerr terms. However using the parameters of Table 2, $\frac{1}{\phi_{+}}$ is roughly four orders of magnitude longer than the chosen time scale ($T=\frac{\pi}{\chi}$) of our protocol. For such long time scales, effects due to photon loss, damping and dephasing would render the protocol impractical.

Appendix C Probability of measuring accepted measurement strings

In this section, we derive the analytic expression for output states obtained from phase estimation protocols and their corresponding probabilities.

C.1 Phase estimation measurement operator of 1 round

From Fig. 2, the operator describing the evolution of a single round of phase estimation in terms of the measurement operator set \Big{\{}M_{x}=\frac{D(-\sqrt{2\pi})+(-1)^{x}e^{i\gamma}D(\sqrt{2\pi})}{2}|x=0,1\Big{\}} can be written as,

[TABLE]

After the measurement, the state become,

[TABLE]

where $\text{Pr}_{x}$ is the probability of measuring $x\in\{0,1\}$.

C.2 Phase estimation measurement operator after $M$ rounds

After $M$ arounds of the phase estimation circuit, the following state will be generated:

[TABLE]

Where $N=\prod_{i}\text{Pr}_{x_{i}}$ is the normalization factor. From the product of $M$ sums, if we choose $j$ of them to be positive shifts ($D(\sqrt{\frac{\pi}{2}})$), $M-j$ to be negative shifts ($D(-\sqrt{\frac{\pi}{2}})$), we can produce a peak at $\sqrt{\pi/2}(j-(M-j))=\sqrt{2\pi}(j-M/2)$ in $q$ space, and the produced phase is $\prod_{l\in S_{j}}e^{i(\gamma_{l}+x_{l}\pi)}$ ($S_{j}$ is the set of rounds we choose to be positive shifts). There are $\begin{pmatrix}M\\ m\end{pmatrix}$ ways to choose such an $S_{j}$, by choosing any of such $S_{j}$ we can produce a peak at $\sqrt{2\pi}(j-M/2)$. Thus, $|\Phi(x[M])\rangle$ can be written as:

[TABLE]

where $N\approx\sum_{j=0}^{M}|c_{j}(x[M])|^{2}$ with $c_{j}(x[M])=\sum_{\{S_{j}\}}\prod_{l\in S_{j}}e^{i(\gamma_{l}+x_{l}\pi)}$.

C.3 Probability

Assume we already performed $M-1$ rounds of phase estimation and the measurement results $x_{1},...,x_{M-1}$ and angle parameters we chose (denoted as $\gamma_{1},...\gamma_{M-1}$) are known. Then from Eq. 110, the input state to the $M$th round is (including the qubit):

[TABLE]

The probability of measuring $x_{M}$ is thus given by

[TABLE]

Notice that the transfer probability from one peak to another is almost 0 ($<10^{-35}$) assuming an initial squeezing level of 0.2. Therefore, we can make the following approximation

[TABLE]

The above expression is recursive, hence substituting $Pr_{x_{M-1}}$ into it, we find

[TABLE]

Appendix D Details of numerical simulation

From the depolarizing noise model described in Section V.2, we computed error probabilities for all Pauli errors arising from a single fault at the locations indicated by the blue boxes in Fig. 9. For instance, the probability of an $I\otimes Z$ error at location $L_{1}$ is given by

[TABLE]

where

[TABLE]

during the first round, and

[TABLE]

in later rounds. Similarly,

[TABLE]

for the first round and

[TABLE]

With a Pauli error added at either locations $L_{1}$, $L_{2}$, $L_{3}$ or $L_{4}$, we update the state of the qubit before performing the master equation simulation described in Section V.2. With the analytic error probability for the Pauli error, in addition to the evolution of the state during the master equation, we can compute the total probability of obtaining the output state. As mentioned in Section V.2, such an analysis ignored second order Pauli error events (before and after the controlled displacement gate).