Ising model formulation for highly accurate topological color codes decoding

Quantum error correction is an essential ingredient for reliable quantum computation for theoretically provable quantum speedup. Topological color codes, one of the quantum error correction codes, have an advantage against the surface codes in that all Clifford gates can be implemented transversely. However, the hardness of decoding makes the color codes not suitable as the best candidate for experimentally feasible implementation of quantum error correction. Here we propose an Ising model formulation that enables highly accurate decoding of the color codes. In this formulation, we map stabilizer operators to classical spin variables to represent an error satisfying the syndrome. Then we construct an Ising Hamiltonian that counts the number of errors and formulate the decoding problem as an energy minimization problem of an Ising Hamiltonian, which is solved by simulated annealing. In numerical simulations on the (4.8.8) lattice, we find an error threshold of 10.36(5)% for bit-flip noise model, 18.47(5)% for depolarizing noise model, and 2.90(4)% for phenomenological noise model (bit-flip error is located on each of data and measurement qubits), all of which are higher than the thresholds of existing efficient decoding algorithms. Furthermore, we verify that the achieved logical error rates are almost optimal in the sense that they are almost the same as those obtained by exact optimizations by CPLEX with smaller decoding time in many cases. Since the decoding process has been a bottleneck for performance analysis, the proposed decoding method is useful for further exploration of the possibility of the topological color codes.


I. INTRODUCTION
In recent years, there has been remarkable progress in the development of quantum computers, with the realization of quantum computers of tens to hundreds of qubits.While entering the realm of hard-to-simulate areas for classical computers [1][2][3][4], demonstrating quantum advantage in meaningful problems remains challenging.This is due to the current high noise level of quantum computers, which hinders the execution of complex quantum algorithms, such as Shor's factorization [5], linear system solver [6], and quantum phase estimation [7].To solve this problem, it is essential to realize a faulttolerant quantum computer, which can perform computation while protecting quantum information from errors through quantum error correction (QEC) [8].
The surface codes [9], one of the topological codes [10], are thought to be one of the most promising approaches for the experimental implementation of QEC due to their simple structure, making physical implementation easier, with a relatively high error threshold.Currently, experimental demonstrations of QEC on a single logical qubit level of the surface code are ongoing [11][12][13].Nevertheless, surface codes have a drawback as QEC codes.Specifically, only X, Z, and CNOT gates can be implemented transversally among the Clifford gates [14].Con-sequently, a special treatment is required to implement H and S gate fault-tolerantly [15][16][17], resulting in additional overhead.
The QEC codes that resolve this drawback are the color codes [18,19].A color code has an advantage against the surface codes in that all Clifford gates can be implemented transversally, due to its high symmetry of the stabilizer operators [20].Despite this advantage, color codes are currently not the mainstream of experimental implementations of QEC due to the disadvantages, the difficulty of decoding and its low error threshold under the circuit-level noise model.For the surface codes, there is a known decoding algorithm, minimumweight perfect matching algorithm (MWPM) [21], which can be executed in polynomial time and with high accuracy.On the other hand, in the case of color codes, such a good decoding algorithm is still missing.The optimal threshold of a color code on the (4.8.8) lattice is estimated to be 10.9% [22] under the bit-flip noise model by a Monte Carlo simulation of the corresponding statistical mechanical model.However, the performances of efficient decoders, such as the renormalization group decoder [23], the restriction decoder with MWPM [24], and the restriction decoder with union-find [24], resulting in a threshold of 8.7%, 10.2%, and 9.8% (analytical estimate), respectively.Although the integer program decoder [20] exhibits good accuracy with a threshold of 10.6%, its decoding process incurs exponential computational time.Regarding the depolarizing noise model, the optimal threshold is estimated to be 18.9% [25].Nonetheless, the neural-network decoder [26] has a comparatively low accuracy of 17.5%.Moreover, for the phenomenolog-ical noise model, which takes into account errors that occur during syndrome measurement, prior studies have shown that the graph matching decoder has a poor accuracy of 2.08% [27].Also, the integer program decoder [20] has better accuracy at 3.05%, but requires exponential time for decoding.Although it has been demonstrated that tensor network decoders can achieve a relatively high threshold for the bit-flip and depolarizing noise model on the (6.6.6)color codes [28,29], it is unclear how high of a threshold value can be obtained for the phenomenological noise model.The threshold for (6.6.6)color codes under phenomenological noise was estimated to be 4.8(3)% [30,31] by mapping the problem onto a statistical-mechanical three-dimensional disordered Ising lattice gauge theory, but it is not the value obtained by decoders.The difficulty in decoding implies a challenge in evaluating its performance, thus obstructing the exploration of methods to overcome the issue of a low error threshold under the circuit-level noise model, which is another drawback of the color codes.
In this paper, we propose a Monte Carlo type decoder for color codes.Monte Carlo type methods have been applied for decoding topological codes so far [32,33].However, in these approaches, an Ising spin variable is assigned for each error, and then syndrome constraints are imposed as penalty terms, leading to relatively low accuracy of the error correction.In the proposed scheme, we instead assign an Ising spin variable for each stabilizer operator.Then we construct an Ising Hamiltonian that counts the number of errors via the multi-body Ising interactions and formulates a decoding algorithm by minimizing the energy of the Ising Hamiltonian, which we solve by using simulated annealing (SA) [34].Specifically, when syndrome measurement is perfect, assigning an Ising spin variable for each stabilizer generator allows us to represent all error patterns.However, when errors are also located on measurement qubits, it becomes not so straightforward.In this situation, we introduce time-like stabilizer generators and assign an Ising spin variable for each of these.As a result, we can represent all patterns of both data errors and measurement errors and then construct an Ising Hamiltonian that counts the number of errors appropriately.The syndrome constraints are imposed on the initial spin configuration and hence satisfied throughout the optimization process.
We perform extensive numerical experiments to examine the performance of the proposed decoder using an open source SA solver and compare its performance with CPLEX [35], which solves the same problems exactly via integer programming.We find the error threshold of 10.36(5)% for the bit-flip noise model, which is almost equivalent to that obtained by the integer program decoder with the highest accuracy in previous studies [20].
For the depolarizing noise model, we find 18.47(5)%, which is higher than the neural-network decoder [26].For the phenomenological noise model, where only bitflip errors occur as data errors and errors are also introduced on the measured syndrome, we find 2.90(4)%, which is again almost equivalent to that of the integer program decoder.While numerical experiments for the circuit-level noise model were not conducted, this decoder can also be applied straightforwardly for the circuit-level noise model.
In terms of decoding time, the proposed method is faster than CPLEX in the cases of the bit-flip noise model and depolarizing noise model.On the other hand, for the phenomenological noise model, the proposed method is faster than CPLEX only for a small code distance.However, we should note that SA employed here can be further optimized for the Ising models formulated here.Furthermore, SA holds the advantage of enabling parallel computation.Therefore, we believe that optimizing and parallelizing SA for this task could lead to decoding with high accuracy while faster than other decoders that achieve comparable accuracy for a wider class of QEC code.
The rest of the paper is organized as follows.In Sec II, we first introduce the definition and characteristics of color codes.Then, we discuss the problems of one of the existing decoders that is closely related to our method, which formulates the decoding problem as a combinatorial optimization problem by assigning variables to each error.In Sec III, we provide a detailed description of our decoding formulation for two types of code capacity noise models: the bit-flip noise model and the depolarizing noise model.In Sec IV, we describe the decoding formulation for the phenomenological noise model.Firstly, we introduce an existing decoding method that formulates the decoding problem as a combinatorial optimization problem by assigning variables to each error, similar to Sec II for the bit-flip noise model.Next, we provide a detailed explanation of our decoding formulation in this noise model.In Sec V, we provide details of the numerical experiments performed and show the logical error rates achieved by our method through Monte Carlo simulations for each noise model.We also report the decoding time measured by the computer and provide a comparison and discussion with CPLEX.Then, Sec VI is devoted to a conclusion.

II. PRELIMINARY A. Color codes
Color codes are topological codes defined on a lattice where each vertex is incident to three edges and adjacent faces are colored with three different colors [18], as shown in Fig. 1.Also, the edges of the lattice are colored with three different colors.Qubits are placed on each vertex of the lattice.For each face, the Z-and X-stabilizer generators are defined as the tensor product of Pauli Z and X operators acting on each qubit located at the vertices included in the face: The code state is defined as the simultaneous +1 eigenspace of stabilizer generators.Specifically, we here use the color codes defined on the (4.8.8) lattice with open boundaries, depicted in Fig. 1.The (4.8.8) lattice is a semi-regular lattice where a square and two octagons meet at every vertex.Similar notations apply to other lattices as well, for example, the lattice denoted as (6.6.6)lattice is a regular lattice where three hexagons meet at every vertex.In terms of transversal gates, a method is proposed to transversally implement all Clifford gates on any 2D color code [36].However, the (4.8.8) color codes are the only 2D color codes that can implement all Clifford gates transversally in a simple way [20].

B. Decoding and combinatorial optimization
For simplicity, let us first consider decoding under the bit-flip noise model.Bit-flip noise model is a noise model in which the Pauli X operator acts on each physical qubit with probability p.This noise model is written by the following map: The problem of finding the error configuration that minimizes the number of errors satisfying the syndrome is formulated as the integer programming problem: i∈f where x i ∈ {0, 1} denotes a binary variable representing the error on the i-th qubit; x i = 1 (= 0) represents the presence (absence) of error.s f ∈ {0, 1} is the syndrome value of the Z-stabilizer generator defined on each face f .The minimum distance decoding is executed by exactly solving the integer programming problem represented by Eqs. ( 4) and ( 5) [20].While it has not been known whether or not the problem of decoding the color codes is NP-hard, in general, decoding problem of QEC codes formulated as an integer programming problem is NP-hard [37], making it unrealistic when the system size is large.
In this paper, we propose a decoding algorithm using SA, which is a heuristic approach for solving combinatorial optimization problems.For example, the integer programming problem represented by Eqs. ( 4) and ( 5) is rewritten as an unconstrained optimization problem by introducing a penalty term [32,33]: where λ is a hyperparameter referred to as the penalty coefficient and its value must be manually set beforehand.
The difficulty in adjusting this penalty coefficient results in poor accuracy.To achieve high accuracy, we avoid this formulation and map the decoding problem to an unconstrained optimization problem, where no penalty term for the syndrome constraints is introduced.

A. Mapping decoding problem to an Ising model
Here, we adopt a different approach to formulate the decoding problem under bit-flip noise as an integer programming problem.Let us use the fact that the Pauli error E can be decomposed as as shown in Ref. [38].Here, T (S) is a pure error, which is a Pauli operator that returns a quantum state with the syndrome S to the code space.G is a stabilizer operator and L is a logical operator.The proof that the Pauli error can be decomposed as shown in Eq. ( 7) is as follows.T (S)E is an operator that acts within the code space, and operators acting within the code space can always be expressed as a product of a stabilizer operator and a logical operator.Therefore, for certain G and L, the relation T (S)E = GL always holds and it is shown that the Pauli error E can be decomposed as Eq. ( 7).Instead of considering the error E as a variable, we consider the stabilizer operator G in this decomposition formula as a variable.Thereby, the error E represented in this manner can always satisfy the syndrome S for an arbitrary stabilizer G.
While the definition of T (S) is not unique, it should be defined in a systematic way for different code distances to simplify the implementation of the decoding algorithm.Specifically, we define T (S) as follows (see Fig. 2).First, we define an error chain that only anti-commutes with the stabilizer generator on each type of face: • The blue face releases error chains from its vertex on the right side through blue edges to the right boundary.((a), (c) of Fig. 2) • The green face releases error chains from its vertex on the left side through green edges to the left boundary.((b), (d), (f) of Fig. 2) • The red square releases error chains from its left lower vertex through red edges to the left bottom, until the bottom edge.((e) of Fig. 2) Then, these error operators, which are defined uniquely for each face operator, are multiplied together for f with the syndrome value s f = 1 to define T (S) for a given syndrome S.An example of T (S) for d = 7 color codes is shown in Fig. 2.Although there are other strategies to choose T (S) [39], these are out of scope in this work because putting in significant effort into determining T (S) can make it difficult to evaluate the performance of SA itself.
In order to define a binary optimization problem, we rewrite the operators T (S), G, L in terms of binary vari- ables {t i (S)}, {g i }, l as follows: where t i (S) is the bit representing the action of T (S) on the i-th qubit, g i is a binary variable indicating whether or not the i-th X-stabilizer generator G i is included, and l is a binary variable indicating whether or not the Xlogical operator L X is included.Note that by changing l and g i , we can represent any operator that commutes with the Z-stabilizer and hence does not change the syndrome value.On the other hand, we can write the error E in terms of the binary variables {x i } as By combining these, we have where B i denotes the set of indices of the stabilizer generators acting on the i-th qubit, and l i is a binary variable by which the logical operator L X is given by L X = i X li i .By changing the binary variables to Ising spin variables via σ j := 1−2g j and J i := (1−2l i •l)(1−2t i (S)) in Eq. ( 12), we obtain the equation so the number of errors i x i is expressed as where N is the number of qubits.Accordingly, we can rewrite the number of errors i x i as an Ising Hamiltonian: where constant factor and term are omitted since they are not important.The Hamiltonian H l corresponds to the three-body Ising model.By minimizing the Hamiltonian H l in Eq. ( 15), it is possible to obtain an error configuration that minimizes the number of errors.It should be noted that H l depends on l, and it is necessary to consider which case, l being 0 or 1, yields the error configuration with the minimum number of errors.The error configuration that minimizes the Hamiltonian is determined for each case with l = 0, 1, and the error configuration achieving smaller energy between the two is the one that provides the minimum number of errors.Thus, we can decode for the bit-flip noise model by solving Therefore, the decoding problem can be formulated as an energy minimization problem of an Ising Hamiltonian.By solving Eq. ( 16) using SA, the optimized g j resulting in the minimum number of errors are obtained.
In the case of the depolarizing noise model, we can consider Z errors in exactly the same way as X errors, since X-stabilizer and Z-stabilizer are symmetric on the color codes.With respect to Y errors, if X and Z errors occur simultaneously on a qubit, such an event should be counted as a Y error.Therefore, the total number of X, Y , and Z errors is given by where z i ∈ {0, 1} is a binary variable that represents the presence or absence of the Z error in the i-th qubit.The mapping between T (S), G, L and binary variables are given as follows: where t X i (S), t Z i (S) are bits representing the action of T (S) on the i-th qubit with respect to X and Z, respectively.g X i and g Z i are binary variables representing whether or not the X-and Z-stabilizer generators G Xi and G Zi are included, respectively.l X and l Z are binary variables representing whether or not the X-and Z-logical operators L X and L Z are included, respectively.Similar to the discussion in the case of bit-flip noise, we can write the error E in terms of binary variables {x i }, {z i } as By combining these, x i and z i are represented as where l X i and l Z i are binary variables by which the logical operators L X and L Z are given by L 23) and ( 24), we convert the binary variables into spin variables (note that σ X j and σ Z j are not Pauli operators).Then, we obtain the equations so Eq. ( 18) is expressed as where N is the number of qubits.Accordingly, the number of errors, namely Eq. ( 18) can be represented as an Ising Hamiltonian: where Here, constant factors and terms are omitted, similar to the previous discussion.The Hamiltonian H l X ,l Z corresponds to the six-body Ising model.By minimizing the Hamiltonian H l X ,l Z in Eq. ( 30), it is possible to obtain an error configuration that minimizes the number of errors.
As in the case of the bit-flip noise model, to obtain the error configuration with the minimum number of errors, we should consider each possible combination of l X and l Z .
Considering the minimization of the Hamiltonian in all cases, the error configuration with the smallest Hamiltonian is the one with the minimum number of errors.Thus, we can decode for the depolarizing noise model by solving min Again, the decoding problem can be formulated as an energy minimization problem of an Ising Hamiltonian.By solving Eq. ( 34) using SA, g X j and g Z j optimized so that the number of errors becomes minimum can be obtained.

IV. DECODING ALGORITHM FOR PHENOMENOLOGICAL NOISE
A. Combinatorial optimization problem for phenomenological noise So far, we have assumed that the error occurs only in the data qubits and that the perfect syndrome is obtained.However, in reality the syndrome measurement is also subject to errors.Here, we extend the proposed decoding method to a phenomenological noise model that introduces errors into the measured syndrome phenomenologically.The phenomenological noise model is a noise model in which bit-flip errors occur on data qubits with a probability p, and measured syndrome also flips with the same probability p.In this situation, it is possible to achieve decoding by repeating the syndrome measurement process a number of times equal to the code distance d.For simplicity, we assume that syndrome measurement at the final round is obtained perfectly.Then, we obtain the error configuration that minimizes the total number of data errors and measurement errors satisfying the syndrome over all time.Note that the syndrome values themselves at each time step reflect the accumulation of data errors through past rounds.The data errors occurred at each time step can be extracted by taking the difference (XOR) of the syndrome values from one time step to the next.
In the following, the cumulative data error occurring in the i-th qubit up to time t is denoted by w (t) i ∈ {0, 1}, the measurement errors that occurred on the face f by r f ∈ {0, 1} at time t.The measured syndrome values at time t and t − 1 can be expressed as The difference s f and s (t−1) f implies data error at time t.Also, the difference between w (t−1) i and w (t) i corresponds to the newly occurring data error at time t, so we write it as x (t) i : Thus, the syndrome condition that should be satisfied is Therefore, the problem of decoding in the phenomenological noise model can be formulated as an integer programming problem with a constraint: As with the noise models we have been dealing with, this constrained combinatorial optimization problem is mapped to an unconstrained optimization problem by changing the variables.

B. Decoder for phenomenological noise model
Even in situations where temporal errors are considered, data errors represented as binary variables can be expressed in a similar form as Eq. ( 12).We can decompose temporal data error vector x i : where each t-th element of the vectors x i , ti (S), g i , and l i , corresponds to variables at time t, e.g.x i .S denotes the syndrome for all time steps.Regarding the pure error, we should consider operations in each time slice in a similar way to the bit-flip noise model.We write a bit that represents the action of a pure error on the i-th qubit at time t as t(t) i (S).With regard to stabilizer operators, the previous approach alone does not work in this case.In the presence of measurement errors, we have to consider not only stabilizer generators defined on each time The time-like stabilizer generator.A data error in a qubit occurs followed by measurement errors in all faces containing that qubit, and then a data error occurs in the same qubit at the next time.
slice but also spatio-temporal error events to cover all transformations that preserve the space-time syndrome values.Here, we call the stabilizer generators acting in the same time slice as "space-like stabilizer generators", and the stabilizer generators that act along the time axis as "time-like stabilizer generators", as described in Fig. 3. Now, in order to prove the decomposition as shown in Eq. ( 41), we provide a proof that any error event including data errors and measurement errors is expressed as a product of a pure error, a space-like stabilizer operator, a time-like stabilizer operator, and a logical operator.To prove this, it is sufficient to show that any error event including data error and measurement error that makes the left side of Eq. ( 38) equal to zero can always expressed as a product of a space-like stabilizer operator, a timelike stabilizer operator, and a logical operator.This is because the quantum state after application of pure errors at each time slice has trivial space-time syndrome values.Firstly, by the definition of time-like stabilizer generators which will be explained in more detail in the following paragraph, any error events including arbitrary patterns of r (t) f for all f and t that make the left side of Eq. ( 38) equal to zero can be expressed using time-like stabilizer generators.Then, for each pattern of r (t) f for all f and t, every patterns of x (t) i for all i and t that make the left side of Eq. ( 38) equal to zero can be expressed using space-like stabilizer operators and logical operators, similar to the discussion of Eq. (7).Therefore, it is proven that any error event including data errors and measurement errors is expressed as a product of a pure error, a space-like stabilizer operator, a time-like stabilizer operator, and a logical operator, which in turn show that the temporal data error vector x i can be decomposed as Eq.(41).
Time-like stabilizer generators are operators that gen-erate a data error at time t in a qubit followed by measurement errors in all faces containing that qubit, and then a data error occurs in the same qubit at time t + 1.
There are n time-like stabilizer generators between each pair of time steps, so there are n(d − 1) of them in total for all time steps.With this simple definition, the order of the interaction of the Hamiltonians is minimized and also it becomes easy to construct the Hamiltonians.Then, by utilizing both time-like and space-like stabilizer generators, the action of the stabilizer operator to the ith qubit at time t can be represented as where g ′ (t) j is a binary variable representing whether or not the space-like stabilizer generator at time t is included, and ḡ(t) i is a binary variable representing whether or not the time-like stabilizer generator between time t and t + 1 is included.B i denotes the set of indices of space-like stabilizer generators acting on the i-th qubit, and this set does not depend on the time t.Note that the logical operator can be defined at an arbitrary one time slice, since its time can be changed by using time-like stabilizers.Therefore, the action of the logical operator to the i-th qubit at time t can be represented as where l is a binary variable representing whether or not the logical operator at a certain time is included, l is a bit that indicates on which qubits and at what time the logical operator acts.Here, l does not need to depend on time t.Then, the data error x (t) i at time t can be written as Also, measurement errors can be represented using timelike stabilizer generators.The measurement error r f is represented as By defining σ ′ (t) i , and ) to transform binary variables into Ising spin variables, the total number of errors is represented as By minimizing Eq. ( 47), it is possible to obtain an error configuration that minimizes the number of errors.Similar to the previous discussion, to obtain the error configuration with the minimum number of errors, we should consider each case, whether l is 1 or 0. Considering the minimization of the Hamiltonian in the two cases, the error configuration with the smaller Hamiltonian is the one with the minimum number of errors.Thus, we can decode for the phenomenological noise model by solving min l (min H l,pheno ) .
Again, the decoding problem can be formulated as an energy minimization problem of an Ising Hamiltonian.
The Hamiltonian H l,pheno corresponds to the eight-body Ising model.By solving Eq. ( 48) with SA, g ′ (t) j and ḡ(t) j optimized so that the number of errors becomes minimum can be obtained.This decoder can also be applied for the circuit-level noise model, with the ability to adjust the weights of data errors and measurement errors based on the probability distribution after error propagation within the circuit.By changing the coefficients of each term in the Ising Hamiltonian, it is possible to adjust the weights because they correspond to the weights of each data error and measurement error.

A. Settings
In order to evaluate the performance of the proposed decoding methods, we perform Monte Carlo simulations to estimate the logical error rates.In the Monte Carlo simulations, errors were generated on the data qubits, and syndrome measurement was performed with and without errors depending on the noise models.Then, the Ising Hamiltonian was constructed from the observed syndrome for each of the noise models and was solved by SA.Specifically, we employed the open-source library OpenJij in Python [40].The solution is used to perform a recovery operation and to see whether or not error correction fails.These Monte Carlo simulations are repeated 10 5 times to estimate the logical error rates accurately.Moreover, in order to validate the accuracy of our method, we performed the same task to estimate the logical error rates by solving the constrained optimization problem exactly by using an integer programming solver, CPLEX.
Let us describe the details of the schedule employed in SA.The initial spin configuration is randomly determined.Also, we employed annealing parameters, such as how to schedule temperatures, how many repetitions of spin updates are taken at each temperature, and how many iterations of the overall process are performed to obtain the best solution of them.In SA, there is a tradeoff relationship between accuracy and computation time.We used the annealing parameters realizing the shortest possible computation time while keeping the accuracy reasonably good compared to that obtained by CPLEX.Specifically, we adopted the annealing parameters as follows.The initial and final values of the inverse temperature are β min = log 2/∆E min and β max = log 100/∆E max , respectively, where ∆E min is a rough lower bound of the energy gap of the Hamiltonians and ∆E max is an upper bound of the energy gap of the Hamiltonians.These values are determined by the coefficients of each term in the Hamiltonians and are the values commonly adopted in OpenJij, respectively.The number of Monte Carlo steps at each inverse temperature is taken to be one, where one Monte Carlo step means to update all spins once.We used the exponential cooling schedule to compute the temperature values at each temperature cycle.The number of temperature cycles and the number of iterations are changed depending on the noise models and solid line data represents our data and the dotted line data represents the data from CPLEX.The threshold obtained by our method is 10.36(5)% decoder with union-find [24].Although our threshold is lower than the optimal threshold of 10.925(5)% [22], it is not a surprising result.While color codes are degenerate codes, the proposed method is based on the concept of minimum distance decoding and does not consider the degeneracy of error configurations, contributing to the slightly lower threshold.Regarding the critical exponent, it belongs to the same universality class as the Ising critical exponent, which characterizes the behavior of the Ising model near its critical point [41].Although our critical exponent ν 0 is nearly consistent with the value of ν 0 =1.463 (6) found for the surface code [41], due to the large impact of statistical errors, it is essentially unclear whether these values are consistent or not.Compared to the logical error rates obtained by CPLEX, the proposed decoder provides almost the same accuracy within statistical errors, which clearly shows that the proposed decoder with SA works well achieving minimum distance decoding.
To compare the decoding time of SA and CPLEX, we show the decoding time for code distances d = 3, 5, 7, 11, 15 with the physical error rates p = 0.01 (a), p = 0.04 (b), p = 0.07 (c), and p = 0.10 (d) in Fig. 5.To make a fair comparison, we compared the computation time using a single core of a CPU for both cases, while our method can be simply parallelized when mul- tiple cores are available.The proposed method with SA succeeded to decode faster than CPLEX when the physical error rate is around the threshold or when d is small.The reason why the proposed method may be slower than CPLEX in cases where the physical error rate is low is that the SA parameters are optimized around the threshold.If we optimize the annealing parameters for each physical error rate, we can decode faster even in cases of low physical error probability, while this is out of scope in this work.Also, we can decode faster by taking advantage of the parallel computation of SA.SA can be parallelized with 100% efficiency for iterations.In addition, each iteration can be parallelized and accelerated using GPUs [42].

C. Depolarizing noise model
Next, we show the logical error rates for d = 3, 5, 7, 11 as functions of a physical error rate p under the depolarizing noise model in Fig. 6.The parameters of SA used in this process are listed in Table II.Again, we fit our data to the formula Eq. ( 50).We found which corresponds closely to the threshold value of 18.6(3)% that we estimated by the integer program decoder.Our threshold is higher than the threshold value of 17.5% found by the neural-network decoder [26].For the same reasons as in the case of the bit-flip noise model, our threshold is less than the value of 18.9% for optimal decoding [25].Moreover, similar to the previous case, the logical error rates of the proposed decoder are almost the same as those obtained by CPLEX.We demonstrate the decoding time measured using a single core of the CPU for code distances d = 3, 5, 7, 11 with the physical error rates p = 0.01 (a), p = 0.07 (b), p = 0.13 (c), and p = 0.18 (d) in Fig. 7. Similar to the case of the bit-flip noise model, our method achieved shorter computation time compared to CPLEX when either the physical error rate is around the threshold, or when d is not large.

D. Phenomenological noise model
Finally, we show the logical error rates for d = 3, 5, 7 as a functions of a physical error rate p under the phenomenological noise model in Fig. 8.In this noise model, bit-flip errors occur on data qubits with a probability p, and measured syndrome also flips with the same probability p.The parameters of SA used in this process are listed in Table III.Again, we fit our data to the form Eq. (50).We found p th = 0.0290 ± 0.0004, (55) which is in good agreement with the threshold values of 3.05(4)% estimated by the integer program decoder in Ref. [20].Compared to other previous studies, our The threshold obtained by our method is 2.90(4)% threshold is higher than the value of 2.08% estimated by graph matching decoder [27].The critical exponent also closely agrees with the value of ν 0 = 1.5(2) previously estimated by the integer program decoder [20].As with the bit-flip and depolarizing noise model, the proposed decoder shows logical error rates that are comparable to those obtained by CPLEX.The decoding time for code distances d = 3, 5, 7 with the physical error rates p = 0.010 (a), p = 0.018 (b), p = 0.024 (c), and p = 0.028 (d) is shown in Fig. 9.In this noise model, we succeeded to decode faster than CPLEX when d is small for almost all physical error rates.While our method is slower than CPLEX with large d, OpenJij can be further optimized in solving the energy minimization problem of the multi-body Ising Hamiltonian.In the case of the phenomenological noise model, the Hamilto-  nian consists of eight bodies.It is believed that we can decode faster by constructing a SA solver that is optimized for the multi-body Ising Hamiltonians instead of using OpenJij.Also, similar to the previous discussion, we can decode faster by tuning the parameters depending on the physical error rate or performing parallel computation.

VI. CONCLUSION
In this paper, we proposed an Ising model formulation for highly accurate color codes decoding, which is solved by SA.The decoding results show that our method can achieve almost the same accuracy as the integer program decoder using CPLEX for all noise models we have employed here if the annealing schedules are appropriately chosen.The decoding time is smaller than CPLEX when the physical error rate is around the threshold or when d is small for the bit-flip and depolarizing noise model.Also, for the phenomenological noise model, we succeeded to decode faster than CPLEX when d is small for almost all physical error rates.
In this work, the number of temperature cycles and the number of iterations in the SA are only optimized.However, the other parameters such as the initial and target inverse temperature, the number of Monte Carlo steps at each inverse temperature, etc., are not optimized but set to be the default values in OpenJij.By setting the target inverse temperature to the Nishimori temperature [43,44], the proposed decoder leads to a performance that is closer to the optimal one, taking into consideration the degeneracy.The proper setting of such annealing parameters has been a long-standing research topic, and several methods have been proposed [45,46], but no optimal method has been established.Therefore, we can decode more accurately and faster by further optimizing the annealing parameters.In this study, we did not optimize the SA solver itself for our purpose, but used the versatile open-source software OpenJij.However, we can further reduce the decoding time by developing an SA solver specific to the Ising model associated with the color code decoding problem and by parallelizing it with multicore CPUs and GPUs.Additionally, instead of selecting parameters that achieve high performance, we can use parameters that result in lower performance but shorter decoding This is because our decoder using SA is in a trade-off relationship between performance and decoding time.
Here we should that it is a rare result that we can solve combinatorial optimization problems through a approach with a smaller time compared to exact solvers [47].The reason is in general, when solving optimization problems SA, the problem needs to be embedded an Ising form, which leads to a loss.our study achieved favorable results with SA, which can be attributed to the compatibility between SA the error correction problem.In our method, the error correction problem is formulated as an Ising problem, and SA can be executed without the loss of problem embedding, resulting in a good performance.
Constructing high-performance decoding methods for the color codes presented here has been a challenging task, unlike the surface codes.However, our proposed decoding method enables highly accurate decoding to be achieved in a shorter amount of time compared to other decoders achieving comparable performance.Although color codes are still suffering from low thresholds in the circuit-level noise model, our method simplifies threshold estimation and makes it easier to improve the performance through trial and error in the architecture design.Recently, QEC on small color codes has been experimentally implemented [48].In such near-term QEC experiments, there is a need to demonstrate that decoders can, in principle, handle correlated errors that often arise when applying e.g., transversal entangling gates.It is reasonable to spend a relatively large amount of time on decoding for this purpose, particularly around the threshold.In this situation, our decoder especially has advantages, as it allows adjusting performance and decoding time within a trade-off by tuning parameters.Furthermore, the decoding method proposed in this paper can be readily applied to any stabilizer code and hence has a great potential to advance the field of QEC, including quantum low-density parity check codes [49], where the decoding problem is highly non-trivial.Further studies can explore the scalability of the method and its effectiveness in real-world scenarios, with the aim of improving the performance of quantum computing systems.

FIG. 2 .
FIG. 2. Definition of T (S) for d = 7 color codes.The figure shows the correspondence between certain faces and the error chains associated with the construction of T (S).(a) The blue octagon.(b) The green octagon.(c) The blue trapezoid other than the bottom.(d) The green trapezoid other than the bottom.(e) The red square.(f) The bottom green trapezoid.

B
. Decoder for the depolarizing noise model So far we have only considered the bit-flip noise model.Next, we extend the previous argument to a more general noise model, the depolarizing noise model.The depolarizing noise model is a noise model in which the Pauli X, Y , and Z operators act on each physical qubit with the probability p x = p y = p z = p/3: 1} at time t, and the measured syndrome values on face f by s (t)

FIG. 3 .
FIG. 3. The space-like stabilizer generator and the time-like stabilizer generator.(a) The space-like stabilizer generator.The four data errors indicate the action of space-like stabilizer generators defined on a blue face.(b)The time-like stabilizer generator.A data error in a qubit occurs followed by measurement errors in all faces containing that qubit, and then a data error occurs in the same qubit at the next time.

FIG. 4 .
FIG. 4. Numerical simulations of logical error rate for the bit-flip noise model for various code distances.solidline data represents our data and the dotted line data represents the data from CPLEX.The threshold obtained by our method is 10.36(5)%

TABLE I .
The parameters of SA used for decoding the bit-flip noise model.

TABLE III .
The parameters of SA used for decoding the phenomenological noise model.