ODMR Simulation

Optically Detected Magnetic Resonance (ODMR) is the standard technique for reading out NV center spin states. This tutorial shows how to simulate ODMR spectra.

CW-ODMR

In continuous-wave ODMR, we measure fluorescence while sweeping the microwave frequency.

import numpy as np
import matplotlib.pyplot as plt
from sim import HamiltonianBuilder, LindbladSolver
from sim.hamiltonian.terms import ZFS, Zeeman, MicrowaveDrive
from sim.states import ground_state, projector_ms

# Parameters
D = 2.87           # GHz
B = 10             # mT (gives ~280 MHz splitting)
mw_power = 1       # MHz Rabi frequency

# Frequency sweep
freqs = np.linspace(2.6, 3.15, 200)  # GHz
fluorescence = []

for f in freqs:
    # Detuning from ms=0 ↔ ms=-1 transition
    # The two transitions are at D ± γB/2π ≈ 2.87 ± 0.28 GHz

    H = HamiltonianBuilder()
    H.add(ZFS(D=D))
    H.add(Zeeman(B=B))
    H.add(MicrowaveDrive(omega=mw_power, detuning=f - D))

    solver = LindbladSolver(H)
    solver.add_t2_dephasing(gamma=5e6)  # T2* = 200 ns

    # Find steady state
    rho_ss = solver.steady_state(t_max=1e-5)

    # Fluorescence ∝ ms=0 population
    # (ms=0 is brighter due to spin-dependent ISC)
    pop_ms0 = np.trace(projector_ms(0) @ rho_ss).real
    fluorescence.append(pop_ms0)

# Normalize and invert (ODMR shows dips)
fluorescence = np.array(fluorescence)
odmr_signal = 1 - 0.3 * (1 - fluorescence)  # 30% contrast

# Plot
plt.figure(figsize=(10, 5))
plt.plot(freqs, odmr_signal)
plt.xlabel("Microwave Frequency (GHz)")
plt.ylabel("Fluorescence (a.u.)")
plt.title("CW-ODMR Spectrum")

# Mark expected resonances
gamma_mhz = 28  # MHz/mT
f_minus = D - gamma_mhz * B / 1000
f_plus = D + gamma_mhz * B / 1000
plt.axvline(f_minus, color='r', linestyle='--', label=f'ms=-1: {f_minus:.3f} GHz')
plt.axvline(f_plus, color='b', linestyle='--', label=f'ms=+1: {f_plus:.3f} GHz')
plt.legend()
plt.show()

Hyperfine Structure

Include N14 hyperfine to see the triplet structure:

from sim.hamiltonian.terms import HyperfineN14

freqs = np.linspace(2.55, 2.65, 300)  # Zoom on ms=-1 transition
fluorescence = []

for f in freqs:
    H = HamiltonianBuilder()
    H.add(ZFS(D=2.87))
    H.add(Zeeman(B=10))
    H.add(HyperfineN14())  # Adds triplet structure
    H.add(MicrowaveDrive(omega=0.5, detuning=f - 2.87))

    solver = LindbladSolver(H)
    solver.add_t2_dephasing(gamma=2e6)  # Narrower linewidth

    rho_ss = solver.steady_state(t_max=2e-5)
    pop_ms0 = np.trace(projector_ms(0) @ rho_ss).real
    fluorescence.append(pop_ms0)

plt.figure(figsize=(10, 5))
plt.plot(freqs * 1000, fluorescence)  # Convert to MHz
plt.xlabel("Microwave Frequency (MHz, relative to D)")
plt.ylabel("ms=0 Population")
plt.title("ODMR with N14 Hyperfine Structure")
plt.show()

You should see three dips separated by ~2.2 MHz (the hyperfine splitting).

Pulsed ODMR

In pulsed ODMR, we apply discrete pulses:

from sim.states import ground_state

freqs = np.linspace(2.55, 2.65, 100)
contrast = []

for f in freqs:
    # 1. Initialize in ms=0
    rho = ground_state()

    # 2. Apply MW π-pulse
    H = HamiltonianBuilder()
    H.add(ZFS(D=2.87))
    H.add(Zeeman(B=10))
    H.add(MicrowaveDrive(omega=10, detuning=f - 2.87))

    solver = LindbladSolver(H)
    solver.add_t2_dephasing(gamma=1e6)

    # π-pulse duration
    t_pi = 50e-9  # 50 ns for 10 MHz Rabi
    result = solver.evolve(rho, (0, t_pi), n_steps=20)
    rho = result.final_state()

    # 3. Readout
    pop_ms0 = np.trace(projector_ms(0) @ rho).real
    contrast.append(1 - pop_ms0)  # Population transfer

plt.figure(figsize=(10, 5))
plt.plot(freqs * 1000, contrast)
plt.xlabel("Microwave Frequency (MHz)")
plt.ylabel("Population Transfer")
plt.title("Pulsed ODMR")
plt.show()

Magnetic Field Sensing

Use ODMR to determine magnetic field:

# Unknown field
B_unknown = 15.3  # mT

# Simulate ODMR
freqs = np.linspace(2.4, 3.4, 500)
signal = []

for f in freqs:
    H = HamiltonianBuilder()
    H.add(ZFS(D=2.87))
    H.add(Zeeman(B=B_unknown))
    H.add(MicrowaveDrive(omega=1, detuning=f - 2.87))

    solver = LindbladSolver(H)
    solver.add_t2_dephasing(gamma=5e6)
    rho_ss = solver.steady_state()
    signal.append(np.trace(projector_ms(0) @ rho_ss).real)

signal = np.array(signal)

# Find dips (minima)
from scipy.signal import find_peaks
peaks, _ = find_peaks(-signal, height=-0.9)

f_dips = freqs[peaks]
print(f"Detected resonances: {f_dips} GHz")

# Calculate field
if len(f_dips) >= 2:
    splitting = abs(f_dips[1] - f_dips[0])  # GHz
    B_measured = splitting / 0.056  # 56 MHz/mT = 0.056 GHz/mT
    print(f"Measured field: {B_measured:.2f} mT")
    print(f"Actual field: {B_unknown:.2f} mT")