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.
Installation uv pip install qutip
Optional packages for additional functionality:
Quantum information processing (circuits, gates)
uv pip install qutip-qip
Quantum trajectory viewer
uv pip install qutip-qtrl
Quick Start 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:
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.
- Time Evolution and Dynamics
Multiple solvers for different scenarios:
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.
- Analysis and Measurement
Compute physical quantities:
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.
- Visualization
Visualize quantum states and dynamics:
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.
- Advanced Methods
Specialized techniques for complex scenarios:
Floquet theory (periodic Hamiltonians)
T = 2 * np.pi / w_drive f_modes, f_energies = floquet_modes(H, T, args) result = fmmesolve(H, psi0, tlist, c_ops, T=T, args=args)
HEOM (non-Markovian, strong coupling)
from qutip.nonmarkov.heom import HEOMSolver, BosonicBath bath = BosonicBath(Q, ck_real, vk_real) hsolver = HEOMSolver(H_sys, [bath], max_depth=5) result = hsolver.run(rho0, tlist)
Permutational invariance (identical particles)
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
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
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
def compute_concurrence(t, psi): rho = ket2dm(psi) if psi.isket else psi return concurrence(rho)
tlist = np.linspace(0, 10, 100) result = mesolve(qeye([2, 2]), psi0, tlist, c_ops)
Compute concurrence for each state
C_t = [concurrence(state.proj()) 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
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 Truncate Hilbert spaces: Use smallest dimension that captures dynamics Choose appropriate solver: sesolve for pure states is faster than mesolve Time-dependent terms: String format (e.g., 'cos(w*t)') is fastest Store only needed data: Use e_ops instead of storing all states Adjust tolerances: Balance accuracy with computation time via Options Parallel trajectories: mcsolve automatically uses multiple CPUs 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