Lindblad Master Equation Solver

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.

---

Theory

The Lindblad Master Equation

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}\)	Population transfer in the stated local basis
Pure dephasing	\(L_k = \sqrt{\gamma_\text{deph}/2}\, \sigma^z_k\)	\(\gamma_\text{deph}\)	Phase randomisation (\(T_2\) decay)

The amplitude channel changes populations and reduces transverse coherences according to its jump-basis convention below. Dephasing destroys off-diagonal coherences without changing populations.

The existing angular-momentum convention is sigma_minus = [[0, 0], [1, 0]]: its local jump maps basis state |0> to |1>. Ground-state identification depends on the Hamiltonian sign; this is not a claim that the channel relaxes every qubit to the computational |0>.

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	Optional positive dense-workspace budget in GiB.

If max_dense_gib is omitted, dense allocation uses SCPN_MAX_DENSE_GIB when set, otherwise the shared default. That default is the smaller of the host's free memory and any cgroup headroom this process still has, so inside a memory-limited container it follows the container's allowance rather than the machine's. 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, initial_density_matrix=None, atol=1e-8, rtol=1e-6) → dict	See below	Evolve through t_max; dt bounds adjacent output spacing. A zero horizon returns the admitted 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.

Every run() checks its active dense allowance, including when operators are cached. Its execution reservation includes the density workspace, integration history and returned scalar histories before their arrays are materialised. Unrepresentable or over-budget histories raise DenseAllocationError; refusal leaves prior cached operators intact. The shared reservation coordinates cooperating calls in this process; it does not guarantee absence of native or external memory pressure.

initial_density_matrix accepts a float64 or complex128 array of shape (2**n, 2**n) in the Hamiltonian's Qiskit little-endian basis order. Admission requires finite values, Hermiticity, unit trace and positive semidefiniteness within absolute 1e-10. Values are copied exactly: there is no trace normalisation, eigenvalue projection or basis reordering, and the caller's array is not mutated. Invalid inputs raise ValueError before building dense operators. The shared admission function is scpn_quantum_control.phase.density_input.validate_density_matrix.

atol and rtol are positive finite SciPy integration tolerances. dt controls output sampling, not the adaptive integrator's internal maximum step or error bound. Invalid grids or tolerances fail with ValueError, and an unsuccessful SciPy integration fails with RuntimeError.

Physical admission constrains the initial matrix. The returned matrix remains the unprojected SciPy solution; tighten integration tolerances when checking trace or positivity at a fixed reference precision. Input validation does not repair numerical drift in the output.

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)
}

When initial_density_matrix is omitted, the legacy pure product seed uses \(R_y(\omega_i \bmod 2\pi)\) rotations with index zero as the most significant tensor factor. That existing seed ordering is preserved; it differs from the Hamiltonian's qubit-index convention. Supply an explicit matrix when comparing routes from the same state. The four result keys and positional arguments are unchanged.

Analytic single-qubit reference

With hbar = 1, the implemented Hamiltonian is H = -omega * Z. Setting omega = -Omega/2 therefore gives H = Omega * Z/2. A plus-state reference has Bloch coordinates X = cos(Omega*t) and Y = sin(Omega*t).

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

Omega, time = 1.3, 0.83
plus = np.array([1, 1], dtype=np.complex128) / np.sqrt(2)
rho = np.outer(plus, plus.conj())
solver = LindbladKuramotoSolver(1, np.zeros((1, 1)), np.array([-Omega / 2]))
result = solver.run(
    time, 0.1, initial_density_matrix=rho, atol=1e-12, rtol=1e-11
)
pauli_x = np.array([[0, 1], [1, 0]])
pauli_y = np.array([[0, -1j], [1j, 0]])
actual = [np.trace(p @ result["rho_final"]).real for p in (pauli_x, pauli_y)]
np.testing.assert_allclose(
    actual, [np.cos(Omega * time), np.sin(Omega * time)], atol=1e-10, rtol=1e-8
)

The direct reference tests compare this density route, Qiskit Trotter evolution and zero-noise sparse action from the same explicit state. A separate noncommuting two-qubit fixture checks norm conservation and first-/second-order product-formula convergence at 4, 8 and 16 repetitions against an independently specified Hamiltonian. Output-grid checks at dt = 0.2, 0.1, 0.05 retain their separate sampling meaning. These local correctness checks are not isolated timings or general stochastic-channel qualification.

---

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

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

---

See Also

• Open-System Hardware Circuits — ancilla Lindblad and MCWF for hardware execution
• Tensor Networks — MPS/DMRG for \(n > 16\)
• Symmetry Sectors — reduce Hilbert space before Lindblad