Estimation of coherent error sources from stabilizer measurements

In the context of measurement-based quantum computation a way of maintaining the coherence of a graph state is to measure its stabilizer operators. Aside from performing quantum error correction, it is possible to exploit the information gained from these measurements to characterize and then counteract a coherent source of errors; that is, to determine all the parameters of an error channel that applies a fixed - but unknown - unitary operation to the physical qubits. Such a channel is generated, e.g., by local stray fields that act on the qubits. We study the case in which each qubit of a given graph state may see a different error channel and we focus on channels given by a rotation on the Bloch sphere around either the x, y or z axis, for which analytical results can be given in a compact form. The possibility of reconstructing the channels at all qubits depends non-trivially on the topology of the graph state. We prove via perturbation methods that the reconstruction process is robust and supplement the analytic results with numerical evidence.


I. INTRODUCTION
In order to construct devices that reliably process quantum information it is necessary to characterize, measure and eventually remove all the known sources of detrimental noise [1][2][3]. It is however inevitable that a certain amount of disturbance will continue to affect the apparatus. Hence it usually will be indispensable to perform quantum error correction (QEC) [4][5][6][7][8] to reduce the error rate to a tolerable amount.
Robust and scalable architectures for quantum information processing require that their elementary building blocks are readily calibrated and easily operated. Better scalability with less manual intervention can be achieved by engineering devices which are capable of autonomously characterizing and correcting the error sources. Also, since these systems will likely become increasingly more integrated and miniaturized [9][10][11][12], it will be difficult to have physical access to their building parts once the apparatus has been assembled. Therefore, it is desirable that these devices are designed and built so that it is possible to gather data while the machine is running, and thus, on-the-fly, achieve good control and feedback on the internal functioning.
In this paper we introduce a method that allows for an automated analysis of the presence of error sources and envision an agent, connected to the quantum device, which performs classical statistical analysis of the data in order to learn which noise processes are occurring on all qubits. A related case study of adapting to a unitary channel acting on an isolated qubit is given in Ref. [13]. Furthermore, as we argued before, the device will necessarily be equipped with a QEC code to protect the logical subspace from noise. Hence we set out to investigate the following natural question: to what extent is it possible to characterize the error sources affecting the apparatus, having access to the syndrome measurements alone? More specifically, we examine the scenario in which the agent collects the syndrome measurements and performs statistical analysis over them to detect and characterize a unitary component in the noise. Therefore we avoid the necessity of adding additional sensors to monitor the performance of the quantum processing device. This approach is independent of the physical implementation of the apparatus and it can be studied through a multitude of different methods, see e.g. Refs. [14][15][16][17][18].
Once the unitary component of the error has been characterized it is possible, at least in principle, to completely correct it through the application of a counter-unitary operation, which undoes the effect of this component of the error channel. As an alternative, we might allow the quantum computation to take place in a "floating frame of reference", in which the logical quantum basis change in time in order to compensate for the rotation created by the unitary error channel.
Hereafter we will study the information that can be recovered from the most restrictive type of QEC. We consider the extreme case of a QEC code in which only one state is protected, i.e., the logical subspace has dimension one. This QEC method can be adopted to produce the resource state for measurement-based quantum computation (MBQC) [19,20]. In MBQC particular classes of quantum states, e.g. the cluster states [21], are universal, meaning that they can be employed to perform arbitrary quantum computations just by applying local measurements. The cluster states are a subset of a more general family of quantum states associated to undirected graphs, which are called graph states. A method for obtaining any graph state consists in measuring a complete set of stabilizer operators (a.k.a. correlators) and in applying the required corrections whenever an "error" is detected. The output of this procedure is a unique quantum state, i.e. the required graph state. Another setting in which this QEC can be employed is for measurementbased quantum repeaters [22][23][24], in which specific graph states are used as resources to perform entanglement purification and entanglement swapping to transmit quantum information over long distances.
The outline of the paper is the following. In Sec. II we give a brief introduction to graph states and set the notation; in Sec. III we describe and motivate the error model we analyse in the paper; in Sec. IV we detail the effect of the unitary error channel over the stabilizer arXiv:1512.07083v1 [quant-ph] 22 Dec 2015 measurements statistics; then in Sec. V we describe a special promise setting in which the problem becomes solvable even with limited resources. The central part of the work begins in Sec. VI, where we provide a general method for reconstructing the error channels in the promise settings; then in Sec. VII we apply the general methods given in Sec. VI to some specific classes of graph states; finally in Sec. VIII we show that our reconstruction methods are resilient to small deviations from the model. We conclude by giving in Sec. IX a completely different approach which allows us to reconstruct the unitary error in complete generality, at the cost of requiring some extra resources.

II. PRELIMINARY CONCEPTS AND NOTATION
A graph state |G is quantum state associated with an undirected graph G = (V, E) having vertices V and edges E = {(a, b) | a, b ∈ V }, which connect pairs of vertices; |G is defined as where |+ := (|0 + |1 )/ √ 2, N := |G| is the number of qubits of the graph state (each vertex represents a qubit) and C-Z a,b denotes the (symmetric) controlled-Z operation between the qubits at vertex a and b. The neighbours of any vertex a ∈ V are given by the vertices in the set N a := {b ∈ V | (a, b) ∈ E}. We also use the notation N a := N a ∪ {a}.
We say that two graph states are locally unitarily equivalent (LU -equivalent) if it is possible to obtain one graph state from the other by applying unitary operations acting on individual qubits.
A standard method that allows us to preserve the coherence of a graph state |G consists in repeatedly measuring the correlator operators ∀a ∈ V . Here X a , Y a and Z a denote the Pauli matrices acting on the qubit in vertex a; we also use the shorthand S a := (X a , Y a , Z a ). All these correlators commute pairwise -thus all the measurements can be performed simultaneously -and form a complete set of observables. We explicitly distinguish the correlator operators K a from the random variables κ a which represent the realized measurement outcomes of K a . The variables κ a assume values κ a ∈ {0, 1}, corresponding respectively to the {+1, −1} eigenvalues of K a . The graph state |G is the unique common eigenvector with eigenvalue equal to +1 for all the correlators; equivalently, this is the state one obtains when κ a = 0, ∀a ∈ V . Hence, when the measurements outcomes of the correlators are given by {κ a } a∈V , the system is projected into the state where, for any quantum gate U , U κa = 1 when κ a = 0 and U κa = U when κ a = 1. The graph state |G is equivalent to |G for the purpose of performing MBQC. The Pauli errors Z a can be accounted for by a redefinition of the logical encoding of the qubits: i.e., when κ a = 1, the sign of the state |1 has to be reversed (that is, we make the association |0 L ⇔ |0 and |1 L ⇔ −|1 ). Thus we can assume, without loss of generality, that the graph state is re-prepared after the measurement of all the correlators in the standard graph state |G . As stated in the introduction, from a quantum error correction perspective these correlators provide a stabilizer code which protects only one word, the state |G . For a thorough introduction to these concepts, see Ref. [25].

III. ERROR MODEL
We consider the following model for a unitary error source that acts on a graph state. Suppose that our graph state is physically implemented with spin-1 2 particles, each representing a qubit with the spin degree of freedom [26]. These particles carry an intrinsic magnetic dipole moment; hence they are susceptible to undergo Larmor precession when immersed in a magnetic field, with angular frequency equal to ω = γB, where B is the magnetic field intensity and γ is the gyro-magnetic ratio of the particle. However, as long as these stray fields are not subject to rapid and stochastic variations in time, it is possible to measure them and then apply a counter-unitary operation which undoes the unwanted precession. Equivalently, if the precession rate of each particle is known, we may let the logical encoding of the qubits become time-dependent, in order to compensate at the logical level the physical precession of the spins of the system. Notice that we focus on unitary processes since, if the error channel is given by a general master equation, the decoherence process can be characterized and mitigated, but not completely corrected.
In summary, we assume that our graph state |G is affected by a space-varying static stray field which has a different strength at each qubit of the state. We further assume that all measurements are performed at discrete time steps ∆t, so that the action of the field on each qubit is given by a SU (2) operation of the form U a = exp(−i γ a ∆t B a · S a /2). We remark that the results we will obtain in this paper will not depend on the specific implementation of the physical qubits, nor on the physical model of the error source; our methods will be applicable whenever the error process acting on the qubits is dominated by a unitary local operation acting on the qubits in the system. Nonetheless, for concreteness we will always refer to the above specified physi-cal error model. Consequently the word field can be assumed to designate any channel acting on single qubits and which is given by a fixed, but unknown, unitary mapping.
Notice that in the special case where the stray field is almost uniform and acts equally on all the qubits in the graph state, then an efficient way to achieve noise resilience is given by encoding the information in decoherence-free subspaces. The properties of these subspaces have been extensively studied elsewhere (see for a review Refs. [27,28]) and thus we will not discuss nor employ them in our work.
In conclusion, at each time step -and on each vertex of the graph -the qubits are affected by the following local operations: where a ∈ V denotes the site at which the field acts, λ a ∈ R is proportional to the strengths of the field (in our physical model λ a = γ∆t |B a |) andn a = (n a x , n a y , n a z ) (with n 2 x +n 2 y +n 2 z = 1) is the normalized vector in R 3 that gives the axis around which the rotation is performed. Observe that the relation n 2 x + n 2 y + n 2 z = 1 allows us to recover n y from (n x , n z ), assuming that its sign is always +1, since this can be enforced by the transformation (λ → −λ,n → −n), which is the identity on SO(3).
Naturally, in any realistic scenario there will be also other uncharted error sources that act on the physical system, causing some decoherence in the system. This further noise will in general render the reconstruction of the fields more difficult, but it will not qualitatively affect our results.

IV. EFFECT OF THE UNITARY ERROR ON STABILIZER MEASUREMENTS
In our setting, we are considering a QEC scheme which aims to preserve the coherence of the graph state by performing the necessary correlator measurements and reading out the corresponding results. Thus, the whole process can be schematized as follows: the graph state is prepared in the standard reference state |G ; we wait for a fixed time ∆t, during which the error channels act on the physical qubits; then we measure the correlators, which project the graph state on some state |G , which depends on the measurement outcomes of the correlators (i.e. on the measured syndromes); we finally use the syndromes to bring back |G into the reference state |G , either physically or by changing the logical encoding of the physical qubits. Implicitly, we are assuming that the state preparation, correlator readout and the final correction can be done very fast, w.r.t. the waiting time ∆t, and almost free of errors. Moreover we will assume for the moment that the stray field is the only error source present, i.e., no other error channel or decoherence process is acting on the qubits. We will lift this assumption in Sec. VIII B, where we will consider the effect of having a decoherence channel acting alongside the stray field.
As a first step in our analysis, we need to find the probability that, for each vertex a ∈ V , the measurement outcome κ a of the correlator K a is equal to 0 or 1, given that the graph state was initialized in a state |G corresponding to a graph G = (V, E), and assuming that a unitary U a has acted on each vertex. The result can be found via a straightforward computation, given in Appendix A. We define the probability difference ∆p a ∈ [−1, +1] as the difference between the probability of getting a "correct" outcome (κ a = 0) and the probability of getting an "error detected" outcome (κ a = 1) ∆p a := p(κ a = 0) − p(κ a = +1) ; then we obtain the following formula for ∆p a : where β a := cos λ a and (λ a , n a x , n b z ) are the parameters that define the action of the field on vertex a as given by Eq. (4).
Notice that from the probability differences alone it is impossible to reconstruct the correct sign of n a x and n a z , as Eq. (6) depends on these only through (n a x ) 2 , (n a z ) 2 , nor the correct value of λ a in the set {±λ a +2kπ | k ∈ Z}, since we can only access β a = cos λ a .
In principle, these ambiguities can be easily resolved. For example, it is sufficient to vary the delay time ∆t to linearly vary the parameters λ a , which can then be discriminated through the functions cos λ a . However, these methods rely on performing some extra operations, not included in the basic stabilizer code, and therefore we will not consider them in the following. In Sec. VII D we will exploit another resource, which is automatically included in MBQC, in order to reconstruct the fields. That is, we will use that some graph states are LU -equivalent even if they correspond to different graphs. LU -equivalent graph states have the same power for performing MBQC, assuming that arbitrary local measurements can be performed on the qubits; nonetheless the information that can be recovered from stabilizer measurements is very different. Hence arises the possibility of exploiting different graph states with equivalent computational power in order to better characterize the error sources.

V. PROMISE SETTINGS
The problem as specified so far is underdetermined, since we need to determine 3N parameters, i.e. the field directions and intensities, from just N parameters, i.e. from the measured probability differences. Thus we need to modify the problem to make it possible to reconstruct the stray fields. We will therefore first consider simplified settings, in which the fields fulfill a condition (a promise) that reduce the number of parameters to be estimated from 3N to N -thus permitting to estimate them through the measurement of the correlators stabilizing the graph state. We will then consider a different scenario in Sec. IX, where we will extend the set of stabilizer measurements in order to collect information on all the 3N parameters.
From now on, we will assume that the field is always globally pointing in a known and well-determined direction, but has unknown site-dependent intensity. Specifically, we investigate the cases in which the field are parallel to one of the Cartesian axes; under these assumptions, Eq. (6) simplifies considerably: the field is always aligned in theẑ direction ∆p a = β a ; 2. X-field [ n x = 1, n z = 0, β a = ? ] the field is always aligned in thex direction 3. Y-field [ n x = 0, n z = 0, β a = ? ] the field is always aligned in theŷ direction

VI. SOLVING METHODS
Here we will present a method that allows us to infer the stray fields, knowing the probability differences ∆p a , in the promise settings specified above. We will show that not in all circumstances it is possible to determine all field intensities β a , and this impossibility depends nontrivially on the connectivity structure of the graph state.
To determine the fields we need to solve a system of N polynomial equations (one for each ∆p a ) in N variables (one for each field intensity β a ) of degree at most d, where d is the maximum number of neighbours of any vertex in the graph. In general, a system of polynomial equations can be solved via computation of the Gröbner basis associated to the polynomials. There exist well-known algorithms that compute Gröbner bases [29,30], but the time complexity in the worst case scenario grows doublyexponentially in N [31]. However, with the method that we will introduce the particular polynomial equations arising from Eqs. (7-9) can be solved efficiently. Notice in particular that the case of Eq. (7) is already trivially solved, as one reads out directly the (cosine of the) intensity of the field from the measured rates. The other two cases, Eqs. (8,9), will be discussed in the following sections.

A. Complex logarithms
Consider the settings in which the stray field is parallel to either theẑ,x orŷ axis. We can write Eqs. (7-9) as b∈N s with s ∈ {x, y, z}, and setting N z a = {a} for the Z-field case, N x a = N a for the X-field case, and N y a = N a ∪ {a} for the Y -field case. They are polynomial equations in the variables β a , in which each polynomial has only two terms, i.e., the monomials appearing on the left hand and on the right hand side. To solve these we take the logarithms of both sides of the equations.
In general the cosines of the fields (β a ) and the probability differences (∆p a ) have support in the interval [−1, 1]. Thus we have to use the logarithm as a complexvalued function; more precisely, we define the logarithm as taking values on the Riemann surface C/ 2πi : where ln(|z|) is the real logarithm (mapping R + in R), C * = C \ {0} and the Riemann surface C/ 2πi is the complex plane modulo identification of numbers that differ by an integer multiple of 2πi; we also denote by = 2πi the equality on this Riemann surface. If we restrict z to be a real number, then (i arg(z)) is in the set {0, iπ}. With this definition the logarithm is a single-valued function and the multiplication rule applies: Notice that on the Riemann surface C/ 2πi the division is not uniquely defined. In particular, when performing a division by m ∈ Z \ {0} we get |m| different valid results:

B. Linear equations
Thus, after taking the logarithm of both sides of Eqs. (7-9), we obtain the following linear equations: The adjacency matrix A of a graph G is a symmetric matrix defined as We also define A s for s ∈ {x, y, z} as where 1 denotes the N × N identity matrix.

C. Multiple solutions
The problem thus has been reduced to solving a linear system of the form A s v = 2πi w, with v := (ln β 1 , . . . , ln β N ) T and w := (ln ∆p 1 , . . . , ln ∆p N ) T . From here, provided that the matrix A s is non-singular, the (logarithm of the cosines of the) fields can be reconstructed via linear algebra. However, in this context even when A s is non-singular it is not formally correct to write the solution as w = (A s ) −1 v, since the system A s v = 2πi w can have multiple solutions. This happens because the entries of A s are taken as elements of C/ 2πi , thus divisions by integer number have to be performed according to Eq. (13). Therefore, in order recover all the valid solutions to the linear system (14) we need to employ Eq. (13) whenever performing a division during the solving algorithm.
We adopt a modified version of Gauss elimination algorithm which allows us to recover these multiple solutions.
Since the entries of A s are all integers, it is possible to bring A s into the form of an upper triangular matrix Q, applying only a sequence of the following elementary Gauss operations on A s : These operations can in fact be used to compute the Hermite normal form [32] of A s . That is, it is always possible to write A s in the form: where P and Q are matrices with integer entries, P unimodular (| det P | = 1) and Q upper-triangular; for more details see, e.g., Ref. [33]. Notice that, P being unimodular, P −1 is also unimodular and has integer entries; moreover, P −1 is obtained via multiplication of the matrices that represent the elementary Gauss operations. Hence the solutions of A s v = 2πi w are in one-to-one correspondence with the solutions of Thus Eqs. (14) have been brought into the form Q v = 2πi w , with w := P −1 w, and Q upper-triangular.
If A s is non-singular, then all the diagonal entries (Q 11 , . . . , Q N N ) of Q are different from 0. Therefore the solutions can be recovered finding the |Q N N | solutions (see Eq. (13)) to the equation Q N N v N = 2πi w N , then substituting the obtained values for v N in Eq. (14), then solving for v N −1 and so forth. That is, having computed the possible values for (v l+1 , . . . , v N ) we obtain v l by solving the equation In total, the number of complex solutions is given by An equivalent way of looking at this operation, which will become useful in the error analysis in Sec. VIII, is to consider the vector w as taking | det A s | different values, and then computing the corresponding solutions for v with standard linear algebra over the complex numbers. Explicitly, we need to solve | det A s | different linear equations in the form Qv = w c1...c N over C, in which the known terms are given by whereê l is the unit vector with a 1 in position l, and In conclusion, we face two possible scenarios: 1. The matrix A s is singular; the system (14) is either inconsistent or underdetermined.
2. The matrix A s is non-singular; the number of complex solutions is given by | det A s |.
Finally, once the solutions v to the system (14) have been computed, the corresponding values of the parameters β a can be obtained via element-wise exponentiation of v.

D. Constraints on the solutions
Now consider w = (ln ∆p 1 , . . . , ln ∆p N ) T ; the components of w are logarithms of real numbers, thus the imaginary part must be an integer multiple of π (i.e. the imaginary component is either 0 or iπ on C/ 2πi ). After the Gauss elimination algorithm, we obtain another vector w = P −1 w whose entries are integer linear combinations of those of w; thus, the entries of w still have an imaginary part which is an integer multiple of π.
Take v l as a complex solution of an equation of the form if Q ll is odd, then among these solutions only one has an imaginary part which is an integer multiple of π; if Q ll is even, then among these solutions only two have imaginary parts which are integer multiples of π. Repeating this argument for all entries of v we get the number of allowed solutions is 2 µ , where µ is the number of even diagonal entries in Q.
The requirement that all the probabilities and the (cosines of the) field strengths must be real thus restricts the allowed solutions to just a small subset of all possible complex solutions. There is still another requirement that should be enforced: the parameters β a and ∆p a have to lie in the range [−1, +1]; this condition further restricts the number of legitimate solutions. Notice that, when we enforce this requirement, the number of acceptable solutions cannot be expressed in terms of properties of the adjacency matrix alone (A s ), since this number does depend on the data actually collected (w).

VII. SOLUTION FOR PARTICULAR CLASSES OF GRAPH STATES
In this section, we will explicitly determine whether a graph state allows a complete reconstruction of the fields acting at each site -assuming that it is parallel to either thex,ŷ orẑ direction -for some relevant classes of graphs. We will consider at the beginning graphs which represent one-dimensional lattices -that is, linear chains. In these graphs each vertex is connected to at most two neighbours. Afterwards, we will show how to extend those results to two-dimensional (and higher dimensional) square lattices. Finally, we will give an application for states used in quantum repeaters and measurement-based approaches to quantum error correction.  (20) in which we introduce the notation . Their determinants can be given as closedform expressions [34]: Therefore, given the solution methods presented in Sec. VI, it is evident that the system is underdetermined when m is even in the case of X-fields, and when m ≡ 2 (mod 3) for Y -fields. In all other cases there is exactly one (complex) solution to the equations.
For a closed linear chain with m vertices the adjacency matrix is given by: and thus A close x and A close y are circulant matrices -where a circulant matrix C is defined by the relation . Also in this case it is possible to give closed-form expressions for their determinants [35], but we have to take into account that m = 1 and m = 2 are special cases:  Table I Open-end linear chain Closed linear chain Adjacency matrix: Adjacency matrix: tridiagonal Toeplitz matrix tribanded circulant matrix The adjacency matrix is not invertible when: X-field: X-field: in which we have introduced the function cos m (θ) to deal with the special cases arising when m = 1 or m = 2: The results above will turn out useful in the next section.

B. Square lattices
Here we describe the method for solving the problem for two dimensional square lattices. The same constructions can be applied straightforwardly in an arbitrary number of dimensions, so in the end we will present the formulas also for the general multidimensional case.
As a first step, we give the expression for the adjacency matrix of a graph representing a two-dimensional square lattice, having m 1 rows in one direction and m 2 columns in the other direction, thus having in total N = m 1 · m 2 vertices. We number the vertices according to a scan-line pattern, i.e. we assign to the vertex in position (j 1 , j 2 ), with j 1 ∈ [m 1 ] a row index and j 2 ∈ [m 2 ] a column index, the value j 2 + (j 1 − 1)m 2 . With this numbering convention, the adjacency matrices A of planar lattices take the form of block-Toeplitz matrices. That is, there are m 1 × m 1 blocks, each one consisting of m 2 × m 2 elements, and the blocks which lie on the same diagonal are all equal; the blocks themselves are Toeplitz matrices.
For a square lattice with open boundaries in both directions, i.e. for a planar graph, the adjacency matrix takes the following tridiagonal block-matrix form: where here a, b are block indices, the blocks are of size m 2 × m 2 , and A open m2 is given by Eq. (20). Now we consider another family of graphs, which are obtained from the above planar square lattice by connecting the vertices of the m 2 -th column to the corresponding vertices in the first column; so we obtain a graph which is topologically a cylinder. Then the adjacency matrix takes the following block-matrix form: where A close m2 is given by Eq. (23). Finally, we consider the family of graphs obtained from the cylindrical graph by connecting the vertices of the m 1 -th row to the corresponding vertices of the first row; so we obtain a graph which is topologically a torus. The adjacency matrix takes the following block-matrix form: with the same conventions as above.
As usual, we have A x = A and A y = A + 1 for all the cases above, where 1 is the m 1 m 2 × m 1 m 2 identity matrix. Now we will show how to obtain the eigenvalues of the above matrices in closed form. Once the eigenvalues are given, it is straightforward to tell whether all information about the fields can be recovered or not, by searching for eigenvalues equal to zero. If these can be found, it means that the adjacency matrix is singular and not all the information can be recovered. The number of zero eigenvalues is the rank defect of the adjacency matrix, and is a measure of the amount of information that one cannot reconstruct.
To solve our problem, we study tridiagonal Toeplitz matrices and circulant matrices with three bands, of size  m × m. A tridiagonal Toeplitz matrix T has the form: and its eigenvalues are given by the expression [34]: By standard spectral methods, it is straightforward to show that this expression is meaningful also when s, t, u are pairwise commuting diagonalisable matrices. Under these assumption, the matrices are simultaneously diagonalisable, and thus it is sufficient to define the necessary operations as acting over the eigenvalues of these matrices as given in the common eigenvector basis. Similarly, a circulant matrix C has the form: for an arbitrary set of coefficients (c 0 , . . . , c m−1 ); to specialise to the case in which there are only three bands, it is sufficient to set c j ≡ 0 for j / ∈ {−1, 0, +1}. The eigenvalues of C are then given by [35]: where ω k := exp(2πik/m) are the roots of unity. The result can be specialised to the case in which s ≡ c 0 , t ≡ c −1 = c 1 , with t = t † = t * and c j = 0 for j / ∈ {−1, 0, +1}, giving the formula: As we argued before, again this formula is valid also when s, t are commuting diagonalisable matrices. A summary of these properties can be found in Table II.
At this point we are ready to compute all the eigenvalues of the adjacency matrices A planar , A cylinder and A torus . The matrix A planar , as given by Eq. (31), is a block tridiagonal Toeplitz matrix, thus its (block) eigenvalues are given by Eq. (35), with the identifications s = A open m2 and t = u = 1 m2 : These "eigenvalues" Λ planar k1 are by themselves m 2 × m 2 tridiagonal Toeplitz matrices, whose eigenvalues can be computed using again formula (35), this time with the identification s = 2 cos(πk 1 /(m 1 + 1)) and t = u = 1. Finally, a (real) eigenvalue of Λ planar An analogous computation can be performed for A cylinder and A torus , resulting in: Naturally, the eigenvalues of A y are equal to the eigenvalues A x plus one, for any adjacency matrix A. A list of the square lattice graph states of size up to 20 × 20 for which the error channel can be reconstructed is given in Fig. 1.

C. Cubic lattices in arbitrary dimension
We conclude this section remarking that the procedure we have outlined carries through any number of dimensions. Therefore we can easily generalize the above results to a general cubic lattice in d = d 1 + d 2 dimensions, assuming closed boundary conditions in d 1 dimensions and open boundary conditions in d 2 dimensions. The eigenvalues of the adjacency matrix of these classes of cubic lattices are given by: where m r is the size of the lattice along dimension r.

D. Application: quantum repeaters
In this section we provide an example of the application of the methods we have developed. We show that in many cases these techniques can be applied to states that are useful resources for measurement-based quantum repeaters and measurement-based QEC.
A quantum repeater is a protocol that allows one to efficiently transmit quantum information over long distances in presence of noise that grows multiplicatively with the distance -which happens, e.g. when transmitting photons along an optical fiber. To achieve this, it is necessary to add several intermediate stations along the transmission line; in these stations entanglement purification and entanglement swapping are applied to an ensemble of noisy Bell pairs in order to establish a single high-fidelity Bell pair between the distant parties [22,23].
In the measurement-based approach to QEC a logical qubit is encoded using a small size graph state as a resource. Then the information encoded in the logical qubit can be teleported to a freshly prepared low-error qubit entangling the graph state with this ancilla qubit and then performing local measurement on all the qubits of the original graph state. Alternatively, both entanglement purification and entanglement swapping can be performed at the same time exploiting again a local measurement pattern on specific graph states [36].
In Refs. [23,24] some examples of graph states are given, up to local unitary (LU ) corrections, that can be used for measurement-based quantum repeaters. Here we will not consider the LU -corrections, since these can always be included in the final read-out measurements; namely, we can effectively perform a LU -correction U a by measuring the (local) observable U † a O a U a in place of the uncorrected observable O a .
These graph-state resources include open-ended and closed linear chains of length 3 and 5; GHZ states with 4 and 6 qubits; and other, more irregular graph states. Since the linear chain cases have been already analysed in Sec. VII A, we focus on the other cases. Table III. Determinants of the adjacency matrices of the fullyconnected graph and of the star graph.
We recall that a GHZ state with N qubits can be represented as a graph state, up to LU -corrections, in at least two ways:  Table III. Another application of graph states is to implement measurement-based versions of QEC. For example the Steane 5-qubit QEC code can be implemented using a graph state representing a closed linear chain with 5 vertices together with a sixth vertex fully connected to all the other vertices. It can be checked that for this graph both the X-field and Y -field cases can be solved. Similarly, more complex graphs can be analysed case by case.
Finally, we remark that there is a very simple test to show in some cases that the X-field setting cannot be reconstructed. This happens whenever any two vertices in the graph have the same set of neighbours. In this circumstance two rows in the adjacency matrix are equal, and thus the matrix itself is singular.

VIII. ERROR PROPAGATION
In this section we address the problem of how uncertainties in the probability differences ∆p a and in the field alignments propagate through the reconstruction method we have given.

A. Effect of uncertainty in the measured probabilities
First we consider the effect of the uncertainty arising from the fact that the correlator measurements are performed only a finite amount of times; thus the probability differences ∆p a are known only up to a finite precision. For the moment however, we still assume that the stray fields are perfectly aligned to one of the axes (x,ŷ orẑ).
Thus, each correlator K a at each vertex a is measured M times, yielding a string of results (κ From these measurements, we can extract an error rate difference ∆R a ∈ [−1, +1], defined as the difference between the frequency of κ a = 0 outcomes and κ a = 1 outcomes: The random variables ∆R a are sampled according to a shifted and rescaled version of the binomial distribution, and therefore they have mean and variance equal to: (45) That is, ∆R a are unbiased estimators of ∆p a . In the following, we will take the limit M 1 and thus assume that R a is approximately distributed according to a normal probability distribution function with very small variance (a delta-like probability distribution function).
We then move to using estimators in logarithmic variables; i.e., we consider the estimator of ln(∆p a ) as given by ω a ≡ w a + δw a := ln(∆R a ), where ω → w when M → ∞, with the definition of w given in Sec. VI C. The logarithm again takes values in C/ 2πi . When the distribution of R a is delta-like and E(∆R a ) = 0 (or, more correctly, E(∆R a ) Var(∆R a )) it is legitimate to approximate the expectation value and variance of w a = ln(R a ) with the following formulas: Then, through the reconstruction procedure given in Sec. VI, it is possible to give estimates of ln(β a ) in terms of estimates of ln(∆p b ). Namely, remembering that there are | det A s | solutions to an equation in the form A s v = 2πi ω (when det(A s ) = 0), we set as estimator of ln(β a ) the value v c a : s ω c denotes the c-th complex solution of the equation A s v = 2πi ω. We recall that this can be achieved by inverting the matrix A s in C and then letting ω assume | det A s | different values ω c , as shown in Eq. (19) (the index c is in a one-to-one correspondence with the indices (c 1 , . . . , c N )), which then get mapped into | det A s | different solutions. Therefore, the error propagation from ω to v is well-behaved if A s is non-singular, and it is legitimate to use the inverse matrix A −1 s , as long as one remembers to apply it to all the vectors ω c .
Thus it is possible to give the mean and variance of the probability distribution function of v, in a neighbourhood of each solution, in terms of the expectation value µ := E[ω] and covariance matrix Σ a,b := Cov(ω a , ω b ) = Cov(ω c a , ω c b ): in which the relation is exact, since the equations are linear. Notice that also the covariance matrix of v c is independent of c. From now on we will omit the index c, since it will always be implied when giving a reconstruction of v depending on ω.
A characterization of error propagation can be given as follows. We have that ω = w + δw represents the measured (perturbed) value of w, with δw w , where · is any chosen vector norm; then the reconstructed solution of the system A s x = w + δw is given by x ≡ v + δv. Then the following bound, arising from standard perturbation analysis, holds: where A s is the operator norm induced by the chosen vector norm and κ(A s ) := A s A −1 s denotes the condition number of matrix A s (which also depends upon the adopted vector norm).
Another measure of the magnitude of the errors that can be adopted is the uncertainty volume in logarithmic variables: where Σ ab := Cov(X a , X b ) is the covariance matrix of set of random variables X. Hence, taking the square root of the determinant of (49), we get Vol(v c ) = det(Σ)/| det A s |. That is, the uncertainty volume around each solution always shrinks by a factor | det A s | Notice, again, that the number of complex solutions is equal to |det A s |, and that the above estimate gives the size of the error region around each valid solution. Thus there is a global conservation of uncertainty volume in our reconstruction method. However, not all complex solutions are physical; namely, as discussed in Sec. VI, we have to restrict ourselves to the case in which the solutions correspond to real values for the stray field intensities. Restricting to these physically relevant cases, the total uncertainty volume indeed can decrease, and it can never increase.

B. Full error analysis
We will now deal with the full error analysis, in which the stray fields satisfy the constraints only approximately. That is, we analyse the case in which the fields are approximately aligned to one of the three Cartesian axeŝ x,ŷ orẑ. In all three situations we have to solve an approximate linear systems of the form where f (x) is a vanishing function for → 0, with a global parameter which bounds the maximum misalignment of the rotation axis from one of the Cartesian axes.
We call x a solution of this perturbed system, and with v a solution of the original unperturbed system. In Appendix B we show that, under some regularity assumptions, the magnitude of the perturbation in the solution is bounded by a continuous function of δw and of the rotation axis misalignments. In principle, one should consider that in any realistic implementation there will be some incoherent noise acting alongside the local unitary rotations induced by the stray fields. Equivalently, we might consider the case in which the stray fields λ a are not fixed vectors, but they vary in time. If these temporal variations are fast and stochastic, then it can be imagined that λ a are random variables drawn according to some probability distribution.
The simplest model, but also the worst case scenario, for this source of incoherent noise is given by the uniformly depolarizing channel where ρ is a single qubit density matrix and 0 ≤ q ≤ 1. It is not necessary to specify in which order the depolarizing channel and the local unitary rotations are applied, since they commute. From Eq. (5) one immediately obtains the difference in probability between the measurement outcomes when both local unitary error sources and depolarizing noise Φ q act on the system: where ∆p a is given by Eq. (6). This noise affects the results, favouring smaller values of β a ; for example, in the Z-field case all β a are just rescaled by a factor (1−q). We point out that, for the X or Y -field cases, for some of the qubits in the graph state the reconstruction algorithm may exhibit a natural resilience to the inaccuracy introduced by the uniformly depolarizing noise. In fact, from the expression in Eq. (55) it can be seen that, to account for the action of the depolarizing noise, the differences in probabilities have to be multiplied by (1−q) for all qubits. Thus the vector log(1 − q) 1 := log(1 − q) (1, 1, . . . , 1) T has to be added to the known term in the linear equation (14), whose entries are the logarithm of the differences in probabilities. By linearity, the effect of the depolarizing noise on w (the vector of the estimators in logarithmic variables) is simply given by the addition of log(1−q)×A −1 s 1. Some of the entries of the vector A −1 s 1 might be zero, implying that the reconstruction of the fields for those entries is unaffected by the depolarizing noise.
We have tested through numerical simulation the performance of our reconstruction method when all the three noise factors considered above are present, i.e., the effect of the finiteness of the number of measurements of the correlators, the misalignment of the field from the promised direction, and the independent action of a depolarizing channel on the qubits. These simulations provide evidence that, indeed, when all the noise factors are small (and the numbers of measurements of each correlator is large) then the accuracy in the reconstruction improves and, in principle, the inaccuracy can be brought arbitrarily close to zero.
We define the reconstruction error E a on any vertex a ∈ V as ]. An excerpt of the reconstruction errors we have obtained, in the case of an open-ended chain with ten vertices, is given in Fig. 2.
A general trend that can be inferred from the numerical data is that -since the correlation among the parameters is not taken into account -the reconstruction error is minimal in the Z-field case, is larger in the X-field case, and is still larger in the Y -field case. This feature is expected, since the degree of the polynomial equations to be solved is minimal in the Z-field case, larger in the X-field case and still larger in the Y -field case. This implies that the cross-correlations between the estimators of the fields are maximal in the Y -field case and they are minimal in the Z-field case.
Another feature that can be inferred from the data is that the precision of the estimation of the fields varies across the vertices of the graph and depends on the global connectivity of the graph itself. Notice that to produce the data shown in Fig. 2 we have assumed that the model of the disturbance is the same for all ten vertices in the chain. Thus for symmetry reasons we expect that, averaging over all possible intensities of the stray field, for a open-ended chain with ten qubits the inaccuracy in the reconstruction of the fields on vertex a and 11 − a are the same.
Also, notice that the reconstruction error decreases approximately as 1/ √ M for small values of M , while for larger values of M the reconstruction error approaches a non-zero asymptotic value, which is due to residual uncertainties introduced by the depolarizing noise and stray field misalignment.
Finally, we remark that at some of the vertices the reconstruction of the fields shows resilience to the depolarizing noise, for the reason given earlier. Explicitly, this happens for qubits on the vertices 3, 4, 7, 8 for the X-field case and for qubits on the vertices 2, 3, 5, 6, 8, 9 for the Yfield case, as can be easily checked computing A −1 s 1 for these two cases. Consequently, it can be seen in Fig. 2 that these qubits are almost insensitive to the magnitude q of the depolarizing noise, provided that q 0.8. When q → 1, then the effect of the finiteness of the number of measurement M renders unreliable the reconstruction of the fields on those qubits.

IX. EXCURSUS: EXTENSION OF THE FAMILY OF GRAPH STATES
A completely different approach to estimating the stray fields consists in adopting different sets of stabilizer measurements in order to collect more information. In this way we do not need to assume that the fields are always aligned in a known direction. Using a different set of stabilizer measurements will produce a different graph state; but this is not an issue, as long as the functionality of this family of graph states is the same of that of the reference graph state.
Explicitly, we consider the following scenario. In the defining equation (2), the correlators are defined in terms of the product of Pauli-X and Pauli-Z matrices. Actually, many other correlators will yield graph states which have the same computational capabilities as the reference graph state. For example, we could as well use correlators in the form Y a b∈Na X b (or in the form Z a b∈Na Y b ). To see that these stabilizers define a graph state which is computationally equivalent to the original one it is enough to observe that these are obtained by a redefinition of the Pauli matrices, (X, Y, Z) → (Y, Z, X) (or (X, Y, Z) → (Z, X, Y )), which is a transformation that preserves the commutation relations between the matrices -and thus all the associated algebraic properties. Correspondingly, the local measurements necessary to perform MBQC have to be modified according to the new Pauli basis.
These considerations apply to any transformation (X, Y, Z) → (r · S,ŝ · S,t · S), in which B ≡ (r,ŝ,t) is any right-oriented orthonormal basis of R 3 . We call these logical Pauli matrices. Correspondingly, Eq. (6) becomes: Thus, using three different graph states, associated to three different logical Pauli bases it is in principle possible to reconstruct all the parameters of the stray fields. Notice that in general there will be many physically acceptable solutions to a system of 3N polynomial equations with 3N unknowns. But then using ≥ 4 logical Pauli bases, from the corresponding graph states one obtains an overdetermined system of N equations with 3N unknowns, which almost certainly admits only one solution; moreover, having more independent equations usually provides a reconstruction which is more robust in presence of small uncertainties. Explicitly, we look for an approximate solution to the equations through numerical Here ∆p a ({β v ,n v }; B) is computed from Eq. (57) and it is a function of the field intensities β v and directionsn v at the vertices in N a = {a} ∪ N a , and of the choice of the logical Pauli basis B (which can be identified with an orthonormal basis in R 3 ). With ∆p (mes) a (B) we denote the value for the probability difference estimated from syndrome measurements for vertex a using the logical Pauli basis B. An example of the result that can be obtained via this method is reported in Fig. 3.

X. CONCLUSIONS
In this paper we have shown how to use the information collected from a set of stabilizer measurements in order to extract information about a unitary error process acting on a graph state. Since the stabilizer measurements are non-local operations, each of them depends on the errors that act on several different sites of the graph state. Therefore the measurement outcomes of these stabilizers are correlated.
We have studied in depth the case in which the local unitary errors are given by rotations aligned to one of the Cartesian axes. This condition is needed to make the problem solvable, i.e., we match the number of parameters to be reconstructed to the number of measured parameters; moreover this restriction allows us to find efficient reconstruction methods.
A surprising result arising from the fact that the stabilizer measurement outcomes are correlated is that in some circumstances it is not possible to reconstruct all the parameters of the unitary error -even under the assumption that we can measure the stabilizers an arbitrary number of times. The possibility of completely reconstructing the error turns out to depend non-trivially on the connectivity structure of the graph state. Specifically, in order to determine the intensity of local unitary rotations around thex axis, the adjacency matrix A associated to a graph state G has to be a full-rank matrix, whereas to reconstruct rotations around theŷ axis the matrix A + 1 has to be full-rank, while rotations around theẑ axis can always be reconstructed.
We have ascertained for some classes of graph states whether it is possible to reconstruct all the errors. Hence we have studied the properties of the adjacency matrices of linear chains (with open and closed boundary conditions), square lattices in arbitrary number of dimensions, and GHZ states (these represented either as fullyconnected graphs or as star graphs).
Finally we have shown that the performed reconstruc-tion is resistant to various kinds of imperfections. In particular we have investigated the effect of a finite number of measurements of stabilizers, of the misalignment of the rotation axis and of the action of a depolarizing noise alongside the action of the unitary channel. In each of these cases, as long as the imperfections are small, the precision in the determination of the unitary channel is only mildly affected.

ACKNOWLEDGMENTS
We acknowledge support from the Austrian Science Fund (FWF) through the SFB FoQuS: F4012, and the Templeton World Charity Foundation grant TWCF0078/AB46.

APPENDICES Appendix A: Derivation of the error probability
In this appendix we derive Eq. (6). This equation gives the difference between the probability of getting κ a = 0 and κ a = 1 when measuring the correlator K a = X a b∈Na Z b on a graph state which is subject to local noise of the form U f ield = a∈V U a = a∈V e −iλana·Sa/2 . We introduce the notation | G := U f ield |G ; observe that K a has spectrum {+1, −1}, therefore there is a unitary mapping V a such that K a = P 0 a − P 1 a , where P j a with j ∈ {0, 1} is the projector on the (−1) j -eigenspace of K a . Hence we get that the difference in the probabilities is given by the expectation value of K a on | G : Now we can compute: where N a = N a ∪ {a}, and in which we have used that K a |G = |G and that U e A U † = e U AU † when U is a unitary operator. Hence: G| K a | G = = G| U † noise K a U noise K a |G = G| e i λa 2n a ·Sa e −i λa Actually, we only need to compute the part proportional to the identity in the operator in the last line of Eq. (A3), since G|X|G = G|Y |G = G|XZ|G = 0 for any non-trivial graph G, i.e. for graphs in which there is no vertex which is completely disconnected from all other vertices (no isolated vertices). Indeed, for the Pauli Z operation on the vertex a we have: in which we have used twice K a |G = |G and the anticommutation relation between Z a and K a . An analogous procedure shows that G|Y a |G = 0, since also Y a and K a anti-commute. Then, for the Pauli X operation on vertex a we have (provided that a has at least one neighbour): Finally, we can evaluate the part proportional to the identity in Eq. (A3) applying the exponentiation formula for Pauli matrices: e iλ(n·S) = 1 cos(λ) + i (n · S) sin(λ) (A6) together with the identity: for α, β ∈ {x, y, z}, obtaining: G| K a | G = = G| 1 a cos 2 λ a 2 + (n a x ) 2 − (n a y ) 2 − (n a z ) 2 sin 2 λ a 2 b∈Na This is exactly the expression given in Eq. (6), once one writes β a := cos λ a .