Introduction to Quantum Optimal Control for Quantum Sensing with Nitrogen-Vacancy Centers in Diamond

Diamond based quantum technology is a fast emerging field with both scientific and technological importance. With the growing knowledge and experience concerning diamond based quantum systems, comes an increased demand for performance. Quantum optimal control (QOC) provides a direct solution to a number of existing challenges as well as a basis for proposed future applications. Together with a swift review of QOC strategies, quantum sensing and other relevant quantum technology applications of nitrogen-vacancy (NV) centers in diamond, we give the necessary background to summarize recent advancements in the field of QOC assisted quantum applications with NV centers in diamond.


B. The Rotating Wave Approximation 24
References 24 The nitrogen-vacancy (NV) center 1,2 is one of the major platforms in the evolving field of quantum technologies. Its remarkable stability, long spin coherence time, and optical properties make it especially attractive for quantum applications. Like any other quantum system, however, it is subject to experimental imperfections and limitations. Quantum optimal control 3 (QOC) aims to improve the efficiency of system manipulation under such constraints. In this introductory review article, we aim to provide an overview of the NV center applications that can be improved with QOC and outline the modus operandi for their implementation. This review has been organised as follows: Section I introduces NV centers in diamond and their most useful properties. Section II reviews the basic quantum sensing techniques, as well as a brief account of possible quantum information and computation applications of NV centers. Section III breaks down the principles and methods of QOC theory. This includes an introduction on the structure of QOC problems and an overview of some of the most common numerical QOC algorithms. Section IV reviews the control techniques discussed in section III applied to the system and methods discussed in section I and II, respectively. This review has been written such that the reader can refer to section I and section III more or less independently. Finally, the Appendices give a description of the relevant Hamiltonians for the NV center spin system.

Examples
Throughout this review, illustrating examples are given in colored text boxes to distinguish them from the general text. Diamonds host a variety of point defects 5,6 i.e. lattice sites where one of the carbon atoms is replaced with a different atom or vacancy. The appearance of color in pure diamond is due to the presence of particular optically active point defects. Ib diamonds, for example, are yellow (see Fig. 2a) due to single substitutional nitrogen impurities. Such fluorescent point defects have a unique spectral signature and are called color-centers. Among them are the silicon-vacancy center 7 , the germanium-vacancy center 8 , the tin-vacancy center 9 , the NV center, and several others which have been extensively studied over the past two decades. 5,10 It is noteworthy that the color centers in diamond emit photoluminescence (PL), which is bright enough to enable observation of single defect sites with a confocal microscope, e.g. in regular patterns of nanostructures (see Fig. 1, and caption for details). To add the capability for spin manipulation, microwaves are applied to the color centers under investigation.  12 The subject of this review, the NV center, is formed when one of the four carbon atoms in the unit cell of the diamond lattice is replaced by a nitrogen atom, accompanied by the formation of a neighboring vacancy site, shown in Fig. 1. NV centers have three known energetically stable states; NV − , NV 0 , and NV + . In this review, we limit our focus to NV − (from here on simply called NV). It is the most promising of the three for applications in quantum metrology and quantum information because of its magneto-optical properties. In the following sections, we briefly describe these useful characteristics of the NV center for the matter of clarity and selfsufficiency with a main focus on optimal control applications. For more details, we refer the reader to more in-depth reviews on NV centers and related applications 2,13,14 .

A. General Properties of NV Centers in Diamond
As a host material, diamond itself has a number of useful properties, including hardness, high refractive index, and non-toxicity. Various techniques can be used to manufacture NV center containing diamond (Fig. 2). The main fabrication methods include high pressure-high temperature (HPHT) synthesis, detonation synthesis and chemical vapor deposition 15,16 (CVD). While detonation synthesis only leads to nanodiamonds (see Fig. 2b). HPHT and CVD produce diamonds of various morphologies including bulk, single crystals. Very recently, up-scaling single crystal CVD diamond to wafer scales (10 cm) has been accomplished via heteroepitaxy 12 (see Fig. 2d). Ion implantation 17 is commonly used to create NV centers in diamond, other methods include doping the crystal during the growth process to incorporate the desired impurity. The availability of a broad range of materials fosters a spectrum of applications including scanning probe diamond tips (Fig. 2c) with single NV centers for nanoscale sensing, [18][19][20][21][22][23] engineered NV ensemble based sensors, 24 wide field imaging with NV center ensembles in bulk diamond, 25 diamonds coupled to cavities [26][27][28][29] as well as biological applications [30][31][32][33] in different environments, temperature ranges and pressures. Despite the many advances in material synthesis techniques, diamonds have inherent noise sources that limit the performance of NV center based applications. In diamonds with exotic geometries like the scanning probes ( Fig. 2) as well as systems with shallow NV centers, even the external environment plays a big role. Accordingly the protection from both inherent and external noise has a high priority for the NV center community to improve performance. Exploring techniques of QOC (section III) promises to make the aforementioned spin systems more robust against unwanted noise. Also important for its applicability is the NV center's electronic structure [34][35][36] (Fig. 3): The vacancy site contains six electrons (three dangling carbon electrons, two nitrogen electrons, and one electron trapped from a nearby donor in the crystal lattice as shown in Fig. 1) which form a spin triplet system with m s = ±1, 0. The electrons are tightly bound to the defect making the NV center a truly atomic sized systems. FIG. 3. NV center electronic structure. The zero-phonon line (ZPL) is represented by red arrows (λ = 637 nm), green arrows signify off-resonant excitation, and the dashed arrows indicate the possible non-radiative decay channels. The spin state lifetimes are given as t LT , and the in-box level scheme shows the Zeeman splitting of the ground-levels under the influence of a static external magnetic field B nv .
As depicted in Fig. 3, the excited and ground state spin triplets are separated by 1.945 eV, corresponding to a purely radiative transition (zero-phonon line [ZPL] in Fig. 4a) at 637 nm. Additionally, intermediate singlet states occur energetically in-between the triplet states. At thermal equilibrium, the spin population in the ground states is governed by the Maxwell-Boltzmann distribution law. 37 Upon irradiation with a green laser (λ = 532 nm), the population in the ground states is transferred to the excited levels through a combination of radiative absorption and non-radiative relaxation processes, involving the excited levels of the NV center and the conduction band of the diamond crystal. The excitation process is followed by direct radiative decay as well as non-radiative decay through the intermediate states. The radiative transitions between the ground and excited states are spin preserving, the decay via the intermediate states however, is not (although it is spin dependent, nonetheless). The decay rate from the excited level m s = ±1 spin states to the intermediate singlet state is comparable to the rate of the direct radiative decay to the ground state. For the excited level m s = 0 spin state, the decay via the singlet state is negligible. This results in a majority of the population being pumped into the m s = 0 ground-level spin state under laser illumination. The transition frequency in absence of an external field between the degenerate m s = ±1 and the m s = 0 spin states in the ground-level is 2.871 GHz. An external magnetic field lifts the degeneracy of the m s = ±1 states as depicted in Fig. 3 and experimentally observed in Fig. 4c (Zeeman splitting). If the magnetic field is sufficiently strong the m s = ±1 states are well separated in frequency allowing for the spin system to be approximated as a qubit, with one of the m s = ±1 states and the m s = 0 state forming the two qubit levels.

Rabi Oscillations
Let us consider an example (also see the example in section III A), approximating the NV center as a two level system with the ground state |0 and excited state |1 . This assumption holds when the m s = ±1 spin states are non-degenerate under the influence of a static magnetic field along the NV axis (the axis connecting the nitrogen atom and the vacancy in the lattice). It is common practice to consider the NV axis as the quantization axis (Fig. 1). The system can be initialised in the m s = 0 ground state with a green laser pulse. Subsequently, a MW field of frequency ω nv is applied to resonantly drive the transition |0 ↔ |1 , leading to a continuous population transfer. The observed coherent oscillation (Rabi flopping) of the population between the spin states is shown in the figure below on the left hand side.
The figure on the right shows the Bloch sphere representation of the NV qubit spin state. The north pole of the Bloch sphere corresponds to the pure |1 state, whereas the south pole indicates the pure |0 state. The Bloch vector (red arrow) represents the spin-state of the system at a given instance,x andŷ indicate the polarisation of the driving field in the rotating frame (see Appendix B). Pulses rotate the state vector with a speed dependent on the driving field strength around an axis dependent on its phase (blue curve). A rotation by 180 • , in this case fully transferring the population from m s = 0 to m s = 1 is referred to as a π-pulse. π/2pulses rotate the state-vector by 90 • . In this case, they would create a superposition state between m s = 0 and m s = 1. The Hamiltonian for this system in the rotating wave approximation (see Appendix B) is given by Eq. (6). If the drive is constantly applied, the system undergoes Rabi oscillations with a frequency Ω that is proportional to the strength of the applied microwave field. The decay of the amplitude of oscillation is owed to the inherent decoherence and decay processes in the spin system.
Following points summarise the properties that render NV centers well suited for quantum technology applications: • The NV center has a magnetically active ground-level spin triplet state with the m s = ±1 states exhibiting Zeeman splitting (Fig. 3) in the presence of external fields. This sensitivity towards external fields is the key to the majority of sensing applications.
• In addition, the transition between the ground-level spin states can be coherently manipulated with microwaves.
• The ground-level spin states also exhibit spin dependent fluorescence due to selection rules and the transition rates mentioned above. It ensures a clear readout contrast between the two qubit levels, which is crucial for quantum applications (see Fig. 4d). Alternatively, other readout schemes based on NV charge state detection are also possible.
• The m s = 0 spin state of the ground-level triplet can be efficiently initialised via an off-resonant laser excitation, this qubit initialisation is a pivotal property for all quantum sensing and quantum computation applications.
• The nuclear spins in the NV center cluster exhibit dynamic nuclear polarisation, 38 opening up the possibil-ity of nuclear spin based quantum computation applications.
• All these properties are exploitable at room as well as cryogenic temperature.
• Additionally, diamond is a bio-compatible host, which is crucial for life science applications.
Various other properties (Fig. 4), e.g. stable emission of single photons, also make these spin systems highly suitable for quantum sensing and quantum information applications.

B. Spin Hamiltonian
Before we explore the vast field of applications of the NV center, it is useful to study the governing Hamiltonian for the spin system in order to understand the theoretical basis of the applications. The NV center spin in the diamond lattice undergoes a variety of interactions because of external magnetic and electric fields, as well as the crystal strain field, and other spins in the crystal lattice. The spin Hamiltonian for the NV center's ground-level spin triplet can be obtained in the basis of the m s = ±1 and m s = 0 spin operators. It can be derived by perturbative expansion of the full Hamiltonian of the system in terms of the spin operators (see Appendix A for more details), such that the thermal degrees of freedom average out. 39,40 Assuming the NV axis to be the quantization axis, the spin Hamiltonian can be written aŝ The first term in the Hamiltonian above is the zero-field term, i.e. the system Hamiltonian in the absence of external fields, whereh is the reduced Planck constant,ˆ S = (Ŝ X ,Ŝ Y ,Ŝ Z ) are the spin operators (see Appendix A), D is the axial and E the non-axial zero-field parameter. At room temperature and ambient pressure D ≈ 2.87 GHz varies by around −80 kHz/K and 1.5kHz/bar with temperature and pressure changes respectively (see sections II A 3 and II A 4). The second term in the Hamiltonian represents the magnetic field interaction, where γ nv = 2π × 28 MHz/mT is the gyromagnetic ratio for the NV center spin and B is the external magnetic field. The electric field interaction with the NV center spin is represented by the third term in the Hamiltonian.
is the effective electric field vector, representing contributions from external electric fields as well as the crystal field, arising from crystal strain and pressure (see section II A 4). δ and δ ⊥ are the axial and transverse coupling constants respectively arising from the symmetry of the crystal structure at the site of the NV center. In comparison to the magnetic field, the electric field has a much weaker interaction with the NV spin, with δ at around 0.17Hz/(Vm −1 ) and δ ⊥ of the order of 10 −3 Hz/(Vm −1 ). The NV center spin interacts with other spin systems in the host crystal and its surroundings. The most dominant interactions are with the 15 N or 14 N nucleus in the NV center and the 13 C nuclear spins in the surrounding diamond lattice. These spin-spin interactions are represented by the hyperfine interaction terms in the Hamiltonian. N i is the hyperfine interaction tensor between the electron and the i th nuclear spin andˆ I i is the nuclear spin operator. There are two main contributions to the hyperfine interaction tensor: First, the isotropic Fermi contact interaction that arises from the interaction of the electron cloud with the nearby nucleus (see Appendix A). Second, the anisotropic magnetic dipole interactions of the NV spin with distant nuclear spins. The magnitude of the latter interaction decays as 1/r 3 with the distance r from the nuclear spin and is comparatively weaker than the Fermi contact interaction terms. On top of the hyperfine coupling, the nuclear spins also interact with the external magnetic field. This is represented as the nuclear Zeeman interaction term, where γ i is the gyromagnetic ratio of the corresponding nuclear spin i. Finally, nuclei like 14 N can exhibit quadrupole splitting; this is represented by the last term in the Hamiltonian. Here Q i represents the quadrupole splitting. As an example, the quadrupole splitting for the 14 N nuclei inside the diamond lattice has been estimated to be around 5 MHz 41 . The reader is advised to refer to Appendix A for a detailed description of these Hamiltonian terms and their origins.

II. QUANTUM TECHNOLOGY APPLICATIONS WITH NV CENTERS
We now proceed with the exploration of the versatility of NV centers in terms of their potential applications in quantum technologies. As quantum sensors, they have been utilised for magnetic and electric field sensing (section II A 1 and II A 2), thermometery (section II A 3), strain/pressure measurements (section II A 4), orientation tracking and more. NV centers have further found application in quantum information related experiments (section II B). In addition, the NV center is also a reliable non-classical source of single photons, and has been utilised in various single photon experiments. 42,43

A. Quantum Sensing Applications
In principle, NV centers possess the necessary properties to be practical quantum sensors. 44 The efficient spin state initialisation/readout and the sensitivity to various physical parameters have resulted in a variety of diamond based sensing applications. As diamonds are simple to handle in terms of logistics, maintenance, and manipulation in comparison to other (atomic) quantum systems, their application has received a lot of attention. The main focus has been on using the quantum nature of the defect center's spin to detect classical physical signals (called classical detection, in contrast to quantum detection which is based on the use of entanglement 44 ). In several sensing applications NV centers have proven to outperform their counterparts. Fig. 5 shows one such example, juxtaposing different scanning magnetometers that are commonly implemented for sensing. A summary of the different measurable parameters with NV center-based sensors is provided in Table I. The sensing methods rely on efficient and coherent spin manipulation, whereby experimental imperfections limit their success. In this section, we briefly discuss different sensing protocols, experimental limitations and other challenges. A review of the current techniques based on quantum optimal control (QOC) (see section III) to overcome/circumvent the described obstacles and efficiently use the NV centers as quantum sensors is given in section IV A. For completeness, we also outline related applications in quantum information and computation.

Magnetic Field Sensing
Magnetic field sensing is the most common application of NV centers as quantum sensors. As illustrated in the full Hamiltonian Eq. (1), the external magnetic field interacts directly with the NV spin. In addition to continuous wave spectroscopy experiments (Fig. 4c), several NMR/EPR based spin manipulation techniques 50 with microwave pulses can be applied to quantify the magnitude of the interaction as well as the directional dependence of the external field to be studied. These techniques pave the way for vector magnetometery, 19 magnetic field imaging at nanoscale, 20  spins, 51 and a number of other applications. All these methods rely on efficient initialisation (laser), manipulation (microwaves) and readout (typically laser, although several other possibilities exist 52,53 , see Barry et al. 13 for details) of the spin system. The most basic sensing protocol is the direct measurement of the Zeeman splitting (see section I B) under the influence of an external magnetic field (Fig. 4c). This technique requires the application of a continuous, off-resonant laser to the spin system. The MW frequency is swept, driving the m s = 0 → m s = ±1 transition when the resonance condition is met. This in turn results in a measurable decrease in the fluorescence from the NV center (see Fig. 4c). Gruber et al. 54 demonstrated the first of these optically detected magnetic resonance (ODMR) experiments with a single NV center, showing that the splitting of resonant peaks (dips) gives information about the external magnetic field. Non-polarized (thus randomly flipping) spins proximal to the NV center create fluctuating magnetic fields at the position of the NV center. The interaction strength of the NV with the fields fluctuating at the transition frequency between m s = 0 → m s = ±1 can be directly inferred from the lifetime (T 1 ) of the spin states. Such relaxometry techniques have in particular been used for life science applications, 25 magnetic noise sensing, 55 and surface and material studies 56 .
Another common method for DC magnetometery is the Ramsey interference pulse technique 57 depicted in Fig. 6. Balasubramanian et al. 58 first experimentally demonstrated this method with NV centers using ultrapure CVD grown single crystal diamond as a host material. In theory, the Ramsey sequence is the spin equivalent of an optical Mach-Zehnder FIG. 6. Ramsey spin manipulation sequence for DC magnetic field sensing. The NV axis is considered to be the quantization axis. The parallel component B , of the external DC magnetic field is the signal under consideration. Resonant MW control pulses are used to apply π 2 -rotations. The first pulse rotates the spin from the ground state into a superposition in the equatorial plane (blue curve). During the free precession time τ the magnetic field component B induces a phase φ (black curve). When the spin is rotated again with a second pulse, the phase can be inferred by projecting the spin along the quantisation axis. The required measurement is done by repeating the process a number of times and reading out whether the NV center is in the ground or excited state after each repetition. The distribution of measurements then gives the projection along the quantisation axis.
interferometer. It is based on measuring the phase induced by external perturbations on a superposition state of the spin system. The "mirrors" and "beam splitters" for this interference experiments are formed by the π 2 -pulses applied at MW frequency (Fig. 6). The induced phase φ is directly proportional to the external magnetic field component B (t): where τ is the free spin precession duration. The phase accumulation arises from the Larmor precession of the spin about the magnetic field vector. Additional pulses can be used to cancel out the phase induced by the DC field, so that the phase information can be used for measuring AC magnetic fields (Fig. 7). This idea is exploited in the spin-echo pulse sequence 59 (Fig. 7). The additional π-pulse refocuses the signal such that the total phase accumulation in this case is When the time period of the AC field is an even multiple of the free precession time τ, the two terms for phase accumulation above are equal and cancel each other (see Fig. 7b). This way the phase induced by slowly varying fields over the entire sequence is zero. Such spin manipulation experiments with NV centers have been independently reported by several authors in the last decade. 58,[60][61][62] The spin-echo protocol with a single refocusing π-pulse is the basic building block of the dynamical decoupling (DD) sequences 63 (Fig. 8) that make use of multiple refocusing πpulses. Decoupling the system from DC or slowly varying  Fig. 6, at the bottom, the red vector represents the state of the system. The first pulse rotates the spin by 90 • from the ground state into a superposition in the equatorial plane (blue curve). During the free precession time τ the parallel magnetic field component B induces a phase φ (black curve). The π-pulse then flips the direction of precession, proceeded by another phase accumulation over a time period of equal length. Finally, a π 2 -pulse rotates the spin again and the magnetic field strength can be inferred by projecting the spin along the quantisation axis. In case of a DC field the phase contributions are of opposite sign for the two free precession periods and hence are subtracted. However, for fields oscillating with frequency ω ac = π/τ, the final projection gives a measurement of its amplitude similar to the Ramsey sequence in Fig. 6. (b) Population in the state |0 after applying spin-echo sequences with different precession times τ. Whenever the frequency of the external field satisfies the condition ω ac = kπ/τ with odd k, the fidelity is reduced. The overall envelope decay of the population depends on the decoherence time T 2 . fields is accompanied by coupling to a given frequency signal (Fig. 7). The general phase accumulation for a series of π pulses over a duration of time t can be estimated as where M is known as the modulation function and reflects whether the phase accumulated over the specific period τ in the sequence is positive or negative (Fig. 8).
In the case of a broadband external field (periodically or FIG. 8. Dynamical decoupling sequence: (a) The protocol consists of a train of π-pulses. The frequency that is measured depends on the spacing between the pulses. The modulation function reflects the state of the NV center at every point of the sequence. By multiplying the signal with the modulation function, one gets the rectified signal, whose integral is proportional to the phase accumulation. (b) The Fourier transform of the modulation function quantifies the weighing function for the detectable frequencies using the given sequence. Please note that the π-pulses are assumed to be instantaneous such that the phase variation during the pulse is negligible (hard pulses).
The peaks (k) in the plot correspond to the different harmonic orders of the filter function. (adapted and modified with permission from Degen et al. 44 ) randomly oscillating), the contribution to the total phase induced by different frequency components depends on the socalled weighing function ( Fig. 8b) for the given pulse sequence. These weighing functions are similar to those of narrow-band frequency filters, i.e. they enhance signals from resonant frequencies and suppress non-resonant external field contributions. 64 Decoupling the spin system from the noisy spin bath in the diamond crystal enhances the coherence time of the spin states. 65 It also improves the sensitivity towards the targeted signal which is limited by the spin dephasing (T * 2 ) and spin decoherence (T 2 ) times. For more information on the nature, usefulness and application of these DD-sequences and filter functions for quantum sensing, the reader may refer to the extensive review on the topic by Degen et al. 44 . QOC (section III) offers an alternative approach to enhance the coherence time for systems with high frequency noise such as shallow NV centers. In addition, the hard pulse approach (i.e. using rectangular pulses) is subjected to cumulative errors with respect to the phase and duration of the pulses, that become significant for longer DD sequences. QOC provides an elegant solution to avoid such errors, which is useful for sensing applications as well as quantum information and computation applications (section II B) that require high fidelity gate operations.
As an indirect application of magnetic field sensing, gyroscopic measurements with NV centers have been performed. 66,67

Electric Field Sensing
The coupling of the magnetic field to the NV spin is much stronger than that of the electric field. Hence, electric field sensing techniques 68 are experimentally more challenging. However, the external electric field directly influences the spin, as discernible from the Hamiltonian in Eq. (1). FIG. 9. Electric field sensing with NV centers: (a) Careful alignment of the external magnetic field B (orthogonal to theẑ-and hence to the NV-axis) is required to measure the effect of the electric field. E ext represents the external electric field and E str the non-axial crystal strain field. E e f f indicates the resulting effective electric field experienced by the spin. (b) The plots show the Stark effect shifting the resonant frequencies between the NV center spin states for different electric field strengths. For details on the experiment see Dolde et al. 69 Comparing their coupling, the magnetic field's effect (represented by γ nv ) is of the order of a few GHz/T 19 , whereas the analogous constant of the electric field E I is of sub-Hz/Vm −1 -order 68 , making the coupling very weak, as indicated in section I B. Nevertheless, electric fields are also detected quantitatively by the phase induced during the pulse sequence. Hence, similar spin manipulation protocols employed for magnetic field sensing are used for detecting electric fields. As the signal is much weaker than for magnetic fields, the challenge lies in separating it from the effects of other strong undesirable interactions. Dolde et al. 69 exploited the interplay between magnetic Zeeman effect, Stark effect and crystal strain field in their carefully devised experiment to suppress the magnetic field effects on the NV spin, making the shifts in the electric field transitions more prominently detectable via spin manipulation protocols ( Fig. 9).

Thermometery
Unlike the electric and magnetic fields, temperature is not directly coupled to the spin system (strictly speaking spinphonon coupling is still temperature dependent, however temperature as a parameter does not directly couple to the spin). Instead, temperature fluctuations have a direct effect on the crystal field as they affect the lattice constant of the crystal. This in turn influences the zero-field splitting parameters discussed in section I B. Acosta et al. 70 reported the dependence of the parameters D and E on temperature, demonstrating that temperatures changes can be directly quantified from ODMR spectra (Fig. 10). They found that the axial zero-field splitting parameter is more sensitive to temperature changes than the non-axial parameter. Neumann et al. 75 utilised the D-Ramsey spin manipulation protocol to suppress the effect of external magnetic fields, paving the way to an enhanced measurement of temperature changes. I. Summary of quantum sensing with NV centers. η B and η E are the sensitivities with respect to the magnetic field B and the electric field E, respectively. σ is the standard deviation in the measurement (related to the detected noise), S is the detected signal, and t m is the measurement time (limited by the coherence time of the system). D and E are the axial and non-axial zero-field parameters, respectively.

Pressure and Strain Sensing
Similar to the temperature, other physical quantities like pressure and strain affect the crystal field. Changes in these quantities can be studied by their coupling to the zero-field splitting parameter.
Doherty et al. 74 reported the effect of hydrostatic pressure on the ZPL and the ground state ODMR of NV centers at room temperature in a hydrostatic pressure medium (Fig. 11). Crystal strain variations lead to changes in the effective electric field as described in section I B. The reader is referred to Maze et al. 19 for a theoretical description of the effect of strain on the NV center's electronic structure. Teissier et al. 76 reported NV spin-strain coupling in diamond cantilever devices.

B. Quantum Information and Computation Applications
The same characteristics that make NV centers a promising sensor, also qualify them in the fields of quantum information and computation. 77 As described in section I A, the NV electron spin interacts with the internal N nuclear spin of the NV center and 13 C nuclear spins in its close proximity (NV center cluster). Abobeih et al. 78 showed that even a much more precise characterisation is possible: They imaged a 27-nuclear spin cluster using a single NV center. Considering the surrounding nuclei, the resulting hyperfine structure 79 (see section I B and Appendix A) forms a more complex quantum mechanical systems than a single electronic qubit. In fact, the majority of quantum computation applications can in principle, be performed with a nuclear spin qubit, using the NV center electron as a quantum bus qubit for initialisation and readout. 41,80 This technique has lead to various NV-based applications such as quantum error correction, 81-83 quantum algorithms, 84,85 quantum simulation, 86,87 fault-tolerant quantum repeaters for long-distance communication, 88 and quantum memory. [89][90][91] Entanglement between optical photons and the NV spin has also been demonstrated experimentally 92 . In rare cases two NV centers might be close enough to each other for dipoledipole interaction 73 . This provides another possible source of entanglement and is further described in Appendix A.
In the proceeding paragraphs, we describe applications mainly regarding the use of NV centers as quantum registers, entangling gates with NV centers, and use of NV center qubits to perform quantum error correction. Quantum registers are at the heart of quantum computation techniques (Fig. 12). They contain the coupled set of qubits which are used to perform quantum computations. The nuclear spins in the NV center cluster exhibit long coherence times at room temperature, which is an essential feature for spin system based quantum registers. Recently, Bradely et al. 93 presented a quantum register based on the electron spin of an NV center together with nine nuclear spins depicted in Fig. 12. To efficiently implement gates in such a large system, they simultaneously used dynamical decoupling based gates on the NV center electron spin and selective phase-controlled driving of the nuclear spins. As a result they saw coherence times from ten seconds up to a minute at room temperature. In a different approach, Neumann et al. 94 demonstrated the working concept of a diamond-based spin register containing two optically addressable NV center spin qubits coupled to each other in an ultrapure and isotopically engineered CVD grown diamond. NV center based spin registers have also been used for quantum simulation of molecules. 86 Similar to the way quantum registers form the underlying structure, entanglement generation enables the connections required for any quantum information technology (entanglement also plays a vital role in several quantum sensing applications, for more information the reader is advised to refer to Degen et al. 44 ). For example, Dolde et al. 95 demonstrated NV-NV spin entanglement at room temperature, while Neumann et al. 96 demonstrated entangled two-and three-particle quantum states with 13 C nuclear spins in an NV center based quantum register. Another essential concept for fault-tolerant quantum computation is quantum error correction (QEC). As mentioned above, NV-based spin registers may provide systems with long coherence times and desired entanglement properties. These are prerequisites to implement standard QEC protocols. Taminiau et al. 82 demonstrated error correction with a spin register of nuclear spins from the NV center cluster, using the NV center electron spin as the control qubit. In another work, Unden et al. 97 demonstrated enhanced quantum sensitivity with metrology techniques based on iterative quantum error correction. Casanova et al. 98 performed a complete set of universal quantum gates with an NV center-based spin register, providing a proof of principle for NV center based quantum computation. To obtain better control over the qubit system, several QOC based methods have been implemented so far. They helped to overcome decoherence effects caused by the environment, and to circumvent experimental limitations. We will review some of these works in section IV B. Furthermore, NV centers have been long known as a non classical photon source 42,99 , which was applied in single and two photon interference experiments. 100,101 The NV center has also been used for Bell's inequality test related experiments. 91

III. OPTIMAL CONTROL THEORY
Without the ability to precisely manipulate quantum systems, researching their properties and applying them for quantum technologies is almost impossible. Quantum Optimal Control (QOC) theory 3,102 improves the shape of dynamical controls (typically electromagnetic field pulses) to achieve a certain goal to maximum precision. The section starts with the details of defining a QOC problem in section III A followed by a description of different numerical optimisation tools in section III B and concluded by a brief discussion of the limits of QOC (section III C). The first part is structured according to the schematic in Fig. 13. In the field of Nuclear Magnetic Resonance (NMR), pulse shaping is used since the 1980s 103 and many of the arguments for pulse shape optimisation 104 equally apply to NV centers, which we will focus on. In many cases, the time scales defining the decay of NV centers are large compared to the control time. In that case, it is sufficient to study closed system dynamics. Indeed, specific open system techniques such as population suppression and the exploitation of useful dissipation processes are often not applicable to the problems considered in this review. 105 Hence, we limit ourselves to a closed system FIG. 13. Schematic drawing of a generic QOC optimisation. The box on the left contains the elements that define a basic QOC problem with blue solid arrows connecting them to the algorithm. The grey box at the center illustrates the optimisation algorithm itself, with the dotted grey arrow indicating its iterative nature. The cost function J is calculated from the controls u i (t) and used to update the controls. In parenthesis the relevant sections in the review paper are indicated where applicable.

A. Defining a Control Problem
The principles of QOC theory derive from early extremisation problems such as Johann Bernoulli's brachistochrone curve problem. 106 Similarly, QOC problems are formulated through a system of equations, which broadly defines three things: First, the system dynamics, i.e. the theoretically obtained description reflecting the system's behaviour, for example given by the Hamiltonian. Alternatively, this first equation might be replaced by a description through the experiment itself. Second and third, the control objectives and control space restrictions. The objectives on the one hand set the goal of the optimisation, like e.g. high fidelity for the transfer to a target state. The control space restrictions on the other hand limit the resources that may be used to reach the desired goal. Together, the three aspects are combined into a so-called control landscape. Each set of controls will result in a different value of the "cost function" J, a measure for how close the system is to reaching the objective. In case of a minimisation (and throughout this review we will always assume minimisations, unless stated otherwise), each valley corresponds to a locally optimal combination of controls. The goal of the optimisation can now be easily defined as reaching the lowest point in the landscape. We will now discuss in more detail these three ingredients of QOC problems as well as the initial guess, stopping criteria and robustness.

System Dynamics
One way to characterise the evolution of a closed quantum system with time dependent controls is through Schrödinger's equation. The system Hamiltonian is usually split into two parts; the drift HamiltonianĤ d , which is constant and cannot be manipulated, and the control HamiltoniansĤ c i which are multiplied with time-dependent coefficients u i (t) called "control pulses". The full Hamiltonian then readŝ Please note that system dynamics for control problems may also be defined through Lindblad-operators and even for non-Markovian dynamics (for a review on open systems QOC see Koch 105 ).

Drift and Control Hamiltonian
As an example, let us consider a NV center, approximated as a qubit with the ground state |0 and excited state |1 . In this simple consideration the goal will be to create a high-fidelity π 2 x -rotation similar to the system in Frank et al. 107 A static magnetic field B is applied in z-direction (the quantisation-/NV-axis) and a circularly polarised microwave field B ⊥ with an amplitude B ⊥ (t), frequency ω mw and phase ϕ is applied orthogonal to B . Let us define the gyromagnetic ratio of the NV center as γ nv . The rotating frame of B ⊥ then gives the Hamiltonian where ∆ = ω nv − ω mw is the detuning, ω nv = B γ nv the NV's resonant frequency, Ω(t) = B ⊥ γ nv is the Rabi frequency andˆ s are the spin operators in the |0 , |1 basis. A derivation of this Hamiltonian can be found in Appendix B. We can easily identify the drift HamiltonianĤ d = ∆ŝ Z . Let us assume, that both the Rabi frequency Ω(t) and the phase of the magnetic field ϕ(t) can be manipu-lated dynamically. The control Hamiltonians may then be identified asĤ c 1 =ŝ X andĤ c 2 =ŝ Y and the control pulses as u 1 (t) = Ω(t) cos ϕ(t) and u 2 (t) = Ω(t) sin ϕ(t).
Once the system has evolved, it is time to test whether the goals have been reached by checking the control objectives. It should be noted at this point that the rotating wave approximation, as presented in Appendix B, is widely used to simplify the NV center's Hamiltonian. While it is useful, when the Rabi frequency is much lower than the NV center's resonant frequency, it can have a detrimental effect on a simulation's accuracy, if the Rabi frequency is of a similar scale as the qubit transition. In fact, Scheuer et al. 108 have shown how the inaccurate use of the RWA can affect the outcomes of optimal control procedures designed for NV centers.
FIG. 14. Example of a QOC landscape. Considering two control parameters, we may represent the cost function J as a surface dependent on the set of controls. The minima correspond to locally optimal control coordinates. Each black path represents a local optimisation starting from a different initial guess.

Control Objective(s)
The cost function (or figure of merit) J defines what is minimised in any QOC problem. This way it describes the goal of the optimisation (terminal cost) and optionally the control limits through penalty-terms (running costs). The terminal costs are determined at the final time of the system's evolution. They quantify for instance the distance between the final state and the desired goal state in the relevant Hilbert space. The running costs are usually related to the restrictions on the control pulses, for example the limited power of a microwave source. In this review, the cost function is defined to be zero, when all objectives are met and to be greater than zero, when they are not met. Note that the running costs have a similar role as the control space restrictions that will be discussed in section III A 3. In the following, we will briefly describe some of the most relevant cost functions found in relation to NV centers in the literature. For more examples of NV center applications, please refer to section II and to section IV for examples specifically combining them with QOC.
• State to state transfer (terminal cost) State to state transfer is the most common optimisation objective and has been used in many papers 107,[109][110][111] . The infidelity is a measure for the distance between two states |φ t and |φ(T ) : If they are equal, it gives zero, if they are orthogonal, it has a value of one. The infidelity can be used directly to define the cost function J state describing the distance between |φ(T ) , the final state of the system at time T , and |φ t , the target state. An alternative way to formulate the transfer is fixing the global phase using J state = 1 − Re{ φ(T )|φ t }.
• Unitary gate optimisation (terminal cost) To measure the distance between the unitary U(T ) produced by the controls and the target gate U t , we define the cost function In the second line, the gate fidelity is defined through the N 0 basis states |ζ i of the initial system and their propagated version |φ i (T ) = U(T ) |ζ i . Similarly to J state , we can also define a global phase dependent version of this cost function, J gate = 1 − 1 N 0 Re{Tr(U † t U(T ))}. Examples for its application can be found in references [112][113][114][115] .
• Sensitivity (terminal cost) In contrast to the previous examples, the sensitivity does not directly contain information about the system. Instead, it quantifies the amount of information about a parameter θ (e.g. the magnetic field) that may be derived from a set of measurements.
The cost function should, however, not contain information about the outcome of the measurement, but rather about its precision. The lower bound of ∆θ is given by the Cramér-Rao bound a value which is inversely proportional to the Fisher information F(θ), calculated by where the second line is specifically related to a discrete number of possible measurement outcomes N x .
One may interpret the Fisher information as the curvature of the logarithmic probability distribution: If it is completely flat, hence giving no information, F(θ 0 ) = 0, if it is strongly peaked, indicating a clear parameter estimate, F(θ 0 ) 0. The corresponding cost function may be defined as Reviews introducing Fisher information in the context of quantum sensing and metrology were written by Degen et al. 44 and Pezze et al. 117 The original paper relating Fisher information and quantum mechanics was published in 1994 by Braunstein and Caves. 118 Applications of Fisher information as a part of optimal control can be found in references. 119,120 • Limited power (running cost) The power of a control pulse is typically calculated as P i = T 0 |u i (t)| 2 dt. To limit P i to a reasonable range P i ∈ [0, P lim ] a penalty term can be introduced which adds a high cost to J, if a certain limit is crossed.
The function κ(P i ) should give very little to no penalty, if the power is within the acceptable range κ(P i P lim ) → 0 and a high penalty, if it is out of range κ(P i P lim ) → ∞. These criteria can be satisfied by a wide variety of functions and it depends on the chosen system. Examples can be found in a number of references. 104,[121][122][123] • Limited bandwidth (running cost) There is a number of ways to limit the bandwidth of the controls. One solution is to gently punish any quickly oscillating solutions through where is some small factor 124 . It should be noted that this expression, does not give strict bounds in terms of bandwidth. An alternative, stricter approach is to punish fast oscillations, only if they lie outside a predefined filter function as described by Schäfer et al. 125 and Kosloff's group. 126,127 A completely different approach is to restrict the basis of the control pulses. This is possible with certain algorithms of the (d)CRAB family including GROUP, and GOAT and will be further discussed in section III B.
There are many more possible terminal costs, each describing a different control problem including partial state transfer, taking into account the full density function, maximising entanglement, 128 or adjusting a certain observable. 129 Similarly, equally many different running costs exist e.g. to avoid populating fast decaying states. 130 Gate Optimisation In the experiment by Frank et al. 107 the control objective was to optimise a unitary defined as the Hadamard gate Hence the cost function may be defined as We can see that this is a good cost function as it is minimal when U(T ) = U t at the final time T (up to a global phase). If we were to also include a bandwidth limitation on the two control pulses u 1 (t) and u 2 (t), we may simply sum up different cost terms. The resulting cost function, where i are some small factors, would be Running costs favour acceptable types of controls, as opposed to physically impossible ones, but if stricter limits are required, control space restrictions might be the more suitable mean of limitation.

Control Space Restrictions
While the running costs (see section III A 2) can only passively punish controls which lie outside the achievable frame, control space restrictions actively change the controls to only allow what is experimentally achievable. One might imagine them as a horizontal squeezing and stretching of the control landscape or as the introduction of hard walls (see Fig. 14), opposed to a vertical distortion induced by running costs.

Restricting the Control Amplitude
As an example, let us consider the amplitude of a control pulse u i (t) that should be restricted to u max i . The pulse could be cut off at the beginning of each iteration according tõ This form ensures maximum exploitation of the amplitude space but is not differentiable, hence it requires that the control pulse is cut off during an extra step (Fig. 13) before the cost function and/or gradient is evaluated. 131 An alternative approach is mapping the control pulse to a restricted subspace using a continuous function. For example by replacing it withũ i (t) = u max i sin(u i (t)). 132 Another common example for the application of mapping are shape functions. They restrict the overall shape of the pulse, which is useful, if e.g. the experiment requires a smoothly rising and falling control pulse with Γ(0) = Γ(T ) = 0, such that u i (t) = Γ(t)u i (t).

Initial Guess and Stopping Criterion
Numerical optimal control techniques are based on iterative algorithms, which require a starting point (called "initial guess") and a clearly defined situation to stop at i.e. the stopping criterion. The optimisation will in most cases find the the closest local minimum to the initial guess (see examples in Fig. 14). Accordingly, it is often helpful to try a number of initial guesses to find which one is closest to the global minimum. The stopping criterion is simpler to define: It might be based on the maximum number of iterations (limited computation time or experimental run time), a measure for convergence or the clear definition of a goal.

Robustness
Usually, there is some discrepancy between the theoretical model and the experiment. In an optimisation, this can be taken into account to ensure that the optimised pulses will work in the presence of such a discrepancy by averaging over cost functions for slightly different systems. Let us consider each system being described by the HamiltonianĤ i , then by taking into account N rob different versions, the cost function becomes

Robustness Against Detuning
The resonance line of the NV center has a finite width and can be described by the normalised distribution f (ω). Off-center NV centers can however still be described with the HamiltonianĤ in Eq. (6) by adjusting the static magnetic field B || and hence the detuning ∆.
One may average the cost over N det different detunings ∆ i to get a robust cost function We can see that this cost function will only reach zero, if all J(Ĥ(∆ i )) are zero, ensuring robustness against the detuning. By including the probability distribution f (B || γ nv ), we ensure that the optimisation favours solutions centered on the average detuning. Due to field inhomogeneities in B ⊥ , the Rabi frequency can also have a finite distribution when considering an ensemble of NV centers. In many optimisations, both detuning and Rabi errors are accounted for simultaneously. 122,133

B. Numerical QOC Algorithms
Once the problem has been defined, an algorithm is required to systematically test possible solutions minimising the cost. 3 In this review, we will only describe numerical optimisation algorithms as they have produced promising results and a variety of packages exist to implement them (see section III B 5 for more details). There are, however, alternative strategies, such as geometrical optimal control 134,135 (GOC) and shortcuts to adiabaticity 136,137 (STA). They usually rely on a deeper analytical analysis of the control problem and hence access a smaller solution space than QOC, but can nevertheless be effective. One example is the direct application of Pontryagins minimum principle (PMP) which falls under the category of GOC. It has been shown to provide time optimal evolution for NV centers. 138,139 Similarly, STA has been used to implement specific gates on NV centers 140 and protect them from decoherence. 141 Section III A started by mentioning the brachistochrone problem, whose solution is usually obtained via an analytical variational approach. In the case of quantum mechanical problems however, one would produce a set of nonlinear equations which, in most cases, cannot be solved analytically. Instead, a variational approach combined with numerical solving was introduced by Konnov and Krotov,142 and Sklarz and Tannor 143 and further adapted by Ohtsuki et al. 144 Since then, different attempts have been made to numerically solve this class of problems. In general, two families of QOC algorithms can be identified: Gradient-free 145 and gradient-based 104 . Gradient-based algorithms determine the derivatives of the cost function with respect to the control pulses to find an improved solution. These methods are usually efficient as they make use of all the available information. Gradient-free methods on the other hand can be applied directly to experiments or to complicated problems, where the gradients are not straightforwardly calculated. In this review, the direct experimental implementation is referred to as closed-loop, while a purely simulation based optimisation is called open-loop. We will start by looking at the working principle of gradient-based optimisation algorithms, before exploring their gradient-free counterparts.

Gradient-based Optimisation
To understand how gradient-based algorithms work, let us first consider the effect of a small change ∆u in some control u(t) on the cost function J(u(t)) (see section III A 2). If the change ∆u is small enough, we can approximate We can now deduce the properties ∆u should have to decrease J. Indeed it enables us to make small changes and update u(t) iteratively. Consider ∆u = − ∂J ∂u(t) , where is a small positive factor. In this case, the new cost function becomes which is smaller than the previous value, implying an optimisation. In order to avoid functional derivatives and iteratively improve the cost, the control function u(t) needs to be split up into time independent control parameters u (k) . According to the simple updating algorithm above, the new control parameters u (k) become More advanced updating algorithms promise faster convergence. Eq. (18) could for example be extended to second-order. 142,146 A popular method approximating the second-order term from the first-order term is the L-BFGS (limited-memory Broyden-Fletcher-Goldfarb-Shanno) quasi-Newton method algorithm. 147,148 In the following, we describe the working principles of the GRadient Ascent Pulse Engineering (GRAPE) algorithm as an example for the whole class of algorithms. It was originally designed for Nuclear Magnetic Resonance (NMR) 104 but has since found various applications with NV centers. For another comprehensive explanation, we refer to Saywell

GRAPE Optimisation of a State Transfer
We start by defining an exemplary cost function J, where J = 0 implies the target state φ t is reached at a time T (also see section III A 2): The initial state is defined as |φ(0) = |φ 0 andT is the time ordering operator. In order to take the derivatives of J, we choose the piece-wise constant control basis, chopping up the control pulses into N small slices u (k) i of width ∆t. This gives a new way to formulate the propagator U(T ): It should be noted that this basis is not the only possible choice but intrinsic to GRAPE (as well as to other gradient based algorithms like Krotov 150 ). We can now reformulate the cost function as We start by calculating the derivatives of φ t |T k U (k) |φ 0 w.r.t. the control parameters.
We have defined the forward propagated state φ (k) and the backward propagated state ξ (k) , which is usually referred to as the adjoint state. They can both be easily calculated by solving Schrödinger's equation. A graphical representation is given in Fig. 15. We applied the chain rule to find the gradient of J. In the case that the control pulses are mapped to a restricted subspace (see section III A 3), the chain rule can be used again. The last thing left to evaluate are the following directional derivatives: All in all, any gradient-based optimisation algorithm relies on calculating the first derivative of the cost function with respect to the control parameters. On top of this specific example of GRAPE, we list below the most commonly used gradient-based QOC algorithms and their natural features (for an illustration of the terms "sequential" and "concurrent" see Fig. 16): • GRAPE 104,148 GRadient Ascent Pulse Engineering concurrently optimises in the piece-wise constant basis.
• Krotov 115,142,143,146,150,151 Krotov's method sequentially optimises one control parameter after the other. It also relies on the piece-wise constant basis.
• GROUP 152 GRadient Optimization Using Parametrisation is based on GRAPE combined with the chain rule and optimises concurrently. Its chopped basis is flexible (also see "CRAB" in section III B 3) but relies on an initial piecewise constant basis.
• GOAT 132 Gradient Optimisation of Analytic conTrols is based on a system of equations of motion obtained by differentiating the full propagator with respect to the control parameters. The parameters are optimised concurrently and its chopped basis is flexible (also see "CRAB" in section III B 3).

Gradient-free Optimisation
In an experiment, the gradients described above cannot be calculated analytically. Certain finite-difference methods help to find them regardless. [154][155][156] However, if the control landscape is not smooth this method might prove inefficient or very costly in terms of measurements. This is where gradientfree optimisation algorithms shine. Even for certain openloop optimisations they can offer an alternative, when their gradient-based counterparts fail: If e.g. the dynamics of a system are significantly more complicated than described in section III B 1, the gradient of the cost function might be hard or impossible to find analytically. One example for such a case is the CRAB algorithm, described below, which was initially introduced to optimise many-body problems using tensor networks to simulate the time dynamics. 145 The first step then is to choose an optimisation basis. In the following, we will focus on the CRAB algorithm (see section III B 3) and consider the basis of trigonometric functions but it should be noted that we could also use Slepians, Chebyshev polynomials or indeed piece-wise constant elements. Broadly following Caneva et al. 121 the expanded control pulses each take the form Each pulse is composed of a sum of N be basis elements. Each basis element is defined by a frequency ω and the control parameters [A n , B n ]. The index n stands for the iteration number. If the number of available basis elements is restricted (i.e. only a certain region of frequency space is accessible), this is called a chopped basis (CB). Especially in the case of bandwidth limitations, only optimising in the accessible restricted control space can be a powerful tool to avoid introducing distorting penalty terms (see section III A 2). Decreasing the number of parameters, also shrinks the size of the search space, potentially making the optimisation a lot more efficient.

CRAB
FIG. 17. Illustration of the chopped random basis. A frequency space is segmented into N be parts. In each part , a frequency ω is randomly selected according to Eq. (29). Two corresponding parameters, A and B , are optimised. They are defined as in Eq. (28) and represented in this plot by black and blue crosses.
The Chopped RAndom Basis algorithm (CRAB) 121,145 is defined by the optimisation of a random choice of basis elements taken from a truncated function space. Intuitively, one might instead chose the basis elements to coincide with the principal harmonics of the pulse. However, Caneva et al. 121 showed that randomness can be surprisingly effective, especially if the energy scales of the system are not fully known. Indeed, a larger function space is covered, if multiple optimizations are done with different randomized bases. The elements that make up the basis of our example in Eq. (28) are defined by the frequencies ω and are illustrated in Fig. 17. These frequencies are chosen according to where ω max is the maximum admissible frequency and = {1, 2, ..., N be } is an index which selects the chunk of the frequency space that ω is chosen from. Let us further choose the random numbers r from an interval [−0.5, 0.5]. Then the bandwidth of the control pulses is automatically limited to [0, ω max ], where a typical choice is ω max = 2πN be /T (T refers to the length of the pulse). By changing the ω max , we can change the bandwidth. Moreover, we can see that the available frequency space has been split into N be regions permitting the optimisation to make use of the entire space. It also conditions the optimisation problem to have clearly distinct control parameters. It should be noted that the number of basis elements should be dependent on the number of degrees of freedom inherent to the system. 144,[157][158][159] During the optimisation the 2N be -dimensional landscape will be followed using any updating algorithm (it could even be gradient-based as in GOAT and GROUP). The most common choice is the gradient-free Nelder-Mead algorithm 160

dCRAB
In the basic version of CRAB, the basis elements are fixed and the local control landscape is explored for all N be frequencies simultaneously. This leads to a restriction in the number of frequencies that can efficiently be optimised. Using the dressed Chopped RAndom Basis algorithm (dCRAB), much fewer basis elements with ω d, need to be optimised at a time (N be (dCRAB) < N be (CRAB)). Instead, when one CRAB routine converges, we move on to ω d+1, . This enables the method to include an arbitrarily large number of bases and to derive the solutions without -whenever no other constraints are present -being trapped by local optima. The extra iterations changing up the basis after each CRAB-run are called superiterations and the index d refers to the d th superiteration. Their effect is illustrated in Fig. 18. If their number is fixed to N SI , the full description of the pulse can be summed up at the end of the optimisation u opt , with all optimised parameters A opt d, , B opt d, as In each superiteration only the parameters with corresponding index d are optimized. By repeatedly changing the basis, dCRAB does not get caught in local minima for most control problems and thus allows to retain this advantageous property of unconstrained control algorithms in a parametrized (e.g. bandwidth limited) setting. Rach et al. 157 explored the improvement from CRAB to dCRAB in detail considering the random Ising model. They found that convergence may be achieved by taking enough parameters to fix the degrees of freedom present in the optimisation problem. The underlying algorithm used for both CRAB and dCRAB performed best for a basis with 10 − 20 parameters. This allowed dCRAB to outperform CRAB as it requires less optimisation parameters per optimisation (i.e. superiteration). All in all, dCRAB, promises faster convergence with respect to CRAB as fewer parameters are optimised in parallel and instead, new basis elements are chosen sequentially. An example for its experimental application to NV centers, among others, can be found in the work of Frank et al. 107 where a Hadamard gate was optimised.

Optimal Control Packages
In the past years, a number of QOC algorithms were implemented in ready-to-use software packages. In this section, we present four of these packages that we deem to be closest to applications with NV centers. An overview over some of their distinguishing features is given in Table II. Nevertheless, more solutions exist. RedCRAB 162,163 is a python based programme, aiming to remotely optimise any experiment or simulation with gradientfree methods. It can be linked to the experiment setup via MATLAB, python, terminal or simple file transfer and is hence very versatile. RedCRAB makes use of the dCRAB alogithm and provides pulse updates. As it does not require any knowledge about the quantum system itself, it is compatible even with more complicated many-body systems and tensor-network simulations. RedCRAB is available from the authors on request. DYNAMO 153 was originally developed as a GRAPE (and Krotov) implementation in MATLAB. It allows the user to choose their own Hamiltonian and dissipator terms as well as one of the available figures of merit. Hence, it combines simulation and optimisation for certain problems dealing with small quantum systems. It allows for the optimisation of robust pulses and includes a large number of examples. The full version is available on github. QuTiP 164,165 is an open source python library for simulating quantum systems. One of its features is a quantum optimal control implementation. As such it offers limited optimisation techniques with GRAPE and CRAB. Conveniently, the optimisation settings are defined with the usual QuTiP structure. The library is available for example via pip or conda. 167 The Krotov package 166 is an open source python library built on top of QuTiP. As such it offers optimisation via Krotov's method. It includes an extended range of settings in comparison to QuTiP's own QOC implementation. The library is available for example via pip and conda. Other QOC packages include Spinach 168 and SIMPSON, 169 which focus on NMR applications, as well as QEngine 170,171 which includes a GROUP implementation designed especially for ultra-cold atom physics. GRAPE was also recently implemented in the GRAPE-Tensorflow python package, 172 using methods known from machine learning to calculate the gradients.

C. Limits of Control: Controllability and the QSL
Whether or not a QOC problem is (approximately) solvable, is not always simple to answer. However, by examining a number of characteristics of the Hamiltonian, some general predictions can be made. First of all, one may ask whether the control objective is in principle reachable. This can be addressed by examining the controllability of the system. 173 The drift and control Hamiltonians define a certain state (and also gate) space that is reachable. A system is called controllable when all states (gates) in the Hilbert space are accessible in finite time. It has been shown that, if the rank of the dynamical Lie algebra generated by the different terms of the Hamiltonian corresponds to the rank of the control space (and fulfills certain symmetry criteria), the system is fully controllable. Alternatively, the question of state-controllability can be examined via a geometric approach based on graph theory, which can be more convenient to check, especially for larger systems. 174 For more information on controllability, please refer to the following books. 173,175 For open quantum systems, the deleterious effect of the environment usually can not be completely canceled and only a subset of the whole set of states (gates) can  105 If the controllability criteria are fulfilled, the question remains whether the controls are complex and energetic enough to navigate the Hilbert space to the specified target. In general, the quantum speed limit (QSL), i.e. the smallest possible control time needed for a system to reach its target, is influenced by two factors. First, the dynamical equation determines how fast the system may change. This is usually quantified by the so-called Schatten p-norm of the dynamical operator (see reference 176 for details). Second, the exact distance between the initial system to the objective needs to be taken into account.
Quantum Speed Limit The minimum time it takes to evolve a system into a target state is mostly dependent on two things: The HamiltonianĤ and the distance between the initial and the target state φ 0 |φ t . For a time-independent Hamiltonian, we obtain the Bhattacharyya-bound: 177,178 This time T QSL is called the quantum speed limit (QSL). It can be interpreted as follows: If the Hamiltonian has a high energy variance calculated on the initial state ∆E = φ 0 |Ĥ 2 |φ 0 − φ 0 |Ĥ |φ 0 2 , then any other state is reached more quickly. It might be more intuitive to consider the case that φ 0 is an eigenstate of H, hence it will never change and as ∆E = 0, the speed limit will go towards infinity. The distance to the target state finally determines the exact time scale.
For a more general and complete picture of the QSL the reader is advised to refer to the following references. 176,178 Similarly to the QSL, the information speed limit (ISL) can also restrict the minimum length of the control pulse. Be-hind this is the idea that the information encoded in the control pulse has to be sufficient to steer the system to the target. For example, in the noiseless case the degrees of freedom in the control (the number of independent frequencies in a bandwidth-limited control field or the number of kicks in a bang-bang control 179 ) should at least reflect the dimension of the system. 158 Note that in the presence of noise in the system or in the controls, more degrees of freedom are required to transmit the same amount of information.

IV. QOC FOR NV CENTERS
We now take a look at various applications of QOC to NV center-based systems including quantum sensing (section IV A), quantum information and quantum computation (section IV B).

A. Quantum Sensing
The sensitivity of NV center-based sensors depends on the coherence times (T 2 ) of the spin states, which in turn are theoretically only strictly limited by their lifetimes. Other limiting factors include efficient and coherent spin manipulation, which is inherently subjected to experimental imperfections and limitations. Several material synthesis techniques 13 have been developed to fabricate NV centers in ultrapure host crystals with minimal noise due to paramagnetic impurities and crystal field induced phenomena (strain) that diminish the coherence time (T 2 ) of the spin qubit. Additionally, dynamical decoupling protocols discussed in section II A 1 have shown promising results enhancing the coherence time. Yet, the currently achieved sensitivities fall short of the theoretical limits set by the spin state lifetimes. A material independent approach, which is the focus of this review, is to design specialized spin manipulation protocols that are optimised for efficiency in consideration of noise and experimental limitations/imperfections. The critical processes in the sensing techniques discussed so far are spin state initialization, state-to-state transfer and spin state readout. All these steps require the system to evolve in a specified way given a certain set of constraints. Broadly speaking, this catches the essence of quantum optimal control (QOC) theory, as described in section III. QOC techniques have shown promising results in fields like NMR 102,180 and atomic physics related experiments. 128,181 In the past decade, a number of interesting results have been obtained by using optimal control for NV spin system.
Häberle et al. 110 showed that quantum limited sensitivities for magnetic field sensing can be achieved using QOC based spectroscopy techniques (Fig. 19). In this work, the team used a single NV qubit system to image nanoscale magnetic fields with scanning probes. The aim of the work was to obtain a wider dynamic range for the spectroscopic microwave (MW) pulses using QOC . An open-loop numerical optimisation technique (GRAPE, see section III B 1 for more details) was used to obtain frequency selective, high bandwidth MW pulses for state transfer. Their results showed photoshot-noise limited sensitivities of around 4.5 µT

√
Hz and a dynamic range of more than 2.2 mT, as well as improved robustness against fluctuations in MW power. Nöbauer et al. 123  proach for open-loop optimisation limiting MW pulses to a certain frequency domain (Fig. 20). They obtained pulses that were robust against MW amplitude variation and frequency detuning. They first demonstrated the working principle employing a single NV using different MW detunings from the resonant frequency. They further concluded that the optimised pulse is ideal for NV ensembles that require the same spin manipulation pulse to be effective for a number of systems with different resonant frequencies and hence detunings as seen in the matrix plot in Fig. 20 shows the comparison of regular rectangular control pulses and the optimised pulses, demonstrating two orders of magnitude enhancement in sensitivity. This is demonstrated in a proof-of-principle experiment with a spin echo sequence for magnetic field sensing. Poggiali et al. 120 demonstrated a different approach (Fig. 21) making use of a modified cost function based on the Fisher information to enhance the sensitivity of NV centers to a wide range of AC magnetic fields (Fig. 21). The Fisher information of the measurement (see section III A 2 specifically Eq. (10)) takes into account the signal of interest as well as the noise. The sensitivity is related to the Fisher information via Eq. (9). A gradient-free minimisation technique (see section III B 2) was used to find the optimal temporal distance between pulses and adjust the initial phase of the signal. Finally, they experimentally demonstrated enhanced sensitivity by up to a factor of two for a single qubit system compared to a CPMG 183 pulse sequences (a type of dynamical decoupling sequence, see section II A 1). Along this line, Müller et al. 119 showed how optimal control-designed frequency filter functions allow a speed-up of the measurement and fast detection of fluctuating signals. In another work, Scheuer et al. 108 demonstrated a novel technique for spin qubit control in the ultra fast driving regime beyond the rotating wave approximation (see Appendix B). When the Rabi frequency of the drive becomes comparable to the transition frequency (in the reference about 30 MHz), the RWA breaks down and the system is no longer described by Eq. (6). To overcome this artificial constraint on the pulse duration, an optimal control pulse considering the counterrotating terms was designed. They used the CRAB algorithm (see section III B 3) to optimise standard πand π 2 -pulses in this regime and experimentally demonstrated Ramsey and spin-echo sensing protocols.
Ziem et al. 184 optimised resonance imaging of 19 F nuclei using GRAPE. The nuclei were part of patterned calcium fluoride on the diamond surface. The optimised pulses significantly improved the robustness of their DD protocol against variations in the driving field in comparison to a standard XY16-N sequence. Recently, Konzelmann et al. 109 showed that QOC pulses can enhance the robustness of temperature sensing for biological applications. They used nanodiamonds in an agarose matrix to demonstrate enhanced signal quality for fast temperature fluctuation measurements in dynamic biological media. For this purpose, a cooperative D-Ramsey pulse sequence, specifically designed for temperature measurements, 185 was optimised using the MATLAB based DYNAMO package that uses a GRAPE optimisation algorithm (see section III B 5 for more details). Here, the typical state-to-state transfer fidelity serves as the cost function for the optimisation (Fig. 22).

B. Quantum Information and Computation
The usefulness of QOC techniques has been demonstrated repeatedly in the field of quantum information and computation. As the basic theoretical protocols and experimental setup need to fulfill similar fundamental demands as for the described sensing techniques, similar optimisation strategies come into play. Dolde et al. 73 realised entanglement of two nuclear spins over a distance of 25 nm by entangling them with one NV center each which in turn were coupled through dipole-dipole interaction. They used GRAPE (see section III B 1) to realise optimised PSWAP and NOT gates which worked despite the hyperfine interactions which interfere with standard controls. The optimised pulses allowed for 20 NOT gate repetitions without significant loss of fidelity, while standard pulses already showed poor performance after a single gate. Waldherr et al. 81 demonstrated three-qubit phase-flip error correction using three nuclear spins and the NV center as an ancilla. While the system was entirely manipulated via the NV center spin, the hyperfine couplings complicated its control. This was solved by applying microwave pulses of two different frequencies simultaneously, finding their shape with a GRAPE-implementation. In order to construct optimal π/2as well as spin-inversion pulses, Frank et al. 107 designed and experimentally implemented closed-loop optimisation for single NV centers (Fig. 23). Their dCRAB implementation autonomously found an optimal solution for the desired goals within the experimental error. Their techniques are translatable to a number of quantum computing applications as π/2-pulses form the building blocks of most common quantum gates. In case of moderate detuning, the closed-loop pulses outperformed any open-loop-generated sequences. More details on the algorithm can be found in section III B 4. Said et al. 112 compared different strategies to entangle FIG. 22. Robust and efficient quantum optimal control of spin probes in a complex (biological) environment: (a) Working principle for the Cooperative D-Ramsey pulse sequence for temperature sensing based on the NV center ground state three level system. The sequence is sliced into 3 segments (seg1, seg2 and seg3), each separated in time by τ/2. ϕ D and ϕ B denote phases accumulated due to change in the zero-field parameter D and due to the external magnetic field, respectively. Optimal control pulses are used for spin projection to include cooperative design in the protocol to cancel out the undesired effects in phase accumulation. (b) The plot shows the measurements of temperature fluctuations over time using this protocol. Refer to Konzelmann et al. 109 for more details. (adapted and modified under the terms of the Creative Commons Attribution 3.0 International License) an NV center with a 13 C nuclear spin, namely sequential pulses, composite pulses, and numerically-optimised pulses (GRAPE). They concluded that the optimised pulse, did not only outperform the others in robustness but was also faster than the composite pulse. Tsurumoto et al. 111 utilised a GRAPE algorithm to optimise a state-to-state transfer pulse applied to a NV center, which was utilised to transfer photon polarization to a nuclear spin qubit. Recently, Chen et al. 186 proposed a QOC technique based on a combination of open-loop and closed-loop optimisation to demonstrate the working principle of a NV-based quantum processor.

V. CONCLUSION
NV centers are a valuable platform for a whole family of quantum applications. As the boundaries of existing implementations of different quantum sensing and quantum computation schemes are being explored, QOC offers a route to transverse these limits. We have reviewed a number of existing applications of QOC to NV center-based systems and provided a recipe for the application of QOC with a special focus on NV centers. It has been shown that QOC methods can increase the precision in qubit control and manipulation, especially in the presence of environmental factors such as surrounding nuclei. These advantages have also been exploited to improve NV-based sensors. Their versatility in combination with QOC shows to be especially effective when the limiting factor is defined by the experimental constraints and/or limited knowledge of the system. Enhanced control also provides a mean to efficiently implement quantum algorithms in NV center-based quantum registers. All combined, QOC enhances the capabilities which are crucial to building a NV center-based quantum computer. QOC methods provide the boost that may lead diamond-based quantum systems into the realm of commercial quantum technologies.

ACKNOWLEDGMENTS
This work has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement N • "765267" (QuSCo), N • "820394" (AsteriQs), and N • "817482" (PASQuanS). This work was partially supported by the Italian PRIN 2017, the QuantERA project QTFLAG and QuantHEP, the DFG via the TWITTER project.

Appendix A: Spin Hamiltonian
In a solid state system, spin interactions are mostly mediated through magnetic fields. In general the spin interaction energies are very small in comparison to electronic interaction energies (Fig. 3)  E 0 represents all the non spin-interaction energies andŝ n are the relevant spin bases (i.e. of different atoms, the NV center itself etc.). The Hamiltonian for such spin interaction energies is called the spin Hamiltonian. For a NV center like system, the spin Hamiltonian can be obtained in terms of spin operators for the corresponding singlet spin state basis: Terms of different order in the expansion A1 represent different kinds of spin interactions. For our purposes, it suffices to consider the terms up to second-order. Consequently, all the different types of interactions can be formalised in the same way through interaction tensors and spin operators. We now look into the details of these different types, for simplicityh is assumed to be equal to 1 in the following equations.
Linear terms: These terms mainly represent the interaction of the spin with external fields, whose origin may also lie within the crystal.
where,ˆ S is the spin operator, B is the magnetic field, and Z is the Zeeman interaction tensor. In the Hamiltonian discussed in section I B the magnetic interaction term is of this nature. The coupling tensor in that case is a diagonal matrix multiplied with the gyromagnetic ratio of the electron.
Bi-linear terms: These terms represent spin-spin interactions such as dipolar coupling, exchange interactions, and hyperfine couplings.
where,ˆ I is the nuclear spin operator vector, and N is the hyperfine coupling tensor. In the Hamiltonian presented in Eq. (1) in section I B the interaction of the NV spin with other spins is bi-linear in nature. In this case, the hyperfine coupling tensor may be simplified and described by the two main contributing interactions; the axial coupling constant N axial and the transverse coupling constant N tran : Note that the Fermi contact interaction term is given by 187 and the dipole interaction term is given by f N is an order of magnitude larger than d N for the 14 N and 15 N nuclei. 187 We may specifically consider the interaction between two different NV center spins. The Hamiltonian for these dipoledipole interactions can be written aŝ whereŜ 1 andŜ 2 are the spin operators for the respective NV centers and r is the spatial vector which joins them. This interaction is relatively weak in comparison to other spin interactions. Nevertheless, NV-NV interactions have been a matter of interest in several works in literature (for example see Dolde et al. 73 ).
Quadratic terms: These terms usually represented the interaction of spin with itself, although the source of the interaction can indirectly be due to an external field. Nuclear quadrupole interaction in NMR and electron zero-field splitting terms are of this nature. They are represented as: where, D is the quadrupole coupling tensor. D is in good approximation a symmetric matrix. The zero-field splitting terms for the Hamiltonian in Eq. (1) in section I B are of this nature. For a NV center this tensor can, to a good approximation, 36 be written in the basis of spin matrices as: where D is the zero-field splitting and E is the non-axial zerofield parameter. For most practical purposes one can consider E ≈ 0. The effect of an electric field on the spin is rather more complicated and for a solid state system it depends on symmetry conditions. For a C 3v symmetry system like the nitrogen-vacancy center, the effect of a linear electric field can be described by the following approximate Hamiltonian: 68,188 H elec = E X −P 11 (Ŝ XŜ Y +Ŝ YŜ X ) + P 15 (Ŝ YŜ Z +Ŝ ZŜ Y + E Y P 11 (Ŝ 2 X −Ŝ 2 Y ) + P 15 (Ŝ XŜ Z +Ŝ ZŜ X ) where, E i are the component of the effective electric field. External strain/pressure is manifested as a crystal strain electric field and is incorporated in the effective electric field. P i j are the components of a third rank coupling tensor 188 . These coupling constants can be determined by electron paramagnetic resonance experiments. In practice P 15 arises from the mixing of the m s = ±1 and m s = 0 spin states 188 , and is negligible with respect to the other coupling terms.

Appendix B: The Rotating Wave Approximation
In the following, the approximate Hamiltonian for a NV qubit is derived using the rotating wave approximation (RWA). It is a simplified version of the second term of Eq. (1). First, let us consider a single NV center in a magnetic field. The magnetic field is composed of two parts: The static component B points along z and a microwave field with amplitude 2B ⊥ oscillating along x with a frequency ω mw and phase ϕ. The magnetic field B can thus be written in a slightly odd form, the usefulness of which will become apparent in the following steps.
The splitting of the energy levels of the NV center ground state due to the static magnetic field allows the description as a two level system only considering |0 and |1 .ˆ σ are the Pauli matrices used to describe the corresponding spin operators which are denoted withˆ s =¯h 2ˆ σ to distinguish them from the three level systemˆ S and genericˆ s spin operators. The resulting Hamiltonian reflects the interaction between the field and the NV spin.
In the rotating frame of the microwave field with the transformation U = exp(iω mwŝZ t/h) this Hamiltonian becomes: Using the following relations the corresponding spin operators may be determined.
Uŝ X U † = cos(ω mw t)ŝ X − sin(ω mw t)ŝ Y Uŝ Y U † = sin(ω mw t)ŝ X + cos(ω mw t)ŝ Y The different terms in the rotating frame, after some algebra and trigonometric identities, can be written as UÂ + U † = B ⊥ (cos((ω mw − ω mw )t + ϕ)ŝ X + sin((ω mw − ω mw )t + ϕ)ŝ Y ), One may notice that the terms in UÂ − U † are of much higher frequency than the rest. Hence, in the rotating wave approximation, it is assumed that they average out and their contribution is negligible. Under this approximation, the Hamiltonian in the rotating frame depends only on the Rabi frequency Ω = γ nv B ⊥ , the detuning ∆ = ω nv − ω mw , with the NV center's resonant frequency ω nv = γ nv B , and the phase ϕ: H RWA ≈H NV−B = ∆ŝ Z + Ω (cos(ϕ)ŝ X + sin(ϕ)ŝ Y ) .