Nitrogen-Vacancy Centers

Color centers in diamond and silicon carbide have been studied for years and demonstrated significant potential in nanoscale sensing of magnetic and electric fields and as long-lived quantum memories. The nitrogen vacancy center is the most well studied point defect in diamond and provides nanoscale resolution, long coherence times at room temperature, all-optical spin state initialization and readout as well as coherent manipulation of spin states with resonant microwave fields. SimOS offers a submodule NV that contains a series of high-level helper functions to simplify simulations of photophysics and nanoscale NMR experiments with NV centers.

System Initialisation

The NV submodule provides an NV-specific class NVSystem to initialise the multipartite quantum system of NV centers. A series of keyword arugments allow to account for various scenarios.

Keyword arguments for NV system initialisation.

Keyword Argument	Default	Description
optics	True	If True, electronic structure of the NV is included.
orbital	False	If True, the excited electronic state is initialised as an orbital doublet.
nitrogen	True	If True, coupling to a nuclear nitrogen spin is included.
synthetic	True	If True, we assume a synthetic (i.e. non-natural) NV center, formed by ion implanation, and the coupled nitrogen spin has a spin \(I= 1/2\). If False, a spin \(I = 1\) is assumed for the coupled nitrogen.
further_spins	[]	Describes further nuclear spin that are coupled to the NV center.

If all keyword-arguments are omitted, the initialised system is representative of an NV center at room temperature and native SimOS code for its construction would be:

# NV Center Electron Spin
S = {'val': 1, 'name':'S', 'type': 'NV-'}
# Nuclear Nitrogen Spin
I = {'val': 1/2, 'name':'I', 'type': 'N'}
# Electronic States
GS = {'val': 0 , 'name':'GS'}
ES = {'val': 0 , 'name':'ES'}
SS = {'val': 0 , 'name':'SS'}
s = sos.System([([(GS, ES), S], SS), N])

The following scheme visualizes the construction of the composite Hilbert space in a stepwise fashion.

NV Hamiltonian

The NVSystem provides users with a method field_hamiltonian() which returns the Hamiltonian of the NV electron spin (\(S=1\)) in the electronic ground and excited state. The ground state spin Hamiltonian

\[\begin{split}\hat{H}_{\mathrm{GS}} = D_{\mathrm{GS}} \left[ \hat{S}_z^2 - \frac{2}{3} \mathbb{I}_3 \right] + \gamma_{\mathrm{NV}}\hat{\vec{S}} \cdot \vec{B} \\ + d_{\mathrm{GS}}^{\parallel} (E_z + \delta_z) + d_{\mathrm{GS}}^{\perp}(E_x + \delta_x) (\hat{S}_y^2 - \hat{S}_x^2) + d_{\mathrm{GS}}^{\perp}(E_y + \delta_y)(\hat{S}_x\hat{S}_y + \hat{S}_y\hat{S}_x)\end{split}\]

contains the ground state zero-field splitting of magnitude \(D_{\mathrm{GS}}\) , the Zeeman interaction of the NV electron spin (gyromagnetic ratio \(\gamma_{\mathrm{NV}}\) ) in an external magnetic field \(\vec{B}\) and the stark interaction with the sum of electric \(\vec{E}\) and strain \(\vec{\delta}\) fields. We distinguish parallel \(d_{\mathrm{GS}}^{\parallel}\) and perpendicular \(d_{\mathrm{GS}}^{\perp}\) electric dipole moments of the NV center ground state.

For the excited state Hamiltonian, we distinguish room- and low-temperature situations. At low temperature the orbital substructure of the excited state is explicitly taken into account and the excited state Hamilonian is

\[\begin{split}\hat{H}_{\mathrm{ES}}^{\mathrm{LT}} = D_{\mathrm{ES}}^{\parallel}\mathbb{I}_2 \otimes \left[ \hat{S}_z^2 - \frac{2}{3} \mathbb{I}_3 \right] + D_{\mathrm{ES}}^{\perp}\left[\hat{\sigma}_z \otimes (\hat{S}_y^2 - \hat{S}_x^2) - \hat{\sigma}_x \otimes (\hat{S}_y\hat{S}_x + \hat{S}_x\hat{S}_y) \right] \\ - \lambda^{\parallel}\hat{\sigma}_y\otimes\hat{S}_z + \lambda^{\perp} \left[\hat{\sigma}_z\otimes(\hat{S}_x\hat{S}_z + \hat{S}_z\hat{S}_x ) - \hat{\sigma}_x\otimes(\hat{S}_y\hat{S}_z + \hat{S}_z\hat{S}_y ) \right] \\ + \gamma_{\mathrm{NV}}\mathbb{I}_2 \otimes (\hat{\vec{S}} \cdot \hat{\vec{B}}) + \gamma_{\lambda} (\hat{\vec{\sigma}}\cdot \vec{B}) \otimes \mathbb{I}_3 \\ + d_{\mathrm{ES}}^{\parallel} (E_z + \delta_z) \mathbb{I}_6+ d_{\mathrm{ES}}^{\perp} (E_x + \delta_x)\hat{\sigma}_z\otimes \mathbb{I}_3 + d_{\mathrm{ES}}^{\perp} (E_y + \delta_y)\hat{\sigma}_x\otimes \mathbb{I}_3\end{split}\]

where \(\hat{\sigma}_i\) (\(i = x, y, z\)) are standard Paul matrices for the excited state orbital operators, \(D_{\mathrm{ES}}^{\parallel}\) and \(D_{\mathrm{ES}}^{\perp}\) are the parallel and perpendicular of the excited state zero field splitting, \(\gamma_{\lambda}\) is the orbital gyromagnetic ratio and \(d_{\mathrm{ES}}^{\parallel}\) and \(d_{\mathrm{ES}}^{\perp}\) are the parallel and perpendicular components of the excited state electric dipole moment. At room temperature, the orbital doublet is not explicitly included and the excited state Hamiltonian simplifies accordingly. The field_hamiltonian() function probes whether the NVSystem was initialised with the orbital substructure to distinguish low- and room-temperature cases.

Since all NV-specific constants are automatically sourced from the SimOS library, the user only has to specify magnetic, electric and strain fields to field_hamiltonian(). The magnetic field is specified with a keyword argument and the method defaults to zero magnetic field if the argument is omitted. Since the excited state dipole moment (i.e. \(d_{\mathrm{ES}}^{\parallel}\) and \(d_{\mathrm{ES}}^{\perp}\) ) are not known, electric and strain fields are not input as fields but rather as frequencies which represent products of electric dipole moment and combined electric and strained fields.

Photophysics

The spin state of the negatively charged NV center can be optically prepared and read out using green light illumination enabling single-spin detection spanning room temperature to cryogenic temperature. To simulate the spin-dependent photoluminescence dynamics of single NV centers, the electronic level structure must be taken into account.

In the simplest case, valid at room temperature, the negatively charged NV center electronic structure may be described by three electronic levels: a triplet ground state, a triplet excited state and a metastable singlet state. Here, we consider off-resonant excitation of the NV center from the ground to the excited state, which is commonly achieved with green (\(\sim\) 532 nm) laser light. The excited state has an optical lifetime of about 16 ns and decays either radiatively by emission of a red photon, or non-radiatively via an intersystem crossing (ISC) to the metastable singlet state. This ISC is strongly spin dependent, resulting in (i) higher fluorescence intensity for the \(m_S = 0\) spin state versus the \(m_S = 1\) spin state and (ii) polarization of the spin state to \(m_S = 0\) throughout multiple excitation-decay cycles.

The simplest simulation of the NV center photophyiscs is purely incoherent, i.e. it does not take into account any coherent spin dynamics and solely simuates the effect of incoherent optical excitation and decay events with the simplified room-temperature model. To perform such a simulation, we construct the representative quantum system as:

NV = sos.NV.NVSystem(nitrogen = False)

The NVSystem provides a method transition_operators() which returns all optical decay rates of the NV center. At room temperature, these can be assumed to be independent of exact temperature and external magnetic and electric fields. A parameter beta allows to indicate the power of the laser illumination relative to the NV center saturation power (saturation for beta = 1).

# c_ops_on transition operators if laser is on, c_ops_off if laser is off
c_ops_on, c_ops_off = NV.transition_operators(beta = 0.2)

We simulate the photoluminescence dynamics during a readout laser pulse for an NV center in the \(m_S = 0\) and \(m_S = 1\) state, respectively, with the code provided below.

pl0 = []
pl1 = []
U = sos.evol( None, dt , c_ops = c_ops_on )
rho0 = (NV.Sp[0]).copy()
rho1 = (NV.GSid*NV.Sp[1]).copy()
for t in tax :
    pl0.append(sos.expect(NV.ESid, rho_0 ) )
    pl1.append(sos.expect(NV.ESid, rho_1 ) )
rho0 = sos . applySuperoperator (U , rho0 )
rho1 = sos . applySuperoperator (U , rho1 )

While this simple and fully incoherent model sufficiently captures NV center photophysics at elevated temperatures, it does not accurately describe effects at cryogenic and intermediate temperature regimes. Here, interplay of spin and orbital dynamics in the excited state may result in fast spin relaxation and vanishing ODMR contrast. Accurate simulations must be performed with a quantum master equation which includes coherent spin dynamics under external magnetic and electric (strain) fields as well as incoherent excitation and decay dynamics accounting for the orbital character of the excited state and phonon-induced hopping. We initialize the representative quantum system and define the system Hamiltonian and collapse operators for a specific temperature T and magnetic and electric (strain) fields Bvec, Evec using the transition_operators() and field_hamiltonian() methods of the NV submodule.

Again, we simulate the photoluminescence dynamics during a readout laser pulse for an NV center in the \(m_S = 0\) and \(m_S = 1\) state, respectively.

pl0 = []
pl1 = []
U  = sos.evol((HGS+HES), dt, c_ops = c_ops_on)
rho0 = (NV.GSid*NV.Sp[0]).copy()
rho1 = (NV.GSid*NV.Sp[1]).copy()
for t in tax:
    pl0.append(sos.expect(NV.ESid, rho_0))
    pl1.append(sos.expect(NV.ESid, rho_1))
    rho0 = sos.applySuperoperator(U,rho0)
    rho1 = sos.applySuperoperator(U,rho1)

The figure below illustrates the results for three different temperature, clearly illustrating vanishing ODMR contrast at intermediate temperature regimes where phonon-mediated hopping interferes with optical spin state preparation and readout.

Note

For a full code example on photophysics simulations you can visit our Virtual Lab and follow the NV.ipynb notebook.

Nanoscale NMR

NV centers have been used to perform nanoscale NMR measurements of proximal nuclei both within the diamond lattice and, when within several nm of the surface, external nuclei in chemical species. The NV center is further capable of detecting the free induction decay of nearby nuclear spins using a measurement technique referred to as weak measurements. Here, the nuclear spin is polarized and rotated with a \(\frac{\pi}{2}\) pulse to initiate precession. Repeated weak measurements on the NV center, consisting of optical spin initialization, XY-8 dynamical decoupling frequency encoding, and optical readout are used to track the nuclear spin precession.

In the simplest case, assuming ideal optical initialization and readout, this measurement protocol can be simulated without considering the photophysics of the NV center. The quantum system is then constructed from the NV center’s electronic spin \(S=1\) which is coupled to a single \(^{13}\) C nuclear spin \(S=1/2\).

C = {"name": "C", "val": 1/2}
NV = sos.NV.NVSystem(further_spins = [C], optics = False, nitrogen = False)

In a next step, we prepare the initial state that represents the polarized and rotated :code:`^{13}`C spin coupled to a polarized NV center and define the system Hamiltonian from the native system operators. We work in the rotating frame for the NV center where the hyperfine interaction is truncated to the secular and pseudo-secular contributions nd utilize the :code:`^{13}`C gyromagnetic ratio of SimOS.

# Initial State
rho0 = sos.state(s, "S[0],C[0.5]")
rho1 = sos.rot(s.Cx,np.pi/2,rho0)
# System Hamiltonian
H0 = sos.yC13*B0*s.Cz + apara*s.Sz*s.Cz + aperp*s.Sx*s.Cx

Finally, we perform the sensing protocol, consisting of a series of \(N\) XY-8 blocks that are interspaced by free evolution periods. Again, we make use of the NV sub-module functionality and utilize the XY8() and meas_NV routines to apply the dipolar decoupling sequence and perform NV center readout. We obtain the free-induction decay of the nuclear spin and, after Fourier transformation, the corresponding spectrum.

N = 500
frf = sos.w2f(B0*sos.yC13)
twait = 0.125/frf
tau = 1/frf/2
store = []
rho = rho1.copy()
for i in range(N):
    rho = sos.rot(s.Sop_y_red,np.pi/2,rho)
    rho = sos.NV.XY8(H0,tau,s,rho,N=16)
    rho = sos.rot(s.Sop_x_red,np.pi/2,rho)
    meas,rho = sos.NV.meas_NV(rho,s)
    store.append(meas)
    rho = sos.evol(H0,twait,rho)
Tsample = 16*tau + twait
foffset = np.round(Tsample*frf)/Tsample

Syntax Reference

class NVSystem(optics=True, orbital=False, nitrogen=True, natural=False, further_spins=[], method='qutip')

A class for initialising quantum systems of NV centers.

field_hamiltonian(Bvec=array([0, 0, 0]), EGS_vec=array([0, 0, 0]), EES_vec=array([0, 0, 0]))

Hamiltonian of the NV center electron spin in the presence of magnetic and electric/strain fields.

Keyword Arguments:

• B_vec – Cartesian magnetic field vector, in Tesla.

• EGS_vec – Cartesian vector representing products of electrin/strain fields and NV electric dipole moment in the GS, in 1/s.

• EES_vec – Cartesian vector representing products of electrin/strain fields and NV electric dipole moment in the ES, in 1/s.

Returns:

Field Hamiltonian of the NV center.

transition_operators(T=298, beta=0.2, Bvec=array([0, 0, 0]), Evec=array([0, 0, 0]))

Collapse operators for laser- and phonon-induced incoherent transitions between electronic levels of the NV center. For a classical rate model of the NV center, only beta must be specified and all other parameters can be omitted.

Keyword Arguments:

• T – Temperature (in Kelvin), defaults to 298 K.

• beta – A parameter from 0-1 to indicate laser power relative to saturation power (saturation for beta = 1).

• Bvec – Cartesian magnetic field vector, in Tesla.

• Evec – Cartesian vector representing products of electrin/strain fields and NV electric dipole moment in the ES, in Hz.

Returns:

Two lists of collapse operators, first list does include laser excitation, second list does not include laser excitation.

Wcp(f, alpha, n, tau)

Filter function of the ideal CPMG sequence.

Parameters:

• f – frequency of signal

• alpha – Phase of signal

• n – number of pulses

:param tau : inter-pulse spacing

Returns:

The filter function.

See https://doi.org/10.1103/RevModPhys.89.035002 page 18

XY8(H, tau, spinsystem, *rho, N=8, phi=0, NV_name='S', c_ops=[], dthet=0)

Simulation of XY8 sequence

.x..y..x..y..y..x..y..x. (. denotes tau/2, x and y are the pi pulse axes)

Parameters:

• H – Hamiltonian for free evolution.

• tau – interpulse delay

• spinsystem – Spinsystem to take the rotation operators out of.

• *rho – Density matrix or state to which the sequence is applied. If none is given, the Propagator will be returned (only for direct method, see below)

• N – Number of pulses (must be multiple of 8)

• phi – change rotation axis, by default 0.

• NV_name – lookup prefix for the rotation operators in spinsystem (will look up ‘[NV_name]x_red’ and corresp. y)

• method – ‘direct’: use evol function and assume no relaxation or timedependence of the Hamiltonian

• method – ‘mesolve’: Pass the collaps operators c_ops to the Master equation solver of Qutip and use this solver for the timeevolution between the pi pulses (slow!).

• args – Arguments for the mesolve routine

• options – Options for the mesolve routine

• dthet – Pulse error in z rotation

References

https://doi.org/10.1016/0022-2364(90)90331-3

auto_pairwise_coupling(spin_system, approx=False, only_to_NV=False, output=True)

Returns the combined Hamiltonian for all pairwise dipolar couplings of a spin system using the tabulated isotropic gyromagnetic ratios of their type and positions.

This routine requires that the dictionaries of the spin members specify their spin type and position.

Parameters:

spin_system (System) – An instance of the System class.

Keyword Arguments:

• approx (bool) – If True, a secular approximation is performed. If False, the full coupling is considered.

• only_to_NV (bool) – If True, only couplings to NV center electronic spins are considered.

• output (bool) – If True, output is printed.

coupl2geom(apara, aperp, angular=False, ynuc=67278750.0, yel=-176085963023.0)

Converts the parallel and perpendicular coupling frequencies as obtained from correlation spectroscopy to r and theta assuming a point dipole model.

Returns [r, theta]

Parameters: apara : parallel coupling frequency aperp : perpendicular coupling frequency

Optional parameters: angular : by default False, assuming frequencies in Hz as input, set to True for inputting angular frequencies ynuc : gyromagnetic ratio for the nuclear spin, by default built-in value for C13 yel : gyromagnetic ratio for the electron spin, by default built-in value for electron

decay_rates(T=298)

Optical decay rates for the NV center at a given temperature. See DOI: 10.1103/PhysRevB.108.085203 for in-depth explanation.

Parameters:

T – Temperature in Kelvin.

Returns:

A dictionary with optical decay rates for fluorescent emission, spin-selective inter-system crossing and decay from the shelving state to the ground state.

exp2cts(x, contrast, ref_cts, ref=0.5)

Calculates photoluminescence counts for a given spin state initialization fidelity.

Parameters:

• x – Expectation value for ms(0) sublevel.

• contrast – ODMR contrast of the NV.

• ref_cts – Photoluminescence counts of the NV center (cts/s) in the ms(0) state.

Returns:

Photoluminescence counts (cts/s)

get_NVaxes(orientation=(0, 0, 1), axisind=(1, 1, 1))

Calculates the NV axis in the laboratory frame (z axis orthogonal to diamond surface) for the specified diamond surface termination and NV axis. Input: - Orietntation: Diamond surface ermination, most commonly (0,0,1) - Axisind: NV axis specification, e.g. (1,1,1) Outout: - np.array specifying the NV axis

lab2NVcoordinates(*args, orientation=(0, 0, 1), axisind=(1, 1, 1), input_coords='cart', output_coords='cart')

Transforms vector from NV coordinate system to laboratory system, arbitrary diamond surface termination and any NV axis

meas_NV(rho, spinsystem, NV_name='S')

Perform measurement operation on NV.

Parameters:

• rho – Density matrix of the state to measure.

• spinsystem – Spinsystem class object containing all operators

• NV_name – NV center spin name. If None, search for all NVs

Returns:

expectation value, rho after measurement.

normalize_data(data, upper_ref, lower_ref)

Normalizes data such that the upper reference correponds to 1, the lower reference trace to 0.

Parameters:

• data – Data trace to normalize

• upper_ref – mean value of the upper reference

• lower_ref – mean value of the lower reference

Returns:

Normalized data trace.

phonon_rates(T, E_perp)

Rates for phonon-mediated hopping between the orbital branches of the NV center excited state. See DOI: 10.1103/PhysRevB.108.085203 for in-depth explanation.

Parameters:

• T – Temperature in Kelvin.

• E_perp – excited state in-plane strain/el. field in Hz.

Returns:

Up- and down rates for phonon-mediated hopping in Hz.