Alterlab Qutip

---
name: alterlab-qutip
description: Simulates open quantum systems with QuTiP, the Quantum Toolbox in Python, solving Lindblad master equations (mesolve), Monte Carlo trajectories (mcsolve), and unitary dynamics (sesolve). Use when studying master-equation or Lindblad dynamics, decoherence, dissipation, quantum optics, cavity QED, or open-system time evolution. NOT for circuit-based quantum computing or hardware execution — for IBM Quantum circuits prefer alterlab-qiskit, for Google Quantum AI or NISQ circuits prefer alterlab-cirq, and for gradient-trained quantum ML prefer alterlab-pennylane. Part of the AlterLab Academic Skills suite.
license: MIT
allowed-tools: Read Write Edit Bash(python:*)
compatibility: No API key required. Runs locally via `uv run python`; requires the qutip Python package.
metadata:
    skill-author: AlterLab
    version: "1.2.0"
    last_updated: "2026-09-23"
---

# QuTiP: Quantum Toolbox in Python

## Overview

QuTiP provides comprehensive tools for simulating and analyzing quantum mechanical systems. It handles both closed (unitary) and open (dissipative) quantum systems with multiple solvers optimized for different scenarios.

## When to Use This Skill

Use this skill when the user wants to:
- Solve Lindblad master equations, quantum trajectories, or Schrödinger dynamics (`mesolve`, `mcsolve`, `sesolve`)
- Model decoherence, dissipation, cavity QED, or driven/periodic systems (Floquet, Bloch-Redfield, HEOM)
- Compute steady states, correlation functions, spectra, entanglement, or phase-space distributions

### Does NOT Trigger

| Scenario | Use Instead |
|----------|-------------|
| Gate-based circuits on IBM hardware, transpilation, Runtime primitives | `alterlab-qiskit` |
| Cirq circuits, Google hardware / Quantum Virtual Machine, XEB/RB | `alterlab-cirq` |
| Training variational circuits with autodiff (quantum ML, differentiable VQE) | `alterlab-pennylane` |
| Symbolic derivation of operators or commutators without numerics | `alterlab-sympy` |

## Installation

```bash
uv pip install qutip
```

This skill targets **QuTiP 5.x** (verified against 5.3.1, current as of 2026-09; Python ≥ 3.11). Key v5 API changes from v4: solver options are a plain `dict` (the old `Options`/`solver.Options` class is gone), stochastic solvers use `heterodyne=True/False` instead of a `noise` code, `progress_bar` takes a string (`'text'`/`'tqdm'`/`''`) rather than `True`, `mcsolve` runs trajectories serially unless you pass `options={"map": "parallel"}`, and the standalone `floquet_modes()` family is deprecated in favour of `FloquetBasis`.

String-format time dependence (`'cos(w*t)'`) is compiled with Cython only when the `runtime-compilation` extra is installed; otherwise QuTiP warns and falls back to slower `eval`:

```bash
uv pip install "qutip[runtime-compilation]"   # Cython + setuptools + filelock
```

Optional packages for additional functionality:

```bash
# Quantum information processing (circuits, gates)
uv pip install qutip-qip

# Quantum optimal control (GRAPE/CRAB pulse optimization)
uv pip install qutip-qtrl
```

## Quick Start

```python
from qutip import *
import numpy as np
import matplotlib.pyplot as plt

# Create quantum state
psi = basis(2, 0)  # |0⟩ state

# Create operator
H = sigmaz()  # Hamiltonian

# Time evolution
tlist = np.linspace(0, 10, 100)
result = sesolve(H, psi, tlist, e_ops=[sigmaz()])

# Plot results
plt.plot(tlist, result.expect[0])
plt.xlabel('Time')
plt.ylabel('⟨σz⟩')
plt.show()
```

## Core Capabilities

### 1. Quantum Objects and States

Create and manipulate quantum states and operators:

```python
# States
psi = basis(N, n)  # Fock state |n⟩
psi = coherent(N, alpha)  # Coherent state |α⟩
rho = thermal_dm(N, n_avg)  # Thermal density matrix

# Operators
a = destroy(N)  # Annihilation operator
H = num(N)  # Number operator
sx, sy, sz = sigmax(), sigmay(), sigmaz()  # Pauli matrices

# Composite systems
psi_AB = tensor(psi_A, psi_B)  # Tensor product
```

**See** `references/core_concepts.md` for comprehensive coverage of quantum objects, states, operators, and tensor products.

### 2. Time Evolution and Dynamics

Multiple solvers for different scenarios:

```python
# Closed systems (unitary evolution)
result = sesolve(H, psi0, tlist, e_ops=[num(N)])

# Open systems (dissipation)
c_ops = [np.sqrt(0.1) * destroy(N)]  # Collapse operators
result = mesolve(H, psi0, tlist, c_ops, e_ops=[num(N)])

# Quantum trajectories (Monte Carlo)
result = mcsolve(H, psi0, tlist, c_ops, ntraj=500, e_ops=[num(N)])
```

**Solver selection guide:**
- `sesolve`: Pure states, unitary evolution
- `mesolve`: Mixed states, dissipation, general open systems
- `mcsolve`: Quantum jumps, photon counting, individual trajectories
- `brmesolve`: Weak system-bath coupling
- `fmmesolve`: Time-periodic Hamiltonians (Floquet)

**See** `references/time_evolution.md` for detailed solver documentation, time-dependent Hamiltonians, and advanced options.

### 3. Analysis and Measurement

Compute physical quantities:

```python
# Expectation values
n_avg = expect(num(N), psi)

# Entropy measures
S = entropy_vn(rho)  # Von Neumann entropy
C = concurrence(rho)  # Entanglement (two qubits)

# Fidelity and distance
F = fidelity(psi1, psi2)
D = tracedist(rho1, rho2)

# Correlation functions
corr = correlation_2op_1t(H, rho0, taulist, c_ops, A, B)
w, S = spectrum_correlation_fft(taulist, corr)

# Steady states
rho_ss = steadystate(H, c_ops)
```

**See** `references/analysis.md` for entropy, fidelity, measurements, correlation functions, and steady state calculations.

### 4. Visualization

Visualize quantum states and dynamics:

```python
# Bloch sphere
b = Bloch()
b.add_states(psi)
b.show()

# Wigner function (phase space)
xvec = np.linspace(-5, 5, 200)
W = wigner(psi, xvec, xvec)
plt.contourf(xvec, xvec, W, 100, cmap='RdBu')

# Fock distribution
plot_fock_distribution(psi)

# Matrix visualization
hinton(rho)  # Hinton diagram
matrix_histogram(H.full())  # 3D bars
```

**See** `references/visualization.md` for Bloch sphere animations, Wigner functions, Q-functions, and matrix visualizations.

### 5. Advanced Methods

Specialized techniques for complex scenarios:

```python
# Floquet theory (periodic Hamiltonians); floquet_modes() is deprecated in v5
T = 2 * np.pi / w_drive
fbasis = FloquetBasis(H, T, args=args)
f_energies = fbasis.e_quasi   # quasi-energies
f_modes = fbasis.mode(0)      # Floquet modes at t = 0
result = fmmesolve(H, psi0, tlist, c_ops, T=T, args=args)

# HEOM (non-Markovian, strong coupling)
from qutip.solver.heom import HEOMSolver, BosonicBath
# Bath correlation as a sum of exponentials: ck_real, vk_real, ck_imag, vk_imag
bath = BosonicBath(Q, ck_real, vk_real, ck_imag, vk_imag)
hsolver = HEOMSolver(H_sys, [bath], max_depth=5)
result = hsolver.run(rho0, tlist)

# Permutational invariance (identical particles)
from qutip.piqs.piqs import dicke, jspin  # not exported at qutip.piqs top level
psi = dicke(N, j, m)  # Dicke states
Jz = jspin(N, 'z')  # Collective operators
```

**See** `references/advanced.md` for Floquet theory, HEOM, permutational invariance, stochastic solvers, superoperators, and performance optimization.

## Common Workflows

### Simulating a Damped Harmonic Oscillator

```python
# System parameters
N = 20  # Hilbert space dimension
omega = 1.0  # Oscillator frequency
kappa = 0.1  # Decay rate

# Hamiltonian and collapse operators
H = omega * num(N)
c_ops = [np.sqrt(kappa) * destroy(N)]

# Initial state
psi0 = coherent(N, 3.0)

# Time evolution
tlist = np.linspace(0, 50, 200)
result = mesolve(H, psi0, tlist, c_ops, e_ops=[num(N)])

# Visualize
plt.plot(tlist, result.expect[0])
plt.xlabel('Time')
plt.ylabel('⟨n⟩')
plt.title('Photon Number Decay')
plt.show()
```

### Two-Qubit Entanglement Dynamics

```python
# Create Bell state
psi0 = bell_state('00')

# Local dephasing on each qubit
gamma = 0.1
c_ops = [
    np.sqrt(gamma) * tensor(sigmaz(), qeye(2)),
    np.sqrt(gamma) * tensor(qeye(2), sigmaz())
]

# Track entanglement
tlist = np.linspace(0, 10, 100)
result = mesolve(qeye([2, 2]), psi0, tlist, c_ops)

# With collapse operators, result.states are density matrices (type 'oper'),
# so pass them straight to concurrence() — do NOT call .proj() on them.
C_t = [concurrence(state) for state in result.states]

plt.plot(tlist, C_t)
plt.xlabel('Time')
plt.ylabel('Concurrence')
plt.title('Entanglement Decay')
plt.show()
```

### Jaynes-Cummings Model

```python
# System parameters
N = 10  # Cavity Fock space
wc = 1.0  # Cavity frequency
wa = 1.0  # Atom frequency
g = 0.05  # Coupling strength

# Operators
a = tensor(destroy(N), qeye(2))  # Cavity
sm = tensor(qeye(N), sigmam())  # Atom

# Hamiltonian (RWA)
H = wc * a.dag() * a + wa * sm.dag() * sm + g * (a.dag() * sm + a * sm.dag())

# Initial state: cavity in coherent state, atom in ground state
psi0 = tensor(coherent(N, 2), basis(2, 0))

# Dissipation
kappa = 0.1  # Cavity decay
gamma = 0.05  # Atomic decay
c_ops = [np.sqrt(kappa) * a, np.sqrt(gamma) * sm]

# Observables
n_cav = a.dag() * a
n_atom = sm.dag() * sm

# Evolve
tlist = np.linspace(0, 50, 200)
result = mesolve(H, psi0, tlist, c_ops, e_ops=[n_cav, n_atom])

# Plot
fig, axes = plt.subplots(2, 1, figsize=(8, 6), sharex=True)
axes[0].plot(tlist, result.expect[0])
axes[0].set_ylabel('⟨n_cavity⟩')
axes[1].plot(tlist, result.expect[1])
axes[1].set_ylabel('⟨n_atom⟩')
axes[1].set_xlabel('Time')
plt.tight_layout()
plt.show()
```

## Tips for Efficient Simulations

1. **Truncate Hilbert spaces**: Use smallest dimension that captures dynamics
2. **Choose appropriate solver**: `sesolve` for pure states is faster than `mesolve`
3. **Time-dependent terms**: String format (e.g., `'cos(w*t)'`) is fastest once compiled (needs `qutip[runtime-compilation]`)
4. **Store only needed data**: Use `e_ops` instead of storing all states
5. **Adjust tolerances**: Balance accuracy with computation time via the `options` dict (`atol`/`rtol`/`nsteps`)
6. **Parallel trajectories**: `mcsolve` is serial by default in v5 — pass `options={"map": "parallel", "num_cpus": 4}` (or `"loky"` / `"mpi"`)
7. **Check convergence**: Vary `ntraj`, Hilbert space size, and tolerances

## Troubleshooting

**Memory issues**: Reduce Hilbert space dimension, use `store_final_state` option, or consider Krylov methods

**Slow simulations**: Use string-based time-dependence, increase tolerances slightly, or try `method='bdf'` for stiff problems

**Numerical instabilities**: Decrease time steps (`nsteps` option), increase tolerances, or check Hamiltonian/operators are properly defined

**Import errors**: Ensure QuTiP is installed correctly; quantum gates require `qutip-qip` package

## References

This skill includes detailed reference documentation:

- **`references/core_concepts.md`**: Quantum objects, states, operators, tensor products, composite systems
- **`references/time_evolution.md`**: All solvers (sesolve, mesolve, mcsolve, brmesolve, etc.), time-dependent Hamiltonians, solver options
- **`references/visualization.md`**: Bloch sphere, Wigner functions, Q-functions, Fock distributions, matrix plots
- **`references/analysis.md`**: Expectation values, entropy, fidelity, entanglement measures, correlation functions, steady states
- **`references/advanced.md`**: Floquet theory, HEOM, permutational invariance, stochastic methods, superoperators, performance tips

## External Resources

- Documentation: https://qutip.readthedocs.io/
- Tutorials: https://qutip.org/qutip-tutorials/
- API Reference: https://qutip.readthedocs.io/en/stable/apidoc/apidoc.html
- GitHub: https://github.com/qutip/qutip

Part of the AlterLab Academic Skills suite.