Optimisation of diamond quantum processors

Diamond is a promising architecture for quantum information processing [1–8] and quantum sensing/microscopy [9–12] at both cryogenic and room temperatures. Optimised diamond quantum processors are crucial building blocks for large-scale diamond quantum computers and the next generation of quantum sensors and microscopes that are enhanced by embedded quantum memories and signal processing [13–15]. To date, diamond quantum processors have been used to implement quantum error correction codes [3, 4, 16], quantum algorithms [17, 18], detection of metallo-protein molecules [19] as well as quantum simulation of the helium hydride cation [5] and topological phase transition of a quantum wire [6]. For technology applications, diamond quantum processors are distinguished from other quantum architectures due to their ability to operate in ambient conditions and with relatively simple microwave, radio-frequency and off-resonant optical control systems [20]. The resulting improvements in complexity, robustness and cost make diamond one of the most flexible and widely applicable quantum technology platforms.

Diamond quantum processors consist of a nitrogen-vacancy (NV) centre with a local cluster of hyperfine-coupled nuclear spins. These coupled nuclear spins include the intrinsic N nuclear spin of the NV centre and isotopic ¹³C lattice impurities. Quantum computations are realised by using the electron spin of the NV centre as a quantum bus that initialises, mediates interactions between, and reads-out the coupled nuclear spins, which act as the physical qubits. Scaling of the diamond quantum architecture requires a mechanism to couple multiple NV centres and their qubit clusters. Coupling of separate NV centres has been demonstrated at cryogenic temperatures using photons [21–23], and at room temperature using magnetic dipole coupling between proximate NV centres [2, 24]. Spin chains [54] and coherent spin transport [25] have also been proposed as coupling mechanisms.

Previous work in optimising the quantum control of either the NV centre's electron spin or coupled nuclear spin qubits has already demonstrated excellent control fidelities, using techniques including dynamical decoupling [7, 8, 26–37] as well as numerical gate pulse-shaping techniques, such as the chopped random basis (CRAB) quantum optimisation algorithm [38, 39] and the gradient ascent pulse engineering (GRAPE) algorithm [40]. Experimental application of the CRAB algorithm to the NV electron spin has yielded ultra-fast single-qubit gates (i.e. beyond the rotating wave approximation) with fidelities of 0.95 ± 0.01 and 0.99 ± 0.016 for π/2 and π pulses, respectively [41]. The GRAPE algorithm has been used to demonstrate single electron spin operations with fidelity F ≈ 0.99 [2] and generate entangled states of three nuclear spins with fidelities exceeding 85% [4]. Moreover, an average single-qubit gate fidelity of 0.999 952 and two-qubit gate fidelity of 0.992 has been reported using composite pulses and a modified GRAPE algorithm, respectively [42]. While some of the demonstrated gate fidelities are impressive, even with a gate fidelity of 0.999, for a fault-tolerant logical qubit to achieve logical error rates comparable to classical computers, effective surface code error correction is anticipated to require up to 10⁴ physical qubits [43, 44]. Therefore, there is a strong motivation to push for further reduction in gate errors.

There are two different optimisation problems to address to improve the performance of diamond quantum processors: improvement of initialisation/readout fidelities and improvement of gate fidelities and speeds. The initialisation/readout fidelities are optimised by selecting the ¹³C nuclear spin lattice sites that have the longest nuclear spin relaxation time during the projective single-shot optical readout process employed in diamond quantum processors [4]. Broadly speaking, the best lattice sites are those whose hyperfine field is well aligned to NV centre's axis. As will be discussed later, overcoming the effects of hyperfine field misalignment and decoherence are crucial to achieve high gate fidelities. As such, there is some correlation between the choices made to optimise initialisation/readout fidelities and gate fidelities. Consequently, for simplicity, in the following we will not discuss the optimisation of initialisation/readout fidelities. Instead, we assume a particular selection of nuclear spin qubits with well-aligned hyperfine fields, and focus on the problem of optimising the gate fidelities due to control errors and speeds.

Optimisation of the processor gate operations requires simultaneous maximisation of gate speed and minimisation of: (1) spurious effects of control fields on qubits other than the target qubit(s) (i.e. cross-talk), (2) the effects of random control field errors (i.e. those caused by fluctuations in amplitude, frequency and phase), (3) the effects of decoherence, and (4) the effects of hyperfine field misalignment. Additional practical requirements are that the gate design: (A) complies with the physical constraints of the control systems, (B) is readily incorporated into a feedback-based optimisation routine that uses measurements to optimise the actual physical processor (and not just models of the processor) and supports updating of the optimisation during operation to adjust for system drifts, and (C) is parameterised so that the degree of convergence to optimisation limits can be deterministically and systematically assessed (in order to support design decisions concerning the costs and benefits of further optimisation) and the dominant error modes can be diagnosed (to improve processor design). As GRAPE involves direct numerical solving of pulses using model systems, it does not readily achieve the practical requirements (B) and (C) and relies heavily on the accuracy of its models [45]. On the other hand, the nature of CRAB allows the integration of (B) but not (C) directly due to its inherent reliance on random numbers, and the number of free parameters required to optimise the control field [39].

We propose a different approach to this optimisation problem. Our approach has three steps. The first step is to generate a complete semi-analytical basis of pulses that comply with (A) and minimises (1) in the approximation where time-ordering is neglected in the simulated quantum evolution. The second step is to find the linear combination of these basis functions that minimise (2). The final step is to further optimise the pulses to account for the effects of time-ordering in the quantum evolution via closed-loop optimisations. We have neglected time-ordering in the quantum evolution to efficiently estimate the initial coefficients of the linear combinations of basis functions required for closed-loop optimisations. Without an initial estimate, optimising the linear combinations of basis functions via closed-loop optimisations is computationally time consuming. This linear minimisation is fast and can include measurements of the processor and therefore readily complies with (B). Furthermore, since the basis is complete, the dimension of the non-trivial basis functions provides a clear parameter to analyse convergence to optimisation limits and interpretation of different error modes, and thus also complies with (C). The principal strategy for both minimising the effects of decoherence (3) and the effects of hyperfine field misalignment (4) is to maximise the gate speed. This minimises the time over which decoherence accumulates and the time under which the nuclear spin experience the hyperfine misalignment during two-qubit gates, while ensuring that the infidelity introduced by control errors is less than that introduced by decoherence and hyperfine field misalignment. We have adopted this strategy to minimise the errors due to decoherence rather than directly targeting the primary source of qubit decoherence (e.g. via dynamical decoupling) because the primary source is the NV centre's electron spin, whose strong interaction with the nuclear spins is required for the implementation of both one- and two-qubit gates. Indeed, owing to the strong interaction, the qubit coherence time 1/T_2,n is limited by the electron spin relaxation time 1/T_1,e [46, 47], which is approximately 2.4 ms at room temperature [49]. Note that our discussion here, and the results of this paper, are in the context of room-temperature operation of diamond quantum processors. At cryogenic temperatures, weakly coupled nuclear spins can instead be employed as qubits and as a result, their decoherence is influenced by other mechanisms [50], at the cost of slower gate speeds.

In this paper, we report (i) a demonstration and analysis of our optimal design approach and (ii) simulation of an optimised diamond quantum processor. In section 2, we first discuss the operating principles of a diamond quantum processor before presenting a complete model of the processor, its gate operations and control system errors. Section 3 demonstrates the generation of gate basis functions, while section 4 demonstrates the optimisation of gate fidelities in the presence of control errors. In section 5, the effects of decoherence on a diamond quantum computer are investigated via master equation simulations and in section 6, we simulated the performance of 3 and 5 qubit quantum Fourier transforms (QFTs) on a diamond quantum processor. Fidelities of QFTs were chosen as a simple performance metric because QFTs are the foundation of many quantum algorithms. Thus, the fidelities of QFTs are basic indicators for the system's performance with more sophisticated algorithms.

2.1. Operating principles of diamond quantum processors

The NV centre is a point defect in diamond consisting of a substitutional nitrogen and an adjacent carbon vacancy [51]. Its electronic structure consists of a ground state spin triplet (³A₂) and an excited state spin triplet (³E) with two intermediate singlet levels (¹E and ¹A₁). There exist spin-selective non-radiative intersystem crossings between the triplet and singlet levels, which lead to initialisation of the electronic spin state upon optical excitation of the centre's ³A₂ → ³E transition as well as read out via the differing fluorescence intensities of the spin states (see reference [51] for further details). In addition to high fidelity optical spin initialisation and readout, the NV centre also has the longest electron spin coherence time of any solid state spin at room temperature (T₂ ≈ 2.4 ms) [49].

Each NV centre is coupled to a register of one or more nuclear spins, which we use as qubits [3, 4]. The quantum register consists of the NV centre's intrinsic nitrogen nuclear spin and nearby ¹³C nuclear spins. Hyperfine coupling between the NV electron spin and the nuclear spins results in a splitting of the electronic and nuclear energy levels. This splitting depends on the particular hyperfine coupling strength between each nucleus and the NV electron spin, and also on the respective electron and nuclear spin states [52, 53] (see figure 1(b)). We choose a register with non-overlapping hyperfine couplings, allowing the use of frequency selectivity to individually address each nuclear spin qubit in the register.

Zoom In Zoom Out Reset image size

Key requirements for universal quantum computation are the initialisation and readout of the qubits, as well as the ability to apply single and two-qubit gate operations. In diamond quantum computing, each of these processes relies on high-fidelity quantum gates on the electron and nuclear spins. Initialisation and readout of a diamond quantum register is performed via a projective, single-shot readout of the nuclear spin qubits. This measurement scheme involves initialising the electron spin, entangling the nuclear spin qubits with the electron spin using a C_nNOT_e gate and then readout of the electron spin [58]. Single-qubit gate operations are realised using radiofrequency (RF) pulses. These pulses correspond to the R_x and R_y gates where they are the rotations about the x and y axes respectively. Other single-qubit gates can be constructed from combinations of these rotations. The intrinsic properties of the NV-nuclear spins system allows direct application of a CZ gate via microwave (MW) pulses. This CZ gate is achieved by performing a selective 2π pulse conditional on the nuclear spin register being in a particular state [4, 56, 57]. The CZ gate can be combined with single-qubit gates to realise any other two-qubit gate.

The splitting in the ³A₂ triplet ground state results in two types of subspaces which we identify as computational subspace and auxiliary subspace (figure 1(b)). The natural computational subspace is either of the m_s = ±1 states as they have non-zero hyperfine interactions, thus allowing the nuclear spin qubits to be individually addressed through frequency selectivity [59]. While the choice of either the m_s = ±1 state as the computational subspace is arbitrary, the m_s = −1 state is more often selected as it requires lower microwave frequencies for qubit gate operations. Single-qubit gates are realised in the computational subspace while a two-qubit CZ gate utilises, but does not occupy, the auxiliary subspace [4, 56, 57].

2.2. Model Hamiltonian

The Hamiltonian H_I of the nuclear spins coupled to the NV centre is

$$\begin{equation}{H}_{I}=\sum _{i} \overrightarrow {S}\cdot {\boldsymbol{A}}_{i}\cdot { \overrightarrow {I}}_{i}-\sum _{i}{\gamma }_{i}\left[{I}_{i,z}{B}_{0}+{I}_{i,x}{B}_{1}\left(t\right)\right]\end{equation} \tag{ 2.1 }$$

where γ_i is the gyromagnetic ratio of the ith nucleus, B₀ is the background static magnetic field aligned with the NV axis, B₁(t) is the applied radio frequency field, A_i is the hyperfine tensor of the ith nucleus with $\overrightarrow {S}$ being the dimensionless electron spin operator and ${ \overrightarrow {I}}_{i}$ is the dimensionless nuclear spin operator of the ith nucleus. For this model, we apply the secular approximation as a very strong magnetic field is applied along the z-axis during the operation of this quantum computer. Therefore, the nuclear spin Hamiltonian in the computational subspace simplifies to

$$\begin{equation}{H}_{I,{m}_{\mathrm{s}}=-1}=-\sum _{i}{ \overrightarrow {A}}_{i,z}\cdot { \overrightarrow {I}}_{i}-\sum _{i}{\gamma }_{i}\left[{I}_{i,z}{B}_{0}+{I}_{i,x}{B}_{1}\left(t\right)\right]\end{equation} \tag{ 2.2 }$$

${ \overrightarrow {A}}_{i,z}$ is now a vector instead of a tensor where

$$\begin{equation}{ \overrightarrow {A}}_{i,z}={A}_{i,xz}\hat{x}+{A}_{i,yz}\hat{y}+{A}_{i,zz}\hat{z}\end{equation} \tag{ 2.3 }$$

We simplify the expression by diagonalising the nuclear spin Hamiltonian in the computational subspace via a rotation of the spin operators about the angles defined by their hyperfine interactions. Assuming that only nuclei with hyperfine fields nearly aligned with the NV axis are chosen, we perform small angle approximations and by undoing a rotation about z-axis for further simplifications, this yields

$$\begin{equation}{H}_{I,{m}_{\mathrm{s}}=-1}=-\sum _{i}{\omega }_{i}{I}_{i,z}-\sum _{i}{I}_{i,x}{\gamma }_{i}{B}_{1}\left(t\right)\end{equation} \tag{ 2.4 }$$

where ω_i is the transition frequency of the ith nucleus. We also transform the Hamiltonian into the interaction picture using the following transformation operator

$$\begin{equation}{H}_{\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{g}\mathrm{l}\mathrm{e}}=\mathcal{T}{H}_{I,{m}_{\mathrm{s}}=-1}{\mathcal{T}}^{-1}-\mathbb{i}\mathcal{T}\enspace \frac{\mathrm{d}}{\mathrm{d}t}\enspace {\mathcal{T}}^{-1}\end{equation} \tag{ 2.5 }$$

where

$$\begin{equation}\mathcal{T}=\mathrm{exp}\left[-\mathbb{i}t\sum _{i}{\omega }_{i}{I}_{i,z}\right]\end{equation} \tag{ 2.6 }$$

Thus, for single-qubit gate operations, the Hamiltonian for the computational subspace in the interaction picture is given by

$$\begin{equation}{H}_{\text{single}}=-\sum _{i}\left[{I}_{i,x}\enspace \mathrm{cos}\enspace {\omega }_{i}t+{I}_{i,y}\enspace \mathrm{sin}\enspace {\omega }_{i}t\right]{\gamma }_{i}{B}_{1}\left(t\right)\end{equation} \tag{ 2.7 }$$

The auxiliary subspace is involved to perform a CZ gate and this is enabled via the NV electron spin. The effective Hamiltonian for two-qubit gate operations is given by

$$\begin{equation}{H}^{\prime }=\frac{{\Delta}}{2}{\sigma }_{z}+\frac{{\Omega}}{2}{\sigma }_{x}{B}_{1}\left(t\right)+\sum _{i}\left({\alpha }_{i}+{\beta }_{i}{\sigma }_{z}\right){I}_{i,z}\end{equation} \tag{ 2.8 }$$

with Δ = D − γ_eB₀ where D denotes the zero-field splitting, ${\Omega}=\sqrt{2}{\gamma }_{\mathrm{e}}$, α_i = −γ_iB₀/2 − ω_i/2 and β_i = −γ_iB₀/2 + ω_i/2. When defining the Hamiltonian above, we ignore the interactions with the m_s = +1 state of the electron spin and the direct interaction between the nuclear spins and the microwave field B₁(t). It is possible to do this because the microwaves are far detuned from these transitions. We have also ignored the interactions of the static magnetic field acting on the nuclear spin when the electron spin is in m_s = 0 state. This interaction term arises from the hyperfine field misalignment and is discussed further in appendix B.2. Likewise to the Hamiltonian for single-qubit gate operations, we transform H' into the interaction picture with the transformation operator ${\mathcal{T}}^{\prime }$ is given by

$$\begin{equation}{\mathcal{T}}^{\prime }=\mathrm{exp}\left[\mathbb{i}t\left(\frac{{\Delta}}{2}{\sigma }_{z}+\sum _{i}\left({\alpha }_{i}+{\beta }_{i}{\sigma }_{z}\right){I}_{i,z}\right)\right]\end{equation} \tag{ 2.9 }$$

The transformed Hamiltonian for two-qubit gate operations is then given by

$$\begin{align}\hfill {H}_{\text{multi}}& =\frac{{\Omega}}{2}{B}_{1}\left(t\right)\left[{\sigma }_{x}\otimes \left(\left\vert 11\right\rangle \left\langle 11\right\vert \enspace \mathrm{cos}\left(\left[{\Delta}+{\beta }_{1}+{\beta }_{2}\right]t\right)+\left\vert 10\right\rangle \left\langle 10\right\vert \enspace \mathrm{cos}\left(\left[{\Delta}+{\beta }_{1}-{\beta }_{2}\right]t\right)\right.\right.\hfill \\ \hfill & \left.\quad +\left\vert 01\right\rangle \left\langle 01\right\vert \enspace \mathrm{cos}\left(\left[{\Delta}-{\beta }_{1}+{\beta }_{2}\right]t\right)\right)\left.+\left\vert 00\right\rangle \left\langle 00\right\vert \enspace \mathrm{cos}\left(\left[{\Delta}-{\beta }_{1}-{\beta }_{2}\right]t\right)\right)\hfill \\ \hfill & \quad -{\sigma }_{y}\otimes \left(\left\vert 11\right\rangle \left\langle 11\right\vert \enspace \mathrm{sin}\left(\left[{\Delta}+{\beta }_{1}+{\beta }_{2}\right]t\right)+\left\vert 10\right\rangle \left\langle 10\right\vert \enspace \mathrm{sin}\left(\left[{\Delta}+{\beta }_{1}-{\beta }_{2}\right]t\right)\right.\hfill \\ \hfill & \quad \left.\left.+\left\vert 01\right\rangle \left\langle 01\right\vert \enspace \mathrm{sin}\left(\left[{\Delta}-{\beta }_{1}+{\beta }_{2}\right]t\right)+\left\vert 00\right\rangle \left\langle 00\right\vert \enspace \mathrm{sin}\left(\left[{\Delta}-{\beta }_{1}-{\beta }_{2}\right]t\right)\right)\right]\hfill \end{align} \tag{ 2.10 }$$

where 0 and 1 are the m_I = −1/2 and +1/2 nuclear spin projections, respectively. We use the notation where the most left entry of a tensor product corresponds to the first qubit, i.e. $\left\vert {q}_{1},{q}_{2},\dots \enspace \right\rangle$

2.3. Control pulses and gate operations

Focusing on single-qubit gate operations within the computational subspace, the applied radio frequency field B₁(t) can be parametrised as a linear combination of oscillating components where

$$\begin{equation}{\gamma }_{i}{B}_{1}\left(t\right)=a\left(t\right)\mathrm{cos}\enspace \omega t+b\left(t\right)\mathrm{sin}\enspace \omega t\end{equation} \tag{ 2.11 }$$

By neglecting the time-ordering in the quantum evolution, we can write the evolution operator for the ith nucleus as

$$\begin{equation}{U}_{i}={\mathrm{e}}^{-\mathbb{i}{I}_{i,x}{X}_{i}}\quad \mathrm{o}\mathrm{r}\quad {U}_{i}={\mathrm{e}}^{-\mathbb{i}{I}_{i,y}{Y}_{i}}\end{equation} \tag{ 2.12 }$$

with

$$\begin{equation}{X}_{i}=-\underset{-\tau /2}{\overset{\tau /2}{\int }}a\left(t\right)\mathrm{cos}\enspace {\omega }_{i}t\enspace \mathrm{cos}\enspace \omega t+b\left(t\right)\mathrm{cos}\enspace {\omega }_{i}t\enspace \mathrm{sin}\enspace \omega t\enspace \mathrm{d}t\end{equation} \tag{ 2.13 }$$

$$\begin{equation}{Y}_{i}=-\underset{-\tau /2}{\overset{\tau /2}{\int }}a\left(t\right)\mathrm{sin}\enspace {\omega }_{i}t\enspace \mathrm{cos}\enspace \omega t+b\left(t\right)\mathrm{sin}\enspace {\omega }_{i}t\enspace \mathrm{sin}\enspace \omega t\enspace \mathrm{d}t\end{equation} \tag{ 2.14 }$$

where τ denotes the gate time while X_i and Y_i parametrises the rotations about the x and y axes respectively, which realises the gate operation. We have neglected time-ordering in the quantum evolution to efficiently estimate the initial coefficients of the linear combinations of basis functions required for closed-loop optimisations. The effects due to time-ordering in the quantum evolution is discussed further in appendix A.

If the jth nucleus is the intended target, then the operations on all other nuclei (i ≠ j) are simply identity operations. Thus we have: for i ≠ j,

$$\begin{equation}{X}_{i}=-\underset{-\tau /2}{\overset{\tau /2}{\int }}a\left(t\right)\mathrm{cos}\enspace {\omega }_{i}t\enspace \mathrm{cos}\enspace {\omega }_{j}t+b\left(t\right)\mathrm{cos}\enspace {\omega }_{i}t\enspace \mathrm{sin}\enspace {\omega }_{j}t\enspace \mathrm{d}t=0\end{equation} \tag{ 2.15 }$$

$$\begin{equation}{Y}_{i}=-\underset{-\tau /2}{\overset{\tau /2}{\int }}a\left(t\right)\mathrm{sin}\enspace {\omega }_{i}t\enspace \mathrm{cos}\enspace {\omega }_{j}t+b\left(t\right)\mathrm{sin}\enspace {\omega }_{i}t\enspace \mathrm{sin}\enspace {\omega }_{j}t\enspace \mathrm{d}t=0\end{equation} \tag{ 2.16 }$$

and for i = j

$$\begin{equation}{X}_{j}=-\underset{-\tau /2}{\overset{\tau /2}{\int }}a\left(t\right)\mathrm{cos}\enspace {\omega }_{j}t\enspace \mathrm{cos}\enspace {\omega }_{j}t+b\left(t\right)\mathrm{cos}\enspace {\omega }_{j}t\enspace \mathrm{sin}\enspace {\omega }_{j}t\enspace \mathrm{d}t={X}_{T}\end{equation} \tag{ 2.17 }$$

$$\begin{equation}{Y}_{j}=-\underset{-\tau /2}{\overset{\tau /2}{\int }}a\left(t\right)\mathrm{sin}\enspace {\omega }_{j}t\enspace \mathrm{cos}\enspace {\omega }_{j}t+b\left(t\right)\mathrm{sin}\enspace {\omega }_{j}t\enspace \mathrm{sin}\enspace {\omega }_{j}t\enspace \mathrm{d}t={Y}_{T}\end{equation} \tag{ 2.18 }$$

We impose the restriction that ω = ω_j where ω_j is the transition frequency of the targeted qubit. Using these gate parametrisations, we introduce the following inverse Fourier transforms, where

$$\begin{equation}a\left(t\right)=\frac{1}{\sqrt{2\pi }}\underset{-\infty }{\overset{\infty }{\int }}a\left(\omega \right){\mathrm{e}}^{\mathbb{i}\omega t}\enspace \mathrm{d}\omega \end{equation} \tag{ 2.19 }$$

$$\begin{equation}b\left(t\right)=\frac{1}{\sqrt{2\pi }}\underset{-\infty }{\overset{\infty }{\int }}b\left(\omega \right){\mathrm{e}}^{\mathbb{i}\omega t}\enspace \mathrm{d}\omega \end{equation} \tag{ 2.20 }$$

As the signal is finite in time domain, we can enforce that a(t) and b(t) are zero outside $t\in \left[-\tau /2,\tau /2\right]$. This enables us to change the limits of the time integral to ±∞ and pass the time integral through the frequency integral. We also use the following identity

$$\begin{equation}\frac{1}{\sqrt{2\pi }}\underset{-\infty }{\overset{\infty }{\int }}{\mathrm{e}}^{-\mathbb{i}\omega t}\enspace \mathrm{d}t=\sqrt{2\pi }\delta \left(\omega \right)\end{equation} \tag{ 2.21 }$$

for further simplifications. For a(t) and b(t) to be real functions with well defined phases, we enforce

$$\begin{equation}{a}^{{\ast}}\left(\omega \right)=a\left(\omega \right)=a\left(-\omega \right)\end{equation} \tag{ 2.22 }$$

$$\begin{equation}{b}^{{\ast}}\left(\omega \right)=b\left(\omega \right)=b\left(-\omega \right)\end{equation} \tag{ 2.23 }$$

Using the above conditions, we arrived at the expressions where: for i ≠ j

$$\begin{equation}{X}_{i}=-\frac{\sqrt{2\pi }}{2}\left[a\left({\omega }_{i}+{\omega }_{j}\right)+a\left({\omega }_{i}-{\omega }_{j}\right)\right]=0\end{equation} \tag{ 2.24 }$$

$$\begin{equation}{Y}_{i}=\frac{\sqrt{2\pi }}{2}\left[b\left({\omega }_{i}+{\omega }_{j}\right)-b\left({\omega }_{i}-{\omega }_{j}\right)\right]=0\end{equation} \tag{ 2.25 }$$

and for i = j

$$\begin{equation}{X}_{j}=-\frac{\sqrt{2\pi }}{2}\left[a\left(2{\omega }_{j}\right)+a\left(0\right)\right]={X}_{T}\end{equation} \tag{ 2.26 }$$

$$\begin{equation}{Y}_{j}=\frac{\sqrt{2\pi }}{2}\left[b\left(2{\omega }_{j}\right)-b\left(0\right)\right]={Y}_{T}\end{equation} \tag{ 2.27 }$$

where X_T and Y_T are the intended angles of rotation on the x-axis and y-axis respectively.

2.4. Statistical model of gate errors

In reality, the control fields have noises and the frequencies of the qubits fluctuate slowly between the computational shots. Thus, the real control field can be written as

$$\begin{equation}{\gamma }_{i}{B}_{1}\left(t\right)=\left(1+{\epsilon}\right)\left[a\left(t\right)\mathrm{cos}\left(\left(\omega +\delta \right)t+\phi \right)+b\left(t\right)\mathrm{sin}\left(\left(\omega +\delta \right)t+\phi \right)\right]\end{equation} \tag{ 2.28 }$$

where , ϕ and δ are free parameters representing the amplitude, phase and frequency noises respectively. Let the target evolution operator for a single qubit gate on the jth nucleus within an N qubit cluster to be

$$\begin{equation}{G}_{T,x}={I}_{1}\otimes \cdots \otimes {I}_{j-1}\otimes {\mathrm{e}}^{-\mathbb{i}\left({I}_{j,x}{X}_{T}\right)}\otimes {I}_{j+1}\otimes \cdots \otimes {I}_{N}\end{equation} \tag{ 2.29 }$$

$$\begin{equation}{G}_{T,y}={I}_{1}\otimes \cdots \otimes {I}_{j-1}\otimes {\mathrm{e}}^{-\mathbb{i}\left({I}_{j,y}{Y}_{T}\right)}\otimes {I}_{j+1}\otimes \cdots \otimes {I}_{N}\end{equation} \tag{ 2.30 }$$

The actual evolution operator is defined as

$$\begin{equation}{G}_{A,x}=\underset{i=1}{\overset{N}{\bigotimes}}{\mathrm{e}}^{-\mathbb{i}\left({I}_{i,x}{X}_{i}\right)}\quad \mathrm{o}\mathrm{r}\quad {G}_{A,y}=\underset{i=1}{\overset{N}{\bigotimes}}{\mathrm{e}}^{-\mathbb{i}\left({I}_{i,y}{Y}_{i}\right)}\end{equation} \tag{ 2.31 }$$

The analytical expression for infidelity can be written as

$$\begin{equation}I=1-\frac{\mathrm{Tr}\left[{G}_{T}^{{\dagger}}{G}_{A}\right]}{\mathrm{Tr}\left[{G}_{T}^{{\dagger}}{G}_{T}\right]}\end{equation} \tag{ 2.32 }$$

This expression for infidelity compares an ideal gate to an experimental gate [60]. It is used as it requires minimal computation when compared to other figures of merit, and is equivalent to the usual performance function implemented in the GRAPE algorithm [40]. Assuming the gates are rotating about the x-axis, the infidelity expression can be written as

$$\begin{align}\hfill I& =1-\frac{1}{{2}^{N}}\left[\prod _{i\ne j}^{N}\mathrm{Tr}\left({\mathrm{e}}^{-\mathbb{i}\left({I}_{i,x}{X}_{i}\right)}\right)\right]\mathrm{Tr}\left[{\left({\mathrm{e}}^{-\mathbb{i}\left({I}_{j,x}{X}_{T}\right)}\right)}^{{\dagger}}{\mathrm{e}}^{-\mathbb{i}\left({I}_{j,x}{X}_{j}\right)}\right]\hfill \\ \hfill & =1-\left[\prod _{i\ne j}^{N}\mathrm{cos}\left(\frac{1}{2}{X}_{i}\right)\right]\left[\mathrm{cos}\left(\frac{{X}_{j}-{X}_{T}}{2}\right)\right]\hfill \end{align} \tag{ 2.33 }$$

Assuming small errors in the regime where X_i ≪ 1 for i ≠ j and $\left({X}_{j}-{X}_{T}\right)\ll 1$, we can then expand the above expression, and keep only the terms up to second order, which gives

$$\begin{equation}I\approx \frac{1}{8}\sum _{i\ne j}^{N}{X}_{i}^{2}+\frac{1}{8}{\delta }_{X}^{2}\end{equation} \tag{ 2.34 }$$

Repeating the same procedure for the rotations about the y-axis gives

$$\begin{equation}I\approx \frac{1}{8}\sum _{i\ne j}^{N}{Y}_{i}^{2}+\frac{1}{8}{\delta }_{Y}^{2}\end{equation} \tag{ 2.35 }$$

where δ_X = X_j − X_T and δ_Y = Y_j − Y_T

The control signal has a general form of

$$\begin{equation}a\left(\omega \right)=\sum _{n}{c}_{n}{a}^{\left(n\right)}\left(\omega \right)\end{equation} \tag{ 2.36 }$$

$$\begin{equation}b\left(\omega \right)=\sum _{n}{d}_{n}{b}^{\left(n\right)}\left(\omega \right)\end{equation} \tag{ 2.37 }$$

As such, the set of linear coefficients c_n and d_n that minimise the infidelity caused by the amplitude, phase and frequency noises can be solved.

We assumed the noise in the qubit frequencies δ, pulse amplitudes $\left(1+{\epsilon}\right)$ where is the fractional error of the pulse amplitudes and phase ϕ are described by Gaussian distributions centred at zero with their respective standard deviations of σ_δ, σ and σ_ϕ. The Gaussian distribution for frequency error is justified by experiments [42, 48]. While there is no clear evidence of either Gaussian or Lorentzian distribution for the amplitude error, we assumed Gaussian distribution and noting that reference [42] observed a Lorentzian distribution for the convolution of power fluctuations and pulse shape errors. The analytical expression for the average infidelity can be written as

$$\begin{align}\hfill \left\langle I\right\rangle & =\underset{-\infty }{\overset{\infty }{\int }}P\left(\delta ;{\sigma }_{\delta }\right)\underset{-\infty }{\overset{\infty }{\int }}P\left({\epsilon};{\sigma }_{{\epsilon}}\right)\underset{-\infty }{\overset{\infty }{\int }}P\left(\phi ;{\sigma }_{\phi }\right)I\left({\omega }_{1},\dots ,{\omega }_{N}\right)\mathrm{d}\delta \enspace \mathrm{d}{\epsilon}\enspace \mathrm{d}\phi \hfill \\ \hfill & =\frac{1}{8}\left[\underset{-\infty }{\overset{\infty }{\int }}p\left(\delta ;{\sigma }_{\delta }\right)\underset{-\infty }{\overset{\infty }{\int }}p\left({\epsilon};{\sigma }_{{\epsilon}}\right)\underset{-\infty }{\overset{\infty }{\int }}p\left(\phi ;{\sigma }_{\phi }\right){\delta }_{X}^{2}\enspace \mathrm{d}\delta \enspace \mathrm{d}{\epsilon}\enspace \mathrm{d}\phi +\sum _{i\ne j}^{N}\underset{-\infty }{\overset{\infty }{\int }}p\left(\delta ;{\sigma }_{\delta }\right)\underset{-\infty }{\overset{\infty }{\int }}p\left({\epsilon};{\sigma }_{{\epsilon}}\right)\right.\hfill \\ \hfill & \quad \left.{\times}\underset{-\infty }{\overset{\infty }{\int }}p\left(\phi ;{\sigma }_{\phi }\right){X}_{i}^{2}\enspace \mathrm{d}\delta \enspace \mathrm{d}{\epsilon}\enspace \mathrm{d}\phi \right]\hfill \end{align} \tag{ 2.38 }$$

A similar infidelity expression can be obtained for rotations about the y-axis. The minimum average gate infidelity is found when

$$\begin{equation}\frac{\partial \left\langle I\right\rangle }{\partial {c}_{n}}=0\end{equation} \tag{ 2.39 }$$

$$\begin{equation}\frac{\partial \left\langle I\right\rangle }{\partial {d}_{n}}=0\end{equation} \tag{ 2.40 }$$

are satisfied for each function in the expansion.

In this paper, we use frequency-shifted sinc functions as an ansatz for our control pulses in the frequency domain. Sinc functions were chosen as they represent pulses of finite duration in the time domain. As previously demonstrated in equations (2.36) and (2.37), the key property is the amplitude at very specific frequencies $\left({\omega }_{i},{\omega }_{j}\right)$. The interference of frequency-shifted sinc functions allow us to cancel the pulse amplitude at certain frequencies while at the same time amplifying other frequencies. The parametrisation of a(ω) and b(ω) is given by

$$\begin{equation}{a}^{\left(n\right)}\left(\omega \right)={f}^{\left(n\right)}\tau \left[\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{c}\left(\frac{\tau }{2}\left(\omega -{\mu }^{\left(n\right)}\right)\right)+\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{c}\left(\frac{\tau }{2}\left(\omega +{\mu }^{\left(n\right)}\right)\right)\right]\end{equation} \tag{ 3.1 }$$

$$\begin{equation}{b}^{\left(n\right)}\left(\omega \right)={g}^{\left(n\right)}\tau \left[\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{c}\left(\frac{\tau }{2}\left(\omega -{\nu }^{\left(n\right)}\right)\right)+\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{c}\left(\frac{\tau }{2}\left(\omega +{\nu }^{\left(n\right)}\right)\right)\right]\end{equation} \tag{ 3.2 }$$

where f⁽ⁿ⁾, g⁽ⁿ⁾ are the pulse amplitudes in the frequency domain, τ is the gate time, μ⁽ⁿ⁾ and ν⁽ⁿ⁾ are the frequency shifts and (n + 1) is the total number of basis functions used in the optimisation procedure.

If the frequency shifts for the nth solution are defined as

$$\begin{equation}{\mu }^{\left(n\right)}=n\left(\frac{2\pi }{\tau }\right)\end{equation} \tag{ 3.3 }$$

$$\begin{equation}{\nu }^{\left(n\right)}=n\left(\frac{2\pi }{\tau }\right)\end{equation} \tag{ 3.4 }$$

Then the above simply becomes a Fourier series, which is known to be a complete basis function. However, as seen in equations (3.1) and (3.2), the a⁽ⁿ⁾(0) and b⁽ⁿ⁾(0) terms would then be

$$\begin{equation}{a}^{\left(n\right)}\left(0\right)=2{f}^{\left(n\right)}\enspace \mathrm{s}\mathrm{i}\mathrm{n}\mathrm{c}\left(n\pi \right)\end{equation} \tag{ 3.5 }$$

$$\begin{equation}{b}^{\left(n\right)}\left(0\right)=2{g}^{\left(n\right)}\enspace \mathrm{s}\mathrm{i}\mathrm{n}\mathrm{c}\left(n\pi \right)\end{equation} \tag{ 3.6 }$$

where for n ⩾ 1, a⁽ⁿ⁾(0) and b⁽ⁿ⁾(0) terms will always be 0. Hence, without any contributions from the a⁽ⁿ⁾(0) and b⁽ⁿ⁾(0) terms, generating an optimal pulse from only the ${a}^{\left(n\right)}\left(2{\omega }_{j}\right)$ and ${b}^{\left(n\right)}\left(2{\omega }_{j}\right)$ terms would require extremely large amplitudes.

One possible method to overcome this complication is to use Kadec's 1/4 theorem, where an additional small frequency shift is introduced to the Fourier series while retaining the completeness of the Fourier series. It was shown that when the additional frequency shift is bounded by a maximum value of $0.25\left(2\pi /\tau \right)$, the inequality will still be able to generate a continuous set of sinc basis functions [61]. In this paper, we have chosen the upper bound to be $0.2\left(2\pi /\tau \right)$ as it enables us to resolve better solutions for the pulse amplitudes (see figure 2).

The frequency shifts are then redefined as

$$\begin{equation}{\mu }^{\left(n\right)}=\left(n+\frac{1}{5}\right)\left(\frac{2\pi }{\tau }\right)\end{equation} \tag{ 3.7 }$$

$$\begin{equation}{\nu }^{\left(n\right)}=\left(n+\frac{1}{5}\right)\left(\frac{2\pi }{\tau }\right)\end{equation} \tag{ 3.8 }$$

As a result, there are many solutions in ${a}^{\left(n\right)}\left(\omega \right)$ and ${b}^{\left(n\right)}\left(\omega \right)$ that we can consider and these solutions form a linear basis for the construction of optimal pulse functions. These optimal pulse functions are found by determining the linear coefficients that minimise the effects of the pulse errors on average.

Zoom In Zoom Out Reset image size

As an example, consider a two-qubit system which consists of ¹⁵N and ¹³C nuclear spins. Their respective hyperfine interactions are given by A_N ≈ 3 MHz [52] and A_C ≈ 0.413 MHz [4]. The background static magnetic field is chosen to be B₀ = 0.62 T [4]. Using these parameters, we demonstrate the optimisation procedure for an X gate targeted at the ¹³C nuclear spin $\left({X}_{2}\right)$. Since we are performing rotations only in the x-axis/y-axis, we assumed there are no mixed signals in the pulse and thus, there are no contributions from the b(ω)/a(ω) components. Hence, the equations that are satisfied by the solutions are given by

$$\begin{equation}{X}_{1}=-\frac{\sqrt{2\pi }}{2}\left[a\left({\omega }_{\mathrm{N}}+{\omega }_{\mathrm{C}}\right)+a\left({\omega }_{\mathrm{N}}-{\omega }_{\mathrm{C}}\right)\right]=0\end{equation} \tag{ 3.9 }$$

$$\begin{equation}{X}_{2}=-\frac{\sqrt{2\pi }}{2}\left[a\left(2{\omega }_{\mathrm{C}}\right)+a\left(0\right)\right]=\pi \end{equation} \tag{ 3.10 }$$

where ω_N = A_N + γ_NB₀ and ω_C = A_C + γ_CB₀. The first step of the optimisation routine is to generate a set of basis functions that minimise the infidelity landscape in the absence of control errors against the pulse amplitude, f⁽ⁿ⁾ and g⁽ⁿ⁾ in the frequency domain as described by equation (2.33). This step minimises crosstalk between the qubits.

The search range for the pulse amplitudes is constrained by the design of our MW/RF system and this corresponds to the shortest gate time that we can perform for our qubit gate operations. In the time domain, the pulse amplitude, a(t) is determined by γ_iB₁(t) where γ_i is the gyromagnetic ratio of the nuclear spin qubits and B₁(t) is the maximum amplitude of the oscillating MW/RF magnetic field. Implementing the typical values used in an experiment, the pulse amplitude a(t) is thus limited to a maximum value of approximately 25 Mrad s⁻¹. The corresponding search range for f⁽ⁿ⁾ and g⁽ⁿ⁾ is then arbitrarily set to ${f}^{ \left(n\right)},{g}^{\left(n\right)}\in \left[-5,5\right]$ and noting that the amplitude in the time domain has an additional factor of $2\sqrt{2\pi }$ from the inverse Fourier transform of the sinc basis functions.

Similarly, the set of basis functions for a two-qubit CZ gate is generated using the same procedure. In this case, the pulse amplitude is limited to approximately 80 Mrad s⁻¹. The search range for f⁽ⁿ⁾ is bounded in the region of ${f}^{\left(n\right)}\in \left[0.5,15\right]$ as the intrinsic infidelity expression used to describe two-qubit gate operations, equation (D.1) is symmetrical. Using f⁽ⁿ⁾ = 0, the optimisation will be stuck in a local minimum and clearly, f⁽ⁿ⁾ = 0 corresponds to no physical pulse. Thus, in order to find sensible solutions, we shifted the initial search boundary by 0.5 to stimulate the optimisation to find another local minimum. The equations that are satisfied by the solutions are given by equations (C.9)–(C.12).

As seen in table E1, the pulse amplitudes are dependent on the angle of rotations and gate time. This is consistent with the formulation of our sinc basis functions where larger amplitudes are expected to generate greater angle of rotations for a given frequency shift and gate time. We are also able to resolve solutions with lower infidelity for smaller angle of rotations at shorter gate times as the pulse amplitudes are smaller. We observed a trend, shown in tables E1 and E2, where at short gate times, higher order basis functions have solutions which correspond to the lower or upper bound of the allowed values. These pulse amplitudes minimise the infidelities within the search range imposed on them as set by physics. While some of the pulse amplitudes correspond to the maximum or minimum allowed values, the critical factor in this optimisation routine is the linear combinations of these generated basis functions as discussed in the following section.

Applying this optimal control method, we optimise the infidelities of an X gate, Hadamard gate and a CZ gate due to the control errors. We approximated the standard deviations of the Gaussian phase, amplitude and frequency noises to be σ_ϕ = 1/2π kHz, σ = 10⁻³, σ_δ = 1 kHz for nuclear spins and σ_δ = 27.5 kHz for the electron spin [48]. Phase noises are excluded from the calculation of average infidelities for a CZ gate as we are performing a 2π pulse. The first order effects due to the phase errors are negligible.

As shown in figure 3, the average infidelities for an X gate, a π/2 rotation about the y-axis and a Hadamard gate fluctuate with a single basis function, since the function parameters depend on the local infidelity landscape during the initial basis function computation. However, the optimised linear combinations of two or more basis functions yield infidelities of approximately 10⁻⁶ for an X gate and a Hadamard gate. Infidelities up to 10⁻⁷ can be achieved for a π/2 rotation about the y-axis. The infidelities are monotonically decreasing with increasing number of basis functions and only 3 basis functions are required for the infidelities to converge to the optimisation limits. Thus, this demonstrates the capability of this optimal control method, and allow us to systematically assess the degree of convergence to the optimisation limits. For a system with two qubits, using more than 3 basis functions does not significantly lower the infidelities. In general, we expect to achieve infidelities on this order (of 10⁻⁶) for any single-qubit operation.

Zoom In Zoom Out Reset image size

Given that this method can achieve low infidelities, the main objective is to shorten the gate time without significantly affect the infidelity, in order to minimise the effects due to decoherence. We analyse the overall amplitudes of the respective linear combinations of the sinc basis functions in order to find the minimum gate time for which the pulse amplitude can still be practically generated. Using estimated constraints described in the previous section, the maximum threshold for the amplitude is set to be log₁₀(25) ≈ 1.40.

Figure 4 depicts the overall amplitudes of the basis functions in the linear combinations for different basis size and gate times. At some gate times, the linear combinations have larger amplitudes when more basis functions are used. Since the optimal control method optimises the linear combinations of the basis functions such that the infidelities are monotonically decreasing, therefore, if the amplitude in the initial basis function is large, then it may result in larger amplitudes for the subsequent basis functions when the basis size is increased. Based on the generated data, we determined that a 0.5  μs X gate with an infidelity of ∼10⁻⁶ can be achieved with just two basis functions (figure 4(b)). Despite being able to perform a 0.5  μs X gate, we need to take into account various gate times required for different types of rotations where slower gates are required for bigger angle of rotations. Thus, on average, a conservative estimate for the fastest single-qubit gate that we are able to perform with ∼10⁻⁶ infidelity is approximately 1  μs.

Zoom In Zoom Out Reset image size

Applying our method now to two-qubit gates, the maximum threshold for the amplitude of a two-qubit CZ gate is set to be log₁₀(80) ≈ 1.9. From figure 5, on average, gate infidelities of 10⁻⁴ ∼ 10⁻⁶ can be achieved with 4 or more basis functions for gate times τ ⩾ 1  μs. The calculated infidelities are monotonically decreasing with increasing basis size for τ ⩾ 1  μs. Monotonicity was not fully demonstrated for τ < 1  μs due to numerical integration errors. These errors are exacerbated for τ < 1  μs because of many linear combinations of basis functions with extremely large amplitudes at short gate times, resulting in a high oscillatory integrand. On average, more than 8 basis functions are required for the convergence of these infidelities to their respective optimisation limits. Based on the optimisations shown in figures 5 and 6, the fastest two-qubit CZ gate that we can perform with an infidelity of ∼10⁻⁶ is 1  μs and it requires 6 linear combinations of basis functions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Additionally, we have analysed the infidelities of a CZ gate for a range of electron spin linewidths that may occur in materials with different densities of ¹³C. Figures 7(a) and (b) depicts the CZ gate infidelities with σ_δ = 2 × 27.5 kHz and 4 × 27.5 kHz, respectively. Meanwhile, figure 7(c) demonstrates convergence of CZ gate infidelities with τ = 2.25  μs occurring at different number of basis functions for different electron spin linewidths, with larger linewidths requiring larger number of basis functions. Overall, we are still able to obtain infidelities less than 10⁻⁴ across the range of linewidths, but at a cost of longer gate time and larger number of basis functions.

Zoom In Zoom Out Reset image size

We now turn our attention to the final step of the optimisation process, which is to incorporate time-ordering in the quantum evolution into our optimisation routine. As a proof of principle demonstration, we have performed this optimisation process on an X gate (see appendix A). Through further optimisations where the quantum evolution is evaluated numerically using a finite-difference method, we found that an X gate can theoretically achieve infidelities of 10⁻⁵ or lower. We believe that with the addition of more advanced numerical optimisation techniques, then higher number of basis functions can be included as well as the treatment of two-qubit gates. Thus, we are confident that we can reduce the control errors to the order of 10⁻⁵ or better for all gates. However, it is not necessary to involve the numerical optimisation techniques if we are performing experiments on the diamond quantum processors as the same feedback-based optimisation routine will occur.

Now, we must assess the implications for neglecting the hyperfine field misalignments when attaining these infidelities. From appendix B.2, the hyperfine field misalignment only affects the two-qubit gate. While we are able to minimise the gate infidelities due to control errors to the order of 10⁻⁵ or 10⁻⁶, the inclusion of the correction terms from the hyperfine field misalignment introduced an additional infidelity on the order of 10⁻³. This additional infidelity arises from the projection of the applied magnetic field due to the rotation of the nuclear spin coordinates in the computational subspace. One way to overcome this effective field is to include radiowave pulses which are resonant with the nuclear spins in the m_s = 0 state in conjunction with the microwave pulses which drive the electron spin transitions. This dual application of radio and microwave pulses is beyond the scope of this work and should be pursued in future work. There is no physical reason why these errors due to hyperfine field misalignment cannot be reduced to the level of our control errors.

As mentioned in section 1, the nuclear spin qubits undergo pure dephasing due to the relaxation of the electron spin. Thus, the coherence time of the qubits is bounded by the relaxation time of the electron spin. Here, we introduce the master equation which is also known as the Lindblad equation [62, 63]

$$\begin{equation}\dot {\rho }\left(t\right)=\mathcal{L}\left(\rho \right)\equiv -\frac{\mathrm{i}}{\hslash }\left[H,\rho \left(t\right)\right]+\sum _{m=1}^{N}\left({L}_{m}\rho \left(t\right){L}_{m}^{{\dagger}}-\frac{1}{2}{L}_{m}^{{\dagger}}{L}_{m}\rho \left(t\right)-\frac{1}{2}\rho \left(t\right){L}_{m}^{{\dagger}}{L}_{m}\right)\end{equation} \tag{ 5.1 }$$

where the sum over m is the summation of the decoherence mechanism over the individual nuclear spins. Instead of solving the differential equations, we can express the Lindblad equation as a vectorized density matrix [64]

$$\begin{equation}\dot {\rho }=\left(\mathbf{H}+\mathbf{G}\right)\rho \end{equation} \tag{ 5.2 }$$

where G is the decoherent part of the Lindblad equation with the form of

$$\begin{equation}\mathbf{G}=\sum _{m=1}^{N}\left({\overline{L}}_{m}\otimes {L}_{m}-\frac{1}{2}\mathbb{I}\otimes \left({L}_{m}^{{\dagger}}{L}_{m}\right)-\frac{1}{2}\left({\overline{L}}_{m}^{{\dagger}}{\overline{L}}_{m}\right)\otimes \mathbb{I}\right)\end{equation} \tag{ 5.3 }$$

and

$$\begin{equation}\mathbf{H}=-\frac{\mathrm{i}}{\hslash }\left(\overline{H}\otimes \mathbb{I}-\mathbb{I}\otimes H\right)\end{equation} \tag{ 5.4 }$$

The overline denotes the complex conjugate, † is the adjoint, H is the Hamiltonian of the qubit system and $\mathbb{I}$ is the 2^N × 2^N identity matrix where N is the number of qubits in the system. For a time dependent Hamiltonian, the Linblad equation is given by

$$\begin{equation}\rho \left(t\right)=\mathrm{exp}\left[\underset{-\tau /2}{\overset{\tau /2}{\int }}\mathbf{H}\enspace \mathrm{d}t+\mathbf{G}t\right]\rho \left(0\right)\end{equation} \tag{ 5.5 }$$

The Lindblad operator L describing the dephasing of the nuclear spin qubits induced by a random electron spin flip can be written as

$$\begin{equation}{L}_{m}=\sqrt{\frac{1}{2{T}_{2}}}{\sigma }_{z,m}\end{equation} \tag{ 5.6 }$$

where the nuclear spin T₂ is defined by the relaxation time T₁ of the electron spin (1.8 ms) [48] and σ_z,m is the Pauli matrix for the z component of each nuclear spin.

To assess the effects of decoherence due to dephasing, we consider an example with a perfect X gate. The aim now is to solve for ρ(t) and calculate its state fidelity defined as

$$\begin{equation}F=\frac{\mathrm{Tr}\left[{\rho }_{\mathrm{I}}^{{\dagger}}\rho \right]}{\mathrm{Tr}\left[{\rho }_{\mathrm{I}}^{{\dagger}}{\rho }_{\mathrm{I}}\right]}\end{equation} \tag{ 5.7 }$$

where ρ_I is the ideal density matrix without the effects of decoherence and ρ is the simulated density matrix of the system.

As shown in figure 8, the errors in both cases have magnitudes of approximately 10⁻³ ∼ 10⁻⁴. These errors are much larger than the errors caused by the effects of frequency, phase and amplitude noises as demonstrated in section 4 and the decoherence errors increase with longer gate times. The decoherence errors also have the same order of magnitude $\left(\sim , 1{0}^{-3}\right)$ as the errors introduced by hyperfine field misalignment on CZ gates (see appendix B.2). Since we cannot remove the electron-nuclear spin coupling as it is required for the selective operations, this leads to the conclusion that for a given electron relaxation time T_1,e, the only solution to improve the gate fidelity is to make the gate operations faster.

Zoom In Zoom Out Reset image size

In this section, we set a benchmark for the optimal performance of a diamond quantum processor by simulating quantum algorithms. The benchmark will provide us with insights into the limits of the processor, which can be used to aid the design and comparison of the device in the near future.

Owing to our expectations that with modifications to the control system, errors due to the hyperfine field misalignment can be reduced to the level of control errors, which are small compared to the effects of decoherence. Hence, building on the results from section 5 and appendix B.2, we can simulate quantum algorithms on the diamond quantum computer by only considering the effects of decoherence. The key metrics for simulation will be the error probability and the total computational time (ignoring initial loading time of the computational control systems). These two metrics are chosen since computational time is the primary resource and the error rate is the key quality of a quantum computer. We have chosen to compute QFTs for simulation as QFT is widely used in quantum algorithms [65]. A further motivation is the similarity to the algorithms used for enhanced quantum sensing using a register of nuclear spin qubits. The initial state was chosen in a way such that an output state of $\left\vert 001\right\rangle$ (3 qubit QFT) and $\left\vert 00\enspace 001\right\rangle$ (5 qubit QFT) will be the only outcome with a probability of 1 (figure 9).

Zoom In Zoom Out Reset image size

The simulations of QFT on the diamond quantum computer are done using equation (5.5) with electron relaxation time of T_1,e ≈ 1.8 ms. These simulations will be iteratively solved for multiples of the fastest single and two-qubit gate times by only considering the effects of decoherence as the effects due to hyperfine field misalignments have the same order of magnitude and only affect CZ gates. The total number of pulses required for 3 qubits (QFT3) and 5 qubits (QFT5) QFTs after the decomposition into rotations about the x and y axis and CZ gate are 75 and 195 for QFT3 and QFT5 respectively (appendix F).

The total computation time on a diamond quantum computer can be broken down into shot time and initialisation/readout time. Shot time can be regarded as the total duration of the pulses required for an experiment. Assuming single-qubit gate times of 1 μs and CZ gate time of 1  μs, the optimal pulse duration for QFT3 and QFT5 are approximately 75  μs and 195  μs, respectively. For the initialisation/readout time (single-shot readout), we are going to apply M number of readout cycles per qubit. Thus, the total time for initialisation/readout is given by

$$\begin{equation}\mathrm{I}\mathrm{n}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{s}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}/\mathrm{R}\mathrm{e}\mathrm{a}\mathrm{d}\mathrm{o}\mathrm{u}\mathrm{t}\enspace \mathrm{T}\mathrm{i}\mathrm{m}\mathrm{e}=n{\times}M{\times}{t}_{\mathrm{c}}\end{equation} \tag{ 6.1 }$$

where n is the number of qubits, M is the number of readout cycles applied per qubit and t_c is the time per cycle. Time per cycle is based on the time of the optical pulse required for readout (t_opt) and the time of the microwave pulse (t_mw) required to perform the CNOT gate for repetitive measurements. These two time quantities are approximately 1 μs each. M is chosen to be 500 as it has the same magnitude as other numbers of repetition which achieve an initialisation fidelity of 0.99 given a specific relative shift of the initialisation threshold [4]. The time required for a single-shot readout is given by

$$\begin{equation}{T}_{\text{QFT}3}=75{\times}1{0}^{-6}+3{\times}500{\times}2{\times}1{0}^{-6}\approx 0.0031\enspace \mathrm{s}\end{equation} \tag{ 6.2 }$$

$$\begin{equation}{T}_{\text{QFT}5}=195{\times}1{0}^{-6}+5{\times}500{\times}2{\times}1{0}^{-6}\approx 0.0052\enspace \mathrm{s}\end{equation} \tag{ 6.3 }$$

As depicted in figure 10(a), simulation of QFT3 is able to achieve higher fidelity (≈0.964) than QFT5 (≈0.855) when they are simulated using the optimal gate time of 1  μs for single-qubit gates and 1  μs for CZ gate as QFT3 has smaller circuit size. Implementing a lesser number of gates will introduce less errors during the evolution of the quantum states, thus an algorithm with a smaller circuit size will achieve greater output state fidelity. The larger decay constant in the fitted model of QFT5 indicates that it is important to perform optimal control on the pulses to obtain the fastest gate possible with the lowest infidelity, as well as optimising our circuit size. Consequently, this will give us the best result when an experiment is performed using an actual diamond quantum computer.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To simulate the performance of a diamond quantum computer over time, we first create a binomial distribution of probabilities with the respective optimal values for QFT3 and QFT5 simulated on a diamond quantum computer. Then, we simulate the probability distribution for a range of shots using Markov chain Monte Carlo method. Markov chain Monte Carlo allows us to approximate the probability distribution of the fidelities through random sampling. With short computation time per shot for QFT3 and QFT5, the simulations of QFT3 and QFT5 on a diamond quantum computer are able to converge to their respective optimal values of 0.964 and 0.855 in less than 10 s (figure 11). These results can be used as a benchmark for comparison with other quantum computing architectures in the future.

In summary, we have presented a complete model of a diamond quantum processor, the gate operations and their implementations. We have developed a semi-analytical optimal control method which theoretically produces fast, high fidelity gates with control errors that are less than those introduced by decoherence and hyperfine field misalignment. This optimal control method has three steps. First, we generate a complete semi-analytical basis of pulses that minimises cross-talk between qubits in absence of time-ordering in the quantum evolution. Second, we generate linear combination of these basis functions to minimise the control field errors. Third, we further optimise the pulses to account for the effects of time-ordering in the quantum evolution. We used frequency-shifted sinc functions as an ansatz for the control pulses. We have demonstrated the first two steps of this optimal control method on a Hadamard gate and a CZ gate with an additional third step on an X gate. The simulated performance of a diamond quantum computer shows promising results where it can perform fast computations with low error probability. Our results will aid the design and the development of diamond quantum computers and enhanced quantum sensors.

A future extension should include more advanced numerical optimisation techniques to optimise the pulses in the presence of time-ordering in the quantum evolution. However, if the optimisation is being performed experimentally, this step may not be required. Instead, we can implement a feedback control system as an optimal control method which tunes the pulse parameters based upon the output of a physical diamond quantum computer. Subsequent work should also include the introduction of a more complex optimal control technique to mitigate the errors of the hyperfine field misalignment through the application of simultaneous radiofrequency and microwave pulses.

We acknowledge support from the Australian Research Council (DP170103098 and DE170100169) as well as the Australian National University, University of Canberra and Charles Sturt University through the Discovery Translation Fund 2.0, managed by ANU Connect Ventures.

In section 2.3, we ignored the time-ordering in the quantum evolution to simplify the unitary operators in terms of rotation about the x-axis or y-axis. In this appendix, we examine the effects due to time-ordering in the quantum evolution on the gate fidelity.

We consider the case where b(t) of equation (2.11) is zero. The nuclear spin Hamiltonian parameterised with the control error parameters δ, and ϕ as described in section 4 is given by

$$\begin{equation}{H}_{\text{single}}=-\sum _{i}a\left(t\right)\left(1+{\epsilon}\right)\left[{I}_{i,x}\enspace \mathrm{cos}\left({\omega }_{i}t\right)\mathrm{cos}\left(\left[\omega +\delta \right]t+\phi \right)+{I}_{i,y}\enspace \mathrm{sin}\left({\omega }_{i}t\right)\mathrm{cos}\left(\left[\omega +\delta \right]t+\phi \right)\right]\end{equation} \tag{ A.1 }$$

This Hamiltonian is the same as the Hamiltonian described in section 2.2. Using a finite-difference method, we approximate the Hamiltonian as being time-independent over the short period of Δt. The quantum evolution of the Hamiltonian can be written as

$$\begin{equation}{U}_{\text{fd}}\left(\tau ,\delta ,{\epsilon},\phi \right)\approx \prod _{m=0}^{N-1}{\mathrm{e}}^{-\mathrm{i}H\left({t}_{m}\right){\Delta}t}\end{equation} \tag{ A.2 }$$

in which the product is a time-ordered product, t_n = nΔt and Δt = τ/N where N is the number of time steps and τ is the gate time.

We use the pulse, a(t) generated from section 4. The functional form of a(t) can be obtained by performing an inverse Fourier transform on a(ω). It is given by

$$\begin{equation}a\left(t\right)=\sum _{n=1}2\sqrt{2\pi }{c}_{n}\enspace \mathrm{cos}\left(\frac{2\pi t\left[\left(n-1\right)\mu \right]}{\tau }\right)\end{equation} \tag{ A.3 }$$

where c_n is the overall amplitudes of the linear combinations for the nth basis functions as estimated by optimising over the control errors and μ = 0.2 is the shift of the basis functions. We further optimise a(t) to account the effects of time-ordering by introducing a change d_n to the coefficients where

$$\begin{equation}a\left(t\right)=\sum _{n=1}2\sqrt{2\pi }\left({c}_{n}+{d}_{n}\right)\mathrm{cos}\left(\frac{2\pi t\left[\left(n-1\right)\mu \right]}{\tau }\right)\end{equation} \tag{ A.4 }$$

The best parameters for d_n are obtained by minimising the average infidelity as given by

$$\begin{equation}\left\langle I\right\rangle =\underset{-\infty }{\overset{\infty }{\int }}P\left(\delta ;{\sigma }_{\delta }\right)\underset{-\infty }{\overset{\infty }{\int }}P\left({\epsilon};{\sigma }_{{\epsilon}}\right)\underset{-\infty }{\overset{\infty }{\int }}P\left(\phi ;{\sigma }_{\phi }\right)\left[1-\frac{\mathrm{Tr}\left({U}_{\text{fd}}^{{\dagger}}\cdot {U}_{\text{perfect}}\right)}{\mathrm{Tr}\left({U}_{\text{fd}}^{{\dagger}}\cdot {U}_{\text{fd}}\right)}\right]\enspace \mathrm{d}\delta \enspace \mathrm{d}{\epsilon}\enspace \mathrm{d}\phi \end{equation} \tag{ A.5 }$$

where the probability distributions P for all control parameters are of Gaussian.

Applying this method to an X gate of 1 μs with 2 basis functions where c₁ ≈ −1.20 and c₂ ≈ 0.79, we found that the infidelities have converged for N > 4000. For the further optimisation of a(t), we performed a simple grid search to obtain the best parameters for d_n. The parameters of d₁ ≈ −0.20 and d₂ ≈ −1.17 resulted in an averaged infidelity of approximately 2.90 × 10⁻⁵, which is comparable to the infidelities obtained in section 4 and is still below the errors introduced by decoherence (figure A1).

Zoom In Zoom Out Reset image size

Calculations for a two-qubit CZ gate are not demonstrated here due to sheer size of the search space considering it requires at least 6 basis functions to achieve an infidelity of approximately 10⁻⁶. It is computationally expensive to perform these calculations without implementing more advanced algorithms for closed-loop optimisations. These calculations will be considered in future extensions.

Since we are able to optimise an X gate with two linear combinations with infidelities below the decoherence limit, there is no reason why we could not achieve comparable infidelities for more linear combinations. Thus, while the time-ordering in the quantum evolution introduces additional infidelities to the gate operations, they can still be optimised to a level below the errors introduced by decoherence.

In section 2.2, we applied the secular approximation and ignored the hyperfine field misalignment to simplify the nuclear spin Hamiltonian of the quantum processor. In this appendix, we examine the validity of these two approximations, identified their leading order corrections and their effects on gate fidelity.

B.1. Secular approximation

In this subsection, we examine the validity of the secular approximation. The full electron–nuclear Hamiltonian can be divided into two terms

$$\begin{equation}H={H}^{\left(a\right)}+{H}^{\left(b\right)}\end{equation} \tag{ B.1 }$$

where

$$\begin{align}\hfill {H}^{\left(a\right)}& =D{S}_{z}^{2}+{\gamma }_{\mathrm{e}}{B}_{0}{S}_{z}+{S}_{z}\left({A}_{zz}{I}_{z}+{A}_{xz}{I}_{x}+{A}_{yz}{I}_{y}\right)-{\gamma }_{n}{B}_{0}{I}_{z}\hfill \\ \hfill {H}^{\left(b\right)}& ={S}_{x}\left({A}_{xx}{I}_{x}+{A}_{yx}{I}_{y}+{A}_{zx}{I}_{z}\right)+{S}_{y}\left({A}_{xy}{I}_{x}+{A}_{yy}{I}_{y}+{A}_{zy}{I}_{z}\right)\hfill \end{align} \tag{ B.2 }$$

H^(a) is the secular component that contains terms with S_z whereas H^(b) is the non-secular component, containing S_x and S_y. For large magnetic fields aligned with the NV axis, the non-secular terms are small compared to the secular terms. In this Hamiltonian, the nuclear spin coordinates are co-aligned with the NV coordinates.

Treating the non-secular terms as a perturbation, we can apply perturbation theory to collapse the full Hamiltonian into effective nuclear spin Hamiltonians corresponding to the different electron spin projections. The nuclear spin Hamiltonian of the m_s = n electron spin projection, correct to second order in the non-secular terms can be written as

$$\begin{equation}{H}_{n}=\left\langle n\left\vert {H}^{\left(a\right)}\right\vert n\right\rangle +\sum _{m\ne n}\frac{\left\langle n\left\vert {H}^{\left(b\right)}\right\vert m\right\rangle \left\langle m\left\vert {H}^{\left(b\right)}\right\vert n\right\rangle }{{E}_{n}-{E}_{m}}\end{equation} \tag{ B.3 }$$

It follows that

$$\begin{equation}{H}_{I,{m}_{\mathrm{s}}=0}=-{\gamma }_{n}{B}_{0}{I}_{z}-\frac{{\gamma }_{\mathrm{e}}{B}_{0}}{{D}^{2}-{\gamma }_{\mathrm{e}}^{2}{B}_{0}^{2}}\left[\left({A}_{yy}{A}_{zx}-{A}_{yx}{A}_{zy}\right){I}_{x}+\left({A}_{xx}{A}_{zy}-{A}_{xy}{A}_{zx}\right){I}_{y}+\left({A}_{xy}^{2}-{A}_{xx}{A}_{yy}\right){I}_{z}\right]\end{equation} \tag{ B.4 }$$

and

$$\begin{align}\hfill {H}_{I,{m}_{\mathrm{s}}=-1}& =D-{\gamma }_{\mathrm{e}}{B}_{0}-\left({A}_{zz}{I}_{z}+{A}_{xz}{I}_{x}+{A}_{yz}{I}_{y}\right)-{\gamma }_{n}{B}_{0}{I}_{z}-\frac{1}{2\left(D-{\gamma }_{\mathrm{e}}{B}_{0}\right)}\hfill \\ \hfill & \quad {\times}\left[\left({A}_{yy}{A}_{zx}-{A}_{yx}{A}_{zy}\right){I}_{x}+\left({A}_{xx}{A}_{zy}-{A}_{xy}{A}_{zx}\right){I}_{y}+\left({A}_{xy}^{2}-{A}_{xx}{A}_{yy}\right){I}_{z}\right]\hfill \end{align} \tag{ B.5 }$$

In equation (B.5), the scalar terms D − γ_eB₀ do not affect the nuclear spin so they can be neglected from the nuclear spin Hamiltonian. However, they do perturb the electron spin transition frequencies. The terms in the square brackets of equations (B.4) and (B.5) are the lowest order corrections to the nuclear spin Hamiltonians. For m_s = −1, the above can be rewritten in terms of ${ \overrightarrow {A}}^{\prime }\cdot \overrightarrow {I}$ where ${ \overrightarrow {A}}^{\prime }$ is the corrected hyperfine field. Thus, the correction terms add a quantitative difference to the values of ${ \overrightarrow {A}}_{i,z}$, but not a functional difference that affects the validity of our model. However, the corrections do represent a change for the m_s = 0 level as there is no equivalent terms used in the model. They are equivalent to an effective magnetic field misaligned with the NV axis that introduces error to our model. The error introduced by the neglect of these terms will be examined in the next subsection alongside with the analysis of hyperfine field misalignments.

B.2. Effects of hyperfine field misalignment

In this subsection, we examine the effects due to the hyperfine field misalignment. After applying the secular approximation, the nuclear spin Hamiltonian (equation (2.1)) in the auxiliary and computational subspace are

$$\begin{equation}{H}_{I,{m}_{\mathrm{s}}=0}=-\sum _{i}{\gamma }_{i}\left[{I}_{i,z}{B}_{0}+{I}_{i,x}{B}_{1}\left(t\right)\right]\end{equation} \tag{ B.6 }$$

$$\begin{equation}{H}_{I,{m}_{\mathrm{s}}=-1}=-\sum _{i}{ \overrightarrow {A}}_{i,z}\cdot { \overrightarrow {I}}_{i}-{\gamma }_{i}\left[{I}_{i,z}{B}_{0}+{I}_{i,x}{B}_{1}\left(t\right)\right]\end{equation} \tag{ B.7 }$$

To diagonalise the nuclear spin Hamiltonian of the computational subspace (i.e. for m_s = −1), we performed a coordinate transform of the nuclear spin operators. We did this by a rotation defined by the angles

$$\begin{equation}\mathrm{tan}\enspace {\phi }_{i}=\frac{{A}_{i,xz}}{{A}_{i,yz}}\end{equation} \tag{ B.8 }$$

$$\begin{equation}\mathrm{tan}\enspace {\theta }_{i}=\frac{{\left({A}_{i,xz}^{2}+{A}_{i,yz}^{2}\right)}^{1/2}}{{A}_{i,zz}+{\gamma }_{i}{B}_{0}}\end{equation} \tag{ B.9 }$$

The transformed Hamiltonians for the auxiliary m_s = 0 and computational m_s = −1 subspaces are then

$$\begin{align}\hfill {H}_{I,{m}_{\mathrm{s}}=0}& =-\sum _{i}{\gamma }_{i}{B}_{0}\left({I}_{i,z}^{\prime }\enspace \mathrm{cos}{\theta }_{i}-{I}_{i,y}^{\prime }\mathrm{sin}{\theta }_{i}\right)-\sum _{i}{\gamma }_{i}{B}_{1}\left(t\right)\left[\left({I}_{i,x}^{\prime },\enspace ,\mathrm{cos},{\phi }_{i}\right.\right.\hfill \\ \hfill & \left.\left.+,{I}_{i,y}^{\prime },\enspace ,\mathrm{sin},{\phi }_{i},\mathrm{cos},\enspace ,{\theta }_{i},+,{I}_{i,z}^{\prime },\enspace ,\mathrm{sin},{\theta }_{i},\mathrm{sin},{\phi }_{i}\right)\right]\hfill \end{align} \tag{ B.10 }$$

$$\begin{equation}{H}_{I,{m}_{\mathrm{s}}=-1}=-\sum _{i}{\omega }_{i}{I}_{i,z}^{\prime }-\sum _{i}{\gamma }_{i}{B}_{1}\left(t\right)\left[\left({I}_{i,x}^{\prime }\enspace \mathrm{cos}\enspace {\phi }_{i}+{I}_{i,y}^{\prime }\enspace \mathrm{sin}\enspace {\phi }_{i}\enspace \mathrm{cos}\enspace {\theta }_{i}+{I}_{i,z}^{\prime }\enspace \mathrm{sin}\enspace {\theta }_{i}\enspace \mathrm{sin}\enspace {\phi }_{i}\right)\right]\end{equation} \tag{ B.11 }$$

where ${\omega }_{i}={\left({\left[{A}_{i,zz}+{\gamma }_{i}{B}_{0}\right]}^{2}+{A}_{i,xz}^{2}+{A}_{i,yz}^{2}\right)}^{1/2}$ is the transition frequency of the ith nucleus in the m_s = −1 subspace. We assume that only nuclei with hyperfine fields nearly aligned to the NV axis are chosen. We then perform small angle approximations on θ_i and simplify the expressions by only including the terms which are zero-order in θ_i. Finally, undoing rotations about the z-axis gives us the zero-order nuclear spin Hamiltonians

$$\begin{equation}{H}_{I,{m}_{\mathrm{s}}=0}^{\left(0\right)}\approx -\sum _{i}{\gamma }_{i}{B}_{0}{I}_{i,z}^{{\prime\prime}}-\sum _{i}{I}_{i,x}^{{\prime\prime}}{\gamma }_{i}{B}_{1}\left(t\right)\end{equation} \tag{ B.12 }$$

$$\begin{equation}{H}_{I,{m}_{\mathrm{s}}=-1}^{\left(0\right)}\approx -\sum _{i}{\omega }_{i}{I}_{i,z}^{{\prime\prime}}-\sum _{i}{I}_{i,x}^{{\prime\prime}}{\gamma }_{i}{B}_{1}\left(t\right)\end{equation} \tag{ B.13 }$$

In the main text, we have dropped the '' signs for the spin operators. The first-order corrections to equations (B.12) and (B.13) are obtained by considering the terms that are first-order in θ_i in equations (B.10) and (B.11). They are given by

$$\begin{equation}{H}_{I,{m}_{\mathrm{s}}=0}^{\left(1\right)}=\sum _{i}{\gamma }_{i}{B}_{0}{\theta }_{i}\left({I}_{i,y}^{\prime }\right)=\sum _{i}{\gamma }_{i}{B}_{0}{\theta }_{i}\left({I}_{i,x}^{{\prime\prime}}\enspace \mathrm{sin}\enspace {\phi }_{i}+{I}_{i,y}^{{\prime\prime}}\enspace \mathrm{cos}\enspace {\phi }_{i}\right)\end{equation} \tag{ B.14 }$$

$$\begin{equation}{H}_{I,{m}_{\mathrm{s}}=-1}^{\left(1\right)}=0\end{equation} \tag{ B.15 }$$

In both m_s = 0 and −1 states, we have ignored the correction term due to the oscillating component of γ_iB₁(t)I'_i,z sin θ_i sin ϕ_i. While the correction term is of first order, within rotating wave approximation, the introduced errors have effects of order greater than one. Therefore, we find that the hyperfine field misalignments only influence the two-qubit gates and the first-order correction for the auxiliary subspace is effectively a static magnetic field in the x and y direction.

Including this leading order correction and the leading order correction from the secular approximation, the new two-qubit gate Hamiltonian is

$$\begin{align}\hfill {H}^{\prime }& =\frac{{\Delta}}{2}{\sigma }_{z}+\frac{{\Omega}}{2}{\sigma }_{x}{B}_{1}\left(t\right)+\sum _{i}\left({\alpha }_{i}+{\beta }_{i}{\sigma }_{z}\right){I}_{i,z}+\sum _{i}\left[\left({\gamma }_{i}{B}_{0}{\theta }_{i}\enspace \mathrm{sin}\enspace {\phi }_{i}-{\chi }_{i}+{\xi }_{i}{\theta }_{i}\enspace \mathrm{sin}\enspace {\phi }_{i}\right){I}_{i,x}\right.\hfill \\ \hfill & \quad \left.+\left({\gamma }_{i}{B}_{0}{\theta }_{i}\enspace \mathrm{cos}\enspace {\phi }_{i}-{\eta }_{i}+{\xi }_{i}{\theta }_{i}\enspace \mathrm{cos}\enspace {\phi }_{i}\right){I}_{i,y}-\left({\chi }_{i}{\theta }_{i}\enspace \mathrm{sin}\enspace {\phi }_{i}+{\eta }_{i}{\theta }_{i}\enspace \mathrm{cos}\enspace {\phi }_{i}+{\xi }_{i}\right){I}_{i,z}\right]\hfill \end{align} \tag{ B.16 }$$

where

$$\begin{equation}{\chi }_{i}=\frac{{\gamma }_{\mathrm{e}}{B}_{0}}{{D}^{2}-{\gamma }_{\mathrm{e}}^{2}{B}_{0}^{2}}\left({A}_{i,yy}{A}_{i,zx}-{A}_{i,yx}{A}_{i,zy}\right)\end{equation} \tag{ B.17 }$$

$$\begin{equation}{\eta }_{i}=\frac{{\gamma }_{\mathrm{e}}{B}_{0}}{{D}^{2}-{\gamma }_{\mathrm{e}}^{2}{B}_{0}^{2}}\left({A}_{i,xx}{A}_{i,zy}-{A}_{i,xy}{A}_{i,zx}\right)\end{equation} \tag{ B.18 }$$

$$\begin{equation}{\xi }_{i}=\frac{{\gamma }_{\mathrm{e}}{B}_{0}}{{D}^{2}-{\gamma }_{\mathrm{e}}^{2}{B}_{0}^{2}}\left({A}_{i,xy}^{2}-{A}_{i,xx}{A}_{i,yy}\right)\enspace \end{equation} \tag{ B.19 }$$

We can estimate the effect of these corrections by considering the rotations that they generate over time when the auxiliary subspace is populated during the implementation of a CZ gate. Taking into account the evolution of electron spin and nuclear spins explicitly, we adopt a simple model as follows to estimate the infidelity introduced by these correction terms. For an ideal unitary operator U_T, we consider a perfect 2π rotation on the electron spin, conditional on the nuclear spin. The unitary operator with the presence of correction terms as described equation (B.16) are given by U_A. None of these corrections are applied to the ¹⁵N nuclear spin as it is perfectly aligned with the NV quantisation axis. Thus, the additional infidelity due to these first-order corrections is given by

$$\begin{equation}\mathrm{I}\mathrm{n}\mathrm{fi}\mathrm{d}\mathrm{e}\mathrm{l}\mathrm{i}\mathrm{t}\mathrm{y}=1-\frac{\mathrm{Tr}\left[{U}_{T}^{{\dagger}}{U}_{A}\right]}{\mathrm{Tr}\left[{U}_{T}^{{\dagger}}{U}_{T}\right]}\end{equation} \tag{ B.20 }$$

where

$$\begin{equation}{U}_{T}=\mathrm{exp}\left[-\frac{\mathbb{i}}{2}\left[2\pi \left({\sigma }_{\mathrm{e},x}{\sigma }_{\mathit{n}}\right)\right]\right]\end{equation} \tag{ B.21 }$$

$$\begin{align}\hfill {U}_{A}& =\mathrm{exp}\left[-\frac{\mathbb{i}}{2}\left[2\pi \left({\sigma }_{\mathrm{e},x}{\sigma }_{\mathit{n}}\right)\right.\right.+\tau \left({\gamma }_{\mathit{n}}{B}_{0}\theta \enspace \mathrm{sin}\enspace \phi -\chi +\xi \theta \enspace \mathrm{sin}\enspace \phi \right)\left({\sigma }_{n}{\sigma }_{x}\right)\hfill \\ \hfill & \quad \left.+,\tau ,\left({\gamma }_{\mathit{n}}{B}_{0}\theta \enspace \mathrm{cos}\enspace \phi -\eta +\xi \theta \enspace \mathrm{cos}\enspace \phi \right),\left({\sigma }_{\mathit{n}}{\sigma }_{y}\right),-,\tau ,\left(\chi \theta \enspace \mathrm{sin}\enspace \phi +\eta \theta \enspace \mathrm{cos}\enspace \phi +\xi \right),\left({\sigma }_{\mathit{n}}{\sigma }_{z}\right)\right]\hfill \\ \hfill {\sigma }_{\mathit{n}}& =\left(\begin{matrix}\hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill \end{matrix}\right)\hfill \end{align} \tag{ B.22 }$$

σ_x,y,z are the Pauli matrices, σ_e and σ_n are the matrices operating on the electron spin and ¹³C nuclear spin, respectively.

For a distant nuclei, we can approximate the hyperfine terms to be A_i,xx = A_i,yy = A_i,zz and they are the leading order terms. We can also approximate A_i,zx = A_i,zy = A_i,nd where ${A}_{i,nd}^{2}={A}_{i,zx}^{2}+{A}_{i,zy}^{2}$. The ¹³C that we have chosen with A_C ≈ 0.413 MHz closely represent a ¹³C lattice site in the S family with their hyperfine parameters given by A_zz = 0.412 MHz and A_nd = 0.060 MHz [66]. The calculated hyperfine field misalignment using B₀ = 0.62 T is approximately 0.0085 rad. We also surveyed a few lattice sites which has similar hyperfine parameters obtained via ab initio method from [66]. Approximating ϕ = π/2 rad and τ = 1  μs to be the gate duration when the auxiliary subspace is populated during the implementation of a CZ gate, the additional infidelities due to hyperfine field misalignments are on the order of 10⁻² ∼ 10⁻³. These additional infidelities are dominated by the hyperfine field misalignments over the secular approximation. The correction terms from the secular approximation χ, η and ξ scale with a factor of approximately 10⁻⁵ rad s⁻¹, compared to the correction terms from the hyperfine field misalignment γ_nB₀, which has a magnitude of 10⁷ rad s⁻¹.

Given that the time evolution operator has the form of

$$\begin{equation}U=\mathrm{exp}\left(-\mathrm{i}\underset{-\tau /2}{\overset{\tau /2}{\int }}{H}_{\text{multi}}\enspace \mathrm{d}t\right)\end{equation} \tag{ C.1 }$$

Assuming there are no mixed signals in a single pulse, the Y components will go to zero. The expression for the evolution operator can be further simplified where

$$\begin{align}\hfill U& =\mathrm{exp}\left(-\mathrm{i}\underset{-\tau /2}{\overset{\tau /2}{\int }}{H}_{11}\enspace \mathrm{d}t\right)\otimes \left\vert 11\right\rangle \left\langle 11\right\vert +\mathrm{exp}\left(-\mathrm{i}\underset{-\tau /2}{\overset{\tau /2}{\int }}{H}_{10}\enspace \mathrm{d}t\right)\otimes \left\vert 10\right\rangle \left\langle 10\right\vert \hfill \\ \hfill & \quad +\mathrm{exp}\left(-\mathrm{i}\underset{-\tau /2}{\overset{\tau /2}{\int }}{H}_{01}\enspace \mathrm{d}t\right)\otimes \left\vert 01\right\rangle \left\langle 01\right\vert +\mathrm{exp}\left(-\mathrm{i}\underset{-\tau /2}{\overset{\tau /2}{\int }}{H}_{00}\enspace \mathrm{d}t\right)\otimes \left\vert 00\right\rangle \left\langle 00\right\vert \hfill \end{align} \tag{ C.2 }$$

where

$$\begin{equation}{H}_{11}=\frac{{\Omega}}{2}{B}_{1}\left(t\right)\left[{\sigma }_{x}\enspace \mathrm{cos}\left(\left[{\Delta}+{\beta }_{1}+{\beta }_{2}\right]t\right)\right]\end{equation} \tag{ C.3 }$$

$$\begin{equation}{H}_{10}=\frac{{\Omega}}{2}{B}_{1}\left(t\right)\left[{\sigma }_{x}\enspace \mathrm{cos}\left(\left[{\Delta}+{\beta }_{1}-{\beta }_{2}\right]t\right)\right]\end{equation} \tag{ C.4 }$$

$$\begin{equation}{H}_{01}=\frac{{\Omega}}{2}{B}_{1}\left(t\right)\left[{\sigma }_{x}\enspace \mathrm{cos}\left(\left[{\Delta}-{\beta }_{1}+{\beta }_{2}\right]t\right)\right]\end{equation} \tag{ C.5 }$$

$$\begin{equation}{H}_{00}=\frac{{\Omega}}{2}{B}_{1}\left(t\right)\left[{\sigma }_{x}\enspace \mathrm{cos}\left(\left[{\Delta}-{\beta }_{1}-{\beta }_{2}\right]t\right)\right]\end{equation} \tag{ C.6 }$$

To get a selective 2π pulse on the $\left\vert 11\right\rangle$ state, we need to enforce

$$\begin{equation}{X}_{1}={\Omega}\underset{-\tau /2}{\overset{\tau /2}{\int }}{B}_{1}\left(t\right)\mathrm{cos}\left[\left({\Delta}+{\beta }_{1}+{\beta }_{2}\right)t\right]\mathrm{d}t=2\pi \end{equation} \tag{ C.7 }$$

and other 3 integrals to be zero. As per the single-qubit gate case, we make the substitution

$$\begin{equation}{\Omega}{B}_{1}\left(t\right)=a\left(t\right)\mathrm{cos}\left(\lambda t\right)+b\left(t\right)\mathrm{sin}\left(\lambda t\right)\end{equation} \tag{ C.8 }$$

where λ = Δ + β₁ + β₂. Thus, we have

$$\begin{equation}{X}_{1}=\frac{\sqrt{2\pi }}{2}\left[a\left({\omega }_{11}+\lambda \right)+a\left({\omega }_{11}-\lambda \right)\right]=2\pi \end{equation} \tag{ C.9 }$$

$$\begin{equation}{X}_{2}=\frac{\sqrt{2\pi }}{2}\left[a\left({\omega }_{10}+\lambda \right)+a\left({\omega }_{10}-\lambda \right)\right]=0\end{equation} \tag{ C.10 }$$

$$\begin{equation}{X}_{3}=\frac{\sqrt{2\pi }}{2}\left[a\left({\omega }_{01}+\lambda \right)+a\left({\omega }_{01}-\lambda \right)\right]=0\end{equation} \tag{ C.11 }$$

$$\begin{equation}{X}_{4}=\frac{\sqrt{2\pi }}{2}\left[a\left({\omega }_{00}+\lambda \right)+a\left({\omega }_{00}-\lambda \right)\right]=0\end{equation} \tag{ C.12 }$$

In accordance with the derivation of infidelity of single-qubit gates, the intrinsic infidelity expression that we used to describe two-qubit gate operations is given by

$$\begin{equation}I=1-\frac{1}{8}\left[-2\enspace \mathrm{cos}\left(\frac{\sqrt{{X}_{1}^{2}}}{2}\right)+2\enspace \mathrm{cos}\left(\frac{\sqrt{{X}_{2}^{2}}}{2}\right)+2\enspace \mathrm{cos}\left(\frac{\sqrt{{X}_{3}^{2}}}{2}\right)+2\enspace \mathrm{cos}\left(\frac{\sqrt{{X}_{4}^{2}}}{2}\right)\right]\end{equation} \tag{ D.1 }$$

Keeping terms up to the second order, the final expression for the infidelity of the CZ gate is

$$\begin{equation}I=\frac{1}{32}\left(\delta {X}_{1}^{2}+{X}_{2}^{2}+{X}_{3}^{2}+{X}_{4}^{2}\right)\end{equation} \tag{ D.2 }$$

where δX₁ = X₁ − 2π. The average infidelity is given by

$$\begin{align}\hfill \left\langle I\right\rangle & =\frac{1}{32}\underset{-\infty }{\overset{\infty }{\int }}p\left(\delta ;{\sigma }_{\delta }\right)\underset{-\infty }{\overset{\infty }{\int }}p\left({\epsilon};{\sigma }_{{\epsilon}}\right)\delta {X}_{1}^{2}\enspace \mathrm{d}{\epsilon}\enspace \mathrm{d}\delta +\frac{1}{32}\underset{-\infty }{\overset{\infty }{\int }}p\left(\delta ;{\sigma }_{\delta }\right)\underset{-\infty }{\overset{\infty }{\int }}p\left({\epsilon};{\sigma }_{{\epsilon}}\right){X}_{2}^{2}\enspace \mathrm{d}{\epsilon}\enspace \mathrm{d}\delta \hfill \\ \hfill & \quad +\frac{1}{32}\underset{-\infty }{\overset{\infty }{\int }}p\left(\delta ;{\sigma }_{\delta }\right)\underset{-\infty }{\overset{\infty }{\int }}p\left({\epsilon};{\sigma }_{{\epsilon}}\right){X}_{3}^{2}\enspace \mathrm{d}{\epsilon}\enspace \mathrm{d}\delta +\frac{1}{32}\underset{-\infty }{\overset{\infty }{\int }}p\left(\delta ;{\sigma }_{\delta }\right)\underset{-\infty }{\overset{\infty }{\int }}p\left({\epsilon};{\sigma }_{{\epsilon}}\right){X}_{4}^{2}\enspace \mathrm{d}{\epsilon}\enspace \mathrm{d}\delta \hfill \end{align} \tag{ D.3 }$$

See (tables E1 and E2).

In general, there are 3 types of gate operations in quantum Fourier transform [65]. They are Hadamard gate, swap gate, and controlled-phase gate. Ignoring the global phase factors, we have

$$\begin{equation}H={R}_{x}\left(\pi \right)\cdot {R}_{y}\left(\frac{\pi }{2}\right)\end{equation} \tag{ F.1 }$$

A swap gate can be constructed from 3 CNOT gates. They can be written as

$$\begin{equation}{\mathrm{S}\mathrm{w}\mathrm{a}\mathrm{p}}_{a,b}={\mathrm{C}}_{a}{\mathrm{N}\mathrm{O}\mathrm{T}}_{b}\cdot {\mathrm{C}}_{b}{\mathrm{N}\mathrm{O}\mathrm{T}}_{a}\cdot {\mathrm{C}}_{a}{\mathrm{N}\mathrm{O}\mathrm{T}}_{b}\end{equation} \tag{ F.2 }$$

where C_aNOT_b is

$$\begin{equation}{\mathrm{C}}_{a}{\mathrm{N}\mathrm{O}\mathrm{T}}_{b}={H}_{b}\cdot {\mathrm{C}}_{a}{\mathrm{Z}}_{b}\cdot {H}_{b}\end{equation} \tag{ F.3 }$$

Here, a and b denotes any two qubits in the system. Phase gate can be written as

$$\begin{equation}\mathrm{P}\mathrm{H}\mathrm{A}\mathrm{S}\mathrm{E}\left(\theta \right)={R}_{y}\left(-\frac{\pi }{2}\right)\cdot {R}_{x}\left(\frac{\theta }{2}\right)\cdot {R}_{y}\left(\frac{\pi }{2}\right)\end{equation} \tag{ F.4 }$$

For a controlled-phase gate, one way to construct them is as follows

$$\begin{align}\hfill {\mathrm{C}}_{a}{\mathrm{P}\mathrm{H}\mathrm{A}\mathrm{S}\mathrm{E}}_{b}\left(\theta \right)& ={R}_{y,a}\left(-\frac{\pi }{2}\right){R}_{x,a}\left(\frac{\theta }{2}\right){R}_{y,a}\left(\frac{\pi }{2}\right){R}_{x,b}\left(-\frac{\pi }{2}\right){R}_{y,b}\left(\frac{-\theta -\pi }{4}\right){\mathrm{C}}_{a}{\mathrm{Z}}_{b}{R}_{y,b}\left(\frac{\theta +\pi }{4}\right)\hfill \\ \hfill & \quad {\times}{R}_{x,b}\left(\frac{\pi }{2}\right){R}_{y,b}\left(-\frac{\pi }{2}\right){R}_{x,b}\left(\frac{-\theta -\pi }{4}\right){\mathrm{C}}_{a}{\mathrm{Z}}_{b}{R}_{x,b}\left(\frac{\theta +\pi }{4}\right){R}_{y,b}\left(\frac{\pi }{2}\right)\hfill \end{align} \tag{ F.5 }$$

where θ is the intended phase. Note that the operations for these gate decompositions are performed from the right to the left. Using these decompositions, the total number of pulses required for QFT3 and QFT5 is 75 and 195 respectively.