States Module

The states module provides functions for creating initial states, density matrices, and measurement projectors.

Initial States

Initial state preparation for NV center simulation.

Provides functions to create common initial states for NV center experiments. The standard initialization is the |g, ms=0, mI=0⟩ state achieved through optical pumping.

Basis Indexing

The 18-dimensional basis follows the convention:

|idx⟩ = |ge⟩ ⊗ |ms⟩ ⊗ |mI⟩

Index = 9 * ge + 3 * ms_idx + mI_idx

where:

ge = 0 (ground) or 1 (excited) ms_idx = 0 (+1), 1 (0), 2 (-1) mI_idx = 0 (+1), 1 (0), 2 (-1)

Common indices:

|g, +1, +1⟩ = 0 |g, +1, 0⟩ = 1 |g, 0, 0⟩ = 4 ← Standard initialization |g, -1, 0⟩ = 7 |e, 0, 0⟩ = 13

sim.states.initial_states.state_index(ge, ms, mI)[source]

Convert quantum numbers to state index.

Parameters:
  • ge (int) – Electronic state: 0 (ground) or 1 (excited)

  • ms (int) – Electron spin: +1, 0, or -1

  • mI (int) – Nuclear spin: +1, 0, or -1

Returns:

Index in 0-17 range

Return type:

int

sim.states.initial_states.ground_state_vector(ms=0, mI=0)[source]

Create pure state vector |g, ms, mI⟩.

Parameters:
  • ms (int) – Electron spin projection: +1, 0, or -1. Default: 0

  • mI (int) – Nuclear spin projection: +1, 0, or -1. Default: 0

Returns:

18-dimensional normalized state vector

Return type:

np.ndarray

Examples

>>> psi = ground_state_vector()  # |g, 0, 0⟩
>>> print(f"Norm: {np.linalg.norm(psi):.1f}")
Norm: 1.0
>>> psi_m1 = ground_state_vector(ms=-1)  # |g, -1, 0⟩
sim.states.initial_states.ground_state(ms=0, mI=0)[source]

Create density matrix for |g, ms, mI⟩⟨g, ms, mI|.

Parameters:
  • ms (int) – Electron spin: +1, 0, or -1. Default: 0

  • mI (int) – Nuclear spin: +1, 0, or -1. Default: 0

Returns:

18×18 density matrix

Return type:

np.ndarray

Examples

>>> rho = ground_state()  # Standard NV initialization
>>> print(f"Trace: {np.trace(rho):.1f}")
Trace: 1.0
sim.states.initial_states.excited_state(ms=0, mI=0)[source]

Create density matrix for |e, ms, mI⟩⟨e, ms, mI|.

Parameters:
  • ms (int) – Electron spin: +1, 0, or -1. Default: 0

  • mI (int) – Nuclear spin: +1, 0, or -1. Default: 0

Returns:

18×18 density matrix

Return type:

np.ndarray

sim.states.initial_states.superposition(coefficients, normalize=True)[source]

Create superposition state from coefficients.

Parameters:
  • coefficients (dict) – Mapping (ge, ms, mI) -> coefficient

  • normalize (bool) – Whether to normalize the state. Default: True

Returns:

18×18 density matrix

Return type:

np.ndarray

Examples

>>> # |ψ⟩ = (|0⟩ + |−1⟩)/√2
>>> coeffs = {
...     (0, 0, 0): 1.0,
...     (0, -1, 0): 1.0
... }
>>> rho = superposition(coeffs)
>>> print(f"Purity: {np.real(np.trace(rho @ rho)):.3f}")
Purity: 1.000
>>> # Mixed nuclear spin
>>> coeffs = {
...     (0, 0, +1): 1/np.sqrt(3),
...     (0, 0, 0): 1/np.sqrt(3),
...     (0, 0, -1): 1/np.sqrt(3)
... }
>>> rho = superposition(coeffs)
sim.states.initial_states.thermal_state(H, temperature)[source]

Create thermal equilibrium state ρ = exp(-βH) / Z.

Parameters:
  • H (np.ndarray) – 18×18 Hamiltonian matrix in rad/s

  • temperature (float) – Temperature in Kelvin

Returns:

18×18 density matrix

Return type:

np.ndarray

Examples

>>> from sim import HamiltonianBuilder
>>> from sim.hamiltonian.terms import ZFS
>>>
>>> H = HamiltonianBuilder().add(ZFS(D=2.87)).build()
>>> rho_300K = thermal_state(H, 300)  # Room temperature
>>> rho_4K = thermal_state(H, 4)      # Cryogenic
sim.states.initial_states.maximally_mixed()[source]

Create maximally mixed state ρ = I/18.

Returns:

18×18 density matrix with equal populations

Return type:

np.ndarray

Notes

This state has maximum entropy S = log(18) and minimum purity Tr(ρ²) = 1/18.

sim.states.initial_states.ms0_superposition(phase=0.0, mI=0)[source]

Create superposition of ms=+1 and ms=-1 states.

|ψ⟩ = (|+1⟩ + e^(iφ)|−1⟩)/√2

Parameters:
  • phase (float) – Relative phase in radians. Default: 0

  • mI (int) – Nuclear spin state. Default: 0

Returns:

18×18 density matrix

Return type:

np.ndarray

sim.states.initial_states.bell_state_electron_nuclear(state_type='phi_plus', ge=0)[source]

Create Bell state between electron and nuclear spin.

Parameters:
  • state_type (str) – Type of Bell state: - “phi_plus”: (|+1,+1⟩ + |−1,−1⟩)/√2 - “phi_minus”: (|+1,+1⟩ − |−1,−1⟩)/√2 - “psi_plus”: (|+1,−1⟩ + |−1,+1⟩)/√2 - “psi_minus”: (|+1,−1⟩ − |−1,+1⟩)/√2

  • ge (int) – Electronic manifold: 0 (ground) or 1 (excited)

Returns:

18×18 density matrix

Return type:

np.ndarray

Common States

from sim.states import (
    ground_state,
    excited_state,
    thermal_state,
    superposition
)

# Ground state |g, ms=0, mI=0⟩
rho0 = ground_state()

# Excited state |e, ms=0, mI=0⟩
rho_e = excited_state()

# Thermal state at temperature T
rho_th = thermal_state(T=300)  # Kelvin

# Superposition state
rho_sup = superposition(
    states=[(0, 0, 0), (0, 1, 0)],  # (g/e, ms, mI)
    weights=[1, 1]
)

Spin Eigenstates

Create states with specific spin projections:

from sim.states import spin_state

# |g, ms=+1, mI=0⟩
rho = spin_state(ms=+1, mI=0, manifold="ground")

# |e, ms=-1, mI=+1⟩
rho = spin_state(ms=-1, mI=+1, manifold="excited")

Density Matrix

Density matrix class for NV center quantum states.

The density matrix ρ describes the quantum state of the system, including mixed states (statistical mixtures) and pure states.

Properties

For a valid density matrix: - Hermitian: ρ = ρ† - Positive semi-definite: all eigenvalues ≥ 0 - Trace one: Tr(ρ) = 1

For a pure state |ψ⟩: ρ = |ψ⟩⟨ψ| with Tr(ρ²) = 1

Basis Convention

The 18-dimensional basis is ordered as:

|g/e⟩ ⊗ |ms⟩ ⊗ |mI⟩

Index mapping:

idx = 9 * ge + 3 * (1 - ms) + (1 - mI)

where ge = 0 (ground) or 1 (excited), ms ∈ {+1, 0, -1}, mI ∈ {+1, 0, -1}.

class sim.states.density_matrix.DensityMatrix(rho=None)[source]

Bases: object

Density matrix for the 18-dimensional NV Hilbert space.

Parameters:

rho (np.ndarray, optional) – 18×18 density matrix. If None, creates zero matrix.

rho

The 18×18 density matrix

Type:

np.ndarray

DIM

Hilbert space dimension (18)

Type:

int

Examples

>>> from sim.states import DensityMatrix, ground_state_vector
>>>
>>> # Create from pure state
>>> psi = ground_state_vector()
>>> dm = DensityMatrix.from_pure_state(psi)
>>> print(f"Purity: {dm.purity():.3f}")
Purity: 1.000
>>>
>>> # Check properties
>>> print(f"Trace: {dm.trace():.3f}")
Trace: 1.000
>>> print(f"Is valid: {dm.is_valid()}")
Is valid: True

Notes

The density matrix is stored as a complex numpy array. All operations return copies to prevent accidental modification.

DIM = 18
__init__(rho=None)[source]

Initialize density matrix.

Parameters:

rho (ndarray | None)

property rho: ndarray

Return copy of the density matrix.

classmethod from_pure_state(psi)[source]

Create density matrix from pure state vector.

Parameters:

psi (np.ndarray) – 18-dimensional state vector

Returns:

ρ = |ψ⟩⟨ψ|

Return type:

DensityMatrix

Examples

>>> psi = np.zeros(18, dtype=complex)
>>> psi[4] = 1.0  # |g, ms=0, mI=0⟩
>>> dm = DensityMatrix.from_pure_state(psi)
>>> dm.purity()
1.0
classmethod from_thermal(H, temperature, zero_point=None)[source]

Create thermal equilibrium state.

Parameters:
  • H (np.ndarray) – 18×18 Hamiltonian matrix in rad/s

  • temperature (float) – Temperature in Kelvin

  • zero_point (float, optional) – Energy zero point. If None, uses ground state energy.

Returns:

ρ = exp(-H/kT) / Z

Return type:

DensityMatrix

Notes

At T → 0, returns the ground state. At T → ∞, returns the maximally mixed state.

trace()[source]

Return trace of density matrix.

Return type:

float

purity()[source]

Return purity Tr(ρ²).

Returns:

Purity in [1/18, 1]. Pure state has purity 1.

Return type:

float

entropy()[source]

Return von Neumann entropy S = -Tr(ρ log ρ).

Returns:

Entropy in natural units (nats). Zero for pure states.

Return type:

float

expectation(operator)[source]

Compute expectation value ⟨O⟩ = Tr(O ρ).

Parameters:

operator (np.ndarray) – 18×18 operator matrix

Returns:

Expectation value

Return type:

complex

population(projector)[source]

Compute population Tr(P ρ) for a projector.

Parameters:

projector (np.ndarray) – 18×18 projector matrix (P² = P)

Returns:

Population in [0, 1]

Return type:

float

is_valid(tol=1e-10)[source]

Check if this is a valid density matrix.

Parameters:

tol (float) – Tolerance for numerical checks

Returns:

True if Hermitian, positive semi-definite, and trace 1

Return type:

bool

fidelity(other)[source]

Compute fidelity F(ρ, σ) between two states.

Parameters:

other (DensityMatrix) – Another density matrix

Returns:

Fidelity in [0, 1]

Return type:

float

Notes

F(ρ, σ) = [Tr(√(√ρ σ √ρ))]²

For pure states: F = |⟨ψ|φ⟩|²

partial_trace(subsystem)[source]

Trace out a subsystem.

Parameters:

subsystem (str) – Which subsystem to trace out: - “electronic”: Trace out |g/e⟩, returns 9×9 - “electron_spin”: Trace out |ms⟩, returns 6×6 - “nuclear_spin”: Trace out |mI⟩, returns 6×6

Returns:

Reduced density matrix

Return type:

np.ndarray

The DensityMatrix class provides utilities for working with quantum states:

from sim.states import DensityMatrix

# Create from array
dm = DensityMatrix(rho_array)

# Properties
print(dm.trace())       # Should be 1
print(dm.purity())      # Tr(ρ²)
print(dm.is_valid())    # Check positivity

# Partial traces
rho_spin = dm.partial_trace_nuclear()    # 6×6
rho_nuclear = dm.partial_trace_spin()    # 6×6
rho_electron = dm.partial_trace_manifold()  # 9×9

Projectors

Projectors and measurement operators for NV center simulation.

Projectors are used to calculate populations and implement measurements. A projector P satisfies P² = P and P† = P.

Population Calculation

The population in a subspace is:

p = Tr(P ρ)

where P is the projector onto that subspace.

Basis Structure

The 18-dimensional basis:

|idx⟩ = |ge⟩ ⊗ |ms⟩ ⊗ |mI⟩

Ground state manifold: indices 0-8 Excited state manifold: indices 9-17

sim.states.projectors.projector_ground()[source]

Projector onto ground state manifold P_g.

Returns:

18×18 projector (indices 0-8)

Return type:

np.ndarray

Examples

>>> P_g = projector_ground()
>>> rho = ground_state()
>>> pop = np.real(np.trace(P_g @ rho))
>>> print(f"Ground population: {pop:.1f}")
Ground population: 1.0
sim.states.projectors.projector_excited()[source]

Projector onto excited state manifold P_e.

Returns:

18×18 projector (indices 9-17)

Return type:

np.ndarray

sim.states.projectors.projector_ms(ms, manifold='both')[source]

Projector onto electron spin state ms.

Parameters:
  • ms (int) – Electron spin projection: +1, 0, or -1

  • manifold (str) – Which electronic manifold: - “ground”: Only ground state - “excited”: Only excited state - “both”: Both manifolds (default)

Returns:

18×18 projector

Return type:

np.ndarray

Examples

>>> P_ms0 = projector_ms(0)  # All ms=0 states
>>> P_ms0_g = projector_ms(0, manifold="ground")  # Only ground
sim.states.projectors.projector_mI(mI, manifold='both')[source]

Projector onto nuclear spin state mI.

Parameters:
  • mI (int) – Nuclear spin projection: +1, 0, or -1

  • manifold (str) – Which electronic manifold: “ground”, “excited”, or “both”

Returns:

18×18 projector

Return type:

np.ndarray

sim.states.projectors.projector_state(ge, ms, mI)[source]

Projector onto single basis state |ge, ms, mI⟩.

Parameters:
  • ge (int) – Electronic state: 0 (ground) or 1 (excited)

  • ms (int) – Electron spin: +1, 0, or -1

  • mI (int) – Nuclear spin: +1, 0, or -1

Returns:

18×18 rank-1 projector

Return type:

np.ndarray

sim.states.projectors.spin_operator_18x18(component, manifold='both')[source]

Electron spin operator in 18×18 space.

Parameters:
  • component (str) – Spin component: “x”, “y”, “z”, “+”, or “-”

  • manifold (str) – Which manifold: “ground”, “excited”, or “both”

Returns:

18×18 spin operator

Return type:

np.ndarray

Examples

>>> Sz = spin_operator_18x18("z")
>>> rho = ground_state(ms=1)
>>> sz = np.real(np.trace(Sz @ rho))
>>> print(f"<Sz> = {sz:.1f}")
<Sz> = 1.0
sim.states.projectors.nuclear_operator_18x18(component, manifold='both')[source]

Nuclear spin operator (N14) in 18×18 space.

Parameters:
  • component (str) – Spin component: “x”, “y”, “z”, “+”, or “-”

  • manifold (str) – Which manifold: “ground”, “excited”, or “both”

Returns:

18×18 nuclear spin operator

Return type:

np.ndarray

sim.states.projectors.population_dict(rho)[source]

Extract all populations from a density matrix.

Parameters:

rho (np.ndarray) – 18×18 density matrix

Returns:

Dictionary with keys like “g_ms0_mI0” and population values

Return type:

dict

Examples

>>> rho = ground_state()
>>> pops = population_dict(rho)
>>> print(f"Ground ms=0: {pops['g_ms0_mI0']:.2f}")
Ground ms=0: 1.00

Spin Projectors

from sim.states import projector_ms, projector_mI

# Project onto ms=0 (ground state)
P_ms0 = projector_ms(0, manifold="ground")

# Project onto ms=±1
P_ms1 = projector_ms(+1, manifold="ground")
P_msm1 = projector_ms(-1, manifold="ground")

# Nuclear spin projectors
P_mI0 = projector_mI(0)

# Measure population
pop = np.trace(P_ms0 @ rho).real

Electronic Manifold Projectors

from sim.states import projector_ground, projector_excited

# Ground state manifold (indices 0-8)
P_g = projector_ground()

# Excited state manifold (indices 9-17)
P_e = projector_excited()

# Excited state population
pop_excited = np.trace(P_e @ rho).real

Combined Projectors

from sim.states import projector_state

# Project onto specific state |g, ms=0, mI=0⟩
P = projector_state(ge="ground", ms=0, mI=0)

# Population in that state
pop = np.trace(P @ rho).real

State Tomography

Reconstruct the density matrix from measurements:

from sim.states import tomography_operators

# Get set of measurement operators
operators = tomography_operators()

# Each operator's expectation gives one data point
expectations = [np.trace(O @ rho).real for O in operators]

# Reconstruct (maximum likelihood)
rho_reconstructed = reconstruct_state(expectations)

Visualization

Tools for visualizing quantum states:

from sim.states import plot_populations, plot_bloch

# Bar plot of spin populations
plot_populations(rho)

# Bloch sphere representation
plot_bloch(rho)

# Density matrix heatmap
import matplotlib.pyplot as plt
plt.imshow(np.abs(rho))
plt.colorbar()