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

QuTiP Analysis and Measurement

Expectation Values

Basic Expectation Values

from qutip import *
import numpy as np

# Single operator
psi = coherent(N, 2)
n_avg = expect(num(N), psi)

# Multiple operators
ops = [num(N), destroy(N), create(N)]
results = expect(ops, psi)  # Returns list

Expectation Values for Density Matrices

# Works with both pure states and density matrices
rho = thermal_dm(N, 2)
n_avg = expect(num(N), rho)

Variance

# Calculate variance of observable
var_n = variance(num(N), psi)

# Manual calculation
var_n = expect(num(N)**2, psi) - expect(num(N), psi)**2

Time-Dependent Expectation Values

# During time evolution
result = mesolve(H, psi0, tlist, c_ops, e_ops=[num(N)])
n_t = result.expect[0]  # Array of ⟨n⟩ at each time

Entropy Measures

Von Neumann Entropy

from qutip import entropy_vn

# Density matrix entropy
rho = thermal_dm(N, 2)
S = entropy_vn(rho)  # Returns S = -Tr(ρ log₂ ρ)

Linear Entropy

from qutip import entropy_linear

# Linear entropy S_L = 1 - Tr(ρ²)
S_L = entropy_linear(rho)

Entanglement Entropy

# For bipartite systems
psi = bell_state('00')
rho = psi.proj()

# Trace out subsystem B to get reduced density matrix
rho_A = ptrace(rho, 0)

# Entanglement entropy
S_ent = entropy_vn(rho_A)

Mutual Information

from qutip import entropy_mutual

# For bipartite state ρ_AB
I = entropy_mutual(rho, [0, 1])  # I(A:B) = S(A) + S(B) - S(AB)

Conditional Entropy

from qutip import entropy_conditional

# S(A|B) = S(AB) - S(B)
S_cond = entropy_conditional(rho, 0)  # Entropy of subsystem 0 given subsystem 1

Fidelity and Distance Measures

State Fidelity

from qutip import fidelity

# Fidelity between two states
psi1 = coherent(N, 2)
psi2 = coherent(N, 2.1)

F = fidelity(psi1, psi2)  # Returns value in [0, 1]

Process Fidelity

from qutip import process_fidelity

# Fidelity between two processes (superoperators)
U_ideal = (-1j * H * t).expm()
U_actual = mesolve(H, basis(N, 0), [0, t], c_ops).states[-1]

F_proc = process_fidelity(U_ideal, U_actual)

Trace Distance

from qutip import tracedist

# Trace distance D = (1/2) Tr|ρ₁ - ρ₂|
rho1 = coherent_dm(N, 2)
rho2 = thermal_dm(N, 2)

D = tracedist(rho1, rho2)  # Returns value in [0, 1]

Hilbert-Schmidt Distance

from qutip import hilbert_dist

# Hilbert-Schmidt distance
D_HS = hilbert_dist(rho1, rho2)

Bures Distance

from qutip import bures_dist

# Bures distance
D_B = bures_dist(rho1, rho2)

Bures Angle

from qutip import bures_angle

# Bures angle
angle = bures_angle(rho1, rho2)

Entanglement Measures

Concurrence

from qutip import concurrence

# For two-qubit states
psi = bell_state('00')
rho = psi.proj()

C = concurrence(rho)  # C = 1 for maximally entangled states

Negativity

from qutip import negativity

# Negativity (partial transpose criterion)
N_ent = negativity(rho, 0)  # Partial transpose w.r.t. subsystem 0

# Logarithmic negativity
from qutip import logarithmic_negativity
E_N = logarithmic_negativity(rho, 0)

Entangling Power

from qutip import entangling_power

# For unitary gates
U = cnot()
E_pow = entangling_power(U)

Purity Measures

Purity

# Purity P = Tr(ρ²)
P = (rho * rho).tr()

# For pure states: P = 1
# For maximally mixed: P = 1/d

Checking State Properties

# Is state pure?
is_pure = abs((rho * rho).tr() - 1.0) < 1e-10

# Is operator Hermitian?
H.isherm

# Is operator unitary?
U.check_isunitary()

Measurement

Projective Measurement

from qutip import measurement

# Measure in computational basis
psi = (basis(2, 0) + basis(2, 1)).unit()

# Perform measurement
result, state_after = measurement.measure(psi, None)  # Random outcome

# Specific measurement operator
M = basis(2, 0).proj()
prob = measurement.measure_povm(psi, [M, qeye(2) - M])

Measurement Statistics

from qutip import measurement_statistics

# Get all possible outcomes and probabilities
outcomes, probabilities = measurement_statistics(psi, [M0, M1])

Observable Measurement

from qutip import measure_observable

# Measure observable and get result + collapsed state
result, state_collapsed = measure_observable(psi, sigmaz())

POVM Measurements

from qutip import measure_povm

# Positive Operator-Valued Measure
E_0 = Qobj([[0.8, 0], [0, 0.2]])
E_1 = Qobj([[0.2, 0], [0, 0.8]])

result, state_after = measure_povm(psi, [E_0, E_1])

Coherence Measures

l1-norm Coherence

from qutip import coherence_l1norm

# l1-norm of off-diagonal elements
C_l1 = coherence_l1norm(rho)

Correlation Functions

Two-Time Correlation

from qutip import correlation_2op_1t, correlation_2op_2t

# Single-time correlation ⟨A(t+τ)B(t)⟩
A = destroy(N)
B = create(N)
taulist = np.linspace(0, 10, 200)

corr = correlation_2op_1t(H, rho0, taulist, c_ops, A, B)

# Two-time correlation ⟨A(t)B(τ)⟩
tlist = np.linspace(0, 10, 100)
corr_2t = correlation_2op_2t(H, rho0, tlist, taulist, c_ops, A, B)

Three-Operator Correlation

from qutip import correlation_3op_1t

# ⟨A(t)B(t+τ)C(t)⟩
C_op = num(N)
corr_3 = correlation_3op_1t(H, rho0, taulist, c_ops, A, B, C_op)

Four-Operator Correlation

from qutip import correlation_4op_1t

# ⟨A(0)B(τ)C(τ)D(0)⟩
D_op = create(N)
corr_4 = correlation_4op_1t(H, rho0, taulist, c_ops, A, B, C_op, D_op)

Spectrum Analysis

FFT Spectrum

from qutip import spectrum_correlation_fft

# Power spectrum from correlation function
w, S = spectrum_correlation_fft(taulist, corr)

Direct Spectrum Calculation

from qutip import spectrum

# Emission/absorption spectrum
wlist = np.linspace(0, 2, 200)
spec = spectrum(H, wlist, c_ops, A, B)

Pseudo-Modes

from qutip import spectrum_pi

# Spectrum with pseudo-mode decomposition
spec_pi = spectrum_pi(H, rho0, wlist, c_ops, A, B)

Steady State Analysis

Finding Steady State

from qutip import steadystate

# Find steady state ∂ρ/∂t = 0
rho_ss = steadystate(H, c_ops)

# Different methods
rho_ss = steadystate(H, c_ops, method='direct')  # Default
rho_ss = steadystate(H, c_ops, method='eigen')   # Eigenvalue
rho_ss = steadystate(H, c_ops, method='svd')     # SVD
rho_ss = steadystate(H, c_ops, method='power')   # Power method

Steady State Properties

# Verify it's steady
L = liouvillian(H, c_ops)
assert (L * operator_to_vector(rho_ss)).norm() < 1e-10

# Compute steady-state expectation values
n_ss = expect(num(N), rho_ss)

Quantum Fisher Information

from qutip import qfisher

# Quantum Fisher information
F_Q = qfisher(rho, num(N))  # w.r.t. generator num(N)

Matrix Analysis

Eigenanalysis

# Eigenvalues and eigenvectors
evals, ekets = H.eigenstates()

# Just eigenvalues
evals = H.eigenenergies()

# Ground state
E0, psi0 = H.groundstate()

Matrix Functions

# Matrix exponential
U = (H * t).expm()

# Matrix logarithm
log_rho = rho.logm()

# Matrix square root
sqrt_rho = rho.sqrtm()

# Matrix power
rho_squared = rho ** 2

Singular Value Decomposition

# SVD of operator
U, S, Vh = H.svd()

Permutations

from qutip import permute

# Permute subsystems
rho_permuted = permute(rho, [1, 0])  # Swap subsystems

Partial Operations

Partial Trace

# Reduce to subsystem
rho_A = ptrace(rho_AB, 0)  # Keep subsystem 0
rho_B = ptrace(rho_AB, 1)  # Keep subsystem 1

# Keep multiple subsystems
rho_AC = ptrace(rho_ABC, [0, 2])  # Keep 0 and 2, trace out 1

Partial Transpose

from qutip import partial_transpose

# Partial transpose (for entanglement detection)
rho_pt = partial_transpose(rho, [0, 1])  # Transpose subsystem 0

# Check if entangled (PPT criterion)
evals = rho_pt.eigenenergies()
is_entangled = any(evals < -1e-10)

Quantum State Tomography

State Reconstruction

from qutip_qip.tomography import state_tomography

# Prepare measurement results
# measurements = ... (experimental data)

# Reconstruct density matrix
rho_reconstructed = state_tomography(measurements, basis='Pauli')

Process Tomography

from qutip_qip.tomography import qpt

# Characterize quantum process
chi = qpt(U_gate, method='lstsq')  # Chi matrix representation

Random Quantum Objects

Useful for testing and Monte Carlo simulations.

# Random state vector
psi_rand = rand_ket(N)

# Random density matrix
rho_rand = rand_dm(N)

# Random Hermitian operator
H_rand = rand_herm(N)

# Random unitary
U_rand = rand_unitary(N)

# With specific properties
rho_rank2 = rand_dm(N, rank=2)  # Rank-2 density matrix
H_sparse = rand_herm(N, density=0.1)  # 10% non-zero elements

Useful Checks

# Check if operator is Hermitian
H.isherm

# Check if state is normalized
abs(psi.norm() - 1.0) < 1e-10

# Check if density matrix is physical
rho.tr() ≈ 1 and all(rho.eigenenergies() >= 0)

# Check if operators commute
commutator(A, B).norm() < 1e-10