claude-scientific-skills/scientific-skills/qutip/references/advanced.md

QuTiP Advanced Features

Floquet Theory

For time-periodic Hamiltonians H(t + T) = H(t).

Floquet Modes and Quasi-Energies

from qutip import *
import numpy as np

# Time-periodic Hamiltonian
w_d = 1.0  # Drive frequency
T = 2 * np.pi / w_d  # Period

H0 = sigmaz()
H1 = sigmax()
H = [H0, [H1, 'cos(w*t)']]
args = {'w': w_d}

# Calculate Floquet modes and quasi-energies
f_modes, f_energies = floquet_modes(H, T, args)

print("Quasi-energies:", f_energies)
print("Floquet modes:", f_modes)

Floquet States at Time t

# Get Floquet state at specific time
t = 1.0
f_states_t = floquet_states(f_modes, f_energies, t)

Floquet State Decomposition

# Decompose initial state in Floquet basis
psi0 = basis(2, 0)
f_coeff = floquet_state_decomposition(f_modes, f_energies, psi0)

Floquet-Markov Master Equation

# Time evolution with dissipation
c_ops = [np.sqrt(0.1) * sigmam()]
tlist = np.linspace(0, 20, 200)

result = fmmesolve(H, psi0, tlist, c_ops, e_ops=[sigmaz()], T=T, args=args)

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

Floquet Tensor

# Floquet tensor (generalized Bloch-Redfield)
A_ops = [[sigmaz(), lambda w: 0.1 * w if w > 0 else 0]]

# Build Floquet tensor
R, U = floquet_markov_mesolve(H, psi0, tlist, A_ops, e_ops=[sigmaz()],
                               T=T, args=args)

Effective Hamiltonian

# Time-averaged effective Hamiltonian
H_eff = floquet_master_equation_steadystate(H, c_ops, T, args)

Hierarchical Equations of Motion (HEOM)

For non-Markovian open quantum systems with strong system-bath coupling.

Basic HEOM Setup

from qutip import heom

# System Hamiltonian
H_sys = sigmaz()

# Bath correlation function (exponential)
Q = sigmax()  # System-bath coupling operator
ck_real = [0.1]  # Coupling strengths
vk_real = [0.5]  # Bath frequencies

# HEOM bath
bath = heom.BosonicBath(Q, ck_real, vk_real)

# Initial state
rho0 = basis(2, 0) * basis(2, 0).dag()

# Create HEOM solver
max_depth = 5
hsolver = heom.HEOMSolver(H_sys, [bath], max_depth=max_depth)

# Time evolution
tlist = np.linspace(0, 10, 100)
result = hsolver.run(rho0, tlist)

# Extract reduced system density matrix
rho_sys = [r.extract_state(0) for r in result.states]

Multiple Baths

# Define multiple baths
bath1 = heom.BosonicBath(sigmax(), [0.1], [0.5])
bath2 = heom.BosonicBath(sigmay(), [0.05], [1.0])

hsolver = heom.HEOMSolver(H_sys, [bath1, bath2], max_depth=5)

Drude-Lorentz Spectral Density

# Common in condensed matter physics
from qutip.nonmarkov.heom import DrudeLorentzBath

lam = 0.1  # Reorganization energy
gamma = 0.5  # Bath cutoff frequency
T = 1.0  # Temperature (in energy units)
Nk = 2  # Number of Matsubara terms

bath = DrudeLorentzBath(Q, lam, gamma, T, Nk)

HEOM Options

options = heom.HEOMSolver.Options(
    nsteps=2000,
    store_states=True,
    rtol=1e-7,
    atol=1e-9
)

hsolver = heom.HEOMSolver(H_sys, [bath], max_depth=5, options=options)

Permutational Invariance

For identical particle systems (e.g., spin ensembles).

Dicke States

from qutip import dicke

# Dicke state |j, m⟩ for N spins
N = 10  # Number of spins
j = N/2  # Total angular momentum
m = 0   # z-component

psi = dicke(N, j, m)

Permutation-Invariant Operators

from qutip.piqs import jspin

# Collective spin operators
N = 10
Jx = jspin(N, 'x')
Jy = jspin(N, 'y')
Jz = jspin(N, 'z')
Jp = jspin(N, '+')
Jm = jspin(N, '-')

PIQS Dynamics

from qutip.piqs import Dicke

# Setup Dicke model
N = 10
emission = 1.0
dephasing = 0.5
pumping = 0.0
collective_emission = 0.0

system = Dicke(N=N, emission=emission, dephasing=dephasing,
               pumping=pumping, collective_emission=collective_emission)

# Initial state
psi0 = dicke(N, N/2, N/2)  # All spins up

# Time evolution
tlist = np.linspace(0, 10, 100)
result = system.solve(psi0, tlist, e_ops=[Jz])

Non-Markovian Monte Carlo

Quantum trajectories with memory effects.

from qutip import nm_mcsolve

# Non-Markovian bath correlation
def bath_correlation(t1, t2):
    tau = abs(t2 - t1)
    return np.exp(-tau / 2.0) * np.cos(tau)

# System setup
H = sigmaz()
c_ops = [sigmax()]
psi0 = basis(2, 0)
tlist = np.linspace(0, 10, 100)

# Solve with memory
result = nm_mcsolve(H, psi0, tlist, c_ops, sc_ops=[],
                     bath_corr=bath_correlation, ntraj=500,
                     e_ops=[sigmaz()])

Stochastic Solvers with Measurements

Continuous Measurement

# Homodyne detection
sc_ops = [np.sqrt(0.1) * destroy(N)]  # Measurement operator

result = ssesolve(H, psi0, tlist, sc_ops=sc_ops,
                   e_ops=[num(N)], ntraj=100,
                   noise=11)  # 11 for homodyne

# Heterodyne detection
result = ssesolve(H, psi0, tlist, sc_ops=sc_ops,
                   e_ops=[num(N)], ntraj=100,
                   noise=12)  # 12 for heterodyne

Photon Counting

# Quantum jump times
result = mcsolve(H, psi0, tlist, c_ops, ntraj=50,
                 options=Options(store_states=True))

# Extract measurement times
for i, jump_times in enumerate(result.col_times):
    print(f"Trajectory {i} jump times: {jump_times}")
    print(f"Which operator: {result.col_which[i]}")

Krylov Subspace Methods

Efficient for large systems.

from qutip import krylovsolve

# Use Krylov solver
result = krylovsolve(H, psi0, tlist, krylov_dim=10, e_ops=[num(N)])

Bloch-Redfield Master Equation

For weak system-bath coupling.

# Bath spectral density
def ohmic_spectrum(w):
    if w >= 0:
        return 0.1 * w  # Ohmic
    else:
        return 0

# Coupling operators and spectra
a_ops = [[sigmax(), ohmic_spectrum]]

# Solve
result = brmesolve(H, psi0, tlist, a_ops, e_ops=[sigmaz()])

Temperature-Dependent Bath

def thermal_spectrum(w):
    # Bose-Einstein distribution
    T = 1.0  # Temperature
    if abs(w) < 1e-10:
        return 0.1 * T
    n_th = 1 / (np.exp(abs(w)/T) - 1)
    if w >= 0:
        return 0.1 * w * (n_th + 1)
    else:
        return 0.1 * abs(w) * n_th

a_ops = [[sigmax(), thermal_spectrum]]
result = brmesolve(H, psi0, tlist, a_ops, e_ops=[sigmaz()])

Superoperators and Quantum Channels

Superoperator Representations

# Liouvillian
L = liouvillian(H, c_ops)

# Convert between representations
from qutip import (spre, spost, sprepost,
                    super_to_choi, choi_to_super,
                    super_to_kraus, kraus_to_super)

# Superoperator forms
L_spre = spre(H)  # Left multiplication
L_spost = spost(H)  # Right multiplication
L_sprepost = sprepost(H, H.dag())

# Choi matrix
choi = super_to_choi(L)

# Kraus operators
kraus = super_to_kraus(L)

Quantum Channels

# Depolarizing channel
p = 0.1  # Error probability
K0 = np.sqrt(1 - 3*p/4) * qeye(2)
K1 = np.sqrt(p/4) * sigmax()
K2 = np.sqrt(p/4) * sigmay()
K3 = np.sqrt(p/4) * sigmaz()

kraus_ops = [K0, K1, K2, K3]
E = kraus_to_super(kraus_ops)

# Apply channel
rho_out = E * operator_to_vector(rho_in)
rho_out = vector_to_operator(rho_out)

Amplitude Damping

# T1 decay
gamma = 0.1
K0 = Qobj([[1, 0], [0, np.sqrt(1 - gamma)]])
K1 = Qobj([[0, np.sqrt(gamma)], [0, 0]])

E_damping = kraus_to_super([K0, K1])

Phase Damping

# T2 dephasing
gamma = 0.1
K0 = Qobj([[1, 0], [0, np.sqrt(1 - gamma/2)]])
K1 = Qobj([[0, 0], [0, np.sqrt(gamma/2)]])

E_dephasing = kraus_to_super([K0, K1])

Quantum Trajectories Analysis

Extract Individual Trajectories

options = Options(store_states=True, store_final_state=False)
result = mcsolve(H, psi0, tlist, c_ops, ntraj=100, options=options)

# Access individual trajectories
for i in range(len(result.states)):
    trajectory = result.states[i]  # List of states for trajectory i
    # Analyze trajectory

Trajectory Statistics

# Mean and standard deviation
result = mcsolve(H, psi0, tlist, c_ops, e_ops=[num(N)], ntraj=500)

n_mean = result.expect[0]
n_std = result.std_expect[0]

# Photon number distribution at final time
final_states = [result.states[i][-1] for i in range(len(result.states))]

Time-Dependent Terms Advanced

QobjEvo

from qutip import QobjEvo

# Time-dependent Hamiltonian with QobjEvo
def drive(t, args):
    return args['A'] * np.exp(-t/args['tau']) * np.sin(args['w'] * t)

H0 = num(N)
H1 = destroy(N) + create(N)
args = {'A': 1.0, 'w': 1.0, 'tau': 5.0}

H_td = QobjEvo([H0, [H1, drive]], args=args)

# Can update args without recreating
H_td.arguments({'A': 2.0, 'w': 1.5, 'tau': 10.0})

Compiled Time-Dependent Terms

# Fastest method (requires Cython)
H = [num(N), [destroy(N) + create(N), 'A * exp(-t/tau) * sin(w*t)']]
args = {'A': 1.0, 'w': 1.0, 'tau': 5.0}

# QuTiP compiles this for speed
result = sesolve(H, psi0, tlist, args=args)

Callback Functions

# Advanced control
def time_dependent_coeff(t, args):
    # Access solver state if needed
    return complex_function(t, args)

H = [H0, [H1, time_dependent_coeff]]

Parallel Processing

Parallel Map

from qutip import parallel_map

# Define task
def simulate(gamma):
    c_ops = [np.sqrt(gamma) * destroy(N)]
    result = mesolve(H, psi0, tlist, c_ops, e_ops=[num(N)])
    return result.expect[0]

# Run in parallel
gamma_values = np.linspace(0, 1, 20)
results = parallel_map(simulate, gamma_values, num_cpus=4)

Serial Map (for debugging)

from qutip import serial_map

# Same interface but runs serially
results = serial_map(simulate, gamma_values)

File I/O

Save/Load Quantum Objects

# Save
H.save('hamiltonian.qu')
psi.save('state.qu')

# Load
H_loaded = qload('hamiltonian.qu')
psi_loaded = qload('state.qu')

Save/Load Results

# Save simulation results
result = mesolve(H, psi0, tlist, c_ops, e_ops=[num(N)])
result.save('simulation.dat')

# Load results
from qutip import Result
loaded_result = Result.load('simulation.dat')

Export to MATLAB

# Export to .mat file
H.matlab_export('hamiltonian.mat', 'H')

Solver Options

Fine-Tuning Solvers

options = Options()

# Integration parameters
options.nsteps = 10000  # Max internal steps
options.rtol = 1e-8     # Relative tolerance
options.atol = 1e-10    # Absolute tolerance

# Method selection
options.method = 'adams'  # Non-stiff (default)
# options.method = 'bdf'  # Stiff problems

# Storage options
options.store_states = True
options.store_final_state = True

# Progress
options.progress_bar = True

# Random number seed (for reproducibility)
options.seeds = 12345

result = mesolve(H, psi0, tlist, c_ops, options=options)

Debugging

# Enable detailed output
options.verbose = True

# Memory tracking
options.num_cpus = 1  # Easier debugging

Performance Tips

1. Use sparse matrices: QuTiP does this automatically
2. Minimize Hilbert space: Truncate when possible
3. Choose right solver:
   • Pure states: sesolve faster than mesolve
   • Stochastic: mcsolve for quantum jumps
   • Periodic: Floquet methods
4. Time-dependent terms: String format fastest
5. Expectation values: Only compute needed observables
6. Parallel trajectories: mcsolve uses all CPUs
7. Krylov methods: For very large systems
8. Memory: Use store_final_state instead of store_states when possible