Robust and efficient quantum optimal control of spin probes in a complex (biological) environment. Towards sensing of fast temperature fluctuations

Spin quantum probes have the potential to measure a wealth of parameters with unprecedented accuracy and spatial resolution. One representative of such quantum sensors, is the negatively charged nitrogen-vacancy center (NV) in diamond. Several applications as a magnetic [1] and electric field [2], pressure [3] and temperature [4–12] probe have been demonstrated. In addition, one can also extract chemical information about samples by nuclear magnetic resonance with ∼Hz spectral and nanometer-scale spatial resolution [13, 14]. These techniques can be modified to measure the temporal dynamics of the quantity under study using NV [15, 16].

The sensor capability of NV can be preserved in nanodiamonds (NDs), even in a biological environment [17], where spin coherence and relaxation times are greatly reduced. As the NV is a spin one system (S = 1), the response to external magnetic fields is anisotropic: tumbling of the host ND will lead to variations in the electron spin resonance (ESR) transitions strength and frequency. Both depend on the orientation of the NV axis (the defect symmetry axis) with respect to the external static magnetic and microwave field. Therefore, proper alignment of the NV spin system is usually mandatory to perform precise quantum measurements.

To counteract variations in the NV spin transitions frequency and its driving strength, one can utilize optimal control theory, which has been used for NV spin control in a variety of applications [18–22]. In general, the strategy is to numerically optimize a pulse to achieve a desired quantum state.

To perform the optimization we use the quantum optimization package DYNAMO [23], which uses the gradient ascent pulse engineering algorithm together with Broyden–Fletcher–Goldfarb–Shanno minimization. DYNAMO can also optimize open quantum system dynamics. By introducing cooperativity between pulses of a given sequence via the quantum state filter method by Braun and Glaser [24], possible errors in phase adjustment can be compensated by the following pulses. Overall this reduces the demand for individual pulses, while still achieving the overarching goal of the entire sequence utilizing less resources like time or energy [24]. To demonstrate the power of the used optimization scheme, we modified the dynamical decoupling sequence D-Ramsey, optimized for temperature sensing with NV [4], to a new sequence, which allows us to monitor local temperature even with tumbling NDs.

The original D-Ramsey involves the allowed ESR transitions $| 0\rangle \leftrightarrow | +1\rangle$ and $| 0\rangle \leftrightarrow | -1\rangle$ of the triplet ground state of the NV center (see figure 1(a)). These two ESR lines can shift due to changes in a variety of internal and external parameters. For example, if the zero field splitting D of the NV changes by local temperature changes [25], they shift in the same direction. A change of the strength of an applied axial magnetic field B₀, shift the ESR lines in opposite directions. The D-Ramsey is designed in a way, that only the center of gravity (COG) of both transitions is measured: hereby the sequence acts as a Ramsey accumulating the phase φ_D = δD · τ during the free evolution interval τ, if the COG is detuned by δD. Slow changes in B₀, introducing a possible phase shift φ_B, are canceled out similar to a Hahn echo sequence. One can separate the D-Ramsey into three pulse segments with individual roles (see figures 1(a) and (b)). 'seg1' initializes quantum sensing, by creating a superposition between $| 0\rangle$ and $| -1\rangle .$ 'seg2' prepares the system to be only sensitive to COG shifts of $| -1\rangle$ and $| +1\rangle ,$ which is accomplished by swapping the populations of $| -1\rangle$ and $| +1\rangle .$ 'seg3' converts phase to state population enabling readout.

Zoom In Zoom Out Reset image size

To understand how we translate the D-Ramsey scheme into a cooperative pulse design, we quickly recall the basic idea of optimal control theory: the dynamics of a quantum system with density matrix ρ can be described by the von Neumann equation:

$$\begin{eqnarray}&&\dot{\rho }(t)=-{\rm{i}}\,\left[\mathop{\underbrace{{H}_{{\rm{drift}}}+\displaystyle {\sum }_{k}{u}_{k}(t){H}_{k}}}\limits_{H},\rho (t)\right],\end{eqnarray} \tag{ 1 }$$

with ${H}_{{\rm{drift}}}$ being the free evolution or drift Hamiltonian, while H_k are the control Hamiltonians that describe the interaction of the system with for example an external microwave field. For magnetic fields, H_k is composed of the spin operators S_x, S_y and S_z. The amplitude u_k(t) at time t of this interaction is typically called control field. The task is to find all u_k(t), such that the system evolves into a predefined final state after time T, represented by the density matrix λ. Therefore, one defines a quantity Φ₀ that can for example be the overlap of the desired final state λ with the actual state ρ evolved under the influence of the control fields: ${{\rm{\Phi }}}_{0}=\langle \lambda | \rho (T)\rangle .$ By maximizing Φ₀, one optimizes the amplitude u_k(t) to the designed functionality. Hereby one of the strengths of optimal control lies in the possibility to find a particular set of u_k(t) that allows robustness against variations of system parameters (e.g. a broad range of ${H}_{{\rm{drift}}}$ and H_k) [26].

To introduce cooperativity using the quantum state filter method, we describe the evolution of the system in the Liouville space. The evolution of a spin system with density matrix ρ_s that is coupled to an environment can be described as:

$$\begin{eqnarray}&&{\rho }_{{\rm{s}}}({t}_{2})=\exp ( {\mathcal L} ({t}_{1},{t}_{2}){\rm{\Delta }}t){\rho }_{{\rm{s}}}({t}_{1}).\end{eqnarray} \tag{ 2 }$$

Hereby ${\mathcal L} =-{\rm{i}}\cdot {\mathcal H} +{\mathscr{D}}$ is the Liouville superoperator, being the sum of the time independent system Hamiltonian ${\mathcal H}$ and the dissipator ${\mathscr{D}}$ in the time interval ${\rm{\Delta }}t={t}_{2}-{t}_{1}.$ The exponential term in equation (2) is typically called a dynamic map. By adding suitable Liouville space projectors P between control sequence segments, the full sequence can be optimized in one go invoking cooperativity. Thereby P projects onto a subspace in the Liouville space of the density matrices, which in our case is any equal superposition of states $| 0\rangle$ and $| -1\rangle$ after 'seg1' and $| 0\rangle$ and $+| 1\rangle$ after 'seg2'. If the state before applying the projector does not belong to this subspace, some purity is irreversibly lost and the fidelity of the final state reduces. This causes DYNAMO to optimize the amplitudes u_k(t) such that the desired subspace is reached before the projector is applied. After the optimization, we can simply remove the projectors and replace them with arbitrary but equally long free evolution intervals between the three segments.

As shown in the section Results and discussion, the introduced strategy allows us to measure the local temperature of a ND, despite its rotation, which would result in incoherent control using the standard sequence. Note, that in principle, one can apply this recipe to any measurement scheme (in particular multi-phase dynamic decoupling sequences).

In this section we define the NV spin Hamiltonians and projectors used to compile an optimal control pulse.

Description of NV spin system

To compile a Coop-D-Ramsey sequence (cooperatively numerically optimized D-Ramsey), one needs to describe the NV quantum system in terms of the system Hamiltonian and the set of parameters, which are varied by the tumbling motion of the NDs. We define ${H}_{{\rm{drift}}}$ and ${H}_{{\rm{c}}}=\displaystyle {\sum }_{k}{u}_{k}(t)\,{H}_{k}$ as:

$$\begin{eqnarray}&&{H}_{{\rm{drift}}}=D\cdot {S}_{z}^{2}+E({S}_{x}^{2}-{S}_{y}^{2})+\displaystyle {\sum }_{k=x,y,z}{\beta }_{0k}{S}_{k}\,,\end{eqnarray} \tag{ 3a }$$

and

$$\begin{eqnarray}&&{H}_{{\rm{c}}}(t)=\displaystyle {\sum }_{k=x,y,z}{\beta }_{1k}(t){S}_{k}.\end{eqnarray} \tag{ 3b }$$

D is the zero field splitting, E is the zero field component introduced by off-axial strain, S_k the spin operators of a spin one system, β_0k static components introduced by an external magnetic field with strength ${\beta }_{0}=\sqrt{\displaystyle {\sum }_{k}{\beta }_{0k}^{2}},$ and β_1k the components of the microwave field.

We distinguish two cases. The first one is a NV with $| E| \ll {\beta }_{0},$ in which the external static field is aligned along the z-axis of the NV spin system: $| {\beta }_{0z}| ={\beta }_{0}.$ We call this bulk case. The second one is a NV with $| E| \gg {\beta }_{0},$ called the ND case. The microwave field is chosen to be an in-plane linear polarized field. Transforming equations (3), into a left-handed and a right-handed rotating reference frame, leads to the overall system Hamiltonian H for the bulk case:

$$\begin{eqnarray}&&H={H}_{{\rm{drift}}}+{u}_{x}(t)\cdot X+{u}_{y}(t)\cdot Y,\end{eqnarray} \tag{ 4a }$$

with

$$\begin{eqnarray}{H}_{{\rm{drift}}}=\left(\begin{array}{ccc}{\rm{\Delta }}D+{\beta }_{0z} & 0 & 0\\ 0 & 0 & 0\\ 0 & 0 & {\rm{\Delta }}D-{\beta }_{0z}\end{array}\right)\end{eqnarray} \tag{ 4b }$$

and

$$\begin{eqnarray}X=\displaystyle \frac{{\beta }_{1}}{2}\left(\begin{array}{ccc}0 & 1 & 0\\ 1 & 0 & 1\\ 0 & 1 & 0\end{array}\right),\,Y=\displaystyle \frac{{\beta }_{1}}{2}\left(\begin{array}{ccc}0 & -{\rm{i}} & 0\\ {\rm{i}} & 0 & {\rm{i}}\\ 0 & -{\rm{i}} & 0\end{array}\right).\end{eqnarray} \tag{ 4c }$$

For equations (4) we make use of the assumption $D\gg {\beta }_{1},$ whereas β₁ is the measured Rabi frequency, which is determined experimentally before pulse compilation. The parameter ΔD denotes the detuning of the zero field splitting D from a reference, which may be the former value of D before local temperature has changed. In later study we introduced a detuning by changing the reference rather than the actual value of D.

Note that the described relations in equations (4) are valid even for a residual microwave field aligned along the z-axis.

For example, if we introduce a cosine modulation linear polarized microwave field along the z-axis with amplitude β_1z, equation (4a) acquires an additional term H_z:

$$\begin{eqnarray}&&{H}_{z}={\beta }_{1z}\cdot \,\cos (Dt)\cdot {S}_{z}.\end{eqnarray} \tag{ 4d }$$

As D is typically 2.87 GHz for NV, and the Rabi amplitude is on the order of several MHz, fast fluctuations introduced by the S_z operator are averaged out. In other words, one has to consider the projection of the driving field onto the xy-plain of the NV spin system only. Tumbling will lead to modulation of the in-plane microwave field amplitude, which we denoted B_xy in figures 2(c) and 3(c). To allow a variation in the driving strength of the microwave field, we introduce the relative amplitude κ. The controls X_l and Y_l for a single set l of the whole parameter set to optimize for, are modified to:

$$\begin{eqnarray}{X}_{l}=\kappa (l)\displaystyle \frac{{\beta }_{1}}{2}\left(\begin{array}{ccc}0 & 1 & 0\\ 1 & 0 & 1\\ 0 & 1 & 0\end{array}\right)\,{\rm{and}}\,{Y}_{l}=\kappa (l)\displaystyle \frac{{\beta }_{1}}{2}\left(\begin{array}{ccc}0 & -{\rm{i}} & 0\\ {\rm{i}} & 0 & {\rm{i}}\\ 0 & -{\rm{i}} & 0\end{array}\right).\end{eqnarray} \tag{ 4e }$$

The variation in κ will be accomplished by changing the signal generator's amplitude.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the ND case the situation becomes more complicated:

$$\begin{eqnarray}&&H={H}_{{\rm{drift}}}+{u}_{x}(t)\cdot [{X}_{\parallel }+{X}_{\perp }]+{u}_{y}(t)\cdot [{Y}_{\parallel }+{Y}_{\perp }]\end{eqnarray} \tag{ 5a }$$

with:

$$\begin{eqnarray}{H}_{{\rm{drift}}}=\left(\begin{array}{ccc}{\rm{\Delta }}D+E & 0 & 0\\ 0 & 0 & 0\\ 0 & 0 & {\rm{\Delta }}D-E\end{array}\right),\end{eqnarray} \tag{ 5b }$$

and:

$$\begin{eqnarray}\begin{array}{l}{X}_{\parallel }=\displaystyle \frac{{\beta }_{1\parallel }\,}{2}\left(\begin{array}{ccc}0 & 1 & 0\\ 1 & 0 & 0\\ 0 & 0 & 0\end{array}\right),{Y}_{\parallel }=\displaystyle \frac{{\beta }_{1\parallel }}{2}\left(\begin{array}{ccc}0 & -{\rm{i}} & 0\\ {\rm{i}} & 0 & 0\\ 0 & 0 & 0\end{array}\right),\\ {X}_{\perp }=\displaystyle \frac{{\beta }_{1\perp }}{2}\left(\begin{array}{ccc}0 & 0 & 0\\ 0 & 0 & 1\\ 0 & 1 & 0\end{array}\right),{Y}_{\perp }=\displaystyle \frac{{\beta }_{1\perp }}{2}\left(\begin{array}{ccc}0 & 0 & 0\\ 0 & 0 & {\rm{i}}\\ 0 & -{\rm{i}} & 0\end{array}\right).\end{array}\end{eqnarray} \tag{ 5c }$$

In equation (5c) ${\beta }_{1\parallel }$ and ${\beta }_{1\perp }$ are again directly the measured Rabi frequency of left and right ESR transitions, respectively. The interpretation of equation (5c) should be understood in the following way: if a strain field is aligned along the microwave field only one transition can be driven. If it is perpendicularly aligned, one can only drive the other transition. Again, for parameter variation equation (5c) is redefined as:

$$\begin{eqnarray}\begin{array}{l}{X}_{l}=\displaystyle \frac{1}{2}\left({\beta }_{1\parallel }\begin{array}{ccc}0 & {\beta }_{1\parallel }\,\kappa (l) & 0\\ \kappa (l) & 0 & {\beta }_{1\perp }\kappa ^{\prime} (l)\\ 0 & {\beta }_{1\perp }\kappa ^{\prime} (l) & 0\end{array}\right),\\ {Y}_{l}=\displaystyle \frac{{\rm{i}}}{2}\left(\begin{array}{ccc}0 & -{\beta }_{1\parallel }\kappa (l) & 0\\ {\beta }_{1\parallel }\,\kappa (l) & 0 & {\beta }_{1\perp }\kappa ^{\prime} (l)\\ 0 & -{\beta }_{1\perp }\kappa ^{\prime} (l) & 0\end{array}\right).\end{array}\end{eqnarray} \tag{ 5d }$$

In case of the fixed ND study $\kappa (l)$ equals $\kappa ^{\prime} (l).$ For the rotating ND $\kappa (l)$ and $\kappa ^{\prime} (l)$ are varied individually, but ${\beta }_{1\parallel }$ is set equal to ${\beta }_{1\perp }.$

Projectors for cooperative pulse design

To introduce cooperativity, we use the projector sets ${P}_{xy}^{(1)}/{P}_{xy}^{(2)}$ (set1) and ${P}_{xy}^{(1b)}/{P}_{xy}^{(2b)}$ (set2) as shown in equations (6) in two sub-ensembles. All P_xy project into a two-dimensional subspace of the Liouville space of density matrices, with the difference, that projectors with the superscript (b) introduce an additional $\tfrac{\pi }{2}$ phase shift. The reason to do so, is that we want to account for the situation, where some phase is picked up during the free evolution intervals with length $\tfrac{\tau }{2}.$ The minimum set to account for all phase accumulations between zero and 2π is covered by both. Note that ${P}_{xy}^{(1b)}$ adds a positive phase, whereas ${P}_{xy}^{(2b)}$ the corresponding negative one. If we would have chosen a positive phase in the latter the final state must be $| +1\rangle$ instead of $| 0\rangle$

$$\begin{eqnarray}\begin{array}{l}{P}_{xy}^{(1)}=\left(\begin{array}{ccccccccc}0 & \, & \, & \, & \, & \, & \, & \, & \,\\ \, & 0 & \, & \, & \, & \, & \, & \, & \,\\ \, & \, & 0 & \, & \, & \, & \, & \, & \,\\ \, & \, & \, & 0 & \, & \, & \, & \, & \,\\ 0.5\, & \, & \, & \, & 0.5 & \, & \, & \, & 0.5\\ \,\, & \, & \, & \, & \, & 1 & \, & 0 & \,\\ \, & \, & \, & \, & \, & \, & 0 & \, & \,\\ \, & \, & \, & \, & \, & 0 & \, & 1 & \,\\ 0.5\, & \, & \, & \, & 0.5 & \, & \, & \, & 0.5\end{array}\right),\\ \,{P}_{xy}^{(2)}=\left(\begin{array}{ccccccccc}0.5 & \, & \, & \, & 0.5 & \, & \, & \, & \,0.5\\ \, & 1 & \, & 0 & \, & \, & \, & \, & \,\\ \, & \, & 0 & \, & \, & \, & \, & \, & \,\\ \, & 0 & \, & 1 & \, & \, & \, & \, & \,\\ \,0.5 & \, & \, & \, & 0.5 & \, & \, & \, & \,0.5\\ \, & \, & \, & \, & \, & 0 & \, & \, & \,\\ \, & \, & \, & \, & \, & \, & 0 & \, & \,\\ \, & \, & \, & \, & \, & \, & \, & 0 & \,\\ \, & \, & \, & \, & \, & \, & \, & \, & 0\end{array}\right),\end{array}\end{eqnarray} \tag{ 6a }$$

$$\begin{eqnarray}\begin{array}{l}{P}_{xy}^{(1b)}=\left(\begin{array}{ccccccccc}0 & \, & \, & \, & \, & \, & \, & \, & \,\\ \, & 0 & \, & \, & \, & \, & \, & \, & \,\\ \, & \, & 0 & \, & \, & \, & \, & \, & \,\\ \, & \, & \, & 0 & \, & \, & \, & \, & \,\\ 0.5\, & \, & \, & \, & 0.5 & \, & \, & \, & 0.5\\ \, & \, & \, & \, & \, & {\rm{i}} & \, & 0 & \,\\ \, & \, & \, & \, & \, & \, & 0 & \, & \,\\ \, & \, & \, & \, & \, & 0 & \, & \mbox{--}{\rm{i}} & \,\\ 0.5\, & \, & \, & \, & 0.5 & \, & \, & \, & 0.5\end{array}\right),\\ {P}_{xy}^{(2b)}=\left(\begin{array}{ccccccccc}0.5 & \, & \, & \, & 0.5 & \, & \, & \, & 0.5\,\\ \, & {\rm{i}} & \, & 0 & \, & \, & \, & \, & \,\\ \, & \, & 0 & \, & \, & \, & \, & \, & \,\\ \, & 0 & \, & \mbox{--}{\rm{i}} & \, & \, & \, & \, & \,\\ \,0.5 & \, & \, & \, & 0.5 & \, & \, & \, & \,0.5\\ \, & \, & \, & \, & \, & 0 & \, & \, & \,\\ \, & \, & \, & \, & \, & \, & 0 & \, & \,\\ \, & \, & \, & \, & \, & \, & \, & 0 & \,\\ \, & \, & \, & \, & \, & \, & \, & \, & 0\end{array}\right).\end{array}\end{eqnarray} \tag{ 6b }$$

To demonstrate that it is mandatory to use both projector sets, we compiled in appendix B one sequence with ${P}_{xy}^{(1)}$ and ${P}_{xy}^{(2)}$ only, and one including all projectors defined in equations (6). Only in the latter case the correct response of the D-Ramsey can be resampled.

As we stated in the introduction part, optimal control theory can be utilized to introduce a certain robustness against parameter variations. For demonstration we use the sequence compiled for later study, which we refer to as bulk case (see also the latter section). Hereby we compare the Coop-D-Ramsey (cooperatively numerically optimized D-Ramsey) to a typical D-Ramsey consisting of conventional π and π/2 pulses. Figure 1(c) shows the population of all three states for different ${H}_{{\rm{drift}}}$ and H_k, for three different situations. The definition of ${H}_{{\rm{drift}}}$ and H_k can be found in equations (4b), (4c), (4e). The upper left plot pictures the population of all three states for a well-set ideal parameter choice in case of a D-Ramsey. The sequence works as desired. The lower plots represent the average over a sample of 63 different sets (for parameters see bulk study and Methods) of ${H}_{{\rm{drift}}}$ and H_k for a D-Ramsey and Coop-D-Ramsey sequence. As one can see, for both, D-Ramsey and Coop-D-Ramsey, the population has a high divergence within the individual segments.

However, the Coop-D-Ramsey converges at the end of each segment to the desired superposition state. To calculate the optimal set of u_k(t) for the Coop-D-Ramsey, we used the projectors, which are described in equations (6). We labeled the projectors between 'seg1' and 'seg2' with ${P}_{xy}^{(1)}$ and ${P}_{xy}^{(1b)}$ and between 'seg2' and 'seg3' ${P}_{xy}^{(2)}$ and ${P}_{xy}^{(2b)}$ (see figure 1(b)). Note that we have to use two sub-ensembles: one with the set ${P}_{xy}^{(1)}$ and ${P}_{xy}^{(2)}$ and the other one with the set ${P}_{xy}^{(1b)}$ and ${P}_{xy}^{(2b)}$ for a successful compilation of the Coop-D-Ramsey sequence (see discussion in appendix B and figure B1). In addition, we design the system to end up in $| 0\rangle ,$ not in $| +1\rangle$ as this is the case for the D-Ramsey.

Next, we compared the performance of a D-Ramsey versus the Coop-D-Ramsey scheme in bulk diamond. To this end, we chose a single NV in bulk diamond and aligned an external weak magnetic field roughly along the NV axis. The splitting of the ESR transitions is 25.7 MHz. The Rabi frequency in the reference configuration (relative amplitude κ = 1) of the energetically lower and higher spin transitions was 3.8 MHz and 3.5 MHz, respectively. We explain the slight difference by the frequency dependence of the microwave power of our apparatus. To determine the maximum spin contrast we performed a spin contrast measurement (T₁-Decay, see figure 2(d) (black)) as described in appendix A or [27]. Then we acquired a Hahn echo to verify, that the coherence time of the NV center is long enough (figure 2(d) (gray)). In the next step we measured the Coop-D-Ramsey and D-Ramsey by varying amplitude and detuning to simulate tumbling and temperature changes. To calculate the sequence, we use the spin system as defined in equations (4a)–(4c), (4e). Amplitude variations are accomplished by varying the relative amplitude κ and detunings by varying ΔD as shown in figure 2(c) and described in former section. We fitted every data set to a harmonic function with an exponential damping term and extracted the peak to peak amplitude (signal contrast). The signal shows an almost constant contrast within the amplitude and detuning variation limit that has been set for compilation (κ: 0.7 to 1.3, ΔD: 0 to 200 kHz, more details in appendix A.2). Interestingly, the sequence has a better performance than the target limit. As one can see in figure 2, the signal contrast of the optimal control sequence is almost twice as large as in the non-optimized case, whereas in case of detuning, the signal contrast is always better within a bandwidth of 1 MHz.

For the ND case we performed a similar analysis. First, we compensated the external magnetic field to satisfy the condition $E\gg {\beta }_{0z}$ as required to use equations (5a)–(5d) for modeling the NV spin system within a highly strained ND. Fitting a continuously driven optically detected magnetic resonance (cwODMR) spectrum as shown in figure C2 reveals $E=\,4.35\,\pm 0.01\,{\rm{MHz}}$ and ${\beta }_{0z}=\,0.53\,\pm 0.03\,{\rm{MHz}}.$ The reference Rabi frequencies on the lower and upper transition are 2.5 MHz and 2.8 MHz, respectively. Again, pulse parameters can be found in appendix A.2.

In contrast to bulk diamond the difference in signal contrast between D-Ramsey and Coop-D-Ramsey is less, but still significant (figure 3). We attribute this difference to the reduction in the hyperfine splitting introduced by the zero field parameter E, which leads to a better result in case of the D-Ramsey. For demonstration we plot the transitions frequencies in figure C1 for E = 5 MHz. With increasing external magnetic field, the splitting of the hyperfine lines develops from a denser packed ND case to the more bulk like behavior with equal hyperfine splitting of 2.16 MHz (see also discussion in appendix C). In the ND study presented, the hyperfine lines gather in the range of several hundred kilohertz (see figure C2), because of the low external magnetic field, whereas in case of the bulk NV the hyperfine lines are equally split by 2.16 MHz. The pulse bandwidth of conventional square pulses is directly proportional to the Rabi driving strength and in our case on the order of several MHz. Now, the more a single spin transition is detuned versus the applied frequency of a microwave control pulse, the less good a desired state is aligned. In case of the ND study the hyperfine lines are much closer within the driving bandwidth of the conventional square pulses used, as in the bulk study. As a consequence, the D-Ramsey signal contrast is different when comparing both studies.

In a next step, we tested the cooperative pulse design on a tumbling ND. Therefore, we incorporated NDs in an agarose matrix, which acts as a cage for the NDs. The latter are spatially fixed but still can change their orientation over time. It turns out that NDs are typically immobilized in the beginning. We used this time window for characterizing the system to extract the Rabi driving strength, the zero field splitting D and E and the NV spin properties. After several hours of analysis, the NDs started to tumble on various timescales with rotational periods ranging from several hours down to milliseconds. Previously McGuiness et al [17] observed tumbling of NDs in cells on similar timescales. We restricted our observations on NDs with slow angular changes compared to the overall sequence length, also because we want the driving strength to be constant during the application of one sequence interval (i.e. in agreement with the quasistatic parameter variation assumed for pulse engineering). In addition, we minimized the residual external magnetic fields by minimizing the linewidth and the COG of a cwODMR spectrum of a highly NV doped ND crystal (see Methods part in [1]), as the COG position depends on the angle an external magnetic field is applied versus the NV axis. This would be interpreted as a temperature shift in case of a slowly tumbling ND (see figure D1).

After compensating for the external magnetic field, we tested the performance of the Coop-D-Ramsey on a ND doped with a single NV containing a ¹⁵N atom. We estimated E to be around $0.2\pm 0.1\,{\rm{MHz}}$ and $| {\beta }_{0z}| =\,0.2\pm 0.1\,{\rm{MHz}}$ in the beginning. The maximum Rabi frequency on the lower and upper transitions are $7.7\pm 0.3\,{\rm{MHz}}$ and $9.0\pm 0.3\,{\rm{MHz}}.$ Different to the former cases we optimized the pulse for every combination of κ and κ' for the two transitions between 0.6 and 1.4 using equation (5d). After compiling a suitable sequence with DYNAMO, we continuously monitored the shift of the zero field splitting via a Coop-D-Ramsey for almost one day. In addition, we interleaved measurements of Rabi oscillations on both electron transitions and the overall temperature of the sample holder via a thermistor (see figure 4).

Zoom In Zoom Out Reset image size

The signal contrast of the Coop-D-Ramsey is around 80% in the beginning. Over almost two hours the zero field splitting D shifted by around 100 kHz. This corresponds to a heating of the confocal cell dish of roughly 1.5 K. This heating was not observed by the thermistor as it is monitoring the environmental temperature only at the sample chamber. During this heat-up phase the signal contrast of the Coop-D-Ramsey stayed almost constant, around 80%. After two hours the contrast dropped. Around 45 min later we reduced the detuning to D by 100 kHz, to stay closer to the defined compilation interval (allowed detuning: 0–200 kHz), also to rule out that the latter is the origin of the signal drop. A closer look on the acquired Rabi data in figures 4(b)–(e) reveals that the FFT shows a significantly larger linewidth resulting from a fast decay of the Rabi oscillations. We interpret this behavior as an averaging of the NV having different orientations to the driving field during tumbling [28]. Furthermore, we acquired the autocorrelation function ${g}^{(2)}(\tau )$ before and after the end of the temperature measurement series (figure E1(a)). Interestingly, the ${g}^{(2)}(\tau )$ measured after the temperature measurements does not approach one, i.e. the value for a totally uncorrelated signal. To find the timescale of the new occurring feature, we performed fluorescence correlation spectroscopy (FCS) from sub-μs to sub-second timescale. The acquired $G(\tau )$ function reveals a bunching feature with a time constant of ${\tau }_{l}=4.0\pm 0.2$ ms. As we do not see a significant increase in width of the point spread function (figure E1) in comparison to the measured point spread function for other fixed NDs, we identify the bunching feature in $G(\tau )$ and the offset in ${g}^{(2)}(\tau )$ to be a consequence of ND rotation.

Finally, we compare the results for the Coop-D-Ramsey to the original D-Ramsey in case of rotating NDs. The D-Ramsey sequence requires by design the NV electronic transitions to be selectively addressable. This condition is not fulfilled by NDs, as their transitions are typical split in the range of ∼kHz to ∼MHz [17, 29, 30]. Nevertheless, we can estimate the maximum signal contrast expected for a rotating ND with the given parameters. To this end we simulated the D-Ramsey for various angles of the microwave field versus the NV axis and averaged the signal. Figure F1 shows the result of this averaging in case of a more bulk or ND like behavior. For the given parameter set (spin transition energy, Rabi amplitude and hyperfine splitting) we would not have been able to observe characteristic modulations corresponding to the shift in D i.e. measure temperature.

There are alternative ways to measure temperature with rotating NDs. In this context one may not be interested in the method's sensitivity only, but also in the accessible bandwidth the individual sensing strategy provides. The zero phonon line energy and the Debye–Waller-factor (DWF) of NV are temperature dependent and allow an all optical approach [11, 12]. Thus, measuring both quantities is very easy to implement, and show a comparable sensitivity to microwave assisted techniques in case of DWF and NVs in NDs. Nevertheless, one has to use a NV ensemble to reach a sensitivity in the range of $\sim 100\,{\rm{mK}}\,{{\rm{Hz}}}^{-1/2}$ [12]. Both methods are compatible with FCS, because thermal information is read out optically via the fluorescent signal. FCS allows the analysis on arbitrary timescales. However, it is questionable, if one can clearly separate the temperature induced correlation signal from features that may origin from, for example, charge state dynamics [31] or rotation, and if the latter do not dominate the desired correlation signal.

The most simple measurement scheme utilizing the NV spin is based on continuous wave excitation and a time resolution of ∼10 μs has been demonstrated [8]. As continuous wave excitation is limited by ${T}_{2}^{* }$ [32], the expected temperature resolution is on the order of ∼0.1 K Hz^−1/2 to ∼1 K Hz^−1/2 for NV ensemble in NDs [5, 8]. However, pulse sequences specifically developed for thermal sensing are ultimately limited by T₂ [4, 5]. They offer a sensitivities of ∼0.1 K Hz^−1/2 with reasonable parameters for single NVs in NDs [4], and can be further improved by higher order decoupling strategies [7].

In this work we continuously readout a single statistically rotating ND. Hereby, every data point in figure 4(f) has been acquired for several minutes, which is a long time period versus the actual rotational correlation time of the ND. The given readout origins from a statistical average of all possible NV angular variations versus the external microwave field. Hereby, it should make less difference, if a single NV or a dense ensemble of N NVs is probed. As sensitivity scales inverse with (relative) signal contrast c_rel and ∼$\sqrt{N}$ [5], a factor of 5 can be gained with $N=\,100$ and ${c}_{{\rm{rel}}}=0.5,$ which we define in our case as the relative signal contrast to the maximum possible signal. Consequently, applying the presented approach to ND ensemble should allow a sensitivity of ≤30 mK Hz^−1/2, using ∼130 mK Hz^−1/2 as a reference for a single NV in ND [4]. All pulsed measurement schemes are naturally expandable to, for example, a so-called stimulated echo sequence. The latter allows to sense correlations on the timescale, which a specific sensing sequence is applied as a lower bound (NV: ∼sub-μs), and the probe spin relaxation time itself as an upper bound (NV: ∼ms). The practicability of such stimulated echoes has been shown in case of NV [15, 16].

Besides, one question arises, namely: which is the best method for robust spin manipulation in pulsed mode. In case of the D-Ramsey scheme one could also replace conventional pulses, by adiabatic passage [33]. For example, chirped pulses, allowing to create a state inversion with high fidelity and bandwidth, or the B₁ insensitive rotation (BIR) pulse, which allows to adjust a superposition state with arbitrary state composition. Here B₁ is the driving microwave field. Even if adiabatic pulses provide a robust way to manipulate the state of a spin system, the pulses itself have to be much longer than the corresponding effective Larmor frequency in the rotating reference frame to fulfill the adiabatic condition. This has to be the case even for the smallest B₁ field amplitude considered. In addition, the BIR pulse scheme can only provide a certain bandwidth in frequency space roughly corresponding to the range of the applied Rabi frequency [33]. Especially in bulk diamond and for a NV with ¹⁴N, the pulse has to have a bandwidth of around 4 MHz at least to drive all hyperfine lines, setting the lower limit of the Rabi driving strength. In the case of NDs, strain usually reduces the hyperfine splitting that shall lead to an increase in the overall performance of the BIR pulses. The problem of overall pulse length still remains.

The used cooperative design gives some advantages over conventional optimal control. In the latter case, each pulse of a sequence would be synthesized as a particular gate, which requires more resources (e.g. time or energy) as synthesis of a particular state transfer. In our cooperative design, not even a particular state transfer is required but just one out of a subset (e.g. any equal superposition state). Then, all different pulses in our sequence can correct errors in state adjustment of other pulses. A particular final state is achieved by cooperative adjustment of all pulses. This introduces an additional degree of freedom, allowing the algorithm to be more efficient concerning the overall length of the sequence [24], providing better results for shorter pulses.

In this work we have chosen an overall pulse length of around $2\,\mu {\rm{s}}$ in case of the rotating ND. In case of a simulated echo, which essentially consist of two identical sequences split by a correlation interval, the shortest echo sequence would be $4\,\mu {\rm{s}}$ long. Consequently, the maximum bandwidth one can access in terms of temporal resolution to measure temperature fluctuations, is about 250 kHz. If we add a free evolution interval on the order of $\tau =0.5\,\mu {\rm{s}},$ one can reach a sensitivity of ∼40 mK Hz^−1/2 in case of a ND ensemble (see equation in main text of [5], $c=0.05,$ $N=100,$ ${T}_{{\rm{coh}}}\hat{=}\tau =\,0.5\,\mu {\rm{s}}$). The overall pulse sequence length would increases to 5 μs, which still allows a bandwidth of 200 kHz. This is comparable to the time resolution reported previously [34], where the resonance frequency at ∼840 kHz of a mechanical oscillator was probed. Therefore, it would be possible to resolve temperature fluctuations with μs-time and potentially even nm-spatial resolution.

This might be of particular importance as recent experiments on living cell using fluorescent probes [35–39] have measured sizeable temperature gradients. Thermodynamic arguments however suggest [40], that changes of the average temperature due to endogenous thermogenesis should we very small, but do not exclude, that measurable temperature fluctuations may still exist. For example, Inomata et al [34] observed a characteristic change of heat generation from a burst-like to a continuous increasing behavior after stimulating brown fat cells with norepinephrine in a well isolated chamber design.

In addition, NDs may be placed on a small structure whose size is sub-micron down to several tens of nanometers, to analyze their thermal dynamics. In difference to scanning probe approaches, where a sharp but microscopic tip is used to probe the local heat dissipation [9, 10, 41], NDs can reach dimensions even below 10 nm and still host NV [29, 30]. As diamond itself has a high thermal conductivity and a low heat capacity, high thermal coupling and low perturbation of the device under study is possible in direct contact. Hence, high spatial resolution can be combined with fast temporal response, enable by the present measurement scheme.

In summary, we have converted the D-Ramsey scheme for sensitive temperature measurements into a more robust Coop-D-Ramsey, which is based on optimal control theory in combination with a cooperative pulse design. We demonstrate superior robustness against variation in the driving strength and resonance mismatching, compared to conventional soft or hard pulses. The recipe we exploited to enable robust spin control is universal and can be used to design other pulse schemes.

In addition, we show that the Coop-D-Ramsey even performs well, when a ND is tumbling within an agarose matrix. Combining all parameters discussed it is reasonable to state, that our work paves the way to measure sub 100 mK temperature fluctuations on microsecond time and nanometer length scale.

We acknowledge the support of the European Commission Marie Curie ETN 'QuSCo' (GA N°765267), the KIST Open Research Program (2E27801) and the German science foundation (SPP 1601). In addition, this work was supported by ERC grant SMel, the Max Planck Society and the Humboldt Foundation and the EU–FET Flagship on Quantum Technologies through the Project ASTERIQS. We acknowledge funding by the Volkswagen Foundation. We thank Thomas Schulte-Herbrüggen and Steffen Glaser for fruitful discussions.

The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript.

Appendix. A.1. Sample preparation

Appendix. A.2 Pulse parameters

Appendix. A.3 Experimental setup

Appendix. A.4 Data analysis

To demonstrate, that it is mandatory to use both projector sets, we compiled one sequence with ${{\boldsymbol{P}}}_{{\boldsymbol{xy}}}^{(1)}$ and ${{\boldsymbol{P}}}_{{\boldsymbol{xy}}}^{(2)}$ only, and one including all projectors given in equations (6).

Zoom In Zoom Out Reset image size

The system was defined as described in equations (4a)–(4c). The shift introduced by an external z-magnetic field has been set to 5 MHz, resulting in a splitting of 10 MHz for the ESR transitions. In addition, DYNAMO shall optimize the pulse response for three different configurations that resample the hyperfine structure of ¹⁴N. Therefore, additional magnetic shifts with 0 MHz and ± 2.16 MHz are introduced. To see how phase accumulates during the evolution interval, the detuning ΔD is set to 200 kHz. Doing so, one would expect a single harmonic modulation between state $| 0\rangle$ and $| +1\rangle$ only with a period of 5 μs matching the adjusted detuning, when varying the time τ of the evolution intervals. Simulating the state dynamics of the model NV for different evolution intervals reveal, that in case of set1 fast oscillations can occur for different hyperfine lines (see figure B1(a)). When using set1 and set2, one can observe a clear oscillation between state $| 0\rangle$ and $| +1\rangle$ with a period of 200 kHz for all hyperfine lines (figure B1(b)). The additional frequency components between 1 MHz and 4 MHz disappeared in the corresponding Fourier transformation (compare figures B1(c) and (d)).

To calculate the ESR transition of the NV center in NDs we use following Hamiltonian:

$$\begin{eqnarray}&&{H}_{{\rm{NV}}}=D\cdot {S}_{z}^{2}+E({S}_{x}^{2}-{S}_{y}^{2})+{\beta }_{0}\cdot {S}_{z}+S\cdot A\cdot I.\end{eqnarray} \tag{ C.1 }$$

Whereas I are the nuclear spin matrices and A the hyperfine interaction tensor [42]. For reasons of simplification we consider the S_zI_z hyperfine coupling only. Then, we calculate the Eigensystem of H_NV and determine the most prominent transitions introduced by the spin operators S_x and S_y in the transformed system. Figure C1(a) shows this calculation for varying β₀ for a zero field splitting E = 5 MHz and ¹⁴N. As one can see the hyperfine lines of the ${m}_{I}=\pm 1$ state overlap for ${\beta }_{0z}\ll E.$ With increasing magnetic field, the lines begin to separate. For ${\beta }_{0z}\geqslant \sim \tfrac{E}{2}$ the hyperfine lines are well separated. In figure C1(b) we plot the corresponding transitions in case of ¹⁵N.

Zoom In Zoom Out Reset image size

As an example, figure C2 shows the cwODMR spectrum of the ND used for the ND case study in the main text. The residual magnetic field along the z-axis correspond to 0.53 ± 0.03 MHz. The zero field parameter E equals $4.35\pm 0.01\,{\rm{MHz}}.$

Zoom In Zoom Out Reset image size

The NV can be used to measure temperature. The reason for the latter lies within the temperature dependence of the zero field splitting D that can be approximated to be linear at ambient conditions. Thereby the shift in D is around $70\mbox{--}80\,{\rm{kHz}}\,{{\rm{K}}}^{-1}$ [25]. For comparison the shift of a single spin transitions of NV is around $2.8\,{\rm{MHz}}\,{{\rm{G}}}^{-1}.$ So already for the Earth magnetic field with ∼0.5 G a line shift of ∼1.4 MHz is expectable and would correspond to a temperature shift around 20 K. Small changes in the sub-G regime may already be interpret as a temperature change. Therefore, the D-Ramsey is designed in a way that magnetic shifts will be partially compensated by measuring the COG of both transitions of NV. If both transitions move in the opposite direction in energy space (for example by a magnetic field aligned along the NV spin axis) the COG does not change. This is not always the case. Changing the orientation of an external magnetic field versus the NV axis lead to a modulation of the COG of both lines (see figures D1(a) and (b) for an external field of 10 G). Therefore, for very slow tumbling NDs its rotation may be interpreted as a change in temperature. We define the maximum temperature uncertainty by the taking the maximum delta in the shift of the COG divided by $74.2\,{\rm{kHz}}\,{{\rm{K}}}^{-1}$ and plot it in dependence of the external magnetic field (see figure D1(c)). For an external field of 1 G an uncertainty in temperature reading of 55 mK is reached.

Zoom In Zoom Out Reset image size

As state in the main text we use an agarose matrix to restrict transversal diffusion. Figure E1 illustrates the behavior of the used NDs in the gel matrix. For a single NV one typically utilizes the correlation function ${g}^{2}(\tau )$ to reveal its single photon emission characteristics. The ${g}^{(2)}(\tau )$ function of the rotating ND in the main text is shown in figure E1(a) (blue dots). Interestingly we observe a characteristic change of ${g}^{(2)}(\tau )$ after several hours of temperature read out (figure E1(a) black dots). A new bunching feature occurs, visible in a biasing of the correlation value for long τ to higher values. Analyzing the fluorescence signal using the so-called $G(\tau )$ function, reveals this bunching to have a characteristic timescale of ${\tau }_{l}=4.0\pm 0.2\,{\rm{ms}}.$

Zoom In Zoom Out Reset image size

As the carrier of information in case of photon emission is the fluorescence intensity and its dynamics, we should see similar intensity fluctuation, if we perform a xy-scan with a suitable high temporal resolution. To this end, we perform a fast xy-scan (figure E1(d)). If one acquires a confocal xy-map with sufficient slow scan speed, the size of the confocal spot still corresponds well to the measured average point spread function of other fluorescent spots (compare white and red cycle in figure E1(c)). Therefore, we interpret the occurring intensity fluctuations to be originated from the rotation of the ND, rather than transversal diffusion, which would also lead to fluctuations of the fluorescence signal.

Zoom In Zoom Out Reset image size