Adaptive Bayesian phase estimation for quantum error correcting codes

Realisation of experiments even on small and medium-scale quantum computers requires an optimisation of several parameters to achieve high-fidelity operations. As the size of the quantum register increases, the characterisation of quantum states becomes more difficult since the number of parameters to be measured grows as well and finding efficient observables in order to estimate the parameters of the model becomes a crucial task. Here we propose a method relying on application of Bayesian inference that can be used to determine systematic, unknown phase shifts of multi-qubit states. This method offers important advantages as compared to Ramsey-type protocols. First, application of Bayesian inference allows the selection of an adaptive basis for the measurements which yields the optimal amount of information about the phase shifts of the state. Secondly, this method can process the outcomes of different observables at the same time. This leads to a substantial decrease in the resources needed for the estimation of phases, speeding up the state characterisation and optimisation in experimental implementations. The proposed Bayesian inference method can be applied in various physical platforms that are currently used as quantum processors.


Introduction
Quantum computers have the potential to solve some computationally hard problems in a more efficient way than classical computers [1]. However, due to coupling with the environment, they are more susceptible than their classical counterparts to dynamical errors that affect the correct behaviour of the algorithms performed [2]. In order to cope with dynamical errors, quantum error correction techniques [3] need to be applied together with a correct initialisation of quantum states that, in general, suffers from different types of noise. These imperfections can often be modelled as irreversible couplings to the environment [4] or as unknown but constant unitary operations appearing due to systematic errors. Due to their constant nature, the latter can be compensated by determining the unknown operations and applying their inverse onto the state. The simplest instance of such systematic errors is given by single-qubit phase shifts which can transform a desired state α |0 + β |1 into α |0 + βe iφ |1 where φ is an unknown but constant phase. Estimates of this phase shift can be obtained by performing Ramsey-type experiments [5,6] and, more recently, by adaptive methods based on application of Bayesian inference [7][8][9][10][11][12][13][14]. These adaptive methods select the measurement to be performed by numerical optimisation of the information gained based on the results obtained in the previous measurements.
The characterisation of multi-qubit states, such as those needed for quantum error correcting codes [15], is a more complex problem. It is well known that quantum state tomography [16] becomes impractical to fully characterise these states since the resources needed scale exponentially with the number of qubits. Additionally, the systematic errors to be corrected can drift slowly over time. Thus, the error estimation needs to be performed on time scales smaller than the drift time in order to correct the errors before the estimates become obsolete. However, as adaptive techniques can take advantage of experimental data collected at each measurement step, they can be successfully employed for increasing the information obtained with each measurement and thus decreasing the amount of resources needed as compared to non-adaptive techniques.
Along this line, in this work, we propose and explore an adaptive method based on the application of Bayesian inference for the characterisation of medium-scale multiqubit states. For concreteness, we focus on the estimation of phases appearing in stabiliser states used in quantum error correction. We examine the efficiency of the method by studying the number of measurements needed and we derive an analytical rule to obtain the optimal measurement to perform at each time. This simple rule does not rely on numerical calculations and ensures the adaptiveness of the method to find the optimal measurement at each step. This renders the protocol particularly suitable for on-chip processing in adaptive control systems.
In addition, we evaluate the efficiency of our adaptive method and compare it to that of the Phase Optimisation Method (PHOM), a non-adaptive phase estimation method developed in [17]. The latter method is based on a generalization of a Ramsey experiment to determine and compensate systematic phases appearing in multi-qubit states and it was proposed and performed experimentally for the seven-qubit quantum error correcting color code (Steane code) [18]. For this reason, we choose to evaluate the efficiency of our adaptive method for the phase characterisation of the states used for the Steane code and we show that our method has an improvement in the efficiency compared to the PHOM.
Remarkably, this adaptive method is not restricted only to the Steane code, but can be used for other multi-qubit states since it only relies on the application of simple single-qubit operations and measurements. As a consequence, the results shown in this paper are applicable to the optimisation of other QEC codes and the compensation of systematic errors appearing in other physical platforms for quantuminformation processing such as, e.g, trapped ions [19][20][21][22], Rydberg atoms [23][24][25] in optical lattices [26][27][28] or tweezer arrays [29,30] or other AMO or solid-state architectures [31][32][33][34][35][36]. This paper is organised as follows: In Section II we introduce the concepts and notation for a one-qubit state phase measurement by using a Ramsey experiment and a Bayesian inference process. In Section III we briefly review basic properties of the Steane code to which we will apply our technique. In Section IV we compare the efficiencies of the PHOM proposed in Ref. [17] and our Bayesian inference method to an intermediate quantum state obtained during the preparation of the logical states of the Steane code since this intermediate state has a less complex structure than the final states of the code. In Section V we generalise the previous results and present a Bayesian inference method to estimate the phase shifts on the fully encoded sevenqubit logical states. Finally, in Section VI we summarise our results, especially the comparison between the results obtained for the efficiency of our method with the method in Ref. [17], and conclude with a brief outlook.

One qubit case
In this Section, we will show how to estimate the unknown phase φ of the following quantum state from a finite set of data obtained from measurements. We suppose we can prepare as many copies of |ψ as needed. Since measurements of theẐ Pauli operator yield no information about φ we will perform measurements of the operator on the XY plane of the Bloch sphere (see Fig. 1 (a)). Thus, the expected value of this operator for the state |ψ is In the following, we will study two different ways in which we can select θ in Eq. (2) in order to obtain the value φ: Ramsey scan and a Bayesian inference process.

Ramsey scan
In order to estimate φ, one can apply a Ramsey-type experiment that can be summarised as follows (see Fig. 1): First we divide the interval [−π, π) in M equidistant points For each of the values θ m we estimate the expected value ofÔ θm by ψ|Ô θm |ψ = (N  (1). The observablesÔ θ are obtained by rotating the observableX around the Z axis by an angle θ. The waŷ O θ is selected depends on the method used to estimate φ. In the Ramsey scan case different values of θ are selected and several measurements are performed for each of these values. This allows a reconstruction of the expected sinusoidal dependence of ψ|Ô θ |ψ with θ that yields an estimate for the phase φ [panel (b)]. In the Bayes case, the outcomes ofÔ θ are used to update the probability distribution P (φ). In this case, the value θ for each measurement is different and it is selected in a way that maximises the information gain per measurement (See Sec 2.2 and Fig. 2). (b) Ramsey scan simulation for the phase estimation of the state |ψ in Eq. (1) with φ = 1. The points represent the valuesÔ θm where the θ m are M = 10 equidistant points in the interval [−π, π). The points are obtained by simulating N = 50 measurements ofÔ θm for each θ m and they are used to fit a cosine whose phase is the estimation obtained for φ.
measuring in theÔ θm basis and N is the total number of measurements. Using the values ψ|Ô θm |ψ obtained, we perform a least squares fit to a function of the form where A is a phase. By comparison of Eqs. (3) and (5), our estimate of φ is given by the fitted parameter A.

Bayesian inference
In this section we discuss how to use Bayes' theorem in order to estimate the value of φ. Bayes' theorem prescribes how to update the prior probability distribution of φ, P (φ), after a measurement M using its likelihood, P (M|φ) (see Fig. 2). The final result is a posterior probability distribution of φ, P (φ|M), with the form up to a normalization factor. If we perform a new measurement we can apply Bayes' theorem again and use the obtained posterior as a prior for the next measurement and so on. For an increasing number of measurements the degree of uncertainty of φ will decrease, allowing us to reach a desired value in the uncertainty of the estimated value of φ. In our case we measure the operatorsÔ θ (with different values of θ for each measurement) with possible outcomes + θ and − θ . The likelihoods of these outcomes for the state as given by Eq. (1) are Assuming no prior knowledge about the value of φ we start with a uniform probability distribution P (φ) = 1/2π as a prior. After N measurements and applying Eq. (6) iteratively, the probability distribution for φ is given by As the number of measurements increases, the posterior probability distribution can be approximated by a normal distribution with decreasing standard deviation.
2.2.1. Efficiency of the parameter learning process. In this section we show how the behaviour of the standard deviation of the distribution in Eq. (8) allows us to choose the value θ after each measurement in order to maximize the information gain of φ. It is expected that as we perform more measurements, the mean value of the probability distribution Eq. (8) gets closer to the true value of φ and its standard deviation decreases. Let us suppose that, after a sufficiently large number n of measurements, the probability distribution P n (φ) can be approximated by a Gaussian with meanφ n and standard deviation σ n . For the measurement n, the probability p ± n of measuring ± when the angle selected is θ n is The probability distribution after having obtained + or − is updated as where p ± n appears due to normalization. These posterior probability distributions will have a standard deviation denoted by σ ± n+1 . We obtain that the average decrease of the variance after performing measurement n is with α n ≡ e σ 2 n sin 2 (φ n − θ n ) 1 − e σ 2 n cos 2 (φ n − θ n ) .
From inspection of Eqs. (11) and (12), we conclude that the maximum decrease on average for the variance is obtained when we select a value of θ n that maximizes the value of α n . This is achieved (see Fig. 3) for (1) with φ = 2. Each step of the inference process consists in performing a measurement of the operatorÔ θn−1 , updating the probability distribution based on the result of the measurement and selecting an optimal angle θ n for the next measurement. 1) One starts with a uniform probability distribution of φ. θ = 0 is selected for the first measurement. 2) If the + outcome is obtained (as assumed here), the probability distribution is updated by multiplying the prior probability distribution by the likelihood of obtaining a + and renormalizing. The maximum of P (φ|1 measurement) is located at φ = 0. 3) The optimal selection of θ will be given by Eq. (13). Thus, we select θ = π/2 for the next measurement. 4) Assume + is obtained again. The probability distribution can be updated again based on this result. This process of measuring, updating and finding the next optimal θ can be performed iteratively. 5) After 10 measurements the probability distribution can be approximated by a normal distribution. 6) After 500 measurements we obtain a normal distribution centered near the value φ = 2 used for the simulation.
For this selection, the value of α n approaches a constant value. Then, it can be proven that the succession in Eq. (11) has the following asymptotic solution: However, as the variance decreases, α n approaches the constant value 1 (except for values close to θ n = φ n ± kπ, k ∈ N) as Fig. 3 shows. This means that after several measurements, the decrease of the variance will be independent on the value of θ we select for the next measurement and σ 2 n will evolve as Similar results can be obtained by performing calculations involving Shannon's entropy [37] for the selection of the optimal measurement. However, this approach requires numerical calculations to obtain the optimal measurement in each step of the iteration. This increases the computational resources needed for in-situ optimisation between Other values ofφ n produce the same plot with a translation on the horizontal axis. As σ 2 n decreases α approaches 1 except for the values θ =φ n ± kπ. experimental measurement runs as compared to the simple rule governed by Eq. (13) of our method. It is worth mentioning that the result in Eq. (15) indicates that our increase in the knowledge of the system satisfies the Standard Quantum Limit (SQL). Although there are methods that aim to obtain better results than the ones given by this limit [8,[38][39][40], these methods rely on having access to, e.g., entangled states, which we will not consider in this work.

Characterisation of multi-qubit states
In quantum error correction the quantum information is encoded in entangled manyqubit systems. This provides protection against noise.
In the following sections we will discuss the estimation of systematic errors appearing in the preparation of the stabiliser states used in quantum error correction. For concreteness, we focus on the states used for the seven-qubit Steane code [18]. This code represents the minimal instance of 2D color codes [41] and it is obtained by restricting a the Hilbert space of seven qubits to the subspace of states which are simultaneous +1 eigenstates of six commuting stabiliser operators S (i) Fig. 4). These stabilisers define a 2-dimensional subspace for this 7qubit system that can be used to encode a logical qubit. Additionally, the logical X and Z operators can be chosen as Similarly, the logical |1 L is defined by Z L |1 L = − |1 L . These states satisfy Z for j = 1, 2, 3 associated with each of the four-qubit plaquettes. The code space is defined as the simultaneous +1 eigenspace of these stabilisers.
In the following sections we develop a method to measure the phases appearing due to systematic errors in the preparation of this class of states. For simplicity, we will first study an intermediate case of the full 7-qubit Steane encoding process (twoplaquette case) to introduce the concepts that will be needed to correct the phases appearing in the fully encoded system (three-plaquette case).

Two-plaquette case
At the start of the preparation of the seven-qubit quantum error correcting code, four-qubit entanglement operations are applied to the first plaquette ( Fig. 5 (a)). This yields the quantum state ⊗7 composed by the superposition of two components in the computational basis which can have a relative phase due to systematic errors. This is equivalent to a single-qubit phase estimation, as the phase can be corrected by rotating one of the four qubits and performing measurements of the S (1) x stabiliser. Therefore, we consider the state |ψ 2 obtained by application of four-qubit entangling operations to the first and second plaquettes (see Fig. 5 This state maximizes the mean value of the X-type stabilisers on the first and second plaquette, S However, systematic phase shifts accumulate during the preparation of |ψ 2 due to experimental errors. The state |ψ 2 containing these unknown phase shifts is x for different rotations on the first qubit. This is similar to the single-qubit phase estimation. (b) Two-plaquette case. After the manipulation of the first and second plaquettes, up to three relative phases can appear in quantum state obtained. To estimate these phases we can perform measurements of S for different rotations on the first, second and fifth qubit.
In order to compensate these relative phase shifts we can apply single qubit Z rotations (see Fig. 5 (b)). For example, by rotating the first, second and fifth qubits we obtain The selection of qubits to be rotated is arbitrary as long as these three rotations do not commute with S (1) x , S (2) x and S (1) x S (2) x , respectively. The expected values of the stabilisers in the state Eq. (19) are given by In order to obtain information about the unknown systematic phases, we can perform measurements of these stabilisers for different values of the rotation angles θ. Once the values of these phases are measured it is possible to perform single-qubit rotations to transform the state |ψ 2 into the desired state |ψ 2 . A way to obtain these values is the Phase Optimisation Method [17]. We propose another method based on application of Bayesian inference. In the following subsections we review the PHOM and introduce our Bayesian protocol.

Phase Optimisation Method
The Phase Optimisation Method introduced in Ref. [17] is given by the following iterative protocol. For concreteness, here we review how it works for the optimisation of the state in Eq. (19).
(i) For each stabiliser, an associated rotation on a qubit i, θ i , is chosen. The selection is arbitrary, but each stabiliser must not commute with its associated qubit rotation. We associate θ 2 with S (1) x , θ 5 with S x and θ 1 with S x S x . (ii) Choose an initial configuration for the set of rotation parameters in a similar way as in the single-qubit case (see Fig. 1 (b)) over its associated angle, θ 2 , in the interval [−π, π] while keeping θ 1 = θ This method gives an estimate of the angles θ 1 , θ 2 and θ 5 that correct the systematic phases φ 1 , φ 2 and φ 3 . The precision of this estimate will become better as the number of measurements used for the scans increases. In this work we also introduce a variation of this PHOM to measure phases that we call the constant cosine PHOM.
4.1.1. Constant cosine PHOM. The constant cosine PHOM is a similar method that also performs scans of the expected values of the stabilisers for different qubit rotations to obtain a correction for systematic phases. This process is more similar to a Ramsey experiment as only one scan of each stabiliser is needed. From Eq. (21) we can see that, if we keep the value θ 2 − θ 1 fixed and vary θ 2 + θ 1 , the mean value of S (1) x will be given by where h is constant for all the measurements since the difference θ 2 − θ 1 is fixed. Thus the angle φ 2 that represents the phase shift to be corrected is given by the value of −2(θ 2 + θ 1 ) for which a maximum in the mean value S

Bayesian inference method
In the following we will introduce and analyse two Bayesian inference methods to measure the phases in the state |ψ 2 of Eq. (19) of the two-plaquette case. These methods are (i) a Bayesian inference method by direct application of the likelihoods, and (ii) a Bayesian inference method using marginal likelihoods. We will first describe this direct Bayesian inference method and then explain the method we propose to improve the PHOM, namely the marginal likelihood method.

Direct Bayesian inference method.
An estimation of the phases using Bayesian inference is performed by measuring the plaquettes and updating the probability distribution based on the results obtained. The likelihoods for each plaquette measurement can be obtained from the expressions for the expected values by where P 1 (± θ |φ) is the likelihood of obtaining a + or a − outcome when measuring S (1) x . For instance, the likelihood for S (1) x is given by Since the likelihoods used are functions of three variables, the obtained probability distribution will be a three-dimensional function. In general, if the number of unknown parameters appearing in the likelihoods of the experiment increases, the probability distribution obtained will be a function of many variables and, therefore, finding the most probable values for the phases and their variance becomes more difficult. This complication can be avoided if in the measurement of each stabiliser we keep one of the cosines constant in a similar way as it is done for the constant cosine PHOM. This yields a likelihood given by (see also Eq. (24)) for the first plaquette, where the value θ 2 − θ 1 is kept constant to ensure that one of the cosines has a constant value given by h 2 . Similar expressions can be obtained for the other stabilisers. This approach yields normal probability distributions defined on two variables, one being h 1 , h 2 or h 3 and the other being φ 1 , φ 2 or φ 3 depending on the stabiliser that is measured. The estimate for each phase is easily obtained from its corresponding probability distribution. We now compare this method and the constant cosine PHOM by simulating them and, by fitting the numerical data we obtained, we find that for both, the scaling of the variance as a function of the number of measurements n is given by σ 2 i,n = 6/n when estimating the single angle φ i (see Fig. 6). However, in order to obtain an estimate of the other two phases, this process must be repeated for each of the other stabilisers. This yields a scaling for the variance σ 2 n that scales as that gives an estimate of the efficiency of the PHOM and the direct Bayesian inference method for the intermediate state |ψ 2 .
In the following, we will introduce and analyze the Bayesian marginal likelihood method, which constitutes an improvement in the efficiency both of the PHOM and of the direct Bayesian inference technique.

Marginal likelihood method.
Let us consider measurements of the first stabiliser whose likelihood is Suppose we are only interested in φ 2 and we perform measurements with values of θ 1 and θ 2 selected randomly. With this selection the cosine containing φ 1 − φ 3 has a completely random argument and it averages to zero. This yields the following marginal likelihood which only depends on φ 2 . This likelihood is similar to the likelihood for the one qubit case of Eq. (7) so we can generalize the analysis for the scaling of the variance for φ 2 (see Eq. (11) and Appendix B) obtaining where andθ 2,n = −2(θ 2,n + θ 1,n ). For small values of σ 2 2,n , α 2,n can be approximated by (see Fig. 7) The quantity α 2,n oscillates between 0 and 1/4 and has a mean value equal to α 2,n = (2 − : Dependence of the factor α i,n on θ i,n for i = 1, 2, 3 (hereφ i,n = 0). α i,n dictates the decrease of the variance step by step. α i,n oscillates between 0 and 0.25 and its mean value is (2 − √ 3)/2 ≈ 0.134, with its maximum at θ i,n =φ i,n ± π/2. By Eq. (14), this is the optimal way of selecting θ n . distribution ofθ 2,n , since we are using random values ofθ 2,n . Taking into consideration Eq. (14) the scaling of the variance is given by A similar derivation can be performed for the likelihoods of the other stabilisers, yielding the same behaviour. These results are obtained for non-adaptive measurements since so far the experimental configuration used is not selected in a way that maximizes the information gained by each measurement (the θs are selected randomly). If one choosesθ i,n =φ i,n + π/2 for i = 1, 2, 3 this causes a displacement h i to appear in each likelihood, e.g., in the likelihood of φ 2 : In order to correct this displacement in the likelihoods a random selection betweeñ θ i,n =φ i,n +π/2 andθ i,n =φ i,n −π/2 can be done. As a consequence the displacements will alternate between the value h i and −h i and its average will be 0. Thus, the same likelihood as for the random angles selection case is obtained, but the factor α i,n improves to 1/4. Thus, the variances scale as σ 2 i,n = 4/n. A pseudocode that summarises the steps presented here is shown in Fig. 8.
Using these simple rules for the selection ofθ i,n ensures that the adaptive way of selecting the parameters of the measurements yields more information than selecting them in a non-adaptive way. Additionally, in each measurement one can use the value of all three stabilisers as opposed to what happens in the PHOM and the Bayesian inference with the constant cosine, in which each measurement only yields information about the scanned stabiliser. Thus the estimate of each phase will still have a variance Input: n copies of the state |ψ Output: List of estimated phases φ = { φ 1 , φ 2 , φ 3 } and corresponding probability distributions where each β j is chosen randomly between {+π/2, −π/2} end for 1 Figure 8: Pseudocode to determine the three phases, φ = {φ 1 , φ 2 , φ 3 }, appearing in the two-plaquette case of Sec. 4 by implementation of the marginal likelihood Bayesian inference method. This pseudocode can also be applied for estimating the seven phases appearing in the three-plaquette case of Sec. 5 after redefining φ, the list of unknown phases and their indices J; P, their probability distributions; θ, the angles for the rotations of the measurements; Π, the likelihoods; S the list of stabilisers; Q the list of qubits to be rotated; the system of equations S according to (A.20) in Appendix A.

Three-plaquette case
The methods seen in the previous sections can be generalized to the more complex case of the entire seven-qubit code. In this case the objective is to measure the 7 phases that can appear in the preparation of the state that represents the logical state |0 L (see Appendix A). To this end, we can measure seven different combinations of stabilisers: S x , S x , S x , S x S (2) x , S x S x , S x S x and S (1) x . The marginal likelihood for the measurement of the first stabiliser is whereθ 2 ≡ −2(θ 1 + θ 2 + θ 3 + θ 4 ). Similar expressions are obtained for the other six combinations of stabilisers (see Appendix A). Following the same process as in the previous sections, the behaviour of the variance and the values of α for each angle φ i at the step n (α i,n ) are now given by (see Appendix B) with where for small values of σ 2 i,n one finds The maximum value of each α i,n is 1/16 for the selectionθ i,n =φ i,n ± π/2. Thus, this rule ensures an adaptive way of selecting the measurements that yields more information than a non-adaptive selection and gives a scaling as for the variance of the estimate for each phase measured.

PHOM simulations
In this section we discuss numerical simulations to obtain the behaviour of the variances for the PHOM with the number of measurements used in the three-plaquette case.
We initially fix the number of iterations I and we choose a number n of copies of the initial quantum state that we can measure. We then choose N different 7-dimensional vectors representing the seven initially unknown phases φ (k) (k = 1, ..., N ). We perform the PHOM following Sec. 4.1 and the expressions (A.4) to (A.11) in Appendix A in order to obtain an estimate of the seven phases φ (k) est . For each of the vectors, we calculate the difference between the estimated phases and the phases chosen initially: φ (k) est − φ (k) . We finally compute an estimate of the variance of the PHOM for a given n and a fixed number of iterations I as where N is a large number to obtain a good estimate of σ 2 n , in our case N = 50000. To perform the scans of each stabiliser we divide the intervals [−π, π] of the angles θs into M = 10 points. Since we have to measure seven stabilisers, the number of measurements per point (mpp) we can perform is mpp = n/(7M I). Figure 9: Behaviour of the variance as a function of the number measurements, n, obtained from simulations. For the PHOM using constant cosines a behaviour of 224/n is obtained. PHOM shows a similar behaviour when using enough iterations for the method to converge. Finally, using marginal Bayes the performance is 16/n. Repeating this process using a different n and I = 1, . . . , 4 yields the data plotted with circles in Fig. 9. A similar process is done for the constant cosine PHOM without performing iterations since they are not needed for this method (see Sec. 4.1). The results obtained for this case are represented by the black circles in Fig. 9. It can be seen that as the number of iterations used for the PHOM increases, the variance with the number of resources decreases until reaching a similar behaviour as that of the constant cosine PHOM. The reason for this is that the PHOM has two limitations, one given by the number of iterations and another by the finite number of measurements used for each scan. If a small number of iterations is used, most of the PHOM simulations do not converge and the result is a wrong estimation for the phases which results in a big variance. As the number of iterations increases, more simulations converge to a correct phase and the variance obtained decreases. After performing enough iterations, all the simulations converge and the only source of error is the number of measurements used for each scan. As this number increases, the statistical error of each scan performed decreases yielding a better estimate of each phase. This behaviour is represented in Fig. 10.
We perform a fit of the relation between the variance and the number of measurements for the constant cosine PHOM and the PHOM when enough iterations are performed for it to converge. This fit yields σ 2 n ≈ 224/n for these cases.

Marginal likelihood Bayes simulations
To obtain the behaviour of the variance with the number of measurements for the marginal likelihood Bayes inference method we apply a generalisation of the method shown in the pseudocode of Fig. 8 with the expressions (A.12) to (A.20) in Appendix A. After performing enough measurements (≈ 100) the probability distribution obtained can be approximated by a normal probability distribution for each of the phases that is used for computing the variances and for estimating the error of this adaptive technique. The results for the variances are shown in Fig. 9 where the data is plotted with blue squares. It is possible to see that the variance decreases as the number of measurements used increases. A numerical fit of the data obtained reveals that the variance decreases as σ 2 n ≈ 16/n, as was expected from the analytical derivation of the method in Sec. 5.

Conclusions and outlook
In this work, we have introduced an adaptive Bayesian method to measure systematic phase shift errors appearing in the experimental preparation of multi-qubit states, in particular, we have used it for correcting the errors appearing in the logical states of the Steane code. This method is capable of finding the experimental configuration that optimises the information gained by each measurement performed. An analytical development of this method that yields a simple rule for this adaptive selection has been shown, thus saving computational power that would be otherwise used in finding numerically the optimal measurement at each step of the process. We compared our method to the Phase Optimisation Method (PHOM), a nonadaptive phase estimation method based on a generalization of a Ramsey experiment for multi-qubit states that was recently realised for the implementation of the Steane code [20]. We simulated both of them to measure quantum phases appearing in the preparation of quantum states needed for the Steane code. The efficiency obtained by simulation of our method is in agreement with the efficiency derived from the theoretical calculations, and shows an improvement by a reduction of the required measurement time by more than one order of magnitude when compared with the efficiency of the PHOM (see Fig. 11).
Furthermore, the method, illustrated for the optimisation of the seven-qubit Steane code, is applicable to other QEC codes and stabiliser states. Additionally, since this method only relies on the application of single-qubit rotations and measurements, it can correct systematic errors appearing in multi-qubit states implemented in different physical platforms. Thus, it has potential application in a variety of systems for quantum information processing such as, e.g., trapped ions, Rydberg atoms in optical lattices or tweezer arrays or other AMO or solid-state architectures.
Appendix A. Three-plaquette case: Quantum states, likelihoods and pseudocode generalisation In this Appendix we provide details on the expressions appearing in the three-plaquette case. For the seven-qubit Steane code the logical 0 is given by Due to experimental errors, phases appear in the preparation of the |0 L state. The state obtained will be then Single-qubit Z rotations can be performed on the seven qubits of the code state |0 L to obtain the following state and likelihoods for each stabiliser combination The expected values of each plaquette can be easily obtained from these expressions. For example, the first plaquette expected value is given by The expected values for the other plaquettes can be obtained from the corresponding likelihood in the same way. The marginal likelihoods are given by As for the generalization of the pseudocode in Fig. 8, it is achieved by changing the previous definitions to the following ones J = 1, ..., 7, x }, Π = {P 2 (± θ |φ 1 ), P 1 (± θ |φ 2 ), P 12 (± θ |φ 3 ), P 3 (± θ |φ 4 ), P 23 (± θ |φ 5 ), P 13 (± θ |φ 6 ), P 123 (± θ |φ 7 )}, where the β i in S are chosen randomly from the set {+π/2, −π/2}.

Appendix B. Analytical study of the scaling of the variance
In this Appendix, we present some details on the scaling of the variance with the number of measurements when applying the Bayesian adaptive method to the singlequbit state of Eq. (1). The expressions obtained are easy to generalize to the two and three-plaquette cases. Let us suppose after n measurements the knowledge about the phase φ is given by a normal distribution with mean value φ n and variance σ After performing a measurement the updated probability distribution is given by where p ± n is the probability at step n of obtaining a + or − in the measurement n + 1. Since at step n it is unknown what the measurement n + 1 will yield, we consider the expected value of the variance after the measurement n + 1, σ 2 n+1 where (σ ± n+1 ) 2 are the variances after obtaining a + or − for the measurement n + 1 given by

(B.12)
Similar calculations can be performed for the two and three plaquette case likelihoods, the only difference being the cosine appearing in the likelihood having amplitude 1/4 and 1/8 respectively and the angleθ i,n being a linear combination of the rotations performed on different qubits at the step n. Taking this into account, for the two plaquette case the decrease in the variance for each φ i is σ 2 i,n+1 − σ 2 i,n = − e −σ 2 i,n σ 4 i,n sin 2 (φ i,n − θ i,n ) 4 − e −σ 2 i,n cos 2 (φ i,n − θ i,n ) (B.13) and for the three plaquette case σ 2 i,n+1 − σ 2 i,n = − e −σ 2 i,n σ 4 i,n sin 2 (φ i,n − θ i,n ) 16 − e −σ 2 i,n cos 2 (φ i,n − θ i,n ) . (B.14)