A memory-assisted decoder for approximate Gottesman-Kitaev-Preskill codes

We propose a quantum error correction protocol for continuous-variable finite-energy, approximate Gottesman-Kitaev-Preskill (GKP) states undergoing small Gaussian random displacement errors, based on the scheme of Glancy and Knill [Phys. Rev. A 73, 012325 (2006)]. We show that information gained from multiple error-syndrome extractions can be used to enhance the protection of GKP encoded qubits from amplitude damping. Furthermore, we show that the expected total displacement error incurred in multiple rounds of error followed by syndrome extraction is bounded by $2\sqrt{\pi}$. By recompiling the syndrome-extraction circuits, we show that all squeezing operations can be subsumed into auxiliary state preparation, reducing them to beamsplitter transformations and quadrature measurements.

Introduction -Encoding and manipulating quantum information in continuous variable (CV) systems [1][2][3] is a promising route to realising a useful quantum computing device. Large scale CV cluster states can be generated on demand [4], and fast, high-quality one-and two-qubit clifford gates are deterministically available [5]. Fault tolerant, measurement-based quantum computation is possible using CV cluster states along with CV measurements and non-Gaussian state injection [6,7], though the levels of squeezing required for fault tolerance are beyond the reach of current experiments [6,[8][9][10][11].
A leading CV quantum computing approach is one proposed by Gottesman, Kitaev and Preskill [12], which is based on the idea of encoding a qubit within an (infinite dimensional) oscillator. Ideal codeword wavefunctions within this paradigm correspond to infinite-energy Dirac combs, and are commonly referred to as GKP states. In practice, this ideal wavefunction is replaced by a finiteenergy approximation, such as a comb of narrow Gaussian peaks modulated by a broad Gaussian envelope.
The appeal of GKP-encoded qubits is that they possess both an intrinsic robustness to physically motivated error channels and natural schemes for error syndrome extraction and correction. The originally proposed GKP error correction protocol was based on CV stabilizer generalisations of the Steane circuits [13] for error syndrome extraction. Subsequently, Glancy and Knill proposed a different error correction method based on a beamsplitter transformation [14], which will provide the basis for the scheme detailed in this paper.
We define an approximate GKP codewords with width ∆ = (∆, κ) as where G Σ (z) = exp (− z 2 2Σ 2 ) and µ ∈ {0, 1} defines the logical basis states. Informally, for the logical 0 (1) state we have a superposition of Gaussians of width ∆ centred at even (odd) values of √ π, with an overall Gaussian envelope of width 1/κ. Better approximations to the ideal GKP state are achieved with smaller values of ∆ and κ, although these also correspond to larger average energy and an apparent increase in experimental difficulty.
Recently, Albert et al. [15] showed that the GKP code outperforms a number of other bosonic codes when states are exposed to Gaussian random displacement errors. We write such an error acting on the stateρ as (2) where the operatorD(α) shifts the state in phase space by Re{α}, Im{α} in the q-and p-quadratures respectively, and the width σ 0 quantifies the extent of the error. Notably, this error model describes amplitude damping that is preceded by an offsetting pre-amplification, and it is therefore highly relevant to many experimental platforms. To date, however, explicit error correction protocols for currently accessible, approximate states are lacking.
In this letter we present a new decoder for the GKP code undergoing Gaussian random displacement errors, based on the Glancy and Knill error recovery scheme and Bayesian estimation. In particular, we extend the scheme to enable enhanced error estimation using multiple syndrome extractions, enabling improved error supression which we show to be useful in extending the lifetime of states with a mean number of bosons as little as ten. We find that many rounds of syndrome extraction without active corrective displacement causes the qubit to drift in phase space by at most ≈ 2 √ π on average in each quadrature. Additionally, we recompile the syndrome extraction circuit in an experimentally friendly way, such that squeezing need only be applied to auxiliary states which can be prepared offline. Syndrome extraction -The GKP syndrome extraction scheme of Glancy and Knill [14] can be broken FIG. 1. The q-quadrature error-syndrome extraction circuit (q-SE), as proposed by Glancy and Knill. The GKP qubit and an auxiliary state are input in the top and bottom modes, respectively. Following beamsplitting and squeezing operations, the error syndrome xm is generated by a q-quadrature measurement on the auxiliary mode.
down into two sequential circuits: q-SE and p-SE, which extract the error syndromes in the q-and p-quadrature respectively (see Fig. 1).
We begin our analysis with an arbitrary input qubit wavefunction Q ∆ (x) that has undergone an unknown displacement error (u, v) This corrupted qubit is input into the top mode of the q-SE circuit, while an auxiliary GKP state |ψ ∆ + is input into the bottom mode. The action of the q-SE circuit can be understood in terms of the two-mode q-quadrature wavefunction. The beamsplitterB π 2 causes an anticlockwise rotation by 45 • , and the subsequent squeezerŜ √ 2 scales the top-mode quadrature by √ 2. The error syndrome is then generated by a q-quadrature measurement of the auxiliary mode.
The p-SE circuit procedes similarly, although there are subtle differences beyond a change of variables q → p, and squeezing in the conjugate direction (Ŝ √ 2 →Ŝ † √ 2 ). Indeed, we note that we must change the auxiliary state 2)-a departure from the Glancy and Knill proposal which we found necessary for sequential q-SE and p-SE circuits to be applied to approximate GKP states (see Supplementary Material).
Following sequential q-SE and p-SE circuits, the qubit transforms to with a high probability when the displacement is small. We bound the probability of this approximation in the next section. Under this approximation, the qubit at the output of the syndrome extraction is a displaced version of the input qubit, where the displacement is a function of the measurement outcomes and the initial error displacements. The quantities θ(x m , u) and θ(p m , v) in equation (4) can interpreted as the total displacement experienced by the qubit after errors and syndrome extraction in the q-and p-quadrature respectively 1 . While an error correcting procedure could involve estimation of the displacement θ(x m , u) and immediately applying a corrective displacement based on this estimate, we note that allowing an uncorrected qubit to undergo a further syndrome extraction process only results in another displaced version of the input qubit. With this in mind, our decoder uses measurement information from multiple rounds of errors and syndrome extraction without intermediate correction to estimate a single corrective displacement to apply. The memoryless and memoryassisted approaches are summarized in Fig. 2.
We now find expressions for multiple rounds of errors and syndrome extractions. Writing the error shifts and measurement results at the start of the h th syndrome extraction round as (u h , v h ) and (x (h) m , p (h) m ), respectively, we define the recursive quantities: The q-quadrature displacement after h syndrome extraction rounds is then given by and the expression for p-quadrature displacement is the same after making the substitutions: The derivation of equation (6) relies on the good approximation that f * step is independent of u k for all values of k (valid for u < √ π/4). If we stop the iterated syndrome extraction at round h, then the failure probability of this approximation is bounded by 1 θ(xm, u) = u 2 − f * step (xm).
a) The standard q-and p-quadrature syndrome extraction circuits (left), which feature squeezing operations on the qubit mode, can be represented by a measurement dependent Kraus operatorK (centre), and recompiled such that squeezing is applied to the p-auxiliary state, with a reinterpreted p-quadrature measurement (right). b) Memoryless error correctionsyndrome measurement results are used to inform a corrective displacement at each step, and no information is carried forward to future correction steps. c) Memory-assisted error correction-no active corrective shift is performed after each syndrome measurement. Instead, all the syndrome measurement results are used together to decode and perform a single corrective displacement after M steps.
Without active corrective shifts at each round, the distance a GKP qubit drifts in phase space after M rounds is bounded by θ(x is a Gaussian random variable with mean 0 and variance σ 2 (see Supplementary Material). It follows that the expected value for total displacement error after any number of syndrome extraction rounds is bounded above by 2 √ π.
In Fig. 2a, we identify a practical simplification to the combined q-SE, p-SE circuit that moves all squeezing operations offline onto auxiliary state preparation. To do so, we use a modified p-SE auxiliary state,Ŝ † √ 2 |ψ ∆ 0 , and reinterpret the p-quadrature measurement (p m → √ 2p m ) (see Supplementary Material).
Bayesian estimation decoder -For a well characterised loss channel (i.e. known value of σ 0 ), we now show how to use the q-and p-SE measurement outcomes to estimate the unwanted displacements.
In the case of a single round of q-SE and p-SE, the probability of measuring x m given a shift error u can be approximated as P(x m |u) ∝ ψ 3 . The corresponding p-quadrature conditional probabilities can be obtained by substituting (u, x m , ∆, κ) → (v, p m , 2∆, κ/2), where the asymmetry between the qand p-quadrature originates from the differing widths of the auxiliary states. For convenience, we set κ = ∆ from this point onward.
For multiple rounds, the probability in the h th round of obtaining the q-quadrature measurement outcome x Since the error displacements are normally distributed random variables with standard deviation σ 0 , we can use Gaussian priors, G σ0 (u h ), along with Bayes' theorem for M total rounds to obtain: The posterior, P M ( u| x m ), is the distribution of the error shifts conditional on a certain set of measurement results, which we note is independent of all p-quadrature measurement results p h m . Under the assumption that σ 0 and ∆ are small compared to √ π, P to order O We note that each of theũ k is straightforwardly computed from the quadrature measurement results, and that the estimate of the total displacement correction is given by The inverse of the covariance matrix (i.e. the precision matrix) for this estimation can be derived explicitly as: with m k,l = max{k, l}.
The q-quadrature displacement θ(x (h) m , U h ) in equation 6 has a correctable part that depends solely on the measurements results and another part that is a random variable with variance V q = Var( M k=1ũ k 2 M −k+1 ) which can be computed using a Neumann series expansion of the inverse covariance matrix. The series converges if σ0 ∆ < 1 2 (which works for all M > 1), giving which implies that V q → σ 2 0 3 very quickly as M grows if σ0 ∆

1.
Numerical results -The aim of a quantum error correction procedure is to protect logical quantum information. Pantaleoni et al. recently reported a method to extract logical information from approximate GKP states [16]. In essence, this reduces a density matrix function describing a CV state to a qubit density matrix (ρ CV →ρ qubit ). Using this, fidelities between qubit density matrices at the input and output of an error correction procedure can be computed in order to benchmark its performance.
In Fig. 3 we present the results of numerical simulations for the specific case of a GKP code with ∆ ≈ 0.22 (corresponding to an average number of bosonsn ≈ 10). For varying numbers of rounds of error channel application followed by syndrome extraction, we benchmark error correction with our memory-assisted decoder against density matrix function simulations of the memoryless decoder and no error correction at all. We observe that error correction is improved by using our memory-assisted decoder instead of the memoryless decoder-a finding similar to that of Vuillot et al. [17] and Noh et al. [9] for Steane-based GKP syndrome extraction. We also see that the quality of the memoryassisted error corrected state is higher than the uncorrected state up to a total error variance of σ 2 tot ≈ 0.06, despite the fact that state contains a low average number of bosons. Furthermore, the fidelity of the input state is approximately preserved through application of the memory-assisted error correction scheme until a critical point (σ 2 tot ≈ 0.06 in Fig. 3). Beyond this point, the failure probability (equation (7)) for our approximation of the total displacement error becomes non-negligible.
Summary and outlook -We have described an explicit protocol for GKP QEC that provides improved protection to Gaussian displacement errors with reduced experimental requirements. Notably, we have shown that GKP states undergo a total displacement (and therefore require a correction) that is approximately bounded by 2 √ π in each quadrature, after multiple rounds of error syndrome extraction. A memory-assisted decoder based on Bayesian estimation was developed to specify the near-optimal corrective displacement, and numerical simulations have shown that it results in significantly enhanced error correction compared to a memoryless decoder when applied to an approximate GKP state with low average number of bosons.
Recent experiments by Campagne-Ibarcq et al. [11] have demonstrated GKP states of a superconducting microwave resonator with (∆, κ) = (0.16, 0.32). This advance suggests that our results might be exploited by near-term experiments to improve implementations of CV error correction. Additionally, the Glancy-Knill syndrome extraction circuit at the core of our method may be better suited to realistic devices than the more studied SUM gate syndrome extraction because the SUM gate requires more gates to implement [18].
The error model that we have considered in this work also describes errors that arise in an injected GKP state 'hopping' on a CV cluster state with finite squeezing [19]. In this situation, the finite squeezing of the cluster state induces a displacement error with width σ 0 . We anticipate that it will therefore be possible to apply our error correction scheme between nodes on a CV cluster state graph, and thereby lower the squeezing threshold for universal fault tolerant quantum computation using CV cluster states and GKP state injection.

DATA AVAILABILITY STATEMENT
This is a theoretical paper and there is no experimental data available beyond the numerical simulation data described in the paper. Kwok Ho Wan performed all numerical experiments, all the authors contributed to the manuscript and the project was supervised by Alex Neville and Steve Kolthammer.
We begin with a graphical illustration with a sydrome extraction step. We define an error shift in the q and p quadrature as (u, v). For example, in the q-SE stage, take an error free (u, v = 0) input state of ψ . The peaks will appear at every integer value of π 2 in the y-axis and every integer value of √ π 2 in the x-axis in figure 2c.
The measurement conditioned displacement in the q-SE can be pictured by taking a line of constant y-values (y = x m ) out from plot in figure 2c and then asking, how much shift is required such that the peak locations are fixed at 2n √ π for logical 0 and (2n + 1) √ π for logical 1, n ∈ Z. Glancy and Knill (GK) [1] suggested an adaptive shift function of the form (for errors |u|, |v| < √ π/2): where rem[z , 4] is the remainder of z when divided by 4, round{z} is the rounding function that includes positive and negative values in its range and domain and SAW is the function defined in FIG. 5. The SAW part of f GK is to deal with errors. Since we have no external errors we can ignore the sawtooth linear part of f GK . We can define f step as: We can look at how the peaks shift in figure 2d. the effect of shift function f step (y). We can see that all the red and blue peaks return to even and odd integers of √ π effectively.
Proof. The input state into the q-quadrature syndrome extraction circuit has the wavefunction: Φ(x, y) ∝ 1 µ=0 n,m∈Z 2 Prior to the conditional shift, the transformation of the beam-splitter, squeezing and measurement outcome at x m results in the change of variables in the wavefunction from: (x, y) → (x + xm √ 2 , −x + xm √ 2 ) . The resulting wavefunction will now be (up to global nor-malisation constant): Regrouping the Gaussian functions with the same width and then competing the square up in the exponent, i.e. complete the square of the exponents of G 1 √ π ] separately, this implies, The natural indices in this sum are 2n + m and 2n − m.
We want to decouple the n and m indices so that the double sum factors out into a product of single sums. If we make the following re-indexing, 2n + m = 4j + c and 2n−m = 4l−c, where j, l ∈ Z 2 and c ∈ {−1, 0, 1, 2}, then we can trade the infinte (n, m) ∈ Z 2 coupled double sum for a decoupled triple sum over c, j, l ∈ {−1, 0, 1, 2}, Z 2 , this implies: Now, we have a product of two separate infinite sums. before we proceed, we need to define the ±1/2 'logical' wavefunctions.
Next, we can expand the wavefunction in c and µ, resulting in 4 terms which are superpositions of Let's name these terms {T i }, so that Φ before shift (x) ∝ i T i . In the limit where the normalization constants are ap- Note that the all terms T i have width change from (∆, κ) → (∆/ √ 2, κ √ 2). Also, T 0 correspounds to the codeword with no logical error and T 1 is a bit flip error of logical 0 to 1 and vice versa, which we refer to as a Pauli error. Finally, T ±1/2 are shifted codewords of ∓ √ π/2 in the envelope and the locations of each individual Gaussian peaks, referred to as non-Pauli errors. Each T i term has coefficients based on ψ , and for a likely measurement result x m , one of these terms dominates.
It is useful to see how the q-SE procedure produces Pauli and non-Pauli error terms and also the width change from (∆, κ) → (∆/ √ 2, κ √ 2), no matter what the measurement conditional shift is. If the p-quadrature syndrome extraction procedure is performed on inputs states and ancilla of (∆, κ), one could show that the output wavefunction will also consist of Pauli and non-Pauli terms, but the width of the output GKP state would be (∆ √ 2, κ/ √ 2), this looks like the inverse width change compared to the q-SE. This suggests, by performing the q-SE with (∆, κ) width qubit and ancilla states and then the p-SE with (∆/ √ 2, κ/ √ 2) width ancilla states, one might be able to reverse this width change induced by the syndrome extraction. The effect of width reduction through syndrome extraction cannot be observed in the ∆ = 0 unphysical states calculation.

Justification for fstep
Looking at the equations T i , we can see a reason for using the f step terms in the shift function. The relative amplitudes between terms are: ψ We can see that if one of the T i is much bigger than the others, then that term dominates. For example, when x m ≈ 2n √ 2π, n ∈ Z, then T 0 dominiates.  1 2 , 1, 3 2 }, this gives you the exact amount to shift as stated in the table I.

q-SE followed by p-SE with errors
Given an input state of ψ ∆,κ α (x), undergoing an error shift in phase space by u, v in the q, p quadratures: by carrying out the syndrome extraction procedures, the output state has an exact decomposition in terms of Pauli/non-Pauli errors and has the same width as the input state (∆, κ).
Definition 2. The Fourier transform wavefunction of an arbitrarily shifted GKP state (α ∈ R + ) is: The output wavefunction of the combined q-SE and p-SE is: We can see from the forms of the L β , with errors (u, v) in phase space, it is more difficult to isolate the dominant L β term. the relative amplitudes of L β depends on xm Hence, we need a way to estimate the probability distribution of the errors u, v given x m , p m , σ 0 , namely P(u|x m ) and P(v|p m ).

Bayesian estimator for q-SE shift
We now seek an optimized estimator for the total displacement after a q-SE step. In particular by using knowledge of our error model, we derive an improvement upon the GK shift function in equation 1.
The probability distribution of measurement outcome x m given error shift u can be approximated by: Note that the probability distribution is proportional to a wavefunction not the absolute square of it here. We shall give the proof to equation 20 now.
Theorem 2. The probability distribution of the measurement result x m in the q-SE circuit given an initial arbitary shift u in the q-quadrature is Proof. The probability distribution can be computed from α T α , with T α from equation 9 to 12, with α ∈ {0, 1, − 1 2 , 1 2 }.
For ∆ 1 √ π , the GKP peak don't overlap sufficiently, hence a valid approximation is: ignoring all the cross terms in the sum. Since all the T α terms can be written as a product of functions in x m and x, i.e. |T α | 2 = |A α (x m )| 2 · |B α (x)| 2 , the product function of |B α (x)| 2 integrates to 1. Note that, With a bit of reindexing, and The non-Pauli terms gives the logical 1 component and the Pauli terms returns the logical 0 components.
Using the marginal rule, one can compute where P(u) ∝ e − (u) 2 2σ 2 0 . Using Bayes' theorem, From the plot in figure 3, one can see why the GK shift has this structure of a sawtooth offset by a step function. The sawtooth part aims to counteract the effects of the shift u and the integer steps bit arise from the rotational structure of the QEC code.
We shall maximise the P(u|x m ) now: For a particular value of x m and u, given u is small, Maximising the log of max(P(u|x m )) gives an estimator of u, we shall call the estimator:ũ.
Now we have the sawtooth structure, shown in figure  3, extracted from maximising the conditional probability distribution.
Alternatively, if we expand the sums involved in P(u|x m ) and only concentrate on the central rod, we can extract the mean slope of the diagonal "rods" in figure 3. This can be used to form a new shift function that adaptively changes depending on the parameters σ 0 and ∆ to achieve a optimal code for our error model. The mean slope is: For unphysical GKP states, S(∆=0,σ0) 2 = 1 √ 2 , which saturates the slope value of 1 √ 2 provided by the striaght line in the GK shift. Also, the mean square error for the slope estimation is the width of the rods in figure 3. This works out to be: To optimise the shift function more, when the rounding function outputs 3, take the output -1 instead. This will minimise the overall magnitude of the shift function and hence preserves the overall 1/κ width of the physical state more. With such a rounding function we shall call it round*. We define the BE shift function f BE as: The sawtooth term in equation 30 aims to counteract the u 2 − x s term in the argument of the L β .

q-SE & p-SE combined Bayesian estimator shift
Theorem 3. The probability distribution of the measurement results x m , p m after the q-SE followed by p-SE given an initial arbitary shift of u, v in the q,p quadrature is Proof. Similar to the proof of theorem 2. We shall illustrate the sketch of the proof here. The probability distribution can be computed from α L β , with L β from equation 16  Since ψ ∆ 0 (z) and ψ ∆ 1 (z) are approximately orthogonal, we can ignore terms in the expansion of L β 2 that involves as they integrate to approximated zero. Collecting terms, use the fact that |a 0 | 2 + |a 1 | 2 = 1 and the identities from equation 21 and 22, we can arrive at: Theorem 3 implies: with, P(p m |v) differs from P(x m |u) by the transformation of ∆ → 2∆, κ → κ/2, x m → p m and u → v. Since u, v are both a Gaussian random variable with width σ 0 , hence, we can quickly write down the optimal slope for the p-SE step by using the same shfting function as the q-SE stage but substituting ∆ → 2∆.

MULTI-STEP SYNDROME EXTRACTION
We aim to derive a Bayesian estimator for multiple steps (M steps) of Gaussian q-and p-quadratures shift errors (width σ 0 ) followed by GKP syndrome extraction under the assumption: σ 0 ∆ 1. This translates to looking at the behaviour of our estimator for small Gaussian errors with good GKP states. Our intuition of GKP states tells us that under small errors and SE, at least in the single shot evolution of the wavefunction, we only need to keep track of a reference point of the GKP wavefunction with respect to the origin in phase space to perform the correction.
Let's assume prior to the p-measurement in a one-step q-then p-SE, the wavefunction is (Q ∆ (x) being the qubit wavefunction): where x labels the qubit mode and z is the mode to be measured in the momentum measurement. Applying a Fourier transform to the variables (x, z) → (p, p m ), together with a change of dummy variable in the Fourier integrals, w = x 2 + z √ 2 and y = x 2 − z √ 2 , results in: Here we have used, which comes from equation 15, moving the 4 terms from the reindexing inside the argument of Q as displacements. This is valid if we assume the envelope width to be large, i.e ∆ 1. Note that the f * step is the f step function with the rem * function embedded. If the displacements u < √ π/4, then the f * step (x m − √ 2u) ≈ f * step (x m ) due to the nature of the round function. If we define θ(x m , u) as: Then, Using reindexing and hiding the four terms as displacement approximately: (41) If we now inverse Fourier transformW(p, p m ) w.r.t p: We can see that the quantities θ(x m , u), θ(p m , v) represents the refrence frame of a GKP state after errors and syndrome extraction. Suppose we have error shifts first then followed by a round of q-then p-SE with no corrective shift at each individual round. Iterate this M times. Then we can derive sets of recursive relationship fo the reference frame in the M th step. For notational convenience, we denote the error shifts of the start of the h th as (u h , v h ) and measurement results as (x If we define the quantities: then the q-and p-quadrature reference frame of the qubit after h steps of error correction is: This is assuming that f step is independent of u k , which is true (due to the nature of the round function) if: h−1 round and all the displacement shifts of all h rounds can be written down easily analogous to equations 20.
Using Bayes' Theorem to flip this distribution for M steps to get: We can see that the q-and p-quadratures equations have the similar forms and can be inter-changable (q ↔ p) if u ↔ v, x m ↔ p m and (∆, κ) ↔ (2∆, κ/2). So we will focus on the q-quadature version for now. Maximising the logarithm of P (q) M ( u| x m ) gives the Bayesian estimator solutions (ũ k is the estimate for u k ).
Small errors: with F h = √ 2X where M is the Bayesian estimate of the error shift u and Σ is the covariance matrix. If we rewrite the terms in the exponent of max(P  We have now characterised the multivariate Gaussian noise model (multivariate because we have many steps of Gaussian errors) with its covariance matrix. If we look at the final displacement, we have a correctable part that depends solely on the measurements results and another part that is a random variable. We wish to find the uncorrectable part's variance, i.e V q = Var( We can see that V q → σ 2 0 3 , very quickly as M grows. We need to stress that V q is the variance of the Bayesian estimator on the final shift. The errors on the qubit can be estimated more precisely using the measurement results from all rounds of syndrome measurements and processing them accordingly. Now we can finally use the expansion of the covariance matrix combined with −2 1