Self-assembling hybrid diamond–biological quantum devices

Dipolar interactions between electron spins of two adjacent NV-centers provide a possibility for implementing gate operations [5, 6]. Combined with the continuous driving of the decoupling field, the total effective Hamiltonian in a frame rotating with the microwave frequency of the driving field can be written as (see appendix C)

$$\begin{eqnarray}&&H^{\prime} \simeq \mathop{\sum }\limits_{i}\frac{{{b}_{i}}(t)}{2}\;\sigma _{z}^{i}+\mathop{\sum }\limits_{i}\frac{{{\Omega }_{i}}}{2}\;\sigma _{x}^{i}+\mathop{\sum }\limits_{i\gt j}\frac{{{J}_{ij}}}{2}\;\sigma _{z}^{i}\;\sigma _{z}^{j}\;.\end{eqnarray} \tag{ 6 }$$

Herein the first two parts describe the dephasing noise and decoupling for the individual NV-centers as introduced in the previous section, respectively. The last part accounts for the dipolar coupling with ${{J}_{ij}}=2\;{{\xi }_{ij}}{{\mu }_{ij}}$, where ξ_ij depends on the external magnetic field strength $|\vec{B}|$ and its orientation with respect to both the NV symmetry axes and the vector ${{\vec{r}}_{ij}}$ connecting NV-center i and j. In the limit of high magnetic fields ${{\xi }_{ij}}=1/4(1-3\;{{{\rm cos} }^{2}}{{\theta }_{ij}})$, where ${{\theta }_{ij}}=\angle ({{\vec{r}}_{ij}},\vec{B})$ and ${{\mu }_{ij}}={{\mu }_{0}}\gamma _{{\rm el}}^{2}\hbar /(4\pi r_{ij}^{3})$ with the latter being ${{\mu }_{ij}}=52\;{\rm kHz}$ for an NV-center distance ${{r}_{ij}}=10\;{\rm nm}$. μ₀ and γ_el denote the magnetic permeability and the electron gyromagnetic ratio, respectively. It is important to note that the configuration of NDs leads to a random relative orientation of the NV symmetry axes in space, the latter forming a 'natural' quantization axis along the N–V direction by the associated crystal-field splitting of the ground state triplet ($D\sim 2.87\;{\rm GHz}$) (see figure 1(a)). This leads to two important consequences: first, a uniform quantization axis has to be defined by a sufficiently strong external magnetic field to guarantee a uniform dipolar coupling, as illustrated in figure 1(c) together with its optimal 2D-orientation. As shown in the figure, this can be achieved applying a magnetic field $B\;\gtrsim \;0.5\;{\rm T}$ (${{\gamma }_{{\rm e}}}\;B/(2\;\pi )\;\gtrsim \;14\;{\rm GHz}$), in which case the magnetic energy shift dominates over the crystal-field splitting (${{\gamma }_{{\rm e}}}\;B\gg D$). Second, the crystal-field splitting is responsible for an orientation dependent transition frequency in that regime, depending on the angle between the symmetry axis and the external field ϑ_i and scaling as $\propto \;D\;{\rm cos} 2{{\vartheta }_{i}}$. As a consequence, the transition frequencies of individual NV-centers differ by typical values of several $100\;{\rm MHz}$, thus providing the possibility for individual microwave addressing for the values of Ω obtained from figure 1 and at the same time maintaining a uniform dipolar coupling interaction. Moreover these inhomogeneous energy splittings justify the omission of dipolar flip-flop interaction terms in (6), that would not be energy conserving in our setup. As a remark, spin-mixing in the ground state triplet by magnetic fields not aligned with the symmetry axis, reduces both the contrast of optical spin readout and the spin initialization efficiency by optical cycling [41]. This might be overcome by projective readout techniques [42] and spin initialization at low magnetic fields followed by an adiabatic increase of the field amplitude.

Combining the dipolar interaction with decoupling suppresses the environmental coupling, i.e. decoherence, but also part of the gate interaction, a general problem in the application of decoupling sequences. We will consider that effect for two distinct decoupling configurations: ${\mathscr{M}_{1}}$, corresponding to a homogeneous decoupling field ${{\Omega }_{i}}=\Omega$, and ${\mathscr{M}_{2}}$ corresponding to ${{\Omega }_{i}}=-{{\Omega }_{j}}$ for $i\in {\rm neighb}(j)$. These will turn out to provide crucial gate interactions and in addition their combination allows one to restore the original form of the dipolar coupling as will be discussed in section 4.2. The effective interaction in the presence of decoupling is best seen by transforming the original Hamiltonian

$$\begin{eqnarray}&&{{H}_{{\mathscr{M}_{k}}}}=H_{{\mathscr{M}_{k}}}^{0}+{{H}_{{\rm dip}}},\end{eqnarray} \tag{ 7 }$$

wherein ${{H}_{{\rm dip}}}$ denotes the dipolar part and $H_{{\mathscr{M}_{k}}}^{0}$ the continuously driven system in the decoupling configuration ${\mathscr{M}_{k}}$, to an interaction picture with respect to the driving, which in the relevant limit $\Omega \gg {{J}_{ij}}$ results in

$$\begin{eqnarray}&&{{H}_{I,{\mathscr{M}_{k}}}}={{{\rm e}}^{+{\rm i}\,H_{{\mathscr{M}_{k}}}^{0}\,t}}\;{{H}_{{\rm dip}}}\;{{{\rm e}}^{-{\rm i}\,H_{{\mathscr{M}_{k}}}^{0}\,t}}\simeq \frac{1}{2}\;\mathop{\sum }\limits_{(i,j)}{{J}_{ij}}\;S_{{\mathscr{M}_{k}}}^{ij}\;.\end{eqnarray} \tag{ 8 }$$

Here the Hamiltonian has been restricted to nearest neighbor interactions, i.e. (i, j) sums over neighboring NV-centers, and the effective interactions take the form

$$\begin{eqnarray}&&S_{{\mathscr{M}_{1}}}^{ij}=s_{+}^{i}\;s_{-}^{j}+{\rm h}.{\rm c}.\end{eqnarray} \tag{ 9a }$$

$$\begin{eqnarray}&&S_{{\mathscr{M}_{2}}}^{ij}=s_{+}^{i}\;s_{+}^{j}+{\rm h}.{\rm c}.\;,\end{eqnarray} \tag{ 9b }$$

taking into account that the off-resonant contributions are suppressed as the decoupling field and dipolar coupling are well-separated in frequency. ${{s}_{+}}=|+\rangle \langle -|$ and ${{s}_{-}}=|-\rangle \langle +|$ are the ladder operators in the σ_x -eigenbasis $|\pm \rangle =1/\sqrt{2}(\pm \;\left| 1 \right\rangle +\left| 0 \right\rangle )$, related to the ladder operators ${{\sigma }_{\pm }}$ in the standard σ_z -eigenbasis by ${{s}_{\pm }}=U{{\sigma }_{\pm }}{{U}^{\dagger }}$ with $U={\rm exp} (-{\rm i}\pi /4{{\sigma }_{y}})$. These basic interactions are illustrated in the inset of figure 3(a) for a two qubit interaction.

Zoom In Zoom Out Reset image size

The two qubit setting is shown in figures 3(a–b), illustrating the fidelity and time evolution to create a maximally entangled state by a $\pi /2$ and $5\pi /2$ two qubit rotation in the ${\mathscr{M}_{1}}$- coupling manifold, respectively. Note that the required time for the latter case exceeds T₂(Ω = 0) by more than a factor of three, therefore impressively demonstrating the decoupling that allows for a 95% final state fidelity with a decoupling field of $\Omega =1.5\;{\rm MHz}$.

In contrast to these simple two qubit examples, recovering the full dipolar Ising-type interaction form, that is the ${{\sigma }_{z}}\otimes {{\sigma }_{z}}$-type coupling appearing in (6), marks a crucial step towards the realization of a universal set of quantum gates. Its importance is reflected in a wide range of applications such as the creation of cluster states for quantum computation or in quantum simulations, both discussed in the upcoming sections 4.3 and 4.4.

The Ising-type interaction, originally altered and partially suppressed by the decoupling field, cannot be recovered by any local operation out of the reduced decoupled manifold interactions ${{H}_{I,{\mathscr{M}_{k}}}}$. However the decoupled dipolar interaction contributions ${{H}_{I,{\mathscr{M}_{1}}}}$ and ${{H}_{I,{\mathscr{M}_{2}}}}$, corresponding to the two different decoupling field settings, can be added in time in order to recover the full dipolar interaction form based on the observation

$$\begin{eqnarray}&&{{H}_{zz}}={{H}_{I,{\mathscr{M}_{1}}}}+{{H}_{I,{\mathscr{M}_{2}}}}=\mathop{\sum }\limits_{(i,j)}({{J}_{i,j}}/2)\sigma _{z}^{i}\sigma _{z}^{j}\;.\end{eqnarray} \tag{ 10 }$$

In the Trotter (11a ) or Suzuki–Trotter (11b ) formalism this can be performed by

$$\begin{eqnarray}&&{\rm exp} \left( -{\rm i}\;{{H}_{zz}}\;t \right)={{\left[ {{R}_{{\mathscr{M}_{1}}}}\left( t/n \right){{R}_{{\mathscr{M}_{2}}}}\left( t/n \right) \right]}^{n}}+\mathscr{O}({{({{J}_{ij}}t/n)}^{2}})\end{eqnarray} \tag{ 11a }$$

$$\begin{eqnarray}&&{\rm exp} \left( -{\rm i}\;{{H}_{zz}}\;t \right)={{\left[ {{R}_{{\mathscr{M}_{1}}}}\left( t/(2n) \right){{R}_{{\mathscr{M}_{2}}}}\left( t/n \right){{R}_{{\mathscr{M}_{1}}}}\left( t/(2n) \right) \right]}^{n}}+\mathscr{O}({{({{J}_{ij}}t/n)}^{3}})\end{eqnarray} \tag{ 11b }$$

with ${{R}_{{\mathscr{M}_{k}}}}(t)={\rm exp} (-{\rm i}\;{{H}_{I,{\mathscr{M}_{k}}}}\;t)$. The switching frequency between the two types of interactions ${\mathscr{M}_{1}}$ and ${\mathscr{M}_{2}}$, or equivalently the maximal magnitudes of the addition sequence time intervals $t/n$ that still provide sufficiently low errors, depend crucially on the Hamiltonian timescales. For the addition of the effective (interaction frame) Hamiltonians ${{H}_{I,{\mathscr{M}_{k}}}}$ as following out of (10), this timescale is determined by the rather small (kHz) dipolar coupling frequency. In contrast, a time addition directly based on (7) as could be obtained by the two purely global π-phase shifted decoupling configurations ${\mathscr{M}_{1}}$ and ${\mathscr{M}_{3}}$, the latter defined by ${{\Omega }_{i}}\equiv -\Omega$, would have to be based on the much higher (MHz) magnitudes of the decoupling field Rabi frequency $|\Omega |$. Even though this would recover the pure dipolar coupling, at the same time the decoupling field sign inversion on a timescale $\tau \;\lesssim \;1/\Omega$ would effectively average out the decoupling effect. Therefore this is not a viable approach and we will focus on the implementation (10) instead. As the effective Hamiltonians ${{H}_{{\mathscr{M}_{1}}}}$ and ${{H}_{{\mathscr{M}_{2}}}}$ are the result of distinct, non-commuting interaction frames ($[H_{{\mathscr{M}_{1}}}^{0},H_{{\mathscr{M}_{2}}}^{0}]\ne 0$), each pulse has to be carefully adjusted for their implementation and will be discussed in the following.

Recalling the interaction picture definition, its time evolution can be created by the evolution sequence

$$\begin{eqnarray}&&{{R}_{{\mathscr{M}_{k}}}}(t)\equiv {{{\rm e}}^{-{\rm i}\,{{H}_{I,{\mathscr{M}_{k}}}}t}}={{{\rm e}}^{+{\rm i}\,H_{{\mathscr{M}_{k}}}^{0}\,t}}\;{{{\rm e}}^{-{\rm i}\,{{H}_{{\mathscr{M}_{k}}}}\,t}}={{E}_{{\mathscr{M}_{k}}}}(t){{U}_{{\mathscr{M}_{k}}}}(t)\end{eqnarray} \tag{ 12 }$$

with ${{U}_{{\mathscr{M}_{k}}}}(t)={\rm exp} (-{\rm i}\;{{H}_{{\mathscr{M}_{k}}}}\;t)$ and ${{E}_{{\mathscr{M}_{k}}}}(t)={\rm exp} ({\rm i}\;H_{{\mathscr{M}_{k}}}^{0}\;t)$. Three options can be exploited in order to create (12).

• (1)  
  Directly implementing the pulse sequence by applying fast local microwave pulses in order to implement ${{E}_{{\mathscr{M}_{k}}}}(t)$. Such local manipulations are feasible as they can be implemented on much faster timescales as opposed to the rather slow (kHz) dipolar coupling dynamics.
• (2)  
  Adjusting the total time by making use of the vastly different timescales of the decoupling field Ω and the dipolar interaction J_ij both appearing in the Hamiltonian ${{H}_{{\mathscr{M}_{k}}}}$. That is, on short timescales $\Delta t\ll J_{ij}^{-1}$, the dipolar part characterized by kHz frequencies can be neglected compared to the MHz evolution of the decoupling field, such that (12) can be implemented to a good approximation by a simple time adjustment. This is based on two observations, namely that ${\rm exp} ({\rm i}\;\Omega /2{{\sigma }_{x}}\;t)={\rm exp} (-{\rm i}\;\Omega /2\;{{\sigma }_{x}}t^{\prime} )$ with the 'short' time $t^{\prime} =(2\;\pi -{\rm mod} \left[ \Omega \;t,2\;\pi \right])/\Omega$ determined by the fast decoupling field frequency scale,and that on such a timescale ${{E}_{{\mathscr{M}_{k}}}}(t^{\prime} )\simeq {{U}_{{\mathscr{M}_{k}}}}(t^{\prime} )$. Therefore (12) can be implemented by the time adjustment ${{R}_{{\mathscr{M}_{k}}}}(t)\simeq {{U}_{{\mathscr{M}_{k}}}}(t+t^{\prime} )$, applicable as long as the decoupling field is sufficiently strong $\Omega \gg {{J}_{ij}}$.
• (3)  
  Implementing echo pulses in order to remove the ${{E}_{{\mathscr{M}_{k}}}}$ contribution, which changes the interaction picture evolution from the normal one according to (12). The application of local π-pulses ${{S}_{\pi }}={\rm exp} (-{\rm i}(\pi /2){{\sum }_{i}}\sigma _{i}^{k})$ either in the z or y-direction (k = z, y), allows for the implementation of (12) by

  $$\begin{eqnarray}&&{{R}_{{\mathscr{M}_{k}}}}(t)=S_{\pi }^{\dagger }\;{{U}_{{\mathscr{M}_{k}}}}(t/2){{S}_{\pi }}\;{{U}_{{\mathscr{M}_{k}}}}(t/2).\end{eqnarray} \tag{ 13 }$$

This follows out of the invariance of ${{H}_{I,{\mathscr{M}_{k}}}}$ under such an echo pulse sequence, $S_{\pi }^{\dagger }H_{{\mathscr{M}_{k}}}^{0}{{S}_{\pi }}=-H_{{\mathscr{M}_{k}}}^{0}$ and the commutation of both parts $[{{H}_{I,{\mathscr{M}_{k}}}},H_{{\mathscr{M}_{k}}}^{0}]=0$. Originating in the non-commutativity of ${{H}_{I,{\mathscr{M}_{k}}}}$ and $H_{{\mathscr{M}_{k^{\prime} }}}^{0}$ for $k\ne k^{\prime}$ this echo pulse sequence has to be applied to each ${\mathscr{M}_{k}}$ interaction separately and cannot be implemented globally on the total evolution sequence.

Based on the two qubit interaction (10) and the ability for local addressing and manipulations as a result of quantization axis dependent magnetic field level shifts, all other types of one and two qubit interactions can be implemented. This again benefits from the different timescales between the two-qubit dipolar interaction and the local microwave addressing. Therefore, local operations can be treated as instantaneous pulses with respect to the two qubit evolution timescale, which besides local operations allows one to change the interaction type to e.g. $\sigma _{i}^{x}\;\sigma _{j}^{x}={{U}_{y}}\;\sigma _{i}^{z}\;\sigma _{j}^{z}\;U_{y}^{\dagger }$ with ${{U}_{y}}={\rm exp} \left( {\rm i}\;\pi /4(\sigma _{y}^{i}+\sigma _{y}^{j}) \right)$. More complicated Hamiltonians can be implemented based on time additions (11a , 11b ) as will be illustrated for a Heisenberg-chain Hamiltonian in section 4.4.

It is also worth noting that we have limited the previous discussion to next neighbor interactions in the dipolar coupling, and only in that case the form of Hamiltonians (10) and (8) is strictly valid. For short rotation angles and linear arrangements this forms a good approximation as can be seen in the upcoming simulations for the cluster state (figure 4(a)) and Heisenberg simulations (figure 5). However for higher dimensional arrangements and long time evolutions, higher-order couplings have to be removed by combining the time addition sequences with local pulses, which will be discussed in more detail in appendix E.3 and is performed in the 2D cluster state simulation depicted in figure 4(b).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

An interesting application of the concepts developed in the preceding sections is the creation of cluster states. These highly entangled states, defined as the unique eigenstates of the multi-body generators ${{K}^{(i)}}=\sigma _{x}^{i}{{\otimes }_{(j,i)}}\sigma _{z}^{j}$ via the eigenvalue equation ${{K}^{(i)}}\;\left| {{\phi }_{\mathscr{C}}} \right\rangle =\left| {{\phi }_\mathscr{C}} \right\rangle \;\forall i$, allow one to perform any quantum computation operation by purely local measurements on individual qubits [43]. It has been shown [43] that the product of two-body phase gates S applied to a specific initial product state allows for the creation of such an eigenstate, namely ${{\left| \phi \right\rangle }_{\mathscr{C}}}=S\left| ++\cdots + \right\rangle$. The connection to the Ising Hamiltonian (10) that can be realized in the ND system proposed follows if we note that $S\simeq {\rm exp} \left( -{\rm i}\;\pi /4\;{{\sum }_{(i,j)}}\sigma _{z}^{i}\;\sigma _{z}^{j} \right)$ up to local operations, the latter being implementable by fast microwave pulses on the electron spin manifold. Thus, the availability of the interaction creating the cluster resource state, combined with local addressing and the well-separated timescales of the microwave and dipolar interaction, makes the assembled ND system a promising candidate for cluster state computation implementations. Numerical simulations of the creation of one- and two-dimensional cluster states for different numbers of qubits and decoupling field strengths are shown in figure 4. For the one-dimensional case, qubit number dependent fidelities between 80% and 90% and even above can be achieved, close to the noise-free limit imposed by non-next nearest neighbor couplings and errors. As expected, the fidelity decreases with the number of qubits involved, which both makes the state overlap to the perfect reference state more prone to deviations and increases the impact of errors. These fidelities could be further improved by removing non-nearest neighbor interactions upon integration of local pulses into the addition sequence as outlined in appendix E.3 and by decreasing the timesteps in the addition sequence, i.e. increasing n in (11b ). For the numerical simulations of the linear arrangement, a modest complexity of n = 2 has been chosen. Note however, that the timesteps (t/n) in the addition sequence have to be sufficiently large, that is $\Omega (t/n)\gg 1$, in order for the effective Hamiltonian form (8) to be valid, i.e. significantly increasing n requires one to increase the Rabi frequency Ω as well.

Figure 4(b) illustrates the fidelity for the creation of a 2D four qubit cluster state. As in that case, cross-diagonal non-nearest neighbor couplings are significant ($\sim 0.35$ of the nearest neighbor interaction), so they have been removed by a more involved addition sequence as illustrated in the lower part of figure 4(b). Herein the lowest level (i) assumes the creation of (10) as described in section 4.2. Including the off diagonal contributions this results in

$$\begin{eqnarray}&&H_{zz}^{2{\rm d}}=\frac{J}{2}\left[ \mathop{\sum }\limits_{(i,j)}\;\sigma _{z}^{i}\sigma _{z}^{j}+{{\left( \frac{1}{\sqrt{2}} \right)}^{3}}\mathop{\sum }\limits_{\lt i,j\gt }\;2\;S_{{\mathscr{M}_{1}}}^{i,j} \right]\end{eqnarray} \tag{ 14 }$$

with ${{S}_{{\mathscr{M}_{1}}}}$ as defined in (9a) and (i, j) summing over nearest neighbors whereas $\lt i,j\gt$ accounts for off-diagonal couplings. By using that ${{U}_{i}}\;S_{{\mathscr{M}_{1}}}^{ij}\;U_{i}^{\dagger }=-S_{{\mathscr{M}_{1}}}^{ij}$ and ${{U}_{i}}\;\sigma _{z}^{i}\sigma _{z}^{j}\;U_{i}^{\dagger }=-\sigma _{z}^{i}\sigma _{z}^{j}$ for ${{U}_{i}}={\rm exp} (-{\rm i}\pi /2\sigma _{x}^{i})$, a local application of U_i pulses on selected (red marked) qubits allows one to remove cross diagonal couplings by the time-addition step (i) in figure 4(b). As this also removes part of the desired nearest neighbor Ising-type coupling, a second time addition (ii) serves to fix that issue. As a remark, such a combination of the addition sequence with local unitary operations can also serve to adjust and weight individual coupling components J_ij . The fidelity for the cluster state creation approaches 80% for such a sequence and a comparison to the noise-free case reveals that this is clearly limited by the coherence time. Note that the low fidelity for $\Omega \;\lesssim \;0.2\;{\rm MHz}$ in the noise-free (dashed) simulation arises from the effect that the effective Hamiltonian description (8) is only valid for sufficiently large Rabi frequencies, whereas the fidelity saturation value for $\Omega \;\gtrsim \;0.5\;{\rm MHz}$ reveals the error that can be ascribed to the time-addition implementation.

As a second application of the concepts developed in section 4.2, we have illustrated the possibility of simulating a Heisenberg-chain Hamiltonian of the form $H=-1/2\;\sum _{j=1,\mu }^{N}{{J}_{\mu }}\;\sigma _{\mu }^{j}\;\sigma _{\mu }^{j+1}$ with $\mu =\{x,y,z\}$ and ${{J}_{y}}={{J}_{z}}=J$, ${{J}_{x}}=\delta \;J$ in figure 5. This task essentially requires three steps as depicted in figure 5(b): creating the $\sigma _{x}^{j}\;\sigma _{x}^{j+1}$-contribution out of the basic decoupled ${\mathscr{M}_{1}}$ and ${\mathscr{M}_{2}}$ interaction following the concept outlined in section 4.2, and mapping σ_z − σ_x by the unitary local pulse operation ${{U}_{y}}={\rm exp} \left( -{\rm i}{{\sum }_{j}}\;\pi /4\;\sigma _{y}^{j} \right)$. Second, creating the $S_{{\mathscr{M}_{1}}}^{j,j+1}$ contribution, which directly corresponds to an effective decoupled interaction (8) and is therefore obtained in a trivial way. In a last step both contributions are added based on the Trotterization approach (11b ) taking the weighting factor δ into account. The accuracy of such a procedure is illustrated in figure 5(a) for a $\theta =\pi /2$ evolution. Significant improvements can be obtained by increasing the decoupling Rabi frequency. A comparison to noise-free simulations, revealing the effect of Trotterization and non-nearest neighbor errors, shows that the final fidelity is decoherence limited.

Distance variations between adjacent vacancy centers and different orientations of the symmetry axes make it hard in practice to guarantee a uniform coupling. Therefore the coupling coefficient J_ij appearing in ${{H}_{I,{\mathscr{M}_{k}}}}$ (8) and H_zz (10) may be replaced by $J(1+{{\epsilon }_{ij}})$ with _ij describing the systematic error from the optimal case. However, extending the concepts provided in [44] allows one to construct compensation sequences as will be discussed in detail in appendix E, provided that the error ${{\epsilon }_{ij}}\;\lesssim \;0.5$ and that non-nearest neighbor couplings can be efficiently suppressed, as discussed in appendix E.3. We analyzed the compensation method for the two qubit case in figure 3(c) and applied it to a four qubit cluster state in figure E1. Due to the significantly increased time, that is eight and sixteen times the original gate operation, respectively, the region of benefit increases with the decoupling field strength provided that it exceeds a threshold magnitude. As expected the two qubit gate compensation is more efficient providing good results already for a $1.5\;{\rm MHz}$ decoupling field in contrast to the multiqubit counterpart that relies on a more general sequence less efficient in time.

Consider a two level system evolving under the Hamiltonian (in the rotating frame)

$$\begin{eqnarray}&&H=\hbar \frac{b(t)}{2}\;{{\sigma }_{z}}+\hbar \;\frac{\Omega (t)}{2}\;{{\sigma }_{x}}\end{eqnarray} \tag{ B.1 }$$

with b(t) a random, zero-mean fluctuating detuning describing the effect of pure dephasing and originating from the environmental coupling and $\Omega (t)$ the classical control decoupling field. We will assume that the system is initially (${{t}_{0}}=0)$ prepared in the state

$$\begin{eqnarray}&&|{{\psi }_{\phi }}\rangle =\frac{1}{\sqrt{2}}\left( |e\rangle +{{{\rm e}}^{{\rm i}\,\phi }}\;|g\rangle \right)\;,\end{eqnarray} \tag{ B.2 }$$

then, after a time t the probability for still finding the system in that initial state is given by

$$\begin{eqnarray}&&\langle {{\psi }_{\phi }}|\rho |{{\psi }_{\phi }}\rangle =\frac{1}{2}\left( 1+{{{\rm cos} }^{2}}\phi \;{{{\rm e}}^{-{{R}_{x}}(t)t}}+{{{\rm sin} }^{2}}\phi \;{{{\rm e}}^{-1/2\left( {{R}_{x}}(t)+{{R}_{\gamma }}(t) \right)t}} \right)\;,\end{eqnarray} \tag{ B.3 }$$

describing purely the effect of decoherence and not taking the coherent evolution into account (see figure B2 ). The decay rates appearing in appendix B.3 can be expressed as

$$\begin{eqnarray}&&{{R}_{k}}(t)=\frac{1}{2t}\frac{1}{2\pi }\;\int _{-\infty }^{\infty }\;S(\omega )F_{t}^{k}(\omega ){\rm d}\omega \end{eqnarray} \tag{ B.4 }$$

with $S(\omega )=\int _{-\infty }^{\infty }{\rm exp} (-{\rm i}\omega \tau )\langle b(t)b(t+\tau )\rangle {\rm d}\tau$ the noise spectrum and $F_{t}^{k}(\omega )$ a decoupling field dependent filter function (see figure B1). Exact expressions can be found in appendix B.3.

Zoom In Zoom Out Reset image size

In the limit of $t\gg \tau$ and $\Omega \;t\gt 1$ and for a constant Ω, ${{R}_{x}}(t)\simeq 1/2\;S(\Omega )$ and ${{R}_{\gamma }}(t)\simeq 0$ (more precise ${{R}_{\gamma }}(t)t\ll {{R}_{x}}(t)t$), the coherence decay in time (B.3) takes the form

$$\begin{eqnarray}&&\langle {{\psi }_{\phi }}|\rho |{{\psi }_{\phi }}\rangle =\frac{1}{2}\left( 1+{{{\rm cos} }^{2}}\phi \;{{{\rm e}}^{-{{R}_{x}}\,t}}+{{{\rm sin} }^{2}}\phi \;{{{\rm e}}^{-1/2\,{{R}_{x}}\,t}} \right)\quad \quad \quad {\rm with}\ \quad {{R}_{x}}=\frac{1}{2}\;S(\Omega ).\end{eqnarray} \tag{ B.5 }$$

This has a clear interpretation in that the first part (the σ_x eigenstate decay) is related to a population decay whereas the second part (the σ_y eigenstate decay) corresponds to a decay of the corresponding coherences, a process that happens in the Markovian limit $t\gg \tau$—also characterized by a time independent decay rate—with half of the population decay rate. In the non-Markovian limit the decay rate of the coherence contributions is increased by an additional factor ${{R}_{\gamma }}(t)$. Note also that for the case of free induction decay, i.e. $\Omega =0$, ${{R}_{\gamma }}(t)={{R}_{x}}(t)$ and thus both decay contributions are equal as expected by the isotropy of the system.

Following the description in [39, 55] we will derive the general decoherence decay formula (B.3). It is advantageous to evaluate the decoherence behavior in the σ_x eigenbasis $|\pm {{\rangle }_{x}}\equiv |\pm \rangle$ in which the effect of the decoupling field can be interpreted as creating an energy gap suppressing flipping processes by the (off-resonant) detuning fluctuations. Using the Nakajima–Zwanzig projection operator approach [57] allows one to obtain a master equation for the system interacting with a noise bath up to second order in the coupling constant for the evolution under (B.1)[39]

$$\begin{eqnarray}&&\frac{{\rm d}\rho }{{\rm d}t}=-\frac{1}{{{\hbar }^{2}}}\;\int _{0}^{t}{\rm d}t^{\prime} \;\phi (t-t^{\prime} )\left[ S(t),S(t^{\prime} )\rho (t) \right]+{\rm h}.{\rm c}.\end{eqnarray} \tag{ B.6 }$$

with $\phi (t-t^{\prime} )=\langle b(t)b(t^{\prime} )\rangle$, $S(t)=\hbar /2\;{{\tilde{\sigma }}_{z}}$ and ${{\tilde{\sigma }}_{z}}={\rm exp} \left( {\rm i}\int _{0}^{t}\Omega (t^{\prime} )/2\;{\rm d}t^{\prime} \;{{\sigma }_{x}} \right){{\sigma }_{z}}\;{\rm exp} \left( -{\rm i}\int _{0}^{t}\Omega (t^{\prime} )/2\;{\rm d}t^{\prime} \;{{\sigma }_{x}} \right)$ the σ_z operator in the interaction picture with respect to the decoupling field. It turns out that coherences and populations are decoupled in the differential equation expressed in the σ_x -basis states leading to (${{\rho }_{ij}}=\langle i|\rho |j\rangle$),

$$\begin{eqnarray}&&\frac{{\rm d}}{{\rm d}t}\left( \begin{array}{ccccccccccccccc} {{\rho }_{++}} \\ {{\rho }_{--}} \\ \end{array} \right)=-\frac{1}{2}\;{{\gamma }_{1}}(t)\left( \begin{array}{ccccccccccccccc} 1 & -1 \\ -1 & 1 \\ \end{array} \right)\left( \begin{array}{ccccccccccccccc} {{\rho }_{++}} \\ {{\rho }_{--}} \\ \end{array} \right),\end{eqnarray} \tag{ B.7 }$$

wherein (see definitions in section appendix B.3, $U={\rm exp} ({\rm i}\int _{0}^{t}\Omega (\tau ){\rm d}\tau )$, $\mathscr{R}$ denotes the real part)

$$\begin{eqnarray}&&\qquad {{\gamma }_{1}}(t)=\int _{0}^{t}\;{\rm d}t^{\prime} \;\phi (t-t^{\prime} )\mathscr{R}\left[ U(t){{U}^{\dagger }}(t^{\prime} ) \right],\quad \quad \quad {{R}_{x}}(t)=\frac{1}{t}\;\int _{0}^{t}{\rm d}t^{\prime} \;{{\gamma }_{1}}(t^{\prime} ){\rm d}t^{\prime} \end{eqnarray} \tag{ B.8 }$$

leading to the solutions

$$\begin{eqnarray}&&\quad {{\rho }_{++}}(t)=\frac{1}{2}\left[ (2\;\rho _{++}^{0}-1){{{\rm e}}^{-{{R}_{x}}(t)t}}+1 \right]\;,\quad \quad \quad {{\rho }_{--}}(t)=\frac{1}{2}\left[ (2\;\rho _{--}^{0}-1){{{\rm e}}^{-{{R}_{x}}(t)t}}+1 \right]\;.\end{eqnarray} \tag{ B.9 }$$

On the other hand the differential equation for the coherences takes the form

$$\begin{eqnarray}&&\frac{{\rm d}}{{\rm d}t}\left( \begin{array}{ccccccccccccccc} {{\rho }_{-+}}-{{\rho }_{+-}} \\ {{\rho }_{-+}}+{{\rho }_{+-}} \\ \end{array} \right)=-\frac{1}{2}\left( \begin{array}{ccccccccccccccc} {{\gamma }_{1}}(t)+{{\gamma }_{2}}(t) & {\rm i}({{\mu }_{1}}(t)-{{\mu }_{2}}(t)) \\ {\rm i}({{\mu }_{1}}(t)+{{\mu }_{2}}(t)) & {{\gamma }_{1}}(t)-{{\gamma }_{2}}(t) \\ \end{array} \right)\left( \begin{array}{ccccccccccccccc} {{\rho }_{-+}}-{{\rho }_{+-}} \\ {{\rho }_{-+}}+{{\rho }_{+-}} \\ \end{array} \right)\end{eqnarray} \tag{ B.10 }$$

with the additional definitions (herein $\mathscr{R}$ denotes the real and $\mathscr{I}$ the imaginary part)

$$\begin{eqnarray}\begin{array}{llllll} {{\gamma }_{2}}(t) & = & \int _{0}^{t}\;{\rm d}t^{\prime} \;\phi (t-t^{\prime} )\mathscr{R}\left[ U(t)U(t^{\prime} ) \right], & \quad \quad \quad {{R}_{\gamma }}(t) & = & \frac{1}{t}\;\int _{0}^{t}{\rm d}t^{\prime} \;{{\gamma }_{2}}(t^{\prime} ){\rm d}t^{\prime} \\ {{\mu }_{1}}(t) & = & \int _{0}^{t}\;{\rm d}t^{\prime} \;\phi (t-t^{\prime} )\mathscr{I}\left[ U(t){{U}^{\dagger }}(t^{\prime} ) \right], & \quad \quad {{\mu }_{2}}(t) & = & \int _{0}^{t}\;{\rm d}t^{\prime} \;\phi (t-t^{\prime} )\mathscr{I}[U(t)U(t^{\prime} )]. \\ \end{array}\end{eqnarray} \tag{ B.11 }$$

As will be seen later, the combination ${{\rho }_{-+}}-{{\rho }_{+-}}$ is responsible for the decay of coherences in the density matrix description. Moreover it is straightforward to show that the off-diagonal elements (the μ-terms) are related to a coherent evolution whereas the diagonal elements describe the decay terms of the corresponding quantities. Since we are not interested in the coherent evolution it is possible to set those off-diagonal contributions to zero resulting in a description of the envelope decay of the quantities. Therefore one ends up with

$$\begin{eqnarray}&&({{\rho }_{-+}}-{{\rho }_{+-}})(t)={{{\rm e}}^{-1/2({{R}_{x}}(t)+{{R}_{\gamma }}(t))t}}({{\rho }_{-+}}-{{\rho }_{+-}})(0).\end{eqnarray} \tag{ B.12 }$$

Now consider the arbitrary phase state (B.2) that can be expressed in the $|\pm \rangle$ basis as (neglecting a global phase factor)

$$\begin{eqnarray}&&|{{\psi }_{\phi }}\rangle ={\rm cos} (\phi /2)|+\rangle -i\;{\rm sin} (\phi /2)|-\rangle ,\quad {{\rho }_{0}}=|{{\psi }_{\phi }}\rangle \langle {{\psi }_{\phi }}|\;.\end{eqnarray} \tag{ B.13 }$$

Using (B.9) and (B.12), the density matrix after a time evolution of t follows to be

$$\begin{eqnarray}\begin{array}{rcl} \rho (t) & = & |+\rangle \langle +|\;\frac{1}{2}\;\left[ {\rm cos} \phi \;{{{\rm e}}^{-{{R}_{x}}(t)t}}+1 \right]+|-\rangle \langle -|\;\frac{1}{2}\left[ -{\rm cos} \phi \;{{{\rm e}}^{-{{R}_{x}}(t)t}}+1 \right] \\ {} & {} & -(i/2)\left( |-\rangle \langle +|-|+\rangle \langle -| \right){\rm sin} \phi \;{{{\rm e}}^{-1/2({{R}_{x}}(t)+{{R}_{\gamma }}(t))t}}\;. \\ \end{array}\end{eqnarray} \tag{ B.14 }$$

Knowing the density matrix time evolution it is straightforward to calculate the decay rate of the initial state $|{{\psi }_{\phi }}\rangle$ towards the completely mixed state ${{\rho }_{{\rm mixed}}}=1/2(|+\rangle \langle +|+|-\rangle \langle -|)$ and the probability for finding the system in the initial state after a time t

$$\begin{eqnarray}&&{\rm tr}\left( |{{\psi }_{\phi }}\rangle \langle {{\psi }_{\phi }}|\;\rho (t) \right)=\langle {{\psi }_{\phi }}|\rho |{{\psi }_{\phi }}\rangle =\frac{1}{2}\left( 1+{{{\rm cos} }^{2}}\phi \;{{{\rm e}}^{-{{R}_{x}}(t)t}}+{{{\rm sin} }^{2}}\phi \;{{{\rm e}}^{-1/2\left( {{R}_{x}}(t)+{{R}_{\gamma }}(t) \right)t}} \right)\end{eqnarray} \tag{ B.15 }$$

which corresponds exactly to the result given in (B.3).

B.3.1. Decay rate and general filter function formula

$$\begin{eqnarray}&&{{R}_{i}}(t)=\frac{1}{2\;t}\;\frac{1}{2\;\pi }\;\int _{-\infty }^{\infty }\;S(\omega )F_{t}^{i}(\omega ){\rm d}\omega \end{eqnarray} \tag{ B.16 }$$

with

$$\begin{eqnarray}&&S(\omega )=\int _{-\infty }^{\infty }\;{{{\rm e}}^{-{\rm i}\,\omega \,\tau }}\;\langle b(t)b(t+\tau )\rangle \;{\rm d}\tau \end{eqnarray} \tag{ B.17 }$$

and

$$\begin{eqnarray*}&&\quad F_{t}^{x}(\omega )={{\left| \int _{0}^{t}\;{{{\rm e}}^{-{\rm i}\,\omega \,t^{\prime} }}\;{\rm cos} \left( \int _{0}^{t^{\prime} }\;{\rm d}\tau \;\Omega (\tau ) \right){\rm d}t^{\prime} \right|}^{2}}+{{\left| \int _{0}^{t}\;{{{\rm e}}^{-{\rm i}\,\omega \,t^{\prime} }}\;{\rm sin} \left( \int _{0}^{t^{\prime} }\;{\rm d}\tau \;\Omega (\tau ) \right){\rm d}t^{\prime} \right|}^{2}}\end{eqnarray*}$$

$$\begin{eqnarray*}&&\quad F_{t}^{\gamma }(\omega )={{\left| \int _{0}^{t}\;{{{\rm e}}^{-{\rm i}\,\omega \,t^{\prime} }}\;{\rm cos} \left( \int _{0}^{t^{\prime} }\;{\rm d}\tau \;\Omega (\tau ) \right){\rm d}t^{\prime} \right|}^{2}}-{{\left| \int _{0}^{t}\;{{{\rm e}}^{-{\rm i}\,\omega \,t^{\prime} }}\;{\rm sin} \left( \int _{0}^{t^{\prime} }\;{\rm d}\tau \;\Omega (\tau ) \right){\rm d}t^{\prime} \right|}^{2}}\;.\end{eqnarray*}$$

B.3.2. Specific filter functions

B.3.3. Spectrum for the Ornstein Uhlenbeck process

$$\begin{eqnarray}&&\langle b(t)b(0)\rangle ={{b}^{2}}\;{{{\rm e}}^{-\frac{|t|}{\tau }}},\quad S(\omega )=\frac{2\;{{b}^{2}}\;\tau }{1+{{\omega }^{2}}\;{{\tau }^{2}}}\;\;.\end{eqnarray} \tag{ B.21 }$$

Here we will briefly discuss the coherence and relaxation timescales in the uncorrelated noise limit $t\gg \tau$, particularly relevant for the high frequency surface noise component. In that case the level transition rate, separated by an energy gap ω₀ and induced by the noise component k, is given by [22, 40]

$$\begin{eqnarray}&&r=\frac{1}{4}\;{{S}_{k}}({{\omega }_{0}})\end{eqnarray} \tag{ B.22 }$$

with ${{S}_{k}}(\omega )$ the corresponding noise spectrum. Here k denotes a noise component capable of inducing transitions, i.e. for relaxation processes in the undressed system it can be either the x or y component. Different components can be added and the final decay rates be obtained by means of classical rate equations [28, 33]. This situation is illustrated in figure B3 (a) for the T₁ relaxation and in figure B3(b) for the dressed state decay in a three level system as the NV-center gound state triplet. This latter case corresponds to the situation appearing in continuous dynamical decoupling.

Zoom In Zoom Out Reset image size

For the configuration (b) this leads to a rate equation for the state populations of the form

$$\begin{eqnarray}&&\frac{{\rm d}}{{\rm d}t}\left( \begin{array}{ccccccccccccccc} {{p}_{|+\rangle }} \\ {{p}_{|-\rangle }} \\ {{p}_{|-1\rangle }} \\ \end{array} \right)=\left( \begin{array}{ccccccccccccccc} -\gamma -{{r}_{0}} & \gamma & {{r}_{0}} \\ \gamma & -\gamma -{{r}_{0}} & {{r}_{0}} \\ {{r}_{0}} & {{r}_{0}} & -2\;{{r}_{0}} \\ \end{array} \right)\left( \begin{array}{ccccccccccccccc} {{p}_{|+\rangle }} \\ {{p}_{|-\rangle }} \\ {{p}_{|-1\rangle }} \\ \end{array} \right)\end{eqnarray} \tag{ B.23 }$$

with $\gamma ={{r}_{1}}+{{r}_{0}}$ the total transition rate between the dressed states induced by the noise z and y-component, respectively, ${{r}_{0}}=1/4\;{{S}_{x}}({{\omega }_{0}})=1/4\;{{S}_{y}}({{\omega }_{0}})$ and ${{r}_{1}}=1/4\;{{S}_{z}}(\Omega )$. Here ω₀ represents the transition frequency of the bare states and r₁ corresponds to the pure dephasing rate characterized by the spectral noise at the microwave Rabi frequency Ω. With ${{p}_{|+\rangle }}(0)=1$ this does lead to the solution

$$\begin{eqnarray}&&{{p}_{|+\rangle }}(t)=\frac{1}{3}+\frac{1}{6}\left[ {{{\rm e}}^{-3{{r}_{0}}\,t}}+3\;{{{\rm e}}^{-(2\gamma +{{r}_{0}})t}} \right]\simeq \frac{1}{2}\left[ 1+{{{\rm e}}^{-(2\gamma +{{r}_{0}})t}} \right]\end{eqnarray} \tag{ B.24 }$$

wherein the last approximation is valid for short enough times based on the realistic assumption that ${{r}_{0}}\ll {{r}_{1}}$. This allows one to identify the (decoupled) coherence time [28]

$$\begin{eqnarray}&&T_{2}^{\Omega }={{\left[ 2\;\gamma +{{r}_{0}} \right]}^{-1}}={{\left[ 2\;{{r}_{1}}+3\;{{r}_{0}} \right]}^{-1}}\end{eqnarray} \tag{ B.25 }$$

leading to (5). Note that the two level analogue can be obtained by retaining only γ in the above equations. An analogue calculation can be performed for the T₁ relaxation illustrated in figure B3(a), which leads to

$$\begin{eqnarray}&&{{T}_{1}}={{\left[ 6\;{{r}_{0}} \right]}^{-1}}\end{eqnarray} \tag{ B.26 }$$

(or ${{T}_{1}}={{[4{{r}_{0}}]}^{-1}}$ for a two level configuration). Comparing (B.26) to (B.25), and noting that r₁ is suppressed with increasing decoupling amplitudes Ω, the maximally achievable coherence time in a decoupled system is given by [22]

$$\begin{eqnarray}&&T_{2}^{\Omega }\leqslant 2\;{{T}_{1}}\;.\end{eqnarray} \tag{ B.27 }$$

An additional complication arises from the hyperfine structure associated with the nuclear spin of the nitrogen atom involved in the vacancy center (spin I = 1 for N-14 or spin $I=1/2$ for the N-15 isotope). Other nuclear spins will be neglected in the description, assuming that they can be efficiently described in the spin bath decoherence framework. The ground state level structure in the principal axis frame (analogue to (C.2)) is then given by [59, 60]

$$\begin{eqnarray}&&\quad {{H}_{{\rm gs}}}={{H}_{0}}+{{H}_{{\rm nucl}}}+{{H}_{{\rm hyp}}}\quad \quad \quad {\rm with}\;\;{{H}_{{\rm nucl}}}=P\;I_{z}^{^{\prime} 2}-{{\gamma }_{I}}\;\vec{B}\vec{I}^{\prime} ,\quad {{H}_{{\rm hyp}}}\simeq {{A}_{\parallel }}\;S_{z}^{^{\prime} }\;I_{z}^{^{\prime} }\end{eqnarray} \tag{ C.8 }$$

with H₀ the pure electron spin-Hamiltonian given by (C.2) and ${{H}_{{\rm nucl}}}$, ${{H}_{{\rm hyp}}}$ the (nitrogen) nuclear spin part and the hyperfine coupling interaction, respectively. Herein the quadrupole splitting $P=-4.95\;{\rm MHz}$, the gyromagnetic ratio ${{\gamma }_{I}}\simeq 3.1\;{\rm MHz}\;{{{\rm T}}^{-1}}$ and the coupling constant $A\simeq -2.16\;{\rm MHz}$. For the regime of our proposal (B = 0.5 T), this leads to a hyperfine-structure level scheme as depicted in figure C1. With the selection rule $\Delta {{m}_{I}}=0$ it becomes clear that neighboring transitions between two different electron spin states differ by $\Delta \omega \sim 2.2\;{\rm MHz}$. In order to avoid off-resonant microwave transitions that would lead to an imperfect decoupling and a different form of the effective dipolar coupling ${{H}_{I,{\mathscr{M}_{k}}}}$ as is most easily verified in the limit of large detunings, several strategies can be followed. (i) The most simple one consists of eliminating the hyperfine structure by sufficiently large microwave driving fields $\Omega \gg |{{A}_{\parallel }}|$. That way, the hyperfine coupling is suppressed, or analogously all possible nuclear states are excited simultaneously in the strong field limit, and there is thus no need to distinguish between individual hyperfine states. (ii) A second option, though not optimal here due to its insufficient noise decoupling strength, would be the regime of weak microwave driving $\Omega \ll |{{A}_{\parallel }}|$, in which only a single nuclear spin state transition is excited. (iii) A third possibility to achieve an effective two level system consists of polarizing the nuclear spin prior to the actual experiment. Together with the nuclear spin selection rules this results in a single possible transition and has been demonstrated as well at room [60–62] and at low temperatures [42, 63].

Zoom In Zoom Out Reset image size

Here we will outline the calculation of the noise spectrum arising from a bath of surface spins, following the approach developed in [22, 64]. This allows one to identify the relevant noise parameters, namely the root mean square noise amplitude and the noise correlation time, that have been used in the discussion of section 3. We will focus on a fluorine (nuclear spin) terminated surface and briefly discuss the influence of potential electronic surface spins at the end of this appendix.

Terminating the ND surface by fluorine, oxygen or hydrogen / hydroxyl groups replaces the electron spins associated with the sp2-hybridized orbitals by much weaker nuclear spins [29] and therefore reduces the dipolar interaction strength (see table D1 ) between surface spins by a factor of 10⁻⁶ (as given by the square ratio of the corresponding gyromagnetic ratios.) Additionally the coupling to the central spin, the electron spin of the NV center, reduces by a factor of 10⁻³. Recalling that the pure dephasing noise responsible for the decoherence mechanism can be described by a fluctuating magnetic field and that the strength and fluctuation rate depends on the (hyperfine) coupling to the central spin and the surface spin flip-flop rate, respectively, a significant improvement can be obtained by terminating the surface purely by nuclear spins. As an illustrative example we will concentrate on the fluorine terminated surface (nuclear spin 1/2), having the additional advantage of not affecting the charge state of the NV-center and forming a lattice with nearest neighbor distance $2.5\;$Å [35].

Table D1.  Magnitude for different coupling constants.

Mutual electron spin coupling:	Electron spin—fluorine	Mutual nuclear
	nuclear spin coupling:	spin fluorine coupling:
${{c}_{{\rm el},{\rm el}}}=\left( \frac{{{\mu }_{0}}\;\gamma _{{\rm el}}^{2}\;\hbar }{4\pi } \right)$	${{c}_{{\rm el},{\rm fl}}}=\left( \frac{{{\mu }_{0}}\;{{\gamma }_{{\rm el}}}{{\gamma }_{{\rm fl}}}\;\hbar }{4\;\pi } \right)$	${{c}_{{\rm fl},{\rm fl}}}=\left( \frac{{{\mu }_{0}}\;\gamma _{{\rm fl}}^{2}\;\hbar }{4\;\pi } \right)$
$52\;{\rm MHz}\;{{({\rm nm})}^{3}}$	$74.4\;{\rm kHz}\;{{({\rm nm})}^{3}}$	$106.3\;{\rm Hz}\;{{({\rm nm})}^{3}}$

Such a coupled system can be described by

$$\begin{eqnarray}&&H=H_{0}^{{\rm NV}}+H_{0}^{{\rm sf}}+H_{\operatorname{int}}^{{\rm sf}}+{{H}_{{\rm hf}}}\end{eqnarray} \tag{ D.1 }$$

with $H_{0}^{{\rm NV}}$ the NV-center energy Hamiltonian as given in (C.1), $H_{0}^{{\rm sf}}$ the magnetic field splitting of the surface nuclear spins, H_int^sf the dipolar coupling between surface nuclear spins and ${{H}_{{\rm hf}}}$ the hyperfine coupling of the surface spins to the central spin (the NV-center electron spin):

$$\begin{eqnarray}\begin{array}{rcl} H_{0}^{{\rm sf}} & = & \mathop{\sum }\limits_{i}\frac{{{\gamma }_{n}}\;B}{2}\;\sigma _{i}^{z} \\ H_{\operatorname{int}}^{{\rm sf}} & \simeq & \mathop{\sum }\limits_{i\gt j}\left( \frac{{{\mu }_{0}}}{4\;\pi }\;\frac{\gamma _{n}^{2}\;\hbar }{r_{ij}^{3}} \right)(3{{{\rm cos} }^{2}}{{\theta }_{ij}}-1)\left( -\sigma _{i}^{z}\;\sigma _{j}^{z}+\frac{1}{2}\left[ \sigma _{i}^{x}\;\sigma _{j}^{x}+\sigma _{i}^{y}\sigma _{j}^{y} \right] \right) \\ {{H}_{{\rm hf}}} & \simeq & \mathop{\sum }\limits_{i}\left( \frac{{{\mu }_{0}}}{4\;\pi }\;\frac{{{\gamma }_{{\rm el}}}\;{{\gamma }_{n}}\;\hbar }{r_{i}^{3}} \right)(1-3\;{{{\rm cos} }^{2}}{{\theta }_{i}})\;\;{{S}^{z}}\;\sigma _{z}^{i}. \\ \end{array}\end{eqnarray} \tag{ D.2 }$$

Herein S denotes the spin-1 operator of the vacancy center, σ_k the spin 1/2 operators of the surface nuclear spins, r_ij the distance between spins i and j, r_i that between surface spin i and the central spin and θ_ij and θ_i the angles between the vector connecting the two coupled spins involved and the external magnetic field. As discussed in appendix C we will assume that the quantization axis of the NV-center is essentially determined by the external magnetic field and not by the symmetry axis of the vacancy center. Alternatively one might assume that the external magnetic field is parallel to the NV symmetry axis. In the two level approximation, used for the driven system in appendix C, one can substitute ${{S}^{z}}\to 1/2({{s}_{z}}\pm 1)$ with s_z the spin 1/2 operator of the quasi-spin. Herein only secular contributions have been retained, and second order nuclear spin flip-flop processes [65–67] are orders of magnitude smaller than the direct ones for the dense spin bath considered and therefore can be neglected.

Calculating the noise spectrum and the corresponding noise parameters is in general intractable for a large number of surface spins due to exponentially increasing computational resources. As early as in 1962 Klauder and Anderson showed by very general arguments that a Lorentzian noise distribution has to be expected in the case of dipolar couplings to a spin bath [68]. Since then various approaches, mean-field and exact approaches in certain limits, have been invented to infer the noise properties and to calculate the decoherence decay rate, ranging from the simple Ornstein–Uhlenbeck model [69] to cluster expansion methods [65, 67, 70, 71]. Here we used the mean-field method described in [22, 64] which corresponds to the lowest order of the cluster expansion methods, the so called pair-correlation approximation, corrected by a mean-field broadening using the method of moments [72]. In the pair-correlation approximation it is assumed that each of the flipping processes in the dipolar coupling (D.2) is independent of all other processes, i.e. the Hilbert-space is substituted by a pair Hilbert space $(ij)$ wherein (${{H}_{ij}}$ follows directly from out of (D.2))

$$\begin{eqnarray}&&H_{\operatorname{int}}^{{\rm sf}}=\mathop{\sum }\limits_{i\gt j}{{\tilde{H}}_{ij}}=\mathop{\sum }\limits_{(ij)}{{\tilde{H}}_{(ij)}}\quad {\rm with}\;\;[{{H}_{ij}},{{H}_{kl}}]=0\;\;\forall ij,kl\;.\end{eqnarray} \tag{ D.3 }$$

This leads to a discrete spectrum of δ-peaks at the different pair-induced transition frequencies and can be justified as long as correlations between those pairs can be neglected, that is as long as the evolution time from an initial thermal state is short enough such that $1-{\rm exp} (-q\;N_{{\rm flip}}^{2}/N)$ is small [67] (with ${{N}_{{\rm flip}}}$ the number of flipped bath spins during the considered time, q the number of nearest neighbors and N the total number of bath spins.) Higher orders would lead to additional frequency peaks as well as couplings of the already existing peaks, i.e. a finite lifetime broadening which is taken into account by using a mean-field type approach based on the theory of moments. From (D.2) and (B.1) the operator for the effective field on the NV-center follows to be

$$\begin{eqnarray}&&\hat{b}=\mathop{\sum }\limits_{i}\left( \frac{{{\mu }_{0}}}{4\;\pi }\;\frac{{{\gamma }_{{\rm el}}}\;{{\gamma }_{n}}\;\hbar }{r_{i}^{3}} \right)(1-3\;{{{\rm cos} }^{2}}{{\theta }_{i}})\sigma _{i}^{z}\end{eqnarray} \tag{ D.4 }$$

and the noise spectrum can be calculated by (see section B.1)

$$\begin{eqnarray}&&S(\omega )=\int _{-\infty }^{\infty }\;{\rm d}t\;\langle \hat{b}(t)\hat{b}(0)\rangle \;{{{\rm e}}^{{\rm i}\omega \,t}}\;\;,\end{eqnarray} \tag{ D.5 }$$

where in good approximation the initial bath states can be assumed to be uncorrelated with the NV-center and at room temperature given by $1/{{2}^{N}}$ [69]. Following the calculation presented in [22] this leads to

$$\begin{eqnarray}&&S(\omega )=2\;\pi \;\mathop{\sum }\limits_{i\lt j}\frac{b_{ij}^{2}\;\Delta _{ij}^{2}}{b_{ij}^{2}+\Delta _{ij}^{2}}\;\;\frac{1}{\sqrt{2\;\pi \;\sigma _{ij}^{2}}}\left( {\rm exp} \left[ -\frac{{{(\omega -{{E}_{ij}})}^{2}}}{2\;\sigma _{ij}^{2}} \right]+{\rm exp} \left[ -\frac{{{(\omega +{{E}_{ij}})}^{2}}}{2\;\sigma _{ij}^{2}} \right] \right)\end{eqnarray} \tag{ D.6 }$$

with

$$\begin{eqnarray}\begin{array}{rcl} {{b}_{ij}} & = & \left( \frac{{{\mu }_{0}}}{4\;\pi }\;\frac{\gamma _{n}^{2}\;\hbar }{r_{ij}^{3}} \right)(3{{{\rm cos} }^{2}}{{\theta }_{ij}}-1) \\ {{\Delta }_{ij}} & = & \frac{1}{4}\left( {{A}_{i}}-{{A}_{j}} \right)\quad {\rm with}\quad {{A}_{i}}=2\left( \frac{{{\mu }_{0}}}{4\;\pi }\;\frac{{{\gamma }_{{\rm el}}}\;{{\gamma }_{n}}\;\hbar }{r_{i}^{3}} \right)(1-3\;{{{\rm cos} }^{2}}{{\theta }_{i}}) \\ \end{array}\end{eqnarray} \tag{ D.7 }$$

and ${{E}_{ij}}=2\sqrt{b_{ij}^{2}+\Delta _{ij}^{2}}$ and σ_ij the mean-field broadening used in replacing the original delta-peaks by Gaussian peaks leading to the same second moment as that obtained by the moment theory

$$\begin{eqnarray}&&\sigma _{ij}^{2}=\frac{b_{ij}^{2}+\Delta _{ij}^{2}}{4\;\Delta _{ij}^{2}\;b_{ij}^{2}}\;\mathop{\sum }\limits_{k\ne i,j}\left( b_{ik}^{2}\;A_{i}^{2}+b_{jk}^{2}\;A_{j}^{2} \right).\end{eqnarray} \tag{ D.8 }$$

The spectrum (D.6) has been numerically evaluated by randomly distributing nuclear surface spins on a sphere with equal inter-spin distance and placing the NV-center in the center of the sphere, i.e. for a sphere radius $r=5\;{\rm nm}$ a number of 5026 (which corresponds to a distance of $2.5\;$Å) nuclear spins were distributed on the surface. Assuming an equal quantization direction for both the NV-center and the nuclear spins (corresponding to a strong magnetic field or to a magnetic field parallel to the NV symmetry axis) the noise spectrum was calculated and subsequently fitted to a Lorentzian in order to obtain the noise amplitude b and correlation time τ defined in (B.21). Depending on the number of surface spins an average over different random positions was performed.

We performed the same calculation as well for different numbers of electron surface spins (see figure D1 ) on diamond with radius $5\;{\rm nm}$. In that case the decoherence parameters are very poor, e.g. ${{T}_{2}}=0.1\;\mu {\rm s}$, $\tau =21.7\;ns$ and $b=3.6\;{\rm MHz}$ for 30 surface spins and getting worse on further increasing the number of spins. For such a situation decoupling fields in the $\Omega \sim 1/\tau \sim {\rm GHz}$ range are required that, despite still allowing simple two qubit gates (provided the Ω stability can be achieved), are intractable for multi-qubit applications that require a certain level of individual addressing. As the validity of the mean-field approach for fast flipping electron spins is debatable, we performed an exact numerical T₂-calculation for smaller electron spin numbers (up to 10) showing that the mean-field values seem at least to be correct in their order of magnitude.

Zoom In Zoom Out Reset image size

Based on the noise parameters as obtained in section 3 and (B.21), the numerical simulations have been performed assuming an Ornstein–Uhlenbeck (Gaussian, Markovian and stationary) noise process [26] and its numerical implementation as described in [6, 73].

For such a process the evolution of the random field is governed by the Langevin equation

$$\begin{eqnarray}&&\frac{{\rm d}b(t)}{{\rm d}t}=-\frac{1}{\tau }\;b(t)+\sqrt{c}\;\;\Gamma (t)\end{eqnarray} \tag{ D.9 }$$

with $\Gamma (t)$ a Gaussian zero-mean noise [6] ($\langle \Gamma (t)\rangle =0$, $\langle \Gamma (t)\Gamma (t^{\prime} )\rangle =\delta (t-t^{\prime} )$). The parameter c appearing in (D.9) is related to the noise correlation time τ and field variance b² as defined in (1) by ${{b}^{2}}=c\;\tau /2$. Exact updating formulas can be obtained for the differential equation (D.9) [73]

$$\begin{eqnarray}&&b(t+\Delta t)=b(t)\mu +{{\xi }_{x}}\;{{n}_{1}}\quad {\rm with}\quad \mu ={\rm exp} (-\Delta t/\tau ),\qquad \xi _{x}^{2}=(c\;\tau /2)(1-{{\mu }^{2}}),\end{eqnarray} \tag{ D.10 }$$

and n₁ a unit normal random number. Introducing (D.10) into the corresponding Hamiltonians directly allows for a numeric integration of the associated equations of motion, followed by an averaging process that reveals the correct influence of decoherence.

The two qubit gate interaction corresponds essentially to a single qubit rotation in a more complicated two qubit manifold ${\mathscr{M}_{1}}$ and ${\mathscr{M}_{2}}$, respectively. Introducing the (unknown) systematic error this results, assuming $\Omega \gg J$, in the following Hamiltonians

$$\begin{eqnarray}&&{{H}_{I,{\mathscr{M}_{k}}}}\simeq \frac{J}{2}(1+\epsilon )\;S_{{\mathscr{M}_{k}}}^{1,2}\end{eqnarray} \tag{ E.1 }$$

which can also be written as ${{H}_{I,{\mathscr{M}_{k}}}}\simeq (J/2)(1+\epsilon )\sigma _{x}^{{\mathscr{M}_{k}}}$ with $\sigma _{l}^{{\mathscr{M}_{k}}}$ the σ_l operation defined in the manifolds $\{|+-\rangle ,|-+\rangle \}$ and $\{|++\rangle ,|--\rangle \}$ for k = 1 and k = 2, respectively (see figure 3(a) in the main text).

References [44, 74] provide a method to remove the systematic error contribution based on noting that multiples of $2\;\pi$ pulses with respect to the ideal error-less gate end up with pure contributions which can, together with the property $\sigma _{\phi }^{{\mathscr{M}_{k}}}+\sigma _{-\phi }^{{\mathscr{M}_{k}}}=2\;{\rm cos} \phi \;\sigma _{x}^{{\mathscr{M}_{k}}}$ ($\sigma _{\phi }^{{\mathscr{M}_{k}}}:={\rm cos} \phi \;\sigma _{x}^{{\mathscr{M}_{k}}}+{\rm sin} \phi \;\sigma _{y}^{{\mathscr{M}_{k}}}$), be used to create a compensation cycle by means of a Suzuki–Trotter time addition. A compensated gate is then obtained by the sequence (herein the compensation sequence is marked by the square brackets, figure E1 (b))

$$\begin{eqnarray}&&\left[ M_{\phi }^{\{\epsilon \}}(\pi )M_{3\phi }^{\{\epsilon \}}(2\;\pi )M_{\phi }^{\{\epsilon \}}(\pi ) \right]\;{{M}^{\{\epsilon \}}}(\theta )={{M}^{\{\epsilon =0\}}}(\theta )+\mathscr{O}({{\epsilon }^{3}})\end{eqnarray} \tag{ E.2 }$$

with

$$\begin{eqnarray}&&{{M}^{\{\epsilon \}}}(\theta )={\rm exp} \left( -{\rm i}(\theta /2)(1+\epsilon )\sigma _{x}^{{\mathscr{M}_{k}}} \right)={\rm exp} \left( -{\rm i}\;{{H}_{I,{\mathscr{M}_{k}}}}\;\theta /J \right),\quad 2\;\pi \;{\rm cos} \phi =-\theta /2\end{eqnarray} \tag{ E.3 }$$

and $M_{\phi }^{\epsilon }(\theta )={{T}_{\phi }}\;{{M}^{\{\epsilon \}}}(\theta )T_{\phi }^{\dagger }$ with ${{T}_{\phi }}={\rm exp} (-{\rm i}\phi /2\;\sigma _{z}^{{\mathscr{M}_{k}}})$. Note that

$$\begin{eqnarray}&&\sigma _{z}^{{\mathscr{M}_{1|2}}}=\left\{ 1/2(\sigma _{x}^{1}\mp \sigma _{x}^{2}),\sigma _{x}^{1},\sigma _{x}^{2} \right\}\end{eqnarray} \tag{ E.4 }$$

with the upper and lower sign referring to ${\mathscr{M}_{1}}$ and ${\mathscr{M}_{2}}$, respectively, and the last two options having additional contributions that however do not affect the gate manifold.

Zoom In Zoom Out Reset image size

Assuming systematic errors _ij in the dipolar coupling term resulting in the Hamiltonian

$$\begin{eqnarray}&&H_{zz}^{\{\epsilon \}}=\frac{J}{2}\;\mathop{\sum }\limits_{\frac{i\gt j}{i\in {\rm neighb}(j)}}(1+{{\epsilon }_{ij}})\sigma _{i}^{z}\;\sigma _{j}^{z}\end{eqnarray} \tag{ E.5 }$$

leads, besides the desired contribution, to additional _ij dependent terms in the defective evolution operator

$$\begin{eqnarray}&&\quad {{M}^{\{\epsilon \}}}(\theta )=\mathop{\prod }\limits_{(i,j)}M_{i,j}^{\epsilon }\quad {\rm with}\;\;M_{ij}^{\epsilon }(\theta )={\rm exp} \left( -{\rm i}\;\frac{\theta }{2}\;\sigma _{z}^{i}\;\sigma _{z}^{j} \right)\;{\rm exp} \left( -{\rm i}\;\frac{\theta }{2}\;{{\epsilon }_{ij}}\;\sigma _{z}^{i}\;\sigma _{z}^{j} \right).\end{eqnarray} \tag{ E.6 }$$

A (compensation) contribution only involving the systematic error part can be obtained by any multiple of $2\;\pi$-pulses $M_{ij}^{\epsilon }(n\;2\pi )={{(-1)}^{n}}\;{\rm exp} \left( -{\rm i}\;n\pi \;{{\epsilon }_{ij}}\;\sigma _{z}^{i}\;\sigma _{z}^{j} \right)$ with $n\in \mathbb{Z}$. This property forms the basis for creating a compensation sequence provided that the relatively fixed phase, limited by the $2\pi$ condition, can be controlled. This latter prerequisite is achieved using that

$$\begin{eqnarray}&&\left( \sigma _{\phi }^{i}+\sigma _{-\phi }^{i} \right)\sigma _{z}^{j}={{T}_{\phi }}\left( \sigma _{z}^{i}\;\sigma _{z}^{j} \right)T_{\phi }^{\dagger }+{{T}_{-\phi }}\left( \sigma _{z}^{i}\;\sigma _{z}^{j} \right)T_{-\phi }^{\dagger }=2\;{\rm cos} \phi \;\sigma _{z}^{i}\;\sigma _{z}^{j}\end{eqnarray} \tag{ E.7 }$$

with ${{T}_{\phi }}={\rm exp} (-{\rm i}\;\phi /2\;\sigma _{x}^{i})$ and ${{\sigma }_{\phi }}={\rm cos} \phi \;{{\sigma }_{z}}-{\rm sin} \phi \;{{\sigma }_{y}}$. Performing a Suzuki–Trotter addition of the two components ${{\sigma }_{\phi }}$ and ${{\sigma }_{-\phi }}$ therefore allows for the construction of a compensated sequence

$$\begin{eqnarray}&&\quad \left[ M_{ij,\phi }^{\epsilon }(n\;2\;\pi )M_{ij,-\phi }^{\epsilon }(n\;4\;\pi )M_{ij,\phi }^{\epsilon }(n\;2\;\pi ) \right]\;M_{ij}^{\epsilon }(\theta )={\rm exp} \left( -{\rm i}\;\frac{\theta }{2}\;\sigma _{z}^{i}\;\sigma _{z}^{j} \right)+\mathscr{O}(\epsilon _{ij}^{3})\end{eqnarray} \tag{ E.8 }$$

with $4\;n\;\pi \;{\rm cos} \phi =-\theta /2$ and $M_{ij,\phi }^{\epsilon }(\theta )={{T}_{\phi }}\;M_{ij}^{\epsilon }(\theta )T_{\phi }^{\dagger }$. Taking into account that the error scales as $\mathscr{O}(n\;\epsilon _{ij}^{3})$ the parameter n should be as small as possible, leading to the obvious choice n = 1. This can then be straightforwardly extended to the total evolution ${{M}^{\{\epsilon \}}}$ as defined in (E.6), exhibiting the same form as (E.8) and ${{T}_{\phi }}={\rm exp} \left( -{\rm i}\;\phi /2\;{{\sum }_{i|i\in {{\mathscr{N}}^{c}}}}\sigma _{i}^{x} \right)$ and ${\mathscr{N}^{c}}$ the space of non-neighboring qubits. As an example the method is illustrated for the creation of a four qubit cluster state in figure E1(a), showing that high decoupling fields are required in the presence of decoherence in order to benefit from the compensation despite the significant increased gate time that is problematic in terms of decoherence.

As outlined above, smaller values for n lead to an improved scaling of the compensation method, i.e. to a smaller remaining error contribution and in the case when decoherence is relevant, to a significant reduction of the compensated gate time. Favorably, the single and two particle case allows for the construction of a compensation sequence with $n=1/2$ (E.2), which is beyond the $n2\pi$ analysis outlined above in that the independent contributions survive in the first stage. However this concept cannot be extended to multiple particles; a detailed analysis shows that for such a $n=1/2$ sequence to be valid, the rotation operator ${{T}_{\phi }}$ has to commute with the remaining -independent contribution involving $\phi \to -\phi$, which is valid for an odd but not for an even commutation, the latter case relevant for the multiparticle situation

$$\begin{eqnarray}&&{{{\rm e}}^{-{\rm i}\pi /2{{\sigma }_{\phi }}}}\;{{{\rm e}}^{-{\rm i}\phi {{\sigma }_{x}}}}={{{\rm e}}^{+{\rm i}\phi {{\sigma }_{x}}}}\;{{{\rm e}}^{-{\rm i}\pi /2{{\sigma }_{\phi }}}}\leftrightarrow \;{{{\rm e}}^{-{\rm i}\pi /2\left[ \sigma _{z}^{1}\sigma _{\phi }^{2}+\sigma _{\phi }^{2}\sigma _{z}^{3} \right]}}\;{{{\rm e}}^{-{\rm i}\phi \sigma _{x}^{2}}}={{{\rm e}}^{-{\rm i}\phi \sigma _{x}^{2}}}\;{{{\rm e}}^{-{\rm i}\pi /2\left[ \sigma _{z}^{1}\sigma _{\phi }^{2}+\sigma _{\phi }^{2}\,\sigma _{z}^{3} \right]}}.\end{eqnarray} \tag{ E.9 }$$

The Hamiltonians (E.1) and (E.5) used to describe the compensation mechanism have been defined up to next-nearest neighbors, i.e. higher order contributions were neglected so far. However, whereas those contributions give only a small correction for simple $\pi /2$ multiqubit pulses, as can be seen from the numerical results in the main text, they are crucial in the compensation sequence due to the prolonged pulse time. One has to distinguish between two types of higher order couplings: even couplings defined as next-nearest neighbor and third, fifth, ...nearest neighbor couplings and odd couplings as second, fourth, ...nearest and diagonal couplings. Whereas the first ones have always the form of a ${{\sigma }_{z}}\otimes {{\sigma }_{z}}$ coupling and are automatically compensated for by the next-nearest neighbor based compensation sequence, the second class consists just of the reduced manifold ${\mathscr{M}_{k}}$ couplings and is not removed by the compensation sequence. This failure for odd couplings originates both from the fact that the compensation pulse acts on both qubits involved such that (E.7) is not fulfilled any more, and from the fact that the manifold interaction does not commute with the ${{\sigma }_{z}}{{\sigma }_{z}}$ type contributions such that a splitting of the evolution in a perfect and defective part as in (E.6) is not possible any more.

A solution to that problem consists of eliminating odd order couplings from the beginning by adding different contributions in time. Examples for that are given for a two qubit state in figure 4(b) in the main text and in figure E2 for second nearest neighbor couplings. With the evolutions obtained that way as a starting point higher orders can be removed in a concatenated way benefiting from the different evolution timescales dependent on the qubit distance. For the multiqubit compensation sequence (E.8) it is sufficient to remove odd order couplings up to the second order (e.g. up to couplings of qubit one to five), because higher order effects are too weak to make significant contributions.

Zoom In Zoom Out Reset image size