Superoperators and Open Systems

def liouvillian(H: Qobj, c_ops: list = None) -> Qobj:
    """
    Create Liouvillian superoperator for master equation evolution.
    
    Parameters:
    - H: System Hamiltonian
    - c_ops: List of collapse operators
    
    Returns:
    - Qobj: Liouvillian superoperator L
    """

def lindblad_dissipator(a: Qobj, b: Qobj = None) -> Qobj:
    """
    Create Lindblad dissipator superoperator.
    
    Parameters:
    - a: Collapse operator
    - b: Second operator (defaults to a.dag())
    
    Returns:
    - Qobj: Lindblad dissipator D[a,b]ρ = aρb† - ½{a†b,ρ}
    """

def spre(A: Qobj) -> Qobj:
    """
    Create superoperator for left multiplication A·ρ.
    
    Parameters:
    - A: Operator for left multiplication
    
    Returns:
    - Qobj: Superoperator representing A·ρ
    """

def spost(A: Qobj) -> Qobj:
    """
    Create superoperator for right multiplication ρ·A.
    
    Parameters:
    - A: Operator for right multiplication
    
    Returns:
    - Qobj: Superoperator representing ρ·A
    """

def sprepost(A: Qobj, B: Qobj) -> Qobj:
    """
    Create superoperator for A·ρ·B.
    
    Parameters:
    - A: Left multiplication operator
    - B: Right multiplication operator
    
    Returns:
    - Qobj: Superoperator representing A·ρ·B
    """

def operator_to_vector(op: Qobj) -> Qobj:
    """
    Convert operator to vector representation (vectorization).
    
    Parameters:
    - op: Operator to vectorize
    
    Returns:
    - Qobj: Vector representation of operator
    """

def vector_to_operator(vec: Qobj) -> Qobj:
    """
    Convert vector back to operator representation.
    
    Parameters:
    - vec: Vector representation
    
    Returns:
    - Qobj: Operator reconstructed from vector
    """

def stack_columns(op: Qobj) -> Qobj:
    """
    Stack columns of operator matrix into vector.
    
    Parameters:
    - op: Operator to stack
    
    Returns:
    - Qobj: Column-stacked vector
    """

def unstack_columns(vec: Qobj, shape: tuple) -> Qobj:
    """
    Unstack vector back into operator matrix.
    
    Parameters:
    - vec: Stacked vector
    - shape: Target operator shape
    
    Returns:
    - Qobj: Reconstructed operator
    """

def reshuffle(obj: Qobj) -> Qobj:
    """
    Reshuffle superoperator indices (Choi to chi representation).
    
    Parameters:
    - obj: Superoperator to reshuffle
    
    Returns:
    - Qobj: Reshuffled superoperator
    """

def to_choi(kraus_ops: list) -> Qobj:
    """
    Convert Kraus operators to Choi matrix representation.
    
    Parameters:
    - kraus_ops: List of Kraus operators
    
    Returns:
    - Qobj: Choi matrix representation
    """

def to_kraus(choi: Qobj) -> list:
    """
    Convert Choi matrix to Kraus operator representation.
    
    Parameters:
    - choi: Choi matrix
    
    Returns:
    - list: List of Kraus operators
    """

def to_super(op: Qobj) -> Qobj:
    """
    Convert operator to superoperator representation.
    
    Parameters:
    - op: Operator to convert
    
    Returns:
    - Qobj: Superoperator representation
    """

import qutip as qt
import numpy as np

# Build Liouvillian for damped harmonic oscillator
N = 10
a = qt.destroy(N)
H = a.dag() * a  # Hamiltonian

# Collapse operators
gamma = 0.1
c_ops = [np.sqrt(gamma) * a]  # Damping

# Create Liouvillian
L = qt.liouvillian(H, c_ops)
print(f"Liouvillian shape: {L.shape}")

# Manual Liouvillian construction
L_manual = -1j * (qt.spre(H) - qt.spost(H))  # Hamiltonian part
for c in c_ops:
    L_manual += qt.lindblad_dissipator(c)  # Add dissipation

# Verify they're the same
diff = (L - L_manual).norm()
print(f"Manual vs automatic Liouvillian difference: {diff:.2e}")

# Superoperator operations
sigma_x = qt.sigmax()
sigma_z = qt.sigmaz()

# Left and right multiplication superoperators
S_left = qt.spre(sigma_x)   # σₓ·ρ
S_right = qt.spost(sigma_z) # ρ·σᵤ
S_both = qt.sprepost(sigma_x, sigma_z)  # σₓ·ρ·σᵤ

# Test on a density matrix
rho = qt.ket2dm(qt.basis(2, 0))  # |0⟩⟨0|

# Apply superoperators
rho_left = S_left * qt.operator_to_vector(rho)
rho_left = qt.vector_to_operator(rho_left)
print(f"σₓ|0⟩⟨0|: {rho_left}")

# Vector/operator conversion
rho_vec = qt.operator_to_vector(rho)
print(f"Vectorized density matrix shape: {rho_vec.shape}")

rho_reconstructed = qt.vector_to_operator(rho_vec)
print(f"Reconstruction fidelity: {qt.fidelity(rho, rho_reconstructed):.6f}")

# Column stacking (alternative vectorization)
rho_stacked = qt.stack_columns(rho)
rho_unstacked = qt.unstack_columns(rho_stacked, rho.shape)
print(f"Stack/unstack fidelity: {qt.fidelity(rho, rho_unstacked):.6f}")

# Example: Amplitude damping channel
def amplitude_damping_kraus(gamma):
    """Create Kraus operators for amplitude damping."""
    K0 = qt.Qobj([[1, 0], [0, np.sqrt(1-gamma)]])
    K1 = qt.Qobj([[0, np.sqrt(gamma)], [0, 0]])
    return [K0, K1]

# Convert to different representations
gamma = 0.3
kraus_ops = amplitude_damping_kraus(gamma)

# Kraus to Choi matrix
choi = qt.to_choi(kraus_ops)
print(f"Choi matrix:\n{choi}")

# Choi back to Kraus
kraus_reconstructed = qt.to_kraus(choi)
print(f"Number of reconstructed Kraus ops: {len(kraus_reconstructed)}")

# Apply channel using different representations
psi = qt.basis(2, 1)  # |1⟩ state
rho_in = psi * psi.dag()

# Using Kraus operators
rho_out_kraus = sum(K * rho_in * K.dag() for K in kraus_ops)

# Using Choi matrix (more complex, via vectorization)
rho_vec = qt.operator_to_vector(rho_in)
# For Choi representation: need to apply proper transformation
# This is a simplified example - full Choi application is more involved

print(f"Output state after amplitude damping: {rho_out_kraus}")
print(f"Excited state population: {qt.expect(qt.num(2), rho_out_kraus):.3f}")

# Superoperator for general channel
def depolarizing_channel(p):
    """Create depolarizing channel superoperator."""
    I = qt.qeye(2)
    X = qt.sigmax()
    Y = qt.sigmay() 
    Z = qt.sigmaz()
    
    # Kraus operators
    K0 = np.sqrt(1 - 3*p/4) * I
    K1 = np.sqrt(p/4) * X
    K2 = np.sqrt(p/4) * Y
    K3 = np.sqrt(p/4) * Z
    
    # Build superoperator
    S = (1 - 3*p/4) * qt.sprepost(I, I)
    S += (p/4) * (qt.sprepost(X, X) + qt.sprepost(Y, Y) + qt.sprepost(Z, Z))
    
    return S

# Test depolarizing channel
p = 0.1
S_depol = depolarizing_channel(p)

# Apply to pure state
psi_pure = (qt.basis(2,0) + qt.basis(2,1)).unit()
rho_pure = psi_pure * psi_pure.dag()
rho_vec = qt.operator_to_vector(rho_pure)
rho_depol_vec = S_depol * rho_vec
rho_depol = qt.vector_to_operator(rho_depol_vec)

print(f"Initial purity: {(rho_pure**2).tr():.3f}")
print(f"After depolarizing: {(rho_depol**2).tr():.3f}")

# Superoperator functions return Qobj instances with type='super'
# representing linear transformations on the space of operators