Hamiltonian Module¶

The Hamiltonian module provides a modular builder pattern for constructing NV center Hamiltonians from individual physical terms.

Hamiltonian Builder¶

Hamiltonian builder for NV center simulation.

The builder combines multiple Hamiltonian terms into a total Hamiltonian and provides analysis methods.

Architecture¶

The builder follows the Composite pattern:

HamiltonianBuilder ├── ZFS(D=2.87) ├── Zeeman(B=10) ├── Hyperfine(…) └── …

H_total(t) = Σ term.build(t)

Units¶

Input: Practical units (GHz, mT, …)
Output: The matrix H is in rad/s (energy with ℏ=1)
Eigenvalues can be retrieved in GHz via eigenvalues_ghz()

Example

>>> from sim import HamiltonianBuilder
>>> from sim.hamiltonian.terms import ZFS, Zeeman
>>>
>>> H = HamiltonianBuilder()
>>> H.add(ZFS(D=2.87))        # 2.87 GHz ZFS
>>> H.add(Zeeman(B=10))       # 10 mT magnetic field
>>>
>>> # Matrix at t=0
>>> H_matrix = H.build(t=0.0)
>>> print(H_matrix.shape)
(18, 18)
>>>
>>> # Eigenvalues in GHz
>>> eigvals = H.eigenvalues_ghz()
>>> print(f"Levels: {np.unique(np.round(eigvals, 2))}")

class sim.hamiltonian.builder.HamiltonianBuilder[source]¶

Bases: object

Builder for constructing the NV center Hamiltonian.

Collects individual terms (ZFS, Zeeman, Hyperfine, …) and combines them into the total matrix.

terms¶

List of added HamiltonianTerm objects

Type:: list

add(term)[source]¶

Add a term

Parameters:: term (HamiltonianTerm)
Return type:: HamiltonianBuilder

build(t=0.0)[source]¶

Build the 18×18 matrix at time t

Parameters:: t (float)
Return type:: ndarray

eigenvalues(t=0.0)[source]¶

Compute eigenvalues in rad/s

Parameters:: t (float)
Return type:: ndarray

eigenvalues_ghz(t=0.0)[source]¶

Compute eigenvalues in GHz

Parameters:: t (float)
Return type:: ndarray

Examples

Simple NV center with ZFS and magnetic field:

>>> H = HamiltonianBuilder()
>>> H.add(ZFS(D=2.87))
>>> H.add(Zeeman(B=10))
>>> print(len(H))
2
>>> print(H)
HamiltonianBuilder([ZFS, Zeeman])

Analyze eigenvalues:

>>> eigvals = H.eigenvalues_ghz()
>>> print(f"Ground state: {min(eigvals):.3f} GHz")
>>> print(f"Highest state: {max(eigvals):.3f} GHz")

Method chaining:

>>> H = HamiltonianBuilder().add(ZFS()).add(Zeeman(B=5))

DIM = 18¶

__init__()[source]¶: Initialize an empty builder.

add(term)[source]¶

Add a Hamiltonian term.

Parameters:: term (HamiltonianTerm) – A term object (e.g., ZFS, Zeeman, Hyperfine)
Returns:: self, for method chaining
Return type:: HamiltonianBuilder
Raises:: TypeError – If term is not a HamiltonianTerm

Examples

>>> H = HamiltonianBuilder()
>>> H.add(ZFS(D=2.87)).add(Zeeman(B=10))
HamiltonianBuilder([ZFS, Zeeman])

clear()[source]¶

Remove all terms.

Returns:: self, for method chaining
Return type:: HamiltonianBuilder

property terms: List[HamiltonianTerm]¶

List of added terms (copy).

Returns:: Copy of the term list
Return type:: list

property is_time_dependent: bool¶

Whether any term is time-dependent.

Returns:: True if at least one term is time-dependent
Return type:: bool

build(t=0.0)[source]¶

Build the total Hamiltonian matrix at time t.

H_total(t) = Σ_i term_i.build(t)

Parameters:: t (float) – Time in seconds
Returns:: 18×18 complex Hermitian matrix in rad/s
Return type:: np.ndarray

Notes

The unit rad/s corresponds to energy with ℏ=1. To convert to Hz: f = ω / (2π) To GHz: f_GHz = ω / (2π × 10⁹)

eigenvalues(t=0.0)[source]¶

Compute eigenvalues of the Hamiltonian.

Parameters:: t (float) – Time in seconds
Returns:: 18 eigenvalues in rad/s, sorted ascending
Return type:: np.ndarray

Notes

Uses np.linalg.eigvalsh (for Hermitian matrices). The eigenvalues are real.

eigenvalues_ghz(t=0.0)[source]¶

Compute eigenvalues in GHz.

Parameters:: t (float) – Time in seconds
Returns:: 18 eigenvalues in GHz, sorted ascending
Return type:: np.ndarray

Examples

>>> H = HamiltonianBuilder().add(ZFS(D=2.87))
>>> eigvals = H.eigenvalues_ghz()
>>> np.unique(np.round(eigvals, 2))
array([0.  , 2.87])

eigenvalues_mhz(t=0.0)[source]¶

Compute eigenvalues in MHz.

Parameters:: t (float) – Time in seconds
Returns:: 18 eigenvalues in MHz, sorted ascending
Return type:: np.ndarray

eigenstates(t=0.0)[source]¶

Compute eigenvalues and eigenvectors.

Parameters:

t (float) – Time in seconds

Returns:

eigenvalues (np.ndarray) – 18 eigenvalues in rad/s
eigenvectors (np.ndarray) – 18×18 matrix, columns are eigenvectors

Return type:

tuple

Usage¶

from sim import HamiltonianBuilder
from sim.hamiltonian.terms import ZFS, Zeeman, HyperfineN14

# Create builder
H = HamiltonianBuilder()

# Add terms
H.add(ZFS(D=2.87))           # Zero-field splitting
H.add(Zeeman(B=10))          # 10 mT magnetic field
H.add(HyperfineN14())        # N14 hyperfine structure

# Build 18×18 matrix at time t
H_matrix = H.build(t=0)

# Check if time-dependent
print(H.is_time_dependent)   # False for static terms

# Get eigenvalues
eigvals = H.eigenvalues_ghz()

Hamiltonian Terms¶

All Hamiltonian terms inherit from HamiltonianTerm and implement a build(t) method returning an 18×18 Hermitian matrix.

Zero-Field Splitting¶

Zero-Field Splitting (ZFS) term for the NV center Hamiltonian.

Physical Background¶

The NV center has two unpaired electrons forming a triplet state (S=1). The spin-spin interaction between these electrons leads to Zero-Field Splitting (ZFS) - an energy splitting even without an external magnetic field.

The ZFS Hamiltonian is:

H_ZFS = D·Sz² + E·(Sx² - Sy²)

where:

D ≈ 2.87 GHz: Axial parameter (splitting ms=0 ↔ ms=±1)
E: Transverse parameter (strain-induced, lifts ms=±1 degeneracy)

Energy Levels¶

Without E-term (E=0):

E(ms=0) = 0
E(ms=±1) = D = 2.87 GHz (degenerate)

With E-term (E≠0):

E(ms=0) = 0
E(ms=+1) = D + E
E(ms=-1) = D - E

→ Splitting of 2E between ms=±1

Typical Values¶

D_GS = 2.87 GHz (ground state at room temperature)

D_ES = 1.42 GHz (excited state)

E = 0 to ~10 MHz (depends on strain/defects)

Temperature Dependence¶

D(T) ≈ D₀ + α·(T - 300K) with α ≈ -74 kHz/K

References

[1] Doherty et al., Physics Reports 528, 1-45 (2013) [2] Maze et al., New J. Phys. 13, 025025 (2011)

class sim.hamiltonian.terms.zfs.ZFS(D=2.87, E=0.0, name=None)[source]¶

Bases: HamiltonianTerm

Zero-Field Splitting Hamiltonian term.

Implements H_ZFS = D·Sz² + E·(Sx² - Sy²)

Parameters:

D (float) – Axial ZFS parameter in GHz. Default: 2.87 GHz (NV ground state at room temperature)
E (float) – Transverse ZFS parameter in GHz. Default: 0 (no strain splitting)
name (str, optional) – Custom name for this term

D¶

Axial parameter in GHz

Type:: float

E¶

Transverse parameter in GHz

Type:: float

Examples

Standard NV center (no strain):

>>> zfs = ZFS(D=2.87)
>>> H = zfs.build()
>>> eigvals = np.linalg.eigvalsh(H) / (2*np.pi*1e9)  # in GHz
>>> np.unique(np.round(eigvals, 2))
array([0.  , 2.87])

With strain splitting (E=5 MHz):

>>> zfs = ZFS(D=2.87, E=0.005)
>>> H = zfs.build()
>>> eigvals = np.linalg.eigvalsh(H) / (2*np.pi*1e9)
>>> # ms=±1 are now split at 2.865 and 2.875 GHz

Notes

The matrix is 18×18, but ZFS only acts on the electron spin. The 6-fold degeneracy at E(ms=0)=0 comes from:

2× g/e states

3× nuclear spin mI states

The 12-fold degeneracy at E(ms=±1)=D comes from:

2× ms=+1 and ms=-1
2× g/e states
3× nuclear spin mI states

__init__(D=2.87, E=0.0, name=None)[source]¶

Initialize the Hamiltonian term.

Parameters:

name (str, optional) – Custom name. If not provided, the class name is used.
D (float)
E (float)

property is_time_dependent: bool¶: ZFS is time-independent (without external modulation).

build(t=0.0)[source]¶

Build the 18×18 ZFS matrix.

Parameters:: t (float) – Time in seconds (ignored, since time-independent)
Returns:: 18×18 Hermitian matrix in rad/s
Return type:: np.ndarray

from sim.hamiltonian.terms import ZFS

# Standard NV center
zfs = ZFS(D=2.87)  # GHz

# With strain splitting
zfs = ZFS(D=2.87, E=0.005)  # 5 MHz E-term

Zeeman Effect¶

Zeeman term for the NV center Hamiltonian.

Physical Background¶

The Zeeman effect describes the interaction of a magnetic moment with an external magnetic field. For the electron spin:

H_Zeeman = γ_e · (B_x·S_x + B_y·S_y + B_z·S_z)
= γ_e · B⃗ · S⃗

where:

γ_e = g_e · μ_B / ℏ ≈ 1.76×10¹¹ rad/(s·T): Gyromagnetic ratio
g_e ≈ 2.0028: Electron g-factor
μ_B: Bohr magneton
B⃗: Magnetic field vector
S⃗: Spin operator vector

Energy Levels¶

For a field B along the z-axis (NV axis):

E(ms) = γ_e · B_z · ms

The splitting between ms=+1 and ms=-1 is:: ΔE = 2·γ_e·B_z = 2·g_e·μ_B·B_z/ℏ
In practical units:: Δf ≈ 28 MHz/mT · B[mT]

Example: At B=10 mT, Δf ≈ 280 MHz between ms=+1 and ms=-1.

Interplay with ZFS¶

In the full picture with ZFS (D=2.87 GHz):

At B=0: ms=0 at 0 GHz, ms=±1 at 2.87 GHz (degenerate)
At B>0: ms=±1 split to 2.87 ± γB

This allows selective excitation of ms=0↔ms=+1 or ms=0↔ms=-1 at different microwave frequencies.

ODMR (Optically Detected Magnetic Resonance)¶

The Zeeman splitting is the basis for NV magnetometry:

ODMR spectrum shows dips at 2.87 ± γB GHz
From the splitting, B can be determined

References

[1] Doherty et al., Physics Reports 528, 1-45 (2013) [2] Schirhagl et al., Annu. Rev. Phys. Chem. 65, 83-105 (2014)

class sim.hamiltonian.terms.zeeman.Zeeman(B=0.0, g=2.0028, name=None)[source]¶

Bases: HamiltonianTerm

Zeeman Hamiltonian for the electron spin in a magnetic field.

Implements H_Zeeman = γ_e · (B_x·S_x + B_y·S_y + B_z·S_z)

Parameters:

B (float, list, tuple, or np.ndarray) – Magnetic field vector [B_x, B_y, B_z] in mT. If a single float is given, B_z=B, B_x=B_y=0 is assumed.
g (float) – Electron g-factor. Default: 2.0028 (NV center)
name (str, optional) – Custom name

B¶

Magnetic field vector [B_x, B_y, B_z] in mT

Type:: np.ndarray

g¶

g-factor

Type:: float

Examples

Field along the NV axis (z):

>>> zeeman = Zeeman(B=10)  # 10 mT in z-direction
>>> zeeman = Zeeman(B=[0, 0, 10])  # equivalent
>>> print(zeeman.splitting_mhz())
560.6  # MHz splitting between ms=+1 and ms=-1

Arbitrary field direction:

>>> zeeman = Zeeman(B=[5, 0, 10])  # tilted field
>>> H = zeeman.build()

Together with ZFS:

>>> from sim import HamiltonianBuilder
>>> from sim.hamiltonian.terms import ZFS, Zeeman
>>> H = HamiltonianBuilder()
>>> H.add(ZFS(D=2.87))
>>> H.add(Zeeman(B=10))
>>> eigvals = H.eigenvalues_ghz()
>>> # Shows splitting of ms=±1 levels

Notes

The typical Zeeman splitting is ~28 MHz/mT.

More precisely: Δf = 2·g·μ_B·B/h ≈ 2·2.0028·9.274e-24·B / 6.626e-34: ≈ 56.06 MHz/mT · B[mT]

I.e., the splitting between ms=+1 and ms=-1 is ~56 MHz per mT. Per ms level (e.g., ms=0 to ms=+1) it is ~28 MHz/mT.

__init__(B=0.0, g=2.0028, name=None)[source]¶

Initialize the Hamiltonian term.

Parameters:

name (str, optional) – Custom name. If not provided, the class name is used.
B (float | List[float] | Tuple[float, ...] | ndarray)
g (float)

property is_time_dependent: bool¶: Static Zeeman field is time-independent.

build(t=0.0)[source]¶

Build the 18×18 Zeeman matrix.

Parameters:: t (float) – Time in seconds (ignored for static field)
Returns:: 18×18 Hermitian matrix in rad/s
Return type:: np.ndarray

splitting_mhz()[source]¶

Compute the Zeeman splitting in MHz.

The splitting is defined as the energy difference between ms=+1 and ms=-1 for the B_z field.

Δf = 2·γ·B_z / (2π)

Returns:: Splitting in MHz
Return type:: float

Examples

>>> Zeeman(B=10).splitting_mhz()
560.6...  # ~56 MHz/mT × 10 mT

larmor_frequency_mhz()[source]¶

Compute the Larmor frequency in MHz.

The Larmor frequency is the precession frequency of the spin in the magnetic field: f_L = γ·B / (2π)

Returns:: Larmor frequency for |B| in MHz
Return type:: float

from sim.hamiltonian.terms import Zeeman

# Field along NV axis
zeeman = Zeeman(B=10)  # 10 mT

# Arbitrary field direction
zeeman = Zeeman(B=[5, 0, 10])  # mT

# Get splitting
print(zeeman.splitting_mhz())  # ~560 MHz

N14 Hyperfine¶

N14 hyperfine coupling for NV center.

The nitrogen-14 nucleus (I=1) at the vacancy site couples to the electron spin through hyperfine interaction and has a nuclear quadrupole moment.

Physics¶

The N14 hyperfine Hamiltonian consists of two parts:

Hyperfine coupling (electron-nuclear spin interaction):

H_hf = A_∥ Sz Iz + A_⊥ (Sx Ix + Sy Iy)

where A_∥ ≈ -2.14 MHz and A_⊥ ≈ -2.70 MHz for NV centers.
Quadrupole coupling (nuclear electric quadrupole):

H_Q = P (Iz² - 2/3)

where P ≈ -5.0 MHz is the quadrupole parameter.

The total Hamiltonian is:

H_N14 = A_∥ Sz Iz + A_⊥ (Sx Ix + Sy Iy) + P (Iz² - 2/3)

Matrix Dimension¶

The hyperfine coupling acts on the space |ms⟩ ⊗ |mI⟩ which has dimension 3 × 3 = 9. This is then extended to the full 18-dimensional space by tensor product with the electronic g/e identity.

Literature Values¶

A_∥ = -2.14 MHz (parallel hyperfine)
A_⊥ = -2.70 MHz (perpendicular hyperfine)
P = -5.0 MHz (quadrupole parameter)

References

[1] Doherty et al., Physics Reports 528, 1-45 (2013) [2] Smeltzer et al., New J. Phys. 13, 025021 (2011)

class sim.hamiltonian.terms.hyperfine_n14.HyperfineN14(A_parallel=-2.14, A_perp=-2.7, P_quad=-5.0, name=None)[source]¶

Bases: HamiltonianTerm

N14 hyperfine coupling with quadrupole interaction.

Parameters:

A_parallel (float) – Parallel hyperfine coupling A_∥ in MHz. Default: -2.14 MHz
A_perp (float) – Perpendicular hyperfine coupling A_⊥ in MHz. Default: -2.70 MHz
P_quad (float) – Quadrupole parameter P in MHz. Default: -5.0 MHz
name (str, optional) – Custom name for the term

A_parallel¶

Parallel coupling in MHz

Type:: float

A_perp¶

Perpendicular coupling in MHz

Type:: float

P_quad¶

Quadrupole parameter in MHz

Type:: float

Examples

>>> from sim import HamiltonianBuilder
>>> from sim.hamiltonian.terms import ZFS, HyperfineN14
>>>
>>> H = HamiltonianBuilder()
>>> H.add(ZFS(D=2.87))
>>> H.add(HyperfineN14())  # Default N14 parameters
>>>
>>> # Check eigenvalues show hyperfine structure
>>> eigvals = H.eigenvalues_mhz()
>>> print(f"Hyperfine splitting visible in eigenvalues")

Custom parameters:

>>> hf = HyperfineN14(A_parallel=-2.2, A_perp=-2.8, P_quad=-4.9)
>>> H_matrix = hf.build()
>>> print(f"Matrix shape: {H_matrix.shape}")
(18, 18)

Notes

The hyperfine structure splits each ms level into three mI levels. The quadrupole term causes asymmetric splitting of the mI = ±1 levels.

The sign convention follows the standard where negative values indicate antiferromagnetic coupling.

__init__(A_parallel=-2.14, A_perp=-2.7, P_quad=-5.0, name=None)[source]¶

Initialize the Hamiltonian term.

Parameters:

name (str, optional) – Custom name. If not provided, the class name is used.
A_parallel (float)
A_perp (float)
P_quad (float)

property is_time_dependent: bool¶: N14 hyperfine coupling is time-independent.

build(t=0.0)[source]¶

Build the 18×18 N14 hyperfine Hamiltonian.

Parameters:: t (float) – Time in seconds (not used, hyperfine is static)
Returns:: 18×18 complex Hermitian matrix in rad/s
Return type:: np.ndarray

Notes

The result is cached since the Hamiltonian is time-independent. The 9×9 matrix is extended to 18×18 by acting identically on both ground and excited electronic states.

hyperfine_tensor()[source]¶

Return the 3×3 hyperfine tensor in the NV frame.

Returns:: 3×3 diagonal tensor with [A_⊥, A_⊥, A_∥] in MHz
Return type:: np.ndarray

Notes

The NV axis is along z, so: - A_xx = A_yy = A_⊥ (perpendicular) - A_zz = A_∥ (parallel)

splitting_mhz(ms=0)[source]¶

Calculate the N14 hyperfine splitting for a given ms level.

Parameters:: ms (int) – Electron spin projection: +1, 0, or -1
Returns:: Dictionary with keys ‘mI_+1’, ‘mI_0’, ‘mI_-1’ containing the energy shifts in MHz for each nuclear spin state
Return type:: dict

Examples

>>> hf = HyperfineN14()
>>> split = hf.splitting_mhz(ms=0)
>>> print(f"mI=0 shift: {split['mI_0']:.2f} MHz")

from sim.hamiltonian.terms import HyperfineN14

# Default literature values
hf = HyperfineN14()

# Custom parameters
hf = HyperfineN14(
    A_parallel=-2.14,  # MHz
    A_perp=-2.70,      # MHz
    P_quad=-5.0        # MHz
)

C13 Hyperfine¶

C13 hyperfine coupling for NV center.

Carbon-13 nuclei (I=1/2) in the diamond lattice couple to the NV electron spin through dipolar hyperfine interaction. The coupling strength depends on the position of the C13 nucleus relative to the NV center.

Physics¶

The C13 hyperfine Hamiltonian is:

H = S · A · I

where A is the hyperfine tensor, S is the electron spin (S=1), and I is the C13 nuclear spin (I=1/2).

The hyperfine tensor has two contributions:

Isotropic (Fermi contact): A_iso = (2μ₀/3) γₑ γ_C ℏ² |ψ(r)|²

This is typically small for C13 in diamond (not at NV site).
Dipolar (anisotropic): A_dip = (μ₀/4π) (γₑ γ_C ℏ²/r³) (I - 3r̂r̂ᵀ)

This is the dominant contribution for C13 in the lattice.

Matrix Structure¶

The C13 hyperfine couples electron spin (3D) to C13 nuclear spin (2D), giving a 6×6 matrix. This is embedded in the 18D space as block-diagonal over the N14 nuclear spin states.

Hilbert space structure:: |ψ⟩ = |g/e⟩ ⊗ |ms⟩ ⊗ |mI_N14⟩ ⊗ |mI_C13⟩

For computational efficiency within the standard 18D basis, we treat C13 coupling in a block-diagonal approximation.

Typical Coupling Strengths¶

Nearest neighbor C13 (~1.5 Å): A_zz ~ 100-200 kHz
Second neighbor (~2.5 Å): A_zz ~ 30-50 kHz
Further away (>5 Å): A_zz < 10 kHz

References

[1] Childress et al., Science 314, 281 (2006) [2] Maze et al., Nature 455, 644 (2008)

class sim.hamiltonian.terms.hyperfine_c13.HyperfineC13(r_vec=None, A_iso=0.0, A_tensor=None, name=None)[source]¶

Bases: HamiltonianTerm

C13 hyperfine coupling (dipolar + isotropic).

Parameters:

r_vec (array_like) – Position vector of C13 nucleus relative to NV center, in nm. Should be [rx, ry, rz] in the NV coordinate frame.
A_iso (float) – Isotropic hyperfine coupling in kHz. Default: 0 (pure dipolar)
A_tensor (array_like, optional) – Custom 3×3 hyperfine tensor in kHz. If provided, r_vec is ignored.
name (str, optional) – Custom name for the term

r_vec¶

Position vector in nm

Type:: np.ndarray or None

A_iso¶

Isotropic coupling in kHz

Type:: float

A_tensor¶

Full 3×3 hyperfine tensor in kHz

Type:: np.ndarray

Examples

>>> from sim import HamiltonianBuilder
>>> from sim.hamiltonian.terms import ZFS, HyperfineC13
>>>
>>> # C13 at 0.25 nm along the NV axis
>>> H = HamiltonianBuilder()
>>> H.add(ZFS(D=2.87))
>>> H.add(HyperfineC13(r_vec=[0, 0, 0.25]))
>>>
>>> # Check the coupling strength
>>> hf_c13 = HyperfineC13(r_vec=[0.154, 0, 0])
>>> print(f"A_zz = {hf_c13.A_zz_khz():.1f} kHz")

Custom tensor:

>>> A = np.array([[10, 0, 0], [0, 10, 0], [0, 0, 50]])  # kHz
>>> hf_c13 = HyperfineC13(A_tensor=A)

Notes

The C13 hyperfine interaction is typically much weaker than the N14 hyperfine (~2 MHz), with typical values in the 10-200 kHz range.

The block-diagonal approximation assumes that the C13-N14 cross-terms are negligible (<0.1%), which is valid for most positions.

__init__(r_vec=None, A_iso=0.0, A_tensor=None, name=None)[source]¶

Initialize the Hamiltonian term.

Parameters:

name (str, optional) – Custom name. If not provided, the class name is used.
r_vec (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None)
A_iso (float)
A_tensor (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None)

property is_time_dependent: bool¶: C13 hyperfine coupling is time-independent.

build(t=0.0)[source]¶

Build the 18×18 C13 hyperfine Hamiltonian.

The 6×6 matrix is embedded in the 18D space as block-diagonal over the N14 nuclear spin states (mI_N14 = +1, 0, -1).

Parameters:: t (float) – Time in seconds (not used, hyperfine is static)
Returns:: 18×18 complex Hermitian matrix in rad/s
Return type:: np.ndarray

Notes

The structure is:: H_18 = I_2 ⊗ (H_6 ⊕ H_6 ⊕ H_6) for N14 mI = +1, 0, -1

But since we’re in the 18D basis without explicit C13, we extend appropriately.

A_zz_khz()[source]¶

Return the secular (A_zz) component in kHz.

Returns:: A_zz component of the hyperfine tensor in kHz
Return type:: float

distance_nm()[source]¶

Return the C13-NV distance in nm.

Returns:: Distance in nm, or None if A_tensor was provided directly
Return type:: float or None

from sim.hamiltonian.terms import HyperfineC13

# C13 at specific position (in Angstroms)
hf_c13 = HyperfineC13(position=[1.54, 0, 0])

Microwave Drive¶

Microwave drive for NV center spin control.

Microwave fields at ~2.87 GHz drive transitions between the ms=0 and ms=±1 spin states. This is the primary control mechanism for NV spin manipulation.

Physics¶

In the rotating frame (rotating wave approximation), the MW drive Hamiltonian is:

H_MW = (Ω/2) · [S₊ e^(-iφ) + S₋ e^(+iφ)] + Δ · Sz

where: - Ω(t) is the Rabi frequency (proportional to MW amplitude) - φ(t) is the phase of the MW field - Δ is the detuning from resonance - S± = Sx ± i·Sy are the raising/lowering operators

The Rabi frequency is related to the MW magnetic field by:

Ω = γ_e · B_MW

where γ_e ≈ 28 GHz/T is the electron gyromagnetic ratio.

Typical Values¶

Rabi frequency: 1-50 MHz (depends on MW power)
π-pulse duration: 10-500 ns
Detuning: 0 for resonant driving

References

[1] Childress & Hanson, MRS Bulletin 38, 134 (2013) [2] Dréau et al., Phys. Rev. B 84, 195204 (2011)

class sim.hamiltonian.terms.mw_drive.MicrowaveDrive(omega=0.0, phase=0.0, detuning=0.0, manifold='ground', name=None)[source]¶

Bases: HamiltonianTerm

Time-dependent microwave drive in rotating frame.

Parameters:

omega (float or callable) – Rabi frequency in MHz. Can be: - float: constant Rabi frequency - callable: omega(t) returning Rabi frequency at time t
phase (float or callable) – Phase in radians. Can be: - float: constant phase - callable: phase(t) returning phase at time t
detuning (float) – Frequency detuning from resonance in MHz. Default: 0
manifold (str) – Which electronic manifold: “ground”, “excited”, or “both”. Default: “ground”
name (str, optional) – Custom name for the term

omega¶

Rabi frequency parameter

Type:: float or callable

phase¶

Phase parameter

Type:: float or callable

detuning¶

Detuning in MHz

Type:: float

manifold¶

Target manifold

Type:: str

Examples

>>> from sim import HamiltonianBuilder
>>> from sim.hamiltonian.terms import ZFS, MicrowaveDrive
>>>
>>> # Constant drive
>>> H = HamiltonianBuilder()
>>> H.add(ZFS(D=2.87))
>>> H.add(MicrowaveDrive(omega=10, phase=0))  # 10 MHz Rabi frequency
>>>
>>> # Time-dependent Rabi frequency (pulsed)
>>> def rabi_pulse(t):
...     if 0 <= t < 100e-9:  # 100 ns pulse
...         return 10.0  # 10 MHz
...     return 0.0
>>>
>>> H.add(MicrowaveDrive(omega=rabi_pulse))
>>>
>>> # Gaussian pulse shape
>>> def gaussian_rabi(t):
...     t0, sigma = 50e-9, 20e-9
...     return 20 * np.exp(-(t - t0)**2 / (2 * sigma**2))
>>>
>>> H.add(MicrowaveDrive(omega=gaussian_rabi))

Notes

The MW drive is always considered time-dependent even if omega and phase are constants, because it represents an oscillating field.

The rotating wave approximation (RWA) assumes near-resonant driving. For large detuning, the RWA breaks down.

__init__(omega=0.0, phase=0.0, detuning=0.0, manifold='ground', name=None)[source]¶

Initialize the Hamiltonian term.

Parameters:

name (str, optional) – Custom name. If not provided, the class name is used.
omega (float | Callable)
phase (float | Callable)
detuning (float)
manifold (str)

property is_time_dependent: bool¶: MW drive is always time-dependent.

build(t=0.0)[source]¶

Build the 18×18 MW drive Hamiltonian at time t.

Parameters:: t (float) – Time in seconds
Returns:: 18×18 complex matrix in rad/s
Return type:: np.ndarray

Notes

The MW drive typically only acts on the ground state manifold in standard ODMR experiments.

pi_time()[source]¶

Calculate π-pulse duration for constant Rabi frequency.

Returns:: π-pulse duration in seconds, or None if Rabi freq is time-dependent
Return type:: float or None

from sim.hamiltonian.terms import MicrowaveDrive

# Constant drive (10 MHz Rabi frequency)
mw = MicrowaveDrive(omega=10, detuning=0, phase=0)

# Time-dependent Rabi frequency
def rabi(t):
    return 10 if 0 < t < 50e-9 else 0

mw = MicrowaveDrive(omega=rabi)

Optical Drive¶

Optical coupling between ground and excited states.

Laser light couples the NV ground state (³A₂) to the excited state (³E) through electric dipole transitions. This is used for optical initialization, readout, and coherent optical control.

Physics¶

The optical Hamiltonian in the rotating frame is:

H_opt = (Ω_L/2) · [|g⟩⟨e| + |e⟩⟨g|]

where Ω_L is the optical Rabi frequency, proportional to the laser electric field amplitude.

For spin-conserving transitions (Δms = 0):

H_opt = (Ω_L/2) · Σ_{ms,mI} [|g,ms,mI⟩⟨e,ms,mI| + h.c.]

The optical Rabi frequency is related to laser parameters by:

Ω_L = d · E_0 / ℏ

where d is the dipole moment and E_0 is the electric field amplitude.

Typical Values¶

Zero-phonon line: 637 nm (1.945 eV)
Green excitation: 532 nm
Dipole moment: d ≈ 1.6 × 10⁻²⁹ C·m
Typical Rabi frequencies: 1-100 MHz

References

[1] Doherty et al., Physics Reports 528, 1-45 (2013) [2] Robledo et al., Nature 477, 574 (2011)

class sim.hamiltonian.terms.optical.OpticalCoupling(omega=0.0, detuning=0.0, spin_conserving=True, name=None)[source]¶

Bases: HamiltonianTerm

Optical coupling between ground and excited states.

Parameters:

omega (float or callable) – Optical Rabi frequency in MHz. Can be: - float: constant Rabi frequency - callable: omega(t) returning Rabi frequency at time t
detuning (float) – Detuning from the optical transition in GHz. Default: 0
spin_conserving (bool) – If True, only spin-conserving transitions (Δms = 0). Default: True
name (str, optional) – Custom name for the term

omega¶

Optical Rabi frequency parameter

Type:: float or callable

detuning¶

Optical detuning in GHz

Type:: float

spin_conserving¶

Transition selection rule

Type:: bool

Examples

>>> from sim import HamiltonianBuilder
>>> from sim.hamiltonian.terms import ZFS, OpticalCoupling
>>>
>>> # Constant optical drive
>>> H = HamiltonianBuilder()
>>> H.add(ZFS(D=2.87))
>>> H.add(OpticalCoupling(omega=50))  # 50 MHz Rabi frequency
>>>
>>> # Pulsed optical excitation
>>> def laser_pulse(t):
...     if 0 <= t < 10e-9:  # 10 ns pulse
...         return 100.0  # 100 MHz
...     return 0.0
>>>
>>> H.add(OpticalCoupling(omega=laser_pulse))

Notes

The optical coupling creates coherent superpositions of ground and excited states. In practice, spontaneous emission quickly destroys these superpositions (lifetime ~12 ns).

For optical pumping (initialization), the Lindblad dissipation terms are more important than the coherent coupling.

__init__(omega=0.0, detuning=0.0, spin_conserving=True, name=None)[source]¶

Initialize the Hamiltonian term.

Parameters:

name (str, optional) – Custom name. If not provided, the class name is used.
omega (float | Callable)
detuning (float)
spin_conserving (bool)

property is_time_dependent: bool¶: Optical coupling is always time-dependent.

build(t=0.0)[source]¶

Build the 18×18 optical coupling Hamiltonian at time t.

Parameters:: t (float) – Time in seconds
Returns:: 18×18 complex matrix in rad/s
Return type:: np.ndarray

Notes

The optical coupling connects the ground state manifold (indices 0-8) to the excited state manifold (indices 9-17).

For spin-conserving transitions:: |g, ms, mI⟩ ↔ |e, ms, mI⟩
The coupling matrix has off-diagonal blocks:: H[0:9, 9:18] = (Ω/2) · I_9 H[9:18, 0:9] = (Ω/2) · I_9 (conjugate)

rabi_period()[source]¶

Calculate Rabi oscillation period for constant Rabi frequency.

Returns:: Period in seconds, or None if Rabi freq is time-dependent
Return type:: float or None

Stark Effect¶

Stark effect for NV center.

The NV center exhibits both DC and AC Stark effects when exposed to electric fields. The quadratic Stark effect shifts the spin levels through coupling to the orbital states.

Physics¶

The Stark Hamiltonian is:

H_stark = d_⊥ · E_⊥ · (S_x² - S_y²) + d_∥ · E_∥ · S_z²

where: - d_⊥ ≈ 17 Hz/(V/m) is the perpendicular Stark coefficient - d_∥ ≈ 3.5 Hz/(V/m) is the parallel Stark coefficient - E_⊥, E_∥ are the electric field components

The full tensor form is:

H_stark = Σ_{ij} D_{ij} E_j {S_i, S_i’}/2

Typical Values¶

d_⊥ = 17 ± 3 Hz/(V/m)
d_∥ = 3.5 ± 0.8 Hz/(V/m)

For typical lab fields of 10⁶ V/m, this gives shifts of ~10-20 kHz.

References

[1] van Oort & Glasbeek, Chemical Physics Letters 168, 529 (1990) [2] Dolde et al., Phys. Rev. Lett. 112, 097603 (2014)

class sim.hamiltonian.terms.stark.Stark(d_perp=0.017, d_par=0.0035, E_field=(0, 0, 0), E_field_func=None, name=None)[source]¶

Bases: HamiltonianTerm

Stark effect (DC and AC electric field coupling).

Parameters:

d_perp (float) – Perpendicular Stark coefficient in Hz/(V/m). Default: 17e-3
d_par (float) – Parallel Stark coefficient in Hz/(V/m). Default: 3.5e-3
E_field (array_like) – Static electric field vector [Ex, Ey, Ez] in V/m. Default: (0,0,0)
E_field_func (callable, optional) – Time-dependent E(t) function returning [Ex, Ey, Ez]. If provided, makes the term time-dependent.
name (str, optional) – Custom name for the term

d_perp¶

Perpendicular Stark coefficient in Hz/(V/m)

Type:: float

d_par¶

Parallel Stark coefficient in Hz/(V/m)

Type:: float

E_field¶

Static electric field in V/m

Type:: np.ndarray

Examples

>>> from sim import HamiltonianBuilder
>>> from sim.hamiltonian.terms import ZFS, Stark
>>>
>>> # Static electric field along z
>>> H = HamiltonianBuilder()
>>> H.add(ZFS(D=2.87))
>>> H.add(Stark(E_field=[0, 0, 1e6]))  # 1 MV/m along z
>>>
>>> # AC Stark effect (from optical field)
>>> def ac_field(t):
...     return [1e6 * np.sin(2 * np.pi * 1e9 * t), 0, 0]
>>>
>>> H.add(Stark(E_field_func=ac_field))

Notes

The Stark effect is typically small compared to ZFS and Zeeman effects unless very strong electric fields are applied.

The AC Stark shift from optical fields can be significant when the laser is detuned from the optical transition.

__init__(d_perp=0.017, d_par=0.0035, E_field=(0, 0, 0), E_field_func=None, name=None)[source]¶

Initialize the Hamiltonian term.

Parameters:

name (str, optional) – Custom name. If not provided, the class name is used.
d_perp (float)
d_par (float)
E_field (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str])
E_field_func (Callable | None)

property is_time_dependent: bool¶: True if E_field_func is provided.

build(t=0.0)[source]¶

Build the 18×18 Stark Hamiltonian at time t.

Parameters:: t (float) – Time in seconds
Returns:: 18×18 complex Hermitian matrix in rad/s
Return type:: np.ndarray

shift_khz(E_par=0, E_perp=0)[source]¶

Calculate the Stark shift for given field components.

Parameters:

E_par (float) – Parallel electric field in V/m
E_perp (float) – Perpendicular electric field in V/m

Returns:

(shift_par, shift_perp) in kHz

Return type:

tuple

Strain¶

Strain-induced splitting for NV center.

Strain in the diamond lattice modifies the zero-field splitting through spin-strain coupling. This can be caused by static stress (e.g., from implantation damage) or dynamic strain (e.g., acoustic waves).

Physics¶

The strain Hamiltonian is:

H_strain = Σ_{ij} ε_{ij} · {S_i, S_j}/2

where ε is the strain tensor and {A, B} = AB + BA is the anticommutator.

For axially symmetric strain along the NV axis (z), this simplifies to:

H_strain = E_x (S_x² - S_y²) + E_y (S_x S_y + S_y S_x)

where E_x and E_y are the transverse strain components in MHz.

This is similar to the E-term in ZFS but can be time-dependent.

Strain Tensor¶

The full 3×3 strain tensor has 6 independent components:

ε = [[ε_xx, ε_xy, ε_xz],: [ε_xy, ε_yy, ε_yz], [ε_xz, ε_yz, ε_zz]]

Only certain combinations couple to the NV spin: - E₁ ~ ε_xx - ε_yy (transverse) - E₂ ~ 2ε_xy (transverse) - Shifts D ~ ε_zz (axial)

References

[1] Doherty et al., Physics Reports 528, 1-45 (2013) [2] MacQuarrie et al., Phys. Rev. Lett. 111, 227602 (2013)

class sim.hamiltonian.terms.strain.Strain(Ex=0.0, Ey=0.0, epsilon=None, time_func=None, name=None)[source]¶

Bases: HamiltonianTerm

Strain-induced splitting (static or time-dependent).

Parameters:

Ex (float) – Transverse strain component E_x in MHz. Default: 0
Ey (float) – Transverse strain component E_y in MHz. Default: 0
epsilon (array_like, optional) – Full 3×3 strain tensor (dimensionless). Overrides Ex, Ey.
time_func (callable, optional) – Function time_func(t) returning (Ex, Ey) or 3×3 tensor at time t. If provided, makes the term time-dependent.
name (str, optional) – Custom name for the term

Ex¶

Static E_x component in MHz

Type:: float

Ey¶

Static E_y component in MHz

Type:: float

epsilon¶

Full strain tensor if provided

Type:: np.ndarray or None

is_time_dependent¶

True if time_func is provided

Type:: bool

Examples

>>> from sim import HamiltonianBuilder
>>> from sim.hamiltonian.terms import ZFS, Strain
>>>
>>> # Static strain
>>> H = HamiltonianBuilder()
>>> H.add(ZFS(D=2.87))
>>> H.add(Strain(Ex=10, Ey=5))  # 10 MHz, 5 MHz transverse strain
>>>
>>> # Time-dependent strain (e.g., acoustic wave)
>>> def acoustic_strain(t):
...     freq = 1e9  # 1 GHz acoustic wave
...     amplitude = 5  # 5 MHz
...     return (amplitude * np.sin(2 * np.pi * freq * t), 0)
>>>
>>> H.add(Strain(time_func=acoustic_strain))

Notes

The strain coupling is analogous to the E-term in ZFS: - Ex corresponds to the E_x splitting - Ey corresponds to additional symmetry breaking

Typical strain values in NV experiments: - Bulk diamond: < 1 MHz - Near implantation damage: 1-50 MHz - Under mechanical stress: 10-100 MHz

__init__(Ex=0.0, Ey=0.0, epsilon=None, time_func=None, name=None)[source]¶

Initialize the Hamiltonian term.

Parameters:

name (str, optional) – Custom name. If not provided, the class name is used.
Ex (float)
Ey (float)
epsilon (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None)
time_func (Callable | None)

property is_time_dependent: bool¶: True if time_func is provided.

build(t=0.0)[source]¶

Build the 18×18 strain Hamiltonian at time t.

Parameters:: t (float) – Time in seconds
Returns:: 18×18 complex Hermitian matrix in rad/s
Return type:: np.ndarray

Jahn-Teller¶

Jahn-Teller coupling in the NV center excited state.

The excited state (³E) of the NV center is orbitally degenerate and couples to vibrational modes through the dynamic Jahn-Teller effect. This causes mixing between the orbital states and affects the optical properties.

Physics¶

The Jahn-Teller Hamiltonian describes electron-phonon coupling:

H_JT = ω(a†a + 1/2) + λ(a + a†) · O_orbital

where: - ω is the phonon mode frequency - a, a† are phonon creation/annihilation operators - λ is the electron-phonon coupling strength - O_orbital is an orbital operator

For the E-symmetry modes coupling to the ³E state:

H_JT = ωx(ax†ax + 1/2) + ωy(ay†ay + 1/2)

λx(ax + ax†)(Sx² - Sy²)

λy(ay + ay†){Sx, Sy}

Typical Values¶

Phonon frequency: ω ~ 65 meV (~ 16 THz)
Coupling strength: λ ~ 0.5 (dimensionless)
Huang-Rhys factor: S ~ λ²/2

Simplified Model¶

For the 18D Hilbert space, we use an effective model that captures the main effects without explicit phonon states. This is valid when: - Temperature is low (kT << ℏω) - Phonon dynamics are fast compared to spin dynamics

References

[1] Fu et al., Phys. Rev. Lett. 103, 256404 (2009) [2] Abtew et al., Phys. Rev. Lett. 107, 146403 (2011)

class sim.hamiltonian.terms.jahn_teller.JahnTeller(omega_x=1000.0, omega_y=1200.0, lambda_x=0.05, lambda_y=0.04, temperature=300.0, name=None)[source]¶

Bases: HamiltonianTerm

Jahn-Teller vibronic coupling in excited state (simplified model).

This is a simplified effective model that captures the main effects of Jahn-Teller coupling without explicit phonon states.

Parameters:

omega_x (float) – Phonon mode frequency (x-mode) in MHz. Default: 1000 (1 GHz)
omega_y (float) – Phonon mode frequency (y-mode) in MHz. Default: 1200 (1.2 GHz)
lambda_x (float) – Coupling strength (x-mode), dimensionless. Default: 0.05
lambda_y (float) – Coupling strength (y-mode), dimensionless. Default: 0.04
temperature (float) – Temperature in Kelvin. Default: 300 K
name (str, optional) – Custom name for the term

omega_x, omega_y

Phonon frequencies in MHz

Type:: float

lambda_x, lambda_y

Coupling strengths (dimensionless)

Type:: float

temperature¶

Temperature in Kelvin

Type:: float

Examples

>>> from sim import HamiltonianBuilder
>>> from sim.hamiltonian.terms import ZFS, JahnTeller
>>>
>>> H = HamiltonianBuilder()
>>> H.add(ZFS(D=2.87))
>>> H.add(JahnTeller())  # Default parameters
>>>
>>> # Low temperature (cryogenic)
>>> H.add(JahnTeller(temperature=4.0))

Notes

The Jahn-Teller effect is primarily important for: - Understanding the excited state fine structure - Calculating optical transition rates - Modeling temperature dependence

At room temperature, the JT effect averages out due to fast phonon dynamics. At low temperature, it can cause strain-like splitting in the excited state.

This simplified model only acts on the excited state manifold.

__init__(omega_x=1000.0, omega_y=1200.0, lambda_x=0.05, lambda_y=0.04, temperature=300.0, name=None)[source]¶

Initialize the Hamiltonian term.

Parameters:

name (str, optional) – Custom name. If not provided, the class name is used.
omega_x (float)
omega_y (float)
lambda_x (float)
lambda_y (float)
temperature (float)

property is_time_dependent: bool¶: Jahn-Teller coupling is time-independent (effective model).

build(t=0.0)[source]¶

Build the 18×18 Jahn-Teller Hamiltonian.

Parameters:: t (float) – Time in seconds (not used)
Returns:: 18×18 complex Hermitian matrix in rad/s
Return type:: np.ndarray

Notes

The JT effect only acts on the excited state manifold (indices 9-17). The ground state is not affected.

huang_rhys_factor()[source]¶

Calculate the Huang-Rhys factors for both modes.

Returns:: (S_x, S_y) Huang-Rhys factors (dimensionless)
Return type:: tuple

Notes

The Huang-Rhys factor S = λ²/2 characterizes the strength of electron-phonon coupling. S ~ 0.5 corresponds to moderate coupling where the Stokes shift equals the phonon energy.

Creating Custom Terms¶

To create a custom Hamiltonian term:

from sim.hamiltonian.terms.base import HamiltonianTerm
import numpy as np

class CustomTerm(HamiltonianTerm):
    def __init__(self, parameter, name=None):
        super().__init__(name=name)
        self.parameter = parameter

    @property
    def is_time_dependent(self):
        return False

    def build(self, t=0.0):
        # Return 18×18 Hermitian matrix
        H = np.zeros((18, 18), dtype=np.complex128)
        # ... fill matrix ...
        return H