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

We present an optimized scheme for nanoscale measurements of temperature in a complex environment using the nitrogen-vacancy center in nanodiamonds. To this end we combine a Ramsey measurement utilized to temperature determination with advanced optimal control theory. We test our new design on single nitrogen-vacancy centers in bulk diamond and fixed nanodiamonds, achieving better readout signal than with common soft or hard microwave control pulses. We demonstrate temperature readout using rotating nanodiamonds in an agarose matrix. Our method opens the way to measure temperature fluctuations in complex biological environment. The used principle however, is universal and not restricted to temperature sensing.


INTRODUCTION
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][5][6][7] probe have been demonstrated. In addition, one can also extract chemical information about samples by nuclear magnetic resonance (NMR) with ~Hz spectral and nanometer-scale spatial resolution. 8,9 These techniques can be modified to measure the temporal dynamics of the quantity under study using NV. 10,11 The sensor capability of NV can be preserved in nanodiamonds (NDs), even in a biological environment, 12 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 in a manifold of variations for NV. [13][14][15][16][17] Thereby one can numerically optimize a pulse to achieve the desired quantum state.
To perform the optimization we use the quantum optimization package DYNAMO, 18 which uses the gradient ascent pulse engineering (GRAPE) algorithm together with Broyden-Fletcher-Goldfarb-Shanno (BFGS) 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, 19 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. 19 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 nanodiamonds. The sequence is sliced into 3 segments (seg1, seg2 and seg3), and each segment is separated equally by /2. and denote quantum phase accumulation due to changes in or the external magnetic field. (B): Before and after seg2 projectors are applied to introduce a cooperative design in the pulse compilation via optimal control. (C): Example of state evolution within the sequence for different drift Hamiltonians and different scenarios. The upper and lower left plots show the state evolution under square pulses with fixed system parameters and a various amplitude and detuned ensemble, respectively. The lower right shows the state evolution under a Coop-D-Ramsey sequence optimized for this ensemble. Thereby, for the lower plots, the solid lines present the average state of the ensemble. The corresponding colored areas are the spreading of state evolutions in terms of one sigma. 0 and 1 are the applied external magnetic field vector and microwave field, respectively.
These two ESR lines can shift due to changes in many internal and external parameters. For example, if the zero field splitting of the NV changes by local temperature changes, 20 they shift in the same direction. A change of the strength of an applied axial magnetic field 0 , 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 = • during the free evolution interval , if the COG is detuned by . Slow changes in 0 , introducing a possible phase shift , are canceled out similar to a Hahn-Echo sequence.
One can separate the D-Ramsey into three pulse segments with individual roles (see Figure 1A and B). 'seg1' initializes quantum sensing, by creating a superposition between |0⟩ and | − 1⟩.
'seg2' prepares the system to be only sensitive to COG shifts of | − 1⟩ and | + 1⟩, which is accomplished by swapping the populations of | − 1⟩ and | + 1⟩.'seg3' converts phase to state population enabling readout.
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: with drift being the free evolution or drift Hamiltonian, while are the control Hamiltonians that describe the interaction of the system with for example an external microwave field. For magnetic fields, is composed of the spin operators , and . The amplitude ( ) at time of this interaction is typically called control field. The task is to find all ( ), such that the system evolves into a predefined final state after time , represented by the density matrix .
Therefore, one defines a quantity 0 that can for example be the overlap of the desired final state with the actual state evolved under the influence of the control fields: 0 = ⟨ | ( )⟩. By maximizing 0 , one optimizes the amplitude ( ) to the designed functionality. Hereby one of the strengths of optimal control lies in the possibility to find a particular set of ( ) that allows robustness against variations of system parameters (e.g. a broad range of drift and ). 21 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 that is coupled to an environment can be described as: Hereby ℒ = −i • ℋ + is the Liouville superoperator, being the sum of the time independent system Hamiltonian ℋ and the dissipator in the time interval Δ = 2 − 1 . The exponential term in eq (2) is typically called a dynamic map . In a different perspective, can also be understood as a projector. By replacing the dynamic map of the system, evolving with ℋ between each segment with a suitable projector, the full sequence can be compiled at once invoking cooperativity. Thereby projects on a subspace in the Liouville space of the density matrices, which in our case is any equal superposition of states |0⟩ and | − 1⟩ after 'seg1' and |0⟩ and +|1⟩ 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 ( ) such that the desired subspace is reached before the projector is applied.
After optimization, we can simply cut the overall sequence between the segments and add arbitrary but equally long free evolution interval in between.
As shown in the following section, 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).

RESULTS AND DISCUSSION
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 referee to as bulk case (see also Supporting Information (SI)). Hereby we compare the Coop-D-Ramsey (cooperatively numerically optimized D-Ramsey) to a typical D-Ramsey consisting of conventional π and π/2 pulses. Figure 1C shows However, the Coop-R-Ramsey converges at the end of each segment to the desired superposition state. To calculate the optimal set of ( ) for the Coop-D-Ramsey, we used the projectors, which are described in eqs (S4a,b) in the SI. We labeled the projectors between 'seg1' and 'seg2' with (1) and (1 ) and between 'seg2' and 'seg3' (2) and (2 ) (see Figure 1B). Note that we have to use two sub-ensembles: one with the set (1) and (2) and the other one with the set (1 ) and (2 ) for a successful compilation of the Coop-D-Ramsey sequence (see discussion in the SI and Figure S1). In addition, we design the system to end up in |0⟩, not in | + 1⟩ as this is the case for the D-Ramsey.
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 frequencies 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 (T1-Decay, see Figure 2D (black) ) as described elsewhere (see also description in Materials and Methods). 22 Then we acquired a Hahn-Echo also to verify that the coherence time of the NV center is long enough ( Figure 2D (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 eqs (S2a,b,c,d). Amplitude variations are accomplished by varying the relative amplitude and detunings by varying ∆ (see also Figure 2C). The parameters for pulse compilation are given in the Materials and Methods part. 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, Δ : 0 to 200kHz). Interestingly, the 7 sequence shows 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 (see Figure 3). First, we compensated the external magnetic field to satisfy the condition ≫ 0 as required to use eqs (S3a,b,c,d) for modeling the NV spin system within a ND. Fitting a continuously driven optical detected magnetic resonance (cwODMR) spectrum ( Figure S3) reveals = 4.35 ± 0.01 MHz and 0 = 0.53 ± 0.03 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 the Materials and Methods part. In contrast to bulk diamond the difference in signal contrast between D-Ramsey and Coop-D-Ramsey is less, but still significant. We attribute this difference to the reduction in the hyperfine splitting introduced by the zero-field parameter : as the hyperfine lines are closer in frequency they are within the driving bandwidth of the conventional square pulses (see also Figure S2).
In a next step, we tested the cooperative pulse design on a tumbling nanodiamond. 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 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). The contrast of the Coop-D-Ramsey is around 80% in the beginning. Over almost two hours the zero field splitting shifted by around 100 kHz. This 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 by 100 kHz, to stay closer to the defined compilation interval (allowed detuning: 0-200kHz), also to rule out that the latter is the origin of the signal drop. A closer look on the acquired Rabi data in Figure 4E 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 23 (compare Figure 4B with C and D ).
Furthermore, we acquired the autocorrelation function (2) ( ) before and after the end of the temperature measurement series (see Figure S5). Interestingly, the (2) ( ) measured after the temperature measurements does not approach one, i.e. the value for a totally uncorrelated signal.
To There are alternative ways to measure temperature with rotating NDs. The simplest measurement scheme is using continuous wave, which, however, is limited by 2 * . 26 In addition, the time resolution of the method is estimated to be on the order of several tens of ms. 12 In case of the D-Ramsey scheme one could also replace conventional pulses, by adiabatic passage. 27 For example, chirped pulses, allowing to create a state inversion with high fidelity and bandwidth, or the 1 insensitive rotation (BIR) pulse, which allows to adjust a superposition state with arbitrary state composition. Here 1 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 1 field amplitude considered. In addition, the BIR pulse scheme can only provide a certain bandwidth in frequency space roughly in the range of the applied Rabi frequency. 27 Especially in bulk diamond and for a NV with 14 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, 19 providing better results for shorter pulses.
In this work we have chosen an overall pulse length of around 2 μ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 μs long. The maximum bandwidth one can access in terms of temporal resolution to measure temperature fluctuations, is thus around 250 kHz. This is comparable to the time resolution reported previously 28 where the resonance frequency at ~840 kHz of a mechanical oscillator was probed. Therefore, it should 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 [29][30][31][32][33] have measured sizeable temperature gradients. Thermodynamic arguments however suggest, 34 that average temperature changes by endogenous thermogenesis should we very small, but do not exclude, that measurable temperature fluctuations may still exist. For example, Inomata et al. 28 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, 35,36 NDs can reach dimensions even below 10 nm and still host NV. 24,25 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.

CONCLUSION
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 nanodiamond is tumbling within an agarose matrix. Ultimately our work paves the way to measure sub 100 mK temperature fluctuations on microsecond time and nanometer length scale.

Sample preparation.
Bulk study: We glue a bulk diamond containing single deep NV on a printed circuit board as sample holder. The latter is also used to apply microwave excitation via a 50 µm thin copper wire spanned over the diamond.
ND study: For fixed nanodiamonds we spin coat a ND sample ("M19-S11c", SN: AA00M7, Diamond Nanotechnologies Inc.) onto a plasma cleaned cover glass glued on a PCB board. A 50 µm thin copper wire is spanned over the glass and connected to the PCB board for microwave excitation.
For tumbling nanodiamonds we modified the above approach by adding a sample chamber on top.
We glued a confocal cell culture dish without the glass slide at the bottom onto the cover glass after spanning the wire. The sample is a 25 μl of ND stock solution ("M19-S11c", SN: AA00M7, Diamond Nanotechnologies Inc.) mixed with 25 mg agarose (target concentration 5 wt%, Agarose fluorescence response of NV is depending on its spin state: The first ~300 ns after turning on the laser give a higher fluorescence signal if being in |0⟩ and lower if being in | ± 1⟩. We perform two spin contrast measurement: the first one is done without any microwave excitation between initialization and readout, giving the reference value for being in |0⟩. In a second run we depopulate |0⟩, using a linear chirp microwave pulse and extract the reference value for being in | ± 1⟩. In the bulk and fixed ND case we determine the spin contrast in a separate measurement. For the rotating nanodiamond study the spin contrast is continuously monitored during the running temperature readout. Further we use DYNAMO to simulate the state evolution of individual NV when varying the free evolution time for a (Coop-)D-Ramsey. As we do not include decoherence effects at this stage, and as we want to compare simulated and experimental results, we also measure a Hahn Echo for the individual NVs. The Hahn Echo is then fit to a stretched exponential function, which parameters are used as an exponential damping term to fit experimental data. In case of the rotating ND we keep the parameters of the stretched exponential as a free fit parameter.
To extract the correlation time of rotating nanodiamond we use a mono exponential decay. Figure   illustrations have been accomplished using Inkscape.

Description of spin system for pulse compilation.
To compile a Coop-D-Ramsey sequence, 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 nanodiamonds. We define drift and c = ∑ ( ) (see main text) as: and is the zero field splitting, is the zero field component introduced by off-axial strain, the spin operators of a spin one system, 0 static components introduced by an external magnetic field with strength 0 = √∑ 0 2 , and 1 the components of the microwave field.
We distinguish two cases. The first case is a NV with | | ≪ 0 , when the external static field shall be aligned along the -axis: | 0 | = 0 . We call this bulk case. And second, an NV with | | ≫ 0 , assigned as the ND case. The microwave field is chosen to be an in-plane linear polarized field with angle to the -axis ( 1 ( ) = 1 cos( + ) , 1 ( ) = 1 sin( + ) and 1 ( ) = 0 ). Transforming eqs (S1a) and (S1b), into a left-handed and a right-handed rotating reference frame, leads to the overall system Hamiltonian for the bulk case using the assumption In real experimental situation 1 will be the measured Rabi frequency that is determined before pulse compilation.
In the ND case the situation becomes more complicated:

= drift + ( ) • [cos( ) ∥ + sin( ) ⊥ ] + ( ) • [cos( ) ∥ + sin( ) ⊥ ]
with: and: Again, for parameter variation eq (S3c) is redefined as: ). The interpretation of eq (S3c) can also 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 perpendicular aligned, one can only drive the other transition.
Projectors to for cooperative pulse design. To introduce cooperativity, we use the projector sets To demonstrate that it is mandatory to use both projector sets, we compile one sequence with (1) and (2) only, and one including all projectors.
The system was defined as described in eq (S2a,b,c). The shift introduced by an external zmagnetic 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 diff erent configurations that resample the hyperfine structure of 14 Figure S1A). When using set1 and set2, one can observe a clear oscillation between state |0⟩ and | + 1⟩ with a period of 200 kHz for all hyperfine lines ( Figure S1B). The additional frequency components between 1 MHz and 4 MHz disappeared in the corresponding Fourier transformation (compare Figure S1C and D). Only Projector (1) and (2) are used to compile the pulse. (B): Compilation includes an additional sub-ensemble applying (1 ) and (2 ) instead of Projector (1) and (2) .  Then, we calculate the Eigensystem of NV and determine the most prominent transitions introduced by the spin operators and in the transformed system. Figure S2 shows this calculation for varying 0 for a zero field splitting =5 MHz and 14 N. As one can see the hyperfine lines of the = ±1 state overlap for 0 ≪ . With increasing magnetic field the lines begin to separate. For 0 ≥~2 the hyperfine lines are well separated. In Figure S2B we plot the corresponding transitions in case of 15 N.
As an example Figure S3 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 equals 4.35 ± 0.01 MHz. 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 Figure S4A and S4B for an external field of 10 G). For very slow tumbling NDs its rotation may therefore 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 kHz K −1 and plot it in dependence of the external magnetic Analysis of rotating ND. As state in the main text we use an agarose matrix to restrict transversal diffusion. Figure S5 illustrates the behavior of the used NDs in the gel matrix. For a single NV one typically utilizes the correlation function 2 ( ) to reveal its single photon emission characteristics.
The (2) ( ) function of the rotating ND in the main text is shown in Figure S5A (blue dots).
Interestingly we observe a characteristic change of (2) ( ) after several hours of temperature read out ( Figure S5A 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 ( ) function, reveals this bunching to have a characteristic time scale of = 4.0 ± 0.2 ms. 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 S5D). 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 S5C). 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.
Simulating the D-Ramsey or a tumbling nanodiamond.