Dynamics Module¶

The dynamics module provides solvers for open quantum system evolution using the Lindblad master equation.

Lindblad Solver¶

Lindblad master equation solver for NV center dynamics.

Solves the Lindblad master equation for open quantum systems:

dρ/dt = -i[H(t), ρ] + Σ_k D[L_k](ρ)

where the Lindblad dissipator is:

D[L](ρ) = L ρ L† - ½{L†L, ρ}

The solver uses scipy.integrate.solve_ivp with adaptive step size control (RK45 method by default).

Vectorization¶

For numerical integration, the density matrix is vectorized:

ρ (18×18) → ρ_vec (324)

The Lindblad equation becomes:

dρ_vec/dt = L ρ_vec

where L is the Liouvillian superoperator (324×324).

For time-dependent Hamiltonians, L is recomputed at each time step.

class sim.dynamics.lindblad.EvolutionResult(times, rho_t, success=True)[source]¶

Bases: object

Container for time evolution results.

Parameters:

times (ndarray)
rho_t (List[ndarray])
success (bool)

times¶

Array of time points in seconds

Type:: np.ndarray

rho_t¶

List of density matrices at each time point

Type:: list of np.ndarray

success¶

Whether the integration succeeded

Type:: bool

expectation(operator)[source]¶

Compute ⟨O⟩(t) for all times

Parameters:: operator (ndarray)
Return type:: ndarray

population(projector)[source]¶

Compute population for all times

Parameters:: projector (ndarray)
Return type:: ndarray

final_state()[source]¶

Return the final density matrix

Return type:: ndarray

times: ndarray¶

rho_t: List[ndarray]¶

success: bool = True¶

expectation(operator)[source]¶

Compute expectation value ⟨O⟩(t) for all times.

Parameters:: operator (np.ndarray) – 18×18 operator matrix
Returns:: Array of expectation values (complex)
Return type:: np.ndarray

population(projector)[source]¶

Compute population Tr(P ρ) for all times.

Parameters:: projector (np.ndarray) – 18×18 projector matrix
Returns:: Array of populations (real, in [0,1])
Return type:: np.ndarray

final_state()[source]¶

Return the final density matrix.

Return type:: ndarray

purity()[source]¶

Compute purity Tr(ρ²) for all times.

Return type:: ndarray

__init__(times, rho_t, success=True)¶

Parameters:

times (ndarray)
rho_t (List[ndarray])
success (bool)

Return type:

None

class sim.dynamics.lindblad.LindbladSolver(hamiltonian)[source]¶

Bases: object

Lindblad master equation solver.

Parameters:: hamiltonian (HamiltonianBuilder) – The system Hamiltonian (can be time-dependent)

hamiltonian¶

The Hamiltonian builder

Type:: HamiltonianBuilder

dissipators¶

List of Lindblad operators

Type:: list

Examples

>>> from sim import HamiltonianBuilder
>>> from sim.hamiltonian.terms import ZFS, Zeeman, MicrowaveDrive
>>> from sim.dynamics import LindbladSolver
>>> from sim.states import ground_state
>>>
>>> # Build time-dependent Hamiltonian
>>> H = HamiltonianBuilder()
>>> H.add(ZFS(D=2.87))
>>> H.add(Zeeman(B=10))
>>> H.add(MicrowaveDrive(omega=10))  # 10 MHz Rabi
>>>
>>> # Setup solver with dissipation
>>> solver = LindbladSolver(H)
>>> solver.add_t1_relaxation(gamma=1e3)    # T1 = 1 ms
>>> solver.add_t2_dephasing(gamma=1e6)     # T2* = 1 μs
>>>
>>> # Evolve from ground state
>>> rho0 = ground_state()
>>> result = solver.evolve(rho0, t_span=(0, 1e-6), n_steps=100)
>>>
>>> # Analyze
>>> from sim.states import projector_ms
>>> pop_ms0 = result.population(projector_ms(0))
>>> print(f"Final ms=0 population: {pop_ms0[-1]:.3f}")

DIM = 18¶

__init__(hamiltonian)[source]¶

Parameters:: hamiltonian (HamiltonianBuilder)

property dissipators: List[ndarray]¶: List of Lindblad operators.

add_dissipator(L)[source]¶

Add a Lindblad operator.

Parameters:: L (np.ndarray) – 18×18 Lindblad operator (already includes √γ)

add_t1_relaxation(gamma, manifold='ground')[source]¶

Add T1 spin-lattice relaxation.

Parameters:

gamma (float) – Relaxation rate 1/T1 in Hz
manifold (str) – “ground”, “excited”, or “both”

add_t2_dephasing(gamma, manifold='ground')[source]¶

Add T2* pure dephasing.

Parameters:

gamma (float) – Dephasing rate 1/T2* in Hz
manifold (str) – “ground”, “excited”, or “both”

add_optical_decay(gamma)[source]¶

Add optical spontaneous emission.

Parameters:: gamma (float) – Emission rate in Hz (~8×10⁷ for NV)

clear_dissipators()[source]¶: Remove all dissipation channels.

evolve(rho0, t_span, n_steps=100, method='RK45', rtol=1e-08, atol=1e-10)[source]¶

Evolve the density matrix in time.

Parameters:

rho0 (np.ndarray) – Initial density matrix (18×18)
t_span (tuple) – (t_start, t_end) in seconds
n_steps (int) – Number of output time points. Default: 100
method (str) – ODE solver method. Default: “RK45”
rtol (float) – Relative tolerance. Default: 1e-8
atol (float) – Absolute tolerance. Default: 1e-10

Returns:

Container with times and density matrices

Return type:

EvolutionResult

steady_state(t_max=0.001, tol=1e-10)[source]¶

Find the steady state by long-time evolution.

Parameters:

t_max (float) – Maximum evolution time in seconds
tol (float) – Convergence tolerance

Returns:

Steady state density matrix, or None if not converged

Return type:

np.ndarray or None

Basic Usage¶

from sim import HamiltonianBuilder, LindbladSolver
from sim.hamiltonian.terms import ZFS, Zeeman
from sim.states import ground_state

# Build Hamiltonian
H = HamiltonianBuilder()
H.add(ZFS(D=2.87))
H.add(Zeeman(B=10))

# Create solver
solver = LindbladSolver(H)

# Add dissipation channels
solver.add_t1_relaxation(gamma=1e3)    # T1 = 1 ms
solver.add_t2_dephasing(gamma=1e6)     # T2* = 1 μs
solver.add_optical_decay(gamma=8e7)    # τ = 12 ns

# Initial state
rho0 = ground_state()

# Evolve
result = solver.evolve(
    rho0,
    t_span=(0, 1e-6),
    n_steps=100,
    method="RK45"
)

Evolution Result¶

The evolve() method returns an EvolutionResult object:

# Access time points
times = result.times  # Array of time values

# Access density matrices
rho_t = result.rho_t  # List of 18×18 matrices

# Get final state
rho_final = result.final_state()

# Compute expectation values
from sim.core.operators import Sz
exp_Sz = result.expectation(Sz)

# Compute populations
from sim.states import projector_ms
pop_ms0 = result.population(projector_ms(0))

# Compute purity
purity = result.purity()

Steady State¶

Find the asymptotic steady state:

rho_ss = solver.steady_state(t_max=1e-3, tol=1e-10)

Dissipation Operators¶

Dissipation operators for Lindblad master equation.

Lindblad operators describe irreversible processes like: - T1 relaxation (spin-lattice relaxation) - T2* dephasing (inhomogeneous broadening) - T2 pure dephasing (homogeneous) - Optical decay (spontaneous emission)

Lindblad Form¶

Each dissipation channel is described by a Lindblad operator L_k that appears in the master equation as:

D[ρ] = L ρ L† - ½{L†L, ρ}

The rate is incorporated into the operator: L = √γ × (bare operator)

Time Scales for NV Centers¶

T1 (spin): 1-10 ms at room temperature
T2* (inhomogeneous): 1-10 μs
T2 (Hahn echo): 100-500 μs
Optical T1: ~12 ns (excited state lifetime)

sim.dynamics.dissipation.t1_operators(gamma, manifold='ground')[source]¶

Create T1 relaxation (spin-lattice) Lindblad operators.

Models incoherent transitions between spin states, driving the system toward thermal equilibrium.

Parameters:

gamma (float) – Relaxation rate 1/T1 in Hz
manifold (str) – Which manifold: “ground”, “excited”, or “both”

Returns:

List of 18×18 Lindblad operators [L_+, L_-]

Return type:

list of np.ndarray

Notes

For spin-1 systems, we use:: L_+ = √γ · S_+ (excitation) L_- = √γ · S_- (relaxation)

At T=0, only L_- contributes (pure decay). At finite T, both contribute with rates determined by detailed balance.

Examples

>>> ops = t1_operators(gamma=1e3)  # T1 = 1 ms
>>> len(ops)
2

sim.dynamics.dissipation.t2_dephasing_operator(gamma, manifold='ground')[source]¶

Create T2* dephasing Lindblad operator.

Models pure dephasing that destroys coherences without changing populations.

Parameters:

gamma (float) – Dephasing rate 1/T2* in Hz
manifold (str) – Which manifold: “ground”, “excited”, or “both”

Returns:

18×18 Lindblad operator

Return type:

np.ndarray

Notes

The dephasing operator is L = √γ · Sz.

This causes decay of off-diagonal elements (coherences) at rate γ, while diagonal elements (populations) are unchanged.

The total dephasing rate is:: 1/T2* = 1/(2T1) + 1/T_φ

where T_φ is the pure dephasing time.

Examples

>>> L = t2_dephasing_operator(gamma=1e7)  # T2* = 100 ns
>>> L.shape
(18, 18)

sim.dynamics.dissipation.optical_decay_operators(gamma)[source]¶

Create optical decay (spontaneous emission) operators.

Models transitions from excited to ground state with photon emission. Spin-conserving transitions.

Parameters:: gamma (float) – Spontaneous emission rate in Hz (~1/12 ns ≈ 8×10⁷ Hz)
Returns:: List of 9 Lindblad operators (one per |e,ms,mI⟩ → |g,ms,mI⟩)
Return type:: list of np.ndarray

Notes

Each spin state decays independently:: L_{ms,mI} = √γ · |g,ms,mI⟩⟨e,ms,mI|
The total emission rate from the excited state is:: Γ = 9γ (summed over all channels)

For NV: τ_rad ≈ 12 ns → γ ≈ 8×10⁷ Hz

Examples

>>> ops = optical_decay_operators(gamma=8e7)
>>> len(ops)
9

sim.dynamics.dissipation.nuclear_relaxation_operators(gamma, manifold='ground')[source]¶

Create nuclear spin T1 relaxation operators.

Models incoherent transitions of the N14 nuclear spin.

Parameters:

gamma (float) – Nuclear relaxation rate in Hz (typically ~100 Hz)
manifold (str) – Which manifold: “ground”, “excited”, or “both”

Returns:

List of 18×18 Lindblad operators [L_+, L_-]

Return type:

list of np.ndarray

Notes

Nuclear T1 is typically very long (~seconds) but can be accelerated by certain processes.

sim.dynamics.dissipation.intersystem_crossing_operators(gamma_isc, gamma_singlet)[source]¶

Create intersystem crossing (ISC) operators.

Models the non-radiative decay pathway through singlet states that enables optical spin polarization.

Parameters:

gamma_isc (float) – ISC rate from excited state (spin-dependent)
gamma_singlet (float) – Decay rate from singlet back to ground state

Returns:

Lindblad operators for ISC process

Return type:

list of np.ndarray

Notes

This is a simplified model. The full ISC involves: - Decay from ³E (ms=±1) to singlet states - Decay from singlet to ³A₂ (preferentially ms=0)

This is what enables optical spin polarization.

T₁ Relaxation¶

Spin-lattice relaxation toward thermal equilibrium:

solver.add_t1_relaxation(gamma=1e3, manifold="ground")

Parameters:

gamma: Relaxation rate 1/T₁ in Hz
manifold: “ground”, “excited”, or “both”

T₂* Dephasing¶

Pure dephasing destroying coherences:

solver.add_t2_dephasing(gamma=1e6, manifold="ground")

Optical Decay¶

Spontaneous emission from excited state:

solver.add_optical_decay(gamma=8e7)  # τ ≈ 12 ns

Custom Dissipators¶

Add any Lindblad operator:

import numpy as np

# Create 18×18 Lindblad operator
L = np.zeros((18, 18), dtype=np.complex128)
L[0, 9] = np.sqrt(gamma)  # |g,+1,+1⟩ ← |e,+1,+1⟩

solver.add_dissipator(L)

Solver Options¶

Integration Methods¶

Available methods (via scipy.integrate.solve_ivp):

"RK45": Default, 4th-order Runge-Kutta (adaptive)
"RK23": 2nd-order Runge-Kutta
"DOP853": 8th-order Dormand-Prince
"Radau": Implicit Runge-Kutta (for stiff problems)
"BDF": Backward differentiation (for stiff problems)

Tolerances¶

result = solver.evolve(
    rho0,
    t_span=(0, 1e-6),
    rtol=1e-8,   # Relative tolerance
    atol=1e-10   # Absolute tolerance
)

For faster but less accurate simulations, increase tolerances.

Time-Dependent Hamiltonians¶

The solver handles time-dependent Hamiltonians automatically:

from sim.hamiltonian.terms import MicrowaveDrive

# Pulsed microwave
def pulse(t):
    if 0 < t < 50e-9:  # 50 ns π-pulse
        return 10  # 10 MHz Rabi
    return 0

H = HamiltonianBuilder()
H.add(ZFS(D=2.87))
H.add(MicrowaveDrive(omega=pulse))

solver = LindbladSolver(H)
# The Hamiltonian is rebuilt at each integration step