In this example we simulate ambiguous resonances in multipulse quantum sensing with nitrogen vacancy (NV) centers in diamonds. This notebook closely follows the work from “Ambiguous Resonances in Multipulse Quantum Sensing with Nitrogen Vacancy Centers” available at:
Tsunaki, A. Singh, K. Volkova, S. Trofimov, T. Pregnolato, T. Schröder, & B. Naydenov. (2025). Ambiguous resonances in multipulse quantum sensing with nitrogen-vacancy centers. Physical Review A, 111(2), 022606. doi: 10.1103/PhysRevA.111.022606.
As this is a quite complex problem, physics discussions will be kept to a minimum here in favor of a focus to the simulation. For a detailed explanation of the problem, please refer to the original paper.
[1]:
import numpy as np
from qutip import tensor, jmat, qeye
from quaccatoo import NV, XY8, Analysis
1. Field Misalignment and \(^{15}\text{N}\) Coupling
We begin simulating a XY8-2 sequence with a misaligned magnetic field, as in Sec. IV A of the paper. Initially, we define the NV system with the experimental parameters of the paper and the XY8 pulse sequence.
[2]:
qsys_1 = NV(
N=15,
B0 = 39.4,
units_B0='mT',
theta=2.6,
units_angles='deg'
)
w1 = 20
XY8_2_15N = XY8(
free_duration = np.linspace(.25, .36, 100),
system = qsys_1,
M=2,
pi_pulse_duration = 1/2/w1,
h1 = w1*qsys_1.MW_h1,
# MW frequency for the ms=0 --> ms=+1 state transition
pulse_params = {'f_pulse': qsys_1.MW_freqs[1]},
# Number of time steps in each MW pulse
time_steps=100,
)
XY8_2_15N.plot_pulses(tau=.1, figsize=(12, 4))
This leads to a \(\pi\)-pulse duration shorter than \(1/(2\omega_1)\), resulting from a hyperfine enhancement effect.
Given the defined system with the experimental parameters, the simplest dynamical decoupling sequence for suppressing the noise from the environment is the Hahn echo technique. The sequence is composed by two free evolutions of duration \(\tau\), with an intermediate \(\pi\)-pulse to the electron spin. In addition, \(\pi/2\)-pulses are included before and after the sequence to drive the electron spin from and to the quantization axis. In NMR, the last projective pulse is not required though, as the measurement is performed perpendicular to the quantization axis. The intermediate \(\pi\)-pulse causes a refocusing of static or slow dephasings, such as those arising from the 14N nuclear spin. This allows the coupled 13C nuclear spin to be probed by the NV.
For performing a Hahn echo simulation in the software component, we use the Hahn class:
Here we can observe the 16 \(\pi\)-pulses plus the initial and final \(\pi/2\) pulses. Now to run the experiment as in Fig. 3 (d):
[3]:
XY8_2_15N.run()
Analysis(XY8_2_15N).plot_results(title=r'XY8-2 $B_0=39$ mT $\theta_0=2.6^\circ$')
Increasing the XY8 order. First for \(M=4\):
[4]:
XY8_4_15N = XY8(
free_duration = np.linspace(.25, .36, 100),
system = qsys_1,
M=4,
pi_pulse_duration = 1/2/w1,
h1 = w1*qsys_1.MW_h1,
pulse_params = {'f_pulse': qsys_1.MW_freqs[1]},
)
XY8_4_15N.run()
Analysis(XY8_4_15N).plot_results(title=r'XY8-4 $B_0=39$ mT $\theta_0=2.6^\circ$')
\(M=6\):
[5]:
XY8_6_15N = XY8(
free_duration = np.linspace(.25, .36, 100),
system = qsys_1,
M=6,
pi_pulse_duration = 1/2/w1,
h1 = w1*qsys_1.MW_h1,
pulse_params = {'f_pulse': qsys_1.MW_freqs[1]},
)
XY8_6_15N.run()
Analysis(XY8_6_15N).plot_results(title=r'XY8-6 $B_0=39$ mT $\theta_0=2.6^\circ$')
\(M=8\):
[6]:
XY8_8_15N = XY8(
free_duration = np.linspace(.25, .36, 100),
system = qsys_1,
M=8,
pi_pulse_duration = 1/2/w1,
h1 = w1*qsys_1.MW_h1,
pulse_params = {'f_pulse': qsys_1.MW_freqs[1]},
)
XY8_8_15N.run()
Analysis(XY8_8_15N).plot_results(title=r'XY8-8 $B_0=39$ mT $\theta_0=2.6^\circ$')
Now, changing the field for XY8-2 as in Fig. 3 (c). First, \(B_0=33\) mT and \(\theta_0 = 5.1^\circ\):
[7]:
qsys_2 = NV(
N=15,
B0 = 33,
units_B0='mT',
theta=5.1,
units_angles='deg'
)
XY8_15N_B1 = XY8(
free_duration = np.linspace(.22, .45, 100),
system = qsys_2,
M=2,
pi_pulse_duration = 1/2/w1,
h1 = w1*qsys_2.MW_h1,
pulse_params = {'f_pulse': qsys_2.MW_freqs[1]}
)
XY8_15N_B1.run()
Analysis(XY8_15N_B1).plot_results(title=r'XY8-2 $B_0=33$ mT $\theta_0=5.1^\circ$')
\(B_0=26\) mT and \(\theta_0 = 9.3^\circ\):
[8]:
qsys_3 = NV(
N=15,
B0 = 26,
units_B0='mT',
theta = 9.3,
units_angles='deg'
)
XY8_15N_B2 = XY8(
free_duration = np.linspace(.22, .45, 100),
system = qsys_3,
M=2,
pi_pulse_duration = 1/2/w1,
h1 = w1*qsys_3.MW_h1,
pulse_params = {'f_pulse': qsys_3.MW_freqs[1]}
)
XY8_15N_B2.run()
Analysis(XY8_15N_B2).plot_results(title=r'XY8-2 $B_0=26$ mT $\theta_0=9.3^\circ$')
\(B_0=17\) mT and \(\theta_0 = 17^\circ\):
[9]:
qsys_4 = NV(
N=15,
B0 = 17,
units_B0='mT',
theta=17,
units_angles='deg'
)
XY8_15N_B3 = XY8(
free_duration = np.linspace(.25, .45, 100),
system = qsys_4,
M=2,
pi_pulse_duration = 1/2/w1,
h1 = w1*qsys_4.MW_h1,
pulse_params = {'f_pulse': qsys_4.MW_freqs[1]}
)
XY8_15N_B3.run()
Analysis(XY8_15N_B3).plot_results(title=r'XY8-2 $B_0=17$ mT $\theta_0=17^\circ$')
If we now include a pulse length error, as in Fig. 6, we have:
[10]:
XY8_15N_B3_errors = XY8(
free_duration = np.linspace(.25, .45, 100),
system = qsys_4,
M=2,
pi_pulse_duration = 1/2/w1*.8,
h1 = w1*qsys_4.MW_h1,
pulse_params = {'f_pulse': qsys_4.MW_freqs[1]}
)
XY8_15N_B3_errors.run()
Analysis(XY8_15N_B3_errors).plot_results(title=r'XY8-2 $B_0=17$ mT $\theta_0=17^\circ$, pulse lenght=0.8 $t_{\pi}$')
2. 13C Coupling
As in Sec. IV B, we need to define the coupling Hamiltonian \(\hat{H}_2\) in the same way as Eq. 3 of the paper. We begin defining the coupling tensor shown in Appendix D and the Zeeman Hamiltonian.
[11]:
A = np.empty((3, 3), dtype=object)
A[0, 0] = -0.25
A[0, 1] = -1.85
A[0, 2] = -.49
A[1, 0] = A[0, 1]
A[1, 1] = 0
A[1, 2] = .01
A[2, 0] = A[0, 2]
A[2, 1] = A[1, 2]
A[2, 2] = 1.01
HhfC = (
A[0,0]*tensor(jmat(1,'x'), qeye(3), jmat(1/2,'x')) + A[0,1]*tensor(jmat(1,'x'), qeye(3), jmat(1/2,'y')) + A[0,2]*tensor(jmat(1,'x'), qeye(3), jmat(1/2,'z'))
+ A[1,0]*tensor(jmat(1,'y'), qeye(3), jmat(1/2,'x')) + A[1,1]*tensor(jmat(1,'y'), qeye(3),jmat(1/2,'y')) + A[1,2]*tensor(jmat(1,'y'), qeye(3), jmat(1/2,'z'))
+ A[2,0]*tensor(jmat(1,'z'), qeye(3), jmat(1/2,'x')) + A[2,1]*tensor(jmat(1,'z'), qeye(3), jmat(1/2,'y')) + A[2,2]*tensor(jmat(1,'z'), qeye(3), jmat(1/2,'z'))
)
def Hzc(B0, theta):
return 10.705e-3*B0*tensor(qeye(3), qeye(3), np.cos(theta*np.pi/180)*jmat(1/2,'z') + np.sin(theta*np.pi/180)*jmat(1/2,'x'))
Now we define the NV Hamiltonian and use the add_spin method to introduce the 13C nuclear spin.
[12]:
qsys_5 = NV(
N=14,
B0 = 18,
units_B0='mT',
theta=2.34,
units_angles='deg',
)
qsys_5.add_spin(HhfC + Hzc(qsys_5.B0, qsys_5.theta))
XY8_13C_B1 = XY8(
free_duration = np.linspace(.4, .85, 100),
system = qsys_5,
M=2,
pi_pulse_duration = 1/2/w1,
h1 = w1*qsys_5.MW_h1,
pulse_params = {'f_pulse': qsys_5.MW_freqs[1]}
)
XY8_13C_B1.run()
Analysis(XY8_13C_B1).plot_results(title=r'XY8-2 $^{13}$C, $B_0=18$ mT, $\theta_0=2.34^\circ$')
Increasing the field as in Fig. 4 (a) for \(B_0=31\) mT and \(\theta_0 = 1.35^\circ\):
[13]:
qsys_6 = NV(
N=14,
B0 = 31,
units_B0='mT',
theta=1.35,
units_angles='deg',
)
qsys_6.add_spin(HhfC + Hzc(qsys_6.B0, qsys_6.theta))
XY8_13C_B2 = XY8(
free_duration = np.linspace(.4, .85, 100),
system = qsys_6,
M=2,
pi_pulse_duration = 1/2/w1,
h1 = w1*qsys_6.MW_h1,
pulse_params = {'f_pulse': qsys_6.MW_freqs[1]}
)
XY8_13C_B2.run()
Analysis(XY8_13C_B2).plot_results(title=r'XY8-2 $^{13}$C, $B_0=31$ mT, $\theta_0=1.35^\circ$')
\(B_0=40\) mT and \(\theta_0 = 1.14^\circ\):
[14]:
qsys_7 = NV(
N=14,
B0 = 40,
units_B0='mT',
theta=1.14,
units_angles='deg',
)
qsys_7.add_spin(HhfC + Hzc(qsys_7.B0, qsys_7.theta))
XY8_13C_B3 = XY8(
free_duration = np.linspace(.4, .85, 100),
system = qsys_7,
M=2,
pi_pulse_duration = 1/2/w1,
h1 = w1*qsys_7.MW_h1,
pulse_params = {'f_pulse': qsys_7.MW_freqs[1]}
)
XY8_13C_B3.run()
Analysis(XY8_13C_B3).plot_results(title=r'XY8-2 $^{13}$C, $B_0=40$ mT, $\theta_0=1.14^\circ$')
Now, we keep the field constant and increase the order, as in Fig. 4(b). \(M=4\):
[15]:
XY8_4_13C = XY8(
free_duration = np.linspace(.4, .85, 100),
system = qsys_7,
M=4,
pi_pulse_duration = 1/2/w1,
h1 = w1*qsys_7.MW_h1,
pulse_params = {'f_pulse': qsys_7.MW_freqs[1]}
)
XY8_4_13C.run()
Analysis(XY8_4_13C).plot_results(title=r'XY8-4 $^{13}$C, $B_0=40$ mT, $\theta_0=1.14^\circ$')
\(M=6\):
[16]:
XY8_6_13C = XY8(
free_duration = np.linspace(.4, .85, 100),
system = qsys_7,
M=6,
pi_pulse_duration = 1/2/w1,
h1 = w1*qsys_7.MW_h1,
pulse_params = {'f_pulse': qsys_7.MW_freqs[1]}
)
XY8_6_13C.run()
Analysis(XY8_6_13C).plot_results(title=r'XY8-6 $^{13}$C, $B_0=40$ mT, $\theta_0=1.14^\circ$')
\(M=8\):
[17]:
XY8_8_13C = XY8(
free_duration = np.linspace(.4, .85, 100),
system = qsys_7,
M=8,
pi_pulse_duration = 1/2/w1,
h1 = w1*qsys_7.MW_h1,
pulse_params = {'f_pulse': qsys_7.MW_freqs[1]}
)
XY8_8_13C.run()
Analysis(XY8_8_13C).plot_results(title=r'XY8-8 $^{13}$C, $B_0=40$ mT, $\theta_0=1.14^\circ$')
3. RXY8 with H2
As discussed in Sec. IV C, the finite length of the pulses can lead to ambiguous resonances. Here, we define the NV Hamiltonian as usual and a time-dependent magnetic field \(B_2(t)\) as in Eq. 4.
[2]:
qsys_8 = NV(
N=15,
B0 = 40,
units_B0='mT'
)
w1 = 20
f2 = 5.5
B2z = 0.3
def B2(t, B2z=B2z, f2=f2):
return B2z*np.sin(f2*t)
H2 = [tensor(jmat(1,'z'), qeye(2)), B2]
XY8_12_H2 = XY8(
free_duration = np.linspace(.06, 0.17, 100),
system = qsys_8,
M=12,
pi_pulse_duration = 1/2/w1,
h1 = w1*qsys_8.MW_h1,
pulse_params = {'f_pulse': qsys_8.MW_freqs[0]},
H2 = H2
)
# the number of CPUs for the parallel computation can be specified with the num_cpus keyword argument
XY8_12_H2.run()
Analysis(XY8_12_H2).plot_results()
To show that these resonances are indeed from the finite pulse duration, we now simulate it with delta pulses.
[3]:
XY8_12_H2_delta = XY8(
free_duration = np.linspace(.06, 0.17, 100),
system = qsys_8,
M=12,
pi_pulse_duration = 0,
Rx = qsys_8.MW_Rx[0],
Ry = qsys_8.MW_Ry[0],
H2 = H2
)
XY8_12_H2_delta.run()
Analysis(XY8_12_H2_delta).plot_results()
These ambiguous resonances can be supressed with phase randomization, for that we use the XY8 class with the RXY8 parameter set to True.
[4]:
RXY8_12_H2 = XY8(
free_duration = np.linspace(.06, 0.17, 100),
system = qsys_8,
M=12,
pi_pulse_duration = 1/2/w1,
h1 = w1*qsys_8.MW_h1,
pulse_params = {'f_pulse': qsys_8.MW_freqs[0]},
H2 = H2,
RXY8=True
)
RXY8_12_H2.run()
Analysis(RXY8_12_H2).plot_results()
