Lindblad Master Equation Solver

July 13, 2026 · View on GitHub

scpn_quantum_control.phase.lindblad

Open-system dynamics for the Kuramoto-XY Hamiltonian via the Lindblad master equation. Solves for the full density matrix under amplitude damping and dephasing channels.

Caveat: Full density matrix evolution scales as $O(4^n)$ in memory. For $n > 10$ , consider the MCWF method or MPS/DMRG instead.

A closed quantum system evolves unitarily: $d|\psi\rangle/dt = -iH|\psi\rangle$ . Real systems interact with their environment. The Lindblad equation is the most general Markovian master equation that preserves trace, Hermiticity, and positivity of the density matrix $\rho$ :

$\frac{d\rho}{dt} = -i[H, \rho] + \sum_k \left( L_k \rho L_k^\dagger - \frac{1}{2}\{L_k^\dagger L_k, \rho\} \right)$

The first term generates coherent (unitary) evolution. The second term — the dissipator — describes irreversible coupling to the environment through Lindblad operators $L_k$ .

Channels in This Module

Two physical channels are implemented, parameterised per qubit:

Channel	Lindblad operator	Rate	Physical meaning
Amplitude damping	$L_k = \sqrt{\gamma_\text{amp}}\, \sigma^-_k$	$\gamma_\text{amp}$	Energy relaxation ( $T_1$ decay)
Pure dephasing	$L_k = \sqrt{\gamma_\text{deph}/2}\, \sigma^z_k$	$\gamma_\text{deph}$	Phase randomisation ( $T_2$ decay)

For the Kuramoto-XY system, amplitude damping destroys synchronisation by relaxing excitations toward the ground state. Dephasing destroys off-diagonal coherences without changing populations.

The XY Hamiltonian

$H = -\sum_{i<j} K_{ij}(X_i X_j + Y_i Y_j) - \sum_i \omega_i Z_i$

where $K_{ij}$ is the coupling matrix (typically exponentially decaying with distance) and $\omega_i$ are natural frequencies.

Order Parameter and Purity

Kuramoto order parameter $R$ : extracted from single-qubit Pauli expectations $\langle X_k \rangle$ , $\langle Y_k \rangle$ via the density matrix. Quantifies synchronisation ( $R=1$ perfect sync, $R \to 0$ incoherent).
Purity $\text{Tr}(\rho^2)$ : $1 $for a pure state, \$ 1/d$ for maximally mixed. Decreases under dissipation.

API Reference

`LindbladKuramotoSolver`

from scpn_quantum_control.phase.lindblad import LindbladKuramotoSolver

Constructor

LindbladKuramotoSolver(
    n_oscillators: int,
    K_coupling: np.ndarray,       # shape (n, n)
    omega_natural: np.ndarray,    # shape (n,)
    gamma_amp: float = 0.0,       # amplitude damping rate
    gamma_deph: float = 0.0,      # dephasing rate
    *,
    max_dense_gib: float | None = None,
)

Parameters:

Parameter	Type	Description
`n_oscillators`	`int`	Positive number of qubits and oscillators.
`K_coupling`	`ndarray (n, n)`	Finite real symmetric coupling matrix. The diagonal is discarded.
`omega_natural`	`ndarray (n,)`	Finite real natural frequencies ordered like the rows of `K_coupling`.
`gamma_amp`	`float`	Finite non-negative amplitude-damping rate per qubit. $\gamma = 0$ disables damping.
`gamma_deph`	`float`	Finite non-negative pure-dephasing rate per qubit. $\gamma = 0$ disables dephasing.
`max_dense_gib`	`float	None`

If max_dense_gib is omitted, dense allocation uses SCPN_MAX_DENSE_GIB when set, otherwise the shared host-aware default. build() estimates the simultaneous Hamiltonian, density-matrix, work-array, and channel-operator footprint and raises DenseAllocationError before an over-budget allocation.

Methods

Method	Signature	Returns	Description
`build()`	`(*, max_dense_gib=None) → None`	—	Build and cache the Hamiltonian and channel operators under the active dense budget.
`run()`	`(t_max, dt, method="RK45", *, max_dense_gib=None) → dict`	See below	Evolve through `t_max`; `dt` bounds adjacent output spacing. A zero horizon returns the initial state without SciPy.
`order_parameter()`	`(rho) → float`	Kuramoto $R$	Return the mean transverse-expectation magnitude.
`purity()`	`(rho) → float`	$\text{Tr}(\rho^2)$	Return density-matrix purity.

The run() budget argument is a build-time override. If the solver is already built, its cached operators are reused. Invalid grids fail with ValueError, and an unsuccessful SciPy integration fails with RuntimeError.

`run()` Return Value

{
    "times": np.ndarray,      # shape (n_samples,)
    "R": np.ndarray,          # Kuramoto R at each sample, shape (n_samples,)
    "purity": np.ndarray,     # Tr(ρ²) at each sample, shape (n_samples,)
    "rho_final": np.ndarray,  # final density matrix, shape (dim, dim)
}

The initial density matrix is a pure product state obtained by applying one $R_y(\omega_i \bmod 2\pi)$ rotation to each qubit of the all-zero state.

Tutorial: Open-System Kuramoto Synchronisation

Step 1: Set Up the System

import numpy as np
from scpn_quantum_control.phase.lindblad import LindbladKuramotoSolver

# 4-oscillator chain with exponentially decaying coupling
n = 4
K = 0.45 * np.exp(-0.3 * np.abs(np.subtract.outer(range(n), range(n))))
np.fill_diagonal(K, 0.0)
omega = np.linspace(0.8, 1.2, n)

Step 2: Closed-System Baseline

solver_closed = LindbladKuramotoSolver(n, K, omega, gamma_amp=0.0, gamma_deph=0.0)
result_closed = solver_closed.run(t_max=2.0, dt=0.05)

print(f"Closed system — R: {result_closed['R'][0]:.3f} → {result_closed['R'][-1]:.3f}")
print(f"Purity: {result_closed['purity'][-1]:.6f}")  # should be 1.000000

Step 3: Add Dissipation

solver_open = LindbladKuramotoSolver(n, K, omega, gamma_amp=0.05, gamma_deph=0.02)
result_open = solver_open.run(t_max=2.0, dt=0.05)

print(f"Open system — R: {result_open['R'][0]:.3f} → {result_open['R'][-1]:.3f}")
print(f"Purity: {result_open['purity'][0]:.3f} → {result_open['purity'][-1]:.3f}")

Step 4: Compare

import matplotlib.pyplot as plt

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))

ax1.plot(result_closed['times'], result_closed['R'], label='Closed')
ax1.plot(result_open['times'], result_open['R'], label='Open (γ=0.05)')
ax1.set_xlabel('Time')
ax1.set_ylabel('R')
ax1.legend()
ax1.set_title('Synchronisation Order Parameter')

ax2.plot(result_closed['times'], result_closed['purity'], label='Closed')
ax2.plot(result_open['times'], result_open['purity'], label='Open')
ax2.set_xlabel('Time')
ax2.set_ylabel('Tr(ρ²)')
ax2.legend()
ax2.set_title('Purity')

plt.tight_layout()
plt.savefig('lindblad_comparison.png', dpi=150)

Examples

Strong Damping Kills Synchronisation

solver_strong = LindbladKuramotoSolver(n, K, omega, gamma_amp=0.5)
result_strong = solver_strong.run(t_max=5.0, dt=0.1)
print(f"R(T=5) = {result_strong['R'][-1]:.4f}")  # → near 0
print(f"Purity(T=5) = {result_strong['purity'][-1]:.4f}")  # → near 1/2^n

Dephasing Only (No Energy Relaxation)

solver_deph = LindbladKuramotoSolver(n, K, omega, gamma_amp=0.0, gamma_deph=0.1)
result_deph = solver_deph.run(t_max=2.0, dt=0.05)
# Populations unchanged, but coherences decay

Verify Density Matrix Properties

rho = result_open['rho_final']
assert np.allclose(np.trace(rho), 1.0), "Trace not preserved"
assert np.allclose(rho, rho.conj().T), "Not Hermitian"
eigenvalues = np.linalg.eigvalsh(rho)
assert np.all(eigenvalues >= -1e-12), "Not positive semidefinite"

Differentiable Objective Evidence

Bounded open-system objective rows are available through scpn_quantum_control.phase.open_system_objectives. The suite evaluates small Kuramoto-XY Lindblad objectives through LindbladKuramotoSolver.run() and certifies the final density matrix before accepting the objective row:

from scpn_quantum_control.phase import run_open_system_objective_suite


suite = run_open_system_objective_suite(backends=("lindblad_density",))
record = suite.records[0]
print(record.gradient)
print(record.invariant_certificate)

The trainable parameters are bounded scalar coupling and damping scales. The recorded gradient is a deterministic central finite difference, so it is useful for local objective diagnostics and reviewer replay, not an adjoint Lindblad gradient or a provider/hardware gradient. The committed evidence artifact is data/differentiable_phase_qnode/open_system_objective_evidence_20260709.json; regenerate it with scpn-bench open-system-objective-evidence --no-diff or scripts/export_open_system_objective_evidence.py.

Comparison with Other Tools

Feature	This module	QuTiP `mesolve`	MISTIQS
Lindblad equation	Yes	Yes	No
Hamiltonian	Kuramoto-XY (built-in)	Any (user-supplied)	TFIM only
Coupling matrix	Arbitrary $K_{ij}$	Any	Nearest-neighbour
Solver	`scipy.solve_ivp` (RK45)	Internal ODE solver	—
GPU	No	No (QuTiP 5: CuPy)	No
Output	$R(t)$ , purity, $\rho$	Arbitrary expect	—

When to use this module: You want open-system dynamics for the Kuramoto-XY Hamiltonian with the SCPN coupling matrix $K_{nm}$ , integrated with the rest of the scpn-quantum-control pipeline.

When to use QuTiP: You need arbitrary Hamiltonians, Floquet theory, stochastic Schrödinger equation, or QuTiP's extensive toolbox. Our module is not a QuTiP replacement — it is a specialised solver for one Hamiltonian.

Scaling

$n$	Hilbert space dim	Density matrix size	Memory (complex128)
4	16	16 × 16	4 KB
8	256	256 × 256	1 MB
10	1,024	1,024 × 1,024	16 MB
12	4,096	4,096 × 4,096	256 MB
14	16,384	16,384 × 16,384	4 GB

Beyond $n = 12$ , wall-time becomes the bottleneck (the RHS evaluation at each time step is $O(d^2)$ where $d = 2^n$ ). For larger systems, use the MCWF method (phase/tensor_jump.py) which evolves state vectors instead of density matrices.

References

Lindblad, G. "On the generators of quantum dynamical semigroups." Commun. Math. Phys. 48, 119–130 (1976).
Gorini, V., Kossakowski, A. & Sudarshan, E. C. G. "Completely positive dynamical semigroups of N-level systems." J. Math. Phys. 17, 821 (1976).
Ameri, V. et al. "Mutual information as an order parameter for quantum synchronization." PRA 91, 012301 (2015).
Giorgi, G. L. et al. "Quantum correlations and mutual synchronization." PRA 85, 052101 (2012).