In the search for scalable, fault-tolerant quantum computing, distributed quantum computers are promising candidates. These systems can be realized in large-scale quantum networks or condensed onto a single chip with closely situated nodes. We present a framework for numerical simulations of a memory channel using the distributed toric surface code, where each data qubit of the code is part of a separate node, and the error-detection performance depends on the quality of four-qubit Greenberger–Horne–Zeilinger (GHZ) states generated between the nodes. We quantitatively investigate the effect of memory decoherence and evaluate the advantage of GHZ creation protocols tailored to the level of decoherence. We do this by applying our framework for the particular case of color centers in diamond, employing models developed from experimental characterization of nitrogen-vacancy centers. For diamond color centers, coherence times during entanglement generation are orders of magnitude lower than coherence times of idling qubits. These coherence times represent a limiting factor for applications, but previous surface code simulations did not treat them as such. Introducing limiting coherence times as a prominent noise factor makes it imperative to integrate realistic operation times into simulations and incorporate strategies for operation scheduling. Our model predicts error probability thresholds for gate and measurement reduced by at least a factor of three compared to prior work with more idealized noise models. We also find a threshold of $ 4 \xd7 10 2$ in the ratio between the entanglement generation and the decoherence rates, setting a benchmark for experimental progress.

## I. INTRODUCTION

A *distributed quantum computer*^{1,2} realizes a large-scale processing system by using entanglement to link smaller quantum processing units. For example, the sub-units may be elements of a photonic chip or form the nodes of a quantum network on a larger scale.^{3} Fault tolerance is naturally achieved by establishing the connectivity according to the architecture of a topological error-correction code. The distributed approach provides advantages in terms of scalability but is limited by the availability of high-quality entanglement. This is because entangled states are required for both inter-node operations and for detecting local errors with the error-correction code.^{4–9}

In this paper, we focus on systems that are capable of generating remote two-qubit entanglement between pairs of connected nodes. There exist several physical systems suitable for generating this type of entanglement with optical interfaces.^{10} Examples of this are ion traps, neutral atoms, and *color centers* in the diamond lattice, such as nitrogen-vacancy (NV),^{11–22} silicon-vacancy (SiV),^{23–27} or tin-vacancy (SnV)^{28–30} centers. As a concrete example, we investigate the distributed surface code with hardware-specific noise based on color centers, also known as *defect centers*. We emphasize that the obtained insights are more general and that our simulation tools allow for implementing error models based on general hardware implementations.

Diamond color centers host long-lived electron spins that exhibit coherent optical transitions, enabling their use as a *communication* qubit. This qubit is used to create entanglement with other nodes and can address up to tens of proximal nuclear spins and or other electron spins occurring in the host material.^{15,19,27} These nuclear spins can be used as local processing qubits—i.e., as a local memory to store and manipulate quantum states.^{16} Hereafter, in the context of diamond color centers, the term “memory qubits” specifically refers to such spins.

Physical systems suitable for distributed quantum computing can be operated as *fault-tolerant* quantum computers by employing a subset of their memory qubits as data qubits of an *error-correction code*, such as the toric surface code.^{31,32} The principle behind error-correction codes is that many physical qubits hold a smaller number of *logical* states, and unwanted errors can be detected and corrected by measuring the *stabilizer operators* of the code. For such a system, fault-tolerance goes hand-in-hand with the existence of *thresholds* for local sources of error: if one manages to keep the error sources below their threshold values, one can make the logical error rate arbitrarily small by increasing the dimension of the error-correction code.

The toric code has a depolarizing *phenomenological* error probability threshold of approximately 10% to 11%.^{33} This error model assumes that all qubits of the code are part of the same quantum system, stabilizer measurements can be carried out perfectly, and the qubits experience depolarizing noise in between the stabilizer measurement rounds. A more precise analysis with a *circuit-level* error model yields error probability thresholds between 0.9% and 0.95%.^{34} In this model, the stabilizer measurement circuit is simulated with noisy gates and measurements, and it is implicitly assumed that the connectivity of the system allows direct two-qubit gates between adjacent qubits of the code topology. Therefore, this error model corresponds to a *monolithic* architecture.

If one wants to implement the toric code in a network setting, where every data qubit of the code is part of a separate network node, the stabilizer operators can be measured with the aid of four-qubit Greenberger–Horne–Zeilinger (GHZ) states. These GHZ states can be created by *fusing* three or more Bell pairs created between the involved nodes. Nickerson *et al.* analyzed the distributed toric code in this setting.^{34,35} They included protocols with a relatively large number of *entanglement distillation* steps that create high-quality GHZ states from imperfect Bell pairs. They found^{34} thresholds for the local qubit operations between 0.6% and 0.82%—i.e., slightly below the monolithic thresholds. In their threshold calculations, Nickerson *et al.* do not explicitly consider circuit operation times and do not include qubit memory *decoherence* during entanglement creation—i.e., the notion that the quality of the code's data qubits decays over time. However, in current physical systems of interest, decoherence during entanglement creation typically constitutes a large source of error. For state-of-the-art NV centers, coherence times during optical Bell pair generation are one to two orders of magnitude lower than estimated by Nickerson *et al.*^{36–38} The influence of this decoherence is further increased by the reality that success probabilities per optical Bell pair generation attempt currently fall significantly short of unity.^{15,39}

Therefore, next to the errors in operations and in entangled states considered in Refs. 34 and 35, decoherence of quantum states over time emerges as the third primary source of noise for accurate assessment of distributed quantum computing systems. The influence of memory qubit decoherence during entanglement creation can be captured with the *link efficiency* $ \eta link *$.^{21} This parameter quantifies the average number of entangled pairs that can be generated within the coherence times.

To investigate the influence of the coherence times, we develop a time-tracking simulator and implement realistic operation durations. Additionally, considering the pivotal role of the operation order in this new scenario, we formulate a strategy for scheduling operations. We find that, with realistic operation and coherence times, the thresholds with the GHZ generation protocols of Refs.34 and 35 disappear. We investigate the quantitative impact of memory decoherence and optimize over GHZ generation protocols with less distillation that can overcome this. For a range of different coherence times during entanglement generation, we find two-qubit gate error and measurement error probability thresholds for diamond color centers up to 0.24%. We find that fault-tolerance is reachable with $ \eta link * \u2248 4 \xd7 10 2$. This improves on the prior results of $ \eta link * = 2 \xd7 10 5$ for the idealized time scale estimates of Nickerson *et al.*^{34} However, this link efficiency is still above the state-of-the-art hardware^{21} reaching up to $ \eta link * \u2248 10$.

In the remainder of this paper, Sec. II describes GHZ creation and distillation protocols necessary for the distributed surface code. Consequently, in Sec. III, we present the full cycle of stabilizer measurements of the surface code. In Sec. IV, we describe error models that allow us to investigate a specific hardware implementation in the distributed surface code setting: diamond color centers. Finally, in Sec. V, we investigate the parameter regimes necessary for fault tolerance with these error models.

## II. GHZ GENERATION PROTOCOLS

As mentioned in Sec. I, the stabilizer operators of a distributed quantum error-correcting code can be measured by consuming GHZ states. In the following, we discuss protocols that create GHZ states by combining Bell pairs. For each GHZ protocol, we identify two parameters. The first one is the minimum number of Bell pairs *k* required to create the GHZ state. This number indicates the amount of distillation taking place in the protocol. The second one is the maximum number of qubits per node *q* necessary to generate the GHZ state. We summarize prior work in Sec. II A. In Sec. II B, we discuss our method for generating GHZ protocols.

### A. Prior GHZ protocols

There is a plethora of prior work considering the generation and purification of GHZ states.^{40–53} Here, we focus on protocols that combine Bell pairs into a four-qubit GHZ state and discuss seven of them.

First, we consider two protocols that we used in an earlier study:^{21} the Plain (*k* = 3, *q* = 2) and Modicum (*k* = 4, *q* = 2) protocols. These protocols were designed to create a GHZ state with no distillation or only a single distillation step. The Plain protocol is the simplest protocol for creating a GHZ state from Bell pairs; it fuses three Bell pairs into a four-qubit GHZ state without any distillation. The Modicum protocol uses a fourth Bell pair to perform one round of distillation on the GHZ state.

### B. Dynamic program to generate GHZ protocols

In this section, we present a method for optimizing GHZ creation with realistic noise models. We focus on creating GHZ states of the form $ | GHZ n \u27e9 = ( | 0 \u27e9 \u2297 n + | 1 \u27e9 \u2297 n ) / 2$, where $ n \u2208 { 2 , 3 , 4}$ represents the number of parties. We call *n* the weight of the GHZ state. For convenience, we use the notation $ | GHZ 2 \u27e9$ to describe a Bell state.

We use a *dynamic program* to optimize over the space of GHZ protocols. This program generates GHZ protocols by using two main operations: *fusing* Bell pairs to create GHZ states, and *distilling* or *purifying* Bell or GHZ states by consuming other ancillary Bell pairs or GHZ states. Figure 1 depicts the two building blocks.

Distillation involves the use of an ancillary state to non-locally measure a stabilizer of the main Bell or GHZ state.^{54} In this process, local control-Pauli gates between ancillary and main state qubits are followed by individual measurements of the ancillary state qubits in the Pauli-*X* basis. Obtaining an even number of −1 measurement outcomes marks a successful distillation attempt. If distillation fails, the post-measurement state is discarded and (part of) the protocol has to be carried out again.

Fusion is executed to create GHZ states out of Bell pairs and to create a larger GHZ state. A state of the form $ | GHZ n 1 \u27e9$ can be fused with a state $ | GHZ n 2 \u27e9$ by applying a CNOT gate between one qubit of both states and measuring out the target qubit in the Pauli-*Z* basis. Obtaining a +1 measurement outcome results in the state $ | GHZ n 1 + n 2 \u2212 1 \u27e9$. A −1 measurement outcome leads to the same state after local Pauli-*X* corrections.

In Algorithm 1, we present a schematic, pseudo-code version of the dynamic program we used to generate and evaluate GHZ protocols. This algorithm is an expanded version of the dynamic program in Ref. 55. In this algorithm, each protocol is created with either a fusion or a distillation operation that combines two smaller GHZ protocols encountered earlier in the search. The protocols created in this fashion can be depicted with a directed *binary tree* graph. An example graph is given on the left side of Fig. 2. For the distillation steps in the binary tree diagrams, we consume the state on the right to distill the state on the left.

Require: $ n m$: number of qubits of final GHZ state |

$ k m$: minimum number of Bell pairs used |

V: set with model parameters used |

$ N b$: protocols stored in buffer per step |

$ N so$: Monte Carlo shots used per protocol |

for $ { ( n , k ) \u2009 | \u2009 2 \u2264 n \u2264 n m , \u2009 n \u2212 n m + k m \u2265 k \u2265 n \u2212 1}$ do |

# Try all non-local measurement combinations. |

for g $\u2208$ {stabilizers of $ | GHZ n \u27e9$} do |

$ n \u2032$ ← weight of g |

for $ k \u2032 \u2208 [ n \u2032 \u2212 1 , \u2009 k \u2212 n + 1 ]$ do |

for $ ( p 1 , p 2 ) \u2208 [ 1 , N b ] \xd7 [ 1 , N b ]$ do |

$ P 1$ ← protocol p_{1} in buffer at $ ( n , k \u2212 k \u2032 )$ |

$ P 2$ ← protocol p_{2} in buffer at $ ( n \u2032 , k \u2032 )$ |

Construct binary tree protocol $ P new$ that |

measures g on $ P 1$ by consuming $ P 2$ |

Construct protocol recipe $ R new$ and evaluate |

quality over $ N so$ iterations times using V |

Store protocol if average performance |

is better than worst protocol in buffer |

# Try all fusion combinations. |

for $ n 2 \u2208 [ 2 , n \u2212 1 ]$ do |

$ n 1 \u2190 n \u2212 n 2 + 1$ |

for $ k 2 \u2208 [ n 2 \u2212 1 , k \u2212 n + 1 ]$ do |

$ k 1 \u2190 k \u2212 k 2$ |

for $ ( p 1 , p 2 ) \u2208 [ 1 , N b ] \xd7 [ 1 , N b ]$ do |

$ P 1$ ← protocol p_{1} in buffer at (n_{1}, k_{1}) |

$ P 2$ ← protocol p_{2} in buffer at (n_{2}, k_{2}) |

for $ ( i , j ) \u2208 [ 1 , n 1 ] \xd7 [ 1 , n 2 ]$ do |

Construct binary tree protocol $ P new$ by |

fusing $ P 1$ at qubit i and $ P 2$ at qubit j |

Construct protocol recipe $ R new$ and evaluate |

quality over $ N so$ iterations using V |

Store protocol if average performance |

is better than worst protocol in buffer |

Require: $ n m$: number of qubits of final GHZ state |

$ k m$: minimum number of Bell pairs used |

V: set with model parameters used |

$ N b$: protocols stored in buffer per step |

$ N so$: Monte Carlo shots used per protocol |

for $ { ( n , k ) \u2009 | \u2009 2 \u2264 n \u2264 n m , \u2009 n \u2212 n m + k m \u2265 k \u2265 n \u2212 1}$ do |

# Try all non-local measurement combinations. |

for g $\u2208$ {stabilizers of $ | GHZ n \u27e9$} do |

$ n \u2032$ ← weight of g |

for $ k \u2032 \u2208 [ n \u2032 \u2212 1 , \u2009 k \u2212 n + 1 ]$ do |

for $ ( p 1 , p 2 ) \u2208 [ 1 , N b ] \xd7 [ 1 , N b ]$ do |

$ P 1$ ← protocol p_{1} in buffer at $ ( n , k \u2212 k \u2032 )$ |

$ P 2$ ← protocol p_{2} in buffer at $ ( n \u2032 , k \u2032 )$ |

Construct binary tree protocol $ P new$ that |

measures g on $ P 1$ by consuming $ P 2$ |

Construct protocol recipe $ R new$ and evaluate |

quality over $ N so$ iterations times using V |

Store protocol if average performance |

is better than worst protocol in buffer |

# Try all fusion combinations. |

for $ n 2 \u2208 [ 2 , n \u2212 1 ]$ do |

$ n 1 \u2190 n \u2212 n 2 + 1$ |

for $ k 2 \u2208 [ n 2 \u2212 1 , k \u2212 n + 1 ]$ do |

$ k 1 \u2190 k \u2212 k 2$ |

for $ ( p 1 , p 2 ) \u2208 [ 1 , N b ] \xd7 [ 1 , N b ]$ do |

$ P 1$ ← protocol p_{1} in buffer at (n_{1}, k_{1}) |

$ P 2$ ← protocol p_{2} in buffer at (n_{2}, k_{2}) |

for $ ( i , j ) \u2208 [ 1 , n 1 ] \xd7 [ 1 , n 2 ]$ do |

Construct binary tree protocol $ P new$ by |

fusing $ P 1$ at qubit i and $ P 2$ at qubit j |

Construct protocol recipe $ R new$ and evaluate |

quality over $ N so$ iterations using V |

Store protocol if average performance |

is better than worst protocol in buffer |

Each binary tree corresponds to multiple inequivalent protocols depending on the time ordering of the steps. We define a *protocol recipe* as a set of instructions for implementing the protocol. The recipe includes the ordering of operations and state generation. An example of a protocol recipe can be seen on the right side of Fig. 2. This step was not required in previous research on distributed surface codes, as the noise models used in previous research did not include memory decoherence. Without a notion of time, the execution order of the tree's branches is irrelevant.

As can be seen in Fig. 2, the conversion to a protocol recipe contains SWAP gates. These gates are required to match the connectivity constraints of our example hardware model—see Sec. IV for more details. The SWAP gates should, therefore, not be considered as fundamental elements of these protocols and can be circumvented or neutralized in hardware systems with more operational freedom. We implement SWAP gates as three CNOT gates.

Although we did not optimize over the conversion from binary tree to protocol recipe, we considered two heuristics to limit the influence of decoherence and of SWAP gates. To limit decoherence, we prioritize creating larger branches of the tree. Here, a branch is defined as an element of the binary tree (i.e., an operation in the GHZ protocol) including all elements that (in)directly point toward it. The size of the branch is the number of elements it contains. Because, generally speaking, creating small branches is faster than creating large branches, this heuristic aims to minimize waiting times for completed intermediate branches of the GHZ protocol.

The SWAP gate count can be limited by making sure a Bell pair that is consumed in a distillation step is the last state to be generated. This prevents the protocol from having to swap this state back and forth between the memory. For this reason, if two branches have equal size, we prioritize creating the left branch over the right one.

In constructing the protocol recipe, we first use these heuristics to determine the order in which the elementary Bell pairs are generated—i.e., the leaves of the binary tree. By following this order, we then check for each Bell pair if other Bell pairs in non-overlapping network nodes can be generated simultaneously. Here, we prioritize based on proximity in the binary tree. We include instructions for distillation, fusion, and SWAP operations at the earliest possible point of execution. This approach gives rise to a unique conversion from binary tree to protocol recipe. More detailed descriptions of the protocol recipe construction and execution procedure can be found in Ref. 56 and the supplementary documents of the repository of Ref. 57.

While our dynamic program explores a large class of protocols, not all of the seven protocols that we introduced in Sec. II A can be generated. This is because, to suppress calculation time, the program limits distillation steps to operations that use an ancillary entangled state to non-locally measure a stabilizer operator of the main state, in a sequential manner. The protocols Refined, Expedient, and Stringent, however, make use of the so-called *double selection* distillation block^{58} that does not directly fit into this stabilizer distillation framework.^{59–61}

## III. DISTRIBUTED TORIC CODE

In this section, we discuss the steps of a distributed toric code and our approach for its simulation.

### A. The toric surface code

In the toric surface code,^{31,32} data qubits are placed on the edges of an *L* × *L* lattice with periodic boundary conditions. It encodes two logical qubit states. The stabilizers of the code come in two forms: the product of *X* operators on the four qubits surrounding every vertex of the lattice, and the product of *Z* operators on the four qubits surrounding every face (or plaquette) of the code.

We consider a network topology with node connectivity matching the connectivity of the toric code lattice. We present a schematic impression in Fig. 3. Each data qubit of the toric code is placed in a separate network node—e.g., in a separate diamond color center. The nodes have access to local memory qubits to create and store entangled links between them. Entangled links can be used to create four-qubit GHZ states, as described in Sec. II, which are then consumed to measure the stabilizers of the code. Figure 4 shows a depiction of the procedure.

The outcomes of the stabilizer measurements, known as the *error syndrome*, are fed to a decoder to estimate the underlying error pattern. Here, we consider an implementation by Hu^{62,63} of the *Union-Find*^{64} error decoder.

We point out that we simulate the toric surface code as a logical quantum *memory* and do not consider its resilience against, e.g., logical operations or operations required for initializing logical information. This means we restrict the study to the code's ability to protect general logical states.

We opted for the toric surface code over the planar surface code (i.e., the surface code with boundaries) because, on top of the weight-4 stabilizer operators in the bulk, the planar code has lower-weight stabilizers at its boundaries. For distributed implementations, measuring these lower-weight stabilizers requires additional entanglement protocols. This makes simulating the planar code more complicated. Studies reveal that the introduction of boundaries typically has a limited effect on the code's threshold values, yet it is likely to result in a slower suppression of the logical error rates below these thresholds.^{65}

### B. Distributed toric code simulation

We split the simulation of the distributed toric code into two levels: simulation of the toric code's stabilizer measurements with the aid of GHZ states and simulation of *L* rounds of stabilizer measurements of the code itself.

The first level characterizes the stabilizer measurement. To this aim, we use Monte Carlo simulation to construct the (average) *Choi state* associated with using the protocol's GHZ state to measure the plaquette or star stabilizer operator on half of the maximally entangled state. Exploiting channel-state duality, the Choi state from each iteration is converted to a *superoperator* describing the stabilizer measurement channel. A formal introduction on channel characterization with a maximally entangled state can be found in Ref. 66. The superoperator construction is described in detail in Appendix C.

The second level is a toric code simulator that takes as noise model the average superoperator obtained in the first level. Following previous research,^{34} we consider a stabilizer measurement cycle consisting of two rounds of plaquette stabilizer measurements and two rounds of star stabilizer measurements. This is because the constraint that each network node only has one single communication qubit in our example hardware model makes it impossible to simultaneously generate entanglement for overlapping stabilizers—see Sec. IV for more details. By splitting the full cycle up into four rounds, each defect center becomes part of exactly one stabilizer measurement per round. This process is schematically depicted in Fig. 5. We note that a different hardware model could require, or benefit from, a different scheduling of the stabilizer measurements as the one used here.

Due to entanglement distillation steps in the GHZ creation protocol, GHZ generation is probabilistic. To fix the duration of the rounds we impose a “GHZ cycle time” $ t GHZ$: if GHZ generation is not successful within this time, it aborts. In that case, the corresponding stabilizer operator cannot be measured. This information could be given to the decoder in the form of an erasure symbol. However, to leverage existing decoders, we opt to duplicate the last measured value of the stabilizer. This choice is suboptimal and better thresholds could be expected for decoders that can handle erasures and noisy measurements. The GHZ cycle time is a free parameter that we explore in Sec. V. In Sec. 5 of Appendix B, we describe a heuristic-driven approach for selecting a suitable GHZ cycle time.

To model GHZ generation failures, at the first level, we construct two separate average superoperators per stabilizer type: a *successful* superoperator $ S \xaf success$ for iterations where the GHZ state is created within $ t GHZ$, and an *unsuccessful* superoperator $ S \xaf fail$ for iterations where the GHZ could not be created. Both superoperators incorporate the influence of decoherence up to the cycle time on the code's data qubits.

## IV. DIAMOND COLOR CENTER MODEL

In this section, we present an overview of the noise models and model parameters used in combination with the distributed surface code structure of Sec. III A. The noise models are based on the nitrogen-vacancy (NV) center.^{16,19,21} More information about the experimental characterization can be found in Appendix A. More details about the models can be found in Appendix B. The parameter values can be found in Table I.

. | State-of-the-art (Refs. 16 and 21) . | Figure 6 . | Figure 7 . | Figure 8 . |
---|---|---|---|---|

Bell pair model input | ||||

Protocol | Single-click (Ref. 68) | Double-click (Ref. 69) | ||

$ F prep$ | 0.99 (Refs. 39 and 71) | 0.999 | ||

$ p EE$ | 0.04 (Ref. 15) | $ p EE ( f \varphi )$ | 0.01 | |

μ | 0.9 (Refs. 15 and 19) | 0.95 | ||

λ | 0.984 (Refs. 15 and 70) | 1 | ||

$ \eta ph$ | 0.0046 (Ref. 72) | 0.4472 | $ \eta ph ( f \eta )$ | |

Bell pair model output | ||||

$ p link$ | 0.0001 | 0.1 | $ p link ( f \eta )$ | |

$ F link$ | 0.8966 | $ F link ( f \varphi )$ | 0.9526 | |

Operation durations | ||||

$ t link$ | $ 6 \xd7 10 \u2212 6$ s | |||

$ t meas$ | $ 4 \xd7 10 \u2212 6$ s | |||

$ t X , Y e$ | $ 0.14 \xd7 10 \u2212 6$ s | |||

$ t X , Y n$ | $ 1.0 \xd7 10 \u2212 3$ s | |||

$ t Z , H e$ | $ 0.1 \xd7 10 \u2212 6$ s | |||

$ t Z , H n$ | $ 0.5 \xd7 10 \u2212 3$ s | |||

$ t C Z , C X , C i Y$ | $ 0.5 \xd7 10 \u2212 3$ s | |||

$ t SWAP$ | $ 1.5 \xd7 10 \u2212 3$ s | |||

Decoherence | ||||

$ T 1 idle n$ | 300 s | |||

$ T 1 link n$ | 0.03 s (Ref. 19) | 0.3 s | $ 0.03 f dec$ s | 0.3 s |

$ T 1 idle e$ | 300 s | |||

$ T 2 idle n$ | 10 s | |||

$ T 2 link n$ | 0.0075 s (Ref. 19) | 0.075 s | $ 0.0075 f dec$ s | 0.075 s |

$ T 2 idle e$ | 1.0 s | |||

$ t pulse$ | $ 1.0 \xd7 10 \u2212 3$ s | |||

$ n DD$ | 500 | 18 | Eq. (B7) | |

Link efficiency—see Eq. (1) | ||||

$ \eta link *$ | $ 2 \xd7 10 \u2212 1$ | $ 2 \xd7 10 3$ | $ 200 f dec$ | $ \eta link * ( f \eta )$ |

Operation noise | ||||

$ p g$ | 0.01 | Threshold | ||

$ p m$ | 0.01 |

. | State-of-the-art (Refs. 16 and 21) . | Figure 6 . | Figure 7 . | Figure 8 . |
---|---|---|---|---|

Bell pair model input | ||||

Protocol | Single-click (Ref. 68) | Double-click (Ref. 69) | ||

$ F prep$ | 0.99 (Refs. 39 and 71) | 0.999 | ||

$ p EE$ | 0.04 (Ref. 15) | $ p EE ( f \varphi )$ | 0.01 | |

μ | 0.9 (Refs. 15 and 19) | 0.95 | ||

λ | 0.984 (Refs. 15 and 70) | 1 | ||

$ \eta ph$ | 0.0046 (Ref. 72) | 0.4472 | $ \eta ph ( f \eta )$ | |

Bell pair model output | ||||

$ p link$ | 0.0001 | 0.1 | $ p link ( f \eta )$ | |

$ F link$ | 0.8966 | $ F link ( f \varphi )$ | 0.9526 | |

Operation durations | ||||

$ t link$ | $ 6 \xd7 10 \u2212 6$ s | |||

$ t meas$ | $ 4 \xd7 10 \u2212 6$ s | |||

$ t X , Y e$ | $ 0.14 \xd7 10 \u2212 6$ s | |||

$ t X , Y n$ | $ 1.0 \xd7 10 \u2212 3$ s | |||

$ t Z , H e$ | $ 0.1 \xd7 10 \u2212 6$ s | |||

$ t Z , H n$ | $ 0.5 \xd7 10 \u2212 3$ s | |||

$ t C Z , C X , C i Y$ | $ 0.5 \xd7 10 \u2212 3$ s | |||

$ t SWAP$ | $ 1.5 \xd7 10 \u2212 3$ s | |||

Decoherence | ||||

$ T 1 idle n$ | 300 s | |||

$ T 1 link n$ | 0.03 s (Ref. 19) | 0.3 s | $ 0.03 f dec$ s | 0.3 s |

$ T 1 idle e$ | 300 s | |||

$ T 2 idle n$ | 10 s | |||

$ T 2 link n$ | 0.0075 s (Ref. 19) | 0.075 s | $ 0.0075 f dec$ s | 0.075 s |

$ T 2 idle e$ | 1.0 s | |||

$ t pulse$ | $ 1.0 \xd7 10 \u2212 3$ s | |||

$ n DD$ | 500 | 18 | Eq. (B7) | |

Link efficiency—see Eq. (1) | ||||

$ \eta link *$ | $ 2 \xd7 10 \u2212 1$ | $ 2 \xd7 10 3$ | $ 200 f dec$ | $ \eta link * ( f \eta )$ |

Operation noise | ||||

$ p g$ | 0.01 | Threshold | ||

$ p m$ | 0.01 |

In our model, qubits undergo both *generalized amplitude damping* and *phase damping* noise,^{67} with separate *T*_{1} and *T*_{2} times. The generalized amplitude damping channel decoheres to the maximally mixed state $ ( | 0 \u27e9 \u27e8 0 | + | 1 \u27e9 \u27e8 1 | ) / 2$. Decoherence of the electron qubit is governed by $ T 1 idle e$ and $ T 2 idle e$ coherence times. For the memory qubits, we use different coherence times: $ T 1 idle n$ and $ T 2 idle n$ for when the node is idling vs $ T 1 link n$ and $ T 2 link n$ during optical entanglement creation. This is because the required operations on the electron qubit typically induce additional decoherence on the memory qubits.^{19,21} The Kraus operators of both channels can be found in Sec. 2 of Appendix B.

To mitigate decoherence from quasi-static noise processes in diamond color center experiments, *dynamical decoupling* is typically interleaved into gate sequences on both the electron and nuclear spins.^{16} Here, the coherence times that we consider are also those achieved for NV center spin registers with dynamical decoupling.^{16} Consequently, in our numerical models, gate operations must be performed only between two consecutive dynamical decoupling pulses—i.e., at the refocusing point of the qubit spins involved in the operations. We define the center-to-center time of consecutive refocusing points as $ t DD = t pulse + 2 n DD t link$, where $ t pulse$ is the time of a *π*-pulse, $ t link$ is the duration of a single Bell pair generation attempt, and $ n DD$ is a fixed number of Bell pair generation attempts. In Sec. 3 of Appendix B, we discuss how $ n DD$ is optimized. In our model, we assume that all memory qubits of a node are decoupled synchronously.

We assume that each diamond color center only possesses a single communication qubit. Within each node, we further assume that measurements are only possible on its communication qubit, and local (i.e., intra-node) two-qubit gates always require the communication qubit to be the control qubit. These requirements mean we have to use SWAP gates to measure the state of a memory qubit or to use a memory qubit as the control of a two-qubit gate, as can be seen in Fig. 2. Finally, we assume that a Bell pair between two centers can only be generated between the communication qubits of the two centers. We note that, in general, one could design (combinations of) diamond color center nodes with multiple communication qubits. Although this limits the number of SWAP gates required for distillation, this gives rise to extra requirements on the memory robustness of the communication qubits.

The generation of Bell pairs between two color centers is modeled as a probabilistic process with a success probability $ p link$ and time $ t link$ per attempt. We constructed an analytic model to calculate $ p link$ and the Bell pair density matrix after success, both for the *single-click* (i.e., the *single-photon*)^{68} entanglement protocol and for the *double-click* (i.e., the *two-photon*)^{69} entanglement protocol. The following five noise sources are included in this analytic model: the preparation error of the initial spin-photon state $ F prep$, the probability of an excitation error $ p EE$, a parameter *λ* based on the standard deviation of the phase uncertainty due to the path-length difference between the two arms (all modeled as dephasing channels on the involved qubits), the photon indistinguishability *μ* for each Bell state measurement (i.e., Hong–Ou–Mandel visibility, modeled with altered measurement operators^{70}), and the total photon detection probability $ \eta ph$ in each arm (modeled with an amplitude damping channel).

For the double-click protocol, we assume the phase uncertainty to not be relevant and set *λ* = 1. The parameters $ F prep , \u2009 p EE$, and *μ* affect the fidelity $ F link$ of the Bell pair state of the double-click protocol, whereas the parameter $ \eta ph$ limits the success probability $ p link$ of a single entanglement attempt. For the single-click protocol, the fidelity $ F link$ is additionally influenced by $ \eta ph$ and *λ*, and $ p link$ depends on the indistinguishability *μ*. A full description of the density matrices and success probabilities of both entanglement protocols can be found in Sec. 1 of Appendix B.

We assume that all operations take a finite amount of time. The time durations can be found in Table I. We neglect the influence of classical communication times, as we consider distances between network nodes to be relatively small, but include synchronization of network nodes when classical communication is required. Furthermore, we assume single-qubit gates to be noiseless, while noise in two-qubit gates is modeled with a depolarizing channel with probability $ p g$ (see Sec. 4 of Appendix B). We model imperfect measurements by flipping the measurement outcome with probability $ p m$.

## V. RESULTS

In this section, we investigate the sensitivity of the distributed toric surface code performance with respect to several physical parameters of the diamond color center hardware model. In particular, we investigate the influence of two-qubit gate and measurement noise, the entanglement success probability, the coherence times during entanglement generation, and the quality of the generated Bell pairs on the noise thresholds. The threshold values $ p th$ are determined with fits of the logical success rates vs the lattice size (*L*) and local two-qubit gate and measurement error rate ( $ p g = p m$). The details of the fitting procedure can be found in Appendix D.

### A. Current state-of-the-art parameter set

Let us first consider a parameter set inspired by state-of-the-art NV center hardware (see the first column of Table I). The operation times in this set are based on typical time scales in nitrogen-vacancy centers with a natural ^{13}C concentration^{16,19,21}—see Appendix A for more details. The Bell pair parameter values are a collection of the best parameters in current NV center literature—see the first column of Table I for relevant citations. In the following, we explicitly refer to this parameter set as the “state-of-the-art” parameter set. As discussed in more detail in Appendix A, for this set, the single-click entanglement protocol outperforms the double-click protocol.

We did not find a noise threshold for the state-of-the-art parameter set—neither with existing GHZ protocols (see Sec. II A) nor when optimization over GHZ protocols (see Sec. II B). This finding is consistent with earlier investigations with a simplified NV center model.^{21} We identify two main limitations. The first one is the link efficiency: in this regime, the average entanglement generation times are longer than coherence times during entanglement generation—i.e., $ \eta link * < 1$. On top of that, the Bell pair fidelity is relatively low. A low Bell pair fidelity requires complex distillation protocols to achieve high-quality GHZ states. This, in turn, magnifies the impact of decoherence.

### B. Near-term parameter sets

As expected, further experimental progress and improved fidelities are required for fault-tolerant quantum computation. In the remainder of this section, we characterize two key parameters that drive the code performance in this regime. These findings can be used to guide future hardware development. Specifically, we investigate the effect of improving the Bell pair fidelity and the link efficiency.

#### 1. Sensitivity to Bell pair fidelity

First, we investigate the influence of the Bell pair fidelity by using a near-future setting parameter set—see the second column in Table I. Compared to the state-of-the-art parameter set of Sec. V A, in this set coherence times during entanglement creation and the photon detection probability are one and two orders of magnitude higher, respectively. The double-click entanglement protocol now gives rise to the best combination of entanglement success probability and Bell pair fidelity, as explained in more detail in Appendix A. This means that these near-future parameters allow for an increase in the link efficiency by four orders of magnitude compared to the state-of-the-art parameter set of Sec. V A—see Eqs. (1), (B2), and (B4).

In our Bell pair model, several parameters contribute to the infidelity of the Bell pair states similarly—i.e., through the parameter $\varphi $ of Eq. (B3) that captures all dephasing noise of the model. To investigate the sensitivity of the performance with respect to the Bell pair fidelity, we vary the influence of dephasing by scaling the probability of double excitation probability and off-resonant excitation errors. These are considered one of the leading error sources in present experiments.^{73} We show the results in Fig. 6. In this figure, the bottom of the horizontal axis indicates the Bell pair fidelity; the top indicates the corresponding excitation error probability.

We find $ p g = p m$ thresholds between $ p th = 0.0066 ( 25 )$% for $ F link \u2248 0.78$ and $ p th = 0.193 ( 4 ) %$ for $ F link \u2248 0.96$. Interestingly, the minimum fidelity for which we find a threshold, $ F link \u2248 0.78$, is lower than the state-of-the-art Bell pair fidelity demonstrated with both the single-click and double-click protocols.^{39,73} This is possible because the link efficiency allows performing several distillation steps.

We find different optimal protocols as a function of the Bell pair fidelity. In particular, we find that the optimal protocols require more distillation steps as we reduce the Bell pair fidelity, ranging from *k* = 12 for $ F link \u2248 0.78$ to *k* = 7 for $ F link \u2248 0.96$. We find lower thresholds as we decrease the Bell pair fidelity since the more complex distillation protocols amplify the effect of decoherence and require more gates. Furthermore, since existing GHZ creation protocols either have a small number ( $ k \u2264 8$) or many ( $ k \u2265 16$) distillation steps, we can understand why the new protocols with $ k \u2208 { 10 , 11 , 12}$ outperform them in this regime.

#### 2. Sensitivity to the link efficiency

Second, we investigate the influence of the link efficiency for near-future parameter values. In particular, we make use of a Bell pair fidelity $ F link \u2248 0.95$—close to the highest value in Sec. V B 1—and we investigate two options for varying the link efficiency.

First, we vary the link efficiency by varying the coherence times during entanglement generation. For this investigation, which we report in Fig. 7, we use the parameter set of the third column of Table I. In this set, we use a high photon detection probability $ \eta ph = 0.4472$, leading to an optical entanglement success probability of $ p link = 0.1$. The $ p g = p m$ threshold values vary between $ p th = 0.020 ( 8 )$% with coherence times corresponding to $ \eta link * = 4 \xd7 10 2$ and $ p th = 0.240 ( 9 )$% for coherence times corresponding to $ \eta link * = 2 \xd7 10 5$. For coherence times corresponding to $ \eta link * = 3 \xd7 10 2$ and lower, we did not find thresholds. At the other end of the spectrum, we evaluate the thresholds in an idealized modular scenario: in particular, in the absence of decoherence and with perfect SWAP gates (the last two points on the horizontal axis of Fig. 7). The last point corresponds to a scenario similar to the one analyzed in Ref. 34. We report a similar threshold value. For the Stringent protocol, the difference of $ p th = 0.775$% reported in Ref. 34 and $ p th = 0.601 ( 29 )$% found here can be attributed to the choice of the error-syndrome decoder and a reduced number of syndrome measurement cycles.

We now verify that, in this regime, similar thresholds can instead be found by varying the link efficiency via the entanglement generation rate. Specifically, we vary the entanglement success probability by adjusting the total photon detection probability $ \eta ph$. For this investigation, which we report in Fig. 8, we use the parameter set in the fourth column of Table I. This set contains coherence times during entanglement generation that are ten times higher than the state-of-the-art coherence times of Sec. V A.^{19} We find $ p g = p m$ thresholds between $ p th = 0.035 ( 4 )$% for a photon detection probability corresponding to $ \eta link * \u2248 4.7 \xd7 10 2$ and $ p th = 0.181 ( 4 )$% for a photon detection probability corresponding to $ \eta link * = 2 \xd7 10 3$. At photon detection probability corresponding to $ \eta link * \u2248 3.6 \xd7 10 2$ and lower, we are not able to find threshold values.

The second investigation gives rise to a similar required link efficiency ( $ \eta link * \u2248 4.7 \xd7 10 2$) as the first investigation ( $ \eta link \u2248 4 \xd7 10 2$). The small difference can be attributed to the slightly larger influence of the idling coherence time $ T 2 idle n$ in a scenario with a smaller entanglement rate. This shows that the link efficiency captures the key trade-off between cycle duration and decoherence rate, even when experimental overhead such as dynamical decoupling is accounted for.

We find that the parameter set used determines which GHZ protocol works the best. However, for a large range of parameters close to the state-of-the-art set, one protocol with *k* = 7 performs the best. We call this protocol Septimum and detail it in Fig. 2. In particular, this protocol is (one of) the best-performing protocol(s) at $ F link \u2248 0.96$ in Fig. 6, in the range $ 5 \xd7 10 2 \u2272 \eta link * \u2272 2 \xd7 10 3$ in Fig. 7, and in the range $ 6.3 \xd7 10 2 \u2272 \eta link * \u2272 1.1 \xd7 10 3$ in Fig. 8. We identify four additional well-performing protocols found with our dynamic program in Appendix E.

### C. GHZ cycle time sensitivity

In the following, we investigate the sensitivity of threshold values to the GHZ cycle time and the associated GHZ completion probability for our diamond defect center hardware model. We present the results in Fig. 9. In this figure, we see a clear dependence of the optimal GHZ completion probability on protocol complexity. In particular, protocols that take longer to finish (i.e., protocols with more distillation steps) peak at lower GHZ completion probabilities than those that finish faster, due to their increased susceptibility to decoherence. We see that for a protocol with relatively small *k*, GHZ cycle times that correspond to GHZ completion probabilities between 99.2% and 99.8% give rise to the highest threshold values in the parameter regimes considered here, whereas protocols with a large *k* peak at GHZ completion probabilities between approximately 92.5% and 98.5%.

We notice that, for some GHZ protocols, noise thresholds are found at relatively low GHZ completion probabilities of 90% and lower. This behavior can be directly attributed to the decoder heralding failures in the GHZ generation: as mentioned in Sec. III B, we utilize the stabilizer outcome from the previous time layer if a GHZ protocol does not finish within the GHZ cycle time, as opposed to naively performing the stabilizer measurement with the state produced by an unfinished GHZ protocol.

The results in Fig. 9 show that thresholds can strongly vary on the GHZ cycle time. For computational reasons, except for the results in this subsection, we do not optimize over the GHZ cycle time. Instead, we use a heuristic method to select this time based on the *k* value of each protocol. We describe this method in Sec. 5 of Appendix B.

## VI. DISCUSSION: FEASIBILITY OF PARAMETER SETS BELOW THRESHOLD

In Sec. V, we observed that the state-of-the-art parameter set is above the threshold. We identified two apparent drivers for this behavior: the Bell pair fidelity and the link efficiency. The sensitivity investigation shows that with a high link efficiency, the requirements on the Bell pair fidelity are modest, while even with a high Bell pair fidelity, a high link efficiency is still necessary.

Let us first discuss the experimental feasibility of the minimum link efficiency $ \eta link * \u2273 4 \xd7 10 2$ found in Fig. 8. First of all, the link efficiency can be increased by either increasing the coherence times of the data and memory qubits, or by increasing the entanglement success probability—or by a combination of both. In Sec. V, we found thresholds with a high success probability ( $ p link = 0.1$) and a modest increase in the coherence times. However, we also found that with high coherence times during entanglement generation (ten times higher than the state-of-the-art^{19}) and Bell pair fidelities of $ F link \u2248 0.95$, the total photon detection probability needs to fulfill $ \eta ph \u2273 0.19$ (Fig. 8). This is a factor 50 above the state-of-the-art parameter value.

The total photon detection probability is the product of multiple probabilities (see Sec. 1 of Appendix B). Present network experiments utilizing NV centers are particularly limited by two of these: the probability of emitting in the zero-phonon-line (ZPL, $ \u2248 3 %$)^{11} and the total collection efficiency ( $ \u2248 10 % \u2212 13 %$).^{19,39} Both values are expected to increase by Purcell enhancement of the NV center emission, for example with the use of optical microcavities^{13,18} or nanophotonic devices.^{74,75} However, even with such devices, the feasibility remains unclear. For microcavities, predicted ZPL emission and collection probabilities are of the order 10 to 46%^{13,18} and 49%,^{71} respectively. Moreover, the successful integration of Purcell-enhanced and optically coherent NV centers in nanophotonic devices remains an open research challenge due to the detrimental effects of surface charges.^{22}

This realization has led to an increased interest in other color centers in the diamond lattice, as, e.g., SiV and SnV defect centers.^{27} These centers have higher intrinsic emission probabilities into the ZPL—for SnV centers this is reportedly in the area of 60%,^{76} whereas SiV centers approximately emit 70 to 80% into the ZPL.^{77} Additionally, the inversion symmetry of SnV and SiV centers makes them less susceptible to proximal charges, facilitating integration into nanophotonic devices. Nanophotonic structures offer advantages over microcavities, such as stronger cooperativities enabled by the small mode volumes,^{23} and reduced sensitivity to vibrations of the cryostat hosting the emitter.^{18} A disadvantage of SnV and SiV centers over NV centers is the fact that they need to be operated at lower temperatures^{24} or under high strain^{78,79} to achieve similar coherence times.

Additionally, these alternative trajectories provide opportunities for “direct” GHZ generation schemes, where a GHZ state is created without Bell pair fusion.^{80} Contrary to the photon-emission-based Bell pair generation with NV centers, these direct GHZ state generation schemes could be based on the transmission or reflection of photons. Since, for nodes with a single communication qubit, SWAP gates are unavoidable when performing fusion, getting rid of SWAP gates during GHZ state generation could relax the requirements for, e.g., the link efficiency and the photon detection probability.

## VII. CONCLUSION

In this paper, we investigated the influence of decoherence and other noise sources on a fault-tolerant distributed quantum memory channel with the toric code. For this, we developed an open-source package that optimizes GHZ distillation for distributed stabilizer measurements and quantifies the impact of realistic noise sources.^{57} The GHZ protocols found with this package are compatible with a second open-source package that calculates logical error rates of the (distributed) surface code.^{63}

We focused our attention on a specific set of noise models inspired by diamond defect centers. We first observed that state-of-the-art nitrogen-vacancy center hardware does not yet satisfy the thresholds. A parameter-sensitivity analysis shows that the main driver of the performance is the link efficiency, giving a benchmark for future experimental efforts. The photon detection probability of state-of-the-art hardware appears to represent the main challenge for operating the surface code below threshold. Sufficient photon detection probabilities could be achieved with the help of Purcell enhancement of NV center emission, or using other color centers such as silicon-vacancy centers or tin-vacancy centers. The use of other color centers also presents opportunities for schemes that directly generate GHZ states between the communication qubits of more than two nodes—i.e., without fusing Bell pairs.^{80}

With our detailed noise models, we found threshold values up to 0.24%. This is three to four times lower than prior thresholds found with less-detailed models. Similarly, the optimal distillation protocols have a small number of distillation steps compared to prior work. For a large parameter regime of parameters, a protocol consuming a minimum of seven Bell pairs was optimal. Its experimental demonstration would be an important step for showing the feasibility of this approach for scalable quantum computation.

We performed a thorough optimization of GHZ distillation protocols. However, further improvements in other elements of the distributed architecture could partially bridge the gap with the performance of monolithic architectures. For instance, the surface code decoder could model unfinished GHZ distribution rounds as erasure noise. The conversion of binary trees to protocol recipes can be optimized and Bell pair distribution could be scheduled dynamically. On top of that, since our software allows for the implementation of general hardware models, further research could focus on analyzing and understanding a broad range of physical systems in the distributed context. In addition to exploring alternative hardware systems, it would be intriguing to implement a more in-depth model of the system's micro-architecture.

## ACKNOWLEDGMENTS

The authors would like to thank Hans Beukers, Johannes Borregaard, Kenneth Goodenough, Fenglei Gu, Ronald Hanson, Sophie Hermans, Stacey Jeffery, Yves van Montfort, Jaco Morits, Siddhant Singh, and Erwin van Zwet for helpful discussions and feedback. The authors gratefully acknowledge support from the joint research program “Modular quantum computers” by Fujitsu Limited and Delft University of Technology, co-funded by the Netherlands Enterprise Agency under project number PPS2007. This work was supported by the Netherlands Organization for Scientific Research (NWO/OCW), as part of the Quantum Software Consortium Program under Project 024.003.037/3368. This work was supported by the Netherlands Organisation for Scientific Research (NWO/OCW) through a Vidi grant. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 852410). D.E. was partially supported by the JST Moonshot R&D program under Grant JPMJMS226C. The authors thank SURF (www.surf.nl) for the support in using the National Supercomputer Snellius.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

Sébastian de Bone, David Elkouss and Paul Möller designed and conceptualized the study. Conor E. Bradley and Tim H. Taminiau designed, built, and conducted experiments to characterize color centers. Conor E. Bradley and Tim H. Taminiau curated the hardware parameter sets. Sébastian de Bone built the software to generate and evaluate GHZ protocols. Sébastian de Bone and Paul Möller built the software to carry out distributed surface code simulations. Sébastian de Bone carried out and analyzed the numerical simulations. The manuscript is written by Sébastian de Bone with input and revision from all authors. David Elkouss supervised the project.

**Sébastian de Bone:** Conceptualization (equal); Data curation (lead); Formal analysis (lead); Methodology (equal); Software (equal); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead). **Paul Möller:** Conceptualization (equal); Methodology (equal); Software (equal). **Conor E. Bradley:** Resources (lead); Writing – review & editing (supporting). **Tim H. Taminiau:** Resources (supporting); Supervision (supporting); Writing – review & editing (supporting). **David Elkouss:** Conceptualization (equal); Formal analysis (supporting); Supervision (lead); Writing – original draft (supporting); Writing – review & editing (supporting).

## DATA AVAILABILITY

The data that support the findings of this study are openly available in 4TU.ResearchData, Ref. 57.

### APPENDIX A: EXPERIMENTAL CHARACTERIZATION

The noise models and parameter values we use in this paper are based on experimental observations by Bradley *et al.*^{16,21} and Pompili *et al.*^{19} Experiments were performed on NV centers at temperatures of 3.7 and 4 K. Time scales for gates are based on microwave and radio frequency pulses for single-qubit gates on the electron and ^{13}C qubits, respectively. Time scales for two-qubit gates are based on phase-controlled radio frequency driving of the carbon spins interleaved with dynamical decoupling sequences on the electron state, following the scheme described by Bradley *et al.*^{12,16}

Characterization of the decoherence model and the associated coherence times was achieved by monitoring the drop in probability of detecting Pauli matrix eigenstates at *t* > 0 after preparing the state at *t* = 0. For the *T*_{1} relaxation time, the $ | 0 \u27e9$ and $ | 1 \u27e9$ states were used, and the expectation value $ \u27e8 Z \u27e9$ of a measurement in the Pauli-*Z* basis was determined after several delay times—the coherence time then directly follows from the observed exponential drop in the expectation value. With the $ | + \u27e9$ and $ | \u2212 \u27e9$ states and the expectation $ \u27e8 X \u27e9$ for measurement in the Pauli-*X* basis, and with $ | + i \u27e9$ and $ | \u2212 i \u27e9$ and $ \u27e8 Y \u27e9$, the *T*_{2} coherence time of our model could be determined. For these four states, the observed exponential decay $ T dec$ corresponds to the *T*_{2} time of our model via $ 1 / T 2 = 1 / T 1 \u2212 2 / T dec$.

The full model is based on the so-called NV^{−} state, a spin-1 electron qubit with spin projection quantum number $ m s \u2208 { \u2212 1 , 0 , 1}$. Stochastic ionization can convert the NV^{−} state to the NV^{0} state with $ m s \u2208 { \u2212 1 / 2 , 1 / 2}$. Under the present understanding, this spin state is no longer usable as a qubit due to a fast orbital relaxation process.^{81} Since the electron-spin state is used to control the ^{13}C spins, ionization accordingly dephases the ^{13}C states. As such, the coherence times of the NV center memory are currently limited by these ionization effects. Reference 21 proposes a method to mitigate ionization-induced memory dephasing by actively recharging the NV^{0} state. They show ionization and recharging can be performed with minimal loss in fidelity in the state of a ^{13}C spin in an isotopically purified device. This marks an important avenue for future research. In our model, we do not specifically include ionization, but simply absorb its influence in the ^{13}C coherence times.

Furthermore, Ref. 21 showed that isotopically engineered nitrogen-vacancy devices with a reduced ^{13}C memory qubit concentration (0.01%) are able to store a quantum state over 10^{5} entanglement attempt repetitions. Compared to samples with natural ^{13}C abundance (1.1%), the memory qubits have a weaker coupling to their color centers. While this increases coherence times during entanglement generation by several orders of magnitude, it also leads to longer time scales for carrying out gates on the memory qubits and gates between the communication and memory qubits. In natural abundance devices, a quantum state can typically be stored over 10^{3} entanglement attempt repetitions,^{19} while demonstrated single qubit ^{13}C gates are typically $ \u2248 13$ times faster and two-qubit gates are typically $ \u2248 50$ times faster than the isotopically purified samples.

Nonetheless, in this paper, we assume that the diamond color centers contain a natural abundance of carbon memory qubits (1.1%). Thus, we trade-off lower coherence times during entanglement generation for faster gates. This choice was made because it is believed that in future systems entanglement success rates are required to be several orders of magnitude higher than current state-of-the-art. In those regimes, fewer entanglement attempts are required and the influence of decoherence during entanglement generation becomes smaller. It is believed that one then benefits more from having the faster operations that samples with natural concentrations of ^{13}C atoms offer.

On top of that, Ref. 21 found an idling carbon coherence time of $ T 2 idle n = 10$ s, which is too low for running the GHZ creation protocols described in this paper when considering the corresponding gate speeds. This time—comparable to those achieved in natural abundance devices^{16}—was limited predominantly by other impurities in the diamond, but the expected linear scaling of coherence time with isotopic concentration remains to be demonstrated in future work. In regimes with high entanglement success rates, the time duration of the GHZ creation protocols is almost solely described by the duration of the two-qubit C*Z*, CNOT, C*iY*, and SWAP gates, which are $ \u2248 50$ slower for the isotopically purified samples. Thus, the operation time of any protocol also takes $ \u2248 50$ time longer with isotopically purified samples. Equivalent performance for isotopically purified samples would require that $ T 2 idle n$ also increases by a factor of 50—i.e., from 10 to 500 s.

Further research into isotopically engineered diamond defect centers is ongoing. This will likely lead to a better understanding of the trade-off between the ^{13}C concentration and the associated operation times and decoherence rates.

### APPENDIX B: SIMULATION MODELS AND SETTINGS

##### 1. Bell pair state

^{68}or double-click

^{69}Bell pair entanglement protocol. We denote the

*bright state*population coefficient for the single-click protocol as

*α*. Phase uncertainty originating from a path length difference between the two involved parties is modeled as a dephasing channel on one of the photonic states before the Bell state measurement, with fidelity

^{70,71}

*I*

_{0}and

*I*

_{1}are modified Bessel functions of the zeroth and first order, and $ \sigma ( \phi )$ is the standard deviation of the phase instability. We only include phase uncertainty for the single-click protocol, for which a phase instability of 14.3° corresponds to $ \lambda = 0.984$.

^{15}

The parameter $ \eta ph$ describes the total photon detection probability per excitation. It can be considered as the product of the total collection efficiency (the transmissivity between defect center and detector multiplied by the detector efficiency) with the probability that a photon is emitted in the detection (time) window and in the zero-phonon line. The photon detection probability only influences the success probability of the entanglement protocol.

We note that, in our model, we neglect photon detector *dark counts*. This is because we consider regimes in which dark counts are negligible.

The parameter $ p EE$ describes the probability of an excitation error during a heralded entanglement generation event. This can occur because an extra photon was emitted during the excitation pulse (a phenomenon known as *double excitation*) or as a result of exciting the dark state. We assume that these excitation errors give rise to a dephasing channel on one of the qubits of the Bell pair.

For NV centers, the double excitation probability is estimated between 4% and 7%.^{15,19} The off-resonant excitations are typically assumed to be negligible. This is because the polarization of the light pulse only leads to a weak driving field on transitions close to the main (bright state) transition, and other transitions are sufficiently far off-resonant.

For other systems, the situation might be different. Here, one can design the excitation pulse to induce a *π* transition on the main transition and a full $ 2 \pi $ rotation on the closest unwanted transition. Tiurev *et al.* created a model that, based on the energy difference between the main transition and the closest unwanted transition, allows one to estimate $ p EE$—i.e., both the double excitation probability and the probability of exciting this unwanted transition.^{82,83} Their model shows that, typically, the larger the energy difference is between the two transitions, the smaller $ p EE$ becomes.

*μ*is the photon indistinguishability per Bell pair measurement.

^{70}The parameter $\varphi $ contains all dephasing contributions of the model:

Our single-click model combines different elements from NV center single-click models for large-scale networks.^{15,19,37,70–72} Our model differs from these models in that we neglect dark counts in the photon detectors. More elaborate versions of the single-click and double-click models can be found in Refs. 84 and 85.

For the parameter value sets used in this paper, as presented in Table I of the main text, we use the single-click protocol to give an impression of the optimal Bell pair fidelity and success probability with state-of-the-art NV center parameter values. This is the parameter set identified as $ F prep = 0.99 , \u2009 p EE = 0.04 , \u2009 \mu = 0.9 , \u2009 \lambda = 0.984$, and $ \eta ph = 0.0046$. In the simulations performed with near-future parameter values, i.e., the set based on $ F prep = 0.999 , \u2009 p EE = 0.01 , \u2009 \mu = 0.95$, *λ* = 1, and $ \eta ph = 0.4472$, we use the double-click protocol model to generate the Bell pair states.

This is because, for the state-of-the-art parameter set, the single-click protocol is the best option. As can be seen in Eqs. (B2) and (B4), for a specific set of parameter values, the double-click protocol has a fixed Bell pair fidelity and success probability. For the single-click protocol, on the other hand, the bright space population parameter *α* allows one to trade in a higher fidelity for a lower success probability, and vice versa. This gives us the option to select *α* such that the success probability is the same for both protocols. If, for that value of *α*, the state generated by the single-click protocol is better than the state generated with the double-click protocol, one could argue that the single-click protocol is the best choice.

This is the case for the state-of-the-art parameter set, as we get, with our double-click model, $ p link ( dc ) \u2248 1 \xd7 10 \u2212 5$ and $ F link ( dc ) \u2248 0.852$. With our single-click model, using $ \alpha \u2248 0.001 \u2009 15$, we get a similar success probability, but a state with higher fidelity: $ F link ( sc ) \u2248 0.905$. We note that setting *α* this low is typically not possible in practical situations. However, with the state-of-the-art parameter set, also for higher values of *α* the single-click model produces better success probabilities and fidelities than the double-click model. In Table I, we have used a higher value for *α*, leading to one order of magnitude higher $ p link ( sc )$, and slightly lower fidelity.

For the near-future parameters, the double-click model becomes favorable with $ p link ( dc ) \u2248 0.1$ and $ F link ( dc ) \u2248 0.953$. Setting *α* to get the same success probability with single-click now only leads to $ F link ( sc ) \u2248 0.873$. Technically, for this parameter set, it is possible to reach slightly higher fidelities with single-click than with double-click, but this gives rise to very low success probabilities—i.e., success probabilities unusable for the coherence times considered in this paper.

##### 2. Decoherence

We describe decoherence noise channels $ N noise ( \rho ) = \u2211 i = 1 \kappa K ( i ) \rho ( K ( i ) ) \u2020$ on a general state *ρ* in terms of their *κ* Kraus operators $ { K ( i )} i = 1 \kappa $, with $ \u2211 i = 1 \kappa K ( i ) ( K ( i ) ) \u2020 = I$.

*t*is the time and

*T*

_{1}is the coherence time]:

*T*

_{2}is the coherence time]:

##### 3. Dynamical decoupling sequence length

*i*th and

*j*th attempt. Furthermore, $ t n \u2032 ( i , j ) = \u2308 max ( i , j ) / ( 2 n \u2032 ) \u2309 ( 2 n \u2032 t link + t pulse )$ is the effective time of performing the required entanglement attempts in this specific scenario, where $ \u2308 max ( i , j ) / ( 2 n \u2032 ) \u2309$ describes how many DD sequences are required for these attempts and $ 2 n \u2032 t link + t pulse$ describes the time of one DD sequence. To solve Eq. (B7) in a practical setting, it suffices to take a large number for $A$ instead of letting it go to infinity. Because finding $ n DD$ in this way only minimizes the waiting and refocusing time during entanglement generation in two nodes, and not the waiting time during the other operations, this process does not lead to the optimal $ n DD$. We found that it does, however, give rise to values for $ n DD$ that, typically, produce good results.

##### 4. Gate and measurement noise

*X*,

*Y*,

*Z*), the Hadamard gate, the CNOT (C

*X*) gate, the C

*Z*gate, and the C

*iY*gate. These are not the native gates of the NV center, but their true gate set can be compiled into the gate set used here without additional costs in terms of two-qubit gates.

^{20}Noise on two-qubit gates are modeled with a depolarizing noise model:

Measurements are restricted to measuring (single-qubit) electron qubits in the Pauli-*Z* basis. Measuring in the Pauli-*X* basis is achieved with an additional Hadamard gate. Measurement errors are modeled by a probability $ p m$ that the measurement projects onto the opposite eigenstate of the measured operator.

##### 5. GHZ cycle time settings

In Sec. V C, we discuss the influence of the GHZ cycle time $ t GHZ$ on the surface code threshold value. In this section, we discuss the heuristic method that we use for finding a suitable GHZ cycle time for a specific protocol at a specific set of error probabilities.

As discussed earlier, protocols with more distillation steps *k* take (on average) longer to finish. This means that, compared to a protocol with a smaller *k*, they require a longer $ t GHZ$ to reach the same GHZ completion probability $ p GHZ$. However, because decoherence plays a larger role at a higher $ t GHZ$, we typically see that the threshold values obtained with protocols with higher *k* peak at a lower GHZ completion probability. This becomes clear in Figs. 9 and 10, where we plot how the threshold values change when using different GHZ cycle times for four protocols with varying *k*.

In an attempt to limit calculation times during our GHZ protocol search, we made use of a heuristic-driven method to estimate the optimal $ t GHZ$ at a certain error probability configuration. Using the protocol-specific parameter *k*, we first identify an adequate GHZ completion probability as $ p GHZ ( aim ) = ( 100.2 \u2212 k / 10 ) %$. On top of that, we use prior knowledge for a rough estimate $ p th ( est )$ of the value of the threshold. We then determine the distribution of the protocol's duration at the error probability $ p th ( est )$, by running it without a GHZ cycle time. Finally, using this distribution, we determine $ t GHZ$ by selecting a time at which *at least* a fraction $ p GHZ ( aim )$ of the iterations will finish.

### APPENDIX C: DETAILS REGARDING THE CALCULATION OF THE SUPEROPERATOR

##### 1. Superoperator calculation

In this section, we describe how we calculate the superoperator that we use in the surface code simulations. Separately calculating a superoperator has the advantage that it breaks up the process of calculating GHZ state creation from the threshold simulations: this drastically decreases the complexity of the full calculation.

Following earlier work by Nickerson *et al.*, we assume that only Pauli errors occur on the data qubits during the toric code simulations.^{34,35} This simplifies the simulation, as every stabilizer measurement now deterministically measures either +1 or −1, and measurement results can be calculated by simply considering commutativity between Pauli errors and the stabilizer operators. In most situations, the stochastic Pauli error model can be considered a good approximation for coherent errors described by continuous rotations.^{86} On top of that, since the nuclear spin qubits (i.e., the memory qubits) of NV centers have no states to leak to, it is believed that a depolarizing channel (i.e., Pauli noise) is a good approximation for noise on these qubits.

Our characterization of the toric code stabilizer measurements is carried out with density matrix calculations that *do* include more general errors. To align these calculations with the toric code calculations themselves, the stabilizer measurement channel is twirled over the Pauli group.^{87–89} This makes sure the superoperators describing the channel only contain Pauli errors. Each superoperator is constructed via the channel's Choi state—i.e., by using the GHZ state created by the concerning protocol to non-locally perform the stabilizer measurement on half of the maximally entangled state.

*P*that lead to orthogonal states $ | \Psi \xb1 ( m ) \u27e9$, i.e., we only use

_{m}*P*that make sure we have

_{m}We see that this procedure gives us the probabilities required to construct the superoperator of the channel. If $ \rho Choi$ is constructed by preparing the post-measurement state according to a +1 measurement outcome, the coefficients $ p + ( m )$ give rise to the Pauli errors $ P m \u2208 E$ without a measurement error on the stabilizer measurement, whereas the coefficients $ p \u2212 ( m )$ describe Pauli errors *P _{m}* accompanied with a measurement error. If $ \rho Choi$ is constructed with a −1 measurement outcome, the role of $ p + ( m )$ and $ p \u2212 ( m )$ is inverted, but the parameter values themselves are the same.

The stabilizer fidelity is defined as the coefficient $ p s ( m )$ corresponding to $ P m = I \u2297 4 \u2297 I \u2297 4$ (i.e., no errors on the data qubits) and no stabilizer measurement error. In our search for well-performing GHZ creation protocols, a good reason for comparing two protocols by using the stabilizer fidelity over, e.g., the GHZ state fidelity is the fact that the surface code data qubits undergo more decoherence for protocols that take longer to finish. This aspect of the optimization problem is not taken into account if we just use the GHZ fidelity to compare the protocols.

##### 2. Convergence of the average superoperator

The construction of a superoperator that we use in the surface code simulator requires averaging over a large number of Monte Carlo simulations. In this section, we investigate the convergence of the average superoperator over an increasing number of Monte Carlo samples. In Fig. 11, we calculate the average Choi state $ \rho \xaf Choi$ over $ 3 \xd7 10 5$ Monte Carlo iterations and calculate the *trace distance* of this state with the average Choi state $ \rho \xaf Choi ( i )$ after a smaller number of *i* iterations. As explained in more detail below, this figure suggests that, after 10^{5} Monte Carlo samples, errors in the average superoperator elements are on the order of $ 10 \u2212 4$.

For $ s \u2208 { + , \u2212}$ and $ P m \u2208 E$, the superoperator used in threshold simulations is calculated as the average $ S \xaf = { p \xaf s ( m )} s , m$ of the individual superoperators $ S = { p s ( m )} s , m$ calculated with the method of Sec. 1 of Appendix C. Alternatively, this average superoperator can be constructed by calculating the average Choi state $ \rho \xaf Choi$, as defined in Sec. 1 of Appendix C, and using Eq. (C6) to calculate $ { p \xaf s ( m )} s , m$ from this average Choi state.

*ρ*and $ \rho \u2032$ is defined as

^{67}

*i*, $ T ( \rho \xaf Choi , \rho \xaf Choi ( i ) )$ provides an upper bound on the difference between the average superoperator $ S \xaf ( i )$ after iteration

*i*and the average superoperator $ S \xaf$ after the full number of iterations. This is because, for the difference in the elements of $ S \xaf$ and $ S \xaf ( i )$, the following holds:

*i*iterations.

In our simulations, as shown schematically in Fig. 5, the superoperator $ { p \xaf s ( m )} s , m$ is subsequently used in a second level of Monte Carlo simulations that emulates the operation of the surface code with this specific superoperator. In Appendix D, we discuss how, for a specific $ { p \xaf s ( m )} s , m$, the statistical uncertainty in the Monte Carlo simulations of the surface code leads to uncertainty in the calculated threshold value.

### APPENDIX D: FITTING PROCEDURE FOR DETERMINING THRESHOLDS

In this appendix, we describe the fitting procedure we use to calculate toric surface code thresholds as well as the associated uncertainties in the threshold values. To perform the fitting procedure described below, we make use of the optimize.curve_fit function of the SciPy package for Python.^{90}

##### 1. Regression model

Calculating threshold values involves varying one or more of the error probabilities. For each combination of error probabilities *p*, we calculate an average superoperator using the methods described in Appendix C. At a specific *p*, we then use Monte Carlo simulations to calculate the logical success rate *r* of the toric surface code for multiple lattice sizes *L*.

*p*,

_{i}*L*) by

_{i}*r*. We make use of $ n C$ to describe the total number of input combinations: $ { ( p i , L i )} i = 1 n C$. For a single (

_{i}*p*,

_{i}*L*) combination, the logical success rate is defined as the number of error-correction iterations

_{i}*M*that do not induce a logical error divided by the full number of error-correction iterations

_{i}*N*. In the context of this paper,

_{i}*N*can be considered as the number of Monte Carlo iterations used for surface code calculations for a certain (

_{i}*p*,

_{i}*L*) and the exact hardware configuration used. We assume that the uncertainty in the observed logical success rates is described by the binomial distribution. This means that the standard deviation can be estimated via

_{i}*et al.*, we fit the logical success rates with the following model:

^{91}

*estimates*$ { r i \u0302} i$ of the logical success rates for all input combinations $ { ( p i , L i )} i$. For a certain (

*p*,

_{i}*L*), the

_{i}*residual*$ \epsilon \u0302 i$ is defined as the difference between the observed logical success rate and the estimated value: $ \epsilon \u0302 i = r i \u2212 r \u0302 i$. Values for the seven fitting parameters $ a \u0302 , \u2009 b \u0302 , \u2009 c \u0302 , \u2009 d \u0302 , \u2009 p \u0302 th , \u2009 \kappa \u0302$, and $ \zeta \u0302$ are found by identifying their (local) minimum with respect to the sum

*Q*of the “weighted” squared residuals. This sum is defined in the following way:

*σ*) are given less priority in the least-squares fit than data points with low uncertainty.

_{i}##### 2. Weighted least-squares fitting procedure

To understand how the confidence intervals in the values of the fitting parameters are determined, we delve a bit deeper into how one could determine fitting parameters for a non-linear regression like Eq. (D2). In line with convention, we denote our input configuration as $ X i = ( p i , L i )$, and we write $ r \u0302 i$ as $ r \u0302 i = f ( \beta \u0302 , X i )$. Here, the function *f* is the function of Eq. (D2) and $ \beta \u0302 = ( a \u0302 , b \u0302 , c \u0302 , d \u0302 , p \u0302 th , \kappa \u0302 , \zeta \u0302 )$ describes the converged values for the fitting parameters after the optimization. We define $ n P = 7$ as the number of fitting parameters of the model. Furthermore, we use the notation $ \beta ( t )$ to indicate the fitting parameter values at a certain step *t* during the optimization. We can now write $ r i = f ( \beta ( t ) , X i ) + \epsilon i ( t )$ for each *t*. Here, the residuals also contain the superscript (*t*) to denote that they depend on the exact values of $ \beta ( t )$.

*Gauss–Newton algorithm*.

^{92,93}This method is a variant of Newton's method for finding the minimum of a non-linear function. We start with a guess $ \beta ( 1 )$ for the fitting parameter values. These values are then iteratively updated by using the fact that we want to minimize the parameter

*Q*of Eq. (D3) until they converge. To go from a certain $ \beta ( t )$ to a new improved version $ \beta ( t + 1 )$, one makes use of $ r i = f ( \beta ( t + 1 ) , X i ) + \epsilon i ( t + 1 )$ to construct a new estimation in terms of the old fitting parameter values $ \beta ( t )$. More explicitly, $ f ( \beta ( t + 1 ) , X i )$ is Taylor expanded around $ \beta ( t )$, and the second and higher order terms are neglected. This allows one to write $ r i = f ( \beta ( t + 1 ) , X i ) + \epsilon i ( t + 1 )$ as a matrix equation—i.e., with each of the $ n C$ input and output combinations of

*X*and

_{i}*r*in a separate row of this equation. To this end, we put the fitting parameter values at step

_{i}*t*of the optimization process into a $ n P \xd7 1$ column vector

**and $ \u03f5 ( t )$, respectively. Finally, we define $ \Delta \beta ( t + 1 )$ as $ \Delta \beta ( t + 1 ) \u2261 \beta ( t + 1 ) \u2212 \beta ( t )$. The full (Taylor expanded) version of $ r i = f ( \beta ( t + 1 ) , X i ) + \epsilon i ( t + 1 )$ can now be rewritten as**

*X*^{92,93}

*f*with respect to the fitting parameters, evaluated at $ \beta ( t )$ and the inputs $ { X i} i$. This matrix is also known as the

*Jacobian*matrix.

*Q*from Eq. (D3) can also be rewritten as a matrix product as follows:

*r*as $ \Sigma \u2261 diag ( \sigma 1 2 , \sigma 2 2 , \u2026 , \sigma n C 2 )$. More generally, one can use the full covariance matrix for $\Sigma $ if this information is available.

_{i}*Q*with respect to $ \Delta \beta ( t + 1 )$, one sets $ \u2202 Q / \u2202 \Delta \beta ( t + 1 )$ to zero. This results in an expression that can be solved for $ \Delta \beta ( t + 1 )$:

^{93}

##### 3. Uncertainty in fitting parameter values

*true*values, which we indicate with $\beta $. One can use another Taylor expansion to express $ \beta \u0302$ in terms of the true set of values $\beta $:

Strictly speaking, in the limit $ n C \u2192 \u221e$, we have $ \beta \u0302 = \beta $ and $ \beta \u0302$ does not have a distribution, since $\beta $ contains constant values. For finite $ n C$, however, we can argue that $ \beta \u0302$ *does* have a distribution, and we use the last term of Eq. (D8) to estimate the uncertainty in the fitting parameters $ \beta \u0302$. This estimation again involves the assumption that $ J \u0302$ is a constant matrix and does not depend on the fitting parameters.

**, the variance of**

*A***is given by $ Var ( AY ) = A \u2009 Var ( Y ) A T$. Together with the assumption that $ Var ( \u03f5 )$ can be estimated by $ Var ( \u03f5 ) = \Sigma $, the covariance matrix of $ \beta \u0302$ can be expressed as**

*AY*^{92,93}

*reduced chi-squared*metric, which is defined as

*p*values, we would typically find fits with $ \chi \nu 2 > 1$. In those situations, we scaled up each

*σ*with $ \chi \nu $, leading to

_{i}*t*-distribution:

*ν*and a confidence interval of $ I ci = 95$%, we have $ t ci \u2248 1.96$. Smaller values of

*ν*lead to $ t ci$ values that are slightly bigger. In Fig. 12, we see an example of a threshold plot with a 95% confidence interval in the value found for the threshold fitting parameter.

### Appendix E: Selection of best-performing GHZ generation protocols

At the end of Sec. V B 2, we discuss the Septimum protocol (depicted in Fig. 2): a GHZ generation protocol found with the dynamic program of Sec. II B that gives rise to the highest thresholds for the simulation parameters used in several segments of Figs. 6–8. In this appendix, we identify four additional GHZ generation protocols that perform the best in multiple segments of the figures in Sec. V B: the protocols Sextimum, Decimum, Undecum, and Duodecum. We depict the timed binary trees of these four protocols in Fig. 13 and provide more information on their performance in Table II.

. | k
. | q
. | Figure 6 . | Figure 7 . | Figure 8 . |
---|---|---|---|---|---|

Sextimum | 6 | 3 | $ [ 4 \xd7 10 2 , 5 \xd7 10 2 ]$ | $ 4.7 \xd7 10 2$ | |

Septimum | 7 | 3 | 0.96 | $ [ 5 \xd7 10 2 , 2 \xd7 10 3 ]$ | $ [ 6.3 \xd7 10 2 , 1.1 \xd7 10 3 ]$ |

Decimum | 10 | 3 | $ [ 0.9 , 0.93 ]$ | $ [ 2 \xd7 10 3 , 2 \xd7 10 5 ]$ | $ 2 \xd7 10 3$ |

Undecum | 11 | 3 | $ [ 0.85 , 0.9 ]$ | ||

Duodecum | 12 | 4 | $ [ 0.78 , 0.82 ]$ |

. | k
. | q
. | Figure 6 . | Figure 7 . | Figure 8 . |
---|---|---|---|---|---|

Sextimum | 6 | 3 | $ [ 4 \xd7 10 2 , 5 \xd7 10 2 ]$ | $ 4.7 \xd7 10 2$ | |

Septimum | 7 | 3 | 0.96 | $ [ 5 \xd7 10 2 , 2 \xd7 10 3 ]$ | $ [ 6.3 \xd7 10 2 , 1.1 \xd7 10 3 ]$ |

Decimum | 10 | 3 | $ [ 0.9 , 0.93 ]$ | $ [ 2 \xd7 10 3 , 2 \xd7 10 5 ]$ | $ 2 \xd7 10 3$ |

Undecum | 11 | 3 | $ [ 0.85 , 0.9 ]$ | ||

Duodecum | 12 | 4 | $ [ 0.78 , 0.82 ]$ |

## REFERENCES

*Lectures on Topological Codes and Quantum Computation*

*Quantum Information Theory*

*Quantum Computation and Quantum Information*

*Advanced Data Analysis & Modelling in Chemical Engineering*