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

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–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 control of a rotationally diffusing NV system and find that the accumulation of geometric phases, due to rotation [20–25], is the critical phenomenon in the application of the NV centre as a bio-marker/magnetometer.

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

$$\begin{equation} H = \frac{1}{\hbar}D S_{z'}^2 +\gamma \boldsymbol{B}\cdot \boldsymbol{S}, \end{equation} \noindent \tag{ 1 }$$

where γ is the gyro-magnetic ratio of the NV (γ = gμ_B/ℏ, g ≈ 2). The first term is the zero-field 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.

Zoom In Zoom Out Reset image size

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–21, 23–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

$$\begin{eqnarray} &&\Phi=\Phi^{\mathrm{G}}+\Phi^{\mathrm{D}}=\int_P \cos \theta\, \mathrm{d}\phi+\gamma\int_0^{\tau} {\hat z}^{\prime}(t)\cdot {\bf B}(t) \mathrm{d}t, \nonumber \\ \end{eqnarray} \tag{ 2 }$$

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, ${\hat z}^{\prime }$. 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 $B_{\perp ^{\prime }}=\sqrt {B^2_{x^{\prime }}+B^2_{y^{\prime }}} {\ll } D$ and $B_{{\perp }^{\prime }}^2/D {\ll } \cal {R}$ where ${\cal R}=\sqrt {B_{z^{\prime }}^2+{\cal E}_{\perp ^{\prime }}^2}$ and $\cal {E}_{\perp ^{\prime }}$ 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 and $E_{\pm }=D{\pm } \cal {R}$.

3.1. Motion of the crystal

A spherical crystal is considered, whose rotation follows Brownian motion, according to the Langevin equation

$$\begin{equation} I\dot{\boldsymbol{\omega}} = -\gamma_{\mathrm{d}} \boldsymbol{\omega} + \boldsymbol{\tau}(t), \end{equation} \tag{ 3 }$$

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³, 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

$$\begin{equation} k_{\mathrm{d}} = \frac{k_{\mathrm{B}} T}{8 \pi r^3 \eta}, \end{equation} \tag{ 4 }$$

where k_BT is the thermal energy. Thus, in room-temperature water, r = 5 nm (r = 50 nm) gives k_d ≈ 1 rad²μs⁻¹ (k_d ≈ 1 rad² ms⁻¹), see red dashed curve in figure 2.

Zoom In Zoom Out Reset image size

3.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 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)_z' → ψ^(m=1)_z' 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 ψ⁽⁰⁾_z', the population of ψ⁽⁰⁾_z' for a single NV as a function of time is then

$$\begin{eqnarray} &&P_0 = [1+\cos\,(\Omega_{\mathrm{R}}t \sin \theta)]/2, \end{eqnarray} \tag{ 5 }$$

where Ω_R = gμ_BB_R/ℏ 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:

$$\begin{eqnarray} S_0(t)&=&\frac{1}{4}\int_0^{\pi} \mathrm{d}\theta \sin \theta \left[1+\cos \left(\Omega_{\mathrm{R}} t \sin \theta\right)\right] \nonumber \\ &=& \frac{1}{4}\left[2+\pi H_{-1} \left(\Omega_{\mathrm{R}}t\right)\right], \end{eqnarray} \tag{ 6 }$$

where H_α is the Struve function. S₀(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.14.

Zoom In Zoom Out Reset image size

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 θ₁ to the z-axis, is first optically pumped into ψ⁽⁰⁾_z'. A 'π/2' pulse is then applied, tuned to the ψ⁽⁰⁾_z' → ψ⁽¹⁾_z' transition, producing a spin rotation of Θ₁ = a(π/2)sin θ₁. The system is allowed to evolve for time τ, during which time the crystal rotates through some trajectory and the ψ⁽¹⁾_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 ±Θ₂ ≡ ± a(π/2)sin θ₂, where θ₂ 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₀ of the ψ⁽⁰⁾_z' state can then be written down explicitly as

$$\begin{eqnarray} &&P_0=\frac{1}{2}\left[1+\cos \Theta_1 \cos \Theta_2 \pm \sin\Theta_1 \sin \Theta_2 \cos \Phi\right]. \end{eqnarray} \noindent \tag{ 7 }$$

The final, normalized, signal is the ensemble average of this value S(t) = 〈P₀〉, 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).

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

$$\begin{eqnarray} H_0& =&\frac{1}{\hbar}D S_{z'}^2 \nonumber \\ &= &D\hbar\left( \begin{array}{@{}ccc@{}} \displaystyle \cos^2\theta +\frac{\sin^2\theta}{2} & \displaystyle \frac{\mathrm{e}^{-\mathrm{i}\phi}\cos\theta\sin\theta}{\sqrt{2}} & \displaystyle \frac{\mathrm{e}^{-2\mathrm{i}\phi}\sin^2 \theta}{2} \\ \displaystyle\frac{\mathrm{e}^{\mathrm{i}\phi}\cos\theta\sin\theta}{\sqrt{2}} & \sin^2\theta & -\displaystyle \frac{\mathrm{e}^{-\mathrm{i}\phi}\cos\theta\sin\theta}{\sqrt{2}} \\ \displaystyle\frac{\mathrm{e}^{\mathrm{2i}\phi}\sin^2 \theta}{2} & - \displaystyle \frac{\mathrm{e}^{\mathrm{i}\phi}\cos\theta\sin\theta}{\sqrt{2}} & \cos^2\theta +\displaystyle \frac{\sin^2\theta}{2} \end{array}\right)_z\!, \end{eqnarray} \noindent \tag{ 8 }$$

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

$$\begin{equation} \Phi_m=\int_0^t \left[\frac{\mathrm{d} \psi_{z^{\prime}}^{(m)}(t)}{\mathrm{d}t} \right]^{*} \psi_{z^{\prime}}^{(m)}(t)\mathrm{ d}t. \end{equation} \noindent \tag{ 9 }$$

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

$$\begin{eqnarray*} \psi_{z^{\prime}}^{(1)}&=& \left(\begin{array}{@{}c@{}} \displaystyle\mathrm{e}^{-\mathrm{i}\phi}\cos^2\frac{\theta}{2}\\ \displaystyle\frac{1}{\sqrt{2}} \sin\theta\\ \displaystyle\mathrm{e}^{\mathrm{i}\phi}\sin^2\frac{\theta}{2}\end{array}\right)_z\!, \quad \psi_{z^{\prime}}^{(0)}= \left(\begin{array}{@{}c@{}} \displaystyle-\frac{1}{\sqrt{2}} \mathrm{e}^{-\mathrm{i}\phi}\sin\theta\\ \cos\theta \\ \displaystyle\frac{1}{\sqrt{2}} \mathrm{e}^{\mathrm{i}\phi}\sin\theta\end{array}\right)_z\!,\\ \psi_{z^{\prime}}^{(-1)}&=& \left(\begin{array}{@{}c@{}} \displaystyle\mathrm{e}^{-\mathrm{i}\phi}\sin^2\frac{\theta}{2}\\ \displaystyle-\frac{1}{\sqrt{2}} \sin\theta\\ \displaystyle\mathrm{e}^{\mathrm{i}\phi}\cos^2\frac{\theta}{2}\end{array}\right)_z\!. \end{eqnarray*} \noindent$$

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 ψ⁽⁰⁾_z' and ψ⁽¹⁾_z', due to the crystal's rotation then evolves according to $\dot \Phi = \dot \phi \cos \theta$, where the phase evolution of each state is given by ${\dot \Phi }_m=m{\dot \Phi }$ [20, 21, 23–25].

Given a particular initial orientation θ₁ the probability density of a particular final orientation θ₂ and geometric phase Φ after time t is p(θ₂,Φ;t|θ₁). The ensemble-averaged signal can then be expressed as the integral of the final ground state population over all possible θ₁, θ₂ and Φ, weighted by the probability distribution function and the probability of the NV starting with orientation θ₁:

$$\begin{eqnarray} S(t)&=&\frac{1}{2}\int \mathrm{d} \theta_1\, \mathrm{d} \theta_2 \,\mathrm{d} \Phi \sin \theta_1 p\left(\theta_2 , \Phi; t |\theta_1\right) P_0\left(\theta_1,\theta_2,\Phi \right) , \nonumber \\ & & \end{eqnarray} \tag{ 10 }$$

where P₀(θ₁,θ₂,Φ) is given by equation (7). The Fokker–Planck equation for the evolution of $p\left (\theta _2 , \Phi ; t |\theta _1\right )$ is

$$\begin{eqnarray} &&\fl \frac{1}{k_{\mathrm{d}}} \frac{\mathrm{d}}{\mathrm{d}t} p(\theta_2,\Phi,t|\theta_1) =- \frac{\partial}{\partial \theta_2}\left [ \sin \theta_2 \frac{\partial}{\partial \theta_2}\left ( \frac{p\left(\theta_2 , \Phi; t |\theta_1\right)}{\sin\theta_2}\right )\right ] + \cot^2\theta_2 \frac{\partial^2}{\partial \Phi^2}p\left(\theta_2 , \Phi; t|\theta_1\right).\nonumber\\ \end{eqnarray} \tag{ 11 }$$

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

$$\begin{equation} p(\theta_2,\Phi,t=0|\theta_1) = \delta(\Phi)\delta(\theta_2 - \theta_1)\end{equation} \tag{ 12 }$$

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, $\Delta k_{\mathrm {d}} \sqrt T {\approx } \frac {1}{\alpha \sqrt {N \tau _{\mathrm {d}}}}$, 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*₂ + k_d, figure 4 (red dashed curve) shows the sensitivity as a function of crystal radius, with an assumed NV density of 10¹⁶ cm⁻³ and a bulk T*₂ value of 10 μs. In assuming a linear decoherence response cooperative decoherence effects have been ignored and it has been assumed that T*₂ does not change with rotation. In figure 2 measured values (black squares) for T₂ are plotted for several cases [12, 19, 30–32]. Assuming that for these systems T*₂ ≈ T₂, in figure 4, we plot the sensitivity for these specific cases (red squares).

Zoom In Zoom Out Reset image size

Up to this point we have been working within the adiabatic approximation. In the context of the NV centre, the adiabatic criterion is

$$\begin{eqnarray} &&\frac{1}{2}\left(\dot{\theta}^2+\dot{\phi}^2\sin^2\theta\right) \ll D^2. \end{eqnarray} \tag{ 13 }$$

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

$$\begin{eqnarray} &&D\gg\frac{k_{\mathrm{B}}T}{I}=\frac{15k_{\mathrm{B}}T}{8\pi \rho r^3}=15 k_{\mathrm{d}}, \end{eqnarray} \tag{ 14 }$$

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²/(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

$$\begin{eqnarray} &&\langle k \rangle=\frac{2k_{\mathrm{d}}}{1+t_{\mathrm{d}}^2D^2}. \end{eqnarray} \noindent \tag{ 15 }$$

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 ψ⁽⁰⁾_z', are plotted in figure 5(a). They show excellent agreement with our analytical expression:

$$\begin{equation} S_{\mathrm{m}}(t) = \frac{1}{3} + \frac{2}{3}\exp(-t/t_{\mathrm{m}}), \end{equation} \tag{ 16 }$$

where t_m = 1/(3〈k〉) is the population mixing time.

Zoom In Zoom Out Reset image size

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₁) 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 ψ⁽⁰⁾_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.

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_{\perp ^{\prime }}=\sqrt {B^2_{x^{\prime }}+B^2_{y^{\prime }}} \ll D$ and $B_{{\perp }^{\prime }}^2/D \ll \cal {R}$ [28, 29], under the secular approximation equation (1) becomes

$$\begin{eqnarray} &&H \approx \frac{1}{\hbar}D S_{z'}^2 + \gamma(B_z\cos\theta+ B_x\sin\theta\cos\phi)S_{z'}, \end{eqnarray} \tag{ 17 }$$

it is only the field component parallel to z' (the NV axis) which affects the phase.

6.1. 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 z- and 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

$$\begin{eqnarray} &&\fl \frac{1}{k_{\mathrm{d}}}\frac{\mathrm{d}}{\mathrm{d}t} p(\theta_2,\phi,\Phi,t|\theta_1) =- \frac{\partial}{\partial \theta_2}\left [\!\sin \theta_2 \frac{\partial}{\partial \theta_2}\left ( \frac{p(\theta_2,\phi,\Phi,t|\theta_1)}{\sin\theta_2}\right )\right ] \nonumber \\ &&+ \frac{1}{\sin^2\theta_2} \left (\frac{\partial}{\partial \phi}+\cos\theta_2 \frac{\partial}{\partial \Phi} \right )^2 p(\theta_2,\phi,\Phi,t|\theta_1) \nonumber \\ &&- \left (\frac{\gamma B_z}{k_{\mathrm{d}}} \cos\theta_2 + \frac{\gamma B_x}{k_{\mathrm{d}}} \sin\theta_2\cos\phi \right )\frac{\partial}{\partial \Phi} p(\theta_2,\phi,\Phi,t|\theta_1). \end{eqnarray} \tag{ 18 }$$

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

$$\begin{equation} \Delta B \sqrt T \approx \frac{1}{\gamma \alpha\sqrt{N \tau_{\mathrm{d}}}}, \end{equation} \tag{ 19 }$$

where 1/τ_d ≈ 1/T*₂ + 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.

6.2. 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²〉 [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²〉 = B²_fc ≪ 1/(γt_c)² which is isotropic, the Fokker–Planck equation is

$$\begin{eqnarray} &&\fl \frac{1}{k_{\mathrm{d}}} \frac{\mathrm{d}p(\theta_2,\Phi,t|\theta_1)}{\mathrm{d}t} =- \frac{\partial}{\partial \theta_2}\left [ \sin \theta_2 \frac{\partial}{\partial \theta_2}\left ( \frac{p\left(\theta_2 , \Phi; t|\theta_1\right)}{\sin\theta_2}\right )\right ]\nonumber\\ &&+ \left(\cot^2\theta_2+\frac{\gamma^2 B_{\mathrm{fc}}^2t_{\mathrm{c}}}{6k_{\mathrm{d}}}\right) \frac{\partial^2 p\left(\theta_2 , \Phi; t|\theta_1\right)}{\partial \Phi^2}, \end{eqnarray} \noindent \tag{ 20 }$$

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/(\gamma B_{\mathrm {fc}} \sqrt {t_{\mathrm {c}}})^2$. The implication of this result is that the decoherence rate is linear in B²_fct_c.

6.3. 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 ($\frac {\pi }{2} - \frac {\tau }{2} - \pi - \frac {\tau }{2} - \frac {\pi }{2}$) effectively rectifies an ac magnetic field of period 2τ, allowing its amplitude to be detected while masking fields oscillating at different frequencies. A spin-echo 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

$$\begin{equation} \tau_{n\pi}^3 \approx \frac{(n+1)^2}{4k_{\mathrm{B}}}\left (\frac{1}{\gamma B_{\mathrm{dc}}}\right )^2, \end{equation} \noindent \tag{ 21 }$$

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 ψ⁽⁰⁾_z' is given by

$$\begin{eqnarray} &&\fl P_0(\theta,\Phi_{\mathrm{field}}) = \frac{1}{2}+ \frac{1}{2}\cos^2 \Theta_1 \cos^2 \Theta_2 \pm \frac{1}{4} \sin^2 \Theta_1(1-\cos^2 \Theta_2) \cos(2\Phi_{\mathrm{ac}}), \end{eqnarray} \noindent \tag{ 22 }$$

where Θ₁ = a(π/2)sin θ and Θ₂ = 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 $\frac {1}{2} \sin \theta$:

$$\begin{eqnarray} &&S =\int_0^{\pi}\mathrm{d} \theta \frac{\sin\theta}{2} \int_0^{2\pi} \mathrm{d}\Phi_{\mathrm{field}}\frac{P_0(\theta,\Phi_{\mathrm{field}})}{2\pi}. \end{eqnarray} \tag{ 23 }$$

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₂ + 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.

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₂ and T*₂ 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].

This research was supported by the Australian Research Council Centre of Excellence for Quantum Computation and Communication Technology (CE110001027). LCLH was supported under an Australian Research Council Professorial Fellowship (DP0770715).

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. The populations of the spin sublevels after some time interval δt is

$$\begin{equation} P_1^{\prime}=P_1-k\delta t P_1+k \delta t P_0 +2\sqrt{P_1P_0}\sqrt{k \delta t}\cos\left(\Phi_0-\Phi_1+\delta \phi\right), \end{equation} \tag{ A.1 }$$

$$\begin{eqnarray} &&\fl P_0^{\prime}=P_0-2k \delta t P_0+k\delta t P_1+k \delta t P_{-1} -2\sqrt{P_1P_0}\sqrt{k \delta t }\cos\left(\Phi_0-\Phi_1+\delta \phi\right)\nonumber \\ &&-2\sqrt{P_{-1}P_0}\sqrt{k \delta t}\cos\left(\Phi_0-\Phi_{-1}+\delta \phi\right), \end{eqnarray} \tag{ A.2 }$$

$$\begin{equation} P_{-1}^{\prime}=P_{-1}-k\delta t P_{-1}+k \delta t P_0 +2\sqrt{P_{-1}P_0}\sqrt{k \delta t}\cos\left(\Phi_0-\Phi_{-1}-\delta \phi\right), \end{equation} \tag{ A.3 }$$

where

$$\begin{eqnarray} &&k=\frac{\delta \theta^2}{D^2\delta t^3}\left[1-\cos\,(\delta tD\right)]. \end{eqnarray} \noindent \tag{ A.4 }$$

Modelling the crystal's Brownian rotation as a random walk of small uncorrelated angular steps of duration δt and angular change δθ on timescales much longer than δt the ensemble average of ${\dot P}_{\alpha }$ can be determined

$$\begin{equation} \frac{\mathrm{d}\langle P_1\rangle}{\mathrm{d}t}=-\langle k \rangle \langle P_1 \rangle+\langle k \rangle \langle P_{-1} \rangle, \end{equation} \tag{ A.5 }$$

$$\begin{equation} \frac{\mathrm{d}\langle P_0\rangle}{\mathrm{d}t}=\langle k \rangle \langle P_1 \rangle+\langle k \rangle \langle P_{-1} \rangle-2\langle k\rangle\langle P_0 \rangle, \end{equation} \tag{ A.6 }$$

$$\begin{equation} \frac{\mathrm{d}\langle P_{-1}\rangle}{\mathrm{d}t}=-\langle k \rangle \langle P_{-1} \rangle+\langle k \rangle \langle P_{1} \rangle. \end{equation} \noindent \tag{ A.7 }$$

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

$$\begin{eqnarray} &&\langle{\boldsymbol{\omega}}(t)\cdot {\boldsymbol{\omega}}(0) \rangle = \langle{\boldsymbol{\omega}}(0)\cdot {\boldsymbol{\omega}}(0) \rangle\mathrm{ e}^{-t/t_{\mathrm{d}}}, \end{eqnarray} \noindent \tag{ 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

$$\begin{eqnarray} &&p(\delta t)=\frac{\delta t}{t^2_{\mathrm{d}}}\mathrm{e}^{-\delta t/t_{\mathrm{d}}}. \end{eqnarray} \tag{ A.9 }$$

We do not need to know the exact distribution of angular velocities, merely that 〈(δθ)²/(δt)²〉 = 2/3〈ω·ω〉 = 2k_d/t_d. The ensemble average population mixing rate is then given by

$$\begin{equation} \langle k \rangle=\frac{1}{D^2}\left\langle\frac{\delta \theta^2 }{\delta t^2}\right\rangle\left\langle \frac{1}{\delta t}\left(1- \cos\,(\delta t D)\right)\right\rangle \end{equation} \tag{ A.10 }$$

$$\begin{equation} \hphantom{\langle k \rangle}=\frac{2k_{\mathrm{d}}}{1+t_{\mathrm{d}}^2D^2}. \end{equation} \tag{ A.11 }$$