Experimental system
The sample used in this work is a type-IIa electronic-grade synthetic diamond (Element Six) with a natural abundance of 13C impurities. The NV centre is at the focus of a solid immersion lens encircled by an antenna for microwave (mw) frequency control. All experiments are performed at room temperature in ambient conditions. A permanent magnet was aligned to the NV symmetry axis using pulsed electron-spin resonance experiments and positioned to create a magnetic-field strength of 338 G. The magnetic-field strength was chosen to minimize nuclear-qubit gate durations and angular errors. Further details of the field alignment and simulations to determine the field strength are provided in Supplementary Section VII.
Green (532-nm) laser pulses of 2 μs were used to (re)initialize the electron spin and charge state through optical pumping, and shorter 300-ns pulses were used to measure the spin-state photoluminescence contrast. The synchronization of the optical and mw signals was achieved using two different configurations. The first used two arbitrary waveform generators, one (Tektronix AWG520) dedicated to optical control and the other (Tektronix AWG7102), for mw control. The second configuration used a Swabian Instruments PulseStreamer 8/2 for both optical and mw control. Additional details are provided in Supplementary Section I. Electron gate errors were quantified using bootstrap tomography of pulses50 (Supplementary Section II).
DD
The Hamiltonian governing the central spin electron interacting with L nuclear qubits is given by
$$H={\mathbb{1}}\otimes \frac{{\omega }_{{\rm{Lar}}}}{2}\mathop{\sum }\limits_{\ell =1}^{L}{\sigma }_{z}^{(\ell )}+\frac{{Z}_{e}}{2}\otimes \mathop{\sum }\limits_{\ell =1}^{L}({A}_{| | }^{(\ell )}{\sigma }_{z}^{(\ell )}+{A}_{\perp }^{(\ell )}{\sigma }_{x}^{(\ell )}),$$
(3)
where ωLar is the nuclear Larmor frequency; \({Z}_{e}={s}_{0}\left\vert 0\right\rangle \left\langle 0\right\vert +{s}_{1}\left\vert 1\right\rangle \left\langle 1\right\vert\) is the electron-spin operator, where sj are the two electron-spin projections chosen as the computational basis (s0 = 0 and s1 = −1 for this work); and \({A}_{\!\parallel,\!\perp }^{(\ell )}\) are the parallel and perpendicular hyperfine couplings between the electron and the ℓth nuclear qubit. This can be rewritten as40
$$H=\sum _{j\in \{0,1\}}\left\vert j\right\rangle {\left\langle j\right\vert }_{e}\otimes \mathop{\sum }\limits_{\ell }^{L}{H}_{j}^{(\ell )},$$
(4)
where each \({H}_{j}^{(\ell )}\) is given by
$${H}_{j}^{(\ell )}=\frac{{\omega }_{L}+{s}_{\!j}{A}_{| | }^{(\ell )}}{2}{\sigma }_{z}^{(\ell )}+\frac{{s}_{\!j}{A}_{\perp }^{(\ell )}}{2}{\sigma }_{x}^{(\ell )}.$$
(5)
The notation \({\sigma }_{i}^{(\ell )}\) in equation (5) means the ith Pauli matrix on the ℓth component of the L-nuclear-qubit Hilbert space and the identity on all other components. This form of the Hamiltonian highlights how the electron-state conditions are different with unique dynamics for each nuclear qubit. This is further made apparent by the free evolution operator Uf(t) for the system:
$${U}_{f}(t)=\sum _{j\in \{0,1\}}\left\vert j\right\rangle {\left\langle j\right\vert }_{e}\mathop{\bigotimes }\limits_{\ell }^{L}\exp \left(-{\rm{i}}t{H}_{j}^{(\ell )}\right),$$
(6)
from which each \(\exp (-{\rm{i}}t{H}_{j}^{(\ell )})\) term can be viewed as a rotation operator acting on the ℓth nuclear qubit. Note a subtle shift in notation from equation (5) to (6), where each index ℓ no longer implies identity operators on the other qubits, and each two-dimensional \({H}_{j}^{(\ell )}\) can be viewed as acting on a distinct subspace. Additional details and derivations are provided in Supplementary Section III.
The free evolution periods of DD sequences leverage equation (6) to control the rotational effects of each nuclear qubit, as well as extend the electron coherence time. The net unitary operator UDD from performing a time-symmetric DD sequence of unit-pulse time t and with N repeats is given by
$${U}_{{\rm{DD}}}=\sum _{j\in \{0,1\}}\left\vert j\right\rangle {\langle j\vert }_{e}\mathop{\bigotimes }\limits_{\ell }^{L}{R}_{{\hat{\mathbf{n}}}_{j}^{(\ell )}(t)}(N{\phi }^{(\ell )}(t)),$$
(7)
where R is a spin-1/2 rotation operator about the axis \({\hat{\mathbf{n}}}_{j}^{(\ell )}\) and by an angle of Nϕ(ℓ) for the ℓth nuclear qubit. Supplementary Section V provides details in calculating each rotation operator based on the hyperfine couplings of the register. This formulation highlights the conditional nature of each nuclear qubit’s rotation depending on the electron state \({\left\vert j\right\rangle }_{e}\).
Resonant X-axis control of a target nuclear qubit is achieved with the proper choice of unit-pulse time tm that creates \({({\hat{\mathbf{n}}}_{0}\cdot {\hat{\mathbf{n}}}_{1})}^{(\ell )}({t}_{m})=\pm 1\). Such a choice of tm occurs periodically and is given by
$${t}_{m}^{(\ell )}=\frac{4\uppi m}{{\omega }_{0}^{(\ell )}+{\omega }_{1}^{(\ell )}},$$
(8)
for \(m\in {{\mathbb{Z}}}^{+}\) and \({\omega }_{j}^{(\ell )}=\sqrt{{({s}_{\!j}{A}_{\perp }^{(\ell )})}^{2}+{({\omega }_{L}+{s}_{\!j}{A}_{\parallel }^{(\ell )})}^{2}}\) (ref. 40). For odd m = 2k + 1, \({({\hat{\mathbf{n}}}_{0}\cdot {\hat{\mathbf{n}}}_{1})}^{(\ell )}=-1\) and for even m = 2k, \({({\hat{\mathbf{n}}}_{0}\cdot {\hat{\mathbf{n}}}_{1})}^{(\ell )}=+1\). Here the integer k specifies the DD order, as discussed earlier. When the electron-state-dependent nuclear rotation axes are maximally anti-aligned, or \({({\hat{\mathbf{n}}}_{0}\cdot {\hat{\mathbf{n}}}_{1})}^{(\ell )}=-1\), the nuclear rotations are maximally dependent on the state of the electron. With N set to create the correct rotation angle, the net gate is \({C}_{e}{X}_{\ell }(\pm \pi /2)=\)\(|0\rangle {\langle 0|}_{e}\otimes {X}_{\ell }(\pi /2)+\)\(|1\rangle {\langle 1|}_{e}\otimes {X}_{\ell }(-\pi /2)\) between the electron and target nuclear qubit qℓ. For all other spins, the choice of t is off-resonance, and the resulting rotation is unconditional and about the Z axis. Similarly, when the unit-pulse time is on resonance and \({({\hat{\mathbf{n}}}_{0}\cdot {\hat{\mathbf{n}}}_{1})}^{(\ell )}=+1\), the resulting ℓth nuclear qubit’s rotation is an unconditional X-axis rotation, with all other nuclear rotations being off-resonance and about the Z axis.
For example, when attempting to rotate the first nuclear qubit unconditionally about the X axis by π/2, the net unitary acting on the register would take the form \(U={I}_{e}\otimes {X}_{\uppi /2}\otimes {Z}_{{\theta }_{(2)}}\ldots \otimes {Z}_{{\theta }_{(L)}}\), where each crosstalk rotation angle θ(ℓ) depends on the choice of t and N that were used to achieve the desired Xπ/2 rotation of q1 and the specific hyperfine couplings of the ℓth nuclear qubit. The goal of the parallelized entangling gate is to leverage this crosstalk in such a way that each nuclear qubit can be maximally entangled for a single choice of t and N. Further information on the t and N parameter choices for each nuclear qubit’s gate, together with their experimental verification, is provided in Supplementary Section VIII.
Entanglement metrics
To quantify the bipartite entangling ability of a DD sequence with a particular nuclear qubit qℓ, one can calculate the first Makhlin invariant, which takes the form
$${G}_{1}^{(\ell )}={\left({\cos }^{2}\frac{N{\phi }^{(\ell )}}{2}+{({\hat{\mathbf{n}}}_{0}\cdot {\hat{\mathbf{n}}}_{1})}^{(\ell )}{\sin }^{2}\frac{N{\phi }^{(\ell )}}{2}\right)}^{2},$$
(9)
for time-symmetric DD sequences such as XY8 (ref. 40). This entanglement metric (bounded from 0 to 1) is minimal when bipartite entanglement is maximal. Using this form of \({G}_{1}^{(\ell )}\), it was shown that with the proper choice of N, \({G}_{1}^{(\ell )}=0\) if \({({\hat{\mathbf{n}}}_{0}\cdot {\hat{\mathbf{n}}}_{1})}^{(\ell )} < 0\). Finding the unit-pulse times that satisfies this condition for each target nuclear qubit is the first step in calibrating the parallel entangling gate. Furthermore, to quantify the multipartite entangling ability of a DD sequence with L target nuclear qubits, one can use the M-qubit entangling power:
$${\varepsilon }_{{\rm{p}},M}({U}_{\rm{DD}})={\left(\frac{d}{d+1}\right)}^{M}\mathop{\prod }\limits_{\ell }^{L}(1-{G}_{1}^{(\ell )}),$$
(10)
where M = L + 1 is the total number of qubits targeted for entangling, including the electron, and d = 2 is the dimension of the qubit subspace41. Often, as shown in Fig. 1, the normalized version of this metric is the most useful, without the constant coefficient in front of the product. The normalized metric ranges from 0 (the DD sequence creates no entanglement) to 1 (the DD sequence is a maximal multipartite entangler). Owing to the central spin nature of solid-state defect systems, εp,M(UDD) depends only on each bipartite entanglement invariant \({G}_{1}^{(\ell )}\). Calculating εp,M(UDD) with each of the target nuclear qubits in the range of unit-pulse times that satisfy \({({\hat{\mathbf{n}}}_{0}\cdot {\hat{\mathbf{n}}}_{1})}^{(\ell )} < 0\) reveals the optimum (t, N) combination to generate maximal multipartite entanglement.
We utilize the non-unitary entangling power to account for the impact of residual entanglement generated with non-targeted nuclear spins41. This entanglement metric is derived using the partial trace quantum channel \({\mathcal{E}}\) over non-targeted nuclear qubits. The set of all nuclear qubits is partitioned into a subset that is targeted (size L) and the rest that are not targeted (size, Ltotal − L). A simple approximate form for this entanglement metric is given by
$${\varepsilon }_{{\rm{p}},M}({\mathcal{E}})=\frac{{\varepsilon }_{{\rm{p}},M}({U}_{\rm{DD}})}{2}\left(1+\mathop{\prod }\limits_{{{\ell \in \,{\text{not}}}\atop {\text{targeted}}}}^{{L}_{{\rm{total}}}-L}{G}_{1}^{(\ell )}\right).$$
(11)
The non-unitary entangling power is bounded above by the unitary entangling power; \({\varepsilon }_{{\rm{p}},M}({\mathcal{E}})\le {\varepsilon }_{{\rm{p}},M}({U}_{\rm{DD}})\), with equality holding when no residual entanglement is generated (\({G}_{1}^{(\ell )}=1\) for all non-targeted nuclear qubits)41. Therefore, any residual entanglement generated leads to a Makhlin invariant less than 1, lowering this entanglement metric. This metric is essential when designing parallel entangling gates with subsets of known nuclear qubits, for example, in the case of L = 2 parallel entangling gates in this work. Additional details regarding how these metrics were used to calibrate each parallel entangling gate are provided in Supplementary Section V.
MQCs
In the original MQC circuit proposed in ref. 43 (Fig. 2a), the M-qubit register is first initialized to \({\left\vert 0\right\rangle }^{\otimes M}\). The control (top) qubit qc is then placed into an equal superposition state so that the following CNOT gates create a GHZ state. Once entangled, each qubit’s relative phase is shifted by an equal amount ϕ, yielding \(\left\vert {\,\text{GHZ}\,}_{\phi }^{M}\right\rangle =\frac{1}{\sqrt{2}}({\left\vert 0\right\rangle }^{\otimes M}+{e}^{-iM\phi }{\left\vert 1\right\rangle }^{\otimes M})\). The system is then disentangled back to the original state by reversing the first half of the circuit. The result (before the last Hadamard that projects the control qubit phase onto the measurement axis) is that the control qubit’s phase is amplified based on how many qubits it was entangled with:
$$\left\vert {\psi }_{f}\right\rangle =\frac{1}{\sqrt{2}}(\left\vert 0\right\rangle +{e}^{-{\rm{i}}M\phi }\left\vert 1\right\rangle )\otimes {\left\vert 0\right\rangle }^{\otimes M-1}.$$
(12)
Thus, the final probability of the entire system returning to the initial state is given by
$$P\left({\left\vert 0\right\rangle }^{\otimes M}\right)=\frac{1}{2}(1+\cos (M\phi )),$$
(13)
which crucially carries a frequency equal to the number of qubits in the entangled state.
Experimental MQC considerations
Nuclear phase gates Zϕ were implemented using off-resonant DD sequences. Since such gates are realized for any off-resonant t, the optimum parameters can be chosen strategically. The off-resonance region before the first-order resonances not only offers fast pulse times (t < 1 μs) but also can parallelize the phase gate, turning the sequence into an unconditional M-qubit gate. Experimentally, the finite pulse duration of electron gates sets a lower limit on t; this restriction, in turn, sets a lower bound on the angular resolution Δϕ = ϕN=1 of the phase gate. With Δϕ specified, t was optimized to minimize the angular error for each nuclear qubit in the register. Then, to increase the phase, the unit pulse was repeated N times, leading to ϕ = NΔϕ. The simulated four-qubit process fidelities for a parallelized gate of the form \({I}_{e}\otimes {Z}_{\phi }^{\otimes 3}\) were ~99% for ϕ = π/2. Further details and a table of pulse parameters are provided in Supplementary Section V.
Entangling-gate fidelities
M-qubit state fidelities are calculated according to the trace overlap of the quantum state ρ with the target state ρtarget; FM = tr(ρ × ρtarget). On the basis of the form of the bipartite and sequential entangling gates, when NE is a multiple of four, the target state is \({\left\vert 0\right\rangle }^{\otimes M}\)—same as the initial state. The repeats of the parallel gate were chosen to maximize the overlap with \({\left\vert 0\right\rangle }^{\otimes M}\) as a target state. This simplifies the state fidelities to be given by only a single component of ρ; \({F}_{M}={\left\langle 0\right\vert }^{\otimes M}\rho {\left\vert 0\right\rangle }^{\otimes M}\). The fidelity of this separable target state can be further approximated by independent Z-axis measurements of each qubit:
$${F}_{M}\approx \frac{1}{{2}^{M}}\mathop{\prod }\limits_{\ell =1}^{M}(1+\langle {Z}_{\ell }\rangle ).$$
(14)
This approximation ignores correlations between qubits, which is reasonable since the initial and final states are separable with vanishing pairwise covariances and cumulants. Experimentally, the electron Z-axis projection is measured directly using spin-dependent fluorescence, whereas nuclear qubits are measured using Z-axis tomography (Fig. 3b). Supplementary Section X provides a derivation of equation (14) and additional details regarding the fidelity measurements.
Generality and extensions
Random NV-nuclear-qubit registers were generated by uniformly sampling nuclear spin positions within a spherical volume surrounding an NV centre. The radius of this sphere was set to 2.3 nm, which encapsulates approximately 100 nuclear spins at natural abundance (1.1%). Nuclei at the surface of this sphere contribute to the spin bath. From the location of each nuclear qubit, the hyperfine matrix was calculated using the dipole–dipole interaction. The boundary between the strongly and weakly coupled qubits was set by the inhomogeneous linewidth of the electron-spin transitions, \(\sqrt{2}/\pi {T}_{2}^{* }\approx 200\) kHz in this work. If a register contained any strongly coupled nuclear qubits, the case was not considered further (Fig. 5a, red region). Approximately 64% of the randomly generated registers contained at least one nuclear qubit with a hyperfine component larger than this cut-off. The remaining 36% of registers (Fig. 5a) were evaluated for parallel entanglement. We further applied lower bounds to separate addressable nuclei from the spin bath: ∣A∥∣ > 15 kHz (based on the location of the spin-bath resonance) and A⊥ > 10 kHz (so that a sufficiently small N can address the qubit). With these cut-offs, the average number of weakly coupled, addressable nuclear qubits per register is 5.7, with a standard deviation of 2.5 at natural 13C concentration.
For each of these registers, we searched for parallel entangling gates following the algorithm in ref. 41. Additional details are provided in Supplementary Sections V and VI. Because a large number of registers were generated to improve statistical significance, a conservative parallel entangling gate search was used. Specifically, a maximum of N ≤ 50 and a minimum non-unitary entangling power of \({\varepsilon }_{{\rm{p}},M}({\mathcal{E}})\ge 0.8\) were imposed. Hence, these results represent a lower bound on the available gates. Supplementary Section VI provides additional simulations of gate durations, comparisons with k = 2 and k = 3 sequential two-qubit gates, and the infidelity arising from residual entanglement.