Nanoscale magnetometry through quantum control of nitrogen-vacancy centres in rotationally diffusing nanodiamonds

The confluence of quantum physics and biology is driving a new generation of quantum-based sensing and imaging technology capable of harnessing the power of quantum effects to provide tools to understand the fundamental processes of life. One of the most promising systems in this area is the nitrogen-vacancy centre in diamond - a natural spin qubit which remarkably has all the right attributes for nanoscale sensing in ambient biological conditions. Typically the nitrogen-vacancy qubits are fixed in tightly controlled/isolated experimental conditions. In this work quantum control principles of nitrogen-vacancy magnetometry are developed for a randomly diffusing diamond nanocrystal. We find that the accumulation of geometric phases, due to the rotation of the nanodiamond plays a crucial role in the application of a diffusing nanodiamond as a bio-label and magnetometer. Specifically, we show that a freely diffusing nanodiamond can offer real-time information about local magnetic fields and its own rotational behaviour, beyond continuous optically detected magnetic resonance monitoring, in parallel with operation as a fluorescent biomarker.


Introduction
The diamond nitrogen-vacancy (NV) defect shows great promise as a biological imaging tool in two distinct capacities. Firstly, its optical properties and biocompatibility make it an excellent fluorescent biolabel [1][2][3][4][5][6][7][8]. It is extremely bright, highly photostable [9], and it can be imaged with spectacular spatial resolution [10] through sub-diffraction localization techniques like stimulated emission depletion [11]. Secondly, optical spin polarization and measurement and microwave control, the NV defect's electron spin can be used to detect magnetic fields with nT Hz −1/2 [12,13] sensitivity, with potentially nanometre spatial resolution [14,15]. The prospect of using the NV centre for nanoscopic MRI has significant promise for a range of problems in biology, including: detecting operation of individual neuronal ion channels [16] and atom-scale magnetic resonance imaging [17,18].
An interesting application is the combination of the NV as a nano-scale biomarker and magnetometer. This would be a super-optical probe capable of reporting motional information simultaneously with atomic level magnetic fields and their fluctuations over arbitrary long times. Such a system would be able to give the precise location of some attached ligand as it moves through a cell, as well as information about its molecular environment by detecting nanoscale magnetic fields. However, nearly all demonstrations of NV magnetometry have assumed the NV spin quantization axis is effectively stationary with respect to the lab-frame control fields. In order to employ the NV centre as a combined nanoscale bio-marker/magnetometer one must fully understand the effect of random rotational motion on the nanomagnetometry and control protocols. In [19] NV nanodiamond nanomagnetometry in living cells was demonstrated. In these experiments, the measurements were conducted on a endocytosed nanodiamond which moderated the motional timescales, facilitating quantum measurement and allowing information about the orientation to be determined over long periods. In this paper we explore the quantum The magnetic field direction of the Rabi pulse defines the z-axis. z is the instantaneous direction of the NV axis, defined with respect to the lab-frame unprimed coordinate system by θ , φ.
control of a rotationally diffusing NV system and find that the accumulation of geometric phases, due to rotation [20][21][22][23][24][25], is the critical phenomenon in the application of the NV centre as a bio-marker/magnetometer.

The diamond nitrogen-vacancy system
The NV defect has a spin triplet ground state with a 2.88 GHz zero-field splitting between the m = 0 state and the degenerate m = ±1 states as shown in figure 1(a). Optical excitation at 532 nm can pump the defect into the m = 0 state, and allows the population of the ground state to be read, since the m = 0 state has a higher fluorescence than the m = ±1 states. The effective Hamiltonian of the NV ground state triplet in the presence of a magnetic field B (lab frame coordinate system defined by z, see figure 1(b)), ignoring crystal asymmetries and hyperfine effects, is where γ is the gyro-magnetic ratio of the NV (γ = gµ B /h, g ≈ 2). The first term is the zerofield splitting of the NV system itself, where D = 2.88 GHz is the zero-field splitting strength. It is this term which makes the crystal's orientation crucial. It defines a quantization direction z (see figure 1(b)) which lies along the axis connecting the nitrogen atom to its adjacent vacancy. The second term is the usual Zeeman splitting interaction with a magnetic field. Previously it has been shown that the NV centre can be used to measure dc and ac magnetic fields [12,14] and fluctuating (fc) fields [26,27] and rotation [19][20][21][23][24][25]. These protocols work by optically polarizing the NV m = 0 spin state. A π/2 microwave pulse then establishes a coherent superposition of the m = 0 and 1 states, and the system evolves under the influence of any magnetic fields by the Zeeman effect. After an evolution time, τ , the rotation and strength of the magnetic field manifests itself as a phase, relative to m = 0, of the m = 1 state. In the adiabatic limit this phase is where G is the accumulated geometric phase, due to rotation along some path, P, defined by θ and φ (see figure 1(b)) [21,23] and D is the accumulated dynamical phase due to the change in the projection of the magnetic field, B onto the NV axis,ẑ . A second π/2 pulse converts this phase into a population difference, which is read out directly as a change in fluorescence. It is common to apply a static magnetic field to split the m = ±1 states, selecting only one of the m = ±1 states by using a sufficiently weak microwave field resonant with one of the transitions (Rabi field). Since a static field would cause different splitting depending on the crystal's orientation [19], we do not take this approach. Instead, we assume that the degeneracy of the m = ±1 states is lifted by strain, which is usually the case for nanodiamonds [19]. Specifically, we consider the regime where z + E 2 ⊥ and E ⊥ characterizes the effective non-axial electric-strain field components [28,29]. In this regime, the field spin states are approximate eigenstates of equation (1) with energies E 0 = 0 and E ± = D±R.

Motion of the crystal
A spherical crystal is considered, whose rotation follows Brownian motion, according to the Langevin equation where I is the crystal's moment of intertia, ω is the crystal's angular velocity, τ is a stochastic delta function correlated torque and γ d is a drag coefficient given by Stokes' law, γ d = 8πηr 3 , where r is the crystal's radius and η is the fluid's viscosity. It can be imagined that the crystal rotates in a 'random walk' of discrete small angle steps, whose directions are uncorrelated with one another. This rotational motion is characterized by a single diffusion coefficient k d which, from equation (3) and the equipartition theorem, is given by where k B T is the thermal energy. Thus, in room-temperature water, r = 5 nm (r = 50 nm) gives k d ≈ 1 rad 2 µs −1 (k d ≈ 1 rad 2 ms −1 ), see red dashed curve in figure 2.

Ensemble averaged measurements
In practice, and even in theory, a single measurement, of a quantum system, contains very little information. The quantity which is usually measured is a time averaged signal from a large ensemble of individual measurement runs. Environmental perturbations cause subtly different quantum evolution in each member of the ensemble, leading to a loss of information in an ensemble-averaged signal after a certain evolution time. In the case of the rotationally diffusing nanodiamond, the ensemble's members differ more dramatically. In addition to its own microscopic environment each member of the ensemble has its own starting orientation and rotational trajectory. A series of sequential measurements is made which may involve several NV centres in a single crystal. If the time over which measurements are made is sufficiently small compared with Also shown are the typical bulk timescales associated with the NV spin state: homogeneous broadening (T 2 = 2 ms), inhomogeneous broadening (T * 2 = 10 µs) and zero field splitting (1/D). For nanodiamonds T 2 and T * 2 can be significantly smaller than the bulk values. For comparison we plot (labelled squares) observed values of T 2 [30] (A), [12] (B) and [19,31,32] (C). the rotational diffusion time, 1/k d , then the situation would be comparable to the case of a static crystal. In this work the opposite limit is considered and it is assumed that sufficiently many measurements are made that the final signal is an average over the full theoretical ensemble: all initial orientations and rotational trajectories are accounted for. Experimentally this average can be constructed via the repeated measurement of the same NV centre or over a collection NV centres [19].
For such a scenario consider a continuously applied microwave field, of amplitude B R , tuned to the ψ (m=0) transition. This will cause Rabi oscillations with a frequency which depends on the angle, θ , between the NV axis z and z, the oscillation direction of the microwave magnetic field. Initialized into ψ (0) z , the population of ψ (0) z for a single NV as a function of time is then where R = gµ B B R /h is the Rabi frequency when θ = π/2. Assuming R 1/k d the ensemble-averaged signal is an average over all the possible different Rabi frequencies: where H α is the Struve function. S 0 (t) (solid curve) is shown in figure 3(a), along with the equivalent curve for a static crystal (dashed curve), with θ = π/2. For an ensemble averaged measurement the signal reaches its first minimum at t R ≈ 1.16π at which point S 0 ≈ 0.14. This result enables the determination of the optimal microwave pulse time to achieve a π/2 rotation. The required time is simply read off from the curve. In our system, there is a certain amount of freedom in deciding how long a π/2 pulse should be. We could take it to be t π/2 = π/(2 R ), which would produce a π/2 rotation for an optimally aligned (θ = π/2) NV centre, and a smaller rotation otherwise. Alternatively, we could use t π/2 = 1.16π/(2 R ) which, based on figure 3(a), would optimize the initial amplitude of a Ramsey-type experiment. One of the strengths of Ramsey-type experiments is their robustness to the π/2 pulse time chosen. We do not, therefore, expect this decision to be an important one (as simulations will confirm). For now we keep our discussion general, introducing a parameter, a, to describe the pulse time. As such when we refer to a 'π/2' pulse we mean a pulse of length t π/2 = aπ/(2 R ), where a ≈ 1, which produces a spin rotation of = a(π/2)sin θ . Now consider a Ramsey-type pulse sequence. The NV centre, starting with z at some polar angle θ 1 to the z-axis, is first optically pumped into ψ (0) z . A 'π/2' pulse is then applied, tuned to the ψ (0) z → ψ (1) z transition, producing a spin rotation of 1 = a(π/2)sin θ 1 . The system is allowed to evolve for time τ , during which time the crystal rotates through some trajectory and the ψ (1) z state develops a phase , due to rotation and the presence of magnetic fields. A final 'π/2' pulse is applied. This final pulse may be in phase with the initial pulse, or out of phase by π. It produces a rotation of ± 2 ≡ ±a(π/2)sin θ 2 , where θ 2 is the final angle between z and z and the + (−) refers to the pulse being in (out) of phase with the original pulse. The final population P 0 of the ψ (0) z state can then be written down explicitly as The final, normalized, signal is the ensemble average of this value S(t) = P 0 , where the ensemble average is taken over all possible starting orientations and rotational trajectories. The upper and lower bounds of the signal envelope are produced by the two curves implied by the ± term of equation (7), provided the mean phase vanishes, = 0 (true for the cases considered in this work). This work primarily focuses on the behaviour (typically decay) of the signal envelope, due to rotation and the presence of magnetic fields, rather than the oscillations themselves. Physically, the two curves correspond to the phase (0 or π respectively) of the final 'π/2' pulse. An example of a signal envelop, due to dephasing from geometric phase accumulation (calculated below), is shown in figure 3(b).

Dephasing due to geometric phase accumulation
To see the influence of the geometric phase, we consider how the spin state of the rotating crystal evolves in the absence of external fields. As the crystal rotates, the direction of z changes, so that equation (1) becomes a time-varying Hamiltonian. With respect to the eigenstates of spin projection along the fixed z-axis, the zero-field Hamiltonian is where θ and φ are time dependent. If the rate of rotation is small compared with 2.88 GHz the evolution will be adiabatic, producing no changes in populations of the spin sublevels (with respect to the z quantization axis). Each sublevel, however, develops a geometric phase There is a gauge degree of freedom in defining the phase of the eigenstates [21,23]. We choose to define the eigenstates with respect to the fixed z basis as This convenient choice of gauge means that an explicit phase factor between sublevels corresponds to the absolute phase of a microwave field linearly polarized along the z-axis. The relative geometric phase accumulation, between ψ (0) z and ψ (1) z , due to the crystal's rotation then evolves according to˙ =φ cos θ , where the phase evolution of each state is given bẏ m = m˙ [20,21,[23][24][25].
Given a particular initial orientation θ 1 the probability density of a particular final orientation θ 2 and geometric phase after time t is p(θ 2 , ; t|θ 1 ). The ensemble-averaged signal can then be expressed as the integral of the final ground state population over all possible θ 1 , θ 2 and , weighted by the probability distribution function and the probability of the NV starting with orientation θ 1 : where P 0 (θ 1 , θ 2 , ) is given by equation (7). The Fokker-Planck equation for the evolution of p (θ 2 , ; t|θ 1 ) is The first term describes the probabilistic rotation of the crystal, and the second term describes phase evolution, which deterministically depends on the probabilistic evolution of azimuthal angle φ. The initial condition is and the boundary conditions are periodic in θ and . Numerically solving equation (11) gives the solution shown in figure 3(b), demonstrating that the decaying envelope is relatively insensitive to the 'π/2' pulse duration, a. There is approximately exponential decay of the signal envelope over the timescale of the crystal's rotation, with decay time ≈0.90/k d . The diffusion rate, k d , thus sets an upper bound on the coherence time for the crystal. In the absence of any external fields, measurement of this decoherence rate allows measurement of k d with shot noise-limited sensitivity. Specifically, where τ d is the overall decoherence time, α ≈ 0.01 is a parameter which accounts for the non-zero fluorescence of the m = ±1 states, and finite collection efficiency [17], N is the number of NV centres in the crystal, along all four crystallographic axes and T is the total measurement time. Assuming a linear decoherence response, 1/τ d ≈ 1/T * 2 + k d , figure 4 (red dashed curve) shows the sensitivity as a function of crystal radius, with an assumed NV density of 10 16 cm −3 and a bulk T * 2 value of 10 µs. In assuming a linear decoherence response cooperative decoherence effects have been ignored and it has been assumed that T * 2 does not change with rotation. In figure 2 measured values (black squares) for T 2 are plotted for several cases [12,19,[30][31][32]. Assuming that for these systems T * 2 ≈ T 2 , in figure 4, we plot the sensitivity for these specific cases (red squares).

Dephasing due to population mixing with non-adiabatic evolution
Up to this point we have been working within the adiabatic approximation. In the context of the NV centre, the adiabatic criterion is 1 2 θ 2 +φ 2 sin 2 θ D 2 .
The angular velocity of the nanodiamond is a consequence of thermal motion. Using the equipartition theorem, the adiabatic criterion for the rotationally diffusing nanodiamond is where ρ is the density of diamond. At room temperature this criterion will be satisfied if the crystal radius is much greater than 0.2 nm (see figure 2, dashed red curve). An additional adiabatic criterion, often overlooked, is that the state of the NV centre changes slowly compared to 1/D. The timescale on which the NV state changes is the angular velocity damping time, t d = I /γ d = ρr 2 /(15η). In figure 2(a) t d (green dashed-dotted curve) is plotted as a function of the radius of the crystal, for a crystal in room temperature water. From this it can be seen that t d is smaller than 1/D when the crystal radius is smaller than 20 nm.
The geometric phase accumulation comes from a zeroth order perturbative expansion of the time evolution operator in the ratio of the rotation rate to the 2.88 GHz zero field splitting. Taking the expansion to higher orders enables the population mixing that comes from nearly adiabatic evolution to be evaluated. To evaluate the properties of an ensemble measurement consider an arbitrary initial state of the NV centre, where each state has a probability weighting of P m and phase m . As shown in the appendix the ensemble average population mixing rate is The above provides an expression for the population mixing rate in terms of physical parameters: the zero-field splitting frequency D, the rotational diffusion rate k d and the angular velocity damping time t d . To obtain the above results several assumptions have been made. Most significantly, the diffusive rotation of the nanodiamond has been modelled as a series of discrete steps, using physical arguments to infer the appropriate statistical behaviour of the steps' sizes and durations.
To check the validity of these approximations a numerical Monte Carlo simulation has also been employed. The simulation models the Langevin equation (3), using a torque which fluctuates on a finite but very small timescale, to produce a set of rotational trajectories. For each trajectory the Schrödinger equation is solved numerically. The results, for an NV centre which starts in where t m = 1/(3 k ) is the population mixing time.
The implication of this non-adiabatic behaviour is more than merely a correction to the behaviour expected from geometric phase accumulation. Measuring the population mixing rate (effectively measuring T 1 ) provides another way to gauge the rotational diffusion rate. More over the measurement can be achieved without microwave fields. It is enough to simply initialize the NV into the ψ (0) z state and observe the decay of the fluorescence rate as a function of free evolution time. The blue dotted curves in figures 2 and 4 show the population mixing time, and associated sensitivity, respectively, as a function of the crystal radii.

Static and fluctuating magnetic fields
We now investigate the properties of rotationally diffusing NV centres in the presence of static (dc), fluctuating (fc) and oscillating (ac) magnetic fields. Magnetic fields enter our model by advancing the phase of the NV centre's spin sublevels. In the regime where B ⊥ = B 2 x + B 2 y D and B 2 ⊥ /D R [28,29], under the secular approximation equation (1) becomes it is only the field component parallel to z (the NV axis) which affects the phase.

Static magnetic fields
Since the orientation of z changes randomly with time, a static field effectively behaves as a field which fluctuates on the timescale of rotational diffusion 1/k d [27]. To demonstrate this consider a static magnetic field with components B z and B x along the zand x-axis (without loss of generality, B y = 0). Unlike the geometric phase, the phase evolution due to a magnetic field depends explicitly on the azimuthal angle φ as well as the polar angle θ of z . The Fokker-Planck equation is then The resulting signal envelope (equation (10)) is plotted in figure 5(b) and the decoherence time (defined here as the time taken for the signal envelope to reach 1/e of its initial value) is plotted in figure 5(c). When γ B k d , the shape of the dephasing signal is close to Gaussian, the characteristic shape for slowly fluctuating fields. In this regime the static field effectively fluctuates on the timescales of 1/k d due to the rotation of the crystal. The crossover between Gaussian and exponential decay occurs on the timescale of 1/k d , but the exponential behaviour is hidden since the geometric phase accumulation causes dephasing to occur on at least this timescale. The decay is significantly stronger for fields in the x (or y) directions than in the z direction. This is because the largest Rabi rotations are experienced by the ensemble members for whom θ ≈ π/2 which are most sensitive to fields in the x-y plane.
For weak fields, the sensitivity approaches zero, which can be seen in the vanishing gradient for small B in figure 5(c). To sense small dc fields, a strong static magnetic field could be applied to bring the field sensitivity to within an appropriate range. The nonlinearity of field sensitivity also allows the possibility of vector magnetometry. By applying a dc field in one direction, only dc fields in that direction would be detected, allowing both the magnitude and the direction of dc magnetic fields to be imaged. Given the optimal strength of applied dc field (γ B ≈ k d ), the magnetic field sensitivity is where 1/τ d ≈ 1/T * 2 + k d . The solid black curve/squares in figure 4 shows the dc magnetic field sensitivity as a function of the crystal radii in room temperature water.

Fluctuating magnetic fields
Many biological systems display behaviour that is neither static nor periodic, but stochastically fluctuating. Although the details may vary, a fluctuating field can be characterized by a correlation time t c , and a mean squared field strength B 2 [27]. For the tumbling NV system, a magnetic field correlation time larger than the rotational diffusion time will mean that the NV experiences a field that is essentially static over its coherence time, although statistically distributed over many runs. In this case, the effect will be very similar to that for a dc field.
For the case of a fluctuating field with correlation time t c 1/k d and mean squared field strength B 2 = B 2 fc 1/(γ t c ) 2 which is isotropic, the Fokker-Planck equation is which is similar to that for the geometric phase alone (11). As such the fluctuating field produces exponential decay of the signal, on a timescale of 1/(γ B fc √ t c ) 2 . The implication of this result is that the decoherence rate is linear in B 2 fc t c .

Oscillating magnetic fields
One application of a nanoscale magnetometer is to detect magnetic resonance from a small number of spins, which can be achieved by measuring the ac field produced by their effective precession, or by driving them via an external microwave field. Applying a spin-echo pulse sequence ( π 2 − τ 2 − π − τ 2 − π 2 ) effectively rectifies an ac magnetic field of period 2τ , allowing its amplitude to be detected while masking fields oscillating at different frequencies. A spinecho pulse sequence has the additional advantage of dramatically increasing the dephasing time of an NV magnetometer, by refocusing the divergent phases accrued due to slowly fluctuating magnetic fields.
In detecting magnetic resonance, a strong dc field is usually applied. A tumbling crystal will see such a field as if it were fluctuating on the timescale of the crystal's rotation, causing dephasing on a timescale, τ nπ , given by where n is the number of π -pulses employed in a concatenated spin-echo pulse sequence. Taking the limit in which very little rotation occurs on the timescale of free evolution. The problem is then a matter of simply averaging over the possible crystal orientations and initial phases of the ac field, rather than solving a Fokker-Planck equation. Consider a spin-echo pulse sequence during which the crystal is oriented at some θ and the free evolution time τ matches the ac fields period. The final population of ψ (0) z is given by P 0 (θ, field ) = 1 2 + 1 2 cos 2 1 cos 2 2 ± 1 4 sin 2 1 (1 − cos 2 2 ) cos(2 ac ), where 1 = a(π/2)sin θ and 2 = aπ sin θ are the angles through which the spin is rotated during the π/2 pulses and π -pulse respectively. For an ac field oscillating along the z-axis, the phase acquired by the NV centre before the π pulse is ac = (τ/2)γ B ac cos θ cos field , where field is the phase of the ac field at the start of the period, and the phase acquired after the π pulse is − ac .
In the limit that the strength of the ac field to be detected is much smaller than the applied dc field strength, it can be assumed that the crystal is stationary over the timescale of free evolution, although an average over the full ensemble of possible orientations still must be performed. To find the ensemble-averaged signal, we integrate over all phases of the ac field, field , and all orientations, θ , weighted by 1 2 sin θ: Figure 5(d) shows the signal envelope due to an ac field. As with dc fields, the optimum ac field strength is such that the dephasing due to the ac field matches the dephasing rate from other sources. The optimal sensitivity to an ac field is the same as that given in (19), except with 1/τ d ≈ 1/T 2 + 1/τ nπ + k d . Using a concatenated pulse sequence can increase the decoherence time. Concatenation also reduces the bandwidth of ac field detection, which is important if different ac frequency components are to be distinguished.

Conclusion
We have shown that a freely diffusing diamond nanocrystal can act as a sensitive magnetometer with sensitivity comparable to that of a fixed crystal magnetometer. We have also shown that, through the accumulation of a geometric phase, the fluorescence signal contains information about the crystal's rotational diffusion rate. This information can be gathered without applying any microwave control sequences, and is thus achievable with little modification to the existing experimental setups that use fluorescent colloidal nanocrystals as biological markers. Before the tumbling nanocrystal system can be used in biological applications, further theoretical and experimental work will be required. It will help to understand anisotropic rotation, and to investigate the effect of the fields produced by a specific biological process rather than the idealized fields considered here. Additionally, in this work the effects of rotation on T 2 and T * 2 have been ignored and it has been assumed that the splitting of the m = ±1 states is small compared to the hyperfine splitting. To provide a quantitative comparison with experiment further work will be required to examine both of these effects, particularly for nanodiamonds with a radius less than 100 nm. The most challenging experimental task will be finding protocols to either track or control rapid spatial movement of the crystal, while still achieving satisfactory microwave control of the quantum system within [19].
of the spin sublevels after some time interval δt is P 1 = P 1 − kδt P 1 + kδt P 0 + 2 P 1 P 0 √ kδt cos ( 0 − 1 + δφ) , (A.1) P 0 = P 0 − 2kδt P 0 + kδt P 1 + kδt P −1 − 2 P 1 P 0 √ kδt cos ( 0 − 1 + δφ) To find the ensemble average of k we need to know how δt is distributed. If we multiply the Langevin equation (3) by ω(0) and take the ensemble average we find that the correlation function of the angular velocity decays exponentially over time ω(t) · ω(0) = ω(0) · ω(0) e −t/t d , (A.8) where we have defined the damping time, t d = I /γ d . In order to retain this correlation function in our discrete random walk approximation, the duration of a given step, δt, must be exponentially distributed. At any given instant, the fraction of the ensemble that is experiencing a step of duration δt is then given by the probability density We do not need to know the exact distribution of angular velocities, merely that (δθ) 2 /(δt) 2 = 2/3 ω · ω = 2k d /t d . The ensemble average population mixing rate is then given by