Custom Hamiltonians¶

This notebook covers advanced usage with custom Hamiltonians in QubitOS. You'll learn how to build complex quantum systems, define custom control operators, and optimize multi-qubit gates.

Prerequisites¶

• Familiarity with the GRAPE optimization notebook
• Understanding of Pauli operators and tensor products
• Basic quantum mechanics (Hamiltonians, unitaries)

Topics Covered¶

1. Pauli String Notation - Compact Hamiltonian specification
2. Building Hamiltonians Programmatically - Using NumPy arrays
3. Multi-Qubit Systems - Tensor products and embedding
4. Custom Control Operators - Beyond X and Y drives
5. Drift Hamiltonians - Modeling real qubit physics
6. Advanced Examples - Transmon anharmonicity, ZZ coupling

In [ ]:

# Import required modules
import numpy as np
import matplotlib.pyplot as plt
from qubitos.pulsegen.grape import GrapeOptimizer, GrapeConfig, generate_pulse
from qubitos.pulsegen.hamiltonians import (
    # Pauli matrices
    PAULI_I, PAULI_X, PAULI_Y, PAULI_Z,
    PAULI_MATRICES,
    # Standard gates
    STANDARD_GATES,
    # Functions
    parse_pauli_string,
    pauli_string_to_matrix,
    build_hamiltonian,
    tensor_product,
    get_target_unitary,
    rotation_gate,
    embed_gate,
)

# Plotting setup
plt.style.use('seaborn-v0_8-whitegrid')
plt.rcParams['figure.figsize'] = (10, 6)
plt.rcParams['font.size'] = 12

# Helper function for displaying matrices
def show_matrix(M, title="Matrix", precision=3):
    """Display a matrix with nice formatting."""
    print(f"{title}:")
    print(np.round(M, precision))
    print()

# Import required modules import numpy as np import matplotlib.pyplot as plt from qubitos.pulsegen.grape import GrapeOptimizer, GrapeConfig, generate_pulse from qubitos.pulsegen.hamiltonians import ( # Pauli matrices PAULI_I, PAULI_X, PAULI_Y, PAULI_Z, PAULI_MATRICES, # Standard gates STANDARD_GATES, # Functions parse_pauli_string, pauli_string_to_matrix, build_hamiltonian, tensor_product, get_target_unitary, rotation_gate, embed_gate, ) # Plotting setup plt.style.use('seaborn-v0_8-whitegrid') plt.rcParams['figure.figsize'] = (10, 6) plt.rcParams['font.size'] = 12 # Helper function for displaying matrices def show_matrix(M, title="Matrix", precision=3): """Display a matrix with nice formatting.""" print(f"{title}:") print(np.round(M, precision)) print()

---

Part 1: Pauli Matrices Review¶

The Pauli matrices are the fundamental building blocks for qubit Hamiltonians.

In [ ]:

# Display the Pauli matrices
print("=" * 50)
print("PAULI MATRICES")
print("=" * 50)

for name, matrix in PAULI_MATRICES.items():
    show_matrix(matrix, f"sigma_{name}")

# Display the Pauli matrices print("=" * 50) print("PAULI MATRICES") print("=" * 50) for name, matrix in PAULI_MATRICES.items(): show_matrix(matrix, f"sigma_{name}")

In [ ]:

# Key properties of Pauli matrices
print("Key Properties:")
print()

# Hermitian
print("1. Hermitian (self-adjoint): X^dag = X")
print(f"   X^dag == X: {np.allclose(PAULI_X.conj().T, PAULI_X)}")
print()

# Unitary
print("2. Unitary: X @ X^dag = I")
print(f"   X @ X^dag == I: {np.allclose(PAULI_X @ PAULI_X.conj().T, PAULI_I)}")
print()

# Square to identity
print("3. Square to identity: X^2 = I")
print(f"   X^2 == I: {np.allclose(PAULI_X @ PAULI_X, PAULI_I)}")
print()

# Anticommutation
print("4. Anticommutation: XY + YX = 0")
anticomm = PAULI_X @ PAULI_Y + PAULI_Y @ PAULI_X
print(f"   XY + YX = 0: {np.allclose(anticomm, 0)}")

# Key properties of Pauli matrices print("Key Properties:") print() # Hermitian print("1. Hermitian (self-adjoint): X^dag = X") print(f" X^dag == X: {np.allclose(PAULI_X.conj().T, PAULI_X)}") print() # Unitary print("2. Unitary: X @ X^dag = I") print(f" X @ X^dag == I: {np.allclose(PAULI_X @ PAULI_X.conj().T, PAULI_I)}") print() # Square to identity print("3. Square to identity: X^2 = I") print(f" X^2 == I: {np.allclose(PAULI_X @ PAULI_X, PAULI_I)}") print() # Anticommutation print("4. Anticommutation: XY + YX = 0") anticomm = PAULI_X @ PAULI_Y + PAULI_Y @ PAULI_X print(f" XY + YX = 0: {np.allclose(anticomm, 0)}")

---

Part 2: Pauli String Notation¶

QubitOS provides a convenient string-based notation for specifying Hamiltonians.

Format¶

"coefficient * Operator0 Operator1 ... + coefficient * ..."

Where:

• coefficient is a real number (optional, defaults to 1.0)
• OperatorN is a Pauli operator (I, X, Y, Z) followed by qubit index N
• Terms are separated by + or -

In [ ]:

# Simple examples
print("=" * 50)
print("PAULI STRING EXAMPLES")
print("=" * 50)

# Single qubit, single term
H1 = parse_pauli_string("X0", num_qubits=1)
show_matrix(H1, "H = X0 (X on qubit 0)")

# Single qubit with coefficient
H2 = parse_pauli_string("0.5 * Z0", num_qubits=1)
show_matrix(H2, "H = 0.5 * Z0")

# Multiple terms
H3 = parse_pauli_string("1.0 * X0 + 0.5 * Z0", num_qubits=1)
show_matrix(H3, "H = X0 + 0.5 * Z0")

# Simple examples print("=" * 50) print("PAULI STRING EXAMPLES") print("=" * 50) # Single qubit, single term H1 = parse_pauli_string("X0", num_qubits=1) show_matrix(H1, "H = X0 (X on qubit 0)") # Single qubit with coefficient H2 = parse_pauli_string("0.5 * Z0", num_qubits=1) show_matrix(H2, "H = 0.5 * Z0") # Multiple terms H3 = parse_pauli_string("1.0 * X0 + 0.5 * Z0", num_qubits=1) show_matrix(H3, "H = X0 + 0.5 * Z0")

In [ ]:

# Two-qubit examples
print("=" * 50)
print("TWO-QUBIT PAULI STRINGS")
print("=" * 50)

# X on qubit 0 only (tensor with I on qubit 1)
H4 = parse_pauli_string("X0", num_qubits=2)
show_matrix(H4, "H = X0 (in 2-qubit space)")

# Compare with explicit tensor product
H4_explicit = np.kron(PAULI_X, PAULI_I)
print(f"Matches X x I: {np.allclose(H4, H4_explicit)}")
print()

# ZZ coupling term
H5 = parse_pauli_string("Z0 Z1", num_qubits=2)
show_matrix(H5, "H = Z0 Z1 (ZZ coupling)")

# Two-qubit examples print("=" * 50) print("TWO-QUBIT PAULI STRINGS") print("=" * 50) # X on qubit 0 only (tensor with I on qubit 1) H4 = parse_pauli_string("X0", num_qubits=2) show_matrix(H4, "H = X0 (in 2-qubit space)") # Compare with explicit tensor product H4_explicit = np.kron(PAULI_X, PAULI_I) print(f"Matches X x I: {np.allclose(H4, H4_explicit)}") print() # ZZ coupling term H5 = parse_pauli_string("Z0 Z1", num_qubits=2) show_matrix(H5, "H = Z0 Z1 (ZZ coupling)")

In [ ]:

# Negative coefficients
print("=" * 50)
print("NEGATIVE COEFFICIENTS")
print("=" * 50)

H6 = parse_pauli_string("1.0 * X0 - 0.5 * Y0", num_qubits=1)
show_matrix(H6, "H = X0 - 0.5 * Y0")

# Verify
expected = PAULI_X - 0.5 * PAULI_Y
print(f"Matches expected: {np.allclose(H6, expected)}")

# Negative coefficients print("=" * 50) print("NEGATIVE COEFFICIENTS") print("=" * 50) H6 = parse_pauli_string("1.0 * X0 - 0.5 * Y0", num_qubits=1) show_matrix(H6, "H = X0 - 0.5 * Y0") # Verify expected = PAULI_X - 0.5 * PAULI_Y print(f"Matches expected: {np.allclose(H6, expected)}")

---

Part 3: Building Hamiltonians Programmatically¶

For complex systems, you can build Hamiltonians using NumPy operations.

In [ ]:

# Tensor products for multi-qubit operators
print("=" * 50)
print("TENSOR PRODUCTS")
print("=" * 50)

# X on qubit 0, I on qubit 1: X x I
X0_I1 = tensor_product([PAULI_X, PAULI_I])
show_matrix(X0_I1, "X x I")

# I on qubit 0, X on qubit 1: I x X
I0_X1 = tensor_product([PAULI_I, PAULI_X])
show_matrix(I0_X1, "I x X")

# ZZ coupling: Z x Z
ZZ = tensor_product([PAULI_Z, PAULI_Z])
show_matrix(ZZ, "Z x Z")

# Tensor products for multi-qubit operators print("=" * 50) print("TENSOR PRODUCTS") print("=" * 50) # X on qubit 0, I on qubit 1: X x I X0_I1 = tensor_product([PAULI_X, PAULI_I]) show_matrix(X0_I1, "X x I") # I on qubit 0, X on qubit 1: I x X I0_X1 = tensor_product([PAULI_I, PAULI_X]) show_matrix(I0_X1, "I x X") # ZZ coupling: Z x Z ZZ = tensor_product([PAULI_Z, PAULI_Z]) show_matrix(ZZ, "Z x Z")

In [ ]:

# Building a transmon-like Hamiltonian
print("=" * 50)
print("TRANSMON-LIKE HAMILTONIAN")
print("=" * 50)

# Parameters (in MHz)
omega_q = 5000.0  # Qubit frequency
delta = 200.0     # Anharmonicity (negative for transmon)

# For a single qubit (2-level approximation)
# H_0 = (omega_q / 2) * Z
H_drift = (omega_q / 2) * PAULI_Z

print(f"Qubit frequency: {omega_q} MHz")
print(f"Drift Hamiltonian: (omega_q / 2) * Z")
print(f"Eigenvalues: {np.linalg.eigvalsh(H_drift)}")
print(f"Energy gap: {np.linalg.eigvalsh(H_drift)[1] - np.linalg.eigvalsh(H_drift)[0]} MHz")

# Building a transmon-like Hamiltonian print("=" * 50) print("TRANSMON-LIKE HAMILTONIAN") print("=" * 50) # Parameters (in MHz) omega_q = 5000.0 # Qubit frequency delta = 200.0 # Anharmonicity (negative for transmon) # For a single qubit (2-level approximation) # H_0 = (omega_q / 2) * Z H_drift = (omega_q / 2) * PAULI_Z print(f"Qubit frequency: {omega_q} MHz") print(f"Drift Hamiltonian: (omega_q / 2) * Z") print(f"Eigenvalues: {np.linalg.eigvalsh(H_drift)}") print(f"Energy gap: {np.linalg.eigvalsh(H_drift)[1] - np.linalg.eigvalsh(H_drift)[0]} MHz")

In [ ]:

# Using build_hamiltonian helper
print("=" * 50)
print("BUILD_HAMILTONIAN HELPER")
print("=" * 50)

# Build drift and control Hamiltonians
drift, controls = build_hamiltonian(
    drift="2500 * Z0",  # omega_q/2 * Z
    controls=["X0", "Y0"],  # X and Y drives
    num_qubits=1,
)

show_matrix(drift, "Drift Hamiltonian")
show_matrix(controls[0], "Control 1 (X drive)")
show_matrix(controls[1], "Control 2 (Y drive)")

# Using build_hamiltonian helper print("=" * 50) print("BUILD_HAMILTONIAN HELPER") print("=" * 50) # Build drift and control Hamiltonians drift, controls = build_hamiltonian( drift="2500 * Z0", # omega_q/2 * Z controls=["X0", "Y0"], # X and Y drives num_qubits=1, ) show_matrix(drift, "Drift Hamiltonian") show_matrix(controls[0], "Control 1 (X drive)") show_matrix(controls[1], "Control 2 (Y drive)")

---

Part 4: Multi-Qubit Systems¶

Building Hamiltonians for multiple qubits requires careful handling of tensor products.

In [ ]:

# Two-qubit system with ZZ coupling
print("=" * 50)
print("TWO-QUBIT SYSTEM WITH ZZ COUPLING")
print("=" * 50)

# Parameters
omega_0 = 5000.0  # Qubit 0 frequency (MHz)
omega_1 = 5200.0  # Qubit 1 frequency (MHz)
J = 5.0           # ZZ coupling strength (MHz)

# Build drift Hamiltonian
# H_drift = (omega_0/2) * Z_0 + (omega_1/2) * Z_1 + J * Z_0 Z_1
drift_str = f"{omega_0/2} * Z0 + {omega_1/2} * Z1 + {J} * Z0 Z1"
H_drift_2q = parse_pauli_string(drift_str, num_qubits=2)

print(f"Drift Hamiltonian:")
print(f"  H = {omega_0/2:.0f} * Z0 + {omega_1/2:.0f} * Z1 + {J:.0f} * Z0 Z1")
print()
show_matrix(H_drift_2q, "H_drift")

# Eigenvalues
eigvals = np.linalg.eigvalsh(H_drift_2q)
print("Eigenvalues (MHz):")
for i, e in enumerate(eigvals):
    print(f"  |{i:02b}>: {e:.1f}")

# Two-qubit system with ZZ coupling print("=" * 50) print("TWO-QUBIT SYSTEM WITH ZZ COUPLING") print("=" * 50) # Parameters omega_0 = 5000.0 # Qubit 0 frequency (MHz) omega_1 = 5200.0 # Qubit 1 frequency (MHz) J = 5.0 # ZZ coupling strength (MHz) # Build drift Hamiltonian # H_drift = (omega_0/2) * Z_0 + (omega_1/2) * Z_1 + J * Z_0 Z_1 drift_str = f"{omega_0/2} * Z0 + {omega_1/2} * Z1 + {J} * Z0 Z1" H_drift_2q = parse_pauli_string(drift_str, num_qubits=2) print(f"Drift Hamiltonian:") print(f" H = {omega_0/2:.0f} * Z0 + {omega_1/2:.0f} * Z1 + {J:.0f} * Z0 Z1") print() show_matrix(H_drift_2q, "H_drift") # Eigenvalues eigvals = np.linalg.eigvalsh(H_drift_2q) print("Eigenvalues (MHz):") for i, e in enumerate(eigvals): print(f" |{i:02b}>: {e:.1f}")

In [ ]:

# Control Hamiltonians for two-qubit system
print("=" * 50)
print("TWO-QUBIT CONTROL HAMILTONIANS")
print("=" * 50)

# Typically: X and Y drives on each qubit
_, controls_2q = build_hamiltonian(
    drift=None,
    controls=None,  # Default: X and Y on each qubit
    num_qubits=2,
)

print(f"Number of control Hamiltonians: {len(controls_2q)}")
for i, H_c in enumerate(controls_2q):
    qubit = i // 2
    drive = "X" if i % 2 == 0 else "Y"
    print(f"Control {i}: {drive} on qubit {qubit}")

# Control Hamiltonians for two-qubit system print("=" * 50) print("TWO-QUBIT CONTROL HAMILTONIANS") print("=" * 50) # Typically: X and Y drives on each qubit _, controls_2q = build_hamiltonian( drift=None, controls=None, # Default: X and Y on each qubit num_qubits=2, ) print(f"Number of control Hamiltonians: {len(controls_2q)}") for i, H_c in enumerate(controls_2q): qubit = i // 2 drive = "X" if i % 2 == 0 else "Y" print(f"Control {i}: {drive} on qubit {qubit}")

In [ ]:

# Embedding single-qubit gates in multi-qubit space
print("=" * 50)
print("EMBEDDING GATES")
print("=" * 50)

# X gate on qubit 0 in a 2-qubit system
X_on_q0 = embed_gate(STANDARD_GATES['X'], num_qubits=2, qubit_indices=[0])
show_matrix(X_on_q0, "X gate on qubit 0 (2-qubit space)")

# X gate on qubit 1 in a 2-qubit system
X_on_q1 = embed_gate(STANDARD_GATES['X'], num_qubits=2, qubit_indices=[1])
show_matrix(X_on_q1, "X gate on qubit 1 (2-qubit space)")

# Verify: X0 @ X1 should flip both qubits
XX = X_on_q0 @ X_on_q1
print("X0 @ X1 applied to |00>:")
state_00 = np.array([1, 0, 0, 0])
result = XX @ state_00
print(f"  |00> -> {result}  (should be |11>)")

# Embedding single-qubit gates in multi-qubit space print("=" * 50) print("EMBEDDING GATES") print("=" * 50) # X gate on qubit 0 in a 2-qubit system X_on_q0 = embed_gate(STANDARD_GATES['X'], num_qubits=2, qubit_indices=[0]) show_matrix(X_on_q0, "X gate on qubit 0 (2-qubit space)") # X gate on qubit 1 in a 2-qubit system X_on_q1 = embed_gate(STANDARD_GATES['X'], num_qubits=2, qubit_indices=[1]) show_matrix(X_on_q1, "X gate on qubit 1 (2-qubit space)") # Verify: X0 @ X1 should flip both qubits XX = X_on_q0 @ X_on_q1 print("X0 @ X1 applied to |00>:") state_00 = np.array([1, 0, 0, 0]) result = XX @ state_00 print(f" |00> -> {result} (should be |11>)")

---

Part 5: Custom Unitaries¶

You can optimize for arbitrary target unitaries, not just standard gates.

In [ ]:

# Rotation gates
print("=" * 50)
print("ROTATION GATES")
print("=" * 50)

angles = [np.pi/4, np.pi/2, np.pi, 3*np.pi/2]

for angle in angles:
    Rx = rotation_gate("X", angle)
    print(f"R_X({angle/np.pi:.2f}*pi):")
    print(np.round(Rx, 3))
    print()

# Rotation gates print("=" * 50) print("ROTATION GATES") print("=" * 50) angles = [np.pi/4, np.pi/2, np.pi, 3*np.pi/2] for angle in angles: Rx = rotation_gate("X", angle) print(f"R_X({angle/np.pi:.2f}*pi):") print(np.round(Rx, 3)) print()

In [ ]:

# Optimizing for a rotation gate
print("=" * 50)
print("OPTIMIZING R_X(pi/4)")
print("=" * 50)

# Target: R_X(pi/4)
target_rx = get_target_unitary("RX", num_qubits=1, angle=np.pi/4)
show_matrix(target_rx, "Target: R_X(pi/4)")

# Optimize
config = GrapeConfig(
    num_time_steps=100,
    duration_ns=20.0,
    target_fidelity=0.999,
    max_iterations=500,
    random_seed=42,
)

optimizer = GrapeOptimizer(config)
result = optimizer.optimize(target_rx, num_qubits=1)

print(f"Optimization result:")
print(f"  Fidelity: {result.fidelity:.6f}")
print(f"  Iterations: {result.iterations}")

show_matrix(result.final_unitary, "Achieved unitary")

# Optimizing for a rotation gate print("=" * 50) print("OPTIMIZING R_X(pi/4)") print("=" * 50) # Target: R_X(pi/4) target_rx = get_target_unitary("RX", num_qubits=1, angle=np.pi/4) show_matrix(target_rx, "Target: R_X(pi/4)") # Optimize config = GrapeConfig( num_time_steps=100, duration_ns=20.0, target_fidelity=0.999, max_iterations=500, random_seed=42, ) optimizer = GrapeOptimizer(config) result = optimizer.optimize(target_rx, num_qubits=1) print(f"Optimization result:") print(f" Fidelity: {result.fidelity:.6f}") print(f" Iterations: {result.iterations}") show_matrix(result.final_unitary, "Achieved unitary")

In [ ]:

# Custom arbitrary unitary
print("=" * 50)
print("CUSTOM UNITARY")
print("=" * 50)

# Create a random unitary (for demonstration)
# In practice, you'd define your specific target
def random_unitary(dim, seed=42):
    """Generate a random unitary matrix."""
    rng = np.random.default_rng(seed)
    # Random complex matrix
    A = rng.normal(size=(dim, dim)) + 1j * rng.normal(size=(dim, dim))
    # QR decomposition gives orthogonal matrix
    Q, R = np.linalg.qr(A)
    # Ensure determinant is 1 (SU(n))
    Q = Q / np.linalg.det(Q)**(1/dim)
    return Q

custom_target = random_unitary(2, seed=123)
show_matrix(custom_target, "Custom target unitary")

# Verify it's unitary
print(f"Is unitary: {np.allclose(custom_target @ custom_target.conj().T, np.eye(2))}")
print()

# Optimize for this target
result_custom = optimizer.optimize(custom_target, num_qubits=1)
print(f"Optimization result:")
print(f"  Fidelity: {result_custom.fidelity:.6f}")
print(f"  Iterations: {result_custom.iterations}")

# Custom arbitrary unitary print("=" * 50) print("CUSTOM UNITARY") print("=" * 50) # Create a random unitary (for demonstration) # In practice, you'd define your specific target def random_unitary(dim, seed=42): """Generate a random unitary matrix.""" rng = np.random.default_rng(seed) # Random complex matrix A = rng.normal(size=(dim, dim)) + 1j * rng.normal(size=(dim, dim)) # QR decomposition gives orthogonal matrix Q, R = np.linalg.qr(A) # Ensure determinant is 1 (SU(n)) Q = Q / np.linalg.det(Q)**(1/dim) return Q custom_target = random_unitary(2, seed=123) show_matrix(custom_target, "Custom target unitary") # Verify it's unitary print(f"Is unitary: {np.allclose(custom_target @ custom_target.conj().T, np.eye(2))}") print() # Optimize for this target result_custom = optimizer.optimize(custom_target, num_qubits=1) print(f"Optimization result:") print(f" Fidelity: {result_custom.fidelity:.6f}") print(f" Iterations: {result_custom.iterations}")

---

Part 6: Drift Hamiltonians in Optimization¶

Including realistic drift Hamiltonians is important for accurate pulse design.

In [ ]:

# Effect of detuning on X gate optimization
print("=" * 50)
print("EFFECT OF DETUNING")
print("=" * 50)

detunings = [0, 5, 10, 20, 50]  # MHz
results_by_detuning = {}

config = GrapeConfig(
    num_time_steps=100,
    duration_ns=20.0,
    target_fidelity=0.999,
    max_iterations=500,
    random_seed=42,
)

target_x = get_target_unitary("X")

for delta in detunings:
    drift = delta * PAULI_Z if delta > 0 else None
    optimizer = GrapeOptimizer(config)
    result = optimizer.optimize(target_x, num_qubits=1, drift_hamiltonian=drift)
    results_by_detuning[delta] = result
    print(f"Detuning = {delta:2d} MHz: fidelity = {result.fidelity:.6f}, iterations = {result.iterations}")

# Effect of detuning on X gate optimization print("=" * 50) print("EFFECT OF DETUNING") print("=" * 50) detunings = [0, 5, 10, 20, 50] # MHz results_by_detuning = {} config = GrapeConfig( num_time_steps=100, duration_ns=20.0, target_fidelity=0.999, max_iterations=500, random_seed=42, ) target_x = get_target_unitary("X") for delta in detunings: drift = delta * PAULI_Z if delta > 0 else None optimizer = GrapeOptimizer(config) result = optimizer.optimize(target_x, num_qubits=1, drift_hamiltonian=drift) results_by_detuning[delta] = result print(f"Detuning = {delta:2d} MHz: fidelity = {result.fidelity:.6f}, iterations = {result.iterations}")

In [ ]:

# Visualize pulses for different detunings
fig, axes = plt.subplots(2, 3, figsize=(15, 8))
axes = axes.flatten()

for i, (delta, result) in enumerate(results_by_detuning.items()):
    if i >= len(axes):
        break
    ax = axes[i]
    t = np.linspace(0, 20, len(result.i_envelope))
    ax.plot(t, result.i_envelope, 'b-', label='I', linewidth=1.5)
    ax.plot(t, result.q_envelope, 'r-', label='Q', linewidth=1.5)
    ax.set_xlabel('Time (ns)')
    ax.set_ylabel('Amplitude (MHz)')
    ax.set_title(f'Detuning = {delta} MHz\n(F = {result.fidelity:.4f})')
    ax.legend(loc='upper right')
    ax.axhline(y=0, color='gray', linestyle='-', alpha=0.3)

# Hide unused subplot
axes[-1].axis('off')

plt.suptitle('X Gate Pulses with Different Detunings', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()

# Visualize pulses for different detunings fig, axes = plt.subplots(2, 3, figsize=(15, 8)) axes = axes.flatten() for i, (delta, result) in enumerate(results_by_detuning.items()): if i >= len(axes): break ax = axes[i] t = np.linspace(0, 20, len(result.i_envelope)) ax.plot(t, result.i_envelope, 'b-', label='I', linewidth=1.5) ax.plot(t, result.q_envelope, 'r-', label='Q', linewidth=1.5) ax.set_xlabel('Time (ns)') ax.set_ylabel('Amplitude (MHz)') ax.set_title(f'Detuning = {delta} MHz\n(F = {result.fidelity:.4f})') ax.legend(loc='upper right') ax.axhline(y=0, color='gray', linestyle='-', alpha=0.3) # Hide unused subplot axes[-1].axis('off') plt.suptitle('X Gate Pulses with Different Detunings', fontsize=14, fontweight='bold') plt.tight_layout() plt.show()

Observation¶

Notice how the Q quadrature (red) becomes more prominent with larger detuning. This is because the optimizer compensates for the Z rotation caused by the detuning by applying an appropriate phase correction.

---

Part 7: Two-Qubit Gate Optimization¶

Two-qubit gates are more challenging but follow the same principles.

In [ ]:

# Available two-qubit gates
print("=" * 50)
print("STANDARD TWO-QUBIT GATES")
print("=" * 50)

two_qubit_gates = ['CZ', 'CNOT', 'ISWAP', 'SWAP']

for gate in two_qubit_gates:
    U = STANDARD_GATES[gate]
    show_matrix(U, gate)

# Available two-qubit gates print("=" * 50) print("STANDARD TWO-QUBIT GATES") print("=" * 50) two_qubit_gates = ['CZ', 'CNOT', 'ISWAP', 'SWAP'] for gate in two_qubit_gates: U = STANDARD_GATES[gate] show_matrix(U, gate)

In [ ]:

# Optimize CNOT gate
print("=" * 50)
print("OPTIMIZING CNOT GATE")
print("=" * 50)

config_2q = GrapeConfig(
    num_time_steps=100,
    duration_ns=50.0,  # Longer for two-qubit
    target_fidelity=0.99,
    max_iterations=1000,
    random_seed=42,
)

target_cnot = get_target_unitary("CNOT", num_qubits=2)
optimizer_2q = GrapeOptimizer(config_2q)

print("Optimizing CNOT (this may take a moment)...")
result_cnot = optimizer_2q.optimize(target_cnot, num_qubits=2)

print(f"\nResult:")
print(f"  Fidelity: {result_cnot.fidelity:.6f}")
print(f"  Iterations: {result_cnot.iterations}")
print(f"  Converged: {result_cnot.converged}")

# Optimize CNOT gate print("=" * 50) print("OPTIMIZING CNOT GATE") print("=" * 50) config_2q = GrapeConfig( num_time_steps=100, duration_ns=50.0, # Longer for two-qubit target_fidelity=0.99, max_iterations=1000, random_seed=42, ) target_cnot = get_target_unitary("CNOT", num_qubits=2) optimizer_2q = GrapeOptimizer(config_2q) print("Optimizing CNOT (this may take a moment)...") result_cnot = optimizer_2q.optimize(target_cnot, num_qubits=2) print(f"\nResult:") print(f" Fidelity: {result_cnot.fidelity:.6f}") print(f" Iterations: {result_cnot.iterations}") print(f" Converged: {result_cnot.converged}")

In [ ]:

# Plot CNOT pulse and verify
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Pulse
t = np.linspace(0, 50, len(result_cnot.i_envelope))
axes[0].plot(t, result_cnot.i_envelope, 'b-', label='I', linewidth=1.5)
axes[0].plot(t, result_cnot.q_envelope, 'r-', label='Q', linewidth=1.5)
axes[0].set_xlabel('Time (ns)')
axes[0].set_ylabel('Amplitude (MHz)')
axes[0].set_title('CNOT Gate Pulse')
axes[0].legend()

# Convergence
axes[1].semilogy(1 - np.array(result_cnot.fidelity_history), 'b-', linewidth=1.5)
axes[1].axhline(y=0.01, color='red', linestyle='--', label='Target 99%')
axes[1].set_xlabel('Iteration')
axes[1].set_ylabel('Infidelity (1 - F)')
axes[1].set_title('Convergence')
axes[1].legend()

plt.tight_layout()
plt.show()

# Plot CNOT pulse and verify fig, axes = plt.subplots(1, 2, figsize=(14, 5)) # Pulse t = np.linspace(0, 50, len(result_cnot.i_envelope)) axes[0].plot(t, result_cnot.i_envelope, 'b-', label='I', linewidth=1.5) axes[0].plot(t, result_cnot.q_envelope, 'r-', label='Q', linewidth=1.5) axes[0].set_xlabel('Time (ns)') axes[0].set_ylabel('Amplitude (MHz)') axes[0].set_title('CNOT Gate Pulse') axes[0].legend() # Convergence axes[1].semilogy(1 - np.array(result_cnot.fidelity_history), 'b-', linewidth=1.5) axes[1].axhline(y=0.01, color='red', linestyle='--', label='Target 99%') axes[1].set_xlabel('Iteration') axes[1].set_ylabel('Infidelity (1 - F)') axes[1].set_title('Convergence') axes[1].legend() plt.tight_layout() plt.show()

---

Part 8: Three-Level Systems (Transmon)¶

For higher accuracy with transmons, we can include the third level (leakage).

In [ ]:

# Three-level (qutrit) operators
print("=" * 50)
print("THREE-LEVEL SYSTEM (QUTRIT)")
print("=" * 50)

# Lowering operator
a = np.array([
    [0, 1, 0],
    [0, 0, np.sqrt(2)],
    [0, 0, 0],
], dtype=np.complex128)

# Raising operator
a_dag = a.conj().T

# Number operator
n = a_dag @ a

print("Lowering operator (a):")
print(np.round(a, 3))
print()

print("Number operator (n = a_dag @ a):")
print(np.round(n, 3))

# Three-level (qutrit) operators print("=" * 50) print("THREE-LEVEL SYSTEM (QUTRIT)") print("=" * 50) # Lowering operator a = np.array([ [0, 1, 0], [0, 0, np.sqrt(2)], [0, 0, 0], ], dtype=np.complex128) # Raising operator a_dag = a.conj().T # Number operator n = a_dag @ a print("Lowering operator (a):") print(np.round(a, 3)) print() print("Number operator (n = a_dag @ a):") print(np.round(n, 3))

In [ ]:

# Transmon Hamiltonian with anharmonicity
print("=" * 50)
print("TRANSMON HAMILTONIAN")
print("=" * 50)

# Parameters
omega_q = 5000.0  # Qubit frequency (MHz)
alpha = -200.0    # Anharmonicity (MHz), negative for transmon

# H = omega_q * n + (alpha/2) * n * (n - I)
I3 = np.eye(3, dtype=np.complex128)
H_transmon = omega_q * n + (alpha / 2) * n @ (n - I3)

show_matrix(H_transmon, "Transmon Hamiltonian")

# Eigenvalues
eigvals = np.linalg.eigvalsh(H_transmon)
print("Energy levels (MHz):")
for i, E in enumerate(eigvals):
    print(f"  |{i}>: {E:.1f}")

print(f"\n0-1 transition: {eigvals[1] - eigvals[0]:.1f} MHz")
print(f"1-2 transition: {eigvals[2] - eigvals[1]:.1f} MHz")
print(f"Anharmonicity:  {(eigvals[2] - eigvals[1]) - (eigvals[1] - eigvals[0]):.1f} MHz")

# Transmon Hamiltonian with anharmonicity print("=" * 50) print("TRANSMON HAMILTONIAN") print("=" * 50) # Parameters omega_q = 5000.0 # Qubit frequency (MHz) alpha = -200.0 # Anharmonicity (MHz), negative for transmon # H = omega_q * n + (alpha/2) * n * (n - I) I3 = np.eye(3, dtype=np.complex128) H_transmon = omega_q * n + (alpha / 2) * n @ (n - I3) show_matrix(H_transmon, "Transmon Hamiltonian") # Eigenvalues eigvals = np.linalg.eigvalsh(H_transmon) print("Energy levels (MHz):") for i, E in enumerate(eigvals): print(f" |{i}>: {E:.1f}") print(f"\n0-1 transition: {eigvals[1] - eigvals[0]:.1f} MHz") print(f"1-2 transition: {eigvals[2] - eigvals[1]:.1f} MHz") print(f"Anharmonicity: {(eigvals[2] - eigvals[1]) - (eigvals[1] - eigvals[0]):.1f} MHz")

In [ ]:

# Control Hamiltonian for transmon (capacitive drive)
print("=" * 50)
print("TRANSMON CONTROL")
print("=" * 50)

# Drive: proportional to (a + a_dag) for X and i*(a_dag - a) for Y
H_drive_x = a + a_dag
H_drive_y = 1j * (a_dag - a)

show_matrix(H_drive_x, "X drive: a + a_dag")
show_matrix(H_drive_y.real, "Y drive: i*(a_dag - a) [real part]")

# Control Hamiltonian for transmon (capacitive drive) print("=" * 50) print("TRANSMON CONTROL") print("=" * 50) # Drive: proportional to (a + a_dag) for X and i*(a_dag - a) for Y H_drive_x = a + a_dag H_drive_y = 1j * (a_dag - a) show_matrix(H_drive_x, "X drive: a + a_dag") show_matrix(H_drive_y.real, "Y drive: i*(a_dag - a) [real part]")

In [ ]:

# Define target X gate in 3-level space (should only affect |0> <-> |1>)
print("=" * 50)
print("X GATE IN 3-LEVEL SPACE")
print("=" * 50)

# X gate that swaps |0> and |1> but leaves |2> unchanged
X_gate_3level = np.array([
    [0, 1, 0],
    [1, 0, 0],
    [0, 0, 1],
], dtype=np.complex128)

show_matrix(X_gate_3level, "X gate (3-level)")

# Verify it's unitary
print(f"Is unitary: {np.allclose(X_gate_3level @ X_gate_3level.conj().T, I3)}")

# Define target X gate in 3-level space (should only affect |0> <-> |1>) print("=" * 50) print("X GATE IN 3-LEVEL SPACE") print("=" * 50) # X gate that swaps |0> and |1> but leaves |2> unchanged X_gate_3level = np.array([ [0, 1, 0], [1, 0, 0], [0, 0, 1], ], dtype=np.complex128) show_matrix(X_gate_3level, "X gate (3-level)") # Verify it's unitary print(f"Is unitary: {np.allclose(X_gate_3level @ X_gate_3level.conj().T, I3)}")

Note: The current GRAPE implementation uses 2^n dimensional Hilbert spaces. For 3-level systems, you would need to extend the optimizer. This is shown here for educational purposes - the concept of leakage and how to model it.

---

Part 9: Practical Example - Calibrated Qubit¶

Combining calibration data with custom Hamiltonians for realistic pulses.

In [ ]:

# Simulating a calibrated qubit scenario
print("=" * 50)
print("CALIBRATED QUBIT EXAMPLE")
print("=" * 50)

# Calibration parameters (would come from measurement)
calibration = {
    'qubit_frequency_ghz': 5.123,      # Measured qubit frequency
    'drive_frequency_ghz': 5.120,      # Drive frequency (slightly detuned)
    'anharmonicity_mhz': -210.5,       # Measured anharmonicity  
    'max_amplitude_mhz': 50.0,         # Hardware amplitude limit
    't1_us': 80.0,                     # T1 relaxation time
    't2_us': 60.0,                     # T2 dephasing time
}

print("Calibration parameters:")
for key, value in calibration.items():
    print(f"  {key}: {value}")

# Simulating a calibrated qubit scenario print("=" * 50) print("CALIBRATED QUBIT EXAMPLE") print("=" * 50) # Calibration parameters (would come from measurement) calibration = { 'qubit_frequency_ghz': 5.123, # Measured qubit frequency 'drive_frequency_ghz': 5.120, # Drive frequency (slightly detuned) 'anharmonicity_mhz': -210.5, # Measured anharmonicity 'max_amplitude_mhz': 50.0, # Hardware amplitude limit 't1_us': 80.0, # T1 relaxation time 't2_us': 60.0, # T2 dephasing time } print("Calibration parameters:") for key, value in calibration.items(): print(f" {key}: {value}")

In [ ]:

# Calculate detuning in rotating frame
detuning_mhz = (calibration['qubit_frequency_ghz'] - calibration['drive_frequency_ghz']) * 1000
print(f"\nDetuning in rotating frame: {detuning_mhz:.1f} MHz")

# Build drift Hamiltonian (in rotating frame, just the detuning)
H_drift_calibrated = detuning_mhz * PAULI_Z / 2

# Optimize X gate with calibrated parameters
config_calibrated = GrapeConfig(
    num_time_steps=100,
    duration_ns=20.0,
    target_fidelity=0.9999,
    max_iterations=500,
    max_amplitude=calibration['max_amplitude_mhz'],
    random_seed=42,
)

optimizer_cal = GrapeOptimizer(config_calibrated)
target_x = get_target_unitary("X")

result_calibrated = optimizer_cal.optimize(
    target_x, 
    num_qubits=1,
    drift_hamiltonian=H_drift_calibrated,
)

print(f"\nOptimization with calibrated parameters:")
print(f"  Fidelity: {result_calibrated.fidelity:.6f}")
print(f"  Iterations: {result_calibrated.iterations}")
print(f"  Max I amplitude: {np.max(np.abs(result_calibrated.i_envelope)):.1f} MHz")
print(f"  Max Q amplitude: {np.max(np.abs(result_calibrated.q_envelope)):.1f} MHz")

# Calculate detuning in rotating frame detuning_mhz = (calibration['qubit_frequency_ghz'] - calibration['drive_frequency_ghz']) * 1000 print(f"\nDetuning in rotating frame: {detuning_mhz:.1f} MHz") # Build drift Hamiltonian (in rotating frame, just the detuning) H_drift_calibrated = detuning_mhz * PAULI_Z / 2 # Optimize X gate with calibrated parameters config_calibrated = GrapeConfig( num_time_steps=100, duration_ns=20.0, target_fidelity=0.9999, max_iterations=500, max_amplitude=calibration['max_amplitude_mhz'], random_seed=42, ) optimizer_cal = GrapeOptimizer(config_calibrated) target_x = get_target_unitary("X") result_calibrated = optimizer_cal.optimize( target_x, num_qubits=1, drift_hamiltonian=H_drift_calibrated, ) print(f"\nOptimization with calibrated parameters:") print(f" Fidelity: {result_calibrated.fidelity:.6f}") print(f" Iterations: {result_calibrated.iterations}") print(f" Max I amplitude: {np.max(np.abs(result_calibrated.i_envelope)):.1f} MHz") print(f" Max Q amplitude: {np.max(np.abs(result_calibrated.q_envelope)):.1f} MHz")

In [ ]:

# Visualize the calibrated pulse
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

t = np.linspace(0, 20, len(result_calibrated.i_envelope))

# Time domain
axes[0].plot(t, result_calibrated.i_envelope, 'b-', label='I', linewidth=1.5)
axes[0].plot(t, result_calibrated.q_envelope, 'r-', label='Q', linewidth=1.5)
axes[0].axhline(y=calibration['max_amplitude_mhz'], color='gray', linestyle='--', alpha=0.5)
axes[0].axhline(y=-calibration['max_amplitude_mhz'], color='gray', linestyle='--', alpha=0.5)
axes[0].set_xlabel('Time (ns)')
axes[0].set_ylabel('Amplitude (MHz)')
axes[0].set_title(f'X Gate Pulse (detuning = {detuning_mhz:.1f} MHz)')
axes[0].legend()

# IQ plane
axes[1].plot(result_calibrated.i_envelope, result_calibrated.q_envelope, 'purple', linewidth=0.8)
axes[1].scatter(result_calibrated.i_envelope[0], result_calibrated.q_envelope[0], 
               color='green', s=100, label='Start', zorder=5)
axes[1].scatter(result_calibrated.i_envelope[-1], result_calibrated.q_envelope[-1], 
               color='red', s=100, label='End', zorder=5)
axes[1].set_xlabel('I Amplitude (MHz)')
axes[1].set_ylabel('Q Amplitude (MHz)')
axes[1].set_title('IQ Plane Trajectory')
axes[1].legend()
axes[1].axis('equal')

plt.tight_layout()
plt.show()

# Visualize the calibrated pulse fig, axes = plt.subplots(1, 2, figsize=(14, 5)) t = np.linspace(0, 20, len(result_calibrated.i_envelope)) # Time domain axes[0].plot(t, result_calibrated.i_envelope, 'b-', label='I', linewidth=1.5) axes[0].plot(t, result_calibrated.q_envelope, 'r-', label='Q', linewidth=1.5) axes[0].axhline(y=calibration['max_amplitude_mhz'], color='gray', linestyle='--', alpha=0.5) axes[0].axhline(y=-calibration['max_amplitude_mhz'], color='gray', linestyle='--', alpha=0.5) axes[0].set_xlabel('Time (ns)') axes[0].set_ylabel('Amplitude (MHz)') axes[0].set_title(f'X Gate Pulse (detuning = {detuning_mhz:.1f} MHz)') axes[0].legend() # IQ plane axes[1].plot(result_calibrated.i_envelope, result_calibrated.q_envelope, 'purple', linewidth=0.8) axes[1].scatter(result_calibrated.i_envelope[0], result_calibrated.q_envelope[0], color='green', s=100, label='Start', zorder=5) axes[1].scatter(result_calibrated.i_envelope[-1], result_calibrated.q_envelope[-1], color='red', s=100, label='End', zorder=5) axes[1].set_xlabel('I Amplitude (MHz)') axes[1].set_ylabel('Q Amplitude (MHz)') axes[1].set_title('IQ Plane Trajectory') axes[1].legend() axes[1].axis('equal') plt.tight_layout() plt.show()

---

Summary¶

In this notebook, we covered:

1. Pauli Matrices - The building blocks of qubit Hamiltonians

2. Pauli String Notation - Compact specification like "0.5 * X0 + Z0 Z1"

3. Programmatic Construction - Using tensor products and NumPy

4. Multi-Qubit Systems - Embedding operators in larger Hilbert spaces

5. Custom Unitaries - Rotation gates and arbitrary targets

6. Drift Hamiltonians - Including detuning and coupling in optimization

7. Two-Qubit Gates - CNOT, CZ, and other entangling gates

8. Three-Level Systems - Modeling transmon anharmonicity

9. Practical Example - Using calibration data for realistic pulses

Key Takeaways¶

• Use Pauli strings for quick Hamiltonian specification
• Use programmatic construction for complex or parameterized systems
• Always include drift Hamiltonians for realistic pulse design
• Two-qubit gates require longer durations and more iterations
• Consider leakage to higher levels for high-fidelity transmon control

Next Steps¶

• Integrate with real hardware using the HAL client
• Use CalibrationLoader to fetch measured qubit parameters
• Explore DRAG pulses for leakage suppression