python v4.7

Qutip

Quantum Toolbox in Python for open quantum systems simulation

Quick install

pip install qutip

Background and History

QuTiP (Quantum Toolbox in Python) was started in 2011 by Robert Johansson and Paul Nation while they were researchers at RIKEN in Japan. The project was motivated by the lack of a comprehensive, open-source Python library for simulating open quantum systems. At the time, most quantum simulation tools were written in Matlab or Mathematica, and Johansson and Nation saw an opportunity to build something that leveraged Python’s growing scientific computing ecosystem (NumPy, SciPy, matplotlib). The first public release, QuTiP 1.0, appeared in 2012.

QuTiP was designed around the physics of open quantum systems rather than the gate-circuit model that dominates quantum computing frameworks. Its core focus is time evolution under Hamiltonians with dissipation, using the Lindblad master equation and Monte Carlo trajectory methods. This made it the standard tool for theoretical quantum optics and superconducting qubit research long before the current wave of quantum computing SDKs existed. Researchers modeling cavity QED systems, spin chains, quantum transport, and decoherence dynamics adopted QuTiP widely.

The project reached version 4.0 in 2017, which brought performance improvements and expanded features including Floquet theory support, stochastic solvers, and improved handling of time-dependent Hamiltonians. Development of QuTiP 5.0 began around 2021, with a focus on modernizing the backend through a new data layer that supports sparse, dense, and JAX-backed representations for quantum objects. QuTiP 5 also introduced a more modular solver architecture and improved GPU support.

As of 2025, QuTiP is one of the most cited quantum simulation libraries in the academic literature, with thousands of published papers referencing it. The project is maintained by a community of contributors and governed by an advisory board. It has over 1,700 GitHub stars and an active mailing list. While QuTiP is not a quantum computing framework in the circuit sense, it occupies an essential role in the broader quantum software ecosystem as the primary tool for physics-level simulation of quantum devices, noise processes, and quantum control protocols. The qutip-qip subpackage adds basic quantum information processing capabilities, including circuit simulation and gate decomposition.

QuTiP (Quantum Toolbox in Python) is the standard Python library for simulating open quantum systems. Its central data type is the Qobj class, which represents kets, bras, density matrices, and operators in a unified way with automatic dimension tracking. Tensor products, partial traces, and operator arithmetic are all expressed naturally through Qobj methods, making it straightforward to build composite systems from subsystems.

Time evolution is handled by mesolve (Lindblad master equation for open systems) and sesolve (Schrodinger equation for closed systems). Both accept a time-dependent Hamiltonian and a list of collapse operators, and return a Result object containing state trajectories at each time step. The Lindblad master equation approach makes QuTiP particularly well suited for modeling dissipation, decoherence, and driven-dissipative systems. The Jaynes-Cummings model (a two-level atom coupled to a cavity mode) is a canonical example that QuTiP handles with only a few lines of code.

Visualization tools cover Bloch sphere rendering, Wigner function plots for continuous-variable states, and Husimi Q functions. QuTiP is widely used in quantum optics research, superconducting qubit theory, and quantum information courses because it operates at the physics level (Hamiltonians, collapse operators, observables) rather than the gate-circuit level. It is hardware-agnostic and focused entirely on numerical simulation.

Key Imports

import qutip as qt
from qutip import (
    basis, Qobj, ket2dm, tensor, destroy, create, sigmax, sigmay, sigmaz,
    qeye, mesolve, sesolve, mcsolve, Options, expect, fidelity,
    wigner, hinton, plot_bloch_sphere
)

State construction: basis, Qobj, ket2dm, tensor for building kets, density matrices, and composite systems. Operators: destroy, create (bosonic ladder operators), sigmax, sigmay, sigmaz (Pauli matrices), qeye (identity). Solvers: mesolve (Lindblad master equation), sesolve (Schrodinger equation), mcsolve (Monte Carlo trajectories). Analysis: expect (expectation values), fidelity (state comparison), Options (solver configuration). Visualization: wigner (Wigner function computation), hinton (Hinton diagram for matrices), plot_bloch_sphere.

Core Objects

The Qobj class is the foundation of QuTiP. Every quantum object, whether a ket, bra, operator, or density matrix, is represented as a Qobj instance. QuTiP infers the type from the object’s dimensions: a column vector is a ket, a row vector is a bra, a square matrix is an operator, and a Hermitian positive-semidefinite square matrix can serve as a density matrix.

Dimensions are tracked as a nested list [[row_dims], [col_dims]]. A single qubit ket has dims=[[2],[1]], while a two-qubit operator has dims=[[2,2],[2,2]]. This dimension tracking enables automatic validation of tensor products, partial traces, and operator-state multiplication.

Key creation functions:

  • basis(N, n) creates the Fock state |n> in an N-dimensional Hilbert space. For a qubit, basis(2, 0) is |0> and basis(2, 1) is |1>.
  • coherent(N, alpha) creates a coherent state with amplitude alpha in an N-dimensional truncated Fock space.
  • thermal_dm(N, n_avg) creates a thermal density matrix with mean occupation number n_avg.
  • ket2dm(psi) converts a ket into its corresponding density matrix |psi><psi|.
  • tensor(A, B, ...) computes the tensor product of multiple quantum objects, used for building composite systems from subsystems.

Common Patterns

Two-level system (qubit) under Lindblad dissipation

This example models a spin-1/2 particle in a magnetic field with energy relaxation (T1 decay). The Hamiltonian drives coherent precession while the collapse operator introduces irreversible decay from |1> to |0>.

import numpy as np
import qutip as qt

# Define qubit basis states
psi0 = qt.basis(2, 1)  # Start in excited state |1>

# Hamiltonian: magnetic field along x (drives Rabi oscillations)
omega = 2.0 * np.pi  # Rabi frequency
H = 0.5 * omega * qt.sigmax()

# Collapse operator for T1 decay (spontaneous emission)
gamma = 0.1  # Decay rate (1/T1)
c_ops = [np.sqrt(gamma) * qt.destroy(2)]

# Time array
tlist = np.linspace(0, 20, 500)

# Solve the Lindblad master equation
result = qt.mesolve(H, psi0, tlist, c_ops, e_ops=[qt.sigmaz(), qt.sigmax()])

# result.expect[0] contains <sigma_z>(t), result.expect[1] contains <sigma_x>(t)
import matplotlib.pyplot as plt

plt.figure(figsize=(8, 4))
plt.plot(tlist, result.expect[0], label=r"$\langle\sigma_z\rangle$")
plt.plot(tlist, result.expect[1], label=r"$\langle\sigma_x\rangle$")
plt.xlabel("Time")
plt.ylabel("Expectation value")
plt.legend()
plt.title("Qubit Rabi oscillations with T1 decay")
plt.show()

Jaynes-Cummings model

The Jaynes-Cummings model describes a two-level atom coupled to a single cavity mode. The Hamiltonian includes the free cavity energy, the atomic energy splitting, and the dipole interaction term. Tensor products combine the cavity Hilbert space with the qubit space.

import numpy as np
import qutip as qt

# Parameters
N = 15           # Cavity Fock space truncation
wc = 1.0         # Cavity frequency
wa = 1.0         # Atom frequency (on resonance)
g = 0.1          # Coupling strength

# Cavity operators (tensor with 2-level identity)
a = qt.tensor(qt.destroy(N), qt.qeye(2))

# Atom operators (tensor with cavity identity)
sm = qt.tensor(qt.qeye(N), qt.destroy(2))  # Atomic lowering operator

# Jaynes-Cummings Hamiltonian
H = wc * a.dag() * a + wa * sm.dag() * sm + g * (a.dag() * sm + a * sm.dag())

# Initial state: cavity in vacuum, atom excited
psi0 = qt.tensor(qt.basis(N, 0), qt.basis(2, 1))

# Time evolution
tlist = np.linspace(0, 100, 1000)
result = qt.mesolve(H, psi0, tlist, [], e_ops=[a.dag() * a, sm.dag() * sm])

import matplotlib.pyplot as plt

plt.figure(figsize=(8, 4))
plt.plot(tlist, result.expect[0], label="Cavity photon number")
plt.plot(tlist, result.expect[1], label="Atom excited state population")
plt.xlabel("Time")
plt.ylabel("Expectation value")
plt.legend()
plt.title("Jaynes-Cummings vacuum Rabi oscillations")
plt.show()

Bloch sphere visualization

QuTiP’s Bloch class renders qubit states on the Bloch sphere. You can add state vectors, points, or trajectories.

import qutip as qt

b = qt.Bloch()

# Add several states
up = qt.basis(2, 0)
down = qt.basis(2, 1)
plus = (up + down).unit()
minus = (up - down).unit()

b.add_states([up, down, plus, minus])
b.show()

Fidelity and quantum state overlap

The fidelity function computes the fidelity between two quantum states or density matrices. For pure states, fidelity equals the absolute value of their inner product. For mixed states, it uses the Uhlmann fidelity F(rho, sigma) = (Tr[sqrt(sqrt(rho) sigma sqrt(rho))])^2.

import qutip as qt

# Two pure states
psi1 = qt.basis(2, 0)
psi2 = (qt.basis(2, 0) + qt.basis(2, 1)).unit()

# Fidelity between kets
F = qt.fidelity(psi1, psi2)
print(f"Fidelity: {F:.4f}")  # ~0.7071 for |0> vs |+>

# Fidelity between density matrices
rho = qt.ket2dm(psi1)
sigma = 0.9 * qt.ket2dm(psi1) + 0.1 * qt.ket2dm(qt.basis(2, 1))  # Mixed state
F_mixed = qt.fidelity(rho, sigma)
print(f"Fidelity (mixed): {F_mixed:.4f}")

Wigner function

The Wigner function provides a phase-space representation of quantum states. It is especially useful for visualizing non-classical features such as interference fringes in cat states.

import numpy as np
import qutip as qt
import matplotlib.pyplot as plt

N = 30  # Fock space dimension
alpha = 2.0

# Create a cat state: (|alpha> + |-alpha>) / sqrt(2)
cat = (qt.coherent(N, alpha) + qt.coherent(N, -alpha)).unit()

# Compute Wigner function over phase space
xvec = np.linspace(-5, 5, 200)
W = qt.wigner(cat, xvec, xvec)

plt.figure(figsize=(6, 5))
plt.contourf(xvec, xvec, W, 100, cmap="RdBu_r")
plt.colorbar(label="W(x, p)")
plt.xlabel("x")
plt.ylabel("p")
plt.title("Wigner function of a cat state")
plt.show()

QuTiP vs Gate-Model Frameworks

QuTiP and gate-model frameworks like Qiskit, Cirq, and PennyLane solve fundamentally different problems. QuTiP operates at the physics level: you define Hamiltonians, collapse operators, and observables, then simulate continuous time evolution of quantum states including dissipation and decoherence. Gate-model frameworks operate at the circuit level: you compose sequences of discrete unitary gates and measure in the computational basis. If your goal is to model how a physical qubit decoheres under noise, simulate cavity QED dynamics, or design quantum control pulses, QuTiP is the right tool. If your goal is to implement Grover’s algorithm, run a variational circuit on real hardware, or compile a quantum program to native gates, use a gate-model framework. The qutip-qip subpackage bridges these worlds by providing circuit simulation within QuTiP, letting you model gate-level quantum information processing while retaining access to QuTiP’s noise and decoherence tools.

Learning Resources