Self-assembling hybrid diamond-biological quantum devices

The realization of scalable arrangements of nitrogen vacancy (NV) centers in diamond remains a key challenge on the way towards efficient quantum information processing, quantum simulation and quantum sensing applications. Although technologies based on implanting NV-center in bulk diamond crystals or hybrid device approaches have been developed, they are limited in the achievable spatial resolution and by the intricate technological complexities involved in achieving scalability. We propose and demonstrate a novel approach for creating an arrangement of NV-centers, based on the self-assembling capabilities of biological systems and its beneficial nanometer spatial resolution. Here, a self-assembled protein structure serves as a structural scaffold for surface functionalized nanodiamonds, in this way allowing for the controlled creation of NV-structures on the nanoscale and providing a new avenue towards bridging the bio-nano interface. One-, two- as well as three-dimensional structures are within the scope of biological structural assembling techniques. We realized experimentally the formation of regular structures by interconnecting nanodiamonds using biological protein scaffolds. Based on the achievable NV-center distances of 11nm, we evaluate the expected dipolar coupling interaction with neighboring NV-center as well as the expected decoherence time. Moreover, by exploiting these couplings, we provide a detailed theoretical analysis on the viability of multiqubit quantum operations, suggest the possibility of individual addressing based on the random distribution of the NV intrinsic symmetry axes and address the challenges posed by decoherence and imperfect couplings. We then demonstrate in the last part that our scheme allows for the high-fidelity creation of entanglement, cluster states and quantum simulation applications.

The realization of scalable arrangements of nitrogen vacancy (NV) centers in diamond remains a key challenge on the way towards efficient quantum information processing, quantum simulation and quantum sensing applications.
Although technologies based on implanting NV-center in bulk diamond crystals [1] or hybrid device approaches [2] have been developed, they are limited in the achievable spatial resolution and by the intricate technological complexities involved in achieving scalability. We propose and demonstrate a novel approach for creating an arrangement of NV-centers, based on the self-assembling capabilities of biological systems and its beneficial nanometer spatial resolution [3,4]. Here, a selfassembled protein structure serves as a structural scaffold for surface functionalized nanodiamonds, in this way allowing for the controlled creation of NV-structures on the nanoscale and providing a new avenue towards bridging the bio-nano interface. One-, two-as well as three-dimensional structures [5] are within the scope of biological structural assembling techniques. We realized experimentally the formation of regular structures by interconnecting nanodiamonds using biological protein scaffolds. Based on the achievable NV-center distances of 11nm, we evaluate the expected dipolar coupling interaction with neighboring NV-center as well as the expected decoherence time. Moreover, by exploiting these couplings, we provide a detailed theoretical analysis on the viability of multiqubit quantum operations, suggest the possibility of individual addressing based on the random distribution of the NV intrinsic symmetry axes and address the challenges posed by decoherence and imperfect couplings. We then demonstrate in the last part that our scheme allows for the high-fidelity creation of entanglement, cluster states and quantum simulation applications.
Coupled solid-state spin systems as electron spins in quantum dots, phosphorous donors in silicon and color centers in diamond form promising candidates for the emerging field of quantum technologies [6]. Among those the negatively charged nitrogen-vacancy center (NV − ) in diamond [7], composed of a substitutional nitrogen atom and an adjacent vacancy within the diamond lattice, subject of this work, stands out due to its long coherence time up to milliseconds at room temperature [8]. It forms a discrete atom-like energy level structure within the diamond bandgap, with the ground state described by an electron-spin triplet (Spin-1) as illustrated in figure 2 (a), that allows for the full coherent control by means of e.g. microwave drivings and static magnetic fields [9]. Remarkably, readout and initialization through the excited state can be performed optically, taking benefit of the spin-dependent fluorescence rates and an intersystem crossing, the latter one allowing for the high-fidelity state preparation by optical pumping [9]. Interaction among neighbouring NV-centers can be mediated by the magnetic dipolar coupling of the electronic spins [10,11], yet this requires distances of the order of 10 nm or below to enable coherent interaction strengths that comfortably exceed the decoherence times. However current techniques for the controlled positioning of NV-centers within bulk crystals, i.e. the creation of NV-centers by ion implantation, are limited to several tens of nanometers in position [12]. In contrast, the ability of biological systems for structural self-assembly [3,4] is a powerful tool allowing for the simple and parallel creation of large ordered arrays on nanometer scales, that holds the potential for outperforming the limited resolution and structural complexity achievable of conventional lithography based on serial pattern creation. We propose to combine the self-assembly of biological systems with the guided attachment of surface functionalized nanodiamonds. The biomolecules are used as a structural scaffold to enable the formation of NV-center configurations with high spatial resolution. This can be achieved with tiled motifs such as short DNA strands [13] or membrane forming complexes including SP1 [14], LH1 [15] and TF55β [16], that allow for the creation of one and two dimensional arrays. Going beyond two-dimensional periodic patterns, the method of DNA-origami [5,17], based on folding a large single stranded DNA molecule directed by staple strands, enables the creation of more complex highly controllable structures, ranging from aperiodic arrays to real three dimensional structures. Elaborate knowledge in genetic engineering now allows control over the assembly, structure and topology of biomolecular complexes as well as the attachment of nanoparticles with nanometer precision [18]. Each of these structures is suitable for hosting nanodiamonds, whose size can be controlled down to 4 nm in fabrication [19], and whose chemical attachment and site-specific binding can be directed by surface functionalization and labeling. This allows for the high precision positioning of NV-centers incorporated in such diamonds, which subsequently interact via dipolar couplings, and paves the way for a highly controllable array of interacting quantum systems.
Formation of SP1 nanodiamond structures -To demonstrate the feasibility of interconnecting nanodiamonds (NDs) with biological structures, we present the controllable and repeatable formation of small nanodiamond complexes using an SP1 (Stable Protein 1, see figure 1(c) and Supplement section 1) protein variant and first steps towards the formation of regular arrays of NDs on SP1 arrays. A crucial first step to enable both experiments is the genetic modification of SP1 to fuse 12 graphite-specific binding peptides to the SP1 N-terminus, that permit site specific binding of SP1 to the carbon sp2-hybridization which forms on the nanodiamond surface [20].
Site-specific binding of the SP1 is essential for small NDs (under 10 nm) to form regular structures on an SP1 array (figure 1 d) [21]. As an important result towards the creation of large ordered NDs structures, we achieved the formation of numerous dimers and trimers along with larger ordered structures such as a 7 NDs hexagon as illustrated in figure 1 (a)&(b). Here a monolayer of genetically modified SP1 was formed by the Langmuir-Blodgett method [22] and subsequently combined with a ND solution (5nm in diameter formed by laser ablation). Using both the SP1-template and a diluted ND solution, ordering of the NDs partially filling the SP1 template has been achieved. Such isolated structures are promising candidates for NV-coupling experiments. Figure 1 (a) shows an image of ND structures on the SP1 monolayer and figure 1 (b) is enlarging the ND hexagon. In this SP1/NDs sample, the size of the measured nano-particles is around 5 nm in diameter, therefore matching the expected NDs size. On the other hand, the distance between centers of adjacent particles corresponds to 11nm, which is the distance between the SP1 proteins on an ordered layer. Moreover, electron diffraction on this sample in the relevant areas proved the existence of diamonds on the surface (see Supplementary section 1).
Two references were used to ensure that the combination of the SP1 template and the NDs is essential for achieving the uniform spacing between adjacent NDs. A sample without the SP1 template that contains only adsorbed NDs, in which case only irregular ND aggregates were observed in the TEM measurements, and a sample with the SP1 template but without the NDs in which case we did not find any nano-particles on the sample.
In a different approach, small ND structures were also achieved using larger (average diameter 30 nm created by grinding) NDs. The ND clusters were formed by mixing solutions of NDs and SP1 under ambient conditions. The average number of NDs in such a complex can be controlled by the concentration ratio of SP1 to NDs. Mixing a 1mg/ml ND solution with a 1mg/ml SP1 solution at a ratio of 1:1 mainly leads to the formation of dimers and trimers and for even higher concentrations of SP1 large diamond clusters can be observed (figure 1 (e)). As the NDs are larger than the SP1, the exact structure of the clusters is controlled by the nanodiamond shape as several SP1 will bind the NDs across a surface. To verify that the creation of the clusters is due to the binding with SP1 rather than electrostatic forces we have also followed the same procedure in the absence of SP1. In this case no ND clusters are observed.
Decoherence and dynamical decoupling -A key challenge for quantum computation and simulation with self-assembled arrays of nanodiamonds are the currently achieved relatively short coherence times of a few µs [19] compared to the typical coupling strength of several kHz between NV-centers in different nanodiamonds [11,23]. This is originated in interactions of the NV-center with the very proximate surface spins [24] and external charges, such that substantial improvements can be obtained by means of surface functionalization methods [19,25]. The remaining decoherence after optimization of material design can be further improved by applying dynamical decoupling techniques, that allow to increase the coherence time by several orders of magnitude. These may be implemented either as pulsed schemes [26] or by continuous driving fields [11,[27][28][29]. At this point we would like to point out that beside surface spins, other impurities well-known from bulk diamond such as spins of C-13 (type-IIa diamond) and P1-centers of nitrogen donors (type-Ib diamond), contribute to the effect of decoherence as well. This influence however is much less dominant as their low concentration leads to a rather small and only weakly coupling spin bath [19]. Moreover it can be described and thus decoupled in the same framework as outlined below and consequently we do not expect those additional impurities to change the results presented here. The pure dephasing noise can be modelled as a fluctuating magnetic field energy shift b(t) with zero mean and b(t) b(0) = b 2 e −t/τ (see Supplement section 5). Adding a continuous driving interaction for decoupling, i.e. a resonant microwave coupling Ω on the relevant electron spin transition, the Hamiltonian of the effective two level system in a frame rotating with the laser frequency for a single NV-center is given by H = (b(t)/2) σ z + (Ω/2) σ x , where σ x,z are the usual spin-1/2 Pauli operators. In the relevant Markovian limit t τ and Ω t > 1 the decoherence decay rate R(t) is determined by the noise spectrum S(ω) evaluated at the decoupling frequency, i.e. R(t) ∝ S(Ω) and the more general formalism will be discussed in the Supplement (section 2) [30]. This allows to define an effective T 2 -time of the system resulting in for the considered Lorentzian noise spectrum. Figure 2 (b) compares this formula to numerical noise simulations showing that both the T 2 -scaling and the exponential decay behaviour are in good agreement with numerical calculations. Here the noise spectrum was calculated for a fluorine terminated surface of a spherical diamond with radius r = 5 nm that results in a fluorine nuclear spin-1/2 surface with a nearest neighbour distance of 2.5Å, using the mean-field approach described in [31] (see Supplement, section 4). We thus expect a noise correlation time of τ = 2.5µs and a mean square root amplitude b = 30.2 kHz (τ ∝ 1/(n 3/2 ) and b 2 ∝ n 2 r with n the surface spin density and r the nanodiamond radius), that could be further improved by spin bath polarization and decoupling schemes [25,32]. For those parameters, a decoupling field of Ω = 1.2 MHz leads to an effective coherence time of T 2 4 ms comparable to the electron spin T 1 -time that imposes an upper limit on the dynamical decoupling. In experiments it is possible to achieve much higher decoupling fields up to 300 MHz such that even a small number of possible electron spins can be decoupled (for a discussion of electron spin noise see Supplement, section 4). Combining self-assembled structures with lithographic techniques might help to integrate current structures within the system [33] beneficial especially for high driving fields. Importantly, decoherence due to intensity fluctuations of the decoupling field can be suppressed by using a concatenated decoupling scheme as proposed in [28].
Engineering of interactions, spin gates and quantum simulation -Dipolar interactions between electron spins of two adjacent NV-centers provide a possibility for implementing gate operations [10,11]. 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 Supplement, section 3) 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 ξ ij µ ij , where ξ ij depends on the external magnetic field strength | B| and its orientation with respect to both the NV symmetry axes and the vector r ij connecting NVcenter i and j. In the limit of high magnetic fields ξ ij = 1/4 (1 − 3 cos 2 θ ij ), where θ ij = ∠( r ij , B) and µ ij = µ 0 γ 2 el /(4πr 3 ij ) with the latter being µ ij = 52 kHz for an NV-center distance r ij = 10 nm. µ 0 and γ el denote the magnetic permeability and the electron gyromagnetic ratio, respectively. It is important to note that the configuration of nanodiamonds 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 (c) Dipolar coupling parameter ξ vs the external magnetic field. Red circles denote the average and blue lines the variance arising from the random axis orientation. For large magnetic field the quantization axis is given by the external field and ξ 1/4. The green dashed area indicates the range that can be corrected using compensation methods. All parameters are obtained by averaging over 10 6 random spin orientations. Inset: Optimal magnetic field configuration for a 2D array as also used for the main graph.
splitting of the ground state triplet (D ∼ 2.87 GHz) (see figure 2(a)). This leads to two important consequences: 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 2 (c) together with its optimal 2D-orientation. As shown in the figure, this can be achieved applying a magnetic field B 0.5 T (γ e B/(2 π) 14 GHz), in which case the magnetic energy shift dominates over the crystal-field splitting (γ e B D). Second, the crystal-field splitting is responsible for an orientation dependent transition frequency in that regime, depending on the angle between the symmetry axes and the external field ϑ i and scaling as ∝ D cos 2ϑ i . As a consequence, the transition frequencies of individual NV-centers differ by typical values of several 100 MHz, thus providing the possibility for individual microwave addressing for the values of Ω obtained from figure 2 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 (2), that would not be energy conserving in our setup.
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. Transforming to an interaction picture with respect to the driving in the relevant limit Ω J ij , the effective interaction up to non-nearest neighbours is given by , wherein (i, j) sums over all neighbouring NV-centers and M 1 corresponds to the situation of a homogeneous decoupling field Ω i = Ω, whereas M 2 corresponds to the situation Ω i = −Ω j for i ∈ neighb(j). s + and s − are the ladder operators in the σ x -eigenbasis.
The two qubit setting is shown in figure 3 (a)-(b), illustrating the fidelity and time evolution to create a maximally entangled state by a π/2 and 5π/2 two qubit rotation in the M 1 -coupling manifold, respectively, depicted in the left inset of figure 3 (a). Note that the required time for the latter case exceeds T 2 (Ω=0) by more than a factor of three, therefore impressively demonstrating the decoupling that allows for a 98% final state fidelity with a decoupling field of Ω = 1.2 MHz. In contrast to these simple two qubit examples, recovering the full dipolar interaction form (σ z ⊗ σ z -type coupling, see (2)) is crucial for a wide range of applications ranging from cluster state computation to quantum simulations. Once achieved this allows to create all other types of interactions by merely applying local unitary transformations. The σ i z ⊗ σ j z -type interaction cannot be recovered by any local operation out of the reduced manifold interactions H I,M k ; however one can add up the contributions H I,M1 and H I,M2 in time by using the Trotter or Suzuki-Trotter formalism to implement the total Hamiltonian (see Supplement, section 6): For this scheme to work, it is crucial that the time addition is based on the interaction frame Hamiltonians H I,M k instead of the one defined in (2), because only in that case the different timescales of the decoupling and coupling strength allow to create a high fidelity σ i z ⊗ σ j z gate at the same time preserving the decoherence decoupling effect (see Supplement section 6.1 for more details).
Compensation of systematic errors -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,M k (for the two qubit situation) and H zz may be replaced by J(1 + ij ) with ij describing the systematic error from the optimal case. However, extending the concepts provided in [34] allows to construct compensation sequences (see Supplement, section 7) provided that the error ij 0.5 and that non-nearest neighbour couplings can be efficiently suppressed, as can for example be achieved by approaches like the one discussed in figure 4 (b). 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 S9 (Supplement). 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.2 MHz decoupling field in contrast to the multiqubit counterpart that relies on a more general sequence less efficient in time.
Cluster-state creation and quantum simulation -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 allow to perform any quantum computation operation by purely local measurements on individual qubits [35]. It has been shown [35] 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 |φ C = S| + + · · · + . The connection to the Ising Hamiltonian (3), that can be realized in the nanodiamond system proposed, follows by noting that S exp −i π/4 (i,j) σ i z σ j z up to local operations, the latter being implementable by fast pulses on the electron spin manifold. As can be seen in figure 4 (a) and (b) the combination of decoupling and time addition is capable of creating one and two dimensional cluster states with fidelities well above 90%. As a second application making use of the available interactions, the simulation of a Heisenberg-chain Hamiltonian is illustrated in the Supplement figure S7, an interesting model for the spin dynamics and magnetism in solid state systems.
Summary -In summary we proposed a new method to create scalable arrangements of NV-centers in diamond by exploiting the ability of biological systems for self-assembly along with the precise positioning of surface functionalized nanodiamonds in such structures. We experimentally realized and verified the creation of ordered nanodiamond structures on a protein scaffold, namely on a SP1 monolayer, as well as the SP1-assisted formation of nanodiamond clusters in solution. Based on the achievable NV distances on the nanometer scale we proposed and analyzed theoretically the implementation of single and multiqubit gates and demonstrated its application for the creation of cluster states, thereby addressing the typical problems as the limited coherence time and heterogeneous dipolar coupling strengths. Moderate decoupling fields around 1MHz, well within reach of current experimental setups, allow for the efficient decoupling from surface spin noise with coherence times comparable to the ultimate T 1 limit. Along with significant dipolar couplings of several tens of kHz and the viability of individual addressing, gate fidelities well above 95% can be expected even for multiple qubits and imperfect couplings. We believe that the combination of nanodiamonds with biological systems provides a promising approach towards scalability, overcoming the limitations of current attempts and offering a high level of control in the structure formation.
Acknowledgements: This work was supported by the Alexander von Humboldt Foundation, the BMBF project QuOReP (FK 01BQ1012), the EU STREP project DIAMANT, the ERC Synergy grant BioQ, the GIF project "Non-linear dynamics in ultra-cold trapped ion crystals", the DARPA (QuASAR project), the DFG with the programs SPP 1601/1 "New frontiers in sensitivity for EPR spectroscopy: from biological cells to nanomaterials" and CU 44/3-2 "Single Molecule based Ultra High Density Memory" and the Taiwan Academia Sinica Research Program on Nanoscience and Nanotechnology "Bio-inspired Organic/Inorganic Hybrid Nano Electronic devices". We thank the Peter Brojde Center for Innovative Engineering and Computer Science for support as well as the bwGRiD project for computational resources.
Author contributions: MBP proposed the concept and project; AR and MBP coordinated the project. AA performed the theoretical calculations with input from AR and MBP; FJ provided the nanodiamonds; YN & OS provided the SP1 protein and performed the genetic engineering; GK performed the experiments under the supervision of DP, SY,YP and FJ; All authors discussed the results; AA drafted the paper with guidance and input from AR & MBP and contributions from all the authors. 183-220 (2009 One remarkable aspect of biomolecules is their ability to recognize a wide range of substances with a high degree of specificity. This feature of biomolecules has been extensively explored and a variety of artificial peptide aptamers, that specifically recognize various inorganic materials, have been created by selecting binders from random arrays of amino acids displayed on phages or bacteria (combinatorial biological method) [1][2][3][4][5].
SP1 is a thermally stable protein, originally isolated from poplar trees [6], which self-assembles to an 11 nm ringshape dodecamer (12-mer). The protein is exceptionally stable under extreme conditions, being resistant to proteolysis, high temperatures, organic solvents and high levels of ionic detergent [7,8]. SP1 multivalency allows the display of 12 binding sites upon one protein complex, resulting in a stable and strong binding agent.
Fusion of Carbon Nano Tube (CNT) specific binding peptides to SP1 N-terminus by genetic engineering resulted in an SP1 ring with 12 CNT binding sites, 6 on each side of the ring. This enabled the creation of SP1 variants which tightly bind to CNTs to form a stable SP1/CNT complex [9]. In this work we used the same SP1 variant to attach and order the nanodiamond structures. The carbon sp 2 -hybridization formed on the surface of the nanodiamonds was used to link the nanodiamonds.
Formation of small ordered nanoparticles areas on an SP1 monolayer is performed using the Langmuir-Blodgett method [10]. In this method a trough is filled with a subphase solution that contains the SP1 proteins and by adding glucose to the subphase solution the SP1 floats to the solution-air interface. Then, by slowly reducing the interface surface area, the SP1 are forced to compress to a point where they become a monolayer. The SP1 monolayer is transferred to a substrate and washed with distilled water to remove excess salts from the subphase solution. In figure S1 (a) an AFM scan of a SP1 layer on a silicon chip with a scratched area is shown. By measuring the height difference between the silicon chip's surface (scratched area) to the rest, we verified the existence of a monolayer of height 2nm, which corresponds to the height of the SP1 protein on a Si surface. This shows that at this stage a dense monolayer of SP1 on the surface has been created. In a second stage the substrate with the SP1 array is inserted into a beaker with nanoparticle solution on an orbital shaker, and afterwards it is washed and dried. In addition to nanodiamonds, experiments with gold particles (5nm nanoparticles Sigma Aldrich) have been performed. In that case small ordered areas have been achieved (see figure S1 (b)); however excess salts remained and formed salt crystals that lead to lattice defects in the monolayer. As for nanodiamonds (5nm produced by laser ablation from Ray Techniques Ltd) the excess salts were sufficiently removed by cleaning and we observed small periodic hexagonal arrangements of nanodiamonds for dilute particle solutions; however large scale periodic structures have not been achieved yet (see figure 1 (a)-(b) in the main text).
In order to confirm that the particles (shown in figure 1 of the main text) are indeed nanodiamonds, electron diffraction measurements have been performed on the sample in the relevant areas. The measurements were taken by the Transmission Electron Microscope Tecnai T12 G 2 Spirit. The measured diffraction (1,1,1) lattice parameter results in 2.08Å and the (2,2,0) lattice parameter is 1.26Å. This agrees well with the corresponding diamond lattice parameters known from literature and given by 2.04Å for (1,1,1) and 1.25Å for (2,2,0) [11] (see figure S2).
The formation of the larger (average nanodiamond size 30 nm created by grinding) nanodiamond complexes as shown in figure 1 (e) of the main text, was achieved by mixing 50 µL of 1 mg/mL of the SP1 solution with 50 µL of 1 mg/mL of the 30 nm nanodiamonds solution. To the mixture we added 900 µL of distilled water. 10 µL of the final solution was applied on a silicon chip and scanned by scanning electron microscopy (Extra High Resolution Scanning Electron Microscopy Magellan TM 400L).   The effect of dynamical decoupling or more precise the decoherence decay rate of a decoupled system can be described in terms of a filter function overlap, representing the effect of the decoupling field, with the noise spectrum [12][13][14] (see figure S3). This allows a very illustrative analysis and description of the working principles of dynamical decoupling methods, opens the possibility of measuring bath-coupling spectra by designing appropriate filter functions and allows for the optimization of decoupling sequences. The only restriction arises from the assumption of the weak coupling limit, i.e. by treating the noise influence in a perturbative way. More practically this means that the coherence time must be large compared to the noise bath correlation time, a condition that -dependent on the type of noise-might not be fulfilled for small decoupling fields. However for the parameters considered in this work, in general T 2 > τ , and hence those limitations are not relevant. Moreover the expressions are exact in any limit for the free induction decay and pulsed decoupling schemes provided that the noise is Gaussian as can be easily checked comparing the final expression with the ones derived in a non-perturbative way in e.g. [15]. A summary of the, in many cases somewhat lengthy, formulas will be given in section 2.4.

Decoherence decay rate
Consider a two level system evolving under the Hamiltonian (in the rotating frame) with δ(t) a random, zero-mean fluctuating detuning describing the effect of pure dephasing and originating from environmental coupling and Ω(t) the classical control decoupling field. We will assume that the system is initially (t 0 = 0) prepared in the state whose special cases φ = 0 and φ = π/2 correspond to the σ x (|+ x ) and σ y (|+ y ) eigenstate, respectively. Then, after a time t the probability for still finding the system in that initial state is given by describing purely the effect of decoherence and not taking the coherent evolution into account (see figure S4) (It is important to note that the dephasing term contributes as well to the coherent evolution. For example, in the limit Ω 2 δ 2 the system performs coherent oscillations with a Rabi frequency Ω ef f ω + δ 2 /(2Ω).). The decay rates appearing in (2.3) can be expressed as dτ the noise spectrum and F k t (ω) a decoupling field dependent filter function (see figure S3). Exact expressions can be found in section 2.4. R x (t) corresponds to the decay rate for an initial σ x eigenstate with R x (t) 1/2 S(Ω) for t τ and a constant decoupling field, wherein S(Ω) is the noise spectrum evaluated at the frequency of the decoupling field. The decay rate R y (t) = 1/2 (R x (t) + R γ (t)) corresponds to the decay rate for an initial σ y -eigenstate and is split into two parts. Splitting the decay rate R y (t) into the two contributions R x (t) and R γ (t) allows for a more simple interpretation of the decay behaviour. In the limit of t τ and Ω t 1 the contribution R γ (t) 0 (more precise R γ (t) t R x (t) t) and thus the decay rate of the second decay contribution corresponds to half the decay rate of the first contribution or as a specific example the decay rate of the σ x is twice as large as the decay rate of the σ y eigenstates. Therefore in that limit (t τ , Ω t > 1) and for a constant decoupling field (2.3) takes the form: 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 τ -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 γ (t). Note also that for the case of free induction decay, i.e. Ω = 0, R γ (t) = R x (t) and thus both decay contributions are equal as expected by the isotropy of the system.  S4. Comparison of the decay rate according to (2.3) (red lines) to a numerical simulation of the noise process (blue) for an initial state (2.2) with (a) φ = 0, (b) φ = π/2 and (c) φ = π/4, Ω = 0.5 MHz and modelling the noise as an Ornstein-Uhlenbeck process. Other lines correspond to decay curves for a zero decoupling field and different initial states as indicated in the figure.

Decay rate derivation
Following the description in [12,13] we will derive the general decoherence decay formula (2.3). It is advantageous to evaluate the decoherence behaviour in the σ x eigenbasis |± x ≡ |± 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 [16] allows 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 (2.1) [12] with φ(t − t ) = δ(t) δ(t ) , S(t) = /2σ z and 0 Ω(t )/2 dt σ x 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 (ρ ij = i|ρ|j ), d dt wherein (see definitions in section 2.4, U = exp(i t 0 Ω(τ ) dτ ), R denotes the real part) leading to the solutions On the other hand the differential equation for the coherences takes the form with the additional definitions (herein R denotes the real and I the imaginary part) (2.11) As will be seen later, the combination ρ −+ − ρ +− 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

(2.13)
Now consider the arbitrary phase state (2.2) that can be expressed in the |± basis as (neglecting a global phase factor) leading to the initial density matrix Using eqn. (2.9) and (2.13), the density matrix after a time evolution of t follows to be (2.16) Knowing the density matrix time evolution it is straightforward to calculate the decay rate of the initial state |ψ φ towards the completely mixed state ρ mixed = 1/2 (|+ +| + |− −|) and the probability for finding the system in the initial state after a time t: what corresponds exactly to the result given in eqn. (2.3).
2.3. Population decay for φ = 0 and φ = π/2 and a constant decoupling field In this section the limiting cases of starting in the σ x eigenstate |+ x (φ = 0) and the σ y eigenstate |+ y (φ = π/2), respectively, shall be reviewed for the case of a constant decoupling field Ω(t) = Ω = const., leading to an explanation of the different decay rates in both basis states.
For φ = 0 the initial state is given by (2.15) and the time evolution follows from (2.16) leading to a decay of the initial state (2.17) In the |± basis, that is optimal for the interpretation of the decoupling effect, this corresponds to a population decay with an effective rate determined by the flipping term δ(t) (the pure dephasing in the original state basis) suppressed by the off resonance due to the decoupling field energy splitting Ω. For φ = π/2 the initial state takes the form (2.15) Note that the first contribution already corresponds to the completely mixed state and does indeed not change in time such that the decay is determined by the decay of the σ x basis coherences (2.16) leading to the probability for finding the system in the initial state (2.17) tr (ρ y 0 ρ y (t)) = Referring again to the decoupled decay picture it is clear that the decay of the σ y eigenstates is governed by the decay of the σ x coherences. This interpretation of pure population (φ = 0) and coherence (φ = π/2) decay allows especially for a simple explanation of the decay rate difference in the Markovian limit (t τ , Ω t > 1): In that case the coherence decay induced by the population decay is just half that rate as also follows by noting that in that limit R γ (t) 0 and R x (t) = 1/2 S(0) = R x = const. For a Markovian master equation description that property follows from the conservation of the density matrix trace as can be easily checked analysing the master equation (i.e. the trace of the right hand side should be zero): Specific filter functions:

HAMILTONIAN, PSEUDOSPIN 1/2 DESCRIPTION AND QUANTIZATION AXIS ADJUSTMENT USING AN EXTERNAL MAGNETIC FIELD
Let us first consider a system of NV centers, namely the electron spin-1 ground state triplet manifold ( 3 A) coupled by a dipolar interaction: with the zero-field and external magnetic field contribution [17] Herein D denotes the orientation dependent zero-field splitting tensor, S i the spin-1 operators, B the external magnetic field and γ el the gyromagnetic ratio of the NV center electron spin. In the principal axis system defined by the NV symmetry axis, the zero-field splitting tensor D is diagonal and takes the form: with D = 2.87 GHz and E introducing an additional coupling capable of lifting the | ± 1 state degeneracy. In nanodiamonds the strain dependent E contribution can be nonzero in contrast to the bulk diamond situation [18]. In the principal axis frame (denoted by the primed quantities) H 0 can be rewritten as providing a good choice for either small magnetic fields or in cases where the magnetic field is parallel to the NV center symmetry axis. An arbitrary choice of the reference frame (and therefore the quantization in that case), e.g. choosing a constant laboratory frame for multiple qubits with different symmetry axes, requires to rotate the tensor correspondingly such that it takes a different form for each of the NV's (unless they do have the same orientation). The dipolar coupling term has the form [19] with µ 0 the magnetic permeability and r ij and e ij the distance and unit direction vector between NV centers i and j, respectively. Assuming an equal quantization axis and imposing the secular approximation allows to approximate eqn. (3.5) by with θ ij the angle between e ij and the quantization axis that is given by the external magnetic field direction for B D. Note that the dipolar coupling in the present situation is rather small, e.g. (µ 0 γ 2 el /(4π r 3 ij )) = 2π · 52 kHz for r ij = 10 nm and therefore small inhomogeneous broadening effects (e.g. different strain contributions) as well as different orientations of the NV center symmetry axis as will be outlined below, essentially reduce (3.6) to Creating an arrangement of nanodiamonds has the disadvantage that the symmetry axis of the NV-center, the axis pointing along 'N-V' and forming a natural quantization direction due to its associated crystal-field energy splitting D in that direction, is not controllable, i.e. there exists no common quantization axis in a laboratory frame. This originates from the fact that the symmetry axis of individual NV centers in our setup have an equal probability of S10 pointing in any arbitrary spatial direction. Without or with a small magnetic field γ el B D, the evaluation of the dipolar coupling (3.5) in the principal axis symmetry frames of the centers involved, leads to a coupling that depends as well on the relative orientation of the NV symmetry axes as on the one to the vector connecting the centers, therefore leading to a vast distribution of coupling frequencies inside a multi-qubit array. To solve that disadvantage an applied external magnetic field (γ el B > D), sufficiently strong in the sense that its associated energy shift (γ el B) outperforms the one of the crystal field D, redefines a new and common quantization axis. Considering the limiting case γ el B D (B 1T ) with the magnetic field pointing in the z-direction, allows to rewrite Hamiltonian (3.1) as (neglecting off-resonant coupling terms of the zero-field splitting contribution) with ϑ i the angle between the magnetic field direction z and the symmetry axis of NV center i and H dip given by (3.5)-(3.7) with θ ij the angle of the vector connecting i and j to the z-axis, i.e. the external magnetic field. Two important prerequisites are achieved that way: First the quantization axis is now completely determined by the external magnetic field and not by the symmetry axis any more such that the dipolar coupling will be homogeneous and independent of the individual orientations (for equal θ ij ). Second it provides individual addressability as there exists an orientation dependent (ϑ i ) distribution of transition frequencies differing in the 10 − 100 MHz range. Moreover this property of individual addressing is directly linked to the suppression of exchange flip-flop terms in the dipolar coupling, simplifying the dipolar coupling to the form (3.7). Recalling that the mutual dipolar coupling is of the order of several tens of kHz (∼ 52 kHz for a distance of 10 nm) the much higher differences in the transition frequencies as stated above will make these (in our case non-energy conserving) transitions very unlikely. Thus in the regime of our proposal we can safely describe the coupling as a pure σ z ⊗ σ z -coupling similar to (3.7).
In the more general case of intermediate magnetic fields the effect on the dipolar couplings is more complicated and cannot be simply described in the S x , S y , S z basis manifold. We analyzed this case numerically below for the relevant quasi-particle two-level system. Effectively, recalling that there will be a continuous microwave driving to reduce decoherence, only a reduction to two states of the Spin-1 system is relevant, i.e. S z j = | + 1 +1| + 1/2 (σ z j − 1 j ) or S z j = | − 1 −1| + 1/2 (σ z j + 1 j ) for a zero magnetic field or a magnetic field parallel to the NV center symmetry axis, and therefore the two level reduction allows to rewrite (3.7) as Note that the single flipping terms can be incorporated in the energy Hamiltonian H 0 , will be suppressed anyway when adding a sufficiently strong continuous decoupling driving field or can also be removed exactly by an echo sequence as discussed in section 6. Adding noise and an additional continuous driving for decoupling, the total Hamiltonian in the two level basis can now be written as with b i (t) a fluctuating frequency as will be discussed in section 5. H TLS drive describes the continuous resonant microwave driving for decoupling on the dressed state two level system (3.11) and is given by with Ω i the Rabi frequency of the driving. The dipolar interaction in the secular approximation can be, noting again that only σ z i σ z j terms will in general survive, written as where in the limiting cases of a strong magnetic field, or more general for equal quantization axes of i and j the parameter ξ ij is given by ξ ij = 1/4 (1 − 3 cos 2 θ ij ). For general magnetic fields the parameter ranges are plotted in figure S5 and figure 2 (c) in the main text, showing that around B 0.5 T are required for providing an almost uniform coupling; noting that the variance is even lower there is a high probability that a sufficiently uniform dipolar coupling is already obtained for B 0.2 T . Finally, transforming to a rotating frame with respect to the laser frequency and neglecting counter-rotating terms, (3.10) then reads In the attempt to achieve a uniform dipolar coupling by means of an external magnetic field the orientation of the magnetic field is crucial as can be seen from (3.7) and (3.14). Claiming that the interaction strength should be equal in all spatial directions the maximal dipolar coupling is achieved in a linear chain and in a two dimensional array when the magnetic field is parallel to the chain (ξ ij = −1/2) and orthogonal to the plane (ξ ij = −1/4), respectively. In a three dimensional array this cannot be achieved, because for the only choice that provides equal interactions in all directions, the space diagonal, cos θ ij = 1/ √ 3 and therefore the interaction strength equals zero (performing an addition in time with non-equal couplings in all spatial directions can however circumvent that problem). As a last remark it should be noted that the combination of a strong magnetic field and the condition cos θ ij = 1/ √ 3 could be used to decouple the system from the dipolar interaction, e.g. after a cluster state is achieved and further local operations and measurements are desired subsequently.

NV − -center hyperfine-structure and its influence on the decoupled gate interaction
An additional complication, that has been neglected so far, 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). In the following we will just consider this 'intrinsic' nitrogen nuclear spin in the center itself and assume that other nuclear spins (e.g. I = 1/2 of C-13) only couple weakly to the center and thus can be treated within the framework of the spin-bath decoherence as outlined in the next section and in the main text. This can be considered as a good approximation for a high purity diamond as well as for strong decoupling fields such that the hyperfine coupling is weak compared to the microwave field, i.e. other nuclear spin are sufficiently far-apart from the NV-center. In that framework the ground-state level structure in the principal axis frame (analogue to (3.4)) can be described as [20,21] H gs = H 0 + H nucl + H hyp (3.16) with H 0 the pure electron spin-Hamiltonian given by (3.4) and H nucl , H hyp the (nitrogen) nuclear spin part and the hyperfine coupling interaction, respectively. The nuclear spin Hamiltonian in the presence of a magnetic field B takes the form with the quadrupole splitting P = −4.95 MHz and the gyromagnetic ratio γ I 3.1 MHz/T. For the hyperfine coupling we will neglect exchange ('flip-flop')-terms due to the large energy mismatch between the electron and nuclear spin transitions, such that For the regime of our proposal (B=0.5 T), this leads to a hyperfine-structure level scheme as depicted in figure S6. With the selection rule ∆m I = 0 it becomes clear that neighbouring transitions between two different electron spin states differ by ∆ω ∼ 2.2 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,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 Ω |A |. 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. In that limit, the NV-center is described by the electron-spin states to a good approximation and the rather small differences in the hyperfine transitions of ∼ 2MHz allow to perform such a strong driving by still keeping the microwave power small enough to enable individual addressing (the latter one determined by the rather strong zero-field splitting ∝ 2.8 GHz). We would like to point out that continuous microwave drivings of ∼ 40 MHz have already been successfully implemented in experiments [22]. (ii) A second option would be the regime of weak microwave driving Ω |A |. In that regime only a single nuclear spin state transition is excited, provided that the microwave is tuned to resonance with a single specific transition. Off-resonant excitations can be neglected in such a regime and one ends up again with a perfect two-level system. However for the noise parameters considered (see figure 2 and 3 in the main text), the decoupling in this regime does not lead to optimal results even it might be sufficient for a short π/2-pulse interaction with tolerable fidelity. As a remark, one might also explore the intermediate regime Ω |A | by identifying conditions where the effect of the off-resonant contributions leads to an effective 2π-multiple rotation similar to the technique performed e.g. in [23]. (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. Such polarization schemes have been demonstrated in numerous experiments to date as well at room temperature [21,24,25] (making use of an excited state level anticrossing that allows direct hyperfine exchange interactions) as at low temperatures using projective measurments [26] or the concept of coherent population trapping [27].

NOISE SPECTRUM FOR A FLUORINE-TERMINATED NANODIAMOND SURFACE
The main source of decoherence for the nitrogen-vacancy center in diamond arises from couplings of the associated electron spin to the spin bath of additional impurities [23]. For the case of bulk-diamonds, it is well-known that the dominant source of decoherence can be attributed to C-13 nuclear spins (high purity type IIa diamond) and paramagnetic P1-centers associated with nitrogen donors (effective electron spin-1/2 system; type Ib diamond), that naturally occur within the diamond crystal. For nanodiamonds the situation becomes quite different as a result of the surface proximity. In that case surface spins, associated either with the carbon hybridization or the terminating elements, form a dense spin bath that couples with significant strength to the electron-spin of the NV-center [18]. We will therefore focus on this surface induced dephasing mechanism in the upcoming analysis. However it should be noted that other decoherence processes, in particular the mechanisms known already from the bulk diamond counterparts, can and indeed have been described in the same framework [28,29] presented in the next two sections. Thus these mechanisms will be decoupled along with the surface-spin decoupling as successfully demonstrated in numerous decoupling experiments to date [29][30][31], and, keeping in mind that such mechanisms form much smaller The nuclear spin consists of the N-14 spin associated with the NV center. The angle ϑ between the external field and the NV-center symmetry axis is assumed to be sufficiently small in this illustration, allowing to neglect couplings between different nuclear spin states.
spin-bathes, we do not expect changes in the reported T 2 -times. Here we would also like to point out that even a small number of much stronger coupled electron spins can be decoupled in that framework, as shown at the end of this section and in figure S8. Terminating the nanodiamond surface by fluorine, oxygen or hydrogen / hydroxyl groups replaces the electron spins associated with the sp2 hybridized orbitals by much weaker nuclear spins [32] and therefore reduces the dipole-dipole interaction strength between surface spins by a factor of 10 −6 (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 −3 . 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 neighbour distance 2.5Å [33]. Such a coupled system can be described by with H NV 0 the NV center energy Hamiltonian as given in (3.2), H sf 0 the magnetic field splitting of the surface nuclear spins, H sf int the dipole dipole coupling between surface nuclear spins and H hf the hyperfine coupling of the surface spins to the central spin (the NV center electron spin): 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 and r i the one 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 section 3 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 section 3, one can substitute S z → 1/2 (s z ± 1) with s z the spin 1/2 operator of the quasi-spin. Herein we just retained secular contributions of the dipole dipole coupling and neglected flipping terms in the hyperfine interaction (that would be related to T 1) due to the large energy mismatch between surface spins and the NV electron spins. However, nevertheless those flip flop processes are very unlikely, second order processes might lead to a significant enhancement of the surface spin flips described by H sf int of the order of A 2 i /∆ with ∆ being the typical energy splitting of the NV levels and A i the hyperfine coupling amplitude [34][35][36]. However we do not expect those terms to be significant in the present proposals that requires magnetic fields B 0.5 T what leads to mediated couplings of the order of A 2 i /Ω ∼ 1 Hz and has to be compared to the direct nuclear nuclear dipole coupling, that, in the dense surface arrangement (2.5Å) leads to next neighbour flipping rates of 6.8 kHz and is comparable to the mediated one for a nuclear spin distance of 5 nm. Thus, in the dense bath at large magnetic fields considered here, the decoherence effect arises mainly from the flip-flop interaction of close nuclear spins which is to a good approximation not affected by second order hyperfine processes whereas the decoherence influence of the weakly interacting remote spins, affected by the mediated contribution, can be neglected compared to the first part (despite the fact that it leads to larger differences in the field).  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 case of dipolar couplings to a spin bath [37]. Since then various approaches, mean-field and exact approaches in certain limits, have been invented to infer the noise properties and calculating the decoherence decay rate, ranging from the simple Ornstein-Uhlenbeck model [28] to cluster expansion methods [34,36,38,39]. In here we used the mean-field method described in [15,40] what 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 [41]. In the pair-correlation approximation one assumes that each of the flipping processes in the dipolar coupling (4.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 eqn. (4.2)) with [H ij , H kl ] = 0 ∀ij, kl . 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 − exp(−q N 2 flip /N ) is small [36] (with N flip the number of flipped bath spins during the considered time, q the number of nearest neighbours 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 what is taken into account by using a mean-field type approach based on the theory of moments. From eqn.
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 1 [28]. Following the calculation presented in [15] this leads to (neglecting static contributions that were also neglected in the decoherence discussion in section 2.1).
and σ ij the mean field broadening used in replacing the original delta-peaks by Gaussian peaks leading to the same second moment as the one obtained by the moment theory The spectrum (4.6) was numerically evaluated by randomly distributing nuclear surface spins on a sphere with equal The dephasing due to the nuclear surface spins (blue) is modelled by randomly placing nuclear spins on the nanodiamond surface with the diamond assumed to be a sphere of radius r = 5 nm. The NV center is placed in the center of that sphere with the symmetry axis aligned with the external magnetic field that determines the quantization axis of the surface spins.
inter-spin distance and placing the NV center in the center of the sphere, i.e. for a sphere radius r = 5 nm a number of 5026 (what corresponds to a distance of 2.5Å) nuclear spins was distributed on the surface (see figure S7). Assuming an equal quantization direction for as well 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 (5.3). Depending on the number of surface spins an average over different random positions is performed (what however does not lead to a significant change for several thousand spins considered here). We performed the same calculation as well for different numbers of electron surface spins (see figure S8) on diamond with radius 5 nm. In that case the decoherence parameter are very poor, e.g. T 2 = 0.1µs, τ = 21.7 ns and b = 3.6 MHz for 30 surface spins and getting worse with further increasing the number of spins. For such a situation decoupling fields in the Ω ∼ 1τ ∼ 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. Of course it is questionable to what extend the mean-field spectrum calculation, initially based on a S16 pair-approximation, retains its validity in the fast flipping regime of electron spins. To gain some insight in the electron spin case we performed some exact numerical T 2 calculations for small electron spin numbers (note that for such small numbers the results obtained by (4.6) are not justified because the spectrum obtained that way differs significantly from a Lorentzian and is very sensitive to the position of the surface spins): For 10 surface spins T 2 is already below 0.5 µs such that the correlations times obtained by the spectrum approach, even they might be not exact, seem realistic at least in the order of magnitude. Therefore it is crucial to avoid the presence of surface electron spins for performing reliable quantum gates.

DEPHASING NOISE SIMULATION AND THE ORNSTEIN-UHLENBECK PROCESS
Dephasing of the nitrogen vacancy center spin states is mainly caused by long range dipolar interactions with a spin-bath, consisting e.g. of substitutional nitrogen atoms (P1 centers) and 13 C nuclear spins [29,31] or in the case of nanodiamonds predominately of unpaired surface spins [18,32,42]. Herein a significant energy mismatch of the transition energies prevents flip-flop processes between the vacancy center and the bath spins limiting the influence to a pure dephasing process. However intra-bath flipping processes are not suppressed and occur on the noise correlation timescale τ such that the NV center is influenced by a random configuration of the bath environment what can be modelled in the mean field approximation as a random magnetic field leading to the frequency shift b(t) and therefore Now assuming that the back action of the nitrogen spin on the large bath is small and that there is no coherence between different bath spin evolutions, one can model b(t) as a random Gaussian, Markovian and stationary process what is also known as an Ornstein-Uhlenbeck process [29] and can be specified by the two time steady state correlation with b a measure of the coupling strength and τ the noise correlation time, or alternatively by the expected Lorentzian noise spectrum [43] For such a process the evolution of the random field is governed by the Langevin equation [44] db(t) dt with Γ(t) a Gaussian zero-mean noise [45] ( Γ(t) = 0, Γ(t)Γ(t ) = δ(t−t )). The parameter c in (5.4) is related to the definitions in (5.2) by b 2 = c τ /2 what follows by solving b(t) b(0) by means of (5.4). Exact updating formulas can be obtained for the differential equation (5.4) and also for the time integral of b(t) given by y(t + dt) = y(t) + b(t) dt [44]: and n 1 , n 2 statistically independent unit normal random numbers.
Simulation of the noise process: Simulating the evolution governed by a Hamiltonian of the form (3.15) can be easily performed by numerically solving the Schrödinger equation and using the exact updating formulas (5.5). In here we used a Suzuki-Trotter decomposition splitting the Hamiltonian in the noise H n and remaining part H r , such that the time evolution is given by ( = 1) and the noise part can easily be evaluated using (5.5).

CREATING AN EFFECTIVE σz σz TYPE COUPLING
A full Ising interaction, that is a σ z σ z type coupling in the presence of decoherence is a crucial task towards the way to a universal set of gates, that, once achieved, easily allows to create all other types of interactions by just applying local unitary transformations, straightforward to implement by local microwave pulses on the electron spin transition, e.g. : σ x i σ x j = U x σ z i σ z j U † x or σ y i σ y j = U y σ z i σ z j U † y with U x = exp i π/4 (σ i y + σ j y ) and U y = exp i π/4 (σ i x + σ j x ) . Noting that (see also definitions in the main text)  Oscillating decoupling field in the filter spectrum description for t = 25 µs. The red curve illustrates the noise spectrum (T2 = 13.3 µs, τ = 2.5 µs, b 2 = 30.2 kHz) and the dashed black lines the filter functions for zero decoupling and Ω0 = 0.5 MHz, respectively. The blue, pink and orange lines correspond to filter functions for an oscillating decoupling Ω(t) = Ω0 cos(ωosc t) with ωosc = 1 kHz, ωosc = 300 kHz and ωosc = 1 MHz, respectively, showing that the decoupling effect vanishes for high oscillation frequencies comparable to the noise correlation time.

Heisenberg Chain simulation
As a second application beside the cluster state creation discussed in the main text, we illustrate the possibility of simulating a Heisenberg-chain Hamiltonian of the form H = −1/2 N j=1,µ J µ σ j µ σ j+1 µ with µ = {x, y, z} and J y = J z = J, J x = δ J. This task essentially requires two steps: Creating the individual parts by using the methods described in the previous discussion, namely the σ j x σ j+1 x -contribution and the trivial σ j y σ j+1 contributions and in a subsequent step adding those two contribution in time. The accuracy of such a procedure is illustrated in figure S10 for a θ = π/2 evolution that is limited by non-nearest neighbour interaction as well as the number of Trotter cycles (2-3) such that each time slice dt still provides Ω dt > 2π to guarantee the desired H M k Hamiltonian form.
References [46,47] provide a method to remove the systematic error contribution based on noting that multiples of 2 π pulses with respect to the ideal error-less gate end up with pure contributions which can, together with the property wherein σ φ = cos φ σ x + sin φ σ y , be used to create a compensation cycle by means of a Suzuki-Trotter time addition of the two contributions in (7.3). Restricting this method to the lowest order of the Suzuki-Trotter expansion (in order to keep the pulse sequence simple and additionally because higher orders require more time and therefore end up in a loss of additional fidelity by decoherence), it follows that for 2 π cos φ = −θ/2 and σ φ = cos φ σ x + sin φ σ y .
Instead of giving a detailed derivation of the compensation sequence here we refer to reference [46] and will discuss the similar idea for the multiparticle case in 7.2 as well as the special case of the M φ (π) M 3φ (2π)M φ (π) sequence S22 in section 7.3 with an emphasis on the complication with increasing particle number. Note that the change of the rotation axis φ can be accomplished by using e −i θ/2 σ φ = e −i φ/2 σz e −i θ/2 σx e i φ/2 σz (7 .7) what allows to rewrite the rotation M φ as The application of (7.4) to the evolution of Hamiltonian (7.2) is straightforward by just replacing σ φ , σ x , σ y , σ z appearing in (7.5-7.8 More precise it follows that wherein the last two options have additional contributions that however do not affect the states in the gate manifolds. In summary, a compensated gate is obtained using the sequence (see figure S11) z + O( 3 ij ) . (7.20) Taking into account that the error scales as O(n 3 ij ) the parameter n should be as small as possible, leading to the obvious choice n = 1. It is straightforward to extend the compensation analysis to the complete evolution (7.14) what finally leads to and N c the space of non-neighbouring qubits. As an example the method is illustrated for the creation of a four qubit cluster state in figure S12 (a), showing that high decoupling fields are required in presence of decoherence in order to benefit from the compensation despite the significant increased gate time that is problematic in terms of decoherence.

S25
Now let us return to the case of multiple particles. Herein the same concepts can be applied as in section 7.2 despite now using the pulse sequence (7.22). However it turns out that in such a situation the condition analogue to (7.28) with T φ and N c defined as below eqn. (7.21), is not fulfilled. That this is in general not valid can best be seen by considering the example of three particles in a linear chain in which case what reveals that due to the double commutation with two σ φ contributions there is no net sign change in the phase when commuting exp −i φ σ 2 x in contrast to the situation for a single particle in (7.28). That prevents the application of relation (7.7) and hence σ φ terms cannot be cast in σ −φ as in (7.27). One should note that this issue can still be fixed for the three particle case by choosing T φ in (7.29) such that the phase change is performed on the first and third particle in what case the left and right hand side of that equation would be equal thus allowing to use a compensation cycle of the form (7.23). For more than three particles such a possibility however does not exist. To conclude: Due to an even number of permutations in (7.29) (i.g. two permutations for a linear chain, four for a two dimensional array and six for a three dimensional one) it is in general not possible to extend the n = 1/2 compensation sequence (7.23) to the multiparticle case (see figure S12(b)). Therefore the general concept (7.21) has to be used in that situation, despite being less efficient in both time and compensation range.

Non-next nearest neighbour couplings in the compensation sequence
The Hamiltonians (7.1) and (7.13) used to describe the compensation mechanism have been defined up to nextnearest neighbours, i.e. higher order contributions were neglected so far. However, whereas those contributions give only a small correction for simple π/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 long pulse times of such a sequence higher order contributions cannot be neglected. One has to distinguish between two types of higher order couplings: Even couplings defined as next-nearest neighbour and third, fifth, . . . nearest neighbour couplings and odd couplings as second, fourth, . . . nearest and diagonal couplings. Whereas the first ones have always the form of a σ z σ z coupling and are automatically compensated for by the next-nearest neighbour based compensation sequence, the second class consists just of the reduced manifold M k couplings and is not removed by the compensation sequence. This failure for odd couplings is as well originating from the fact that the compensation pulse acts on both qubits involved such that (7.17) is not fulfilled any more, as on the fact that the manifold interaction does not commute with the σ z σ z type contributions such that a splitting of the evolution in a perfect and defective part as in (7.14) 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 S13 for second nearest neighbour couplings. With the evolutions obtained that way as a starting point higher orders can be removed in a concatenated way taking benefit of the different evolution timescales dependent on the qubit distance. For the multiqubit compensation sequence (7.20) 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 orders effects are too weak for giving significant contributions. S26 + = + + FIG. S13. Removing 2nd nearest neighbour interactions for a linear four qubit configuration. Each of the blocks removes 2nd nearest neighbour couplings; however both blocks are needed to restore the proper next nearest neighbour interaction. Blue lines denote σzσz and red lines M1 interactions and dashed lines denote a negative sign of the corresponding interaction. Red marked qubits denote that the interaction is embedded between a U and U † pulse with U = exp(−iπ/2σx). Couplings higher than second order are neglected in the illustration.