Thresholds for the distributed surface code in the presence of memory decoherence Open Access

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.³ 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.¹⁰ 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.¹⁶ 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%.³³ 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%.³⁴ 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³⁴ 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 $η link *$⁠.²¹ 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 $η link * ≈ 4 × 10 2$⁠. This improves on the prior results of $η link * = 2 × 10 5$ for the idealized time scale estimates of Nickerson et al.³⁴ However, this link efficiency is still above the state-of-the-art hardware²¹ reaching up to $η link * ≈ 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.

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.

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:²¹ 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.

On top of that, we consider five GHZ protocols found by Nickerson et al. in the context of distributed implementations of the toric code: Expedient (k = 22, q = 3) and Stringent (k = 42, q = 3) from Ref. 34, and Basic (k = 8, q = 3), Medium (k = 16, q = 4), and Refined (k = 40, q = 5) from Ref. 35.

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 ⟩ = ( | 0 ⟩ ⊗ n + | 1 ⟩ ⊗ n ) / 2$⁠, where $n ∈ { 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 ⟩$ 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.⁵⁴ 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 ⟩$ can be fused with a state $| GHZ n 2 ⟩$ 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 − 1 ⟩$⁠. 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.

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⁵⁸ that does not directly fit into this stabilizer distillation framework.^59–61 

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

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⁶⁴ 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.⁶⁵ 

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,³⁴ 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 ¯ success$ for iterations where the GHZ state is created within $t GHZ$⁠, and an unsuccessful superoperator $S ¯ 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.

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.

In our model, qubits undergo both generalized amplitude damping and phase damping noise,⁶⁷ with separate T₁ and T₂ times. The generalized amplitude damping channel decoheres to the maximally mixed state $( | 0 ⟩ ⟨ 0 | + | 1 ⟩ ⟨ 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.¹⁶ Here, the coherence times that we consider are also those achieved for NV center spin registers with dynamical decoupling.¹⁶ 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)⁶⁸ entanglement protocol and for the double-click (i.e., the two-photon)⁶⁹ 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⁷⁰), and the total photon detection probability $η 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 , p EE$⁠, and μ affect the fidelity $F link$ of the Bell pair state of the double-click protocol, whereas the parameter $η 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 $η 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$⁠.

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.

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 ¹³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.²¹ 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., $η 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.

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.

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 $ϕ$ 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.⁷³ 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 ≈ 0.78$ and $p th = 0.193 ( 4 ) %$ for $F link ≈ 0.96$⁠. Interestingly, the minimum fidelity for which we find a threshold, $F link ≈ 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 ≈ 0.78$ to k = 7 for $F link ≈ 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 ≤ 8$⁠) or many (⁠ $k ≥ 16$⁠) distillation steps, we can understand why the new protocols with $k ∈ { 10 , 11 , 12 }$ outperform them in this regime.

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 ≈ 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 $η 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 $η link * = 4 × 10 2$ and $p th = 0.240 ( 9 )$% for coherence times corresponding to $η link * = 2 × 10 5$⁠. For coherence times corresponding to $η link * = 3 × 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 $η 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.¹⁹ We find $p g = p m$ thresholds between $p th = 0.035 ( 4 )$% for a photon detection probability corresponding to $η link * ≈ 4.7 × 10 2$ and $p th = 0.181 ( 4 )$% for a photon detection probability corresponding to $η link * = 2 × 10 3$⁠. At photon detection probability corresponding to $η link * ≈ 3.6 × 10 2$ and lower, we are not able to find threshold values.

The second investigation gives rise to a similar required link efficiency (⁠ $η link * ≈ 4.7 × 10 2$⁠) as the first investigation (⁠ $η link ≈ 4 × 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 ≈ 0.96$ in Fig. 6, in the range $5 × 10 2 ≲ η link * ≲ 2 × 10 3$ in Fig. 7, and in the range $6.3 × 10 2 ≲ η link * ≲ 1.1 × 10 3$ in Fig. 8. We identify four additional well-performing protocols found with our dynamic program in  Appendix E.

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.

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 $η link * ≳ 4 × 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¹⁹) and Bell pair fidelities of $F link ≈ 0.95$⁠, the total photon detection probability needs to fulfill $η ph ≳ 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, $≈ 3 %$⁠)¹¹ and the total collection efficiency (⁠ $≈ 10 % − 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%,⁷¹ 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.²² 

This realization has led to an increased interest in other color centers in the diamond lattice, as, e.g., SiV and SnV defect centers.²⁷ These centers have higher intrinsic emission probabilities into the ZPL—for SnV centers this is reportedly in the area of 60%,⁷⁶ whereas SiV centers approximately emit 70 to 80% into the ZPL.⁷⁷ 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,²³ and reduced sensitivity to vibrations of the cryostat hosting the emitter.¹⁸ A disadvantage of SnV and SiV centers over NV centers is the fact that they need to be operated at lower temperatures²⁴ 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.⁸⁰ 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.

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.⁵⁷ 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.⁶³ 

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.⁸⁰ 

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.

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.

The authors have no conflicts to disclose.

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).

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

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.¹⁹ 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 ¹³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₁ relaxation time, the $| 0 ⟩$ and $| 1 ⟩$ states were used, and the expectation value $⟨ Z ⟩$ 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 $| + ⟩$ and $| − ⟩$ states and the expectation $⟨ X ⟩$ for measurement in the Pauli-X basis, and with $| + i ⟩$ and $| − i ⟩$ and $⟨ Y ⟩$⁠, the T₂ coherence time of our model could be determined. For these four states, the observed exponential decay $T dec$ corresponds to the T₂ time of our model via $1 / T 2 = 1 / T 1 − 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 ∈ { − 1 , 0 , 1 }$⁠. Stochastic ionization can convert the NV⁻ state to the NV⁰ state with $m s ∈ { − 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.⁸¹ Since the electron-spin state is used to control the ¹³C spins, ionization accordingly dephases the ¹³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⁰ state. They show ionization and recharging can be performed with minimal loss in fidelity in the state of a ¹³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 ¹³C coherence times.

Furthermore, Ref. 21 showed that isotopically engineered nitrogen-vacancy devices with a reduced ¹³C memory qubit concentration (0.01%) are able to store a quantum state over 10⁵ entanglement attempt repetitions. Compared to samples with natural ¹³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³ entanglement attempt repetitions,¹⁹ while demonstrated single qubit ¹³C gates are typically $≈ 13$ times faster and two-qubit gates are typically $≈ 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 ¹³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¹⁶—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 CZ, CNOT, CiY, and SWAP gates, which are $≈ 50$ slower for the isotopically purified samples. Thus, the operation time of any protocol also takes $≈ 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 ¹³C concentration and the associated operation times and decoherence rates.

The parameter $η 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 π$ 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.

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 , p EE = 0.04 , μ = 0.9 , λ = 0.984$⁠, and $η ph = 0.0046$⁠. In the simulations performed with near-future parameter values, i.e., the set based on $F prep = 0.999 , p EE = 0.01 , μ = 0.95$⁠, λ = 1, and $η 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 ) ≈ 1 × 10 − 5$ and $F link ( dc ) ≈ 0.852$⁠. With our single-click model, using $α ≈ 0.001 15$⁠, we get a similar success probability, but a state with higher fidelity: $F link ( sc ) ≈ 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 ) ≈ 0.1$ and $F link ( dc ) ≈ 0.953$⁠. Setting α to get the same success probability with single-click now only leads to $F link ( sc ) ≈ 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.

We describe decoherence noise channels $N noise ( ρ ) = ∑ i = 1 κ K ( i ) ρ ( K ( i ) ) †$ on a general state ρ in terms of their κ Kraus operators ${ K ( i ) } i = 1 κ$⁠, with $∑ i = 1 κ K ( i ) ( K ( i ) ) † = I$⁠.

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.

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 − 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.

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.⁸⁶ 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.