Fault-tolerant preparation of approximate GKP states

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.

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 faulttolerant phase estimation protocol maintain good error correcting capabilities even in the presence of noise.

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 [1,2]. 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 [3]. Lastly, given a supply of GKP-encoded Pauli eigenstates, universal fault-tolerant quantum computation can be achieved using only Gaussian operations [4].
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 [5][6][7][8][9][10]. However to our knowledge, no clear definitions for fault-tolerantly preparing GKP states 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 [11], 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 A, we briefly review the phase estimation protocol described in [7]. In Section IV B, we obtain a new phase estimation protocol and prove that it is fault-tolerant under our definitions. In Section IV C, 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 A, we provide analytic expressions for the time evolution of the qubitcavity coupling when implementing a controlled-displacement gate. The expressions are derived in Appendix B. In Section V B, 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 [7,12], 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 √ πp and S q = e 2i √ πq , can correct shift errors of size at most √ π 2 . Note that for the displacement operator D(α) = e αa † −α * a , we can write S p = D( √ 2π) and S q = D(i √ 2π). 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 |0 and |1 states are given by 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 state can be written as a squeezed vacuum state |sq with squeezing parameter ∆ which is given by One can then apply a sum of displacements with a Gaussian filter to obtain where C is a normalization coefficient. As was shown in [7], it is natural to have∆ = ∆. 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 [7] 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 √ π 2 with certainty. For example, suppose we have an approximate |0 GKP state with a peak at q = 0 subject to a shift error e −ivp with |v| ≤ √ π 2 . The finite width of the Gaussian peaks will have a non-zero overlap in the region Thus with non-zero probability the state can be decoded to |1 instead of |0 (see Fig. 2 for an illustration). . Peaks centered at even integer multiples of √ π in q space. The peak on the left contains large tails that extend into the region where a shift error is decoded to the logical |1 state. The peak on the right is much narrower. Consequently for some interval δ, the peak on the right will correct shift errors of size 1 . In what follows, we will assume that after preparing the desired ancilla states using some state preparation protocol, the ancillas will be used in a Steane error correction scheme (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 Steane error correction, it is important to prepare approximate |0 and |+ 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 √ π 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 ( √ π 2 − δ, ) p,q correctable for some desired δ with large (so that the probability of correcting a shift √ π 2 − δ is small). Note that the notion of correctable in Definition 3 assumes that only the data and ancilla qubits used in Steane error correction have shift errors. If other operations such as the CNOT gates introduce additional errors, the correctable threshold would be smaller than √ π/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. 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δ occurs on the oscillator Hilbert space. We will say that the protocol is an (m,δ)-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. 1 Using the optimizations considered in [13], using Knill error correction could potentially increase the threshold of √ π 6 to a larger value.
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 [11]). 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 [7]). By performing Steane error correction, 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δ in Definition 5 relates to the largest allowed size of the shift error which occurs during the protocol. This should not be confused with δ in Definitions 1 and 2 which relates to the size of a shift error that can be corrected by an approximate GKP state 2 .
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δ, 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δ on the oscillator Hilbert space, the remaining shift error on the output state will be corrected in a perfect Steane error correction protocol. As we will see in Section IV, due to measurement errors, the phase estimation protocol presented in [7] needs to be modified in order to be a (1,δ) 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 Steane error correction. Therefore we will compare the ( √ π 2 − δ, ) 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 max is computed as Similarly, in p space we have For protocols preparing |+ , 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 A we will give a brief review of the phase estimation protocol presented in [7]. In Section IV B, we will show how the protocol can be modified so that it becomes a (1,δ)-Fault-tolerant state preparation of an approximate GKP state and we will provide the value ofδ.
A. Brief review of the phase estimation protocol for preparing approximate GKP states The phase estimation circuit for preparing an approximate |0 GKP state is given in Fig. 2. The H gate is the Hadamard gate and Λ(e iφ ) = diag(1, e iφ ). 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 eigenvalue 3 e iθ . Additionally, an estimated value for the phase θ is obtained. After computing the phase, the state can be shifted back to an approximate +1 eigenstate of S p . 2 The probability 1 − of correcting a shift of size √ π 2 − δ depends on the location and width of the peaks. 3 The phase θ that is obtained depends on the measurement outcome of the ancilla qubit in each round. Figure 2. Phase estimation circuit for preparing an approximate |0 GKP state. The initial state of the cavity is taken to be the squeezed vacuum state defined in Eq. (2). The circuit effectively projects the input squeezed vacuum state onto an approximate eigenstate of the Sp operator while at the same time estimating the phase of the eigenvalue. The phase can be computed so that the output state can be shifted back to an approximate +1 eigenstate of Sp.
There are a variety of ways to choose the phases φ at each round in the gate Λ(e iφ ). In a non-adaptive phase estimation protocol, half of the values for φ can be chosen to be 0 and the other half can be chosen to be π/2. In an adaptive phase estimation protocol, if the phase to be estimated is θ and assuming no prior knowledge of θ, during the first round the phase can be chosen as φ 1 = 0. For later rounds, it was shown in [7] that the next phases can be chosen as In Eq. (10), the probability P φ (x[m]|θ) is the probability of obtaining measurement outcomes x 1 , · · · , x m when the produced output eigenstate |ψ θ has eigenvalue e iθ . Since the measurement results for each round are independent, the estimated probability is given by By analytically computing the evolution of the state, the exact probability of obtaining the outcomes x[M ] is derived in Appendix C. Given the final measurement record x[M ], the estimated phaseθ is chosen as The shift back correction based on the estimated phaseθ is given by e iθ 2 √ π q . Note that for either the adaptive or non-adaptive protocol, the estimated phase and probabilities are computed using Eqs. (11) and (12). 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 √ π > √ π/6, then the shift error will be uncorrectable and thus the phase estimation protocol will not be fault-tolerant 4 . 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 , · · · , 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 out [m] will have one bit which differs from the bit string x corr [m] that would have been obtained in the fault free case (so 4 Note that |v − w| should be taken modulo √ π since a perfect |0 GKP state in p space has peaks at any integer multiple of √ π. Calculation of acceptance set A m andδ. Consider all output states obtained from an m round fault-free phase estimation protocol of Section IV A. Let A m = ∅ (which we call the acceptance set) and Γ m = ∅. For each output state x[m], do the following (1, δ) fault-tolerant m round phase estimation protocol for measurement readout errors. Consider the acceptance set A m = ∅ andδ 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] / ∈ 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 √ π/6. Our protocol only guarantees protection against a single readout error.  Table I. The first row displays the accepted measurement strings (set A4) for the four round non-adaptive phase estimation protocol. If the protocol does not output an element of A4, 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. The third row gives the probabilities of obtaining each state at the output of the phase estimation protocol.
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. This indicates that the adaptive phase estimation protocol described in Section IV A is not fault-tolerant to single measurement readout errors 5 . However, if the phases are computed using the non-adaptive phase estimation protocol and applying the protocol described above, we find thatδ = √ π 6 − 0.235 ≈ 0.0604. The set A 4 of accepted measurement strings is given in Table I and has |A 4 | = 5. This result indicates that although the adaptive phase estimation protocol outperforms the non-adaptive protocol in the fault-free case ( [7,18]), the adaptive phase estimation protocol is cannot be used in the presence of measurement readout errors. Therefore in a noisy implementation of phase estimation, the non-adaptive protocol should be used. Lastly, for the four round noise free non-adaptive phase estimation protocol, the probability of obtaining an accepted √ 2π) gate resulting in a D( √ 2π) error, the flag qubit will be measured as 1 instead of 0. If the flag is measured as 1, we abort the protocol and start anew.
state was numerically computed to being approximately 0.4828. These values are in close agreement to the probabilities obtained using the analytical expressions found in Appendix C.
Suppose now that for m rounds the set A m is not empty. We then know there exists aδ such that the protocol is a (1,δ) 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( √ 2π) gate will do nothing since the qubit is in the |+ 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( √ 2π) gate can result in incorrectly applying the D( √ 2π) displacement on the cavity. If several damping events occur during the protocol, the output state in p space could potentially be badly deformed. In [7] 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( √ 2π 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 [19] for preparing cat codes. However large cavity displacements were mitigated by interacting the cavity with a three level system |f , |e and |g instead of a qubit. A transition from |f to |e 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( √ 2π) gate which uses a single extra flag qubit [20][21][22][23][24][25]. Consider the circuit in Fig. 3. If a bit flip error occurs before or during the application of the controlled-D( √ 2π) 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 |+ ancilla qubit causes the wrong D( √ 2π) 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 Λ(e iφ ) gate will change the sign of φ. 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 Λ(e iφ ) 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δ in p space (i.e. a shift of the form e −iδ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 < √ π/6. Following Definition 5, we conclude that the protocol described in this section is a (1,δ)-fault-tolerant state preparation of an approximate logical |0 GKP state withδ = √ π/6 − 0.235 ≈ 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.
C. Goodness of approximate |0 states obtained from the noise free phase estimation protocol In Section V B, 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 nonadaptive 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 δ, we computed 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 √ π 2 − δ. Plots for the states 0101 and 1000 which belong to A 4 (see Table I) are given in Fig. 4. The ( √ π 2 − δ, ) 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 √ π 2 − δ, both states can correct the shift error with probability close to 1.

V. CIRCUIT QED IMPLEMENTATION
A. State evolution in the dispersive regime In this section we will describe a direct implementation of the controlled-D( √ 2π) gate. From [26,27], the Hamiltonian describing the coupling between a qubit and a cavity in the dispersive regime is given by Here g is the coupling strength between the qubit and the cavity and ∆ = |ω r − ω a |. 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 φ [27]. Note that a more systematic treatment of the qubit as an anharmonic oscillator leads to an additional term in Eq. (13) given by − K 2 (a † a) 2 which is referred to as the Kerr non-linearity [35]. Hence, in our analysis, we choose the system Hamiltonian to be given by The direct implementation of the controlled-D( √ 2π) gate can be achieved using the drive Hamiltonian where E(t) describes the pulse shape of the drive and ω d is the drive frequency. Thus the total Hamiltonian describing the evolution of the qubit-cavity system during the implementation of the controlled-D( √ 2π) gate is given by Going into the rotating frame of the qubit and the cavity and defining φ ± = φ ± K and ω ± = ±χ, in Appendix B we show that where In deriving Eqs. (18) to (20), we kept only leading order terms in φ ± since the Hamiltonian in Eq. (14) neglects higher order terms in φ and K. We keep only leading order terms since with current experimental parameter values, φ ± is roughly three orders of magnitude smaller than χ. Firstly, notice that the terms R ± (T ) introduce a relative rotation between the |0 and |1 state of the qubit. Neglecting the terms proportional to φ ± , it was shown in [7] that by choosing the total interaction time to be T = π/χ, the relative rotation between |0 and |1 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 ± (T ) does not depend on the pulse shape 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 π/χ to ensure that we eliminate the relative rotation from the leading order terms in Eq. (18). In Appendix 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 1 φ± which (for the parameter values in Table II) is four orders of magnitude larger than π χ . 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 ± 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 φ ± 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 where µ = π 2χ and σ = π 10χ . With this Gaussian pulse we obtain The amplitude of the Gaussian, given by −2.09562 √ 2π T , 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 − π 2 in Fig. 3 is not required. Furthermore, the first term in Eq. (22) is independent of χ. 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 χ will minimize the undesired effects arising from the non-linear dispersive shift and the Kerr term while still producing the desired gate.
The B ± terms of Eq. (20) result in a phase difference between the |0 and |1 qubit states. However, computing the first integrals in Eq. (20) (i.e. terms before φ ± ), 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 ± . Note however that the φ ± terms in B ± will not be cancelled and will thus introduce additional shift errors.

B. 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 B. The simulation is performed in several steps that we now describe. First, the controlled displacement gate is modeled using the following master equatioṅ where H(t) is given by Eq. (16) after going into the rotating frame and D[L]ρ = (2LρL † − L † Lρ − ρL † L)/2. The density matrix corresponds to the joint state of the cavity, ancilla qubit and flag qubit. The parameters κ, γ 1 and γ φ 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 π χ with χ = g 2 ∆ − φ. Both before and after the controlled displacement gate, the cavity is subject to photon loss. 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µs. Similarly, the evolution time after the controlled-displacement gate to the time the ancilla is measured is also chosen to be 0.14µs. Thus during  Table II. Parameter values for the two simulations of the non-adaptive noisy phase estimation protocol. The values for κ are chosen based on state of the art 3D cavities. The parameters in the second column results in a T1 time given by T1 = 2π κ = 10ms. For a resonator frequency fr = 7GHz, the quality factor is given by Q = fr κ/(2π) = 7 × 10 7 . Figure 5. Probability of correcting a shift errors of size √ π 2 − δ for parameters chosen from the second (labeled "Simulation 1") and third (labeled "Simulation 2") column of Table II in addition to the case where no noise is present. The plots are for the states 1000, 0101 and 1111 obtained from the four round fault-tolerant phase estimation protocol described in Section IV.
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 Table III. The second row corresponds to the probabilities of obtaining the output states of A4 for the parameters chosen from the second column of Table II conditioned on all flag measurements being 0. The third row is identical but for the parameters chosen from the third column of Table II. 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. one round, the cavity freely evolves with photon loss for a total of 0.28µs. In addition, we allow all qubit locations to fail with the following depolarizing noise model 1. With probability p, each two-qubit gate is followed by a two-qubit Pauli error drawn uniformly and independently from {I, X, Y, Z} ⊗2 \ {I ⊗ I}.
3. With probability 2p 3 , any single qubit measurement has its outcome flipped. 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 [28], 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 immediately after the first CNOT gate, before both measurement locations and before the phase gate. Using these probabilities, all single fault locations were considered. At each location, all possible Pauli operators based on their associated probabilities were added. Therefore, instances with two or more faults occurring on the qubit Hilbert space before and after the controlled-displacement gate were negelcted. However our analysis is still more complete than previous implementations which considered only measurement errors and errors during the controlled-dispalcement gate.
We performed two different simulations where for each simulation, the chosen parameter values are given in Table II. The parameters chosen for the first simulation (middle column in Table II) are based on current experimental values for 2D and 3D cavities [19,[29][30][31][32][33]. 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 √ π 2 − δ 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 − δ plots. Note that noise during the phase estimation protocol introduces more shift errors in p space. Therefore in Fig. 5, only − δ plots in p space are shown.
For the four round phase estimation protocol, we find that the state 1010 has a similar − δ plot to the one for the state 0101 whereas the state 0010 has a similar − δ plot to the one for the state 1000. The probability of obtaining each state for the two noisy simulations (parameters in Table II) are given in Table III. 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,δ)-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δ 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 A 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 [34,35] 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 [8,10]. 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 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. 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 |+ ancilla state before the controlled displacement gate as shown in Fig. 8. Let |ϕ be the input state part of the cavity Hilbert space. Defining p as the damping rate, after the controlleddisplacement gate, the state of the system can be written as 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 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) Applying the remaining gates in Fig. 8, the final state prior to the measurement of the ancilla qubit is If the measurement of the ancilla is 0, the output state will be If the measurement of the ancilla is 1, the output state will be 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( √ 2π) gate.
1. Implementation of the controlled-D( √ 2π) 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: where and The term φ(a † a) 2 Z corresponds to the non-linear dispersive shift. Note that we have not included the Kerr term − K 2 (a † 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. (B3), we neglected a term of the form λX (where λ = (g/∆) and all higher powers in λ. 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. (B3), we represent the drive pulse E(t) as where ω d =ω r − χ 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. (B1) is given by The right-hand side of Eq. (B5) can be computed using the Suzuki-Trotter decomposition, so that wheret ≡ T /n, ω ± ≡ω r ± χ. Now, the two terms on the right-hand side of Eq. (B94) can be decomposed as where and Hence from the above, we see that we need to compute terms of the form In order to compute the products appearing in Eq. (B10), we will use the Heisenberg equation of motion as a trick to simplify commutation rules. Defining, and with we have that where ∂At ∂t H = e iH t ∂At ∂t e −iH t .
Computing the partial derivative, we have that The terms inside the bracket of Eq. (B16) can be computed by invoking Heisenberg's equation of motion with the creation and annihilation operators. We obtain with where we defined We thus have the following set of coupled differential equations: and The solution to Eqs. (B20) and (B21) is given by and To see this, we take a derivative of Eq. (B22) to obtain Comparing Eqs. (B20) and (B24), we must have that Using Eqs. (B22) and (B23) and defining H te ≡ −kt(φ + 2φ(a † a − 1)), we have that as desired. A similar calculation can be done for a † (t). Hence we have that with Defining H ≡ sE(t k )(a + a † ) and using [H , a † ] = sE(t k ), we obtain e iH t a † e −iH t = isE(t k )t + a † so that Eq. (B31) becomes [a † a, At] = −itE(t k )At(a † − a + itE(t k )).
Next we compute [(a † a) 2 , At]. Again using Eq. (B29), we have that with Using the relations e −sB a † e sB = istE(t k ) + a † , e −sB ae sB = −istE(t k ) + a and Eq. (B34), we obtain Computing the integral and inserting the result in Eq. (B33) we obtain To obtain the final version of Heisenberg's equation of motion, we must now compute e iH t [a † a, At]e −iH t and e iH t [(a † a) 2 , At]e −iH t . For ease of notation, we define We then have Similarly, we have Inserting Eqs. (B27), (B28), (B38) and (B39) into Eq. (B14), we obtain Note that for arbitrary functions f (t) and g(t), a differential equation dx . Using this and the fact that A(0) = I, the solution to Eq. (B40) is given by where Before computing the integral in Eq. (B42), it is important to point out that in the dispersive regime, we are only interested in terms that are leading order in φ. The integrals will require us to compute the inverse of Φ k . Before doing so, let us defineδ and write Then it is straightforward to see that Hence we can truncate to leading order inδ as long as (n − 1) 2δ2 1. In general then, we will have Again keeping only leading order terms in φ, we can approximate all the integrals in Eq. (B42) as when m ≥ 1. We separate terms in M (t) into two parts. The first part consists of terms which are proportional to ω + e ±iΦ k x and e ±iΦ k x which we define as I 1 and the second part are terms proportional to φe ±iΦ k x which we define as I 2 . Hence we have Note that in Eq. (B48) we have kept terms proportional to Φ −2 k that are multiplied by ω + . We can approximate the term (kω which is kω + Φ −1 k and can thus be neglected. Lastly, we also have a term of the form kω + Φ −1 k which is given by Inserting Eq. (B50) into Eq. (B48), we obtain The remaining terms in M (t) are then grouped into the operator I 2 Performing the same approximations as in the calculation of I 1 and defining we obtain Combining Eqs. (B51) and (B54) into Eq. (B42), the final expression for M (t) is given by We know that A(t) = exp M (t) is given by Eq. (B13) (so that A(t) = R −k n D t k R k n ). Since the unitary operator describing the evolution during the controlled-displacement gate is obtained by taking products of A(t) (for different values of k) as in Eq. (B7) and taking the limit where n goes to infinity, only terms proportional to ( T n ) or k( T n ) 2 will be non-zero. The reason is that although many of terms in the exponents of R −k n D t k R k n and R −j n D tj R j n do not commute, summing exponents (arising from the products of these terms) can give a factor of at most O(n 2 ). Now since the terms in the exponents don't commute, from the Baker-Campbell-Hausdorff lemma, we cannot directly sum the terms in the exponents. However products between terms consisting of higher powers of T /n in Eq. (B55) cannot sum to cancel the highest power of n in the denominator and will thus vanish when taking the limit n goes to infinity.
Writing V (0, T ) (Eq. (B94)) as using Eq. (B55) (in addition to the paragraph below it) and the fact that R n n = e −iT (ω+−φa † a)a † a ≡ R + (T ), we can write By combining the terms in Eq. (B57), a great simplification takes place. To see this, notice that where we used Eq. (B43) when going from the third to last line. Lastly notice that φ . Thus keeping terms to leading order in φ (we will also expand terms such as e ±iΦ k T /n to leading order in φ at the end of the calculation), we obtain where we defined 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. (B59), 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. (B60)). Since A k is proportional to T /n and remembering that we will take the limit where n → ∞, using Baker-Campbell-Hausdorff, we have the exact relation The commutator on the right hand side of Eq. (B61) has four terms As we will show, two of these terms vanish when taking the limit n → ∞.
First, we compute with Now, terms such as e − iT n Φj a † e iT n Φj can be computed by defining Using φ δ = 2φ and expanding Eq. (B67) to leading order in φ, we obtain Inserting Eq. (B68) into Eq. (B63), we obtain A similar calculation shows that Next we compute the cross terms. We first have with Expanding the term proportional to a † a to leading order in φ, Eq. (B72) becomes Inserting Eq. (B73) into Eq. (B71), we obtain The last commutator to compute is .

(B78)
Since the product D(A n )D(A n−1 ) · · · D(A 2 )D(A 1 ) will produce sums in the exponents proportional to n 2 , all terms in Eq. (B78) proportional to ( T n ) 3 will vanish in the limit where n → ∞. Hence Eq. (B78) simplifies to where Φ (k−j) is defined as We thus have that the product of the modified displacement operators is given by We now wish to compute sin T n Φ (k−j) D( m A m ), the commutator of the two exponents will be proportional to ( T n ) 3 which will vanish in the limit where n → ∞. Hence, we can commute all terms exp iE(t k )E(t j ) T n 2 sin T n Φ (k−j) to the right hand side of P n . Hence we have Now, taking the limit where n → ∞ of P n , we obtain where R + (T ) = e −iT (ω+−φa † a)a † a , (B86) and Recall that in deriving the expression for V + (0, T ), in several steps of the calculation we expanded to leading order in φ. Hence the expressions obtained in Eqs. (B86) to (B88) are only valid to leading order in φ. Hence to leading order in φ, we have and (B92)

Including the Kerr non-linearity
When including the Kerr non-linearity, the Hamiltonian during the control-displacement gate is now given by H(t) =ω r a † a +ω a Z + χa † aZ − φ(a † a) 2 Z − K 2 (a † a) 2 .
In this case, we have that Hence, we see that the analysis of Appendix B 1 is exactly the same, with φ → φ ± K 2 in R ± (T ), A ± and B ± .

Controlled-displacement gate in the rotating frame
In this section, we consider the same Hamiltonian as in Eq. (B1) but with a drive term of the form H d (t) = E(t)a † e −iω d t + E * (t)ae iω d t .
We will go into the rotating frame of both the qubit and the cavity. We define H s =ω r a † a +ω a Z.
With U (t) = e iH s t and applying the transformation U (t)H(t)U −1 (t) − iU (t) ∂ ∂t U −1 (t) ≡ H R (t) to the Hamiltonian in Eq. (B1) with H d (t) given in Eq. (B95), we obtain H R (t) = χZa † a − φ(a † a) 2 Z + E(t)a † e i(ωr−ω d )t + E * (t)ae −i(ωr−ω d )t (B97) The frequency dependence in the drive term of Eq. (B97) can be eliminated by choosing ω d =ω r . Again, we wish to compute the unitary evolution Using the Suzuki-Trotter decomposition where now we have that ω ± = ±χ. Comparing Eq. (B99) to Eq. (B94), we see that we can follow the same steps as in Appendix B 1 by using the new values for ω ± and replacing E by E * in the conjugate expressions. Doing so, we obtain where R ± (T ) = e −iT ω±a † a (1 ± iT φ(a † a) 2 ), (B101) 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 φ ± , 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 1 φ± . Using Eqs. (B86), (B91) and (B92), we can write (for instance choosing the terms affecting the |0 0| R + (T )D(A + ) = e −iω+T a † a e iω+T φ+(a † a) 2 e I1a † −I *