Skip to content

Quick Start

This tutorial walks you through the most common workflows in SpinDefectSim, from a single spin to a full sensing experiment with an inhomogeneous defect ensemble.


1. Single-spin ODMR

Pick a defect and apply a bias field

import numpy as np
import SpinDefectSim as sds

# VB⁻ in hBN at 2 mT bias field (along z by default)
defect = sds.SpinDefect("vb_minus", B_mT=2.0)

Compute ODMR transition frequencies

# At zero electric field
freqs = defect.transition_frequencies()
print(f"Transitions: {freqs / 1e9} GHz")
# → array([3.403..., 3.517...]) GHz

# With an applied electric field (V/m)
freqs_E = defect.transition_frequencies(E_vec_Vpm=(1e6, 0, 0))

Compute the spin Hamiltonian matrix

H = defect.hamiltonian(E_vec_Vpm=(0, 0, 0))
# H is a (2S+1 × 2S+1) complex128 array in Hz units
eigenvalues, eigenvectors = defect.diagonalize()

Include nuclear spins (hyperfine structure)

H_full = defect.full_hamiltonian()
# Full Hilbert space: electron ⊗ 3 × ¹⁴N for VB⁻
hf_freqs = defect.hyperfine_transitions()
print(f"Hyperfine transitions: {hf_freqs / 1e6} MHz")

2. CW ODMR spectrum

Single-defect spectrum

from SpinDefectSim.spin.spectra import PL_model

f_axis = np.linspace(3.3e9, 3.6e9, 500)
pl = PL_model(
    f_MW_Hz=f_axis,
    transitions_Hz=freqs,
    linewidth_Hz=5e6,   # 5 MHz FWHM
    contrast=0.02,      # 2% CW contrast typical for VB⁻
)

import matplotlib.pyplot as plt
plt.plot(f_axis / 1e9, pl)
plt.xlabel("MW frequency (GHz)")
plt.ylabel("Normalised PL")
plt.title("CW ODMR — single VB⁻ defect")
plt.show()

Ensemble spectrum with inhomogeneous broadening

from SpinDefectSim.spin.spectra import (
    ensemble_transitions_from_Efields,
    ensemble_odmr_spectrum,
)

# 200 defects with random local E-fields (disorder)
E_fields = np.random.randn(200, 3) * 0.5e6   # ±0.5 MV/m spread

sp = defect.spin_params
transitions_list = ensemble_transitions_from_Efields(E_fields, sp)

pl_ens = ensemble_odmr_spectrum(f_axis, transitions_list, fwhm=5e6, contrast=0.02)
plt.plot(f_axis / 1e9, pl_ens)
plt.title("Ensemble ODMR — inhomogeneous broadening")
plt.show()

3. Building an ensemble

Generate defect positions and compute local fields

# Create ensemble: 500 defects within 300 nm radius patch
ens = sds.DefectEnsemble(N_def=500, R_patch=300e-9)
ens.generate_defects(seed=42)

# Uniform gate bias + Yukawa-screened disorder charges
charges = np.array([
    [0, 0, 0, 1.6e-19],    # single electron charge at origin
    [50e-9, 0, 0, -1.6e-19],
])
ens.compute_efields(
    E0_gate=(0, 0, 1e6),   # 1 MV/m uniform z-field
    disorder_xyzq=charges,
)

print(f"E-field range: {ens.E_fields[:, 2].min():.2e}{ens.E_fields[:, 2].max():.2e} V/m")

Random quantization axes

ens.generate_defects(seed=42, quantization_axis="random")
# Each defect now has a random orientation drawn from the sphere

4. Sensing experiments

Wrap the ensemble into a sensing experiment

# "E" → modulate an E-field; "B" → modulate a B-field; "both" → both
exp = ens.to_experiment(sensing="E")

Ramsey (free-induction decay)

tau_s, S_w, S_no, dS, tau_opt, dS_peak = exp.ramsey()

plt.plot(tau_s * 1e9, dS)
plt.axvline(tau_opt * 1e9, ls="--", label=f"τ_opt = {tau_opt*1e9:.0f} ns")
plt.xlabel("Free-precession time τ (ns)")
plt.ylabel("Lock-in signal ΔS")
plt.legend()
plt.title("Ramsey — differential lock-in signal")
plt.show()

Hahn-echo (static-field rejection)

tau_s, S_w, S_no, dS, tau_opt, dS_peak = exp.echo_static()

plt.plot(tau_s * 1e6, dS)
plt.xlabel("τ (µs)")
plt.ylabel("ΔS (echo)")
plt.title("Hahn-echo lock-in")
plt.show()

CW ODMR from the experiment

f_axis = np.linspace(3.35e9, 3.57e9, 400)
pl_w, pl_no, dpl = exp.cw_odmr(f_axis)

fig, axes = plt.subplots(1, 2, figsize=(10, 4))
axes[0].plot(f_axis / 1e9, pl_w, label="with signal")
axes[0].plot(f_axis / 1e9, pl_no, label="reference")
axes[0].legend()
axes[0].set_title("CW ODMR spectra")
axes[1].plot(f_axis / 1e9, dpl)
axes[1].set_title("Lock-in differential")
for ax in axes:
    ax.set_xlabel("MW frequency (GHz)")
plt.tight_layout()
plt.show()

5. SNR and integration time

from SpinDefectSim.sensing.snr import snr, n_avg_for_threshold
from SpinDefectSim.sensing.sequences import RamseySequence

# How many averages to reach SNR = 5?
N_needed = exp.n_avg_to_detect(dS_peak, snr_target=5.0)
print(f"Need {N_needed:.0f} averages")

# Convert to wall-clock time with a Ramsey sequence
seq = RamseySequence()
T_per_shot = seq.total_time(tau_opt)
T_total_s = N_needed * T_per_shot
print(f"Integration time ≈ {T_total_s:.3f} s  ({T_total_s * 1e3:.1f} ms)")

6. Magnetometry

from SpinDefectSim.magnetometry.geometry import DiskGeometry
from SpinDefectSim.magnetometry.magnetometry import MagnetometryExperiment

# Disk-shaped magnetic sample, radius 1 µm
geom = DiskGeometry(radius=1e-6)

# Uniform in-plane magnetization (A, equivalent to surface current)
magnetization = np.ones((100, 100)) * 1e-3   # 1 mA equivalent

mag_exp = MagnetometryExperiment(
    geometry=geom,
    magnetization=magnetization,
    z_defect=50e-9,   # defect 50 nm above sample
)

# Stray B-field map
x_obs = np.linspace(-2e-6, 2e-6, 50)
y_obs = np.linspace(-2e-6, 2e-6, 50)
Bz = mag_exp.B_z_map(x_obs, y_obs)

plt.imshow(Bz * 1e3, extent=[-2, 2, -2, 2], origin="lower", cmap="RdBu_r")
plt.colorbar(label="Bz (mT)")
plt.xlabel("x (µm)")
plt.ylabel("y (µm)")
plt.title("Stray field map above magnetic disk")
plt.show()

7. Electrometry

from SpinDefectSim.electrometry.electrometry import ElectrometryExperiment

# Point charge (e.g. a single electron) trapped at the surface
charges = np.array([[0.0, 0.0, 0.0, 1.6e-19]])  # [x, y, z, q]

elec_exp = ElectrometryExperiment(
    charges_xyzq=charges,
    z_defect=0.34e-9,
    E0_gate=(0, 0, 0),
    screening_model="dual_gate",
    d_gate=15e-9,
)

# Map of ODMR frequency shift across a 200 nm field of view
x_obs = np.linspace(-100e-9, 100e-9, 40)
y_obs = np.linspace(-100e-9, 100e-9, 40)
freq_map = elec_exp.transition_frequency_map(x_obs, y_obs)

# Plot the lower transition frequency
plt.imshow(freq_map[:, :, 0] / 1e9, extent=[-100, 100, -100, 100], origin="lower")
plt.colorbar(label="Transition frequency (GHz)")
plt.xlabel("x (nm)")
plt.ylabel("y (nm)")
plt.title("Electrometry — transition frequency map")
plt.show()

8. Custom defect types

from SpinDefectSim.spin.defects import DefectType
from SpinDefectSim.spin.nuclear import NuclearSpin, axial_A_tensor, GAMMA_14N

# Define a custom spin-1 defect
my_defect_type = DefectType(
    name="my_centre",
    spin=1,
    D0_Hz=2.0e9,
    E0_Hz=30e6,
    d_perp=0.25,
    gamma_Hz_T=28e9,
    nuclear_spins=[
        NuclearSpin(
            spin=1,
            A_tensor_Hz=axial_A_tensor(A_zz_Hz=-2e6, A_perp_Hz=-2.5e6),
            gamma_Hz_T=GAMMA_14N,
            label="14N",
            quadrupole_Hz=-4.95e6,
        )
    ],
    notes="My custom spin-1 defect with one 14N",
)

defect = sds.SpinDefect(my_defect_type, B_mT=5.0)
print(defect.transition_frequencies() / 1e9)

9. Parameter sweeps

from SpinDefectSim.analysis.sweep import ParameterSweep

sweep = ParameterSweep(N_def=300, seed=0)

# Sweep over gate field strength and defect density
def run_one(E_gate_Vpm, N_def):
    ens = sweep.make_ensemble(N_def=N_def)
    ens.generate_defects()
    ens.compute_efields(E0_gate=(0, 0, E_gate_Vpm))
    exp = ens.to_experiment(sensing="E")
    _, _, _, dS, tau_opt, dS_peak = exp.ramsey()
    return {"dS_peak": dS_peak, "tau_opt_ns": tau_opt * 1e9}

results = sweep.sweep(
    run_one,
    E_gate_Vpm=[1e6, 2e6, 5e6],
    N_def=[100, 300, 600],
)
# results is a list of dicts with every combination
for r in results:
    print(r)

Next steps