QuTiP - Quantum Dynamics

QuTiP is designed to be a flexible and efficient framework for quantum mechanics. It allows for easy creation of quantum objects and provides powerful solvers to track their evolution in both closed and open environments.

When to Use

• Simulating the time evolution of a quantum system (Schrödinger or Master equation)
• Calculating steady states of open quantum systems (e.g., a cavity under drive and dissipation)
• Analyzing entanglement, Wigner functions, and quantum correlations
• Studying light-matter interactions (Jaynes-Cummings model, Rabi oscillations)
• Pulse sequence optimization and quantum control
• Calculating spectrums and multi-time correlation functions
• Parallelizing quantum simulations across multiple CPUs

Reference Documentation

Official docs: https://qutip.org/docs/latest/
Tutorials: https://qutip.org/tutorials.html
Search patterns: qutip.Qobj, qutip.mesolve, qutip.wigner, qutip.expect, qutip.basis

Core Principles

The Qobj (Quantum Object)

The central data structure. A Qobj can represent a state (ket or bra), an operator, or a superoperator. It automatically handles dimensions and validates operations (e.g., preventing the addition of a ket and an operator).

Composite Systems

QuTiP uses the Kronecker product (tensor) to represent systems composed of multiple sub-systems (e.g., a qubit coupled to a resonator).

Solver Workflow

1. Define the Hamiltonian (H)
2. Define collapse operators (C_n) for dissipation
3. Set initial state (ρ₀)
4. Define time sequence (t)
5. Run the solver (mesolve, sesolve, etc.)

Quick Reference

Installation

pip install qutip

Standard Imports

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

Basic Pattern - Rabi Oscillations

import qutip as qt
import numpy as np

# 1. Define Hamiltonian: H = h_bar * omega * sigma_x / 2
omega = 1.0 * 2 * np.pi
H = 0.5 * omega * qt.sigmax()

# 2. Initial state: Spin down (|1>)
psi0 = qt.basis(2, 1)

# 3. Time points
times = np.linspace(0.0, 10.0, 100)

# 4. Solve Schrödinger Equation
result = qt.sesolve(H, psi0, times, [qt.sigmaz()])

# 5. Get expectation values
sz_expt = result.expect[0]

Critical Rules

✅ DO

• Use basis(N, i) - Use built-in functions to create states instead of manual arrays
• Check Qobj Dimensions - Always verify .dims when working with composite systems to ensure tensor products are correct
• Use expect() - Let solvers calculate expectation values directly to save memory and time
• Specify c_ops for Open Systems - Use mesolve with collapse operators to model decoherence
• Vectorize Hamiltonians - If your Hamiltonian is time-dependent, use the [H0, [H1, 'sin(w*t)']] list format for speed
• Visualize with Wigner - Use qt.wigner(state, xvec, yvec) for a deep look into the quantum nature of the state

❌ DON'T

• Extract data manually - Avoid state.data.toarray() unless absolutely necessary; QuTiP's Qobj is optimized for sparse math
• Ignore solver warnings - Numerical instabilities often mean the time step is too large or the Hilbert space is truncated too early
• Mix Hilbert Spaces - Don't perform operations between objects with different dimensions without proper tensor products
• Forget Truncation - For infinite systems (like oscillators), ensure the Fock space size (N) is large enough for convergence

Anti-Patterns (NEVER)

import qutip as qt

# ❌ BAD: Creating states using raw numpy arrays (no metadata)
psi = qt.Qobj([[1], [0]])

# ✅ GOOD: Use built-in constructors
psi = qt.basis(2, 0) # Ket |0>

# ❌ BAD: Manual Kronecker product
H_total = np.kron(H_atom, np.eye(N_cav))

# ✅ GOOD: Use qt.tensor
H_total = qt.tensor(H_atom, qt.qeye(N_cav))

# ❌ BAD: Calculating expectation values in a Python loop
# (Extremely slow and memory intensive)

# ✅ GOOD: Pass operators to the solver
result = qt.mesolve(H, psi0, times, c_ops, [qt.sigmaz(), qt.sigmax()])

Quantum Objects (qt.Qobj)

Creation and Properties

# Kets and Bras
psi = qt.basis(5, 0) # Ket |0> in 5-level system
bra = psi.dag()      # Adjoint (Bra)

# Operators
a = qt.destroy(10) # Annihilation operator (size 10)
n = a.dag() * a    # Number operator
sz = qt.sigmaz()   # Pauli Z

# Properties
print(psi.type)    # 'ket'
print(sz.dims)     # [[2], [2]]
print(n.eigenstates()) # Get eigenvalues and vectors

Composite Systems (Tensors)

# Qubit (2) + Cavity (10)
sz = qt.sigmaz()
a = qt.destroy(10)

# Operator acting on qubit only
sz_total = qt.tensor(sz, qt.qeye(10))

# Interaction term: sigma_x * (a + a_dag)
H_int = qt.tensor(qt.sigmax(), a + a.dag())

Dynamics Solvers (qt.mesolve, qt.sesolve)

Master Equation (Open Systems)

# Parameters
kappa = 0.1 # Cavity decay rate
H = qt.tensor(qt.sigmax(), qt.qeye(10)) # Hamiltonian
c_ops = [np.sqrt(kappa) * qt.tensor(qt.qeye(2), a)] # Decay list
psi0 = qt.tensor(qt.basis(2, 0), qt.basis(10, 0))
times = np.linspace(0, 50, 200)

# Solve
result = qt.mesolve(H, psi0, times, c_ops, [qt.tensor(qt.sigmaz(), qt.qeye(10))])

# Access expectation values
exp_z = result.expect[0]

Time-Dependent Hamiltonians

# Using string-based syntax (fastest)
H0 = qt.sigmax()
H1 = qt.sigmaz()
args = {'w': 0.5}
H_t = [H0, [H1, 'sin(w * t)']]

result = qt.sesolve(H_t, qt.basis(2, 0), times, args=args)

Steady State and Correlation

# Steady state of a driven-dissipative system
rho_ss = qt.steadystate(H, c_ops)

# Expectation value in steady state
avg_n = qt.expect(a.dag() * a, rho_ss)

# Emission Spectrum
w_list = np.linspace(0, 10, 100)
spec = qt.spectrum(H, w_list, c_ops, a.dag(), a)

Visualization

Wigner and Bloch Sphere

# Wigner function for a Cat state
cat = (qt.coherent(20, 2) + qt.coherent(20, -2)).unit()
xvec = np.linspace(-5, 5, 100)
W = qt.wigner(cat, xvec, xvec)

# Bloch Sphere
b = qt.Bloch()
b.add_states(qt.basis(2, 0))
b.add_states((qt.basis(2,0) + qt.basis(2,1)).unit())
# b.show()

Practical Workflows

1. Jaynes-Cummings Model (Qubit-Cavity)

def simulate_jc(g, kappa, n_max, times):
    # Operators
    a = qt.tensor(qt.qeye(2), qt.destroy(n_max))
    sm = qt.tensor(qt.destroy(2), qt.qeye(n_max))
    
    # Hamiltonian
    H = g * (sm.dag() * a + sm * a.dag())
    
    # Dissipation
    c_ops = [np.sqrt(kappa) * a]
    
    # Initial: Excited atom, 0 photons
    psi0 = qt.tensor(qt.basis(2, 0), qt.basis(n_max, 0))
    
    res = qt.mesolve(H, psi0, times, c_ops, [sm.dag()*sm, a.dag()*a])
    return res.expect

2. Quantum State Tomography Analysis

def fidelity_check(state_ideal, state_measured):
    """Calculates fidelity between two quantum states."""
    # Works for both kets and density matrices
    return qt.fidelity(state_ideal, state_measured)

3. Calculating G2 Correlation Function

def calculate_g2(H, c_ops, a_op, times):
    """Second-order correlation function g2(tau)."""
    # Use built-in correlation function
    g2_tau = qt.correlation_2op_1t(H, None, times, c_ops, a_op.dag(), a_op)
    # Normalize by steady state intensity
    rho_ss = qt.steadystate(H, c_ops)
    n_ss = qt.expect(a_op.dag() * a_op, rho_ss)
    return np.real(g2_tau) / (n_ss**2)

Performance Optimization

Sparse vs Dense

QuTiP uses Scipy sparse matrices by default. For very small systems (N < 4), converting to dense might be faster, but for most N > 10, sparse is mandatory.

Parallel Map

Use qt.parallel_map for parameter sweeps.

def task(omega):
    H = 0.5 * omega * qt.sigmax()
    return qt.mesolve(H, psi0, times, c_ops, [qt.sigmaz()]).expect[0]

omegas = np.linspace(0, 10, 20)
results = qt.parallel_map(task, omegas)

Common Pitfalls and Solutions

The "Dimensions Mismatch" Error

When tensoring systems, the order must be consistent.

# ❌ Problem: H is (2x10) but Op is (10x2)
# ✅ Solution: Check dims and ensure consistent order in tensor()
print(H.dims) # [[2, 10], [2, 10]]

Insufficient Fock Space

If your photon number distribution hits the upper limit of your N, your results are invalid.

# ✅ Solution: Check populations
occ = qt.expect(qt.destroy(N).dag() * qt.destroy(N), final_state)
# If occ > 0.8 * N, increase N.

Hermitian Hamiltonian

Solvers assume H is Hermitian. If your time-dependent code results in a non-Hermitian H, sesolve will produce unphysical states (norm ≠ 1).

# ✅ Solution: Check norm
if not np.allclose(result.states[-1].norm(), 1.0):
    print("Warning: Norm not preserved!")

Best Practices

1. Always use built-in constructors (basis, coherent, etc.) instead of manual array creation
2. Verify dimensions with .dims before performing tensor products
3. Pass operators to solvers for expectation values rather than calculating them manually
4. Use appropriate Fock space truncation (N) - check that populations don't saturate
5. For time-dependent Hamiltonians, use the string-based format for best performance
6. Monitor solver warnings - they often indicate numerical issues
7. Use sparse matrices (default) for systems with N > 10
8. Check state normalization after time evolution, especially for custom solvers
9. Use parallel_map for parameter sweeps to leverage multiple CPUs
10. Visualize quantum states with Wigner functions and Bloch spheres for deeper insight

QuTiP is the standard tool for exploring the non-intuitive world of quantum dynamics. It bridges the gap between theoretical equations and numerical experimentation, providing physicists with a high-performance lab in a Python script.