Introduction

The past years have seen rapid advances in controlling small groups of qubits encoded in atomic or atom-like quantum memories. An important question now concerns the development of architectures to efficiently combine these memories into large-scale systems capable of general-purpose quantum computing,1,2,3 quantum simulation,4 and metrology near the quantum limit.5 A promising approach is entangling the atomic qubits with optical links to generate cluster states, which can perform general-purpose quantum computing with adaptive single-qubit measurements.6 A key challenge is to produce this cluster state fast enough to allow computation and error correction within a coherence time of memories.

Here, we show that percolation of heralded entanglement allows us to create arbitrarily large cluster states. Our cluster state-based architecture does not require reconfiguration of physical qubits by the result of the entanglement success unlike conventional approaches.7 Instead, the optical network reconfigures the connectivity. Meanwhile, the concept of percolation greatly relaxes the high success probability required in the previous proposals.1,2,3 The process is fast enough for implementation with the device parameters demonstrated; one does not need high cooperativity cavities and ancilla single photons. It requires ‘feed-forward’, operations conditioned on the previous measurement results. In the absence of errors, the missing bonds can be compensated with constant overhead.8,9,10,11 In the presence of errors, the scheme can be adapted to be fault-tolerant.12,13 Our approach also provides tolerance for site imperfections, and we show trade-offs between reaching percolation threshold and site-loss-resilience. Combined with transparent nodes implemented on nanophotonic platform, long-range connections are possible, reducing the percolation threshold. We derive a theoretical lower bound on the minimum time required to percolate in any lattice, and found that the proposed lattices are only a factor of \(1.6 \sim 3\) above this limit.

We focus on nitrogen vacancy (NV) centers in diamond14 and propose the control sequence to map the physical properties to cluster state quantum computation. NV centers have favorable properties as quantum memories. The NV\({}^{-}\) charge state has a robust optical transition for heralded entanglement between remote NV centers15,16 and a long electronic spin (\(S\)=1) coherence time approaching 1 second.17 Recently, single-qubit gates with fidelities up to 99\(\%\) were achieved with optimal control techniques.18 Electronic spins of NV centers can be coupled with nearby nuclear spins, which have coherence times exceeding 1 minute.19 In addition, they can be coupled with integrated nanophotonic devices.20

Our approach possibly applies to a number of physical systems, including atomic gases,21 ion traps,22 semiconductor quantum dots,23 and rare earth ions.24

Results

Protocol

Figure 1 illustrates the percolation approach to generate cluster states with NV centers. We work in the framework of cluster states where nodes represent qubits in the superposition state \((\left|0\right\rangle +\left|1\right\rangle )/\sqrt{2}\) and bonds represent the controlled-Z gate. Consider a square lattice where edges exist with probability \(p\) (Fig. 1(a–c)). The computational power of the cluster state corresponding to the graph is related to the size of the largest connected component (LCC) (shown in red). When \(p\;<\;0.5\), the graph produces small disconnected islands. For a lattice with \(N\) nodes, the size of the LCC is \(O(\mathrm{log}(N))\).25 Single-qubit measurements on such a cluster state can be efficiently simulated on classical computers. When the bond probability \(p\) exceeds the critical value \({p}_{c}=0.5\), called the percolation threshold of the lattice, there is a sudden transition in the size of the LCC. On crossing this boundary, the size of the LCC jumps to \(\Theta (N)\) and the lattice is ‘percolated’. A large percolated lattice can be ‘renormalized’ with single-qubit measurements to obtain a perfect lattice.8 Renormalization consumes a constant fraction of qubits8 and requires classical computation (pathfinding).26,27 Resulting perfect lattice is a resource for universal quantum computation. Thus, the percolation transition is accompanied by a sudden transition in computational power; adaptive single-qubit measurements on the cluster state enable universal quantum computing.11

Fig. 1
figure 1

Cluster state generation by percolation. a, b Transition in the size of the largest connected component (LCC) with increasing bond probability. Spheres and lines represent nodes and bonds, respectively, and the red spheres are in LCC. When the bond probability (\(p\)) goes above the percolation threshold (\({p}_{c}\)), the size of the LCC suddenly increases. Corresponding graph states change from being classically simulable to a resource for universal quantum computation. c Expanded view of a. d Physical implementation of nodes and bonds with NV centers in diamond. Probabilistic Bell measurement is attempted on two nearest neighbor electronic spins (blue spheres). Conditioned on two single-photon detection events, the two electron spins are entangled (Bell states). Hyperfine interaction entangles electronic spins and nuclear spins (\({}^{15}\)N) through controlled-Z gates. Measurement of electronic spins in transverse direction projects remote nuclear spins onto an entangled state (entanglement swapping).

Figure 1d shows the physical implementation of bond creation (entanglement) with NVs. The nuclear spins (red spheres) are ‘client qubits’ that store entanglement. They are coupled to NV electronic spins (‘broker qubits’), that can be entangled remotely by photon-mediated Bell measurements. In each time step, we attempt to create one bond (entanglement) at each node. We consider the Barrett–Kok protocol28 on the broker qubits of neighboring nodes. If the probabilistic Bell measurement succeeds, the electron spins in each node are entangled. This entanglement is transferred to the nuclear spins with the entanglement swapping procedure, as illustrated in Fig. 1(d)1. The whole cycle from initialization to entanglement swapping is assumed to be approximately t0 = 5 µs based on experimental demonstrations16 (See Supplemental Material, which includes a discussion of emitters coupled with nanophotonic circuits,29 ultrasmall mode volume cavities,30,31 waveguide-integrated single-photon detectors,32 low jitter detectors,33 effects of frequency mismatch between two emitters,34,35 spectral diffusion,36,37 nitrogen implantation,38,39 and charge state initialization40,41).

If the Bell measurement fails, one needs to restore the nuclear spins as before the trial. In addition, the electronic spins should be initialized to spin ground state without disrupting the nuclear spin. For the first problem, we just wait for the nuclear spin and electronic spin to be decoupled, which happens after a time period of the hyperfine interaction. This time period serves as a global clock and can be synchronized across the whole graph only using the \({}^{15}\)N nuclear spin, which is the host atom of NV centers. For the electron spin initialization, optical pumping technique cannot be used as in many experiments,15,16 because hyperfine coupling persists during a random amount of time when the electronic spin is in ms = 1 states. Instead, we readout the spin-state via ms = 0 optical transition where the electronic spin decouples from the nuclear spin. If the spin is in the other state, one can use a fast microwave pulse to put it in the spin ground state (See Supplemental Material).

We assume the Barrett–Kok protocol because it does not require ancilla single photons or high cooperativity cavities, unlike reflection-based protocols.1 Furthermore, unlike single-photon detection protocols,42 photon loss does not degrade the entanglement fidelity \(\left\langle \Psi \right|\hat{\rho }\left|\Psi \right\rangle\) where \(\Psi\) is the Bell state, and \(\hat{\rho }\) is density matrix operator of the system after successful heralding. High entanglement fidelity is important for reducing the error correction overhead in a fault-tolerant architecture. High fidelity comes at the price of low bond success probability,28 that can be overcome in the percolation based architecture.

Physical unit

The physical unit cell could be very small, on the order of tens of microns, so that the entire lattice may be integrated on a chip, as illustrated in Fig. 2. Each node in the architecture requires an atomic memory and a \(1\times d\) switch, where \(d\) is the number of nearest neighbors in the underlying lattice, i.e., the degree of the lattice. Each bond in the lattice requires waveguides, a beam-splitter, and two detectors for the Bell measurement.

Fig. 2
figure 2

Physical implementation of the proposed architecture. A unit cell consists of an atomic memory, a \(1\times 4\) switch, a 50/50 beam-splitter, waveguides and four single-photon detectors. Single photons emitted from the atomic memory are coupled to the waveguide and directed to the switch. The switch chooses one of the nearest neighbor nodes to be entangled with, and single photons are interfered using a 50/50 beam-splitter. Single-photon detectors detect interfered photons leaving electronic spins entangled.

At each time step, a state-selective optical \(\pi\)-pulse entangles a photonic mode with the electronic spin. Photonic modes coupled to neighboring electronic spins undergo probabilistic linear optic Bell measurement. Successful Bell measurement corresponds to the creation of a bond in the lattice.

Let us define \({n}_{{\rm{a}}tt}\) as the total number of entanglement trials per node. Then each bond is attempted \({n}_{{\rm{a}}tt}/d\) times. If a bond is created before \({n}_{{\rm{a}}tt}/d\) allocated trials, the node idles on the rest of the attempts for that bond. Each switch only needs to be flipped \(d-1\) times during cluster generation, and hence the switching time is negligible compared with the whole process. For example, electro-optic modulators can switch at sub-nanosecond time scales, but the time spent on each bond is on the order of milliseconds as we will see.

The probability of successfully heralding entanglement of two NV centers is \({p}_{0}={\eta }^{2}/2\),28 where \(\eta\) is the efficiency of emitting, transmitting, and detecting the photon entangled with the electronic spin (zero phonon line) from the NV-excited state. Table 1 summarizes \({p}_{0}\) for three representative types of emitter-photon interfaces: low-efficiency interfaces with \({p}_{0}=5\times 1{0}^{-5}\) representative of today’s state-of the art circular gratings or solid immersion lenses (SILs),15,16,43 medium-efficiency interfaces with \({p}_{0}=2\times 1{0}^{-4}\) for NV centers coupled to diamond waveguides,20 and high-efficiency \({p}_{0}=5\times 1{0}^{-2}\) for nanocavity-coupled NV centers.44 For all three coupling mechanisms, we assumed coupling efficiencies that are realistic today (See Supplemental Material). After \({n}_{{\rm{a}}tt}/d\) entanglement attempts with a nearest neighbor, the probability of having generated a bond is \(p=1-{(1-{p}_{0})}^{{n}_{{\rm{a}}tt}/d}\).

Table 1 Bell measurement success probability (\({p}_{0}\)), Bond trial time (\({t}_{0}\)), and readout time for three different coupling schemes.

Percolation thresholds

We evaluate the growth of clusters using the Newman-Ziff algorithm with nine million nodes. Figure 3a plots the fraction of nodes that are in the largest cluster component (\({f}_{{\rm{L}}CC}\)), as a function of time from the start of the protocol for the three values of \({p}_{0}\), assuming t0 = 5 µs, with an underlying square lattice. Initially, \({f}_{{\rm{L}}CC}\) is O(\(\mathrm{log}(N)/N\))25 where \(N\) is the total number of nodes in the lattice. As the bond success probability passes the bond percolation threshold (\({p}_{c}\)), \({f}_{{\rm{L}}CC}\) rapidly rises approaching \(\Theta (1)\). For a degree \(d\) lattice, the bond probability after time \(t\) is \(p=1-{(1-{p}_{0})}^{t/{t}_{0}d}\). The time required to obtain a resource for universal quantum computation is \({t}_{c}={t}_{0}d\mathrm{ln}(1-{p}_{c})/\mathrm{ln}(1-{p}_{0})\), which is shown with vertical dashed lines in the figure. The transition becomes sharper as the number of nodes in the lattice (\(N\)) increases.

Fig. 3
figure 3

Size of the largest connected component vs time for a different values of \({p}_{0}\), the Bell measurement success probability in one attempt and b different underlying lattice geometries. A square lattice is used in a and \({p}_{0}=0.02\, \%\) is used for b.

In all collection schemes, the bond success probability exceeds the percolation threshold within 1 second, which is a conservative estimate of the nuclear spin coherence time in nanostructured, integrated systems based on the experimentally demonstrated coherence time, 1 minute.19 These simulations reveal that an arbitrarily large cluster can be generated even with free space couplings.

Degree of the lattice and imperfection

It is known that higher degree lattices have a lower percolation threshold. However, \({t}_{c}\) is nearly the same for the honeycomb (\(d=3\)), the square (\(d=4\)) and the triangular (\(d=6\)) lattices (Fig. 3b). This is because increasing \(d\) lowers the bond percolation threshold, but it also decreases the number of entanglement attempts per bond, \({n}_{{\rm{a}}tt}/d\); a single broker qubit per NV requires entanglement attempts to proceed serially. Increasing \(d\) would in fact substantially lower \({t}_{c}\) if each site contained multiple broker qubits that could be entangled simultaneously. We do not explore the possibility of multi-NV nodes owing to a lack of studies.45 Ion trap systems are presently more mature for this purpose.2

Let’s consider the most general case where we can attempt Bell measurements on any pair of NVs at any time step. What is the minimum time, \({t}_{c}^{(LB)}\), required to obtain a resource for universal quantum computation over all lattice geometries? The bond probability after time \(t\) is \(p=1-{(1-{p}_{0})}^{t/{t}_{0}d}\). For percolation, \(p\ge {p}_{c}\), i.e., \(t\ge {t}_{0}d\cdot \mathrm{ln}(1-{p}_{c})/\mathrm{ln}(1-{p}_{0})\). For a degree \(d\) lattice, \({p}_{c}\ge 1/(d-1)\),46 with equality for a degree-\(d\) Bethe lattice (infinite tree with fixed degree at each node). This leads to \({t}_{c}\ge {t}_{0}d\cdot \mathrm{ln}(1-1/(d-1))/\mathrm{ln}(1-{p}_{0})\). \({t}_{c}\) is minimized as \(d\to \infty\) in which case we obtain \({t}_{c}^{(LB)}=-{t}_{0}/\mathrm{ln}(1-{p}_{0})\). \({t}_{c}^{(LB)}\) is plotted as a black dashed line in Fig. 3b. The lattice corresponding to the lower bound is the infinite-degree Bethe lattice, and such lattices are not a resource for universal quantum computing.47 Meanwhile, we find that the simple 2D lattices with nearest neighbor connectivity are only a factor \(\,3\) above this limit and are resources for universal quantum computing.

Practically, a scalable architecture should be able to tolerate non-functional sites. For example, trapping in a metastable state, a far-detuned transition, failed charge initialization greatly reduced the system performance in the recent state-of-the art experiment.48 Even if all faulty nodes and the bonds connected to it are removed, the lattice can retain enough bonds for a percolated lattice (Fig. 4a, inset). The problem maps to site-bond percolation. We define the site yield \(q\) as the fraction of functional nodes. Figure 4a plots the minimum time required to obtain a percolated lattice as a function of \(q\), assuming NVs coupled to diamond waveguides (\({p}_{0}=0.02 \%\)). In general, a reduced site yield can be compensated with a larger bond probability, which would require a longer time (more attempts). The site percolation threshold, \({q}_{c}\), corresponds to the minimum possible site yield for percolation with all bonds having succeeded (\(p=1\)). The triangular (\({q}_{c}=0.5\)) performs better than the square lattice (\({q}_{c}\approx 0.593\)) that outperforms the honeycomb lattice (\({q}_{c}\approx 0.697\)).

Fig. 4
figure 4

a The minimum time required to obtain a percolated lattice with sub-unity site yield. b \({f}_{{\rm{L}}CC}\) as a function of time for different values of \(\epsilon\) in the transparent node architecture. The insets show the bonds that can be attempted in a square lattice if the sites marked with crosses are inactive/transparent. \({p}_{0}=0.02 \%\).

The architecture that we have discussed thus far only allows for nearest neighbor interaction. Adding long-range connections shown in the inset of Fig. 4b can decrease the threshold time and increase tolerance to imperfect site yield. Such an architecture can be implemented by replacing the \(1\times 4\) switch in Fig. 2 with a \(5\times 5\) optical switch, which is depicted in Supplementary Fig. S1. Seven Mach-Zehnder interferometers (MZIs) with two phase shifters can implement switching between the set of input and output modes (Fig. S1c). The MZI arrays allow a node to be ‘transparent’. If a node is transparent, photons pass through the node to the next adjacent node for interference (Fig. 4b inset). By turning multiple adjacent nodes transparent, this bypass enables long-range entanglement in a planar architecture.

Interestingly, this transparent node architecture can be used to reduce percolation threshold by randomly turning a fraction \(1-\epsilon\) of the nodes transparent. Figure 4b plots \({f}_{{\rm{L}}CC}\) vs. time. As \(\epsilon\) decreases, the maximum possible value of \({f}_{{\rm{L}}CC}\) is reduced from one to \(\epsilon\) because only a fraction \(\epsilon\) nodes are active. However, reduced \(\epsilon\) decreases \({t}_{c}\) because the transparent nodes increase the effective dimensionality of the lattice. Therefore, there is an optimum value of \(\epsilon\) that maximizes LCC for a given time. We numerically found a minimum \({p}_{c}\) of \(0.33\) with transparent nodes, which is achieved when \(1/N\ll \epsilon \ll 1\), i.e., \(\epsilon \to 0\) but the number of non-transparent nodes in the lattice is still \(\Theta (N)\). Faulty sites can be incorporated into the fraction of transparent nodes as long as the yield exceeds \(1/N\).

Small size system

We used nine million qubits (nodes) in the simulations above in order to show the limiting behavior of the cluster creation. However, the same qualitative behavior is also observed in smaller systems. In Fig. 5a, we plot the mean (line) and the standard deviation (shade) of \({f}_{LCC}\) as a function of time for a \(5\times 5\) (green) and \(10\times 10\) (red) square lattices (\({p}_{0}=0.02 \%\)). In the simulation, we generated 300 lattice instances at each time with a periodic boundary condition and calculated the statistical mean and standard deviation. The transition from small disconnected islands to a large cluster on the order of the lattice size becomes sharper as the size of the system increases, but even for a \(5\times 5\) qubit system, there is a clearly visible transition near the percolation threshold (dashed line). Because of the small system size, there is statistical variation in \({f}_{LCC}\) and the shaded regions represent one standard deviation. As expected, the relative variation become smaller in larger systems.

Fig. 5
figure 5

a Size of LCC as a function of time in a \(5\times 5\) (green) and \(10\times 10\) (red) square grid with \({p}_{0}=0.02 \%\) (corresponding to NVs coupled to diamond waveguides). The lines represent the average size of the LCC and the shaded regions represent one standard deviation. The dashed line is the percolation threshold. b The combination of bond probability and measurement error probability achievable for different values of \({p}_{0}\). In order to achieve fault tolerance, the curve should have a portion inside the shaded regions which correspond to adaptive and non-adaptive measurement schemes from ref. 51.

Fault tolerance

It is possible to obtain a regular lattice from the percolated lattices in Fig. 3 with a constant overhead by finding crossing paths and using single-qubit measurements (renormalization).8,9,10,11 The single-qubit measurements used to obtain the regular lattice may bring additional errors, but we can adapt our architecture to the Raussendorf lattice.12,13 The Raussendorf lattice is a 3D lattice with degree-4 (Fig. S3a) and a means of translating surface code error correction into the cluster state model of quantum computation. The Raussendorf lattice can be constructed in a 2+1D architecture49 where qubits are arranged in 2D, and an additional dimension is constructed in time (Fig. S3 b)1. Then, one needs only two layers of 2D lattices to implement a 2+1D lattice because memories are reused after measurement (Fig. S3c). This architecture requires a small number of waveguide-crossings, which can be implemented with low loss.50 Alternatively, they can be absorbed in the optical switches.

Recent results have evaluated the fault tolerance in the Raussendorf lattice with non-deterministic entangling gates.26,51 Following,51 fault tolerance requires the bond probability (\(p\)) and measurement error to be in the two shaded regions in Fig. 5b, which correspond to adaptive and non-adaptive measurement schemes; in the non-adaptive scheme, qubits are measured in the \(X\)-basis regardless of the failure of bonds; in the adaptive scheme, one of the qubits connected to the missing bond is measured in the \(Z\)-basis transforming bond-loss to site-loss. The error threshold of each scheme vanishes at \(p \sim 0.935\) and \(p \sim 0.855\), respectively. These correspond to the site percolation threshold of the Raussendorf lattice (\({q}_{c} \sim 0.75\)) when lost bonds are transformed to missing sites (see more details in ref. 51).

We assume that the single-qubit measurement error probability increases with time as \(1-{e}^{-t/{t}_{coh}}\), where we use \({t}_{coh}=1\) sec. At \(t=0\), the measurement error probability and bond probability (\(p\)) are both zero. Both of these probabilities increase as we spend more time attempting entanglement generation, resulting in the curves shown in Fig. 5b depending on the entanglement success probability per attempt (\({p}_{0}\)). Assuming a bond trial time \({t}_{0}=5\)\(\mu\)s (Table 1), we find that a value of \({p}_{0}\approx 4\times 1{0}^{-3}\) is required to meet the fault tolerance threshold, which corresponds to an NV to detector coupling efficiency of \(\sim 9\)\(\%\).

As an example, we estimated the resource overhead for fault-tolerant factorization of 2000-bit numbers using Shor’s algorithm (Fig. 6). Various methods from refs. 52,53,54,55,56,57,58,59,60,61,62,63,64 are used for the calculation. Reference 54 in particular, has the least resource overhead even with the additional overhead from space-time layout and routing. With the assumptions on the error model described in the Supplementary Discussion, 2000-bit Shor’s algorithm requires 2.9 billion (280 million) qubits and 2.2 hours (43 mins) of computation given a physical error rate \(p=1{0}^{-3}\,(1{0}^{-5})\). Alternatively, 57 million (3.8 million) physical qubits can run the algorithm in a day (half a day) as a result of space-time trade-off. Essentially, a large number of high purity \(T\)-gates for modular exponentiation in Shor’s algorithm require a formidable amount of physical qubits. On the other hand, a lower bound is set at 12.2 million (1.46 million) physical qubits for 6000 logical qubits (2028 (243) physical qubits per a logical qubit) for distance-25 (8) surface code.

Fig. 6
figure 6

Number of physical qubits vs. time for factoring 2000-bit numbers with Shor’s algorithm with a physical error rate \(p=1{0}^{-3}\) for a and \(p=1{0}^{-5}\) for b. Darker lines assume a bond trial time t0 = 1 µs, and lighter lines denote \({t}_{0}=100\) ns. See Supplementary Discussion for detailed information and calculation. Results marked with red lines use (double defect) braiding qubits52 with two-step 15-to-1 distillation for high purity \(T\)-gate creation. Green lines show results with the lattice surgery qubits62 and catalyzed-\(2T\) factories.53 Windowed arithmetic and autoCCZ factories dramatically reduce the resource overhead (blue).54 The last result incorporates space-time layout64 implying that the improvement is even larger. Results are terminated on the left-hand side by either measurement time (red and green) or surface code cycle (blue) and on the right hand side by logical error rates.

Discussion

We assumed here that both bit and phase-flip probabilities are the same. However, our results may be significantly improved using recent work on tailoring the surface code65 if the noise is biased. We consider other important properties of NVs in the Supplementary Notes. For example, NVs can be ionized under strong optical pulses, but their charge state can be recovered by optical repumping. We designed the microwave and optical pulse sequences so that the nuclear spin is not disturbed by failed trials and subsequent spin-state initialization. However, precise control of the microwave and optical transitions presents technical challenges and mark an area of active research.18,66 Future work should also refine the error model to explicitly include contributions from the higher order internal dynamics of NVs.

Though there have been huge progress in reducing the resource overhead in fault-tolerant quantum computing, it still requires a large number of physical qubits and long computation time owing to slow clock cycles. In this aspect, noisy intermediate scale quantum (NISQ) technology67 can be a more near to medium term path for our proposal. Two observations are favorable to the NISQ direction; the proposed system well behaves with only 25 qubits as shown in Fig. 5a; the LCC size scales linearly with the total number of qubits in the supercritical regime. Thus, the quantum simulation on this architecture quickly reaches thermodynamic behavior.68

We emphasize that the architecture is compatible with any quantum memories with enough coherence time. Especially, silicon vacancy centers and other group IV emitters such as GeV, SnV, or PbV showed promising emission properties.69,70 For example, silicon vacancy centers have higher Debye–Waller factor of \(\sim 0.7\), with narrow inhomogeneous distribution of \(\sim 51\) GHz.71 Emission wavelength of emitters can be matched by strain tuning with tuning range upto \(\sim 440\) GHz.72 The spin coherence time of \(10\) ms has been shown at dilution fridge temperature.73